Romberg's method

Method

Using

h n = 1 2 n ( b − a )

the method can be inductively defined by

R ( 0 , 0 ) = h 1 ( f ( a ) + f ( b ) ) R ( n , 0 ) = 1 2 R ( n − 1 , 0 ) + h n ∑ k = 1 2 n − 1 f ( a + ( 2 k − 1 ) h n ) R ( n , m ) = R ( n , m − 1 ) + 1 4 m − 1 ( R ( n , m − 1 ) − R ( n − 1 , m − 1 ) )

or

R ( n , m ) = 1 4 m − 1 ( 4 m R ( n , m − 1 ) − R ( n − 1 , m − 1 ) )

where n ≥ m and m ≥ 1 . In big O notation, the error for R(n, m) is (Mysovskikh 2002):

O ( h n 2 m + 2 ) .

The zeroeth extrapolation, R(n, 0), is equivalent to the trapezoidal rule with 2ⁿ + 1 points; the first extrapolation, R(n, 1), is equivalent to Simpson's rule with 2ⁿ + 1 points. The second extrapolation, R(n, 2), is equivalent to Boole's rule with 2ⁿ + 1 points. Further extrapolations differ from Newton Cotes formulas. In particular further Romberg extrapolations expand on Boole's rule in very slight ways, modifying weights into ratios similar as in Boole's rule. In contrast, further Newton Cotes methods produce increasingly differing weights, eventually leading to large positive and negative weights. This is indicative of how large degree interpolating polynomial Newton Cotes methods fail to converge for many integrals, while Romberg integration is more stable.

When function evaluations are expensive, it may be preferable to replace the polynomial interpolation of Richardson with the rational interpolation proposed by Bulirsch & Stoer (1967).

A geometric example

To estimate the area under a curve the trapezoid rule is applied first to one-piece, then two, then four, and so on.

After trapezoid rule estimates are obtained Richardson extrapolation is applied

For the first iteration the two piece and one piece estimates are used in the formula (4 X (more accurate) - (less accurate))/3 The same formula is then used to compare the four piece and the two piece estimate, and likewise for the higher estimates

For the second iteration the values of the first iteration are used in the formula (16(more accurate)-less accurate)/15

The third iteration uses the next power of 4: (64 (More accurate) - less accurate)/63 on the values derived by the second iteration.

The pattern is continued until there is one estimate.

MA stands for more accurate, LA stands for less accurate

Example

As an example, the Gaussian function is integrated from 0 to 1, i.e. the error function erf(1) ≈ 0.842700792949715. The triangular array is calculated row by row and calculation is terminated if the two last entries in the last row differ less than 10⁻⁸.

0.77174333 0.82526296 0.84310283 0.83836778 0.84273605 0.84271160 0.84161922 0.84270304 0.84270083 0.84270066 0.84243051 0.84270093 0.84270079 0.84270079 0.84270079

The result in the lower right corner of the triangular array is accurate to the digits shown. It is remarkable that this result is derived from the less accurate approximations obtained by the trapezium rule in the first column of the triangular array.

Implementation

Here is an example of a computer implementation of the Romberg method (in the C programming language).