>I think the error message is appropriate. You've requested a relerr of
>1d-8, but the actual relerr is about 1d-5. What >did you expect qags to do
>in this case?
There are three issues here.
1. I initially accepted the default nominal rel error goal
which ignored both the approximate value of the
expected result and the number of significant figures
desired.
2. the actual rel error of slatec returned results is usually
of the order of 10^(-6) times the requested rel error precision
(ie, much more accurate).
3. given a requested precision, one can avoid the slatec error
message by increasing the limit value.
-------------------------------------------------
here we calculate the true value first using integrate and float.
(%i1) (display2d:false,fpprintprec:6)$
(%i2) g : (x-2)^2 * sin (4000*x)$
(%i3) integrate(g,x,2,3);
(%o3) (4000*sin(12000)-7999999*cos(12000))/32000000000-cos(8000)/32000000000
(%i4) tval : float(%);
(%o4) -1.586246E-4
proceeding very roughly,
for six figure precision, we need abs error <= 10^(-10)
or rel error <= (abs err)/ value = 6.7 10^(-7) =appr 10^(-6)
so use quad_qags directly:
(%i5) qval : quad_qags(g,x,2,3,epsrel=1.0e-6,limit=400)[1];
***MESSAGE FROM ROUTINE DQAGS IN LIBRARY SLATEC.
***INFORMATIVE MESSAGE, PROG CONTINUES, TRACEBACK REQUESTED
* ABNORMAL RETURN
* ERROR NUMBER = 1
*
***END OF MESSAGE
(%o5) -1.586246E-4 <---- answer
(%i6) abs(qval - tval);
(%o6) 7.375285E-17 <------abs error
(%i7) %/abs(tval);
(%o7) 4.64952E-13 <------rel error - note how small
so increase limit value:
(%i8) qval : quad_qags(g,x,2,3,epsrel=1.0e-6,limit=600)[1];
no message from slatec
(%o8) -1.586246E-4 <---- answer
(%i9) abs(qval-tval);
(%o9) 5.838429E-17 <-------abs error
(%i10) %/abs(tval);
(%o10) 3.680657E-13 <------- rel error note - how small
likewise with nint:
(%i11) load(nint);
(%o11) "c:/work2/nint.mac"
(%i12) nint(g,x,2,3,epsrel=1.0e-6,limit=600);
(%o12) -1.586246E-4
The actual relative error of slatec returned result
is typically 10^(-6) to 10^(-7) smaller than
the requested value, so we can get good results
using here epsrel=1.0e-2 (but still have to
increase limit to avoid slatec message 1).
(%i14) nint(g,x,2,3,epsrel=1.0e-2,limit=400);
(%o14) -1.586245E-4
---------------------------------------------------
>
>For such highly oscillatory integrals some other approach could be used
>like quad_qawo. It says:
>
>quad_qawo((x-2)^2,x, 2,3,4000,'sin);
> [- 1.586246466526855e-4, 0.0, 25, 0]
>
>As accurate, and much less expensive too! Double win!
---------------------------------------------------------------
I completely agree a double win, but a lot more work to
incorporate this efficiency bonus into a wrapper for quadpack
routines.
1. The nintegrate function must detect the highly oscillatory nature
of the integrand.
2. Only sometimes is the oscillation due simply to a factor of
sin(...something ), which would also have to be identified
to include the correct sin or cos expected by quad_qawo
which is designed to calculate fourier series coefficients.
3. The only quadpack function really suitable for highly
oscillating integrands in general is quad_qag, in which
you have to supply the 'key' as an additional argument,
a number from 1 to 6 to select the order of the method.
So using quad_qag, you need to both detect rapid oscillations,
and decide on the 'key' value.
My initial approach is to sacrifice efficiency for a reasonably
accurate result for down and dirty use, by a student or engineer,
who often is not interested in high precision, but in four or six
figure accuracy.
One approach is to have the user input a requested rel error
precision goal (if not the default), and inside nint we could
use our knowledge that usually quad_qags returns a result
that has a lot more rel error precision than the nominal
value input to quad_qags.
Thanks for your feedback.
Ted