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.

<!-- Richardson Extrapolation demo — compute a derivative with a forward difference, then combine F(h) and F(h/2) to cancel the leading error term. -->
<div class="controls">
  <label>
    Function
    <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>
    Step count (rows) <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>Raw F(h) F(h)</th>
      <th>|error| raw</th>
      <th>Extrapolated R</th>
      <th>|error| extrap.</th>
    </tr>
  </thead>
  <tbody id="tbody"></tbody>
</table>
<p class="note" id="note"></p>

/* Table and control styles */
* { 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; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

/* Functions with exact derivatives for comparison */
const FUNCS = {
  sin: { f: x => Math.sin(x), df: x => Math.cos(x), label: 'sin(x)' },
  exp: { f: x => Math.exp(x), df: x => Math.exp(x), label: 'exp(x)' },
  ln:  { f: x => Math.log(x), df: x => 1 / x,       label: 'ln(x)' },
  x3:  { f: x => x * x * x,  df: x => 3 * x * x,    label: 'x³' }
};

const X0 = 1; // evaluate at x = 1

/* First-order forward difference — error is O(h) */
function fwdDiff(fn, x, h) {
  return (fn(x + h) - fn(x)) / h;
}

/* Richardson extrapolation — combines F(h) and F(h/2) to get O(h^2) accuracy */
function richardson(fn, x, h) {
  const Fh  = fwdDiff(fn, x, h);
  const Fh2 = fwdDiff(fn, x, h / 2);
  return 2 * Fh2 - Fh; // = 2*F(h/2) - F(h), cancels the O(h) error term
}

const STEPS = [0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625];
const stepIn  = document.getElementById('stepIn');
const stepVal = document.getElementById('stepVal');
const tbody   = document.getElementById('tbody');
const note    = document.getElementById('note');
const funcSel = document.getElementById('funcSel');

function fmt6(v) { return v.toExponential(2); }

function update() {
  const key   = funcSel.value;
  const fn    = FUNCS[key];
  const maxRow = parseInt(stepIn.value, 10); // number of rows = chosen step count
  const h     = STEPS[maxRow - 1];
  stepVal.textContent = h.toFixed(6).replace(/0+$/, '');

  const exact = fn.df(X0);
  tbody.innerHTML = '';

  /* fill one row per halving of h */
  for (let i = 0; i < maxRow; i++) {
    const hi   = STEPS[i];
    const raw  = fwdDiff(fn.f, X0, hi);
    const rich = richardson(fn.f, X0, hi);
    const eRaw  = Math.abs(raw  - exact);
    const eRich = Math.abs(rich - exact);

    const tr = document.createElement('tr');
    if (i === maxRow - 1) tr.className = 'highlight';
    tr.innerHTML =
      `<td>${hi}</td>` +
      `<td>${raw.toFixed(6)}</td>` +
      `<td class="err-raw">${fmt6(eRaw)}</td>` +
      `<td>${rich.toFixed(6)}</td>` +
      `<td class="err-rich">${fmt6(eRich)}</td>`;
    tbody.appendChild(tr);
  }

  // error ratio between last two extrapolated estimates
  if (maxRow >= 2) {
    const hi0  = STEPS[maxRow - 2];
    const hi1  = STEPS[maxRow - 1];
    const eR0  = Math.abs(richardson(fn.f, X0, hi0) - exact);
    const eR1  = Math.abs(richardson(fn.f, X0, hi1) - exact);
    const ratio = eR0 > 0 && eR1 > 0 ? (eR0 / eR1).toFixed(1) : '—';
    note.textContent = i18nFmt('Each time h halves, the extrapolated error shrinks by ~%ratio%x (ideal: ~4x for O(h²)).', { ratio });
  } else {
    note.textContent = '';
  }
}

stepIn.addEventListener('input', update);
funcSel.addEventListener('change', update);
update();

<!-- Richardson Extrapolation demo — compute a derivative with a forward difference, then combine F(h) and F(h/2) to cancel the leading error term. -->
<div class="controls">
  <label>
    Function
    <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>
    Step count (rows) <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>Raw F(h) F(h)</th>
      <th>|error| raw</th>
      <th>Extrapolated R</th>
      <th>|error| extrap.</th>
    </tr>
  </thead>
  <tbody id="tbody"></tbody>
</table>
<p class="note" id="note"></p>

/* Table and control styles */
* { 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; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

/* Functions with exact derivatives for comparison */
const FUNCS = {
  sin: { f: x => Math.sin(x), df: x => Math.cos(x), label: 'sin(x)' },
  exp: { f: x => Math.exp(x), df: x => Math.exp(x), label: 'exp(x)' },
  ln:  { f: x => Math.log(x), df: x => 1 / x,       label: 'ln(x)' },
  x3:  { f: x => x * x * x,  df: x => 3 * x * x,    label: 'x³' }
};

const X0 = 1; // evaluate at x = 1

/* First-order forward difference — error is O(h) */
function fwdDiff(fn, x, h) {
  return (fn(x + h) - fn(x)) / h;
}

/* Richardson extrapolation — combines F(h) and F(h/2) to get O(h^2) accuracy */
function richardson(fn, x, h) {
  const Fh  = fwdDiff(fn, x, h);
  const Fh2 = fwdDiff(fn, x, h / 2);
  return 2 * Fh2 - Fh; // = 2*F(h/2) - F(h), cancels the O(h) error term
}

const STEPS = [0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625];
const stepIn  = document.getElementById('stepIn');
const stepVal = document.getElementById('stepVal');
const tbody   = document.getElementById('tbody');
const note    = document.getElementById('note');
const funcSel = document.getElementById('funcSel');

function fmt6(v) { return v.toExponential(2); }

function update() {
  const key   = funcSel.value;
  const fn    = FUNCS[key];
  const maxRow = parseInt(stepIn.value, 10); // number of rows = chosen step count
  const h     = STEPS[maxRow - 1];
  stepVal.textContent = h.toFixed(6).replace(/0+$/, '');

  const exact = fn.df(X0);
  tbody.innerHTML = '';

  /* fill one row per halving of h */
  for (let i = 0; i < maxRow; i++) {
    const hi   = STEPS[i];
    const raw  = fwdDiff(fn.f, X0, hi);
    const rich = richardson(fn.f, X0, hi);
    const eRaw  = Math.abs(raw  - exact);
    const eRich = Math.abs(rich - exact);

    const tr = document.createElement('tr');
    if (i === maxRow - 1) tr.className = 'highlight';
    tr.innerHTML =
      `<td>${hi}</td>` +
      `<td>${raw.toFixed(6)}</td>` +
      `<td class="err-raw">${fmt6(eRaw)}</td>` +
      `<td>${rich.toFixed(6)}</td>` +
      `<td class="err-rich">${fmt6(eRich)}</td>`;
    tbody.appendChild(tr);
  }

  // error ratio between last two extrapolated estimates
  if (maxRow >= 2) {
    const hi0  = STEPS[maxRow - 2];
    const hi1  = STEPS[maxRow - 1];
    const eR0  = Math.abs(richardson(fn.f, X0, hi0) - exact);
    const eR1  = Math.abs(richardson(fn.f, X0, hi1) - exact);
    const ratio = eR0 > 0 && eR1 > 0 ? (eR0 / eR1).toFixed(1) : '—';
    note.textContent = i18nFmt('Each time h halves, the extrapolated error shrinks by ~%ratio%x (ideal: ~4x for O(h²)).', { ratio });
  } else {
    note.textContent = '';
  }
}

stepIn.addEventListener('input', update);
funcSel.addEventListener('change', update);
update();

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.