Simpson's rule area and volume integrals using buildq & lambda methods
Subject: Simpson's rule area and volume integrals using buildq & lambda methods
From: Edwin Woollett
Date: Fri, 9 Dec 2011 15:26:19 -0800
On Dec. 9, I wrote:
-------------------------------------
>Since I am a complete novice using buildq with lambda methods,
>there are surely more efficient ways to do this, which I would be
>glad to hear about.
------------------------------------
An obvious thing to try is to call area when trying
to invoke volume, although it is not obvious that this
is more efficient than the present volume.
volume1 calls area, but I need to worry about interference
between some of the symbols already used by area.
Also, this last work showed that I need to call the
one dimensional code simp1 (say), instead of simp,
which Maxima tried to use as an evaluation flag.
-----------------------------------------------------------------
simp1(expr,var,a,b,n):=
block ([h : abs (b-a)/n,jj,
sf : buildq([expr,var],lambda([var],expr)) ],
if not evenp(n) then return("n must be an even integer"),
float((h/3)*( sf(a) + sf(b) +
4*sum( sf(a + (2*jj-1)*h), jj, 1, n/2 ) +
2*sum( sf(a + 2*jj*h), jj, 1, (n/2) - 1 ) )))$
/* area = same and volume = same */
volume1(ex1,u1L,v1L,w1L,n1L) :=
block([%%u1,u11,u21,%%v1,v11,v21,%%w1,w11,w21,nu1,nv1,nw1,%uu1,
exu : buildq([%%u1,ex1,v11,v21,%%v1,nv1,w11,
w21,%%w1,nw1 ],
lambda([u01],
block([%exuv,%v11,%v21,%w11,%w21],
%exuv:subst(u01,%%u1,ex1),
%v11:subst(u01,%%u1,v11),
%v21:subst(u01,%%u1,v21),
%w11:subst(u01,%%u1,w11),
%w21:subst(u01,%%u1,w21),
area(%exuv,[%%v1,%v11,%v21],[%%w1,%w11,%w21],[nv1,nw1]))))],
[%%u1,u11,u21] : u1L,
[%%v1,v11,v21] : v1L,
[%%w1,w11,w21] : w1L,
[nu1,nv1,nw1] : n1L,
simp1(exu(%uu1),%uu1,u11,u21,nu1))$
-------------------------------------------------
/*
(%i1) load(nintmd);
(%o1) "c:/work2/nintmd.mac"
(%i2) float(integrate(integrate(integrate(x*y^2*z^3,z,0,1),y,0,1),x,0,1));
(%o2) 0.041666666666667
(%i3) volume(x*y^2*z^3,[x,0,1],[y,0,1],[z,0,1],[4,4,4]);
(%o3) 0.041666666666667
(%i4) volume1(x*y^2*z^3,[x,0,1],[y,0,1],[z,0,1],[4,4,4]);
(%o4) 0.041666666666667
(%i5) float(integrate(integrate(integrate(x*y^2*z^3,z,0,y),y,0,x),x,0,1));
(%o5) 0.003968253968254
(%i6) volume(x*y^2*z^3,[x,0,1],[y,0,x],[z,0,y],[16,16,16]);
(%o6) 0.003969543636397
(%i7) volume1(x*y^2*z^3,[x,0,1],[y,0,x],[z,0,y],[16,16,16]);
(%o7) 0.003969543636397
*/
--------------------------------
Ted Woollett