Richardson Extrapolation

Section C.1 Richardson Extrapolation

There are many approximation procedures in which one first picks a step size \(h\) and then generates an approximation \(A(h)\) to some desired quantity \(\cA\text{.}\) For example, \(\cA\) might be the value of some integral \(\int_a^b f(x)\,\dee{x}\text{.}\) For the trapezoidal rule with \(n\) steps, \(\De x=\tfrac{b-a}{n}\) plays the role of the step size. Often the order of the error generated by the procedure is known. This means

\begin{equation*} \cA=A(h)+Kh^k+K_1h^{k+1}+K_2h^{k+2}+\ \cdots \tag{E1} \end{equation*}

with \(k\) being some known constant, called the order of the error, and \(K,\ K_1,\ K_2,\ \cdots\) being some other (usually unknown) constants. If \(A(h)\) is the approximation to \(\cA=\int_a^b f(x)\,\dee{x}\) produced by the trapezoidal rule with \(\De x=h\text{,}\) then \(k=2\text{.}\) If Simpson’s rule is used, \(k=4\text{.}\)

Let’s first suppose that \(h\) is small enough that the terms \(K'h^{k+1}+K''h^{k+2}+\ \cdots\) in (E1) are small enough

Typically, we don’t have access to, and don’t care about, the exact error. We only care about its order of magnitude. So if \(h\) is small enough that \(K_1h^{k+1}+K_2h^{k+2}+\ \cdots\) is a factor of at least, for example, one hundred smaller than \(Kh^k\text{,}\) then dropping \(K_1h^{k+1}+K_2h^{k+2}+\ \cdots\) would not bother us at all.

that dropping them has essentially no impact. This would give

\begin{equation*} \cA=A(h)+Kh^k \tag{E2} \end{equation*}

Imagine that we know \(k\text{,}\) but that we do not know \(A\) or \(K\text{,}\) and think of (E2) as an equation that the unknowns \(\cA\) and \(K\) have to solve. It may look like we have one equation in the two unknowns \(K\text{,}\) \(\cA\text{,}\) but that is not the case. The reason is that (E2) is (essentially) true for all (sufficiently small) choices of \(h\text{.}\) If we pick some \(h\text{,}\) say \(h_1\text{,}\) and use the algorithm to determine \(A(h_1)\) then (E2), with \(h\) replaced by \(h_1\text{,}\) gives one equation in the two unknowns \(\cA\) and \(K\text{,}\) and if we then pick some different \(h\text{,}\) say \(h_2\text{,}\) and use the algorithm a second time to determine \(A(h_2)\) then (E2), with \(h\) replaced by \(h_2\text{,}\) gives a second equation in the two unknowns \(\cA\) and \(K\text{.}\) The two equations will then determine both \(\cA\) and \(K\text{.}\)

To be more concrete, suppose that we have picked some specific value of \(h\text{,}\) and have chosen \(h_1=h\) and \(h_2=\tfrac{h}{2}\text{,}\) and that we have evaluated \(A(h)\) and \(A(h/2)\text{.}\) Then the two equations are

\begin{align*} \cA&=A(h)+Kh^k \tag{E3a}\\ \cA&=A(h/2)+K\big(\tfrac{h}{2}\big)^k \tag{E3b} \end{align*}

It is now easy to solve for both \(K\) and \(\cA\text{.}\) To get \(K\text{,}\) just subtract (E3b) from (E3a).

\begin{align*} (\textrm{E3a})-(\textrm{E3b}):&\quad 0=A(h)-A(h/2) +\big(1-\tfrac{1}{2^k}\big)Kh^k\\ \implies&\quad K=\frac{A(h/2)-A(h)}{[1-2^{-k}]h^k} \tag{E4a} \end{align*}

To get \(\cA\) multiply (E3b) by \(2^k\) and then subtract (E3a).

\begin{align*} 2^k(\textrm{E3b})-(\textrm{E3a}):&\quad [2^k-1]\cA=2^kA(h/2)-A(h)\\ \implies&\quad \cA=\frac{2^kA(h/2)-A(h)}{2^k-1} \tag{E4b} \end{align*}

The generation of a “new improved” approximation for \(\cA\) from two \(A(h)\)’s with different values of \(h\) is called Richardson

Richardson extrapolation was introduced by the Englishman Lewis Fry Richardson (1881--1953) in 1911.

Extrapolation. Here is a summary

Equation C.1.1. Richardson extrapolation.

Let \(A(h)\) be a step size \(h\) approximation to \(\cA\text{.}\) If

\begin{equation*} \cA=A(h)+K h^k \end{equation*}

then

\begin{equation*} K=\frac{A(h/2)-A(h)}{[1-2^{-k}]h^k}\qquad \cA=\frac{2^kA(h/2)-A(h)}{2^k-1} \end{equation*}

This works very well since, by computing \(A(h)\) for two different \(h\)’s, we can remove the biggest error term in (E1), and so get a much more precise approximation to \(\cA\) for little additional work.

Example C.1.2. Richardson extrapolation with the trapezoidal rule.

Applying the trapezoidal rule (1.11.6) to the integral \(\cA=\int_0^1\frac{4}{1+x^2}\dee{x}\) with step sizes \(\frac{1}{8}\) and \(\frac{1}{16}\) (i.e. with \(n=8\) and \(n=16\)) gives, with \(h=\frac{1}{8}\text{,}\)

\begin{equation*} A(h)=3.1389884945\qquad A(h/2)=3.1409416120 \end{equation*}

So (E4b), with \(k=2\text{,}\) gives us the “new improved” approximation

\begin{equation*} \frac{2^2\times 3.1409416120 -3.1389884945}{2^2-1}=3.1415926512 \end{equation*}

We saw in Example 1.11.3 that \(\int_0^1\frac{4}{1+x^2}\dee{x}=\pi\text{,}\) so this new approximation really is “improved”:

\(A(1/8)\) agrees with \(\pi\) to two decimal places,
\(A(1/16)\) agrees with \(\pi\) to three decimal places and
the new approximation agrees with \(\pi\) to eight decimal places.

Beware that (E3b), namely \(\cA=A(h/2)+K\big(\tfrac{h}{2}\big)^k\text{,}\) is saying that \(K\big(\frac{h}{2}\big)^k\) is (approximately) the error in \(A(h/2)\text{,}\) not the error in \(\cA\text{.}\) You cannot get an “even more improved” approximation by using (E4a) to compute \(K\) and then adding \(K\big(\frac{h}{2}\big)^k\) to the “new improved” \(\cA\) of (E4b) — doing so just gives \(\cA+K\big(\tfrac{h}{2}\big)^k\text{,}\) not a more accurate \(\cA\text{.}\)

Example C.1.3. Example 1.11.16 revisited.

Suppose again that we wish to use Simpson’s rule (1.11.9) to evaluate \(\int_0^1 e^{-x^2}\,\dee{x}\) to within an accuracy of \(10^{-6}\text{,}\) but that we do not need the degree of certainty provided by Example 1.11.16. Observe that we need (approximately) that \(|K|h^4 \lt 10^{-6}\text{,}\) so if we can estimate \(K\) (using our Richardson trick) then we can estimate the required \(h\text{.}\) A commonly used strategy, based on this observation, is to

first apply Simpson’s rule twice with some relatively small number of steps and
then use (E4a), with \(k=4\text{,}\) to estimate \(K\) and
then use the condition \(|K| h^k\le 10^{-6}\) to determine, approximately, the number of steps required
and finally apply Simpson’s rule with the number of steps just determined.

Let’s implement this strategy. First we estimate \(K\) by applying Simpson’s rule with step sizes \(\tfrac{1}{4}\) and \(\tfrac{1}{8}\text{.}\) Writing \(\tfrac{1}{4}=h'\text{,}\) we get

\begin{equation*} A(h')=0.74685538 % 0.746855379790987 \qquad A(h'/2)=0.74682612 %0.746826120527467 \end{equation*}

so that (E4a), with \(k=4\) and \(h\) replaced by \(h'\text{,}\) yields

\begin{equation*} K=\frac{0.74682612 - 0.74685538}{[1-2^{-4}](1/4)^4} %=-\frac{0.000031211}{(1/4)^4} =-7.990\times 10^{-3} \end{equation*}

We want to use a step size \(h\) obeying

\begin{equation*} |K|h^4\le 10^{-6} \iff 7.990\times 10^{-3} h^4\le 10^{-6} \iff h \le\root{4}\of{\frac{1}{7990}} =\frac{1}{9.45} \end{equation*}

like, for example, \(h=\tfrac{1}{10}\text{.}\) Applying Simpson’s rule with \(h=\tfrac{1}{10}\) gives

\begin{equation*} A(1/10) = 0.74682495 \end{equation*}

The exact answer, to eight decimal places, is \(0.74682413\) so the error in \(A(1/10)\) is indeed just under \(10^{-6}\text{.}\)

Suppose now that we change our minds. We want an accuracy of \(10^{-12}\text{,}\) rather than \(10^{-6}\text{.}\) We have already estimated \(K\text{.}\) So now we want to use a step size \(h\) obeying

\begin{align*} |K|h^4\le 10^{-12} &\iff 7.99\times 10^{-3} h^4\le 10^{-12}\\ &\iff h \le\root{4}\of{\frac{1}{7.99\times 10^9}} =\frac{1}{299.0} \end{align*}

like, for example, \(h=\tfrac{1}{300}\text{.}\) Applying Simpson’s rule with \(h=\tfrac{1}{300}\) gives, to fourteen decimal places,

\begin{equation*} A(1/300) = 0.74682413281344 \end{equation*}

The exact answer, to fourteen decimal places, is \(0.74682413281243\) so the error in \(A(1/300)\) is indeed just over \(10^{-12}\text{.}\)

Prev Top Next