Application Center - Maplesoft
  Submit your work
Application Center Applications Romberg Algorithm for Integration Preview

App Preview:

Romberg Algorithm for Integration

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application

 

Image 

Romberg Algorithm for Integration 

Jay Pedersen  

Explanation 

This estimates the integral value over interval [a,b] of input function f(x). 

 

Routine: romberg
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 

Algorithm 

> 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:
 

Examples 

> # Estimate the integral of e^(-x^2)/2 from [0,1]
f := x -> exp(-x^2/2);
val := romberg(f,0,1,4,false);
f := proc (x) options operator, arrow; exp(-1/2*x^2) end proc
val := .8556243918

> # Estimate integral of 2^x over [0,4], display table
f := x -> 2^x;
val := romberg(f,0,4,4,true);
f := proc (x) options operator, arrow; 2^x end proc
34.0000000000
25.0000000000 22.0000000000
22.5000000000 21.6666666700 21.6444444500
21.8566017200 21.6421356300 21.6405002300 21.6404376300
21.6945506600 21.6405336400 21.6404268400 21.6404256800 21.6404256300
val := 21.64042563
 

Restrictions 

The input function must be Riemann integrable. 

Notes 

(1) Uses the trapezoid rule (with uniform spacing) to generate
     the first column.  Uses properties of the error formula for
     the trapezoid rule integration (which is defined using the
     second derivative of the function).
(2) To check if the method is working for a function; check table
     entries (set print_table = true); to see if the following holds:  

               Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mi( ≅ Typesetting:-mrow(Typesetting:-msup(Typesetting:-mn(  

References 

>  (1) Numerical Mathematics and Computing, 5th Edition
      Ward Cheney and David Kincaid.  ISBN 0-534-38993-7.
      Section 5.3, Romberg Algorithm, Pages 223-224.
 (2) Dr. Steve From - University of Nebraska at Omaha
                                  mathematics professor; lecture notes
                                  from Numerical Methods course.

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.
 

Image