Reading Time: 8 minutes

Key Takeaways

  • Truncation error is the foundation of all code verification — understanding Taylor series derivation explains why MMS works and why convergence rates match theoretical predictions.
  • The Method of Manufactured Solutions (MMS) is now automated — frameworks like MOOSE and FEniCSx use SymPy to derive source terms symbolically, eliminating the hand-calculation errors that plagued early verification efforts.
  • You only need 3–4 mesh refinements — once observed order-of-accuracy matches theoretical prediction, further refinement is wasted effort. The real question is: when is verification “good enough”?
  • Leading-order error terms are diagnostic — the structure of the truncation error expansion tells you exactly which derivative approximation is failing, guiding your scheme selection.

What You’re Actually Measuring

Before diving into equations, let’s clarify what we’re actually verifying. In the simulation community, verification and validation are often confused — but they answer fundamentally different questions:

  • Verification asks: “Are we solving the equations correctly?” It’s about internal consistency. Does the numerical implementation faithfully solve the mathematical model, regardless of whether that model describes reality?
  • Validation asks: “Are we solving the right equations?” It compares simulation results against physical experiments or benchmark data.

If Post 182’s V&V framework overview introduced you to the big picture, this article drills into the mathematical machinery that makes verification rigorous. Think of it as the companion piece to our convergence studies guide, which covers the empirical side. Here, we’ll explain why convergence studies work and show you how to automate them.

Mental framework: the truncation-error compass. Imagine truncation error as your navigation tool. The Taylor series expansion doesn’t just give you a bound — it reveals the direction your scheme is biased. First-order asymmetry in forward differences? That’s your compass needle pointing to the dominant error derivative.

Truncation Error: The Mathematical Foundation

The most common source of numerical error in finite difference, finite element, and finite volume methods is truncation error — the difference between the exact differential operator and its discrete approximation.

Here’s the thing most textbooks don’t emphasize enough: truncation error derivation via Taylor series isn’t an academic exercise. It’s the practical foundation for every grid refinement study you’ll ever run. Let’s see why.

Consider the backward Euler time derivative:

$$\frac{u(t) – u(t – \Delta t)}{\Delta t} = \frac{du}{dt} + \frac{1}{2} \frac{d2u}{dt2} \Delta t + O(\Delta t^2)$$

The leading error term is positive and proportional to Δt. For forward Euler, the sign flips. For a central difference in space:

$$\frac{u(x + \Delta x) – u(x – \Delta x)}{2\Delta x} = \frac{du}{dx} + O(\Delta \Delta x^2)$$

Notice the symmetry — the odd-order terms cancel out, giving second-order accuracy. The even/odd symmetry of the Taylor expansion is why central differences outperform forward/backward schemes.

A Practical Comparison

Scheme Formula Truncation Error Order
Forward difference (u(x+Δx) - u(x)) / Δx −½ u'' Δx 1st
Backward difference (u(x) - u(x-Δx)) / Δx −½ u'' Δx 1st
Central difference (u(x+Δx) - u(x-Δx)) / 2Δx +(1/6) u''' Δx² 2nd

The table above isn’t just comparison — it’s a diagnostic tool. If your convergence study shows first-order when second-order was expected, the leading-order error term tells you exactly what’s wrong.

The Truncation Error Expression

For a general discretization, the truncation error takes the form:

$$\tau = C h^r$$

where h is the discretization parameter (mesh size or time step) and r is the convergence rate. This expression is not just an asymptotic bound — it’s the actual quantity you verify empirically. Every grid refinement study is fundamentally an attempt to confirm that the theoretical truncation error expression matches the measured error.

This connection between the theoretical τ and the empirical error measurement is what makes code verification rigorous rather than hand-wavy.

From Truncation Error to Convergence

Now that we have the truncation error expression, let’s connect it to the big theorem of numerical analysis: the Lax-Richtmyer equivalence theorem.

The theorem states, in plain English:

For linear, well-posed problems, a consistent, stable discretization is convergent.

Let’s unpack each term at a practitioner level:

Consistent means the truncation error goes to zero as h → 0. If your scheme has τ = O(h²), it’s consistent.

Stable means errors don’t grow unboundedly. The numerical solution stays bounded relative to the initial data. You verify stability through von Neumann analysis, energy methods, or practical CFL condition checks.

Convergent means the numerical solution approaches the exact solution as h → 0.

The theorem’s proof sketch (which you should know intuitively, not memorize):

  1. Consistency ensures the discrete operator approximates the continuous operator.
  2. Stability bounds the propagation of errors through each time step.
  3. Together, they guarantee that the total error — truncation error accumulated over N = T/Δt steps — remains bounded and converges to zero.

The Lax-Richtmyer theorem explains why stability requirements exist. You can’t just reduce h arbitrarily; you must also ensure your scheme is stable. For explicit time integration, that’s the CFL condition. For implicit schemes, you’re usually safe, but the iteration tolerance becomes the new stability concern.

The Method of Manufactured Solutions

Now for the practical payoff: how do you verify that your code solves the equations correctly? The most rigorous approach is the Method of Manufactured Solutions (MMS).

Salari and Park’s 2000 Sandia report established MMS as the industry standard, and it now has over 555 citations. The reason? MMS works for any PDE — linear, nonlinear, coupled, multiphysics — whereas analytical benchmarks only exist for simple test cases.

The MMS Workflow

The beauty of MMS is that it flips the verification problem on its head:

  1. Choose a manufactured solution u_manufactured(x, t) — an arbitrary smooth function
  2. Substitute into the PDE to derive the forcing/source term that makes u_manufactured an exact solution
  3. Derive initial and boundary conditions from u_manufactured
  4. Run the simulation with these modified inputs
  5. Compare the computed solution against u_manufactured

If the numerical solution matches the manufactured solution within the expected error bounds, your code is verified.

Let’s walk through this with concrete Python code. We’ll use two of the most widely adopted open-source frameworks: MOOSE and FEniCSx.

MOOSE MMS: Symbolic Automation

MOOSE’s mms module wraps SymPy to derive forcing functions automatically. Here’s how you set up a spatial convergence study for a diffusion equation:

import mms

# Define the PDE and manufactured solution
fs, ss = mms.evaluate("-div(grad(u))", "sin(2*pi*x)*sin(2*pi*y)")

# Print forcing function for MOOSE input file
mms.print_fparser(fs)

# Print exact solution and forcing function as MOOSE hit syntax
mms.print_hit(fs, "force")
mms.print_hit(ss, "exact")

The output tells you exactly what to put in your .i input file:

8*pi^2*sin(2*x*pi)*sin(2*y*pi)
[force]
  type = ParsedFunction
  expression = '8*pi^2*sin(2*x*pi)*sin(2*y*pi)'
[]
[exact]
  type = ParsedFunction
  expression = 'sin(2*x*pi)*sin(2*pi*y)'
[]

This symbolic derivation is crucial. For a simple 1D diffusion equation, you could derive the source term by hand. For Navier-Stokes or elasticity with coupled terms, symbolic computation isn’t optional — it’s the only way to avoid errors.

The MOOSE input file then looks like this:

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 8
  ny = 8
[]

[Kernels]
  [diff]
    type = ADDiffusion
    variable = u
  []
  [force]
    type = BodyForce
    variable = u
    function = force
  []
[]

[BCs]
  [all]
    type = FunctionDirichletBC
    variable = u
    function = exact
    boundary = 'left right top bottom'
[]
[]

[Postprocessors]
  [error]
    type = ElementL2Error
    function = exact
    variable = u
[]
[]

Then you automate the convergence study:

import mms

# Run 4 levels of refinement for both 1st and 2nd order elements
df1 = mms.run_spatial("diffusion_mms.i", 4, console=False)
df2 = mms.run_spatial("diffusion_mms.i", 4, "Variables/u/order=SECOND")

fig = mms.ConvergencePlot(xlabel="Element Size ($h$)", ylabel="$L_2$ Error")
fig.plot(df1, label="1st Order")
fig.plot(df2, label="2nd Order")
fig.save("convergence_plot.png")

On a log-log plot, the slope of each line gives you the observed convergence rate. For first-order elements, the slope should approach 2. For second-order elements, it should approach 3.

FEniCSx: Pure Python Convergence Studies

FEniCSx (the successor to FEniCS/Dolfin) provides an equally powerful Python interface. Here’s how you compute error norms and convergence rates:

from dolfinx import default_scalar_type
from dolfinx.fem import (
    Expression, Function, functionspace,
    assemble_scalar, dirichletbc, form,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dx, grad, inner
from mpi4py import MPI

import ufl
import numpy as np

def u_ex(mod):
    return lambda x: mod.cos(2 * mod.pi * x[0]) * mod.cos(2 * mod.pi * x[1])

u_numpy = u_ex(np)
u_ufl = u_ex(ufl)

def solve_poisson(N=10, degree=1):
    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    x = SpatialCoordinate(mesh)
    f = -div(grad(u_ufl(x)))
    V = functionspace(mesh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v)) * dx
    L = f * v * dx
    u_bc = Function(V)
    u_bc.interpolate(u_numpy)
    facets = locate_entities_boundary(
        mesh, mesh.topology_dim - 1, lambda x: np.full(x.shape[1], True)
    )
    dofs = locate_dofs_topological(V, mesh.topology_dim - 1, facets)
    bcs = [dirichletbc(u_bc, dofs)]
    problem = LinearProblem(
        a, L, bcs=bcs,
        petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )
    return problem.solve(), u_ufl(x)

The key insight here is reliable error norm computation. When the error is small, direct computation (u_ex - uh)^2 can suffer from round-off errors because you’re subtracting two nearly equal numbers. The FEniCSx tutorial recommends interpolating both solutions into a higher-order function space first:

def error_L2(uh, u_ex, degree_raise=3):
    degree = uh.function_space.ufl_element().degree
    family = uh.function_space.ufl_element().family_name
    mesh = uh.function_space.mesh
    
    # Create higher-order space for accurate subtraction
    W = functionspace(mesh, (family, degree + degree_raise))
    u_W = Function(W)
    u_W.interpolate(uh)
    
    u_ex_W = Function(W)
    u_ex_W.interpolate(u_ex)
    
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
    
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_global = mesh.comm.allreduce(assemble_scalar(error), op=MPI.SUM)
    return np.sqrt(error_global)

Then you run the convergence study:

Ns = [4, 8, 16, 32, 64]
Es = np.zeros(len(Ns))
hs = np.zeros(len(Ns))

for i, N in enumerate(Ns):
    uh, u_ex = solve_poisson(N, degree=1)
    Es[i] = error_L2(uh, u_numpy)
    hs[i] = 1.0 / Ns[i]
    print(f"h: {hs[i]:.2e} Error: {Es[i]:.2e}")

# Compute observed convergence rates
rates = np.log(Es[1:] / Es[:-1]) / np.log(hs[1:] / hs[:-1])
print(f"Rates: {rates}")

The output for first-order elements shows rates approaching 2:

Rates: [1.61 1.89 1.97 1.99]

For second-order elements, the rates approach 3. This is the empirical manifestation of the theoretical truncation error analysis we discussed earlier.

Practical Convergence Studies: Interpreting the Results

You’ve run your grid refinement study. You have your log-log plot. The slopes are close to theoretical values. But what does this actually tell you?

The Observed Order Formula

For mesh sizes h_i and h_{i-1} with corresponding errors E_i and E_{i-1}, the observed order is:

$$p \approx \frac{\ln(E_{i-1} / E_i)}{\ln(h_{i-1} / h_i)}$$

If the grid refinement ratio is r ≈ 2 (common in studies), this simplifies to:

$$p \approx \log_2(E_{i-1} / E_i)$$

When your observed order approaches the theoretical prediction, you’ve confirmed that your code operates in the asymptotic regime. This is the verification milestone.

Degree-Raising: A Practical Trick

For L2 error norms, direct computation can hide the true convergence rate due to round-off. The degree-raising technique (interpolating into a space one degree higher before subtraction) is the standard workaround in FEniCSx. Without it, you might report second-order convergence when the code is actually converging at third-order.

When to Stop Verifying

Here’s the practical reality: verification is expensive. Each mesh refinement doubles (or quadruples for 2D) your computational cost. You can’t refine infinitely. So when do you stop?

The field has converged on three practical heuristics:

1. Observed order matches theoretical prediction. When you have 3+ mesh refinements and the observed order converges to the theoretical value within tolerance (typically ±0.1), your code is verified in the asymptotic regime. No further refinement is needed for verification purposes.

2. Grid Convergence Index (GCI) falls below ~1%. The GCI quantifies the uncertainty band from mesh refinement. When GCI < 1% of your quantity of interest, the numerical error is negligible relative to modeling uncertainty.

3. Numerical error falls below physical uncertainty. If your material properties have 5% uncertainty, refining until numerical error is 0.001% is wasteful. The error estimate τ = C h^r tells you where the practical floor is.

A Practical Decision Checklist

  • [ ] Run 3–4 mesh refinements with refinement ratio r ≈ 2
  • [ ] Compute observed order p from log-log slope
  • [ ] Verify p matches theoretical prediction within ±0.1
  • [ ] Compute GCI for the finest mesh
  • [ ] Compare GCI to modeling uncertainty (typically 1–5%)
  • [ ] Stop when: p is confirmed AND GCI < 1% of QoI

If you’re running MMS for a production code, this checklist is sufficient. For safety-critical applications (like nuclear reactor simulation per the MPACT manual), additional checks for each term in the governing equation are required.

Summary

The journey from truncation error to code verification follows a clear path:

  1. Taylor series reveals the theoretical truncation error expression τ = C h^r
  2. Lax-Richtmyer theorem guarantees convergence when the scheme is consistent and stable
  3. MMS provides the practical workflow to verify that your code achieves the theoretical convergence rate
  4. Grid refinement studies confirm that observed rates match theoretical predictions

The key insight that ties this together: truncation error isn’t just a theoretical concept — it’s the measurable quantity your grid refinement studies confirm. When the observed order matches theoretical prediction across 3+ refinements, your code is verified.

For next steps, check our V&V framework overview for the big-picture context and our convergence studies guide for practical mesh quality recommendations that complement these verification workflows.

References

  • Salari, K. & Park, K. (2000). “Code Verification by the Method of Manufactured Solutions.” Sandia National Laboratories Report SAND2000-0949. Source
  • Oberkampf, W. & Roy, C. (2010). “Verification and Validation in Scientific Computing.” Cambridge University Press.
  • Langtangen, H.P. “Truncation Error Analysis for Finite Difference Methods.” Numerical Methods for PDEs documentation.
  • MOOSE MMS documentation. MOOSE framework
  • FEniCSx convergence tutorial. Dolfinx tutorial
  • Kindo, T. “Verify Simulations with the Method of Manufactured Solutions.” COMSOL Blog. Source