Subject: how to "fool" Gauss-Kronrod 21-point rule
From: Edwin Woollett
Date: Fri, 13 Apr 2012 14:45:26 -0700
On April 13, 2012, Raymond Toy wrote:
-------------------------------
>Not exactly sure what you're looking for but here is one example where
>quad_qags gives no error or warning messages, but the result is totally
>wrong:
>
> quad_qags(4^(-16)/((x - %pi/4)^2 + 16^(-16)),x,0,1);
> [1.4602196559157534e-6, 1.8681728113217153e-6, 693, 5]
-----------------------------------------------
Thanks for the example, which my present quadpack wrapper
does successfully:
-----------------------------------------------------
(%i1) load(quad);
load nint_util.mac
load quad1d.mac
load quad2d.mac
(%o1) "c:/work2/quad.mac"
(%i2) quad(4^(-16)/((x - %pi/4)^2 + 16^(-16)),x,0,1);
(%o2) 3.1415926
(%i3) nargL;
(%o3) [qag,1/(4294967296*((x-%pi/4)^2+1/18446744073709551616)),x,0.0,1.0,3,
limit = 800]
(%i4) noutL;
(%o4) [qag,3.1415926,9.28343605E-9,2201,0]
(%i5) mdefint(4^(-16)/((x - %pi/4)^2 + 16^(-16)),x,0,1);
(%o5) atan(1073741824*%pi)-atan(1073741824*%pi-4294967296)
(%i6) float(%);
(%o6) 3.1415927
---------------------------------------------------
I am getting to the point of looking for more obscure
examples for testing of the code.
The call to quad above (with no general options)
sends the expression to a competition between
qags and qag. For each method the epsrel criterion
is tried, and if it doesn't work, the epsabs criterion
is tried. For this case quad_qag with the epsrel
criterion (the default) succeeds, and neither of the
quad_qags tries succeeds.
Ted