FA-86921 / Physics integrator stability / Open access

Complex eigenvalues are tested with the real-axis step criterion · case 01

Undamped oscillators are called stable although explicit Euler always gains energy on them.

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

ROOT CAUSE

Amplification is replaced by |dt*lambda|/2, the interval test valid only for negative real eigenvalues.

VERIFIED REPAIR

Use the modulus of 1+dt*lambda in the complex plane.

Unsuccessful approach: Dropping the imaginary part of dt*lambda hides the oscillatory growth.

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 / 2), complex(-a / 2, -r / 2)]
    amp = max(abs(dt * lam) / 2 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, 1.0, 0.5], ['marginal', 1.0]], [[30.639, 3.582, 3.661, 0.4], ['unstable', 1.399912]], [[32.4, 3.038, 3.326, 0.2], ['unstable', 1.098925]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[11.082, 0.828, 5.406, 0.8], ['unstable', 2.083893]], [[38.284, 3.135, 4.579, 0.2], ['unstable', 1.093778]], [[42.221, 2.418, 7.819, 0.01], ['stable', 0.984586]]], [[[50.892, 1.977, 6.092, 0.8], ['unstable', 3.874242]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 2.92, 3.001, 0.1], ['marginal', 1.0]], [[2.124, 2.111, 2.83, 0.8], ['stable', 0.755952]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[1.0, 1.0, 3.0, 0.1], ['stable', 0.961803]], [[8.452, 3.883, 5.856, 0.4], ['stable', 0.863147]]], [[[7.364, 1.8, 2.173, 0.1], ['stable', 0.959265]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[14.137, 0.782, 5.45, 1.0], ['unstable', 3.479755]], [[33.332, 2.078, 0.851, 0.05], ['unstable', 1.009765]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[34.941, 2.954, 0.853, 0.2], ['unstable', 1.189699]], [[26.203, 2.441, 7.403, 0.4], ['unstable', 1.226546]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]]], [[[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 3.727, 0.555, 0.2], ['marginal', 1.0]], [[59.202, 0.85, 4.717, 0.1], ['unstable', 1.068435]], [[25.042, 3.597, 6.989, 0.1], ['stable', 0.935584]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[32.965, 3.721, 3.681, 0.1], ['stable', 0.99482]], [[28.743, 3.175, 4.507, 0.8], ['unstable', 2.378706]]], [[[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[36.647, 2.751, 2.545, 0.01], ['stable', 0.996033]], [[22.607, 3.348, 6.875, 0.2], ['stable', 0.92704]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[14.174, 2.424, 0.649, 0.8], ['unstable', 2.127938]], [[1.0, 1.0, 3.0, 0.1], ['stable', 0.961803]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[29.176, 2.254, 7.792, 0.01], ['stable', 0.983222]]]]
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	`['stable', 0.25]`	`['marginal', 1.0]`	Failed
case 1	`['stable', 0.584931]`	`['unstable', 1.399912]`	Failed
case 2	`['stable', 0.326572]`	`['unstable', 1.098925]`	Failed
case 3	`['stable', 0.025]`	`['marginal', 1.0]`	Failed
case 4	`['unstable', 1.0004]`	`['unstable', 1.0008]`	Failed
case 5	`['unstable', 1.463369]`	`['unstable', 2.083893]`	Failed
case 6	`['stable', 0.349454]`	`['unstable', 1.093778]`	Failed
case 7	`['stable', 0.020893]`	`['stable', 0.984586]`	Failed

SHA-256 / 740f11e96eaf2ad83da7d4ecd3532ee5b24490ba2b585e86080fe3618d1a3341

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 / 2, r / 2), complex(-a / 2, -r / 2)]
    amp = max(abs(1 + dt * lam.real) 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, 1.0, 0.5], ['marginal', 1.0]], [[30.639, 3.582, 3.661, 0.4], ['unstable', 1.399912]], [[32.4, 3.038, 3.326, 0.2], ['unstable', 1.098925]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[11.082, 0.828, 5.406, 0.8], ['unstable', 2.083893]], [[38.284, 3.135, 4.579, 0.2], ['unstable', 1.093778]], [[42.221, 2.418, 7.819, 0.01], ['stable', 0.984586]]], [[[50.892, 1.977, 6.092, 0.8], ['unstable', 3.874242]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 2.92, 3.001, 0.1], ['marginal', 1.0]], [[2.124, 2.111, 2.83, 0.8], ['stable', 0.755952]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[1.0, 1.0, 3.0, 0.1], ['stable', 0.961803]], [[8.452, 3.883, 5.856, 0.4], ['stable', 0.863147]]], [[[7.364, 1.8, 2.173, 0.1], ['stable', 0.959265]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[14.137, 0.782, 5.45, 1.0], ['unstable', 3.479755]], [[33.332, 2.078, 0.851, 0.05], ['unstable', 1.009765]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[34.941, 2.954, 0.853, 0.2], ['unstable', 1.189699]], [[26.203, 2.441, 7.403, 0.4], ['unstable', 1.226546]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]]], [[[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 3.727, 0.555, 0.2], ['marginal', 1.0]], [[59.202, 0.85, 4.717, 0.1], ['unstable', 1.068435]], [[25.042, 3.597, 6.989, 0.1], ['stable', 0.935584]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[32.965, 3.721, 3.681, 0.1], ['stable', 0.99482]], [[28.743, 3.175, 4.507, 0.8], ['unstable', 2.378706]]], [[[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[36.647, 2.751, 2.545, 0.01], ['stable', 0.996033]], [[22.607, 3.348, 6.875, 0.2], ['stable', 0.92704]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[14.174, 2.424, 0.649, 0.8], ['unstable', 2.127938]], [[1.0, 1.0, 3.0, 0.1], ['stable', 0.961803]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[29.176, 2.254, 7.792, 0.01], ['stable', 0.983222]]]]
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	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 1	`['stable', 0.795589]`	`['unstable', 1.399912]`	Failed
case 2	`['stable', 0.89052]`	`['unstable', 1.098925]`	Failed
case 3	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 4	`['unstable', 1.0008]`	`['unstable', 1.0008]`	Passed
case 5	`['unstable', 1.611594]`	`['unstable', 2.083893]`	Failed
case 6	`['stable', 0.853939]`	`['unstable', 1.093778]`	Failed
case 7	`['stable', 0.983832]`	`['stable', 0.984586]`	Failed

SHA-256 / b403bd752602e6f2b5f903d169695e788dbaf440deea4e1fcdaa9b0bb20b202b

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, 1.0, 0.5], ['marginal', 1.0]], [[30.639, 3.582, 3.661, 0.4], ['unstable', 1.399912]], [[32.4, 3.038, 3.326, 0.2], ['unstable', 1.098925]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[11.082, 0.828, 5.406, 0.8], ['unstable', 2.083893]], [[38.284, 3.135, 4.579, 0.2], ['unstable', 1.093778]], [[42.221, 2.418, 7.819, 0.01], ['stable', 0.984586]]], [[[50.892, 1.977, 6.092, 0.8], ['unstable', 3.874242]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 2.92, 3.001, 0.1], ['marginal', 1.0]], [[2.124, 2.111, 2.83, 0.8], ['stable', 0.755952]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[1.0, 1.0, 3.0, 0.1], ['stable', 0.961803]], [[8.452, 3.883, 5.856, 0.4], ['stable', 0.863147]]], [[[7.364, 1.8, 2.173, 0.1], ['stable', 0.959265]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[14.137, 0.782, 5.45, 1.0], ['unstable', 3.479755]], [[33.332, 2.078, 0.851, 0.05], ['unstable', 1.009765]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[34.941, 2.954, 0.853, 0.2], ['unstable', 1.189699]], [[26.203, 2.441, 7.403, 0.4], ['unstable', 1.226546]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]]], [[[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[0.0, 3.727, 0.555, 0.2], ['marginal', 1.0]], [[59.202, 0.85, 4.717, 0.1], ['unstable', 1.068435]], [[25.042, 3.597, 6.989, 0.1], ['stable', 0.935584]], [[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[32.965, 3.721, 3.681, 0.1], ['stable', 0.99482]], [[28.743, 3.175, 4.507, 0.8], ['unstable', 2.378706]]], [[[0.0, 2.0, 1.0, 0.1], ['marginal', 1.0]], [[36.647, 2.751, 2.545, 0.01], ['stable', 0.996033]], [[22.607, 3.348, 6.875, 0.2], ['stable', 0.92704]], [[0.0, 1.0, 2.0008, 1.0], ['unstable', 1.0008]], [[14.174, 2.424, 0.649, 0.8], ['unstable', 2.127938]], [[1.0, 1.0, 3.0, 0.1], ['stable', 0.961803]], [[0.0, 1.0, 1.0, 0.5], ['marginal', 1.0]], [[29.176, 2.254, 7.792, 0.01], ['stable', 0.983222]]]]
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	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 1	`['unstable', 1.399912]`	`['unstable', 1.399912]`	Passed
case 2	`['unstable', 1.098925]`	`['unstable', 1.098925]`	Passed
case 3	`['marginal', 1.0]`	`['marginal', 1.0]`	Passed
case 4	`['unstable', 1.0008]`	`['unstable', 1.0008]`	Passed
case 5	`['unstable', 2.083893]`	`['unstable', 2.083893]`	Passed
case 6	`['unstable', 1.093778]`	`['unstable', 1.093778]`	Passed
case 7	`['stable', 0.984586]`	`['stable', 0.984586]`	Passed

SHA-256 / 8c490f8c23daa05a95eb44f7ade5b815bfc640cf474e9655a094b936861dddd6

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

Case digest / 5a480e81f8c030237cfbb3c538d21d6e4c5775bb43c648cebad9e0296d67f965