Verifying scientific code

How a piece of simulation code is shown to be right or wrong: independent tests from limits, conservation, and convergence, six deliberately faulty steppers that the tests must catch, and the written contract that says what a verifier checks, in what order, and with what evidence.

Course

Classical mechanics, lesson 15 of 28

Slot

AI2, number 3 of the 9 slots of the path Physics reasoning and scientific agents

Methods

Simulation

Difficulty

Advanced, level 3 of 3

What you’ll learn

Test a simulation against limits, conservation, and convergence with references independent of the code
Catch a deliberately faulty stepper with a named test and explain why other tests miss it
Evaluate a piece of code against a stated claim, with the evidence in a fixed form

Prerequisites

Numerical integrators, Energy conservation.

Before you begin

Python: reading a short function and a unit test

Related

none

Previous

The Kepler orbit

Next

Automated analysis pipelines

With scripts, this page keeps a record in your browser of whether you completed the lesson. Nothing is stored without them.

With scripts, this place suggests what to do next, with the reason, from the progress saved in your browser. Without them, follow the learning path.

A simulation that runs and prints a smooth curve has shown nothing. The earlier lessons of this course built their integrators carefully, and each was compared with an exact solution, but a reader who is given a stepper written by someone else, or by a program, has no such history. The question is how to find out whether it is right, and what “right” is allowed to mean.

The answer is a set of tests that do not trust the code they test. Each compares the result with something the code did not compute: the motion with no force, the exact oscillator, a quantity the physics conserves. This lesson runs five such tests on the steppers of the course, and then on six steppers that are wrong on purpose, each with one fault of the kind that is made in practice: a missing factor, a sign, a clock that does not run, a state that leaks from one step to the next. Every fault is caught by a test with a name, and the tests that miss a fault are as instructive as the one that finds it.

The tests, their order, and the evidence a verifier reports are written down as the Scientific validation contract, in the file docs/validation-contract.md of the repository. The lesson follows the contract and names it where it applies; the multi-agent review of a later lesson and the rubric for code written by an agent use the same contract.

Assumptions

The tests and the program rest on these assumptions. Each one names what goes wrong when it does not hold.

A claim exists. The tests check what a method says about itself: its order of accuracy and whether it keeps the energy bounded. Code with no claim has nothing to check against, and a claim written after seeing the results is a description, not a test.
The oscillator is a fair reference. It has an exact solution, so every error is the stepper’s. It is also a linear system with a force that does not depend on the velocity, so a fault that only matters for a nonlinear or velocity-dependent force can pass every test here; the contract says what to add for a new model.
The tolerances are fixed first. They are constants of the module, chosen from what a correct method does, before any faulty one is run. A tolerance widened until a test passes makes the test measure nothing, as the first exercise shows.
Faults are single. Each faulty stepper has one fault. Real faults come in groups and hide each other; the tests do not become less useful for it, but the pattern of failures becomes harder to read.
Round-off is below the tolerances. The exact tests allow $10^{-12}$ , several orders above the round-off of the numbers used. A method with large intermediate values would need its own tolerance.

What a test may compare with

A test of a simulation has three parts: a claim, a reference, and a tolerance. The reference decides whether the test is worth anything. If the reference is computed by the code under test, or by code written from the same mistaken derivation, the two agree whatever the mistake. A reference is independent when it comes from somewhere the code did not: an exact solution, a limiting case, a law, or the experiment’s own bookkeeping.

The contract orders five tests, from the cheapest and most basic to the ones that need a reference solution and a long run. They are the names the rest of the page uses.

clock. A step of length $\Delta t$ advances the time by $\Delta t$ . The reference is the step the experiment asked for.
purity. A step depends only on its arguments. The same state, taken again after other steps, gives the same result.
limits. Cases with a closed form. With no force the motion is uniform, $x = x_0 + v_0 t$ ; under a constant force $a$ the velocity is $v_0 + a t$ exactly, for every method of the course, because the velocity update of each is exact for a constant acceleration, and the position is $x_0 + v_0 t + \tfrac12 a t^2$ , exactly for a method of order 2 or more. A method of order 1 may be off by $\tfrac12 a\,\Delta t^2$ in each step, and only the excess over that allowance counts.
convergence-order. Halving the step divides the error by $2^p$ for a method of order $p$ . The reference is the exact oscillator of the earlier lessons, and the number is the slope of the error against the step on logarithmic axes.
bounded-energy. For a method that claims a bounded energy, the energy of the oscillator stays within a band that does not grow. A single number cannot tell a band from a trend, so the test compares the width of the band in the first periods of a long run with its width in the last.

The first three need no reference solution at all, which is why they come first: they cost little, and a method that fails them has a fault that the later results cannot be trusted to show. The last two are the ones that find a fault which changes the numbers only slightly.

No test proves a method right. Each asks one question that a right method answers one way. The contract says what a report must therefore contain: for each test the claim, the measured value, the tolerance, and the verdict, and a statement of what was not covered. A faulty method is caught when some test fails; a method that passes has passed those tests, on this model, at these steps, and no more.

The program

The estimates are in pbc.verification.estimates, the tests and the verdict in pbc.verification.checks, and the faulty steppers in pbc.verification.faulty. The listings are read from the files when the site is built. The order of convergence between two steps, and the fitted order over all of them, are the two estimates the convergence test uses:

def observed_orders(steps: ArrayLike, errors: ArrayLike) -> NDArray[np.float64]:
    """The order of convergence seen between each pair of successive steps.

    ``steps`` are decreasing step sizes and ``errors`` the error of the run
    at each. With the error proportional to ``h ** p``, the pair gives
    ``p = log(e_i / e_(i+1)) / log(h_i / h_(i+1))``; for halved steps that is
    the base-2 logarithm of the ratio of the errors. An error that is zero
    or not finite gives ``nan``: nothing was measured.
    """
    h, e = _convergence_data(steps, errors)
    with np.errstate(divide="ignore", invalid="ignore"):
        orders = np.log(e[:-1] / e[1:]) / np.log(h[:-1] / h[1:])
    orders[
        ~np.isfinite(e[:-1]) | ~np.isfinite(e[1:]) | (e[:-1] <= 0) | (e[1:] <= 0)
    ] = np.nan
    return orders

src/pbc/verification/estimates.py, lines 28 to 43, read from the repository when this page was built.

def fitted_order(steps: ArrayLike, errors: ArrayLike) -> float:
    """The slope of the straight line through the points (log h, log error),
    fitted by least squares: one number for the order of a method.

    ``nan`` when any error is zero or not finite.
    """
    h, e = _convergence_data(steps, errors)
    if not np.all(np.isfinite(e)) or np.any(e <= 0):
        return float("nan")
    return float(np.polyfit(np.log(h), np.log(e), 1)[0])

src/pbc/verification/estimates.py, lines 46 to 55, read from the repository when this page was built.

The claim of a method and the evidence of a test are two small records:

@dataclass(frozen=True)
class Claim:
    """What a method says about itself, and so what is checked.

    ``order`` is the order of accuracy; ``bounded_energy`` says whether the
    energy of an oscillator stays within a band for all time.
    """

    order: int
    bounded_energy: bool

src/pbc/verification/checks.py, lines 41 to 50, read from the repository when this page was built.

@dataclass(frozen=True)
class Finding:
    """The evidence of one check: the claim it tests, what was measured,
    against what tolerance, the verdict, and how to run the check again."""

    check: str
    claim: str
    reproduction: str
    passed: bool
    measured: float
    tolerance: float
    detail: str

src/pbc/verification/checks.py, lines 53 to 64, read from the repository when this page was built.

The convergence test measures the largest distance from the exact oscillator over a run, at the times of the experiment and not of the stepper, for a series of steps, and compares the fitted order with the claim:

def check_convergence_order(build: Callable[[], Stepper], claim: Claim) -> Finding:
    """Convergence: halving the step divides the error by 2 ** order."""
    counts = [50, 100, 200, 400, 800]
    errors = np.array([oscillator_error(build, n) for n in counts])
    order = fitted_order(PERIOD / np.array(counts), errors)
    return _finding(
        "convergence-order",
        "halving the step divides the error by 2 ** the claimed order",
        bool(np.isfinite(order)) and abs(order - claim.order) <= ORDER_MARGIN,
        order,
        ORDER_MARGIN,
        f"fitted order {order:.2f}, claimed {claim.order}",
    )

src/pbc/verification/checks.py, lines 189 to 201, read from the repository when this page was built.

The names of the five tests, in the order of the contract:

def list_checks():
    print(", ".join(CHECKS))

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 40 to 41, read from the repository when this page was built.

clock, purity, limits, convergence-order, bounded-energy

Worked examples

The correct steppers pass

Three methods of the course, with what each claims: velocity Verlet is of order 2 and keeps the energy bounded, symplectic Euler is of order 1 and keeps it bounded, and Runge-Kutta 4 is of order 4 and loses energy steadily, so it claims no bound. The cell runs every test that applies to each claim and prints the verdict of each as pass or FAIL, with - where the claim does not include the test.

def verdicts(build, claim):
    findings = {f.check: f for f in verify(build, claim)}
    return findings, {
        check: "-"
        if check not in findings
        else "pass"
        if findings[check].passed
        else "FAIL"
        for check in CHECKS
    }

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 44 to 53, read from the repository when this page was built.

def print_matrix(rows):
    print(f"{'':38}" + "".join(f"{check:>19}" for check in CHECKS))
    for name, marks in rows:
        print(f"{name:38}" + "".join(f"{marks[check]:>19}" for check in CHECKS))

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 56 to 59, read from the repository when this page was built.

def correct_steppers():
    """The verdicts of every check on the correct steppers."""
    correct = {name: verdicts(build, claim) for name, build, claim in CORRECT}
    print_matrix([(name, marks) for name, (_, marks) in correct.items()])
    return correct

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 62 to 66, read from the repository when this page was built.

                                                    clock             purity             limits  convergence-order     bounded-energy
velocity Verlet                                      pass               pass               pass               pass               pass
symplectic Euler                                     pass               pass               pass               pass               pass
Runge-Kutta 4                                        pass               pass               pass               pass                  -

Every applicable test passes, and the tolerances were fixed before this cell ran. The evidence behind a verdict is the finding. For velocity Verlet:

def verlet_findings(correct):
    """What each check measured for velocity Verlet."""
    for finding in correct["velocity Verlet"][0].values():
        print(
            f"{finding.check:18} measured {finding.measured:9.3g}"
            f"  tolerance {finding.tolerance:7.2g}  {finding.detail}"
        )
        print(f"{'':18} claim: {finding.claim}")
        print(f"{'':18} reproduce: {finding.reproduction}")

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 69 to 77, read from the repository when this page was built.

clock              measured         0  tolerance   1e-12  a step of 0.125 s advanced the time by 0.125 s
                   claim: a step of length dt advances the time by dt
                   reproduce: pbc.verification.checks.check_clock(build, claim) with build() returning a fresh stepper
purity             measured         0  tolerance   1e-12  the same step, repeated after other steps, differs by 0.0e+00
                   claim: the result of a step depends on its arguments alone
                   reproduce: pbc.verification.checks.check_purity(build, claim) with build() returning a fresh stepper
limits             measured  7.11e-15  tolerance   1e-12  uniform motion off by 8.9e-16 m, velocity under constant force by 7.1e-15 m/s, position by 3.6e-15 m (allowed 0.0e+00 m)
                   claim: no force gives uniform motion; a constant force gives the exact velocity, and the exact position for a method of order 2 or more
                   reproduce: pbc.verification.checks.check_limits(build, claim) with build() returning a fresh stepper
convergence-order  measured         2  tolerance     0.4  fitted order 2.00, claimed 2
                   claim: halving the step divides the error by 2 ** the claimed order
                   reproduce: pbc.verification.checks.check_convergence_order(build, claim) with build() returning a fresh stepper
bounded-energy     measured         1  tolerance     1.5  the energy band is 9.6e-03 in the first 5 periods and 9.6e-03 in the last 5
                   claim: the energy of the oscillator stays in a band that does not grow
                   reproduce: pbc.verification.checks.check_bounded_energy(build, claim) with build() returning a fresh stepper

The convergence test reports the fitted order, not the table behind it. The table is the same estimate one pair of steps at a time, and is how a reader looks at a result that is close to a tolerance:

counts = np.array([50, 100, 200, 400, 800])
steps = np.pi / counts

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 35 to 36, read from the repository when this page was built.

def verlet_convergence():
    """The error of velocity Verlet on the oscillator for each step, with
    the observed order."""
    errors = [oscillator_error(lambda: verlet_step, n) for n in counts]
    print(" steps per period   step (s)     error (m)   observed order")
    for n, (h, e, p) in zip(counts, convergence_table(steps, errors), strict=True):
        print(f"{n:17d}{h:11.5f}{e:14.3e}{'' if p is None else f'{p:15.2f}'}")

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 80 to 86, read from the repository when this page was built.

 steps per period   step (s)     error (m)   observed order
               50    0.06283     8.827e-03
              100    0.03142     2.205e-03           2.00
              200    0.01571     5.514e-04           2.00
              400    0.00785     1.378e-04           2.00
              800    0.00393     3.446e-05           2.00

The slope on logarithmic axes is the order. The figure puts the correct velocity Verlet beside two faulty steppers that pose as methods of higher order than they have:

def draw_verification_convergence() -> Figure:
    """The error against the step for two correct and two faulty steppers."""
    series = {
        "velocity Verlet": (verlet_step, "#00429d", "-"),
        "Runge-Kutta 4": (rk4_step, "#00734f", "-"),
        "verlet without the half": (verlet_without_the_half, "#93003a", "--"),
        "rk4 with equal weights": (rk4_with_equal_weights, "#b34700", "-."),
    }
    fig, ax = plt.subplots(figsize=(6.4, 4.2))
    for name, (step, color, style) in series.items():
        errors = [oscillator_error(lambda: step, n) for n in counts]  # noqa: B023
        ax.loglog(
            steps,
            errors,
            color=color,
            linestyle=style,
            marker="o",
            markersize=3,
            label=name,
        )
    ax.set_xlabel("time step (s)")
    ax.set_ylabel("largest error over two periods (m)")
    ax.legend(loc="lower right")
    return fig

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 89 to 112, read from the repository when this page was built.

Log-log plot of the error in metres against the time step in seconds. Velocity Verlet and Runge-Kutta 4 lie on straight lines of slope 2 and 4, the second far lower. Verlet without the half has slope 1 and lies above velocity Verlet. Runge-Kutta 4 with equal weights has slope 2 and lies on a line parallel to velocity Verlet, far above Runge-Kutta 4. — Figure 1: simulated Largest error over two periods against the time step for velocity Verlet, Runge-Kutta 4, and two faulty steppers. A fault shows as a slope below the claimed order.

Six faulty steppers

Each stepper below is wrong on purpose, in one way. They live in pbc.verification.faulty, a module that says so in its first line and that nothing which simulates may import. Each entry records which correct method it poses as, so that it is held to that method’s claim, and which tests must catch it. Three of them:

def verlet_without_the_half(
    state: State, acceleration: Acceleration, dt: float
) -> State:
    """Velocity Verlet whose position update has ``a dt**2`` where the
    Taylor expansion has ``a dt**2 / 2``."""
    a = acceleration(state.t, state.x, state.v)
    t = state.t + dt
    x = state.x + state.v * dt + a * dt**2
    a_next = acceleration(t, x, state.v + a * dt)
    return State(t=t, x=x, v=state.v + 0.5 * (a + a_next) * dt)

src/pbc/verification/faulty.py, lines 26 to 35, read from the repository when this page was built.

def verlet_with_a_stopped_clock(
    state: State, acceleration: Acceleration, dt: float
) -> State:
    """Velocity Verlet that returns the time it was given."""
    after = verlet_step(state, acceleration, dt)
    return State(t=state.t, x=after.x, v=after.v)

src/pbc/verification/faulty.py, lines 46 to 51, read from the repository when this page was built.

def verlet_that_remembers_the_acceleration() -> Stepper:
    """Velocity Verlet that keeps the acceleration of the last step it took
    and uses it at the start of the next one, whatever state it is given.

    Run in sequence on a force that does not depend on the velocity it is
    exact, which is why the fault survives every test that only runs a
    simulation."""
    remembered: list = []

    def step(state: State, acceleration: Acceleration, dt: float) -> State:
        a = remembered[0] if remembered else acceleration(state.t, state.x, state.v)
        t = state.t + dt
        x = state.x + state.v * dt + 0.5 * a * dt**2
        a_next = acceleration(t, x, state.v + a * dt)
        remembered[:] = [a_next]
        return State(t=t, x=x, v=state.v + 0.5 * (a + a_next) * dt)

    return step

src/pbc/verification/faulty.py, lines 86 to 103, read from the repository when this page was built.

The six run the same tests as the correct methods, under the claim of the method they pose as:

def gallery():
    """The verdicts of every check on the gallery of faulty steppers."""
    caught = {g.name: verdicts(g.build, g.claim) for g in GALLERY}
    print_matrix([(name, marks) for name, (_, marks) in caught.items()])
    print()
    for g in GALLERY:
        named = ", ".join(g.caught_by)
        print(f"{g.name}: {g.fault}\n    poses as {g.poses_as}; caught by {named}")

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 115 to 122, read from the repository when this page was built.

                                                    clock             purity             limits  convergence-order     bounded-energy
verlet without the half                              pass               pass               FAIL               FAIL               pass
symplectic Euler with the wrong sign                 pass               pass               FAIL               FAIL               FAIL
verlet with a stopped clock                          FAIL               pass               pass               pass               pass
rk4 with equal weights                               pass               pass               pass               FAIL                  -
explicit Euler as symplectic                         pass               pass               pass               pass               FAIL
verlet that remembers                                pass               FAIL               pass               pass               pass

verlet without the half: the position update has a dt² where the Taylor expansion has a dt²/2
    poses as velocity Verlet; caught by limits, convergence-order
symplectic Euler with the wrong sign: the velocity update subtracts the acceleration
    poses as symplectic Euler; caught by limits, convergence-order, bounded-energy
verlet with a stopped clock: the returned state keeps the old time
    poses as velocity Verlet; caught by clock
rk4 with equal weights: the four rates are averaged with equal weights
    poses as Runge-Kutta 4; caught by convergence-order
explicit Euler as symplectic: the position moves with the old velocity: the updates are in the wrong order
    poses as symplectic Euler; caught by bounded-energy
verlet that remembers: the acceleration of the last step is reused, whatever state comes next
    poses as velocity Verlet; caught by purity

Each fault fails the tests the gallery names for it and no others, and no correct method fails any test: the same tests separate the two groups. The pattern matters more than the count.

Verlet without the half has the wrong coefficient on a term that is of second order in the step. Everything it computes is plausible and the energy stays in its band. The order is wrong, which the convergence test sees, and so is the position under a constant force, which the limit test sees without a reference solution of the oscillator.
Symplectic Euler with the wrong sign runs the velocity the wrong way. Three tests fail: it is the fault any test would find, and the limit test finds it first, with no reference solution and a run of twenty steps.
Verlet with a stopped clock is exact in every physical test, because the oscillator does not depend on time. Only the interface test, clock, finds it, and it would corrupt every simulation with a time-dependent force.
Runge-Kutta 4 with equal weights is of order 2, not 4: accurate enough to look like a good method, until the fitted order is compared with the claim.
Explicit Euler as symplectic is a first-order method, as claimed, so the convergence test passes. Only the energy test sees that the updates are in the wrong order. A test of order alone would have approved it.
Verlet that remembers is exact when its steps are taken one after another, which is how every simulation takes them, so every test that only runs a simulation passes. Only purity, which asks for the same step again after others, finds the leaked state.

A test that cannot fail

A test that passes everything is not evidence. The first protection is the gallery: every test has a faulty stepper it must fail. The second is to claim the wrong thing about a correct method. Velocity Verlet claimed to be of order 1 or of order 4 must fail the convergence test, and the true claim must pass:

def claimed_orders():
    """The convergence check of velocity Verlet against three claimed
    orders."""
    for order in (1, 2, 4):
        findings = verify(lambda: verlet_step, Claim(order=order, bounded_energy=True))
        result = next(f for f in findings if f.check == "convergence-order")
        print(
            f"claimed order {order}: {'pass' if result.passed else 'FAIL'}"
            f"  ({result.detail})"
        )

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 125 to 134, read from the repository when this page was built.

claimed order 1: FAIL  (fitted order 2.00, claimed 1)
claimed order 2: pass  (fitted order 2.00, claimed 2)
claimed order 4: FAIL  (fitted order 2.00, claimed 4)

Only the true claim passes. A tolerance is part of the claim; how loose it may be is the subject of the first exercise.

An evaluation of a piece of code against the contract

Given a stepper and a claim, the evaluation is the table of findings, a verdict, and a sentence on what was covered. For a stepper someone else wrote, here the Runge-Kutta 4 with equal weights, claimed to be Runge-Kutta 4, the report is generated from the findings and nothing is judged by eye:

def report(name, build, claim):
    findings = verify(build, claim)
    failed = [f.check for f in findings if not f.passed]
    bound = "claimed" if claim.bounded_energy else "not claimed"
    print(f"subject: {name}; claim: order {claim.order}, energy bound {bound}")
    for f in findings:
        print(f"  {'pass' if f.passed else 'FAIL'}  {f.check:18} {f.detail}")
        print(f"{'':26}claim: {f.claim}")
        print(f"{'':26}reproduce: {f.reproduction}")
    print(
        "verdict: "
        + (
            "claim not supported by " + ", ".join(failed)
            if failed
            else "claim supported"
        )
    )
    print("covered: the oscillator at 50 to 800 steps per period")
    print("not covered: nonlinear forces, velocity-dependent forces, the Kepler orbit")

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 137 to 155, read from the repository when this page was built.

def report_on_equal_weights():
    """Exercise 1: the report on Runge-Kutta 4 with equal weights."""
    rk4_equal_weights = next(g for g in GALLERY if g.name == "rk4 with equal weights")
    report(rk4_equal_weights.name, rk4_equal_weights.build, rk4_equal_weights.claim)

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 158 to 161, read from the repository when this page was built.

subject: rk4 with equal weights; claim: order 4, energy bound not claimed
  pass  clock              a step of 0.125 s advanced the time by 0.125 s
                          claim: a step of length dt advances the time by dt
                          reproduce: pbc.verification.checks.check_clock(build, claim) with build() returning a fresh stepper
  pass  purity             the same step, repeated after other steps, differs by 0.0e+00
                          claim: the result of a step depends on its arguments alone
                          reproduce: pbc.verification.checks.check_purity(build, claim) with build() returning a fresh stepper
  pass  limits             uniform motion off by 8.9e-16 m, velocity under constant force by 7.1e-15 m/s, position by 5.3e-15 m (allowed 0.0e+00 m)
                          claim: no force gives uniform motion; a constant force gives the exact velocity, and the exact position for a method of order 2 or more
                          reproduce: pbc.verification.checks.check_limits(build, claim) with build() returning a fresh stepper
  FAIL  convergence-order  fitted order 2.00, claimed 4
                          claim: halving the step divides the error by 2 ** the claimed order
                          reproduce: pbc.verification.checks.check_convergence_order(build, claim) with build() returning a fresh stepper
verdict: claim not supported by convergence-order
covered: the oscillator at 50 to 800 steps per period
not covered: nonlinear forces, velocity-dependent forces, the Kepler orbit

The last two lines are part of the report. The tests were run on one model, and a claim about another is not supported until it is run there.

Exercises

Exercise 1. A tolerance that is too loose. The convergence test accepts a fitted order within 0.4 of the claim. Compute the fitted orders of velocity Verlet and of Verlet without the half, and judge each against the claim of order 2 with a margin of 0.4 and with a margin of 1.5. What does the loose margin do, and what does that say about choosing a tolerance after seeing the faulty result?

Solution to exercise 1

def margins():
    """Exercise 2: the fitted order of two steppers against two margins."""
    for name, step in (
        ("velocity Verlet", verlet_step),
        ("verlet without the half", verlet_without_the_half),
    ):
        errors = [oscillator_error(lambda: step, n) for n in counts]  # noqa: B023
        order = fitted_order(np.pi / counts, errors)
        marks = {m: "pass" if abs(order - 2) <= m else "FAIL" for m in (0.4, 1.5)}
        print(
            f"{name:26} fitted order {order:.2f}; margin 0.4: {marks[0.4]},"
            f" margin 1.5: {marks[1.5]}"
        )

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 164 to 176, read from the repository when this page was built.

velocity Verlet            fitted order 2.00; margin 0.4: pass, margin 1.5: pass
verlet without the half    fitted order 0.94; margin 0.4: FAIL, margin 1.5: pass

A margin of 1.5 passes both: the fault has made the method first order, and the distance from order 2 is smaller than the margin. The margin has to be chosen from the scatter of a correct method, here a few hundredths, before the faulty one is seen; fixed that way, it fails the fault with room to spare.

Exercise 2. A fault of your own. Write a stepper that is velocity Verlet except that the velocity update uses the acceleration at the start of the step only, and claim order 2 with a bounded energy for it. Predict which tests fail, then run verify.

Solution to exercise 2

def verlet_with_one_acceleration(state, acceleration, dt):
    a = acceleration(state.t, state.x, state.v)
    return State(
        t=state.t + dt,
        x=state.x + state.v * dt + 0.5 * a * dt**2,
        v=state.v + a * dt,
    )

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 179 to 185, read from the repository when this page was built.

def one_acceleration():
    """Exercise 3: the checks on a Verlet step that evaluates the
    acceleration once."""
    findings = verify(
        lambda: verlet_with_one_acceleration, Claim(order=2, bounded_energy=True)
    )
    for f in findings:
        print(f"{'pass' if f.passed else 'FAIL'}  {f.check:18} {f.detail}")

src/pbc/mechanics/experiments/verifying_scientific_code.py, lines 188 to 195, read from the repository when this page was built.

pass  clock              a step of 0.125 s advanced the time by 0.125 s
pass  purity             the same step, repeated after other steps, differs by 0.0e+00
pass  limits             uniform motion off by 8.9e-16 m, velocity under constant force by 7.1e-15 m/s, position by 3.6e-15 m (allowed 0.0e+00 m)
FAIL  convergence-order  fitted order 1.07, claimed 2
FAIL  bounded-energy     the energy band is 2.0e+01 in the first 5 periods and 3.3e+26 in the last 5

The velocity is updated as in explicit Euler, so the method is of order 1 and the convergence test fails. The clock and purity tests pass, and so does the limit test, which only needs a velocity update that is exact for a constant force, as this one is. The energy test fails as well, with a band that grows by many orders of magnitude over the run: a method that is second order in the position and first order in the velocity is not a symplectic method of any order. The three failures are one fault seen three ways, which is the pattern a verifier reads. The linear oscillator would not show a fault that only matters for a force that depends on the velocity; a test on the throw with drag of the earlier lessons would be the next one to add.

Run this experiment

This page was built from commit 6417e2ff90258bd2d2ba21d217eb7666fd24f680 of the repository.

In a clone at that commit, with the project’s Python packages installed, this command runs the experiment of this page:

uv run --locked --no-dev python -m pbc.mechanics.experiments.verifying_scientific_code results

Files, tests, and every command

The files behind it:

site/lessons/mechanics/09-verifying-scientific-code/index.qmd: this page, with every code cell.
src/pbc/verification/estimates.py: code the experiment runs.
src/pbc/verification/checks.py: code the experiment runs.
src/pbc/verification/faulty.py: code the experiment runs.
src/pbc/mechanics/experiments/verifying_scientific_code.py: code the experiment runs.
src/pbc/experiment.py: code the experiment runs.
tests/unit/test_verification.py: tests of that code.
tests/unit/test_verifying_scientific_code_experiment.py: tests of that code.

With Git and uv installed (see the setup page), run:

git clone https://github.com/taspinar/physics-by-construction.git
cd physics-by-construction
git checkout 6417e2ff90258bd2d2ba21d217eb7666fd24f680
uv sync --locked --no-dev
uv run --locked --no-dev python -m pbc.mechanics.experiments.verifying_scientific_code results

The experiment writes every figure of this page into results/, as <label>.png, and every number and table of the page, as the page shows it, into results/results.json. Its numbers match this page to the digits shown, and its figure files are the ones the page shows; digits beyond those can differ between machines, because floating-point libraries differ. Change the code of the experiment and run it again to see what changes.