The question that verification can't dodge
The previous guide drew a hard line between verification and validation: verification asks "am I solving the equations right?", validation asks "am I solving the right equations?". This guide is about the sharpest tool we have for the first question. Suppose you have written a solver for a differential equation and your textbook proved the scheme is second-order accurate — its discretization error should shrink like O(h^2) as the grid spacing h goes to zero. A green test suite and a plot that "looks right" prove nothing about that. You need to watch the error actually fall, and fall at the promised rate.
Here is the trap that catches almost everyone the first time. A subtle bug — a stencil coefficient that should be 1/12 but is 1/10, a boundary node handled one cell off — often still converges. The error still shrinks as you refine, so a single run on a single grid looks fine. What the bug quietly destroys is the rate: your shiny second-order method silently degrades to first order, O(h^1). The answer is not wrong enough to scream, just wrong enough to waste a thousand CPU-hours buying accuracy you never get. A convergence study is the X-ray that exposes this.
Manufacturing a problem you already know the answer to
To measure error you need the truth to compare against, and that is exactly what real problems refuse to give you — if you knew the exact solution to your Navier-Stokes flow, you wouldn't be running a simulation. The method of manufactured solutions (MMS) cuts this knot with a beautifully backwards trick. Instead of picking an equation and hunting for its solution, you pick the solution first — any smooth function you like — and then work out which equation it solves.
Concretely: say your code claims to solve the Poisson equation u''(x) = f(x) on [0,1]. You manufacture a solution by decree — pick u_exact(x) = sin(3 x), say, just because it is smooth and non-trivial. Now plug it into the left side: u_exact''(x) = -9 sin(3 x). So if you set the source term to f(x) = -9 sin(3 x) and the boundary values to match sin(3x) at the ends, then sin(3 x) is the exact, provable solution of that modified problem. You feed f into your solver, get a numerical u_h back, and the error u_h - u_exact is now something you can compute to the last digit, on any grid.
- Choose the answer. Pick a smooth solution by decree, e.g. u_exact(x) = sin(3 x). It must not be accidentally special for your scheme.
- Derive the forcing. Apply the equation's operator to it. For u'' = f that means f(x) = u_exact''(x) = -9 sin(3 x).
- Set boundaries to match. Use u_exact evaluated at the ends as the boundary values, so the manufactured problem is fully consistent.
- Solve and measure. Feed f and the boundaries into your solver to get u_h, then compute the exact error e(h) = ||u_h - u_exact|| on any grid spacing h.
Reading the rate off a refinement table
With an exact error in hand, the order-of-accuracy study is mechanical. Run the same manufactured problem on a sequence of grids, each twice as fine as the last — h, h/2, h/4, h/8 — and record the error e(h) each time. If the method truly converges at order p, then e(h) is roughly C times h^p for some constant C. Halving h should multiply the error by (1/2)^p. For a second-order method that means the error drops by a factor of 4 each time you double the resolution; for first order, only a factor of 2.
You don't have to eyeball it. The observed order is recovered directly from any two consecutive grids: p_observed = log(e(h) / e(h/2)) / log(2). If your theory says 2 and the table prints 1.98, 2.01, 1.99 as the grid refines, you have verified your code — that convergence to the right rate is hard evidence the implementation matches the math. If instead it prints 1.0, you have a first-order bug masquerading as a second-order method, and you have caught it before it cost anyone a real result.
h e(h) ratio p_observed -------------------------------------------- 1/16 3.91e-03 - - 1/32 9.83e-04 3.98 1.99 1/64 2.46e-04 4.00 2.00 1/128 6.15e-05 4.00 2.00 p_observed = log( e(h) / e(h/2) ) / log(2)
Asymptotic, not magic: the honest fine print
The relation e(h) ~ C h^p is asymptotic — it only holds as h gets small enough to enter the asymptotic range, where the leading O(h^p) term dominates all the others. On coarse grids, higher-order terms (h^{p+1}, h^{p+2}, ...) still contribute, so the observed order on your first two grids may be nowhere near p. The fix is simply to refine further and watch the observed order settle toward a constant. If it never settles, or wanders, something is wrong: a bug, a non-smooth solution, or a singularity the theory never assumed.
There is also a floor you will hit if you refine too far, and it is the same floor every guide in this rung keeps returning to. Each grid you double costs more — refining a 3D mesh by 2x multiplies the work by 8 — and meanwhile the absolute error is shrinking toward the level of round-off in floating-point arithmetic. Eventually the discretization error you are trying to measure is buried under accumulated round-off, the table goes ragged, and the observed order collapses. The useful window is the band of grids between "too coarse for the asymptotic range" and "so fine that round-off pollutes the answer."
From one check to a standing guard
A convergence study is not a one-time graduation ritual; it is a test you keep. Once you have a manufactured problem with a known order, freeze it into your test suite as a regression test: every commit re-runs the refinement and asserts the observed order stays near p, say within a tolerance of 0.1. Now the day a colleague "optimizes" the stencil and quietly knocks it down to first order, the build goes red the same afternoon — not six months later when a paper's results won't reproduce. This is how the abstract idea of verification becomes a habit your software actually keeps.
Two honest caveats keep this from becoming overconfidence. First, MMS verifies the code, not the model: it proves you are solving your chosen equations correctly, but says nothing about whether those are the right equations for reality — that is validation against a real benchmark problem or experiment, a separate job. Second, MMS only certifies the regimes you actually exercise. A manufactured solution that is everywhere smooth will never test how your scheme handles a shock or a steep boundary layer; for those you manufacture solutions that contain them, or you accept that part stays unverified. Verification narrows where bugs can hide, but it never proves there are none.