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 P_{4}(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(h^{6}). 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).
Denote the trapezoidal rule as R_{1}(h):
R_{1}(h) = 0.5 h (I_{0} + 2I_{1} + 2I_{2} + ... + 2I_{n-1} + I_{n})
The trapezoidal rule has the global truncation error of order O(h^{2}). If the data set is equally spaced and the number of trapezoids n is even, we can find two estimates for the same integral S_{T}[I(t)]: one is with the step size h (R_{1}(h)) and the other one is with the double step size 2h (R_{1}(2h)). By canceling the O(h^{2}) truncation error, we define a new integration rule:
R_{2}(h) = ( 4R_{1}(h) - R_{1}(2h) ) / 3
The new integration formula R_{2}(h) is in fact the Simpson's rule:
R_{2}(h) = h (I_{0} + 4I_{1} + 2I_{2} + 4I_{3} + 2I_{4} + ... + 2I_{n-2} + 4I_{n-1} + I_{n}) / 3
The Simpson's rule has the global truncation error of order O(h^{4}). If the step size is sufficiently small, the Simpson's rule gives a much better numerical estimate for the integral S_{T}[I(t)] compared to the trapezoidal rule.
Construct the following mapping from R_{k}(h) to R_{k+1}(h):
R_{k+1} = R_{k}(h) + ( R_{k}(h) - R_{k}(2h) ) / ( 4^{k} - 1 )
There are two modifications of the Romberg integration algorithm: one is to construct R_{k}(h) by doubling the step size h and the other one is to construct R_{k}(h) by halving the step size h. The first modification is useful if a limited number of data values (t_{k},I_{k}) is given as input. The second modification is more preferable if the function I = I(t) is known analytically and the data values (t_{k},I_{k}) can be computed for any given step size h. We will use here the first variant of Romberg integration. Given the data values (t_{k},I_{k}) for k = 0,1,...,n, where n = 2^{m-1}, the numerical approximations for the integral S_{T}[I(t)] can be computed by using the trapezoidal rule for the step sizes: h, 2h, 4h, ..., and up to 2^{m-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 | R_{1}(h) | R_{2}(h) | R_{3}(h) | R_{4}(h) |
2h | R_{1}(2h) | R_{2}(2h) | R_{3}(2h) | |
4h | R_{1}(4h) | R_{1}(4h) | ||
8h | R_{1}(8h) |
The first row of this table gives a sequence of higher-order approximations of the integral S_{T}[I(t)]. The approximation R_{k}(h) has the truncation error of O(h^{2k}). 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).