Mathematics Source Library
C & ASM


Romberg's Method

Romberg's Method

The Euler-Maclaurin summation formula relates the integral of a function f(x) over a closed and bounded interval [a,b] , abf( x ) dx, and the composite trapezoidal rule,

Th(f) = h [ f(a)/2 + f(a+h) + ··· + f(b-h) + f(b)/2 ] ,
by
Th(f) = abf( x ) dx + (h2/12) [ f'(b) - f'(a) ] - (h4/720) [ f'''(b) - f'''(a) ]
+ ··· + K h 2p-2 [f(2p-3)(b) - f(2p-3)(a) ] + O(h 2p)
,
where f', f''', and f(2p-3) are the first, third and (p-3)rd derivatives of f and K is a constant.

If the subinterval length is halved, then
Th/2(f) = abf( x ) dx + (h2/4·12) [ f'(b) - f'(a) ] - (h4/16·720) [ f'''(b) - f'''(a) ]
+ ··· + K (h/2) 2p-2 [f(2p-3)(b) - f(2p-3)(a) ] + O(h 2p).
So that
(4Th/2(f) - Th(f))/3 = abf( x ) dx + (h4/2880) [ f'''(b) - f'''(a) ] - (h6/96768) [ f(5)(b) - f(5)(a) ]
+ ··· + K h 2p-2 [f(2p-3)(b) - f(2p-3)(a) ] + O(h 2p)
.

The term has vanished. This process can be continued, each halving of the subinterval length results in a new composite trapezoidal rule estimate of the integral which can be combined with previous estimates to yield an estimate in which the lowest order term involving h vanishes.
The easiest way to combine all the estimates from applications of the trapezoidal rule by halving the length of the subintervals is to arrange the estimates in a column from the coarsest estimate to the finest estimate. To form the second, column take two adjacent values from the first column, subtract the finer estimate from the coarser estimate divide by 3 and add to the finer estimate. The second column is automatically arranged from the coarsest estimate to the finest with one less element. Continue, form the third column by taking two adjacent values from the second column, subtract the finer estimate from the coarser estimate, divide by 15 and add to the finer estimate. The third column is automatically arranged from the coarsest to the finest estimate with one fewer element than in the second column. This process is continued until there is only one element in the last column, this is the estimate of the integral. The numbers which are used the divide the difference of two adjacent elements in the ith column is 4i - 1.

Function List

  • double Rombergs_Integration_Method( double a, double h, double tolerance, int max_cols, double (*f)(double), int *err )

    Integrate, using Romberg's method, the user supplied function (*f)(x) from a to a + h where a is the lower limit of integration, h > 0 is the length of the interval, tolerance is the tolerance, max_cols is the maximum number of columns to be used in the Richardson extrapolation to the limit, and err is set to 0 if iteration was successful and -1 if not.

C Source Code

  • The file, rombergs_method.c, contains a version of
    Rombergs_Integration_Method( ) written in C.