[Maxima] romberg or numerical integration

Jim FitzSimons cherry at getnet.net
Fri Jun 30 12:00:03 CDT 2006

```This is a multi-part message in MIME format.

------=_NextPart_000_0000_01C69C2B.F6030E40
Content-Type: text/plain;
charset="us-ascii"
Content-Transfer-Encoding: 7bit

How does Maxima best do numerical integration?
I tried romberg and I failed.
Attached is an example.

Maxima 5.9.3 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1)
batching #pC:/dfw5/MATH/testrf.mac
(%i3) mrf : bfloat([[1, 2, 0, 1.3110287771461],
[%i, - %i, 0, 1.8540746773014], [0.5, 1, 0, 1.8540746773014],
[%i - 1, %i, 0, 0.79612586584234 - %i 1.2138566698365],
[2, 3, 4, 0.58408284167715], [%i, - %i, 2, 1.0441445654064],
[%i - 1, %i, 1 - %i, 0.93912050218619 - %i 0.53296252018635]])
Warning:  Float to bigfloat conversion of 1.3110287771461
Warning:  Float to bigfloat conversion of 1.8540746773014001
Warning:  Float to bigfloat conversion of 0.5
Warning:  Float to bigfloat conversion of 1.8540746773014001
Warning:  Float to bigfloat conversion of 0.79612586584234002
Warning:  Float to bigfloat conversion of -1.2138566698365001
Warning:  Float to bigfloat conversion of 0.58408284167714997
Warning:  Float to bigfloat conversion of 1.0441445654064001
Warning:  Float to bigfloat conversion of 0.93912050218618992
Warning:  Float to bigfloat conversion of -0.53296252018635004
(%o3) [[1.0b0, 2.0b0, 0.0b0, 1.3110287771461b0],
[%i, - 1.0b0 %i, 0.0b0, 1.8540746773014b0],
[5.0b-1, 1.0b0, 0.0b0, 1.8540746773014b0],
[%i - 1.0b0, %i, 0.0b0, 7.9612586584234b-1 - 1.2138566698365b0 %i],
[2.0b0, 3.0b0, 4.0b0, 5.840828416771499b-1],
[%i, - 1.0b0 %i, 2.0b0, 1.0441445654064b0],
[%i - 1.0b0, %i, 1.0b0 - 1.0b0 %i, 9.391205021861899b-1
- 5.3296252018635b-1 %i]]
(%i4) 				  fpprec : 20
(%o4) 				      20
(%i5) 			       bferrtol : 1.0b-6
(%o5) 				    1.0b-6
(%i6) 				   maxn : 16
(%o6) 				      16
(%i7) 		   bfrf((mrf ) , (mrf ) , (mrf ) ) - (mrf )
5 1	     5 2      5 3	 5 4
(%o7) 			   1.7679644369588548569b-15
2
(%i8) rf(x, y, z) := romberg(u/SQRT((x - 2 u x + u  (1 + x))
2			  2
(y - 2 u y + u  (1 + y)) (z - 2 u z + u  (1 + z))), u, 0,
1)
2
(%o8) rf(x, y, z) := romberg(u/SQRT((x - 2 u x + u  (1 + x))
2			  2
(y - 2 u y + u  (1 + y)) (z - 2 u z + u  (1 + z))), u, 0,
1)
(%i9) 	        bfloat(rf((mrf ) , (mrf ) , (mrf ) ) - (mrf ) )
5 1      5 2      5 3	   5 4
Warning:  Float to bigfloat conversion of 1.0
Warning:  Float to bigfloat conversion of 2.0
Warning:  Float to bigfloat conversion of 1.0
Warning:  Float to bigfloat conversion of 2.0
Warning:  Float to bigfloat conversion of 1.0
Warning:  Float to bigfloat conversion of 2.0
Maxima encountered a Lisp error:

Error in MACSYMA-TOP-LEVEL [or a callee]: ((MTIMES SIMP) 1.0
((MEXPT SIMP)
((|\$sqrt| SIMP
(9
"C:/dfw5/MATH/testrf.mac"
SRC \$RF 9))
((BIGFLOAT SIMP 69)
295147905179352825856 1))
-1)) is not of type (OR
RATIONAL

LISP:FLOAT).

Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

Regards, Jim FitzSimons

------=_NextPart_000_0000_01C69C2B.F6030E40
Content-Type: application/octet-stream;
name="elliptbf2.mac"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="elliptbf2.mac"

/* elliptbf.mac jim fitzsimons 25 january, 1998
for complex numbers */

bfrc(x, y):=3D block([xn,yn,z,w,a,an,pwr4,n,epslon,lamda,sn,s],
yn:expand(bfloat(y)),
if ((imagpart(yn) =3D 0) and (realpart(yn) < 0.b0)) then
(xn:bfloat(x-y),
yn:-yn,
z:yn,
w:sqrt(bfloat(x)/xn))
else
(xn:bfloat(x),
z:yn,
w:1.b0),

a:bfloat((xn+yn+yn)/3.b0),
epslon:cabs(a-xn)/bferrtol,
an:a,
pwr4:1.b0,
n:0,
while (epslon*pwr4 > cabs(an)) do block([],
pwr4:pwr4/4.b0,
lamda:expand(bfloat(2.b0*sqrt(xn)*sqrt(yn)+yn)),
an:expand(bfloat((an+lamda)/4.b0)),
xn:expand(bfloat((xn+lamda)/4.b0)),
yn:expand(bfloat((yn+lamda)/4.b0)),
n:n+1),
/* c2=3D3/10,c3=3D1/7,c4=3D3/8,c5=3D9/22,c6=3D159/208,c7=3D9/8 */
sn:(z-a)*pwr4/an,
s:sn*sn*(3.b0/10.b0+sn*(1.b0/7.b0+sn*(3.b0/8.b0+
sn*(9.b0/22.b0+sn*(159.b0/208.b0+9.b0*sn/8.b0))))),
rectform(expand(bfloat(w*(1.b0+s)/sqrt(an)))) )\$

bfrd(x, y, z):=3D block([xn,yn,zn,a,epslon,an,sigma,power4,n,
xnroot,ynroot,znroot,lamda,xndev,yndev,zndev,ee2,ee3,ee4,ee5,s],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
a:bfloat((xn+yn+3.b0*zn)/5.b0),
epslon:max(cabs(a-xn),cabs(a-yn),cabs(a-zn))/bferrtol,
an:a,
sigma:0.b0,
power4:1.b0,
n:0,
while (power4*epslon > cabs(an)) do block([],
xnroot:bfloat(sqrt(xn)),
ynroot:bfloat(sqrt(yn)),
znroot:bfloat(sqrt(zn)),
lamda:expand(xnroot*ynroot+xnroot*znroot+ynroot*znroot),
sigma:expand(bfloat(sigma+power4/(znroot*(zn+lamda)))),
power4:power4/4.b0,
xn:expand(bfloat((xn+lamda)/4.b0)),
yn:expand(bfloat((yn+lamda)/4.b0)),
zn:expand(bfloat((zn+lamda)/4.b0)),
an:expand(bfloat((an+lamda)/4.b0)),
n:n+1),
/* =
c1=3D-3/14,c2=3D1/6,c3=3D9/88,c4=3D9/22,c5=3D-3/22,c6=3D-9/52,c7=3D3/26 =
*/
xndev:(a-bfloat(x))*power4/an,
yndev:(a-bfloat(y))*power4/an,
zndev:(-xndev-yndev)/3.b0,
ee2:xndev*yndev-6.b0*zndev*zndev,
ee3:(3.b0*xndev*yndev-8.b0*zndev*zndev)*zndev,
ee4:3.b0*(xndev*yndev-zndev*zndev)*zndev*zndev,
ee5:xndev*yndev*zndev*zndev*zndev,
/* s:ee2*(-3.b0/14.b0+9.b0*ee2/88.b0-9.b0*ee3/52.b0),
s2:ee3/6.b0-3.b0*ee4/22.b0+3.b0*ee5/26.b0, */
s:1.b0-3.b0*ee2/14.b0+ee3/6.b0+9.b0*ee2*ee2/88.b0-3.b0*ee4/22.b0
-9.b0*ee2*ee3/52.b0+3.b0*ee5/26.b0
-1.b0*ee2*ee2*ee2/16.b0+3.b0*ee3*ee3/10.b0+3.b0*ee2*ee4/20.b0
+45.b0*ee2*ee2*ee3/272.b0-9.b0*(ee2*ee5+ee3*ee4)/68.b0,
rectform(expand(bfloat(3.b0*sigma+power4*s/(an^(3/2))))) )\$

bfrf(x, y, z):=3D block([xn,yn,zn,a,epslon,an,pwr4,xnroot,ynroot,znroot,
lamda,xndev,yndev,zndev,ee2,ee3,s,n],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
a:bfloat((xn+yn+zn)/3.b0),
epslon:max(cabs(a-xn),cabs(a-yn),cabs(a-zn))/bferrtol,
an:a,
pwr4:1.b0,
n:0,
while ((epslon*pwr4 > cabs(an)) and (n < maxn)) do block([],
xnroot:bfloat(sqrt(xn)),
ynroot:bfloat(sqrt(yn)),
znroot:bfloat(sqrt(zn)),
lamda:expand(xnroot*ynroot+xnroot*znroot+ynroot*znroot),
xn:expand(bfloat((xn+lamda)/4.b0)),
yn:expand(bfloat((yn+lamda)/4.b0)),
zn:expand(bfloat((zn+lamda)/4.b0)),
an:expand(bfloat((an+lamda)/4.b0)),
pwr4:pwr4/4.b0,
n:n+1),
/* c1=3D-1/10,c2=3D1/14,c3=3D1/24,c4=3D-3/44 */
xndev:(a-bfloat(x))*pwr4/an,
yndev:(a-bfloat(y))*pwr4/an,
zndev:-xndev-yndev,
ee2:xndev*yndev-zndev*zndev,
ee3:xndev*yndev*zndev,
s:1.b0-1.b0*ee2/10.b0+ee3/14.b0+ee2*ee2/24.b0-3.b0*ee2*ee3/44.b0,
rectform(expand(bfloat(s/sqrt(an)))) )\$

bfrj(x, y, z, p):=3D block([xn,yn,zn,pn,v,s],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
qn:bfloat(-p),
if ((imagpart(xn) =3D 0) and (realpart(xn) >=3D 0.b0) and
(imagpart(yn) =3D 0) and (realpart(yn) >=3D 0.b0) and
(imagpart(zn) =3D 0) and (realpart(zn) >=3D 0.b0) and
(imagpart(qn) =3D 0) and (realpart(qn) > 0.b0)) then
(v:sort([xn,yn,zn]),
xn:v[1],yn:v[2],zn:v[3],
pn:yn+(zn-yn)*(yn-xn)/(yn+qn),
s:(pn-yn)*bfrj1(xn,yn,zn,pn)-3.b0*bfrf(xn,yn,zn),
s:s+3.b0*sqrt(xn*yn*zn/(xn*zn+pn*qn))*bfrc(xn*zn+pn*qn,pn*qn),  =20
rectform(expand(bfloat(s/(yn+qn)))) )  =20
else
bfrj1(x, y, z, p))\$

bfrj1(x, y, z, p):=3D block([xn,yn,zn,pn,sigma,power4,k,a,epslon,an,
xnroot,ynroot,znroot,pnroot,lamda,dn,en,xndev,yndev,zndev,pndev,
ee2,ee3,ee4,ee5,s],
xn:bfloat(x),
yn:bfloat(y),
zn:bfloat(z),
pn:bfloat(p),
en:expand(bfloat((pn-xn)*(pn-yn)*(pn-zn))),
sigma:0.b0,
power4:1.b0,
k:0,
a:bfloat((xn+yn+zn+pn+pn)/5.b0),
epslon:max(cabs(a-xn),cabs(a-yn),cabs(a-zn),cabs(a-pn))/bferrtol,
an:a,
while (power4*epslon > cabs(an)) do block([],
xnroot:bfloat(sqrt(xn)),
ynroot:bfloat(sqrt(yn)),
znroot:bfloat(sqrt(zn)),
pnroot:bfloat(sqrt(pn)),
lamda:expand(xnroot*ynroot+xnroot*znroot+ynroot*znroot),
dn:expand(bfloat((pnroot+xnroot)*(pnroot+ynroot)*(pnroot+znroot))),
sigma:expand(bfloat(sigma+power4*bfrc(1.b0,1.b0+en/(dn*dn))/dn)),
power4:power4/4.b0,
en:expand(en/64.b0),
xn:expand(bfloat((xn+lamda)/4.b0)),
yn:expand(bfloat((yn+lamda)/4.b0)),
zn:expand(bfloat((zn+lamda)/4.b0)),
pn:expand(bfloat((pn+lamda)/4.b0)),
an:expand(bfloat((an+lamda)/4.b0)),
k:k+1),
/*1-3/14*ee2+1/6*ee3+9/88*ee2^2-3/22*ee4-9/52*ee2*ee3+3/26*e5
-1/16*ee2^3+3/10*ee3^2+3/20*ee2*ee4+45/272*ee2^2*ee3-9/68*(ee2*ee5+ee3*ee=
4) */
xndev:(a-bfloat(x))*power4/an,
yndev:(a-bfloat(y))*power4/an,
zndev:(a-bfloat(z))*power4/an,
pndev:(-xndev-yndev-zndev)/2.b0,
ee2:xndev*yndev+xndev*zndev+yndev*zndev-3.b0*pndev*pndev,
ee3:xndev*yndev*zndev+2.b0*ee2*pndev+4.b0*pndev*pndev*pndev,
ee4:(2.b0*xndev*yndev*zndev+ee2*pndev+3.b0*pndev*pndev*pndev)*pndev,
ee5:xndev*yndev*zndev*pndev*pndev,
/* s:1.b0-3.b0*ee2/14.b0+ee3/6.b0+9.b0*ee2*ee2/88.b0-9.b0*ee2*ee3/52.b0
+3.b0*ee5/26.b0, */=20
s:1.b0-3.b0*ee2/14.b0+ee3/6.b0+9.b0*ee2*ee2/88.b0-3.b0*ee4/22.b0
-9.b0*ee2*ee3/52.b0+3.b0*ee5/26.b0
-1.b0*ee2*ee2*ee2/16.b0+3.b0*ee3*ee3/10.b0+3.b0*ee2*ee4/20.b0
+45.b0*ee2*ee2*ee3/272.b0-9.b0*(ee2*ee5+ee3*ee4)/68.b0,
rectform(expand(bfloat(6.b0*sigma+power4*s/(an^(3/2))))) )\$

bfellf(phi, ak):=3Dblock([s,c,cc,sak,q],
s:sin(phi),
c:cos(phi),
cc:c * c,
sak:s * bfloat(ak),
q:(1.b0 - sak) * (1.b0 + sak),
rectform(bfloat(expand(s * bfrf(cc, q, 1.b0)))) )\$

bfelle(phi, ak):=3Dblock([s,c,cc,sak,q],
s:sin(phi),
c:cos(phi),
cc:c * c,
sak:s * bfloat(ak),
q:(1.b0 - sak) * (1.b0 + sak),
rectform(expand(bfloat(s*(bfrf(cc, q, 1.b0)-sak*sak*bfrd(cc, q, =
1.b0)/3.b0)))) )\$

bfellpi(phi, en, ak):=3Dblock([s,c,cc,sak,q,enss],
s:sin(phi),
enss:en * s * s,
c:cos(phi),
cc:c * c,
sak:s * bfloat(ak),
q:(1.b0 - sak) * (1.b0 + sak),
rectform(expand(bfloat(s*(bfrf(cc, q, 1.b0)-enss*bfrj(cc, q, 1.b0, 1.b0 =
+ enss)=20
/ 3.b0)))) )\$

bfam(u, m):=3Dblock([a,b,c,k:bfloat(sqrt(m)),fac,phi,i,j,dela,
ee2:bferrtol*bferrtol],
array([a,b,c],256),
a[0]:1.b0,
b[0]:bfloat(sqrt((1.b0-k)*(1.b0+k))),
c[0]:k,
if (1.b0-k) > ee2 then
(fac:1.b0,dela:1.b0,i:1,
while dela > ee2 do block([],
fac:fac*2.b0,
a[i]:(a[i-1]+b[i-1])/2.b0,
b[i]:bfloat(sqrt(a[i-1]*b[i-1])),
c[i]:(a[i-1]-b[i-1])/2.b0,
dela:abs(a[i]-a[i-1]),
i:i+1),
phi:u*fac*a[i-1],
for j:i-1 step -1 thru 1 do
(phi:(phi+asin(c[j]/a[j]*sin(phi)))/2.b0))
else
(phi:2.b0*atan(exp(u))-bfloat(%pi/2.b0)),
rectform(expand(bfloat(phi))) )\$

bferrtol:1.b-6;

maxn:16;

------=_NextPart_000_0000_01C69C2B.F6030E40
Content-Type: application/octet-stream;
name="testrf.mac"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="testrf.mac"

/* Use numerical checks from Numerical Algorithms 10(1995)13-26.
B.C. Carlson Numerical computation of real or complex elliptic integrals
page 22 3. Numerical checks */

mrf:bfloat([[1,2,0,1.3110287771461],
[%i,-%i,0,1.8540746773014],
[0.5,1,0,1.8540746773014],
[%i-1,%i,0,0.79612586584234-%i*(1.2138566698365)],
[2,3,4,0.58408284167715],[%i,-%i,2,1.0441445654064],
[%i-1,%i,1-%i,0.93912050218619-%i*(0.53296252018635)]]);

fpprec:20;
bferrtol:1.b-6;
maxn:16;

bfrf(mrf[5][1],mrf[5][2],mrf[5][3])-mrf[5][4];

rf(x,y,z):=3Dromberg =
(u/SQRT((u^2*(x+1)-2*u*x+x)*(u^2*(y+1)-2*u*y+y)*(u^2*(z+1)-2*u*z+z)),u,0,=
1);

bfloat(rf(mrf[5][1],mrf[5][2],mrf[5][3])-mrf[5][4]);

------=_NextPart_000_0000_01C69C2B.F6030E40--

```