› Richardson Extrapolation & Order of Accuracy Simulator Back
Numerical Analysis

Richardson Extrapolation & Order of Accuracy Simulator

From three step sizes and three numerical solutions, get the observed order of accuracy p and the Richardson-extrapolated true value u_ext in real time. Designed for CFD/FEM mesh convergence (V&V), ODE solver verification, and quick order checks on midpoint/trapezoidal/RK methods.

Parameters
Coarse step h₃
Largest step size (h₃ = max)
Medium step h₂
Intermediate step (usually h₃/2)
Fine step h₁
Smallest step size (h₁ = min)
Coarse solution u₃
Numerical solution obtained with h₃
Medium solution u₂
Numerical solution obtained with h₂
Fine solution u₁
Numerical solution obtained with h₁
Results
—
Step ratio r = h₂/h₁
—
Observed order p
—
Richardson value u_ext
—
Fine-step truncation error
—
Asymptotic ratio e₁/e₂
—
Order classification
—
Log-log plot — |error| vs step h (slope = p)

Blue points: errors of the three steps h₁,h₂,h₃. Yellow line: best-fit line through them (slope = observed order p). Green diamond: Richardson value (h → 0).

Solution vs step h (with extrapolated value)
|Error| vs step h (log-log axes)
Theory & Key Formulas

$$u_{exact} \approx u_1 + \frac{u_1-u_2}{r^p-1},\qquad p=\frac{\ln\bigl(|u_2-u_3|/|u_1-u_2|\bigr)}{\ln r}$$

r is the step ratio (h₂/h₁) and p is the observed order of accuracy. u₁,u₂,u₃ are the fine, medium and coarse solutions. When the three are in the asymptotic range, u_ext approximates the true value with 1-2 orders of magnitude smaller error.

$$e_h \approx C\,h^{p},\qquad \frac{e_1}{e_2}=\frac{1}{r^{p}}$$

Discrete error e_h falls like h^p. If the asymptotic ratio e₁/e₂ is close to 1/r^p, you are in the asymptotic range and the observed order is reliable.

$$\mathrm{GCI}=\frac{F_s\,|u_1-u_2|}{(r^p-1)\,|u_1|}$$

Grid Convergence Index (Roache 1998). Safety factor F_s=1.25 (three meshes) gives a reliable error band for the Richardson estimate. Standard practice in CFD V&V reporting.

Richardson Extrapolation and Order of Accuracy

🙋
"Richardson extrapolation" — first time I hear it. Why do I really need three step sizes? Isn't one enough?
🎓
Good question. A single solution u_h contains an unknown error e_h. To know "how far off" and "in which direction", you need at least two step sizes so you can take a difference. And if you also want to estimate the convergence order p as an unknown, you need three. With one point you cannot tell whether the true value is 0.81 or 0.82 or 0.85 — but with three you can back out "ah, this is converging to 0.811 with order p = 2".
🙋
So "observed order" and "theoretical order" are different things? The textbook says "midpoint rule is second-order accurate" — that's something else?
🎓
That is exactly the point. The "theoretical order" is derived on paper (midpoint = 2, RK4 = 4); the "observed order" is what you actually measure by running your code. They only agree when (1) the program has no bugs, (2) the step is fine enough to be in the asymptotic range, and (3) the solution is smooth. If a code that is supposed to be second-order gives an observed order of 0.8, it is almost certainly a boundary condition that has dropped to first order. Richardson extrapolation is the diagnostic that tells you whether your code really delivers the order it claims.
🙋
With the defaults I get p = 1.983 — labelled "2nd order". So is this midpoint rule or RK2?
🎓
Yes — this example is ∫₀¹ 1/(1+x²) dx = π/4 ≈ 0.7854 computed with the midpoint rule. The midpoint rule is second-order in theory (error ∝ h²), and the observed value of 1.983 lands right on 2. The Richardson value u_ext = 0.8111 is closer to the limit than u₁ = 0.8125 — although for this particular integrand symmetry makes the midpoint rule "super-converge" to a value slightly off the true integral. That is itself the moral: even with a clean order, you must verify that u_ext is physically the right limit, by comparing with an analytic solution, experiment, or another solver.
🙋
CFD courses always say "do a mesh convergence study" — in practice, how do I actually do it?
🎓
Three steps. (1) Run the same case on three meshes whose cell counts differ by about 2x each — coarse, medium, fine. (2) Record a quantity of interest (drag coefficient, pressure loss, peak temperature) as u₃, u₂, u₁. (3) Feed them into this tool to get p and the GCI. For a journal paper or conference, you can only claim "mesh-independent results" once p is within ±0.5 of the theoretical order and the GCI is below your error tolerance. Conversely, a single-mesh CFD report that just says "the picture looks nice" has zero V&V credibility. AIAA G-077 and ASME V&V20 are the standardised procedures.
🙋
The notes say "check the asymptotic ratio e₁/e₂ is close to 1/r^p". What if it is far off?
🎓
Then the three points are not yet in the asymptotic range, and the estimated order itself cannot be trusted. Two fixes. (a) Take a fourth, even finer mesh to form a 2:1 pair you trust more. (b) Check for discontinuities or strong singularities in the solution. With shocks or stress concentrations the smoothness of the solution caps the order, and p often plateaus at 1.0 to 1.5. In that case do not force it back to 2 — report "p = 1.2 observed" honestly, that is the correct attitude.

Frequently Asked Questions

From three step sizes h1<h2<h3 with a constant ratio r=h2/h1=h3/h2 and three corresponding numerical solutions u1,u2,u3, the observed order is p = ln(|u2−u3|/|u1−u2|) / ln(r). If it matches the theoretical order p_theory (1 for Euler, 2 for trapezoidal/RK2, 4 for RK4), the code and mesh are in the asymptotic range. A large deviation means the step is still too coarse, there is a bug in the code, or the order has dropped due to boundary conditions. The tool uses the same formula even when r is not exactly constant, as a quick estimate.
The Richardson estimate u_ext = u1 + (u1−u2)/(r^p − 1) removes the leading-order term of the asymptotic error expansion. When the three solutions are in the asymptotic range, the error of u_ext is typically 1-2 orders of magnitude smaller than the original error. The standard check is that the asymptotic error ratio e1/e2 is close to 1/r^p — for r=2 and p=2, the expected value is 0.25. If it is far off, the observed order is not trustworthy. Roache's (1998) Grid Convergence Index (GCI) adds a safety factor Fs=1.25 to convert this into a reliable error bound.
There are four common causes. (1) The step is still too coarse — the solutions are not in the asymptotic range yet. Refine further. (2) Boundary or initial conditions are implemented at low order (for example, a first-order staircase approximation) — this dominates the discrete error and pulls the observed order down to 0.5-1. (3) The solution has singularities or discontinuities (shocks, stress concentrations, re-entrant corners) — the smoothness of the solution caps the order, and the theoretical smooth-case order is not reached. (4) The iterative solver has not converged tightly enough, so iteration error exceeds discrete error — drive the residual at least one decade lower.
Yes — this is the canonical use case. CFD V&V (AIAA G-077, ASME V&V20) requires reporting the observed order of accuracy, the Richardson-extrapolated value, and the GCI from three meshes. The procedure is: (1) run the same case on coarse / medium / fine meshes (roughly 2x cells each), (2) record a quantity of interest (drag coefficient, peak stress, Nu number) as u3, u2, u1, (3) feed them here and verify p and u_ext and the GCI. If p is within ±0.5 of the theoretical order (typically 2 for FVM/FEM), the result is considered mesh-converged.

Real-World Applications

CFD verification & validation (V&V): Drag prediction for aircraft, efficiency of turbomachinery, automotive aerodynamics — commercial CFD reports (ANSYS Fluent, OpenFOAM, STAR-CCM+) almost always include an observed order and GCI. AIAA G-077 (aerospace), ASME V&V20 (fluid/thermal) and the NAFEMS guidelines all require "at least three meshes", "observed order check" and "GCI reporting". A "nice picture" from a single mesh is rejected at most conferences.

Finite element stress analysis: Peak stress at stress concentrations or J-integral values fluctuate strongly with mesh refinement. Run coarse / medium / fine meshes, take Richardson extrapolation to estimate "the true stress". With singularities such as re-entrant corners the observed order can drop below 1 and the extrapolation simply does not converge — in that case report "the finest mesh value and its uncertainty". Many NAFEMS benchmark problems include this kind of convergence study.

ODE solver verification: When you implement a new time integrator (Runge-Kutta, Adams-Bashforth, IMEX), halving Δt should drop the error by 2^p. A bug shows up as an "in-between" observed order (1.8 instead of 3.2, say) and is a strong debugging signal. This is a standard practice in textbook exercises and numerical-experiment sections of papers.

Climate and ocean models: General circulation models (GCM), ocean models (NEMO, MITgcm) also undergo spatial/temporal resolution convergence studies. Because of strong non-linearity and turbulence, the observed order is typically lower than the theoretical one and reaching the asymptotic range needs enormous compute. In practice people stop "when the quantity of interest stops moving in a steady direction" — a good example of where Richardson extrapolation hits its limits.

Common Misconceptions and Pitfalls

The most frequent mistake is to assume "if the three points are close together, the calculation has converged". For example with u₃=1.000, u₂=1.001, u₁=1.0011 the change is 0.001 then 0.0001, and you obtain an observed order of ln(10)/ln(2) ≈ 3.3. Calling that "third order!" is premature — it might just be three points stuck at the same wrong value (bias). Richardson extrapolation tells you "where you are heading and how fast", not "the value is correct". Always cross-check the limit u_ext against analytic solutions, experiments or an alternative solver.

Next, using the p formula on non-uniform step ratios. The p formula here implicitly assumes a constant ratio r = h₂/h₁ = h₃/h₂. With adaptive meshes or hybrid refinement the ratio can be 2 / 1.7 / 2.3 etc., and the observed-order estimate goes wildly off. The proper fix is Roache's (1998) iterative formula for non-uniform meshes — or, even better, deliberately producing a strict r = 2 mesh pair. This tool gives a meaningful number only when r is approximately constant.

Finally, "high observed order means I do not need fine meshes" — a tempting but wrong conclusion. Even with p = 4 (RK4, high-order FEM), the actual error e ≈ C·h^p can be large if the error constant C is large. What matters for engineering decisions is the actual current error, expressed by E_h = (u₁ − u₂)/(r^p − 1) in the GCI formula, not the order alone. The statement "high-order methods reach the same accuracy with coarser meshes" is true only after you have reached the asymptotic range; at the entry of that range high-order methods can be less accurate than first-order ones.

How to Use

  1. Enter three step sizes (h₁ coarse, h₂ medium, h₃ fine) in your numerical method—typically halving or reducing by factor 2–10 for finite difference, finite element, or ODE integration schemes.
  2. Input the corresponding three numerical solutions (u₁, u₂, u₃) obtained at each step size from your discretization (e.g., deflection in beam bending, pressure field in CFD, or integration error in Runge-Kutta).
  3. The simulator calculates step ratio r, observed order of accuracy p, Richardson extrapolated value u_ext, truncation error estimates, and asymptotic ratio e₁/e₂ to classify convergence behavior.

Worked Example

Cantilever beam deflection using finite elements: h₁=0.5m gives u₁=12.4mm; h₂=0.25m gives u₂=11.8mm; h₃=0.125m gives u₃=11.65mm (exact≈11.6mm). Step ratio r=2.0, observed order p≈1.9 (nearly second-order accurate FEM). Richardson extrapolation yields u_ext≈11.58mm, fine-step truncation error≈0.05mm, and asymptotic ratio e₁/e₂≈3.8 confirms p≈1.9 convergence pattern.

Practical Notes

  1. Ensure consistent step ratio r across all pairs (h₁/h₂ and h₂/h₃ ideally equal) to avoid spurious order estimates; non-uniform ratios degrade Richardson accuracy below 0.1p uncertainty.
  2. When p is negative or exceeds theoretical expectation (e.g., FEM should be p≈2), check for insufficient grid refinement, asymptotic region not reached, or nonsmooth solutions contaminating convergence.
  3. Richardson value u_ext provides an improved estimate only when p>1.5; use it to benchmark reference solutions for code validation in structural analysis, aerodynamics, or heat transfer problems.