Richardson Extrapolation

Richardson extrapolation is a technique to mathematically eliminate the leading error term in a numerical approximation. Applied to the central difference formula, each level of extrapolation removes the dominant error term and raises the accuracy by two orders. With just a few function evaluations, we can achieve an approximation far more accurate than what the central difference formula alone provides.

The Central Difference Error Structure

Before deriving the extrapolation, we need to know the exact form of the error in the central difference formula. Expanding $f(x+h)$ and $f(x-h)$ in Taylor series:

$f(x+h) = f(x) + f'(x)h + \frac{f''(x)}{2!}h^2 + \frac{f^{(3)}(x)}{3!}h^3 + \frac{f^{(4)}(x)}{4!}h^4 + \frac{f^{(5)}(x)}{5!}h^5 + O(h^6)$

$f(x-h) = f(x) - f'(x)h + \frac{f''(x)}{2!}h^2 - \frac{f^{(3)}(x)}{3!}h^3 + \frac{f^{(4)}(x)}{4!}h^4 - \frac{f^{(5)}(x)}{5!}h^5 + O(h^6)$

Subtracting the two equations, the even-power terms cancel:

$f(x+h) - f(x-h) = 2f'(x)h + \frac{2f^{(3)}(x)}{3!}h^3 + \frac{2f^{(5)}(x)}{5!}h^5 + O(h^7)$

Dividing by $2h$ :

$D_h = \frac{f(x+h)-f(x-h)}{2h} = f'(x) + \frac{f^{(3)}(x)}{3!}h^2 + \frac{f^{(5)}(x)}{5!}h^4 + O(h^6)$

This shows that $D_h$ equals the true derivative $f'(x)$ plus an infinite series of even-power error terms. The dominant error is $\frac{f^{(3)}(x)}{3!}h^2$ , which is $O(h^2)$ .

First-Level Extrapolation: $D_h^{(1)}$

Let’s derive the first level of Richardson extrapolation, denoted as $D_h^{(1)}$ . The main idea here is to combine two estimates of the derivative to completely cancel out the leading error term—which, for the central difference, is $O(h^2)$ .

We’ll start with the standard central difference formula at a step size $h$ :

$D_h = f'(x) + \frac{f^{(3)}(x)}{3!}h^2 + \frac{f^{(5)}(x)}{5!}h^4 + O(h^6)$

If we halve the step size to $h/2$ , the equation becomes:

$D_{h/2} = f'(x) + \frac{f^{(3)}(x)}{3!}\frac{h^2}{4} + \frac{f^{(5)}(x)}{5!}\frac{h^4}{16} + O(h^6)$

To knock out the $h^2$ term, we can multiply the $D_{h/2}$ equation by $2^2$ (or 4) to match the coefficients, and then subtract $D_h$ :

$4D_{h/2} - D_h = 3f'(x) - \frac{3}{4}\frac{f^{(5)}(x)}{5!}h^4 + O(h^6)$

Finally, dividing by 3 isolates $f'(x)$ and gives us our new estimate, free of the $h^2$ error:

$D_h^{(1)} = \frac{4D_{h/2} - D_h}{3}$

With the $h^2$ term gone, this first-level extrapolation formula is now accurate to $O(h^4)$ .

Second-Level Extrapolation: $D_h^{(2)}$

Now we can take it a step further. Let’s eliminate the new leading error term, $h^4$ , to find $D_h^{(2)}$ .

From the previous derivation, we know the exact error structure of our first-level estimate looks like this:

$D_h^{(1)} = f'(x) - \frac{1}{4}\frac{f^{(5)}(x)}{5!}h^4 + O(h^6)$

If we evaluate this same formula using a step size of $h/2$ , we get:

$D_{h/2}^{(1)} = f'(x) - \frac{1}{4}\frac{f^{(5)}(x)}{5!}\frac{h^4}{16} + O(h^6)$

This time around, we need to multiply $D_{h/2}^{(1)}$ by $2^4$ (or 16) so the $h^4$ coefficients match up perfectly:

$16D_{h/2}^{(1)} = 16f'(x) - \frac{1}{4}\frac{f^{(5)}(x)}{5!}h^4 + O(h^6)$

Subtracting $D_h^{(1)}$ from this equation eliminates the $h^4$ term entirely:

$16D_{h/2}^{(1)} - D_h^{(1)} = 15f'(x) + O(h^6)$

Divide by 15, and we arrive at our second-level extrapolation formula:

$D_h^{(2)} = \frac{16D_{h/2}^{(1)} - D_h^{(1)}}{15}$

By chaining these approximations together, we’ve bumped our accuracy up to $O(h^6)$ .

Generalized Formula

The same elimination pattern continues at every level. To eliminate the leading $h^{2n}$ error term, the general step is:

$D_h^{(n)} = \frac{2^{2n}\,D_{h/2}^{(n-1)} - D_h^{(n-1)}}{2^{2n} - 1}$

Level	Formula	Accuracy
Base	$D_h = \dfrac{f(x+h)-f(x-h)}{2h}$	$O(h^2)$
First	$D_h^{(1)} = \dfrac{4\,D_{h/2} - D_h}{3}$	$O(h^4)$
Second	$D_h^{(2)} = \dfrac{16\,D_{h/2}^{(1)} - D_h^{(1)}}{15}$	$O(h^6)$

Each additional level requires one more central difference evaluation at a halved step size, and raises the accuracy order by 2.

Worked Example

Let’s put the formulas to the test with a concrete problem. We have the function $f(x) = 4x^3 - 9e^{7x}$ , and we want to evaluate its derivative at $x=2.7$ using Richardson extrapolation, specifically finding $D_{0.2}^{(1)}$ and $D_{0.2}^{(2)}$ to 4 significant figures.

Part (a): First-Level Richardson Extrapolation

To find our first-level extrapolation, $D_{0.2}^{(1)}$ , we first need two basic central difference estimates: one at our starting step size $h=0.2$ , and another at half that step size, $h=0.1$ .

Let’s calculate $D_{0.2}$ first:

$D_{0.2} = \frac{f(2.7 + 0.2) - f(2.7 - 0.2)}{2(0.2)} = \frac{f(2.9) - f(2.5)}{0.4}$

Plugging in the numbers gives us $D_{0.2} = -1.384 \times 10^{10}$ (to 4 s.f.).

Now, we halve the step size to find $D_{0.1}$ :

$D_{0.1} = \frac{f(2.7 + 0.1) - f(2.7 - 0.1)}{2(0.1)} = \frac{f(2.8) - f(2.6)}{0.2}$

Which evaluates to $D_{0.1} = -1.103 \times 10^{10}$ (to 4 s.f.).

With both estimates in hand, we can blend them using our first-level extrapolation formula to knock out the $O(h^2)$ error:

$D_{0.2}^{(1)} = \frac{4D_{0.1} - D_{0.2}}{3} = \frac{4(-1.103 \times 10^{10}) - (-1.384 \times 10^{10})}{3}$

$= \frac{-4.412 \times 10^{10} + 1.384 \times 10^{10}}{3} = \frac{-3.028 \times 10^{10}}{3}$

This leaves us with a much sharper estimate:

$D_{0.2}^{(1)} = -1.009 \times 10^{10}$

Part (b): Second-Level Richardson Extrapolation

Moving on to the second level, our goal is to compute $D_{0.2}^{(2)}$ , which requires eliminating the $O(h^4)$ error. The formula for this is:

$D_{0.2}^{(2)} = \frac{16D_{0.1}^{(1)} - D_{0.2}^{(1)}}{15}$

We already figured out $D_{0.2}^{(1)}$ , but we’re missing $D_{0.1}^{(1)}$ . To get that, we have to drill down one step further and compute a raw central difference at $h=0.05$ .

$D_{0.05} = \frac{f(2.7 + 0.05) - f(2.7 - 0.05)}{2(0.05)} = \frac{f(2.75) - f(2.65)}{0.1}$

This comes out to $D_{0.05} = -1.038 \times 10^{10}$ (to 4 s.f.).

Next, we combine $D_{0.05}$ and $D_{0.1}$ to get our missing first-level piece:

$D_{0.1}^{(1)} = \frac{4D_{0.05} - D_{0.1}}{3} = \frac{4(-1.038 \times 10^{10}) - (-1.103 \times 10^{10})}{3}$

$= \frac{-4.152 \times 10^{10} + 1.103 \times 10^{10}}{3} = \frac{-3.049 \times 10^{10}}{3}$

So, $D_{0.1}^{(1)} = -1.016 \times 10^{10}$ .

Finally, we have everything we need. We can plug our two first-level estimates into the second-level formula:

$D_{0.2}^{(2)} = \frac{16(-1.016 \times 10^{10}) - (-1.009 \times 10^{10})}{15}$

$= \frac{-16.256 \times 10^{10} + 1.009 \times 10^{10}}{15} = \frac{-15.247 \times 10^{10}}{15}$

Giving us our highly accurate final answer:

$D_{0.2}^{(2)} = -1.016 \times 10^{10}$

Custom Extrapolation: Non-Standard Step Ratio

Richardson extrapolation isn’t limited to the standard $h$ and $h/2$ pairing. As long as we know the ratio between two step sizes, we can combine their estimates to knock out that leading $h^2$ error term.

Let’s see how this plays out if we decide to replace $h$ with $3h/4$ . We’ll start by recalling our standard central difference expansion:

$D_h = f'(x) + \frac{f^{(3)}(x)}{3!}h^2 + O(h^4)$

If we evaluate this using a step size of $3h/4$ instead, the $h^2$ gets scaled by $(3/4)^2$ , which gives us:

$D_{3h/4} = f'(x) + \frac{f^{(3)}(x)}{3!}\frac{9h^2}{16} + O(h^4)$

Our strategy remains exactly the same: we need to make the $h^2$ coefficients match so we can subtract them away. To turn that $9/16$ back into a 1, we multiply the entire $D_{3h/4}$ equation by its reciprocal, $16/9$ :

$\frac{16}{9}D_{3h/4} = \frac{16}{9}f'(x) + \frac{f^{(3)}(x)}{3!}h^2 + O(h^4)$

Now, if we subtract our original $D_h$ equation from this scaled version, the $h^2$ terms cancel out perfectly:

$\frac{16}{9}D_{3h/4} - D_h = \left(\frac{16}{9} - 1\right)f'(x) + O(h^4) = \frac{7}{9}f'(x) + O(h^4)$

To wrap things up, we isolate $f'(x)$ by dividing both sides by $7/9$ (which is the same as multiplying the numerator and denominator of our fraction by 9). This gives us our custom extrapolation formula:

$D_h^{(1)} = \frac{16D_{3h/4} - 9D_h}{7}$

And just like that, using an arbitrary step ratio, we’ve successfully derived an estimate that completely avoids the $O(h^2)$ error.