[Maxima] simple block for romberg method

[Luigi Marino]

There is a simple block
for romberg method:

trap(n,a,b,f):=block(
h:(b-a)/n,
s:0,
for i:1 thru n-1 step 1 do
(s:s+f(a+i*h)),
AV:(h/2)*(f(a)+2*s+f(b)),
return (float(AV)))$

romb(max_iter,a,b,f):=block(
integ[1,1]:trap(1,a,b,f),
  for i:1 thru max_iter do
    (n:2^i,
    integ[i+1,1]:trap(n,a,b,f),
      for j:2 thru i+1 do
        (k:2+i-j,
        integ[k,j]:(4^(j-1)*integ[k+1,j-1]-integ[k,j-1])/(4^(j-1)-1))),
  integ_value:integ[1,max_iter],
return(integ_value))$

Example:

f(x):=cos(x);

romb(10,0,10,f);
-0.54402111088937

f(x):=sqrt(1+x^3);

romb(12,-1,1,f);
1.952753609576008

Best wishes
Luigi Marino
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20100402/55875e92/attachment-0001.html>