Skip to content

Convergence-study pitfalls

A convergence study can go wrong in ways that make a correct code look broken, or a broken one look fine. These came up again and again when testing public code, both in the codes' own studies and in the studies written to check them.

Comparing at the wrong time

int(t_end / h) can come out one short. In floating point, 2*pi / (2*pi/200) is 199.99999999999997, and int makes it 199. The run stops one step before t_end, its error is then O(h) instead of O(h^4), and an order-4 method looks first order.

Comparing arrays index by index, or taking the maximum over "the last N samples", moves the comparison time with h.

trajectory_order_study takes the whole trajectories and compares every run at one fixed time, failing loudly when a run has no sample there:

"""Convergence-study pitfalls: a run that stops one step short, and comparing at a whole period."""
import math

import numpy as np

from physverify.studies import trajectory_order_study


def rk4_run(n_steps, h):
    """RK4 on x'' = -x; returns the times and states of every step."""
    t, y = [0.0], [np.array([1.0, 0.0])]

    def f(y):
        return np.array([y[1], -y[0]])

    for k in range(n_steps):
        x = y[-1]
        k1 = f(x)
        k2 = f(x + h / 2 * k1)
        k3 = f(x + h / 2 * k2)
        k4 = f(x + h * k3)
        y.append(x + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4))
        t.append((k + 1) * h)
    return np.array(t), np.array(y)


T = 2 * math.pi


def exact(t):
    return np.array([math.cos(t), -math.sin(t)])


ns = (40, 80, 160, 320)
good = [rk4_run(n, T / n) for n in ns]
print("at 3T/4:", trajectory_order_study(good, exact, 0.75 * T, 4).message)

# the classic mistake: n = int(t_end / h) can come out one short in floating point
short = [rk4_run(n - 1 if n == 80 else n, T / n) for n in ns]
r = trajectory_order_study(short, exact, T, 4)
print("one run short:", r.status, "-", r.message)

r = trajectory_order_study(good, exact, T, 4, period=T)
print("at a whole period:", r.status, "-", r.message)

Output:

at 3T/4: observed order 4.00 (expected 4)
one run short: FAIL - the comparison time is not on every grid: run 1 (h = 0.0785) has no sample at t = 6.28319; it ends at t = 6.20464549, 1 step(s) short (a loop over int(t_end / h) steps?)
at a whole period: PASS - observed order 4.00 (expected 4); t = 6.28319 is 1 whole period(s), where phase errors can cancel and show a higher order than the method has: compare at a time that is not a multiple too

Comparing at a whole number of periods

At a whole period, phase errors of an oscillator can cancel, and a method can look more accurate than it is: a symplectic Euler method showed order 2 there in a real study. The last line above is a note, not a failure: compare at a time that is not a multiple of the period too.

Loops that test the time instead of counting steps

while t < t_end: t += h takes one step too many whenever rounding leaves t just below t_end, and the right number otherwise. time_grid_check runs the code with several step sizes and checks the number of steps and the final time:

"""`while t < t_end: t += h` can take one step too many, depending on rounding."""
from physverify.studies import time_grid_check


def times_while_loop(n, t_end=1.0):
    h, t, out = t_end / n, 0.0, [0.0]
    while t < t_end:  # the common loop
        t += h
        out.append(t)
    return out


def times_counted(n, t_end=1.0):
    h = t_end / n
    return [k * h for k in range(n + 1)]  # exactly n steps


ns = [4, 8, 10, 3]
print("while loop:", time_grid_check(times_while_loop, 0.0, 1.0, ns).message)
print("counted:   ", time_grid_check(times_counted, 0.0, 1.0, ns).message)

Output:

while loop: h=0.1: 11 steps ending at t=1.1 (expected 10 ending at 1)
counted:    ends at t=1 with the expected number of steps for every h tested

The fix is to count: for k in range(n): t = (k + 1) * h.

Errors already at round-off

If the coarsest step already gives an error of 1e-13, the error cannot shrink any more and the study reports INCONCLUSIVE or EXACT, not a failure. Use coarser steps. The message says when this happened.

Absolute tolerances

A check that requires "error below 1e-6 at the finest step" depends on the problem's scale and on how fine the finest step is. The observed order does not. Prefer the order; if an absolute tolerance fails while the observed order is right, the code is most likely fine.

Start-up transients

When fitting a decay or relaxation rate from a run that starts from equilibrium, skip the start-up transient. In one real code, fitting from t = 0 gave a viscosity 2% off; skipping the transient gave 0.01%.