Richardson Extrapolation

Plotting $T(h)$ as a function of $h^2$, Figure 3 shows that the behavior $T(h)=I+\alpha h^2$ is evident for even the modest values $N=2, 4, 8$. This plot suggests a very powerful method of finding much better estimates for the integral $I \equiv T(0)$: simply extrapolate the data in the plot to $h=0$. In general, we will not know a priori the analytic derivatives $f'(a)$ and $f'(b)$ and thus the slope $\alpha$ of the approach to the origin in the plot. But, with two values of $T(h)$ we can determine a very good estimate. Algebraically, we simply solve

$$\begin{eqnarray*} T(h) & \approx & T(0) + \alpha h^2\\ T(h/2) & \approx & T(0) + \alpha (h/2)^2, \end{eqnarray*}$$

to find the value

$$T(0) \approx \frac{4}{3} T(h/2) - \frac{1}{3} T(h).$$ (14)

The fourth column of Table 1 shows that the resulting extrapolation formula (16) indeed converges much more quickly than $T(h)$, as $h^4$ in fact.

Figure 3: Linear behavior of the trapezoid rule when viewed as a function of $h^2$.

The extrapolation (16) translates directly into an uniform-interval method. From the definition of the trapezoid rule (13),

$$-\frac{1}{3} T(h) = -\frac{1}{3} h \left( \frac{1}{2} f_0 + f_1 +f_2 + \ldots + f_{N-1} + \frac{1}{2} f_{N} \right)$$

$$\frac{4}{3} T(h/2) = \frac{4}{3} \frac{h}{2} \left( \frac{1}{2} f_0 + f_{1/2} +f_1 + \ldots + f_{(N-1)/2} + \frac{1}{2} f_{N} \right),$$

where $f_{i-1/2} \equiv f(x_{i-1/2})$ with $x_{i-1/2} \equiv a + (i-1/2) h$ denoting sample values at the intermediate points which $T(h/2)$ uses but $T(h)$ does not. Summing the above two formulae then gives an explicit interval method,

$$S(h/2) \equiv \frac{4}{3} T(h/2) - \frac{1}{3} T(h)$$

$$= \frac{h}{6} \left( f_0 + 4 f_{1/2} + 2 f_1 + \ldots + 4 f_{(N-1)/2} + f_{N} \right),$$ (15)

which is none other than the familiar Simpson's rule! For this rule, we denote the fact that the spacing between sample points $x_i$ is now $h/2$ and not $h$ by writing $S(h/2)$. The final column of Table 1 confirms numerically that Simpson's rule evaluated in the traditional way is precisely equivalent to the extrapolation (16).

Understanding this relationship allows us to rapidly determine the error in Simpson's from the error in the trapezoid rule (16):

$$S(h/2)-I = \frac{4}{3} T(h/2) - \frac{1}{3} T(h) -I$$

$$= \frac{4}{3} \left(T(h/2)-I\right) - \frac{1}{3} \left(T(h) -I\right)$$

$$= \frac{4}{3} \left( \frac{1}{12} \left(\frac{h}{2}\right)^2 \left(... ...t(\frac{h}{2}\right)^4 \left( f'''(b)-f'''(a) \right)+ {\mathcal O}(h^6)\right)$$ (16)

$$- \frac{1}{3} \left( \frac{1}{12} h^2 \left( f'(b)-f'(a) \right) - \frac{1}{720} h^4 \left( f'''(b)-f'''(a) \right)+ {\mathcal O}(h^6)\right)$$

$$= \frac{1}{2880} h^4 \left( f'''(b)-f'''(a) \right)+ {\mathcal O}(h^6),$$

where the quadratic terms have cancelled by the construction which led to (16), leaving only the quartic term. The patterns of shifting digits in Table 1 again confirm our computation of the prefactor,

$$\begin{eqnarray*} \frac{1}{2880} \left(f'''(b)-f'''(a)\right) & = & \frac{1}{28... ...c{1}{16}\right) = \frac{90}{46080 \ln 2}\\ & = & 0.002818\ldots \end{eqnarray*}$$

One may proceed in this manner to methods of arbitrarily high order. For instance, extrapolating $T(h)$ to $h^2=0$ as a quadratic function using the information which $T(h)$, $T(h/2)$ and $T(h/4)$ provide results in the sixth-order approach known as Bode's rule.

The idea of extracting an exact analytic result by extrapolating several calculations at finite values of a parameter to the result at zero value is an extremely powerful tool in scientific computing, where it is known as Richardson extrapolation. (The same idea occurs in statistical mechanics where it is often known as finite-size scaling.)