Romberg's method
Encyclopedia
In numerical analysis
Numerical analysis
Numerical analysis is the study of algorithms that use numerical approximation for the problems of mathematical analysis ....

, Romberg's method is used to estimate the definite integral
Integral
Integration is an important concept in mathematics and, together with its inverse, differentiation, is one of the two main operations in calculus...




by applying Richardson extrapolation
Richardson extrapolation
In numerical analysis, Richardson extrapolation is a sequence acceleration method, used to improve the rate of convergence of a sequence. It is named after Lewis Fry Richardson, who introduced the technique in the early 20th century. In the words of Birkhoff and Rota, ".....

  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
Gaussian quadrature
In numerical analysis, a quadrature rule is an approximation of the definite integral of a function, usually stated as a weighted sum of function values at specified points within the domain of integration....

 and Clenshaw–Curtis quadrature are generally more accurate.

The method is named after Werner Romberg (1909–2003), who published the method in 1955.

Method

The method can be defined inductively




or


where




In big O notation
Big O notation
In mathematics, big O notation is used to describe the limiting behavior of a function when the argument tends towards a particular value or infinity, usually in terms of simpler functions. It is a member of a larger family of notations that is called Landau notation, Bachmann-Landau notation, or...

, the error for R(nm) is :


The zeroeth extrapolation, R(n, 0), is equivalent to the trapezoidal rule with 2n + 1 points; the first extrapolation, R(n, 1), is equivalent to Simpson's rule
Simpson's rule
In numerical analysis, Simpson's rule is a method for numerical integration, the numerical approximation of definite integrals. Specifically, it is the following approximation:...

 with 2n + 1 points. The second extrapolation, R(n, 2), is equivalent to Boole's rule
Boole's rule
In mathematics, Boole's rule, named after George Boole, is a method of numerical integration. It approximates an integralby using the values of ƒ at five equally spaced pointsIt is expressed thus Abramowitz and Stegun :...

 with 2n + 1 points. Further extrapolations differ from Newton Cote's 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 .

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's 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.

|-
|Number of pieces>
Trapezoid estimates First iteration second iteration >-
|
(4MA-LA)/3* (16MA-LA)/15 >-
|1
0 (4*480-0)/3 = 640 (16*880-640)/15 =896 >-
|2
480 (4*780-480)/3 = 880 (16*1006.67-880)/15 = 1015.11.. >-
|4
780 (4*950-780)/3 =1006.666.. >-
|8
950
  • 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
Error function
In mathematics, the error function is a special function of sigmoid shape which occurs in probability, statistics and partial differential equations...

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

Implementation

Here is an example of a computer implementation of the Romberg method (In the C programming language), it needs one vector and one variable; it also needs a sub-routine trapeze:


#include
#include
#include
#define MAX 6

int main
{
double s[MAX];
int i,k;
double var ;
for (i = 1; i< MAX; i++)
s[i] = 1;

for (k=1; k< MAX; k++)
{
for (i=1; i <=k; i++)
{
if (i1)
{
var = s[i];
s[i] = Trap(0, 1, pow(2, k-1)); // sub-routine trapeze
} // integrated from 0 and 1
/* pow is the number of subdivisions*/
else
{
s[k]= ( pow(4 , i-1)*s[i-1]-var )/(pow(4, i-1) - 1);

var = s[i];
s[i]= s[k];
}
}

for (i=1; i <=k; i++)
printf (" %f ", s[i]);

printf ("\n");
}

return 0;
}

External links
The source of this article is wikipedia, the free encyclopedia.  The text of this article is licensed under the GFDL.
 
x
OK