Richardson Extrapolation

Every numerical estimate carries a price tag: error. Approximate a derivative by a finite difference and you get the right answer plus a leftover term that shrinks as the step size $h$ goes to zero. The trouble is that shrinking $h$ costs computation, and at some point floating-point rounding wins anyway.

Richardson extrapolation, invented by Lewis Fry Richardson around 1910 and published in 1927, offers a shortcut. Instead of taking one estimate at a small $h$, take two estimates — one at $h$ and one at $h/2$ — and combine them in a ratio that makes the dominant error term cancel exactly. The result is an estimate that behaves as if you had used a much smaller step without actually computing at that step.

The trick is purely algebraic: if the error expands as $c_{1}h + c_{2}h^{2} + \cdots$, then the right linear combination of $F(h)$ and $F(h/2)$ eliminates the $c_{1}h$ term and leaves you with an error of order $h^{2}$ — one whole power better, for free.

Richardson's idea generalizes: apply it again to two second-order estimates and you reach fourth order. This iterated process is the foundation of Romberg integration, one of the most accurate practical quadrature methods available.

Choose a target function and a step size $h$. The demo computes the derivative numerically at $x = 1$ using a simple first-order forward difference at step $h$ and at $h/2$, then combines them via Richardson extrapolation. Compare the errors side by side.

<!-- {{c_intro}} -->
<div class="controls">
  <label>
    {{lbl_func}}
    <select id="funcSel">
      <option value="sin">sin(x)</option>
      <option value="exp">exp(x)</option>
      <option value="ln">ln(x)</option>
      <option value="x3">x³</option>
    </select>
  </label>
  <label>
    {{lbl_step}} <em>h</em>
    <input id="stepIn" type="range" min="1" max="6" value="3" step="1">
    <span id="stepVal"></span>
  </label>
</div>
<table id="tbl">
  <thead>
    <tr>
      <th>h</th>
      <th>{{th_raw}} F(h)</th>
      <th>{{th_err_raw}}</th>
      <th>{{th_rich}}</th>
      <th>{{th_err_rich}}</th>
    </tr>
  </thead>
  <tbody id="tbody"></tbody>
</table>
<p class="note" id="note"></p>

/* {{c_style}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; }
.controls { display: flex; flex-wrap: wrap; gap: .8rem 1.4rem; margin-bottom: 1rem; align-items: flex-end; }
label { display: flex; flex-direction: column; font-size: .85rem; color: #555; gap: .25rem; }
select, input[type=range] { font-size: .95rem; }
table { border-collapse: collapse; width: 100%; font-size: .88rem; }
th { background: #1d3557; color: #fff; padding: .4rem .6rem; text-align: center; }
td { padding: .35rem .6rem; text-align: right; border-bottom: 1px solid #dde3ea; }
td:first-child { text-align: center; font-family: ui-monospace, monospace; }
tr.highlight td { background: #e8f4e8; font-weight: 600; }
.err-raw  { color: #a33; }
.err-rich { color: #080; }
.note { font-size: .85rem; color: #555; margin-top: .6rem; line-height: 1.5; }

// Code not found

Notice that halving $h$ roughly halves the raw error — that is first-order convergence. But the extrapolated estimate cuts the error by roughly a factor of four each time $h$ halves, because the leading error term has been cancelled and only the $h^{2}$ term remains. That jump in the error column is Richardson extrapolation at work.

Start with the forward-difference formula for the derivative of $f$ at $x$:

$$F(h) = \frac{f(x+h) - f(x)}{h} = f'(x) + \frac{h}{2}f''(x) + \frac{h^{2}}{6}f'''(x) + \cdots$$

So $F(h)$ differs from the true $f'(x)$ by $c_{1}h + c_{2}h^{2} + \cdots$ where $c_{1} = \tfrac{1}{2}f''(x)$.

Replace $h$ with $h/2$:

$$F\!\left(\tfrac{h}{2}\right) = f'(x) + \frac{c_{1}h}{2} + \frac{c_{2}h^{2}}{4} + \cdots$$

Now form $2F(h/2) - F(h)$:

$$2F\!\left(\tfrac{h}{2}\right) - F(h) = f'(x) + \left(c_{2} - \frac{c_{2}}{2}\right)h^{2} + \cdots$$

The $c_{1}h$ terms cancel exactly, leaving an error of order $h^{2}$ — one order better than either raw estimate.

The general rule: if $F(h) = L + c_{p}h^{p} + c_{p+1}h^{p+1} + \cdots$ (error starts at order $p$), then

$$R = \frac{2^{p}\,F(h/2) - F(h)}{2^{p} - 1}$$

has error starting at order $p+1$. For $p=1$ the factor is $2^{1} = 2$, giving the formula above.

Romberg's method exploits this iteratively: build a triangular table of estimates, and each diagonal eliminates one more error order. Starting from a trapezoidal rule (which has error $O(h^{2})$, so $p=2$), Romberg's table converges extraordinarily fast — often reaching machine precision in fewer than a dozen function evaluations.

Richardson extrapolation is not just a textbook curiosity — it is embedded in tools scientists and engineers use every day:

• Romberg integration: build a table of trapezoidal-rule estimates and apply Richardson extrapolation at every level; the result converges faster than almost any other general-purpose quadrature rule.
• ODE step-size control: solvers like Runge-Kutta pairs (e.g., Dormand-Prince, the engine behind ode45 in MATLAB and solve_ivp in SciPy) compare a solution at step $h$ with one at $h/2$, use Richardson's formula to estimate the error, and adjust the step accordingly — automatically balancing speed against accuracy.
• Finite-difference derivatives: whenever a gradient or Jacobian is computed by finite differences, Richardson extrapolation can double the effective accuracy of each partial derivative with almost no extra cost.
• Lattice QCD and quantum field theory: numerical simulations on discrete spacetime grids must extrapolate to the continuum limit $a \to 0$ (where $a$ is the lattice spacing); Richardson-style extrapolation drives that analysis.
• Financial Greeks: options pricing libraries compute sensitivities (delta, gamma, vega) by finite differences; Richardson extrapolation makes those sensitivities more stable without shrinking the step into floating-point noise.

The underlying idea — combine two estimates to kill one error order — reappears wherever a computation has a known error expansion in a tunable parameter.

Richardson extrapolation is one of those ideas that feels almost too simple once you see it: the dominant error term in two estimates cancels when you subtract them in the right ratio. No smaller step, no more function evaluations — just algebra.

Yet that algebra unlocks genuine power. It is the engine behind Romberg integration, behind the step-size controllers in virtually every modern ODE solver, and behind the continuum extrapolations in lattice physics. Whenever a numerical result depends on a tunable parameter $h$ and the error is known to expand as a power series in $h$, Richardson's 1910 insight tells you how to do better at almost no cost.

The next time a simulation converges slowly, or a finite-difference gradient looks noisy, remember: two estimates chosen wisely are often worth far more than one estimate computed expensively.