simple block for romberg method



By the way: The hashed arrays in the body of your function makes your code easy to break:

 (%i3) integ[a] : 'b$
 (%i4) f(x):=cos(x);
 (%o4) f(x):=cos(x)
 (%i5) romb(10,0,10,f);
  Array integ has dimension 1; it was called by integ[1,1] #0: romb(max_iter=10,a=0,b=10,f=f)

--Barton

-----maxima-bounces at math.utexas.edu wrote: -----


>you?should?use?apply(f,?...)??or?else?the?function?f?defined?globally?will
>be?used?instead?of?the?f,?the?parameter?to?trap.

>-----?Original?Message?-----

>>?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
>>?_______________________________________________
>>?Maxima?mailing?list
>>?Maxima at math.utexas.edu
>>?http://www.math.utexas.edu/mailman/listinfo/maxima
>_______________________________________________
>Maxima?mailing?list
>Maxima at math.utexas.edu
>http://www.math.utexas.edu/mailman/listinfo/maxima