FA-87056 / Physics integrator stability / Open access

Mixed error norm scales by the initial value instead of the current state · case 01

After the solution decays or grows the tolerance is effectively too tight or too loose.

Verified by executionVariant 1 · 8 checks per implementationDownload source bundle ↓JSON ↗

ROOT CAUSE

The denominator is max(1,|y0|) rather than max(1,|y|).

VERIFIED REPAIR

Scale by the magnitude of the current state.

Unsuccessful approach: Using 1+|y| over-relaxes the tolerance for small states.

Case contract

solve(y0, lam, s, t_end, dt0, tol): integrate y'=-lam*y+s to t_end with the embedded Euler/Heun pair. err=|y_heun-y_euler|/max(1,|y|) at the step start. Accept if err<=tol (propagate Heun); grow h by min(2, 0.9*(tol/err)^(1/2)) (2 when err=0); on reject shrink by max(0.2, 0.9*(tol/err)^(1/2)). Steps are clipped to t_end; at most 500 attempts. Return [y, accepted, rejected, last h] rounded.

Why this case matters

Game and robotics physics loops depend on integrator update order, step control and stabilization terms; a wrong decision point turns a stable simulation into drifting or exploding motion.

1 / The failure

Exit 1

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(y0, lam, s, t_end, dt0, tol):
    def f(y):
        return -lam * y + s
    t, y, h = 0.0, y0, dt0
    acc = rej = 0
    while t < t_end - 1e-12 and acc + rej < 500:
        h = min(h, t_end - t)
        k1 = f(y)
        k2 = f(y + h * k1)
        y_low = y + h * k1
        y_high = y + h / 2 * (k1 + k2)
        err = abs(y_high - y_low) / max(1.0, abs(y0))
        if err <= tol:
            t += h
            y = y_high
            acc += 1
            grow = 2.0 if err == 0 else min(2.0, 0.9 * (tol / err) ** 0.5)
            h = h * grow
        else:
            rej += 1
            h = h * max(0.2, 0.9 * (tol / err) ** 0.5)
    return [round(y, 6), acc, rej, round(h, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[1.2, 2.023, 1.796, 0.5, 0.05, 0.001], [1.001506, 12, 1, 0.012662]], [[-2.569, 4.482, -1.811, 1.5, 0.1, 0.001], [-0.406872, 68, 2, 0.019689]], [[0.355, 0.0, -0.581, 2.0, 0.3, 0.0001], [-0.807, 3, 0, 2.2]], [[-0.952, 1.695, -1.777, 2.0, 0.3, 0.01], [-1.043824, 6, 1, 0.806875]], [[-2.23, 0.286, -1.937, 2.0, 0.1, 0.0001], [-4.208794, 48, 1, 0.056624]], [[-0.696, 4.83, -1.861, 0.5, 0.7, 0.0001], [-0.41309, 64, 3, 0.004672]], [[0.106, 3.87, -0.733, 0.5, 0.7, 0.01], [-0.145157, 7, 2, 0.09505]], [[0.621, 1.843, 0.407, 2.0, 0.3, 0.001], [0.23109, 29, 2, 0.18722]]], [[[-0.207, 3.318, -0.298, 1.0, 0.05, 0.001], [-0.094257, 16, 1, 0.151659]], [[2.09, 3.725, -1.778, 1.0, 0.7, 0.0001], [-0.415384, 205, 4, 0.005852]], [[0.881, 2.451, 0.201, 1.5, 0.7, 0.0001], [0.102256, 121, 3, 0.035311]], [[-0.387, 1.826, -1.624, 1.5, 0.05, 0.01], [-0.855226, 11, 0, 0.239175]], [[1.324, 1.927, -0.437, 2.0, 0.3, 0.001], [-0.193674, 55, 2, 0.105477]], [[1.164, 0.0, 0.608, 0.5, 0.05, 0.0001], [1.468, 4, 0, 0.3]], [[1.539, 0.883, 0.651, 0.5, 0.7, 0.0001], [1.252849, 24, 2, 0.022362]], [[-0.677, 2.921, 1.193, 1.5, 0.7, 0.0001], [0.394822, 149, 4, 0.036265]]], [[[-1.057, 1.633, -1.528, 1.5, 0.3, 0.0001], [-0.946197, 41, 2, 0.02163]], [[0.32, 2.426, -0.767, 1.5, 0.3, 0.01], [-0.297783, 13, 1, 0.315679]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[-0.927, 2.194, -1.702, 2.0, 0.1, 0.01], [-0.779335, 9, 0, 0.051506]], [[-2.368, 1.168, -1.559, 1.0, 0.3, 0.001], [-1.656414, 17, 2, 0.075628]], [[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]], [[-1.862, 0.445, 0.515, 2.0, 0.1, 0.0001], [-0.082609, 92, 1, 0.025555]], [[2.17, 0.606, 1.006, 2.0, 0.7, 0.01], [1.81484, 5, 1, 0.239475]]], [[[-2.511, 1.07, 0.297, 0.5, 0.7, 0.0001], [-1.35563, 46, 3, 0.00122]], [[-0.242, 3.577, -0.4, 1.5, 0.7, 0.001], [-0.112625, 20, 2, 0.288449]], [[-2.172, 0.921, -1.832, 1.0, 0.3, 0.01], [-2.063818, 3, 0, 0.447421]], [[-2.345, 2.942, -1.888, 2.0, 0.7, 0.01], [-0.64795, 21, 2, 0.218574]], [[2.523, 2.59, 0.434, 1.5, 0.05, 0.0001], [0.215993, 189, 2, 0.001913]], [[-1.51, 3.253, -1.883, 1.5, 0.05, 0.0001], [-0.585952, 137, 2, 0.043884]], [[0.764, 2.831, -0.75, 0.5, 0.05, 0.01], [-0.013275, 9, 1, 0.080888]], [[-2.405, 0.0, -1.998, 1.5, 0.3, 0.0001], [-5.402, 3, 0, 1.2]]], [[[-0.374, 3.843, -1.056, 1.5, 0.7, 0.001], [-0.27525, 19, 2, 0.109434]], [[-0.799, 4.136, -0.079, 1.0, 0.7, 0.01], [-0.03315, 15, 2, 0.206325]], [[0.995, 2.051, -1.425, 0.5, 0.05, 0.001], [-0.088635, 27, 1, 0.024806]], [[-2.606, 3.176, -1.717, 1.5, 0.1, 0.001], [-0.55846, 60, 2, 0.071009]], [[1.917, 0.404, 1.543, 1.5, 0.1, 0.01], [2.779594, 6, 0, 0.239185]], [[-0.597, 0.0, -0.974, 2.0, 0.7, 0.01], [-2.545, 2, 0, 2.6]], [[-1.151, 1.35, -1.656, 0.5, 0.1, 0.001], [-1.188015, 5, 0, 0.046485]], [[-0.509, 0.0, 1.966, 1.5, 0.3, 0.001], [2.44, 3, 0, 1.2]]]]
for i, (args, expected) in enumerate(fixtures[N-1]):
    check("case %d" % i, solve(*args), expected)
print(json.dumps({"observations": observations, "passed": all(x["passed"] for x in observations)}, ensure_ascii=False))
raise SystemExit(0 if all(x["passed"] for x in observations) else 1)

Boundary fixture	Actual	Expected	Outcome
case 0	`[1.001514, 11, 1, 0.061887]`	`[1.001506, 12, 1, 0.012662]`	Failed
case 1	`[-0.407071, 49, 2, 0.138838]`	`[-0.406872, 68, 2, 0.019689]`	Failed
case 2	`[-0.807, 3, 0, 2.2]`	`[-0.807, 3, 0, 2.2]`	Passed
case 3	`[-1.043828, 6, 1, 0.817997]`	`[-1.043824, 6, 1, 0.806875]`	Failed
case 4	`[-4.208806, 57, 1, 0.017008]`	`[-4.208794, 48, 1, 0.056624]`	Failed
case 5	`[-0.41309, 64, 3, 0.004672]`	`[-0.41309, 64, 3, 0.004672]`	Passed
case 6	`[-0.145157, 7, 2, 0.09505]`	`[-0.145157, 7, 2, 0.09505]`	Passed
case 7	`[0.23109, 29, 2, 0.18722]`	`[0.23109, 29, 2, 0.18722]`	Passed

SHA-256 / 56c9af50c91c00cbfe6c5686d8ea08f99ce365dc8ab8fb9f8778499924bdf26e

2 / The unsuccessful fix

Exit 1

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(y0, lam, s, t_end, dt0, tol):
    def f(y):
        return -lam * y + s
    t, y, h = 0.0, y0, dt0
    acc = rej = 0
    while t < t_end - 1e-12 and acc + rej < 500:
        h = min(h, t_end - t)
        k1 = f(y)
        k2 = f(y + h * k1)
        y_low = y + h * k1
        y_high = y + h / 2 * (k1 + k2)
        err = abs(y_high - y_low) / (1.0 + abs(y))
        if err <= tol:
            t += h
            y = y_high
            acc += 1
            grow = 2.0 if err == 0 else min(2.0, 0.9 * (tol / err) ** 0.5)
            h = h * grow
        else:
            rej += 1
            h = h * max(0.2, 0.9 * (tol / err) ** 0.5)
    return [round(y, 6), acc, rej, round(h, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[1.2, 2.023, 1.796, 0.5, 0.05, 0.001], [1.001506, 12, 1, 0.012662]], [[-2.569, 4.482, -1.811, 1.5, 0.1, 0.001], [-0.406872, 68, 2, 0.019689]], [[0.355, 0.0, -0.581, 2.0, 0.3, 0.0001], [-0.807, 3, 0, 2.2]], [[-0.952, 1.695, -1.777, 2.0, 0.3, 0.01], [-1.043824, 6, 1, 0.806875]], [[-2.23, 0.286, -1.937, 2.0, 0.1, 0.0001], [-4.208794, 48, 1, 0.056624]], [[-0.696, 4.83, -1.861, 0.5, 0.7, 0.0001], [-0.41309, 64, 3, 0.004672]], [[0.106, 3.87, -0.733, 0.5, 0.7, 0.01], [-0.145157, 7, 2, 0.09505]], [[0.621, 1.843, 0.407, 2.0, 0.3, 0.001], [0.23109, 29, 2, 0.18722]]], [[[-0.207, 3.318, -0.298, 1.0, 0.05, 0.001], [-0.094257, 16, 1, 0.151659]], [[2.09, 3.725, -1.778, 1.0, 0.7, 0.0001], [-0.415384, 205, 4, 0.005852]], [[0.881, 2.451, 0.201, 1.5, 0.7, 0.0001], [0.102256, 121, 3, 0.035311]], [[-0.387, 1.826, -1.624, 1.5, 0.05, 0.01], [-0.855226, 11, 0, 0.239175]], [[1.324, 1.927, -0.437, 2.0, 0.3, 0.001], [-0.193674, 55, 2, 0.105477]], [[1.164, 0.0, 0.608, 0.5, 0.05, 0.0001], [1.468, 4, 0, 0.3]], [[1.539, 0.883, 0.651, 0.5, 0.7, 0.0001], [1.252849, 24, 2, 0.022362]], [[-0.677, 2.921, 1.193, 1.5, 0.7, 0.0001], [0.394822, 149, 4, 0.036265]]], [[[-1.057, 1.633, -1.528, 1.5, 0.3, 0.0001], [-0.946197, 41, 2, 0.02163]], [[0.32, 2.426, -0.767, 1.5, 0.3, 0.01], [-0.297783, 13, 1, 0.315679]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[-0.927, 2.194, -1.702, 2.0, 0.1, 0.01], [-0.779335, 9, 0, 0.051506]], [[-2.368, 1.168, -1.559, 1.0, 0.3, 0.001], [-1.656414, 17, 2, 0.075628]], [[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]], [[-1.862, 0.445, 0.515, 2.0, 0.1, 0.0001], [-0.082609, 92, 1, 0.025555]], [[2.17, 0.606, 1.006, 2.0, 0.7, 0.01], [1.81484, 5, 1, 0.239475]]], [[[-2.511, 1.07, 0.297, 0.5, 0.7, 0.0001], [-1.35563, 46, 3, 0.00122]], [[-0.242, 3.577, -0.4, 1.5, 0.7, 0.001], [-0.112625, 20, 2, 0.288449]], [[-2.172, 0.921, -1.832, 1.0, 0.3, 0.01], [-2.063818, 3, 0, 0.447421]], [[-2.345, 2.942, -1.888, 2.0, 0.7, 0.01], [-0.64795, 21, 2, 0.218574]], [[2.523, 2.59, 0.434, 1.5, 0.05, 0.0001], [0.215993, 189, 2, 0.001913]], [[-1.51, 3.253, -1.883, 1.5, 0.05, 0.0001], [-0.585952, 137, 2, 0.043884]], [[0.764, 2.831, -0.75, 0.5, 0.05, 0.01], [-0.013275, 9, 1, 0.080888]], [[-2.405, 0.0, -1.998, 1.5, 0.3, 0.0001], [-5.402, 3, 0, 1.2]]], [[[-0.374, 3.843, -1.056, 1.5, 0.7, 0.001], [-0.27525, 19, 2, 0.109434]], [[-0.799, 4.136, -0.079, 1.0, 0.7, 0.01], [-0.03315, 15, 2, 0.206325]], [[0.995, 2.051, -1.425, 0.5, 0.05, 0.001], [-0.088635, 27, 1, 0.024806]], [[-2.606, 3.176, -1.717, 1.5, 0.1, 0.001], [-0.55846, 60, 2, 0.071009]], [[1.917, 0.404, 1.543, 1.5, 0.1, 0.01], [2.779594, 6, 0, 0.239185]], [[-0.597, 0.0, -0.974, 2.0, 0.7, 0.01], [-2.545, 2, 0, 2.6]], [[-1.151, 1.35, -1.656, 0.5, 0.1, 0.001], [-1.188015, 5, 0, 0.046485]], [[-0.509, 0.0, 1.966, 1.5, 0.3, 0.001], [2.44, 3, 0, 1.2]]]]
for i, (args, expected) in enumerate(fixtures[N-1]):
    check("case %d" % i, solve(*args), expected)
print(json.dumps({"observations": observations, "passed": all(x["passed"] for x in observations)}, ensure_ascii=False))
raise SystemExit(0 if all(x["passed"] for x in observations) else 1)

Boundary fixture	Actual	Expected	Outcome
case 0	`[1.001644, 9, 0, 0.040046]`	`[1.001506, 12, 1, 0.012662]`	Failed
case 1	`[-0.406909, 54, 2, 0.131273]`	`[-0.406872, 68, 2, 0.019689]`	Failed
case 2	`[-0.807, 3, 0, 2.2]`	`[-0.807, 3, 0, 2.2]`	Passed
case 3	`[-1.042931, 5, 0, 0.725388]`	`[-1.043824, 6, 1, 0.806875]`	Failed
case 4	`[-4.208783, 42, 1, 0.047209]`	`[-4.208794, 48, 1, 0.056624]`	Failed
case 5	`[-0.4131, 52, 3, 0.011228]`	`[-0.41309, 64, 3, 0.004672]`	Failed
case 6	`[-0.144962, 7, 2, 0.057449]`	`[-0.145157, 7, 2, 0.09505]`	Failed
case 7	`[0.231136, 25, 2, 0.207631]`	`[0.23109, 29, 2, 0.18722]`	Failed

SHA-256 / 89732ef248bdf8fb23558bc47443dd7ffd161736787b261d6fb1b20c87dfe0d5

3 / The verified repair

Exit 0

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(y0, lam, s, t_end, dt0, tol):
    def f(y):
        return -lam * y + s
    t, y, h = 0.0, y0, dt0
    acc = rej = 0
    while t < t_end - 1e-12 and acc + rej < 500:
        h = min(h, t_end - t)
        k1 = f(y)
        k2 = f(y + h * k1)
        y_low = y + h * k1
        y_high = y + h / 2 * (k1 + k2)
        err = abs(y_high - y_low) / max(1.0, abs(y))
        if err <= tol:
            t += h
            y = y_high
            acc += 1
            grow = 2.0 if err == 0 else min(2.0, 0.9 * (tol / err) ** 0.5)
            h = h * grow
        else:
            rej += 1
            h = h * max(0.2, 0.9 * (tol / err) ** 0.5)
    return [round(y, 6), acc, rej, round(h, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[1.2, 2.023, 1.796, 0.5, 0.05, 0.001], [1.001506, 12, 1, 0.012662]], [[-2.569, 4.482, -1.811, 1.5, 0.1, 0.001], [-0.406872, 68, 2, 0.019689]], [[0.355, 0.0, -0.581, 2.0, 0.3, 0.0001], [-0.807, 3, 0, 2.2]], [[-0.952, 1.695, -1.777, 2.0, 0.3, 0.01], [-1.043824, 6, 1, 0.806875]], [[-2.23, 0.286, -1.937, 2.0, 0.1, 0.0001], [-4.208794, 48, 1, 0.056624]], [[-0.696, 4.83, -1.861, 0.5, 0.7, 0.0001], [-0.41309, 64, 3, 0.004672]], [[0.106, 3.87, -0.733, 0.5, 0.7, 0.01], [-0.145157, 7, 2, 0.09505]], [[0.621, 1.843, 0.407, 2.0, 0.3, 0.001], [0.23109, 29, 2, 0.18722]]], [[[-0.207, 3.318, -0.298, 1.0, 0.05, 0.001], [-0.094257, 16, 1, 0.151659]], [[2.09, 3.725, -1.778, 1.0, 0.7, 0.0001], [-0.415384, 205, 4, 0.005852]], [[0.881, 2.451, 0.201, 1.5, 0.7, 0.0001], [0.102256, 121, 3, 0.035311]], [[-0.387, 1.826, -1.624, 1.5, 0.05, 0.01], [-0.855226, 11, 0, 0.239175]], [[1.324, 1.927, -0.437, 2.0, 0.3, 0.001], [-0.193674, 55, 2, 0.105477]], [[1.164, 0.0, 0.608, 0.5, 0.05, 0.0001], [1.468, 4, 0, 0.3]], [[1.539, 0.883, 0.651, 0.5, 0.7, 0.0001], [1.252849, 24, 2, 0.022362]], [[-0.677, 2.921, 1.193, 1.5, 0.7, 0.0001], [0.394822, 149, 4, 0.036265]]], [[[-1.057, 1.633, -1.528, 1.5, 0.3, 0.0001], [-0.946197, 41, 2, 0.02163]], [[0.32, 2.426, -0.767, 1.5, 0.3, 0.01], [-0.297783, 13, 1, 0.315679]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[-0.927, 2.194, -1.702, 2.0, 0.1, 0.01], [-0.779335, 9, 0, 0.051506]], [[-2.368, 1.168, -1.559, 1.0, 0.3, 0.001], [-1.656414, 17, 2, 0.075628]], [[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]], [[-1.862, 0.445, 0.515, 2.0, 0.1, 0.0001], [-0.082609, 92, 1, 0.025555]], [[2.17, 0.606, 1.006, 2.0, 0.7, 0.01], [1.81484, 5, 1, 0.239475]]], [[[-2.511, 1.07, 0.297, 0.5, 0.7, 0.0001], [-1.35563, 46, 3, 0.00122]], [[-0.242, 3.577, -0.4, 1.5, 0.7, 0.001], [-0.112625, 20, 2, 0.288449]], [[-2.172, 0.921, -1.832, 1.0, 0.3, 0.01], [-2.063818, 3, 0, 0.447421]], [[-2.345, 2.942, -1.888, 2.0, 0.7, 0.01], [-0.64795, 21, 2, 0.218574]], [[2.523, 2.59, 0.434, 1.5, 0.05, 0.0001], [0.215993, 189, 2, 0.001913]], [[-1.51, 3.253, -1.883, 1.5, 0.05, 0.0001], [-0.585952, 137, 2, 0.043884]], [[0.764, 2.831, -0.75, 0.5, 0.05, 0.01], [-0.013275, 9, 1, 0.080888]], [[-2.405, 0.0, -1.998, 1.5, 0.3, 0.0001], [-5.402, 3, 0, 1.2]]], [[[-0.374, 3.843, -1.056, 1.5, 0.7, 0.001], [-0.27525, 19, 2, 0.109434]], [[-0.799, 4.136, -0.079, 1.0, 0.7, 0.01], [-0.03315, 15, 2, 0.206325]], [[0.995, 2.051, -1.425, 0.5, 0.05, 0.001], [-0.088635, 27, 1, 0.024806]], [[-2.606, 3.176, -1.717, 1.5, 0.1, 0.001], [-0.55846, 60, 2, 0.071009]], [[1.917, 0.404, 1.543, 1.5, 0.1, 0.01], [2.779594, 6, 0, 0.239185]], [[-0.597, 0.0, -0.974, 2.0, 0.7, 0.01], [-2.545, 2, 0, 2.6]], [[-1.151, 1.35, -1.656, 0.5, 0.1, 0.001], [-1.188015, 5, 0, 0.046485]], [[-0.509, 0.0, 1.966, 1.5, 0.3, 0.001], [2.44, 3, 0, 1.2]]]]
for i, (args, expected) in enumerate(fixtures[N-1]):
    check("case %d" % i, solve(*args), expected)
print(json.dumps({"observations": observations, "passed": all(x["passed"] for x in observations)}, ensure_ascii=False))
raise SystemExit(0 if all(x["passed"] for x in observations) else 1)

Boundary fixture	Actual	Expected	Outcome
case 0	`[1.001506, 12, 1, 0.012662]`	`[1.001506, 12, 1, 0.012662]`	Passed
case 1	`[-0.406872, 68, 2, 0.019689]`	`[-0.406872, 68, 2, 0.019689]`	Passed
case 2	`[-0.807, 3, 0, 2.2]`	`[-0.807, 3, 0, 2.2]`	Passed
case 3	`[-1.043824, 6, 1, 0.806875]`	`[-1.043824, 6, 1, 0.806875]`	Passed
case 4	`[-4.208794, 48, 1, 0.056624]`	`[-4.208794, 48, 1, 0.056624]`	Passed
case 5	`[-0.41309, 64, 3, 0.004672]`	`[-0.41309, 64, 3, 0.004672]`	Passed
case 6	`[-0.145157, 7, 2, 0.09505]`	`[-0.145157, 7, 2, 0.09505]`	Passed
case 7	`[0.23109, 29, 2, 0.18722]`	`[0.23109, 29, 2, 0.18722]`	Passed

SHA-256 / 1074f751b0eb499ec9c8b867057772da2a33b376dcc0faedd62dff2f0b478cc0

Verification & scope

A deterministic bounded teaching model with stipulated toy conventions and rounded float output; not a production physics engine or a proof of numerical stability. This reproducer isolates one failure mechanism. Results cover the supplied fixtures. Variants within a family share a test contract and should remain grouped when constructing evaluation splits. Related mechanisms with a shared evaluation_group must also remain together; these controlled models are not independent production incidents.

Observations recorded using Python 3.12.14 at 2026-09-29T14:50:55.286349+00:00.

Case digest / cc5d517dc10836a7f1f52cc8fdb7f6f0d727183ca40213413235144d1f2586eb