Romberg's method
From Wikipedia, the free encyclopedia
In numerical analysis, Romberg's method generates a triangular array consisting of numerical estimates of the definite integral
by using Richardson extrapolation 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.
The method can be defined inductively in this way:
or
where
In big O notation, the error for R(n,m) is:
[edit] Python implementation of Romberg's method
Here is an implementation of Romberg's method in the 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=[]
R.append([ 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] )
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 erf(1). 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.
