Subject: how to "fool" Gauss-Kronrod 21-point rule
From: Edwin Woollett
Date: Sat, 14 Apr 2012 09:59:07 -0700
On April 13, 2012, Raymond Toy wrote:
------------------------------
>It's a numerical issue. e1 is clearly more accurate when x is near %pi/4.
>Compare the value of e1 and e2 for x =.78539816. (If this doesn't cause a
>problem, add a few for digits from %pi/4.).
--------------------------------------
Yes, that is the problem, and I should have looked at that difference in
accuracy of the two expressions.
---------------------------------------------------
(%i1) load(quad);
(%o1) "c:/work2/quad.mac"
(%i2) e1 : 4^(-16)/((x - %pi/4)^2 + 16^(-16));
(%o2) 1/(4294967296*((x-%pi/4)^2+1/18446744073709551616))
(%i3) e2 : expand(e1);
(%o3) 1/(4294967296*x^2-2147483648*%pi*x+268435456*%pi^2+1/4294967296)
(%i4) xL : [0.785,0.7853,0.78539,0.785398]$
(%i5) for xv in xL do
print(" ",xv," ",fcompare(e1,e2,x,xv,20))$
0.785 [1.53645769E-13, 6.40931986E-10]
0.7853 [6.23703981E-13, 2.47483209E-9]
0.78539 [7.500549E-12, 4.79013095E-7]
0.785398 [3.74731578E-10, 0.00215256]
-----------------------------------------------------------
where fcompare is:
----------------------------------
fcompare (ee1,ee2,xx,xval,fppval):=
block([fpprec:fppval, tval, e1f, e2f, relerr1, relerr2],
tval : bfloat (subst (xx = xval, ee1)),
e1f : float (subst (xx = xval, ee1)),
e2f : float (subst (xx = xval, ee2)),
relerr1 : float (abs ((e1f - tval)/tval)),
relerr2 : float (abs ((e2f - tval)/tval)),
[relerr1, relerr2])$
---------------------------------------------
Ted