 Application Center - Maplesoft

# Romberg Algorithm for Integration

You can switch back to the summary page by clicking here. 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);  > # Estimate integral of 2^x over [0,4], display table f := x -> 2^x; val := romberg(f,0,4,4,true); 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 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:  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. 