Question

Write a program to implement the Romberg Algorithm.

Input: f(x), an interval [a,b], N (the number of subintervals on the finest grid on level 0 is 2^N, therefore, N is usualy a small integer)

Output: the triangle generated by Romberg Algorithm.

Answer 1

Parameters:
f                - function to be integrated
a, b           - interval of integration
N              - number of columns to generated (zero-based)
print_table - true or false, true to display generated table

romberg := proc(f::algebraic, a, b, N,print_table)
local R,h,k,row,col;
R := array(0..N,0..N);

# Compute column 0, Trapezoid Rule approximations of
#                   1,2,4,8,..2^N subintervals
h := evalf(b - a);
R[0,0] := evalf(h/2 * (f(a)+f(b)));
for row from 1 to N do;
   h := h/2;
   R[row,0] := evalf(0.5*R[row-1,0] +
                     sum(h*f(a+(2*k-1)*h),k=1..2^(row-1)));
   # Compute [row,1]:[row,row], via Richardson extrapolation
   for col from 1 to row do;
     R[row,col] := ((4^col)*R[row,col-1] - R[row-1,col-1]) /
                       (4^col - 1);
   end do;
end do;

# Display results if requested
if (print_table) then
   for row from 0 to N do;
     for col from 0 to row do;
       printf("%12.10f ",R[row,col]);
     end do;
     printf("\n");
   end do;
end if;

return(R[N,N]);

end proc:

Ask if you have any queries. Thank you