On Thu, Oct 6, 2011 at 2:25 PM, Edwin Woollett <woollett at charter.net> wrote:
>
> 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? I think it just means that quadpack was rather conservative in its
estimate of what the relative error actually was and the limit was too low
for quadpack to refine the estimate any further.
But what if you didn't know the analytical solution and therefore the "true"
error? All you would have is the estimated value. Then what?
Also, how about this:
quad_qags(g,x,2,3)->
[- 1.6023580453445135e-4, .005299669215071409, 8379, 1]
The actual relative error is 1d-2, far above the requested default relative
error.
> 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.
>
I didn't say it would be easy to automate this. :-)
>
> 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.
>
For some quick and dirty hack it's ok.
Ray