FA-87036 / Physics integrator stability / Open access
Step controller omits the error-order exponent on growth · case 01
Accepted steps grow erratically and trigger rejection cascades.
ROOT CAUSE
Growth uses 0.9*(tol/err) instead of the square root for a first-order error estimate.
VERIFIED REPAIR
Use exponent 1/(p+1)=1/2 for the Euler/Heun pair.
Unsuccessful approach: Exponent 1/3 corresponds to a second-order estimator, not this pair.
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(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))
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 = [[[[-2.661, 0.0, 0.567, 1.5, 0.05, 0.001], [-1.8105, 5, 0, 1.5]], [[1.455, 0.0, 0.495, 1.5, 0.3, 0.001], [2.1975, 3, 0, 1.2]], [[-0.678, 2.724, 0.55, 1.5, 0.05, 0.01], [0.185448, 16, 0, 0.164424]], [[1.301, 0.0, 1.346, 0.5, 0.1, 0.0001], [1.974, 3, 0, 0.4]], [[1.761, 1.735, -1.932, 1.5, 0.05, 0.01], [-0.898458, 21, 0, 0.089026]], [[1.155, 1.454, 1.298, 1.0, 0.7, 0.001], [0.954176, 14, 2, 0.10605]], [[0.621, 1.843, 0.407, 2.0, 0.3, 0.001], [0.23109, 29, 2, 0.18722]], [[1.121, 4.954, -0.6, 0.5, 0.3, 0.01], [-0.014719, 14, 2, 0.068659]]], [[[2.09, 3.725, -1.778, 1.0, 0.7, 0.0001], [-0.415384, 205, 4, 0.005852]], [[-0.876, 0.0, -0.17, 0.5, 0.05, 0.0001], [-0.961, 4, 0, 0.3]], [[-2.532, 0.766, -0.213, 0.5, 0.3, 0.001], [-1.815013, 9, 1, 0.056945]], [[-2.2, 0.0, -0.933, 1.0, 0.3, 0.001], [-3.133, 3, 0, 0.2]], [[2.431, 0.0, 0.71, 1.0, 0.1, 0.01], [3.141, 4, 0, 0.6]], [[2.544, 3.611, 0.037, 1.5, 0.05, 0.0001], [0.021528, 217, 2, 0.032019]], [[-2.813, 4.676, -0.613, 1.0, 0.1, 0.01], [-0.158154, 23, 1, 0.132111]], [[0.366, 0.277, -0.958, 2.0, 0.3, 0.0001], [-1.260745, 75, 2, 0.009378]]], [[[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]], [[-1.429, 4.357, 1.585, 2.0, 0.7, 0.001], [0.363343, 71, 3, 0.130298]], [[-2.296, 0.0, -1.609, 1.5, 0.3, 0.001], [-4.7095, 3, 0, 1.2]], [[1.267, 2.693, -0.039, 0.5, 0.7, 0.0001], [0.318914, 87, 3, 0.008112]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[-1.302, 4.391, -1.764, 1.5, 0.3, 0.01], [-0.404267, 18, 2, 0.40767]], [[-2.618, 0.0, -1.013, 1.0, 0.1, 0.001], [-3.631, 4, 0, 0.6]], [[1.412, 4.795, 1.541, 1.0, 0.3, 0.001], [0.330607, 50, 3, 0.043844]]], [[[-2.405, 0.0, -1.998, 1.5, 0.3, 0.0001], [-5.402, 3, 0, 1.2]], [[-1.257, 0.0, -1.478, 0.5, 0.1, 0.01], [-1.996, 3, 0, 0.4]], [[0.67, 3.55, -0.067, 1.0, 0.7, 0.01], [0.002935, 14, 2, 0.005435]], [[-0.836, 1.36, -1.381, 1.0, 0.05, 0.0001], [-0.969367, 34, 1, 0.042607]], [[-0.556, 1.399, -1.805, 0.5, 0.05, 0.0001], [-0.925415, 41, 1, 0.007374]], [[-2.5, 0.0, 0.727, 1.0, 0.05, 0.0001], [-1.773, 5, 0, 0.5]], [[0.551, 2.318, 0.186, 1.5, 0.7, 0.0001], [0.094812, 92, 3, 0.032361]], [[-1.416, 3.263, -0.863, 0.5, 0.05, 0.0001], [-0.489781, 93, 2, 0.002341]]], [[[1.902, 0.0, 1.399, 0.5, 0.1, 0.001], [2.6015, 3, 0, 0.4]], [[-0.195, 0.0, -1.909, 2.0, 0.1, 0.001], [-4.013, 5, 0, 1.0]], [[1.292, 4.788, 1.763, 1.5, 0.1, 0.01], [0.37002, 19, 1, 0.240185]], [[-0.799, 4.136, -0.079, 1.0, 0.7, 0.01], [-0.03315, 15, 2, 0.206325]], [[0.854, 0.751, 1.114, 0.5, 0.05, 0.01], [1.050601, 4, 0, 0.243552]], [[-0.597, 0.0, -0.974, 2.0, 0.7, 0.01], [-2.545, 2, 0, 2.6]], [[-2.168, 3.971, -0.843, 1.5, 0.1, 0.01], [-0.218747, 23, 1, 0.239273]], [[1.604, 3.262, 1.026, 1.0, 0.1, 0.0001], [0.363958, 141, 2, 0.004622]]]]
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.8105, 5, 0, 1.5] | [-1.8105, 5, 0, 1.5] | Passed |
| case 1 | [2.1975, 3, 0, 1.2] | [2.1975, 3, 0, 1.2] | Passed |
| case 2 | [0.184312, 15, 0, 0.006089] | [0.185448, 16, 0, 0.164424] | Failed |
| case 3 | [1.974, 3, 0, 0.4] | [1.974, 3, 0, 0.4] | Passed |
| case 4 | [-0.898084, 20, 4, 0.077825] | [-0.898458, 21, 0, 0.089026] | Failed |
| case 5 | [0.954211, 13, 2, 0.147149] | [0.954176, 14, 2, 0.10605] | Failed |
| case 6 | [0.23113, 27, 3, 0.235747] | [0.23109, 29, 2, 0.18722] | Failed |
| case 7 | [-0.014271, 13, 3, 0.080243] | [-0.014719, 14, 2, 0.068659] | Failed |
SHA-256 / 73cc9b6a066ad171eaad219911252bda968c7f07267a9f9da8cc403b0141ff54
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) / 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) ** (1 / 3))
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 = [[[[-2.661, 0.0, 0.567, 1.5, 0.05, 0.001], [-1.8105, 5, 0, 1.5]], [[1.455, 0.0, 0.495, 1.5, 0.3, 0.001], [2.1975, 3, 0, 1.2]], [[-0.678, 2.724, 0.55, 1.5, 0.05, 0.01], [0.185448, 16, 0, 0.164424]], [[1.301, 0.0, 1.346, 0.5, 0.1, 0.0001], [1.974, 3, 0, 0.4]], [[1.761, 1.735, -1.932, 1.5, 0.05, 0.01], [-0.898458, 21, 0, 0.089026]], [[1.155, 1.454, 1.298, 1.0, 0.7, 0.001], [0.954176, 14, 2, 0.10605]], [[0.621, 1.843, 0.407, 2.0, 0.3, 0.001], [0.23109, 29, 2, 0.18722]], [[1.121, 4.954, -0.6, 0.5, 0.3, 0.01], [-0.014719, 14, 2, 0.068659]]], [[[2.09, 3.725, -1.778, 1.0, 0.7, 0.0001], [-0.415384, 205, 4, 0.005852]], [[-0.876, 0.0, -0.17, 0.5, 0.05, 0.0001], [-0.961, 4, 0, 0.3]], [[-2.532, 0.766, -0.213, 0.5, 0.3, 0.001], [-1.815013, 9, 1, 0.056945]], [[-2.2, 0.0, -0.933, 1.0, 0.3, 0.001], [-3.133, 3, 0, 0.2]], [[2.431, 0.0, 0.71, 1.0, 0.1, 0.01], [3.141, 4, 0, 0.6]], [[2.544, 3.611, 0.037, 1.5, 0.05, 0.0001], [0.021528, 217, 2, 0.032019]], [[-2.813, 4.676, -0.613, 1.0, 0.1, 0.01], [-0.158154, 23, 1, 0.132111]], [[0.366, 0.277, -0.958, 2.0, 0.3, 0.0001], [-1.260745, 75, 2, 0.009378]]], [[[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]], [[-1.429, 4.357, 1.585, 2.0, 0.7, 0.001], [0.363343, 71, 3, 0.130298]], [[-2.296, 0.0, -1.609, 1.5, 0.3, 0.001], [-4.7095, 3, 0, 1.2]], [[1.267, 2.693, -0.039, 0.5, 0.7, 0.0001], [0.318914, 87, 3, 0.008112]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[-1.302, 4.391, -1.764, 1.5, 0.3, 0.01], [-0.404267, 18, 2, 0.40767]], [[-2.618, 0.0, -1.013, 1.0, 0.1, 0.001], [-3.631, 4, 0, 0.6]], [[1.412, 4.795, 1.541, 1.0, 0.3, 0.001], [0.330607, 50, 3, 0.043844]]], [[[-2.405, 0.0, -1.998, 1.5, 0.3, 0.0001], [-5.402, 3, 0, 1.2]], [[-1.257, 0.0, -1.478, 0.5, 0.1, 0.01], [-1.996, 3, 0, 0.4]], [[0.67, 3.55, -0.067, 1.0, 0.7, 0.01], [0.002935, 14, 2, 0.005435]], [[-0.836, 1.36, -1.381, 1.0, 0.05, 0.0001], [-0.969367, 34, 1, 0.042607]], [[-0.556, 1.399, -1.805, 0.5, 0.05, 0.0001], [-0.925415, 41, 1, 0.007374]], [[-2.5, 0.0, 0.727, 1.0, 0.05, 0.0001], [-1.773, 5, 0, 0.5]], [[0.551, 2.318, 0.186, 1.5, 0.7, 0.0001], [0.094812, 92, 3, 0.032361]], [[-1.416, 3.263, -0.863, 0.5, 0.05, 0.0001], [-0.489781, 93, 2, 0.002341]]], [[[1.902, 0.0, 1.399, 0.5, 0.1, 0.001], [2.6015, 3, 0, 0.4]], [[-0.195, 0.0, -1.909, 2.0, 0.1, 0.001], [-4.013, 5, 0, 1.0]], [[1.292, 4.788, 1.763, 1.5, 0.1, 0.01], [0.37002, 19, 1, 0.240185]], [[-0.799, 4.136, -0.079, 1.0, 0.7, 0.01], [-0.03315, 15, 2, 0.206325]], [[0.854, 0.751, 1.114, 0.5, 0.05, 0.01], [1.050601, 4, 0, 0.243552]], [[-0.597, 0.0, -0.974, 2.0, 0.7, 0.01], [-2.545, 2, 0, 2.6]], [[-2.168, 3.971, -0.843, 1.5, 0.1, 0.01], [-0.218747, 23, 1, 0.239273]], [[1.604, 3.262, 1.026, 1.0, 0.1, 0.0001], [0.363958, 141, 2, 0.004622]]]]
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.8105, 5, 0, 1.5] | [-1.8105, 5, 0, 1.5] | Passed |
| case 1 | [2.1975, 3, 0, 1.2] | [2.1975, 3, 0, 1.2] | Passed |
| case 2 | [0.185817, 17, 0, 0.219708] | [0.185448, 16, 0, 0.164424] | Failed |
| case 3 | [1.974, 3, 0, 0.4] | [1.974, 3, 0, 0.4] | Passed |
| case 4 | [-0.898778, 22, 0, 0.118055] | [-0.898458, 21, 0, 0.089026] | Failed |
| case 5 | [0.954152, 15, 2, 0.08185] | [0.954176, 14, 2, 0.10605] | Failed |
| case 6 | [0.231058, 32, 2, 0.000582] | [0.23109, 29, 2, 0.18722] | Failed |
| case 7 | [-0.01506, 15, 2, 0.060544] | [-0.014719, 14, 2, 0.068659] | Failed |
SHA-256 / 8c513a4ea958cbd81524eb78f0a0bd3f57bdb6f39d6244d7784a756241abbae3
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 = [[[[-2.661, 0.0, 0.567, 1.5, 0.05, 0.001], [-1.8105, 5, 0, 1.5]], [[1.455, 0.0, 0.495, 1.5, 0.3, 0.001], [2.1975, 3, 0, 1.2]], [[-0.678, 2.724, 0.55, 1.5, 0.05, 0.01], [0.185448, 16, 0, 0.164424]], [[1.301, 0.0, 1.346, 0.5, 0.1, 0.0001], [1.974, 3, 0, 0.4]], [[1.761, 1.735, -1.932, 1.5, 0.05, 0.01], [-0.898458, 21, 0, 0.089026]], [[1.155, 1.454, 1.298, 1.0, 0.7, 0.001], [0.954176, 14, 2, 0.10605]], [[0.621, 1.843, 0.407, 2.0, 0.3, 0.001], [0.23109, 29, 2, 0.18722]], [[1.121, 4.954, -0.6, 0.5, 0.3, 0.01], [-0.014719, 14, 2, 0.068659]]], [[[2.09, 3.725, -1.778, 1.0, 0.7, 0.0001], [-0.415384, 205, 4, 0.005852]], [[-0.876, 0.0, -0.17, 0.5, 0.05, 0.0001], [-0.961, 4, 0, 0.3]], [[-2.532, 0.766, -0.213, 0.5, 0.3, 0.001], [-1.815013, 9, 1, 0.056945]], [[-2.2, 0.0, -0.933, 1.0, 0.3, 0.001], [-3.133, 3, 0, 0.2]], [[2.431, 0.0, 0.71, 1.0, 0.1, 0.01], [3.141, 4, 0, 0.6]], [[2.544, 3.611, 0.037, 1.5, 0.05, 0.0001], [0.021528, 217, 2, 0.032019]], [[-2.813, 4.676, -0.613, 1.0, 0.1, 0.01], [-0.158154, 23, 1, 0.132111]], [[0.366, 0.277, -0.958, 2.0, 0.3, 0.0001], [-1.260745, 75, 2, 0.009378]]], [[[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]], [[-1.429, 4.357, 1.585, 2.0, 0.7, 0.001], [0.363343, 71, 3, 0.130298]], [[-2.296, 0.0, -1.609, 1.5, 0.3, 0.001], [-4.7095, 3, 0, 1.2]], [[1.267, 2.693, -0.039, 0.5, 0.7, 0.0001], [0.318914, 87, 3, 0.008112]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[-1.302, 4.391, -1.764, 1.5, 0.3, 0.01], [-0.404267, 18, 2, 0.40767]], [[-2.618, 0.0, -1.013, 1.0, 0.1, 0.001], [-3.631, 4, 0, 0.6]], [[1.412, 4.795, 1.541, 1.0, 0.3, 0.001], [0.330607, 50, 3, 0.043844]]], [[[-2.405, 0.0, -1.998, 1.5, 0.3, 0.0001], [-5.402, 3, 0, 1.2]], [[-1.257, 0.0, -1.478, 0.5, 0.1, 0.01], [-1.996, 3, 0, 0.4]], [[0.67, 3.55, -0.067, 1.0, 0.7, 0.01], [0.002935, 14, 2, 0.005435]], [[-0.836, 1.36, -1.381, 1.0, 0.05, 0.0001], [-0.969367, 34, 1, 0.042607]], [[-0.556, 1.399, -1.805, 0.5, 0.05, 0.0001], [-0.925415, 41, 1, 0.007374]], [[-2.5, 0.0, 0.727, 1.0, 0.05, 0.0001], [-1.773, 5, 0, 0.5]], [[0.551, 2.318, 0.186, 1.5, 0.7, 0.0001], [0.094812, 92, 3, 0.032361]], [[-1.416, 3.263, -0.863, 0.5, 0.05, 0.0001], [-0.489781, 93, 2, 0.002341]]], [[[1.902, 0.0, 1.399, 0.5, 0.1, 0.001], [2.6015, 3, 0, 0.4]], [[-0.195, 0.0, -1.909, 2.0, 0.1, 0.001], [-4.013, 5, 0, 1.0]], [[1.292, 4.788, 1.763, 1.5, 0.1, 0.01], [0.37002, 19, 1, 0.240185]], [[-0.799, 4.136, -0.079, 1.0, 0.7, 0.01], [-0.03315, 15, 2, 0.206325]], [[0.854, 0.751, 1.114, 0.5, 0.05, 0.01], [1.050601, 4, 0, 0.243552]], [[-0.597, 0.0, -0.974, 2.0, 0.7, 0.01], [-2.545, 2, 0, 2.6]], [[-2.168, 3.971, -0.843, 1.5, 0.1, 0.01], [-0.218747, 23, 1, 0.239273]], [[1.604, 3.262, 1.026, 1.0, 0.1, 0.0001], [0.363958, 141, 2, 0.004622]]]]
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.8105, 5, 0, 1.5] | [-1.8105, 5, 0, 1.5] | Passed |
| case 1 | [2.1975, 3, 0, 1.2] | [2.1975, 3, 0, 1.2] | Passed |
| case 2 | [0.185448, 16, 0, 0.164424] | [0.185448, 16, 0, 0.164424] | Passed |
| case 3 | [1.974, 3, 0, 0.4] | [1.974, 3, 0, 0.4] | Passed |
| case 4 | [-0.898458, 21, 0, 0.089026] | [-0.898458, 21, 0, 0.089026] | Passed |
| case 5 | [0.954176, 14, 2, 0.10605] | [0.954176, 14, 2, 0.10605] | Passed |
| case 6 | [0.23109, 29, 2, 0.18722] | [0.23109, 29, 2, 0.18722] | Passed |
| case 7 | [-0.014719, 14, 2, 0.068659] | [-0.014719, 14, 2, 0.068659] | Passed |
SHA-256 / 77269e56c86491657308c5e1426bd96346a6fe435ef73ddb2232daa9d247088e
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.102290+00:00.
Case digest / 1ab105353603aa665cf88b149bf01c133ac3aab6323112093725a1117be4803f