JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

Richardson Extrapolation

If you know exactly how an approximation's error shrinks with the step size, you can combine two crude estimates to cancel the leading error and leap to a far more accurate one — almost for free. Richardson extrapolation is that trick, and it turns the trapezoidal rule into the formidable Romberg method.

The leftover-error idea

In the previous guide we met the central difference D(h) = (f(x+h) - f(x-h)) / (2h) as an approximation to f'(x), and we saw its truncation error: a Taylor expansion shows D(h) = f'(x) + c2*h^2 + c4*h^4 + ..., where c2 = f'''(x)/6 and the higher terms ride on f''''', f''''''', and so on. The headline is that the error is O(h^2): halving h cuts it by four. But look closer at the structure. The error is not some shapeless blur — it is a known polynomial in h, dominated by a single leading term c2*h^2 whose coefficient we may not know but whose shape we do.

That single fact is the whole opening. If the leading error behaves exactly like h^2, then I can deliberately compute two estimates, one at step h and one at step h/2, and play them against each other. Each is wrong, but they are wrong in a predictable, proportional way — and two wrongs whose ratio you know can be combined to make a right. This is the heart of Richardson extrapolation: do not fight the leading error, measure it and subtract it. The same idea applies to any method whose order of accuracy is a known power of h, which is most of numerical analysis.

Two wrongs make a right: the cancellation

Write A for the true value we want, A(h) for the estimate at step h, and suppose A(h) = A + c*h^p + (higher-order terms), where the order p is known — p = 2 for the central difference. Halving the step gives A(h/2) = A + c*(h/2)^p + ... = A + c*h^p / 2^p + ... . Now the leading error of A(h) is exactly 2^p times the leading error of A(h/2). So form the combination that kills that term: take 2^p copies of the finer estimate, subtract one copy of the coarser, and divide by (2^p - 1). The unknown c cancels cleanly, A survives, and what remains of the error is the next term — now O(h^{p+2}) for an even-power series like the central difference.

Improved = ( 2^p * A(h/2) - A(h) ) / ( 2^p - 1 )

# central difference: error is O(h^2), so p = 2, 2^p = 4
Improved = ( 4 * D(h/2) - D(h) ) / 3      # error now O(h^4)
One extrapolation step. The weights 2^p and -1 are chosen so the c*h^p terms exactly cancel; dividing by (2^p - 1) keeps the true value A unchanged. For the central difference this turns two O(h^2) estimates into one O(h^4) estimate.

Make it concrete. Suppose D(h) = 1.040 and D(h/2) = 1.012, both crude guesses at a derivative whose true value is 1.000. Richardson gives (4*1.012 - 1.040)/3 = (4.048 - 1.040)/3 = 3.008/3 = 1.0027 — already much closer than either input, and we spent only the function evaluations we needed for the two estimates anyway. The two errors, +0.040 and +0.012, stood in roughly the 4-to-1 ratio the theory predicts, and the weighted difference squeezed most of that error out.

Climbing the table, one order at a time

Why stop at one step? The improved estimate has its own leading error, now of order p+2, with its own unknown coefficient — the same situation as before, one rung higher. So extrapolate again, this time with the weight chosen for the new order. Build a triangular table: the first column holds raw estimates A(h), A(h/2), A(h/4), ...; each later column is formed by Richardson-combining two neighbours from the column to its left, using the weight 2^q for whatever order q that column kills. Each new column wipes out the next error term, so accuracy can climb O(h^2), O(h^4), O(h^6), ... astonishingly fast — for smooth functions, a handful of halvings reaches near machine precision.

Romberg: Richardson married to the trapezoidal rule

The most famous payoff lives in integration. The composite trapezoidal rule is the humble workhorse — approximate the area under a curve by straight-line trapezoids — and on its own it is only O(h^2) accurate, which sounds mediocre. But for a smooth integrand a deep result (the Euler-Maclaurin formula) says its error is a clean even-power series in h: error = c2*h^2 + c4*h^4 + c6*h^6 + ... . That is precisely the structure Richardson devours. So compute trapezoidal estimates at h, h/2, h/4, ..., then run them up the extrapolation table. The result is Romberg integration, and the very first extrapolation of the trapezoidal rule turns out to be exactly Simpson's rule in disguise — a satisfying hint that we have rediscovered something the next guide treats head-on.

  1. Compute the composite trapezoidal estimate R(0,0) using one interval, then R(1,0), R(2,0), ... by repeatedly halving the step. Reuse the previous evaluations — each halving only adds the new midpoints, so refining is cheap.
  2. Fill in each extrapolated column with R(i,j) = ( 4^j * R(i, j-1) - R(i-1, j-1) ) / ( 4^j - 1 ). Here 4^j = 2^p plays the weight that kills the order-2j error term.
  3. The diagonal entries R(i,i) are your best estimates; watch successive diagonals stop changing in the digits you care about and stop there.

The practical magic is that each trapezoidal refinement reuses every function value from the coarser grid, so the dominant cost is just the new points, while the extrapolation column-building is a few cheap arithmetic combinations. You pay roughly for the finest trapezoidal grid and get accuracy that no trapezoidal grid of any feasible size could reach. For a smooth integrand on a fixed interval, Romberg often nails ten or more correct digits from a table only five or six rows tall.

What it cannot rescue you from

Richardson extrapolation attacks truncation error — the error of the formula itself — and only that. It does nothing about round-off, and for differentiation that distinction is sharp. Recall the step-size trade-off from the previous guide: as h shrinks, truncation error falls but subtractive cancellation in (f(x+h) - f(x-h)) makes round-off error grow, so there is an optimal h below which a smaller step makes things worse. Extrapolating a finite-difference formula with very small h does not escape that floor; it can amplify it, because the weights 4, -1 magnify the noise already polluting the inputs. Use moderate, not tiny, steps when extrapolating a derivative.

There is also an honest limit to how far the table climbs. The extrapolation assumes infinitely many smooth error terms exist, but in floating-point every column is computed from inputs already accurate only to about 16 digits, so once a column reaches that accuracy, the next columns merely shuffle round-off noise — the table stalls, and a later diagonal can even drift away from the true value. The right reading of a Romberg or differentiation table is to watch the diagonal: trust the point where successive entries first agree to the precision you need, and ignore the churn beyond it. The method is a beautiful accelerator, not a perpetual-motion machine.