Romberg's method

From Wikipedia, the free encyclopedia

In numerical analysis, Romberg's method (Romberg 1955) generates a triangular array consisting of numerical estimates of the definite integral

by using Richardson extrapolation (Richardson 1910) repeatedly on the trapezium rule. Romberg's method evaluates the integrand at equally-spaced points. The integrand must have continuous derivatives, though fairly good results may be obtained if only a few derivatives exist. If it is possible to evaluate the integrand at unequally-spaced points, then other methods such as Gaussian quadrature and Clenshaw–Curtis quadrature are generally more accurate.

[edit] Method

The method can be defined inductively in this way:

or

where

In big O notation, the error for R(n,m) is:

The zeroeth extrapolation, R(n,0), is equivalent to the Trapezoidal Rule with n + 2 points; the first extrapolation, R(n,1), is equivalent to Simpson's rule with n + 2 points.

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).

[edit] Python implementation of Romberg's method

Here is an implementation of Romberg's method in Python.

def print_row(lst):
    print ' '.join('%11.8f' % x for x in lst)
 
def romberg(f, a, b, eps = 1E-8):
    """Approximate the definite integral of f from a to b by Romberg's method.
    eps is the desired accuracy."""
    R = [[0.5 * (b - a) * (f(a) + f(b))]]  # R[0][0]
    print_row(R[0])
    n = 1
    while True:
        h = float(b-a)/2**n
        R.append((n+1)*[None])  # Add an empty row.
        R[n][0] = 0.5*R[n-1][0] + h*sum(f(a+(2*k-1)*h) for k in range(1, 2**(n-1)+1)) # for proper limits
        for m in range(1, n+1):
            R[n][m] = R[n][m-1] + (R[n][m-1] - R[n-1][m-1]) / (4**m - 1)
        print_row(R[n])
        if abs(R[n][n-1] - R[n][n]) < eps:
            return R[n][n]
        n += 1
 
from math import *
 
# In this example, the error function erf(1) is evaluated.
print romberg(lambda t: 2/sqrt(pi)*exp(-t*t), 0, 1)

[edit] Example

As an example, the Gaussian function is integrated from 0 to 1, i.e. the Error function . The triangular array is calculated row by row and calculation is terminated if the two last entries in the last row differ less than 1E-8.

 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.

[edit] References

• Richardson, L. F. (1911), “The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, with an Application to the Stresses in a Masonry Dam”, Philosophical Transactions of the Royal Society of London. Series A 210: pp. 307–357, <http://links.jstor.org/sici?sici=0264-3952(1911)210%3C307%3ATAASBF%3E2.0.CO%3B2-J> 
• Romberg, W. (1955), “Vereinfachte numerische Integration”, Norske Videnskabers Selskab Forhandlinger (Trondheim) 28 (7): pp. 30–36 
• Thacher, Jr., Henry C. (1964), “Remark on Algorithm 60: Romberg integration”, Communications of the ACM 7 (7): 420-421, <http://portal.acm.org/citation.cfm?id=364520.364542> 
• Bauer, F.L.; Rutishauser & Stiefel, E. (1963), Metropolis, N. C., et al., ed., “New aspects in numerical quadrature”, Experimental Arithmetic, high-speed computing and mathematics, Proceedings of Symposia in Applied Mathematics (AMS) (no. 15): pp. 199–218 
• Bulirsch, Roland & Stoer, Josef (1967), “Handbook Series Numerical Integration. Numerical quadrature by extrapolation”, Numerische Mathematik 9: 271–278, <http://www-gdz.sub.uni-goettingen.de/cgi-bin/digbib.cgi?PPN362160546_0009> 
• Mysovskikh, I.P. (2002), “Romberg method”, in Hazewinkel, Michiel, Encyclopaedia of Mathematics, Springer-Verlag, ISBN 1-4020-0609-8, <http://eom.springer.de/r/r082570.htm> 

[edit] External links

• ROMBINT -- code for MATLAB (author: Martin Kacenak)
• Romberg's method is implemented in Maxima CAS
• Romberg Integration Source-Code (C++) with documentation
• ROMBERG -- c++ code for romberg integration
• Module for Romberg Integration
• Romberg's method -- plugin for Yacas