In numerical analysis, Romberg's method is used to estimate the definite integral <math display="block"> \int_a^b f(x) \, dx </math> by applying Richardson extrapolation repeatedly on the trapezium rule or the rectangle rule (midpoint rule). The estimates generate a triangular array. Romberg's method is a Newton–Cotes formula – it 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 is named after Werner Romberg, who published the method in 1955.
Method
Using <math display="inline">h_n = \frac{(b-a)}{2^{n+1</math>, the method can be inductively defined by
<math display="block">\begin{align}
R(0,0) &= h_0 (f(a) + f(b)) \\
R(n,0) &= \tfrac{1}{2} R(n{-}1,\,0) + 2h_n \sum_{k=1}^{2^{n-1 f(a + (2k-1)h_{n-1}) \\
R(n,m) &= R(n,\,m{-}1) + \tfrac{1}{4^m-1} (R(n,\,m{-}1) - R(n{-}1,\,m{-}1)) \\
&= \frac{1}{4^m-1} ( 4^m R(n,\,m{-}1) - R(n{-}1,\, m{-}1))
\end{align}</math>
where <math> n \ge m </math> and In big O notation, the error for is: <math> O{\left(h_n^{2m+2}\right)}.</math>
The zeroth extrapolation, , is equivalent to the trapezoidal rule with points; the first extrapolation, , is equivalent to Simpson's rule with points. The second extrapolation, , is equivalent to Boole's rule with points. The 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.
By labelling our <math display="inline">O(h^2)</math> approximations as <math display="inline">A_0\big(\frac{h}{2^n}\big)</math> instead of <math display="inline">R(n,0)</math>, we can perform Richardson extrapolation with the error formula defined below:
<math display="block"> \int_a^b f(x) \, dx = A_0\bigg(\frac{h}{2^n}\bigg)+a_0\bigg(\frac{h}{2^n}\bigg)^{2} + a_1\bigg(\frac{h}{2^n}\bigg)^{4} + a_2\bigg(\frac{h}{2^n}\bigg)^{6} + \cdots </math>
Once we have obtained our <math display="inline">O(h^{2(m+1)})</math> approximations <math display="inline">A_m\big(\frac{h}{2^n}\big)</math>, we can label them as <math display="inline">R(n,m)</math>.
When function evaluations are expensive, it may be preferable to replace the polynomial interpolation of Richardson with the rational interpolation proposed by .
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.
alt=One-piece approximation|thumb|One-piece. Note since it starts and ends at zero, this approximation yields zero area.
alt=Two-piece approximation|thumb|Two-piece
alt=Four-piece approximation|thumb|Four-piece
alt=Eight-piece approximation|thumb|Eight-piece
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 . 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
- The third iteration uses the next power of 4: on the values derived by the second iteration.
- The pattern is continued until there is one estimate.
{| class="wikitable col1right col2right"
|-
! Number of pieces || Trapezoid estimates || First iteration || Second iteration || Third iteration
|-
| || || || ||
|-
|1||0 || || ||
|-
|2||16 || || ||
|-
|4||30 || || ||
|-
|8||39 || || ||
|}
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 10<sup>−8</sup>.
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):
<syntaxhighlight lang="c" line="1">
- include <stdio.h>
- include <math.h>
void print_row(size_t i, double *R) {
printf("R[%2zu] = ", i);
for (size_t j = 0; j <= i; ++j) {
printf("%f ", R[j]);
}
printf("\n");
}
/*
INPUT:
(*f) : pointer to the function to be integrated
a : lower limit
b : upper limit
max_steps: maximum steps of the procedure
acc : desired accuracy
OUTPUT:
Rp[max_steps-1]: approximate value of the integral of the function f for x in [a,b] with accuracy 'acc' and steps 'max_steps'.
- /
double romberg(double (*f)(double), double a, double b, size_t max_steps, double acc)
{
double R1[max_steps], R2[max_steps]; // buffers
double *Rp = &R1[0], *Rc = &R2[0]; // Rp is previous row, Rc is current row
double h = b-a; //step size
Rp[0] = (f(a) + f(b))*h*0.5; // first trapezoidal step
print_row(0, Rp);
for (size_t i = 1; i < max_steps; ++i) {
h /= 2.;
double c = 0;
size_t ep = 1 << (i-1); //2^(n-1)
for (size_t j = 1; j <= ep; ++j) {
c += f(a + (2*j-1) * h);
}
Rc[0] = h*c + .5*Rp[0]; // R(i,0)
for (size_t j = 1; j <= i; ++j) {
double n_k = pow(4, j);
Rc[j] = (n_k*Rc[j-1] - Rp[j-1]) / (n_k-1); // compute R(i,j)
}
// Print ith row of R, R[i,i] is the best estimate so far
print_row(i, Rc);
if (i > 1 && fabs(Rp[i-1]-Rc[i]) < acc) {
return Rc[i];
}
// swap Rn and Rc as we only need the last row
double *rt = Rp;
Rp = Rc;
Rc = rt;
}
return Rp[max_steps-1]; // return our best guess
}
</syntaxhighlight>
Here is an implementation of the Romberg method (in the Python programming language):
<syntaxhighlight lang='python' line='1'>
from numpy import arange
def romberg(f,a,b,maxorder):
"""
Approximates an integral using Romberg integration.
Args:
f: The function to integrate.
a: Lower limit of integration.
b: Upper limit of integration.
maxorder: Maximum recursion depth
Returns:
The approximate value of the integral.
"""
assert a < b and isinstance(maxorder,int)
minstepsize = float( (1./2)**maxorder * (b-a) )
- algorithm is simplified and accelerated by caching
f=memorize(f)
@memorize
def rule(h, order):
if ( order == 1 ):
partial_sum = (f(a) + f(b))/2
partial_sum += sum(map(f, arange(a+h, b-h/2, h)))
return h*partial_sum
- Richardson extrapolation formula
W = 4.**(order-1)
return (W*rule(h, order-1) - rule(2*h, order-1)) / (W - 1)
return rule(minstepsize, maxorder)
def memorize(f):
_cache_ = dict()
def cachedf(*argv):
if not ( argv in _cache_):
_cache_[argv] = f(*argv)
return _cache_[argv]
cachedf._cache_ = _cache_
return cachedf
def main():
from numpy import sin, exp, pi
- Example: integrate sin(x)*exp(x/6)+x over [0,2pi]
a,b,f = 0, 2*pi, lambda x : sin(x)*exp(x/6) + x
exact = 36*(1-exp(pi/3))/37 + 2*pi**2
print("Order Value Rel_Err\n"+'-'*30)
for maxorder in range(1,9):
value = romberg(f, a, b, maxorder)
err = abs(value - exact)/exact
print('%2d'%maxorder, '%.15f'%value, '%1.1e'%err)
main()
</syntaxhighlight>
References
Citations
Bibliography
External links
- ROMBINT – code for MATLAB (author: Martin Kacenak)
- Free online integration tool using Romberg, Fox–Romberg, Gauss–Legendre and other numerical methods
- SciPy implementation of Romberg's method
- Romberg.jl — Julia implementation (supporting arbitrary factorizations, not just <math>2^n+1</math> points)