FA-86916 / Physics integrator stability / Open access

Underdamped eigenvalues use twice the imaginary part · case 01

Oscillatory systems report too much amplification and are flagged unstable early.

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

ROOT CAUSE

The complex branch forgets to halve sqrt(-disc).

VERIFIED REPAIR

Imaginary part is sqrt(4k/m-(c/m)^2)/2.

Unsuccessful approach: Moving the halving to the real part instead leaves the damping rate wrong.

Case contract

solve(k, m, c, dt): classify explicit Euler on x'=v, v'=(-k*x-c*v)/m. Eigenvalues are the roots of s^2+(c/m)s+k/m; amplification is max |1+dt*lambda|. Return ["marginal" if |amp-1|<=1e-9, else "stable" if amp<1, else "unstable", amp 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(k, m, c, dt):
    a = c / m
    b = k / m
    disc = a * a - 4 * b
    if disc >= 0:
        r = math.sqrt(disc)
        lams = [complex((-a + r) / 2, 0), complex((-a - r) / 2, 0)]
    else:
        r = math.sqrt(-disc)
        lams = [complex(-a / 2, r), complex(-a / 2, -r)]
    amp = max(abs(1 + dt * lam) for lam in lams)
    if abs(amp - 1) <= 1e-9:
        label = 'marginal'
    elif amp < 1:
        label = 'stable'
    else:
        label = 'unstable'
    return [label, round(amp, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[33.665, 3.374, 3.698, 0.05], ['stable', 0.984958]], [[17.863, 3.409, 0.738, 0.4], ['unstable', 1.323555]], [[15.545, 2.583, 6.586, 1.0], ['unstable', 2.11387]], [[46.563, 1.397, 1.187, 0.01], ['stable', 0.997415]], [[39.672, 0.379, 2.933, 1.0], ['unstable', 9.896296]], [[0.0, 3.87, 4.121, 1.0], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]]], [[[23.283, 2.64, 7.844, 0.1], ['stable', 0.889422]], [[24.43, 0.273, 3.575, 0.1], ['stable', 0.76508]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[8.452, 3.883, 5.856, 0.4], ['stable', 0.863147]], [[47.405, 3.962, 4.696, 0.05], ['stable', 0.985215]], [[0.0, 3.884, 4.78, 0.2], ['marginal', 1.0]], [[56.015, 2.524, 6.052, 1.0], ['unstable', 4.560172]]], [[[15.636, 2.616, 6.979, 0.8], ['unstable', 1.640448]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[8.52, 3.836, 5.13, 0.4], ['stable', 0.90578]], [[29.423, 0.607, 3.273, 0.1], ['stable', 0.972378]], [[0.0, 3.996, 2.36, 1.0], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[18.776, 1.514, 0.866, 0.01], ['stable', 0.997758]], [[17.34, 2.005, 4.299, 0.01], ['stable', 0.989658]]], [[[9.079, 3.699, 5.728, 0.05], ['stable', 0.963696]], [[15.373, 0.488, 0.687, 0.8], ['unstable', 4.476057]], [[0.0, 0.88, 2.713, 0.01], ['marginal', 1.0]], [[29.305, 3.778, 6.728, 0.2], ['stable', 0.976782]], [[41.131, 0.314, 6.187, 0.4], ['unstable', 3.751925]], [[3.996, 2.055, 0.915, 0.05], ['stable', 0.991261]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]]], [[[49.306, 2.663, 2.331, 0.2], ['unstable', 1.251216]], [[58.722, 3.822, 5.191, 0.8], ['unstable', 3.121945]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.689, 5.168, 0.01], ['marginal', 1.0]], [[39.004, 1.806, 2.569, 0.01], ['stable', 0.993949]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[36.246, 3.029, 2.637, 0.2], ['unstable', 1.142163]], [[13.954, 1.677, 2.862, 1.0], ['unstable', 2.759383]]]]
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	`['unstable', 1.0008]`	`['unstable', 1.0008]`	Passed
case 1	`['unstable', 1.021139]`	`['stable', 0.984958]`	Failed
case 2	`['unstable', 2.064304]`	`['unstable', 1.323555]`	Failed
case 3	`['unstable', 4.200848]`	`['unstable', 2.11387]`	Failed
case 4	`['unstable', 1.002388]`	`['stable', 0.997415]`	Failed
case 5	`['unstable', 19.158456]`	`['unstable', 9.896296]`	Failed
case 6	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 7	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed

SHA-256 / b309332c62cdeb96d263ec73559eda598bb19ad9b6dc16db4856fdcfe10f9fcc

2 / The unsuccessful fix

Exit 1

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(k, m, c, dt):
    a = c / m
    b = k / m
    disc = a * a - 4 * b
    if disc >= 0:
        r = math.sqrt(disc)
        lams = [complex((-a + r) / 2, 0), complex((-a - r) / 2, 0)]
    else:
        r = math.sqrt(-disc)
        lams = [complex(-a, r / 2), complex(-a, -r / 2)]
    amp = max(abs(1 + dt * lam) for lam in lams)
    if abs(amp - 1) <= 1e-9:
        label = 'marginal'
    elif amp < 1:
        label = 'stable'
    else:
        label = 'unstable'
    return [label, round(amp, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[33.665, 3.374, 3.698, 0.05], ['stable', 0.984958]], [[17.863, 3.409, 0.738, 0.4], ['unstable', 1.323555]], [[15.545, 2.583, 6.586, 1.0], ['unstable', 2.11387]], [[46.563, 1.397, 1.187, 0.01], ['stable', 0.997415]], [[39.672, 0.379, 2.933, 1.0], ['unstable', 9.896296]], [[0.0, 3.87, 4.121, 1.0], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]]], [[[23.283, 2.64, 7.844, 0.1], ['stable', 0.889422]], [[24.43, 0.273, 3.575, 0.1], ['stable', 0.76508]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[8.452, 3.883, 5.856, 0.4], ['stable', 0.863147]], [[47.405, 3.962, 4.696, 0.05], ['stable', 0.985215]], [[0.0, 3.884, 4.78, 0.2], ['marginal', 1.0]], [[56.015, 2.524, 6.052, 1.0], ['unstable', 4.560172]]], [[[15.636, 2.616, 6.979, 0.8], ['unstable', 1.640448]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[8.52, 3.836, 5.13, 0.4], ['stable', 0.90578]], [[29.423, 0.607, 3.273, 0.1], ['stable', 0.972378]], [[0.0, 3.996, 2.36, 1.0], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[18.776, 1.514, 0.866, 0.01], ['stable', 0.997758]], [[17.34, 2.005, 4.299, 0.01], ['stable', 0.989658]]], [[[9.079, 3.699, 5.728, 0.05], ['stable', 0.963696]], [[15.373, 0.488, 0.687, 0.8], ['unstable', 4.476057]], [[0.0, 0.88, 2.713, 0.01], ['marginal', 1.0]], [[29.305, 3.778, 6.728, 0.2], ['stable', 0.976782]], [[41.131, 0.314, 6.187, 0.4], ['unstable', 3.751925]], [[3.996, 2.055, 0.915, 0.05], ['stable', 0.991261]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]]], [[[49.306, 2.663, 2.331, 0.2], ['unstable', 1.251216]], [[58.722, 3.822, 5.191, 0.8], ['unstable', 3.121945]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.689, 5.168, 0.01], ['marginal', 1.0]], [[39.004, 1.806, 2.569, 0.01], ['stable', 0.993949]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[36.246, 3.029, 2.637, 0.2], ['unstable', 1.142163]], [[13.954, 1.677, 2.862, 1.0], ['unstable', 2.759383]]]]
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	`['unstable', 1.0008]`	`['unstable', 1.0008]`	Passed
case 1	`['stable', 0.957911]`	`['stable', 0.984958]`	Failed
case 2	`['unstable', 1.292605]`	`['unstable', 1.323555]`	Failed
case 3	`['unstable', 2.606648]`	`['unstable', 2.11387]`	Failed
case 4	`['stable', 0.993174]`	`['stable', 0.997415]`	Failed
case 5	`['unstable', 11.623876]`	`['unstable', 9.896296]`	Failed
case 6	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 7	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed

SHA-256 / 0073a40991c281594fb0dd2dd0575c2b656c34f2c260fdb824d3b71f98a6f5d3

3 / The verified repair

Exit 0

"""Failure Map reference implementation. Python standard library only."""
import json
import math
N = 1
observations = []
def solve(k, m, c, dt):
    a = c / m
    b = k / m
    disc = a * a - 4 * b
    if disc >= 0:
        r = math.sqrt(disc)
        lams = [complex((-a + r) / 2, 0), complex((-a - r) / 2, 0)]
    else:
        r = math.sqrt(-disc)
        lams = [complex(-a / 2, r / 2), complex(-a / 2, -r / 2)]
    amp = max(abs(1 + dt * lam) for lam in lams)
    if abs(amp - 1) <= 1e-9:
        label = 'marginal'
    elif amp < 1:
        label = 'stable'
    else:
        label = 'unstable'
    return [label, round(amp, 6)]
def check(label, actual, expected):
    observations.append({"check": label, "actual": actual, "expected": expected, "passed": actual == expected})
fixtures = [[[[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[33.665, 3.374, 3.698, 0.05], ['stable', 0.984958]], [[17.863, 3.409, 0.738, 0.4], ['unstable', 1.323555]], [[15.545, 2.583, 6.586, 1.0], ['unstable', 2.11387]], [[46.563, 1.397, 1.187, 0.01], ['stable', 0.997415]], [[39.672, 0.379, 2.933, 1.0], ['unstable', 9.896296]], [[0.0, 3.87, 4.121, 1.0], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]]], [[[23.283, 2.64, 7.844, 0.1], ['stable', 0.889422]], [[24.43, 0.273, 3.575, 0.1], ['stable', 0.76508]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[8.452, 3.883, 5.856, 0.4], ['stable', 0.863147]], [[47.405, 3.962, 4.696, 0.05], ['stable', 0.985215]], [[0.0, 3.884, 4.78, 0.2], ['marginal', 1.0]], [[56.015, 2.524, 6.052, 1.0], ['unstable', 4.560172]]], [[[15.636, 2.616, 6.979, 0.8], ['unstable', 1.640448]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[8.52, 3.836, 5.13, 0.4], ['stable', 0.90578]], [[29.423, 0.607, 3.273, 0.1], ['stable', 0.972378]], [[0.0, 3.996, 2.36, 1.0], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[18.776, 1.514, 0.866, 0.01], ['stable', 0.997758]], [[17.34, 2.005, 4.299, 0.01], ['stable', 0.989658]]], [[[9.079, 3.699, 5.728, 0.05], ['stable', 0.963696]], [[15.373, 0.488, 0.687, 0.8], ['unstable', 4.476057]], [[0.0, 0.88, 2.713, 0.01], ['marginal', 1.0]], [[29.305, 3.778, 6.728, 0.2], ['stable', 0.976782]], [[41.131, 0.314, 6.187, 0.4], ['unstable', 3.751925]], [[3.996, 2.055, 0.915, 0.05], ['stable', 0.991261]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]]], [[[49.306, 2.663, 2.331, 0.2], ['unstable', 1.251216]], [[58.722, 3.822, 5.191, 0.8], ['unstable', 3.121945]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.689, 5.168, 0.01], ['marginal', 1.0]], [[39.004, 1.806, 2.569, 0.01], ['stable', 0.993949]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[36.246, 3.029, 2.637, 0.2], ['unstable', 1.142163]], [[13.954, 1.677, 2.862, 1.0], ['unstable', 2.759383]]]]
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	`['unstable', 1.0008]`	`['unstable', 1.0008]`	Passed
case 1	`['stable', 0.984958]`	`['stable', 0.984958]`	Passed
case 2	`['unstable', 1.323555]`	`['unstable', 1.323555]`	Passed
case 3	`['unstable', 2.11387]`	`['unstable', 2.11387]`	Passed
case 4	`['stable', 0.997415]`	`['stable', 0.997415]`	Passed
case 5	`['unstable', 9.896296]`	`['unstable', 9.896296]`	Passed
case 6	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 7	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed

SHA-256 / 16159c1572373a0abc8fa7e3d3d51cc7b9f8f5968b2febd6d9f9585c3988e93e

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

Case digest / 1c327c5a94dd07385e7f27bd0867dc67519900e10d8c74f5ed848dd4b065ef73