FA-87096 / Physics integrator stability / Open access

Backward Euler denominator omits the implicit stiffness term · case 01

Stiff springs at large dt overshoot and gain energy, defeating the implicit scheme.

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

ROOT CAUSE

The dt^2*k/m term is missing from the denominator.

VERIFIED REPAIR

Include 1+dt*c/m+dt^2*k/m from eliminating x_{n+1}.

Unsuccessful approach: Using dt*k/m has the wrong time-step power.

Case contract

solve(x0, v0, k, m, c, dt, steps): backward Euler for x'=v, m v'=-k x-c v in closed form v1=(v-dt*k*x/m)/(1+dt*c/m+dt^2*k/m), x1=x+dt*v1. Energy E=0.5 m v^2+0.5 k x^2; report the peak E_n/E_0 including step 0 (0 when E_0=0). Return [x, v, peak] rounded to 6.

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(x0, v0, k, m, c, dt, steps):
    x, v = x0, v0
    def energy(p, q):
        return 0.5 * m * q * q + 0.5 * k * p * p
    e0 = energy(x, v)
    peak = 1.0 if e0 > 0 else 0.0
    for _ in range(steps):
        v = (v - dt * k * x / m) / (1 + dt * c / m)
        x = x + dt * v
        if e0 > 0:
            peak = max(peak, energy(x, v) / e0)
    return [round(x, 6), round(v, 6), round(peak, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[-1.795, -0.563, 97.971, 1.205, 2.366, 0.01, 5], [-1.618428, 6.062195, 1.0]], [[-1.173, 0.504, 23.099, 0.515, 1.355, 0.01, 5], [-1.077545, 2.77492, 1.0]], [[-0.221, -1.064, 98.514, 1.675, 2.729, 0.05, 6], [-0.014738, 1.136191, 1.0]], [[1.763, 0.622, 84.046, 2.119, 2.812, 0.1, 3], [0.130435, -5.944283, 1.0]], [[0.191, 1.945, 34.85, 2.123, 1.81, 0.5, 2], [0.064229, -0.321773, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.228, -1.699, 31.924, 0.826, 2.927, 0.1, 3], [-0.077028, -0.539607, 1.0]], [[-1.468, -0.813, 5.598, 2.782, 0.907, 0.5, 6], [0.352262, 0.016453, 1.0]]], [[[0.194, -0.92, 25.988, 2.712, 0.886, 0.1, 8], [-0.209628, 0.163927, 1.0]], [[-1.593, 1.89, 51.296, 0.937, 1.405, 0.01, 5], [-1.382702, 5.629739, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[1.808, -0.132, 49.997, 2.669, 1.474, 0.01, 6], [1.731053, -2.081744, 1.0]], [[1.15, -0.015, 47.859, 1.721, 0.258, 0.05, 4], [0.522542, -4.482089, 1.0]], [[1.593, 1.157, 80.887, 1.679, 0.691, 0.05, 6], [-0.29915, -7.152492, 1.0]], [[1.605, 0.764, 71.519, 1.364, 1.033, 0.1, 8], [0.009419, 1.778105, 1.0]], [[-0.519, -0.968, 62.61, 2.039, 0.152, 0.1, 6], [0.217976, 0.56902, 1.0]]], [[[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[-1.201, -0.558, 89.531, 1.586, 0.182, 0.1, 5], [0.386107, -0.013533, 1.0]], [[-1.096, -0.407, 20.232, 0.596, 1.678, 0.1, 9], [0.100413, -0.614384, 1.0]], [[-1.106, -1.223, 60.531, 2.608, 1.36, 0.01, 8], [-1.10939, 0.861284, 1.0]], [[1.841, -0.172, 82.672, 2.727, 1.716, 0.05, 6], [0.034505, -7.469333, 1.0]], [[1.659, 0.2, 31.533, 0.495, 2.184, 0.1, 3], [0.104813, -4.59631, 1.0]], [[-1.724, -0.957, 58.152, 0.412, 2.15, 0.01, 9], [-0.983354, 13.094621, 1.0]], [[-0.643, 1.644, 71.947, 1.526, 2.926, 0.01, 9], [-0.398185, 3.378185, 1.0]]], [[[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.48, 0.698, 50.79, 1.825, 2.065, 0.05, 5], [0.243232, -1.683845, 1.0]], [[-1.762, 1.692, 87.1, 2.431, 2.792, 0.01, 12], [-1.153456, 7.390115, 1.0]], [[1.857, -0.996, 11.568, 1.719, 0.343, 0.5, 9], [-0.005974, -0.045732, 1.0]], [[-0.978, -0.995, 90.294, 1.951, 0.766, 0.1, 7], [0.166803, -1.198449, 1.0]], [[-1.031, -0.971, 69.316, 1.698, 1.184, 0.01, 12], [-0.833114, 3.639745, 1.0]], [[-0.828, -0.185, 67.12, 0.517, 2.786, 0.5, 10], [-0.0, -0.0, 1.0]], [[-1.767, -1.709, 7.751, 0.748, 0.809, 0.5, 5], [0.035361, -0.158519, 1.0]]], [[[1.966, 0.072, 60.969, 1.44, 1.411, 0.05, 13], [-0.558619, 3.747659, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.001, -0.283, 51.629, 2.856, 1.037, 0.01, 3], [-0.007389, -0.277456, 1.0]], [[1.647, -0.33, 15.134, 2.849, 1.27, 0.1, 7], [0.042104, -2.734236, 1.0]], [[0.544, 1.609, 65.688, 1.374, 0.719, 0.01, 10], [0.552951, -1.128422, 1.0]], [[1.48, -1.5, 47.647, 0.462, 0.24, 0.5, 7], [-1.4e-05, 2e-05, 1.0]], [[-0.907, 1.778, 31.733, 1.951, 2.217, 0.1, 10], [0.196352, -0.928892, 1.0]], [[0.318, -1.09, 68.068, 1.546, 0.961, 0.5, 5], [0.000624, -0.001624, 1.0]]]]
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.615333, 6.197261, 1.0]`	`[-1.618428, 6.062195, 1.0]`	Failed
case 1	`[-1.076525, 2.814533, 1.0]`	`[-1.077545, 2.77492, 1.0]`	Failed
case 2	`[0.037159, 1.517624, 1.0]`	`[-0.014738, 1.136191, 1.0]`	Failed
case 3	`[-0.625179, -9.088269, 1.0]`	`[0.130435, -5.944283, 1.0]`	Failed
case 4	`[-0.514147, -1.674846, 1.630464]`	`[0.064229, -0.321773, 1.0]`	Failed
case 5	`[0.0, 0.0, 0.0]`	`[0.0, 0.0, 0.0]`	Passed
case 6	`[-0.193138, -0.792642, 1.0]`	`[-0.077028, -0.539607, 1.0]`	Failed
case 7	`[0.616329, -0.775822, 1.0]`	`[0.352262, 0.016453, 1.0]`	Failed

SHA-256 / e73510ee79ba333b598a85875df1bcf21e56ccc3fd65acd37f52b431bb91e8da

2 / The unsuccessful fix

Exit 1

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(x0, v0, k, m, c, dt, steps):
    x, v = x0, v0
    def energy(p, q):
        return 0.5 * m * q * q + 0.5 * k * p * p
    e0 = energy(x, v)
    peak = 1.0 if e0 > 0 else 0.0
    for _ in range(steps):
        v = (v - dt * k * x / m) / (1 + dt * c / m + dt * k / m)
        x = x + dt * v
        if e0 > 0:
            peak = max(peak, energy(x, v) / e0)
    return [round(x, 6), round(v, 6), round(peak, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[-1.795, -0.563, 97.971, 1.205, 2.366, 0.01, 5], [-1.618428, 6.062195, 1.0]], [[-1.173, 0.504, 23.099, 0.515, 1.355, 0.01, 5], [-1.077545, 2.77492, 1.0]], [[-0.221, -1.064, 98.514, 1.675, 2.729, 0.05, 6], [-0.014738, 1.136191, 1.0]], [[1.763, 0.622, 84.046, 2.119, 2.812, 0.1, 3], [0.130435, -5.944283, 1.0]], [[0.191, 1.945, 34.85, 2.123, 1.81, 0.5, 2], [0.064229, -0.321773, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.228, -1.699, 31.924, 0.826, 2.927, 0.1, 3], [-0.077028, -0.539607, 1.0]], [[-1.468, -0.813, 5.598, 2.782, 0.907, 0.5, 6], [0.352262, 0.016453, 1.0]]], [[[0.194, -0.92, 25.988, 2.712, 0.886, 0.1, 8], [-0.209628, 0.163927, 1.0]], [[-1.593, 1.89, 51.296, 0.937, 1.405, 0.01, 5], [-1.382702, 5.629739, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[1.808, -0.132, 49.997, 2.669, 1.474, 0.01, 6], [1.731053, -2.081744, 1.0]], [[1.15, -0.015, 47.859, 1.721, 0.258, 0.05, 4], [0.522542, -4.482089, 1.0]], [[1.593, 1.157, 80.887, 1.679, 0.691, 0.05, 6], [-0.29915, -7.152492, 1.0]], [[1.605, 0.764, 71.519, 1.364, 1.033, 0.1, 8], [0.009419, 1.778105, 1.0]], [[-0.519, -0.968, 62.61, 2.039, 0.152, 0.1, 6], [0.217976, 0.56902, 1.0]]], [[[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[-1.201, -0.558, 89.531, 1.586, 0.182, 0.1, 5], [0.386107, -0.013533, 1.0]], [[-1.096, -0.407, 20.232, 0.596, 1.678, 0.1, 9], [0.100413, -0.614384, 1.0]], [[-1.106, -1.223, 60.531, 2.608, 1.36, 0.01, 8], [-1.10939, 0.861284, 1.0]], [[1.841, -0.172, 82.672, 2.727, 1.716, 0.05, 6], [0.034505, -7.469333, 1.0]], [[1.659, 0.2, 31.533, 0.495, 2.184, 0.1, 3], [0.104813, -4.59631, 1.0]], [[-1.724, -0.957, 58.152, 0.412, 2.15, 0.01, 9], [-0.983354, 13.094621, 1.0]], [[-0.643, 1.644, 71.947, 1.526, 2.926, 0.01, 9], [-0.398185, 3.378185, 1.0]]], [[[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.48, 0.698, 50.79, 1.825, 2.065, 0.05, 5], [0.243232, -1.683845, 1.0]], [[-1.762, 1.692, 87.1, 2.431, 2.792, 0.01, 12], [-1.153456, 7.390115, 1.0]], [[1.857, -0.996, 11.568, 1.719, 0.343, 0.5, 9], [-0.005974, -0.045732, 1.0]], [[-0.978, -0.995, 90.294, 1.951, 0.766, 0.1, 7], [0.166803, -1.198449, 1.0]], [[-1.031, -0.971, 69.316, 1.698, 1.184, 0.01, 12], [-0.833114, 3.639745, 1.0]], [[-0.828, -0.185, 67.12, 0.517, 2.786, 0.5, 10], [-0.0, -0.0, 1.0]], [[-1.767, -1.709, 7.751, 0.748, 0.809, 0.5, 5], [0.035361, -0.158519, 1.0]]], [[[1.966, 0.072, 60.969, 1.44, 1.411, 0.05, 13], [-0.558619, 3.747659, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.001, -0.283, 51.629, 2.856, 1.037, 0.01, 3], [-0.007389, -0.277456, 1.0]], [[1.647, -0.33, 15.134, 2.849, 1.27, 0.1, 7], [0.042104, -2.734236, 1.0]], [[0.544, 1.609, 65.688, 1.374, 0.719, 0.01, 10], [0.552951, -1.128422, 1.0]], [[1.48, -1.5, 47.647, 0.462, 0.24, 0.5, 7], [-1.4e-05, 2e-05, 1.0]], [[-0.907, 1.778, 31.733, 1.951, 2.217, 0.1, 10], [0.196352, -0.928892, 1.0]], [[0.318, -1.09, 68.068, 1.546, 0.961, 0.5, 5], [0.000624, -0.001624, 1.0]]]]
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.734399, 1.611463, 1.0]`	`[-1.618428, 6.062195, 1.0]`	Failed
case 1	`[-1.128877, 1.003246, 1.0]`	`[-1.077545, 2.77492, 1.0]`	Failed
case 2	`[-0.180176, 0.187271, 1.0]`	`[-0.014738, 1.136191, 1.0]`	Failed
case 3	`[1.340747, -1.453948, 1.0]`	`[0.130435, -5.944283, 1.0]`	Failed
case 4	`[0.122912, -0.175341, 1.0]`	`[0.064229, -0.321773, 1.0]`	Failed
case 5	`[0.0, 0.0, 0.0]`	`[0.0, 0.0, 0.0]`	Passed
case 6	`[0.139979, -0.158905, 1.0]`	`[-0.077028, -0.539607, 1.0]`	Failed
case 7	`[0.057799, 0.202336, 1.0]`	`[0.352262, 0.016453, 1.0]`	Failed

SHA-256 / 0222c1cc2d3028036e94786369d6a90bca629eb60403e9b0d594eae452a7df4a

3 / The verified repair

Exit 0

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(x0, v0, k, m, c, dt, steps):
    x, v = x0, v0
    def energy(p, q):
        return 0.5 * m * q * q + 0.5 * k * p * p
    e0 = energy(x, v)
    peak = 1.0 if e0 > 0 else 0.0
    for _ in range(steps):
        v = (v - dt * k * x / m) / (1 + dt * c / m + dt * dt * k / m)
        x = x + dt * v
        if e0 > 0:
            peak = max(peak, energy(x, v) / e0)
    return [round(x, 6), round(v, 6), round(peak, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[-1.795, -0.563, 97.971, 1.205, 2.366, 0.01, 5], [-1.618428, 6.062195, 1.0]], [[-1.173, 0.504, 23.099, 0.515, 1.355, 0.01, 5], [-1.077545, 2.77492, 1.0]], [[-0.221, -1.064, 98.514, 1.675, 2.729, 0.05, 6], [-0.014738, 1.136191, 1.0]], [[1.763, 0.622, 84.046, 2.119, 2.812, 0.1, 3], [0.130435, -5.944283, 1.0]], [[0.191, 1.945, 34.85, 2.123, 1.81, 0.5, 2], [0.064229, -0.321773, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.228, -1.699, 31.924, 0.826, 2.927, 0.1, 3], [-0.077028, -0.539607, 1.0]], [[-1.468, -0.813, 5.598, 2.782, 0.907, 0.5, 6], [0.352262, 0.016453, 1.0]]], [[[0.194, -0.92, 25.988, 2.712, 0.886, 0.1, 8], [-0.209628, 0.163927, 1.0]], [[-1.593, 1.89, 51.296, 0.937, 1.405, 0.01, 5], [-1.382702, 5.629739, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[1.808, -0.132, 49.997, 2.669, 1.474, 0.01, 6], [1.731053, -2.081744, 1.0]], [[1.15, -0.015, 47.859, 1.721, 0.258, 0.05, 4], [0.522542, -4.482089, 1.0]], [[1.593, 1.157, 80.887, 1.679, 0.691, 0.05, 6], [-0.29915, -7.152492, 1.0]], [[1.605, 0.764, 71.519, 1.364, 1.033, 0.1, 8], [0.009419, 1.778105, 1.0]], [[-0.519, -0.968, 62.61, 2.039, 0.152, 0.1, 6], [0.217976, 0.56902, 1.0]]], [[[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[-1.201, -0.558, 89.531, 1.586, 0.182, 0.1, 5], [0.386107, -0.013533, 1.0]], [[-1.096, -0.407, 20.232, 0.596, 1.678, 0.1, 9], [0.100413, -0.614384, 1.0]], [[-1.106, -1.223, 60.531, 2.608, 1.36, 0.01, 8], [-1.10939, 0.861284, 1.0]], [[1.841, -0.172, 82.672, 2.727, 1.716, 0.05, 6], [0.034505, -7.469333, 1.0]], [[1.659, 0.2, 31.533, 0.495, 2.184, 0.1, 3], [0.104813, -4.59631, 1.0]], [[-1.724, -0.957, 58.152, 0.412, 2.15, 0.01, 9], [-0.983354, 13.094621, 1.0]], [[-0.643, 1.644, 71.947, 1.526, 2.926, 0.01, 9], [-0.398185, 3.378185, 1.0]]], [[[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.48, 0.698, 50.79, 1.825, 2.065, 0.05, 5], [0.243232, -1.683845, 1.0]], [[-1.762, 1.692, 87.1, 2.431, 2.792, 0.01, 12], [-1.153456, 7.390115, 1.0]], [[1.857, -0.996, 11.568, 1.719, 0.343, 0.5, 9], [-0.005974, -0.045732, 1.0]], [[-0.978, -0.995, 90.294, 1.951, 0.766, 0.1, 7], [0.166803, -1.198449, 1.0]], [[-1.031, -0.971, 69.316, 1.698, 1.184, 0.01, 12], [-0.833114, 3.639745, 1.0]], [[-0.828, -0.185, 67.12, 0.517, 2.786, 0.5, 10], [-0.0, -0.0, 1.0]], [[-1.767, -1.709, 7.751, 0.748, 0.809, 0.5, 5], [0.035361, -0.158519, 1.0]]], [[[1.966, 0.072, 60.969, 1.44, 1.411, 0.05, 13], [-0.558619, 3.747659, 1.0]], [[0.0, 0.0, 5.0, 1.0, 0.5, 0.1, 3], [0.0, 0.0, 0.0]], [[0.001, -0.283, 51.629, 2.856, 1.037, 0.01, 3], [-0.007389, -0.277456, 1.0]], [[1.647, -0.33, 15.134, 2.849, 1.27, 0.1, 7], [0.042104, -2.734236, 1.0]], [[0.544, 1.609, 65.688, 1.374, 0.719, 0.01, 10], [0.552951, -1.128422, 1.0]], [[1.48, -1.5, 47.647, 0.462, 0.24, 0.5, 7], [-1.4e-05, 2e-05, 1.0]], [[-0.907, 1.778, 31.733, 1.951, 2.217, 0.1, 10], [0.196352, -0.928892, 1.0]], [[0.318, -1.09, 68.068, 1.546, 0.961, 0.5, 5], [0.000624, -0.001624, 1.0]]]]
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.618428, 6.062195, 1.0]`	`[-1.618428, 6.062195, 1.0]`	Passed
case 1	`[-1.077545, 2.77492, 1.0]`	`[-1.077545, 2.77492, 1.0]`	Passed
case 2	`[-0.014738, 1.136191, 1.0]`	`[-0.014738, 1.136191, 1.0]`	Passed
case 3	`[0.130435, -5.944283, 1.0]`	`[0.130435, -5.944283, 1.0]`	Passed
case 4	`[0.064229, -0.321773, 1.0]`	`[0.064229, -0.321773, 1.0]`	Passed
case 5	`[0.0, 0.0, 0.0]`	`[0.0, 0.0, 0.0]`	Passed
case 6	`[-0.077028, -0.539607, 1.0]`	`[-0.077028, -0.539607, 1.0]`	Passed
case 7	`[0.352262, 0.016453, 1.0]`	`[0.352262, 0.016453, 1.0]`	Passed

SHA-256 / 4855530ecbfc5f4cb74120493ffe4cfe3274888a359887116e79c837e5a0dde2

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

Case digest / 0ce1be13a2806600352060fe46c8458f93a915732ae067b17c1f7f7faca739bc