FA-87061 / Physics integrator stability / Open access

Exact steps never grow the step size · case 01

A constant-slope phase keeps the tiny initial step forever and wastes work.

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

ROOT CAUSE

The err == 0 branch uses growth factor 1.

VERIFIED REPAIR

Use the maximum growth factor when the estimate is exactly zero.

Unsuccessful approach: A growth factor of 1.5 still deviates from the contracted cap.

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 = 1.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.23, 0.286, -1.937, 2.0, 0.1, 0.0001], [-4.208794, 48, 1, 0.056624]], [[-0.603, 0.0, 1.691, 2.0, 0.1, 0.0001], [2.779, 5, 0, 1.0]], [[0.355, 0.0, -0.581, 2.0, 0.3, 0.0001], [-0.807, 3, 0, 2.2]], [[1.314, 3.637, 1.469, 1.0, 0.1, 0.001], [0.428088, 42, 2, 0.054341]], [[-1.56, 0.0, -1.68, 1.5, 0.05, 0.001], [-4.08, 5, 0, 1.5]], [[2.009, 1.907, -1.617, 0.5, 0.3, 0.01], [0.255071, 10, 1, 0.029016]], [[1.996, 1.795, -1.923, 1.0, 0.3, 0.0001], [-0.561722, 156, 3, 0.009874]], [[1.301, 0.0, 1.346, 0.5, 0.1, 0.0001], [1.974, 3, 0, 0.4]]], [[[0.881, 2.451, 0.201, 1.5, 0.7, 0.0001], [0.102256, 121, 3, 0.035311]], [[-0.911, 2.367, 1.851, 1.5, 0.7, 0.0001], [0.733371, 173, 4, 0.019948]], [[-2.2, 0.0, -0.933, 1.0, 0.3, 0.001], [-3.133, 3, 0, 0.2]], [[0.366, 0.277, -0.958, 2.0, 0.3, 0.0001], [-1.260745, 75, 2, 0.009378]], [[0.034, 0.0, -0.536, 1.5, 0.3, 0.001], [-0.77, 3, 0, 1.2]], [[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]], [[-2.161, 4.981, 0.821, 2.0, 0.05, 0.0001], [0.164701, 232, 2, 0.153215]]], [[[-2.605, 4.58, 1.505, 0.5, 0.1, 0.01], [0.029147, 18, 1, 0.021233]], [[-0.022, 0.0, -1.06, 2.0, 0.05, 0.001], [-2.142, 6, 0, 0.9]], [[-1.772, 4.412, -1.322, 0.5, 0.05, 0.0001], [-0.461828, 121, 2, 0.007093]], [[-2.296, 0.0, -1.609, 1.5, 0.3, 0.001], [-4.7095, 3, 0, 1.2]], [[-0.931, 0.0, 0.846, 1.0, 0.05, 0.01], [-0.085, 5, 0, 0.5]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[0.886, 0.0, -0.594, 2.0, 0.1, 0.0001], [-0.302, 5, 0, 1.0]], [[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]]], [[[-2.124, 0.0, -1.31, 0.5, 0.7, 0.01], [-2.779, 1, 0, 1.0]], [[0.551, 2.318, 0.186, 1.5, 0.7, 0.0001], [0.094812, 92, 3, 0.032361]], [[1.022, 0.992, -0.734, 2.0, 0.05, 0.001], [-0.497403, 43, 1, 0.07961]], [[-2.5, 0.0, 0.727, 1.0, 0.05, 0.0001], [-1.773, 5, 0, 0.5]], [[-1.664, 0.608, -1.437, 1.0, 0.05, 0.01], [-1.981366, 5, 0, 0.426033]], [[-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]], [[-2.419, 2.259, -1.694, 2.0, 0.1, 0.0001], [-0.768126, 161, 2, 0.040422]]], [[[2.915, 0.0, 1.496, 2.0, 0.1, 0.0001], [5.907, 5, 0, 1.0]], [[-1.909, 4.509, 1.302, 1.5, 0.05, 0.0001], [0.286193, 222, 2, 0.001448]], [[2.746, 4.491, -1.815, 2.0, 0.05, 0.001], [-0.403609, 86, 2, 0.189861]], [[-2.167, 0.0, 1.842, 2.0, 0.7, 0.0001], [1.517, 2, 0, 2.6]], [[1.902, 0.0, 1.399, 0.5, 0.1, 0.001], [2.6015, 3, 0, 0.4]], [[1.604, 3.262, 1.026, 1.0, 0.1, 0.0001], [0.363958, 141, 2, 0.004622]], [[-1.375, 1.447, 0.137, 2.0, 0.7, 0.01], [0.011153, 17, 2, 0.010327]], [[-1.008, 0.873, 1.67, 1.0, 0.05, 0.0001], [0.69286, 96, 2, 0.013152]]]]
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	`[-4.208794, 48, 1, 0.056624]`	`[-4.208794, 48, 1, 0.056624]`	Passed
case 1	`[2.779, 20, 0, 0.1]`	`[2.779, 5, 0, 1.0]`	Failed
case 2	`[-0.807, 7, 0, 0.2]`	`[-0.807, 3, 0, 2.2]`	Failed
case 3	`[0.428088, 42, 2, 0.054341]`	`[0.428088, 42, 2, 0.054341]`	Passed
case 4	`[-4.08, 30, 0, 0.05]`	`[-4.08, 5, 0, 1.5]`	Failed
case 5	`[0.255071, 10, 1, 0.029016]`	`[0.255071, 10, 1, 0.029016]`	Passed
case 6	`[-0.561722, 156, 3, 0.009874]`	`[-0.561722, 156, 3, 0.009874]`	Passed
case 7	`[1.974, 5, 0, 0.1]`	`[1.974, 3, 0, 0.4]`	Failed

SHA-256 / d5b7824d3324b3249e51127bb3c273ef1b23e70a7325a452740dcde9582d5bdd

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 = 1.5 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.23, 0.286, -1.937, 2.0, 0.1, 0.0001], [-4.208794, 48, 1, 0.056624]], [[-0.603, 0.0, 1.691, 2.0, 0.1, 0.0001], [2.779, 5, 0, 1.0]], [[0.355, 0.0, -0.581, 2.0, 0.3, 0.0001], [-0.807, 3, 0, 2.2]], [[1.314, 3.637, 1.469, 1.0, 0.1, 0.001], [0.428088, 42, 2, 0.054341]], [[-1.56, 0.0, -1.68, 1.5, 0.05, 0.001], [-4.08, 5, 0, 1.5]], [[2.009, 1.907, -1.617, 0.5, 0.3, 0.01], [0.255071, 10, 1, 0.029016]], [[1.996, 1.795, -1.923, 1.0, 0.3, 0.0001], [-0.561722, 156, 3, 0.009874]], [[1.301, 0.0, 1.346, 0.5, 0.1, 0.0001], [1.974, 3, 0, 0.4]]], [[[0.881, 2.451, 0.201, 1.5, 0.7, 0.0001], [0.102256, 121, 3, 0.035311]], [[-0.911, 2.367, 1.851, 1.5, 0.7, 0.0001], [0.733371, 173, 4, 0.019948]], [[-2.2, 0.0, -0.933, 1.0, 0.3, 0.001], [-3.133, 3, 0, 0.2]], [[0.366, 0.277, -0.958, 2.0, 0.3, 0.0001], [-1.260745, 75, 2, 0.009378]], [[0.034, 0.0, -0.536, 1.5, 0.3, 0.001], [-0.77, 3, 0, 1.2]], [[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]], [[-2.161, 4.981, 0.821, 2.0, 0.05, 0.0001], [0.164701, 232, 2, 0.153215]]], [[[-2.605, 4.58, 1.505, 0.5, 0.1, 0.01], [0.029147, 18, 1, 0.021233]], [[-0.022, 0.0, -1.06, 2.0, 0.05, 0.001], [-2.142, 6, 0, 0.9]], [[-1.772, 4.412, -1.322, 0.5, 0.05, 0.0001], [-0.461828, 121, 2, 0.007093]], [[-2.296, 0.0, -1.609, 1.5, 0.3, 0.001], [-4.7095, 3, 0, 1.2]], [[-0.931, 0.0, 0.846, 1.0, 0.05, 0.01], [-0.085, 5, 0, 0.5]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[0.886, 0.0, -0.594, 2.0, 0.1, 0.0001], [-0.302, 5, 0, 1.0]], [[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]]], [[[-2.124, 0.0, -1.31, 0.5, 0.7, 0.01], [-2.779, 1, 0, 1.0]], [[0.551, 2.318, 0.186, 1.5, 0.7, 0.0001], [0.094812, 92, 3, 0.032361]], [[1.022, 0.992, -0.734, 2.0, 0.05, 0.001], [-0.497403, 43, 1, 0.07961]], [[-2.5, 0.0, 0.727, 1.0, 0.05, 0.0001], [-1.773, 5, 0, 0.5]], [[-1.664, 0.608, -1.437, 1.0, 0.05, 0.01], [-1.981366, 5, 0, 0.426033]], [[-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]], [[-2.419, 2.259, -1.694, 2.0, 0.1, 0.0001], [-0.768126, 161, 2, 0.040422]]], [[[2.915, 0.0, 1.496, 2.0, 0.1, 0.0001], [5.907, 5, 0, 1.0]], [[-1.909, 4.509, 1.302, 1.5, 0.05, 0.0001], [0.286193, 222, 2, 0.001448]], [[2.746, 4.491, -1.815, 2.0, 0.05, 0.001], [-0.403609, 86, 2, 0.189861]], [[-2.167, 0.0, 1.842, 2.0, 0.7, 0.0001], [1.517, 2, 0, 2.6]], [[1.902, 0.0, 1.399, 0.5, 0.1, 0.001], [2.6015, 3, 0, 0.4]], [[1.604, 3.262, 1.026, 1.0, 0.1, 0.0001], [0.363958, 141, 2, 0.004622]], [[-1.375, 1.447, 0.137, 2.0, 0.7, 0.01], [0.011153, 17, 2, 0.010327]], [[-1.008, 0.873, 1.67, 1.0, 0.05, 0.0001], [0.69286, 96, 2, 0.013152]]]]
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	`[-4.208794, 48, 1, 0.056624]`	`[-4.208794, 48, 1, 0.056624]`	Passed
case 1	`[2.779, 6, 0, 1.021875]`	`[2.779, 5, 0, 1.0]`	Failed
case 2	`[-0.807, 4, 0, 0.8625]`	`[-0.807, 3, 0, 2.2]`	Failed
case 3	`[0.428088, 42, 2, 0.054341]`	`[0.428088, 42, 2, 0.054341]`	Passed
case 4	`[-4.08, 7, 0, 0.691406]`	`[-4.08, 5, 0, 1.5]`	Failed
case 5	`[0.255071, 10, 1, 0.029016]`	`[0.255071, 10, 1, 0.029016]`	Passed
case 6	`[-0.561722, 156, 3, 0.009874]`	`[-0.561722, 156, 3, 0.009874]`	Passed
case 7	`[1.974, 4, 0, 0.0375]`	`[1.974, 3, 0, 0.4]`	Failed

SHA-256 / 98db207aef38976f4d26e02be0b9c75039e9a4b85f70a2397004c9d4801dff68

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.23, 0.286, -1.937, 2.0, 0.1, 0.0001], [-4.208794, 48, 1, 0.056624]], [[-0.603, 0.0, 1.691, 2.0, 0.1, 0.0001], [2.779, 5, 0, 1.0]], [[0.355, 0.0, -0.581, 2.0, 0.3, 0.0001], [-0.807, 3, 0, 2.2]], [[1.314, 3.637, 1.469, 1.0, 0.1, 0.001], [0.428088, 42, 2, 0.054341]], [[-1.56, 0.0, -1.68, 1.5, 0.05, 0.001], [-4.08, 5, 0, 1.5]], [[2.009, 1.907, -1.617, 0.5, 0.3, 0.01], [0.255071, 10, 1, 0.029016]], [[1.996, 1.795, -1.923, 1.0, 0.3, 0.0001], [-0.561722, 156, 3, 0.009874]], [[1.301, 0.0, 1.346, 0.5, 0.1, 0.0001], [1.974, 3, 0, 0.4]]], [[[0.881, 2.451, 0.201, 1.5, 0.7, 0.0001], [0.102256, 121, 3, 0.035311]], [[-0.911, 2.367, 1.851, 1.5, 0.7, 0.0001], [0.733371, 173, 4, 0.019948]], [[-2.2, 0.0, -0.933, 1.0, 0.3, 0.001], [-3.133, 3, 0, 0.2]], [[0.366, 0.277, -0.958, 2.0, 0.3, 0.0001], [-1.260745, 75, 2, 0.009378]], [[0.034, 0.0, -0.536, 1.5, 0.3, 0.001], [-0.77, 3, 0, 1.2]], [[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]], [[-2.161, 4.981, 0.821, 2.0, 0.05, 0.0001], [0.164701, 232, 2, 0.153215]]], [[[-2.605, 4.58, 1.505, 0.5, 0.1, 0.01], [0.029147, 18, 1, 0.021233]], [[-0.022, 0.0, -1.06, 2.0, 0.05, 0.001], [-2.142, 6, 0, 0.9]], [[-1.772, 4.412, -1.322, 0.5, 0.05, 0.0001], [-0.461828, 121, 2, 0.007093]], [[-2.296, 0.0, -1.609, 1.5, 0.3, 0.001], [-4.7095, 3, 0, 1.2]], [[-0.931, 0.0, 0.846, 1.0, 0.05, 0.01], [-0.085, 5, 0, 0.5]], [[-0.18, 0.0, 0.811, 2.0, 0.05, 0.0001], [1.442, 6, 0, 0.9]], [[0.886, 0.0, -0.594, 2.0, 0.1, 0.0001], [-0.302, 5, 0, 1.0]], [[0.299, 2.489, -1.677, 0.5, 0.3, 0.001], [-0.393348, 24, 2, 0.029906]]], [[[-2.124, 0.0, -1.31, 0.5, 0.7, 0.01], [-2.779, 1, 0, 1.0]], [[0.551, 2.318, 0.186, 1.5, 0.7, 0.0001], [0.094812, 92, 3, 0.032361]], [[1.022, 0.992, -0.734, 2.0, 0.05, 0.001], [-0.497403, 43, 1, 0.07961]], [[-2.5, 0.0, 0.727, 1.0, 0.05, 0.0001], [-1.773, 5, 0, 0.5]], [[-1.664, 0.608, -1.437, 1.0, 0.05, 0.01], [-1.981366, 5, 0, 0.426033]], [[-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]], [[-2.419, 2.259, -1.694, 2.0, 0.1, 0.0001], [-0.768126, 161, 2, 0.040422]]], [[[2.915, 0.0, 1.496, 2.0, 0.1, 0.0001], [5.907, 5, 0, 1.0]], [[-1.909, 4.509, 1.302, 1.5, 0.05, 0.0001], [0.286193, 222, 2, 0.001448]], [[2.746, 4.491, -1.815, 2.0, 0.05, 0.001], [-0.403609, 86, 2, 0.189861]], [[-2.167, 0.0, 1.842, 2.0, 0.7, 0.0001], [1.517, 2, 0, 2.6]], [[1.902, 0.0, 1.399, 0.5, 0.1, 0.001], [2.6015, 3, 0, 0.4]], [[1.604, 3.262, 1.026, 1.0, 0.1, 0.0001], [0.363958, 141, 2, 0.004622]], [[-1.375, 1.447, 0.137, 2.0, 0.7, 0.01], [0.011153, 17, 2, 0.010327]], [[-1.008, 0.873, 1.67, 1.0, 0.05, 0.0001], [0.69286, 96, 2, 0.013152]]]]
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	`[-4.208794, 48, 1, 0.056624]`	`[-4.208794, 48, 1, 0.056624]`	Passed
case 1	`[2.779, 5, 0, 1.0]`	`[2.779, 5, 0, 1.0]`	Passed
case 2	`[-0.807, 3, 0, 2.2]`	`[-0.807, 3, 0, 2.2]`	Passed
case 3	`[0.428088, 42, 2, 0.054341]`	`[0.428088, 42, 2, 0.054341]`	Passed
case 4	`[-4.08, 5, 0, 1.5]`	`[-4.08, 5, 0, 1.5]`	Passed
case 5	`[0.255071, 10, 1, 0.029016]`	`[0.255071, 10, 1, 0.029016]`	Passed
case 6	`[-0.561722, 156, 3, 0.009874]`	`[-0.561722, 156, 3, 0.009874]`	Passed
case 7	`[1.974, 3, 0, 0.4]`	`[1.974, 3, 0, 0.4]`	Passed

SHA-256 / 88f6bd6e72fa51312b1426b3e540f1cf10666e21a9a99a0e6171f657c4c7c75f

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.287598+00:00.

Case digest / b5c466a5404bb0711b6cace5c82eade024930e55d0d2c98b89a5c45a01427d5d