Manufactured solutions

Question: does the code solve these equations, at the order it claims?

You choose a smooth function u(t), and physverify adds the source term that makes u an exact solution of the reference equations. A correct code then reproduces u, and its error shrinks at the method's order. A code that solves slightly different equations, or takes a wrong step, does not.

Example: a half-step typo in the midpoint method

The midpoint method takes half a step to a midpoint, then a full step using the slope there. Writing h instead of h / 2 still gives a method that runs, stays stable and looks plausible on a plot.

"""A midpoint (RK2) solver, correct and with a one-character bug, checked by a manufactured solution."""
import numpy as np

from physverify import ManufacturedODE, mms_order_study

OMEGA = 2.0


def rhs(t, y):  # the code under test: a harmonic oscillator
    return np.array([y[1], -OMEGA**2 * y[0]])


def integrate(source, y0, t_end, n, buggy=False):
    """Midpoint method on y' = rhs(t, y) + source(t)."""
    h, y = t_end / n, np.array(y0, dtype=float)

    def f(t, y):
        return rhs(t, y) + source(t)

    for k in range(n):
        t = k * h
        k1 = f(t, y)
        mid = y + (h if buggy else h / 2) * k1  # bug: a full step instead of half a step
        y = y + h * f(t + h / 2, mid)
    return y


def reference_rhs(t, y):  # written again from the physics, not copied from the code
    return np.array([y[1], -4.0 * y[0]])


mms = ManufacturedODE(lambda t: np.array([np.sin(1.3 * t) + 0.2, 1.3 * np.cos(1.3 * t)]), reference_rhs)
ns = [10, 20, 40, 80, 160]
print("correct:", mms_order_study(integrate, mms, 1.0, ns, 2).message)
print("buggy:  ", mms_order_study(lambda s, y0, T, n: integrate(s, y0, T, n, buggy=True), mms, 1.0, ns, 2).message)

Output:

correct: observed order 2.00 (expected 2)
buggy:   observed order 0.98 is below the expected 2

The buggy version is still close to the right answer at these step sizes, so a comparison with one answer at a loose tolerance can pass it. The order study does not: the error halves when it should quarter.

What you provide

• ManufacturedODE(u, reference_rhs): the manufactured solution u(t) (returning an array) and the reference right-hand side, written from the physics. The source term is computed by complex-step differentiation, so u must accept complex t (use np.sin, not math.sin).
• integrate(source, y0, t_end, n): runs your code for exactly n steps from t = 0 with source(t) added to the right-hand side, and returns the final state.

Write the reference independently

If the reference right-hand side is copied from the code under test, a wrong sign in the code is copied too and the check passes. Write it again from the equations.

Variants

• mms_trajectory_study: the largest error over the whole run, not only at the end (catches errors that cancel by the final time).
• backward_mms_study: integrate from t_end back to 0. A code that evaluates its stages at t + h/2 while stepping backwards passes forwards but drops to first order here.
• ode_checks(...): bundles the standard checks for an ODE solver; see Getting started.