Lecture 3.5:
Recursive integration formulas from Romberg integration
More accurate integration formulas with smaller truncation error can be
obtained by interpolating several data points with higher-order interpolating
polynomials. For example, the fourth-order interpolating polynomial
P4(t) between five data points leads to the Boole's rule
of numerical integration. The Boole's rule has the global truncation error
of order O(h6).
However, the higher-order interpolating polynomials often do not provide
good approximations for integrals because they tend to oscillate wildly between
the samples (polynomial wiggle). As a result, they are seldom used past
Boole's rule. Another popular numerical algorithm is used instead to
reduce the truncation error of numerical integration. This is Romberg
integration based on the Richardson extrapolation algorithm (see
Lecture 3.3).
Example: from trapezoidal to Simpson's rule
Denote the trapezoidal rule as R1(h):
R1(h) = 0.5 h (I0 + 2I1 + 2I2 + ...
+ 2In-1 + In)
The trapezoidal rule has the global truncation error of order O(h2).
If the data set is equally spaced and the number of trapezoids n is even,
we can find two estimates for the same integral ST[I(t)]: one is
with the step size h (R1(h)) and the other one is
with the double step size 2h (R1(2h)). By canceling the
O(h2) truncation error, we define a new integration rule:
R2(h) = ( 4R1(h) - R1(2h) ) / 3
The new integration formula R2(h) is in fact the Simpson's rule:
R2(h) = h (I0 + 4I1 + 2I2 +
4I3 + 2I4 + ... + 2In-2 + 4In-1
+ In) / 3
The Simpson's rule has the global truncation error of order O(h4).
If the step size is sufficiently small, the Simpson's rule gives a much better
numerical estimate for the integral ST[I(t)] compared to
the trapezoidal rule.
Romberg integration algorithm
Construct the following mapping from Rk(h) to Rk+1(h):
Rk+1 = Rk(h) +
( Rk(h) - Rk(2h) ) / ( 4k - 1 )
There are two modifications of the Romberg integration algorithm: one is
to construct Rk(h) by doubling the step size h
and the other one is to construct Rk(h) by halving
the step size h. The first modification
is useful if a limited number of data values (tk,Ik)
is given as input. The second modification is more preferable if the function
I = I(t) is known analytically and the data values
(tk,Ik) can be computed for any given step size
h. We will use here the first variant of Romberg integration.
Given the data values (tk,Ik) for k = 0,1,...,n,
where n = 2m-1, the numerical approximations for the
integral ST[I(t)]
can be computed by using the trapezoidal rule for the step sizes: h,
2h, 4h, ..., and up to 2m-1 h. Then,
the recursive Romberg iterations can be performed and the numerical approximations
can be written in the following table
(similar to the table of Newton's divided differences):
| step size |
trapezoidal
rule |
Simpson's
rule |
Boole's
rule |
third
improvement |
| h |
R1(h) |
R2(h) |
R3(h) |
R4(h) |
| 2h |
R1(2h) |
R2(2h) |
R3(2h) |
|
| 4h |
R1(4h) |
R1(4h) |
|
|
| 8h |
R1(8h) |
|
|
|
The first row of this table gives a sequence of higher-order approximations
of the integral ST[I(t)].
The approximation Rk(h) has the truncation error of
O(h2k). As the order of numerical approximation increases, the
effects of round-off error begin to dominate and no further increase in accuracy
is obtained. The algorithm can be terminated when the relative error falls below
an error tolerance or when a maximum order m has been reached.
More accurate methods of numerical integration are based on Gauss quadrature methods
for orthogonal polynomials such as Legendre, Chebyshev, Laguerre and Hermite polynomials
(optional reading - chapter 7.5 of the main textbook).
- MATLAB codes for Romberg integration
- The algorithm outputs all
computed Romberg approximations up to the order m