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