1import json
 2import numpy as np
 3
 4
 5def transverse_matrix(a=1.0, lambda_p=2.0, lambda_i=2.0, k_p=0.0, k_i=0.0):
 6    """Scalar n=1 instance of the paper's transverse matrix."""
 7    return np.array([[-a - k_p * lambda_p, -1.0],
 8                     [k_i * lambda_i, 0.0]], dtype=float)
 9
10
11def euler_radius(M, h):
12    return float(np.max(np.abs(np.linalg.eigvals(np.eye(M.shape[0]) + h * M))))
13
14
15def simulate(a=1.0, bias_difference=1.0, k_p=0.0, k_i=0.0,
16             h=0.01, steps=5000):
17    """Explicit Euler difference dynamics for two copies and L transverse mode 2."""
18    d, z = 0.0, 0.0
19    ds = np.empty(steps + 1)
20    ds[0] = d
21    for t in range(steps):
22        dd = -a * d + bias_difference - 2.0 * k_p * d - z
23        dz = 2.0 * k_i * d
24        d, z = d + h * dd, z + h * dz
25        ds[t + 1] = d
26    return ds
27
28
29def eig_pairs(M):
30    return [[float(x.real), float(x.imag)] for x in np.linalg.eigvals(M)]
31
32
33def main():
34    a, h = 1.0, 0.01
35    cases = {
36        "uncoupled": (0.0, 0.0),
37        "proportional": (1.0, 0.0),
38        "PI_stable": (1.0, 0.2),
39        "PI_unstable": (1.0, 200.0),
40    }
41    results = {}
42    for name, (kp, ki) in cases.items():
43        d = simulate(k_p=kp, k_i=ki, h=h)
44        M = transverse_matrix(k_p=kp, k_i=ki)
45        results[name] = {
46            "kP": kp, "kI": ki,
47            "spectral_radius": euler_radius(M, h),
48            "max_continuous_real_part": float(np.max(np.linalg.eigvals(M).real)),
49            "continuous_eigenvalues": eig_pairs(M),
50            "final_disagreement": float(d[-1]),
51            "mean_abs_last_500": float(np.mean(np.abs(d[-500:]))),
52            "max_abs": float(np.max(np.abs(d))),
53        }
54
55    scan = []
56    for ki in np.linspace(0.0, 200.0, 201):
57        M = transverse_matrix(k_p=1.0, k_i=float(ki))
58        rho = euler_radius(M, h)
59        d = simulate(k_p=1.0, k_i=float(ki), h=h, steps=3000)
60        scan.append({"kI": float(ki), "rho": rho,
61                     "strictly_stable_by_rho": bool(rho < 1.0),
62                     "tail_abs": float(np.mean(np.abs(d[-300:])))})
63    positive = [x for x in scan if x["kI"] > 0]
64    stable = [x for x in positive if x["strictly_stable_by_rho"]]
65    unstable = [x for x in positive if not x["strictly_stable_by_rho"]]
66    results["scan"] = {
67        "largest_positive_stable_kI": max(x["kI"] for x in stable),
68        "first_positive_unstable_kI": min(x["kI"] for x in unstable),
69        "points": scan,
70    }
71
72    # Reproducible checks of both parts of the claimed signature.
73    assert results["proportional"]["mean_abs_last_500"] > 0.3
74    assert results["PI_stable"]["mean_abs_last_500"] < 0.001
75    assert results["PI_stable"]["spectral_radius"] < 1.0
76    assert results["PI_unstable"]["spectral_radius"] > 1.0
77    assert results["PI_unstable"]["max_abs"] > 1e6
78
79    with open("pi_results.json", "w") as f:
80        json.dump(results, f, indent=2)
81    summary = {k: v for k, v in results.items() if k != "scan"}
82    print(json.dumps(summary, indent=2))
83    print("positive-kI Euler boundary:",
84          results["scan"]["largest_positive_stable_kI"],
85          "to", results["scan"]["first_positive_unstable_kI"])
86
87
88if __name__ == "__main__":
89    main()