1import json
  2import math
  3import numpy as np
  4
  5SEED = 2771
  6UC = 2.0 / (3.0 * math.sqrt(3.0))
  7
  8def F(x, u):
  9    return x - x**3 + u
 10
 11def jac(x):
 12    return 1.0 - 3.0*x*x
 13
 14def equilibria(u):
 15    roots = np.roots([1.0, 0.0, -1.0, -u])
 16    return sorted(float(r.real) for r in roots if abs(r.imag) < 1e-8)
 17
 18def integrate_many(x0, u, T=35.0, damping=1.0):
 19    """Vectorized fixed-step RK4 integration of the scalar flow."""
 20    x = np.asarray(x0, dtype=float).copy()
 21    n = max(1, int(math.ceil(T / 0.04)))
 22    dt = T / n
 23    for _ in range(n):
 24        k1 = damping * F(x, u)
 25        k2 = damping * F(x + 0.5*dt*k1, u)
 26        k3 = damping * F(x + 0.5*dt*k2, u)
 27        k4 = damping * F(x + dt*k3, u)
 28        x += dt*(k1 + 2*k2 + 2*k3 + k4)/6.0
 29    return x
 30
 31def integrate(x0, u, T=35.0, damping=1.0):
 32    return float(integrate_many(np.array([x0]), u, T, damping)[0])
 33
 34def desired_eq(u):
 35    return max(equilibria(u))
 36
 37def estimate_saddle_node():
 38    us = np.linspace(0.0, 0.45, 9001)
 39    counts = np.array([len(equilibria(u)) for u in us])
 40    return float(us[np.where(counts >= 3)[0][-1]])
 41
 42def main():
 43    numeric_uc = estimate_saddle_node()
 44    u_vals = [0.0, 0.10, 0.20, 0.35, 0.38]
 45    stability = []
 46    for u in u_vals:
 47        roots = equilibria(u)
 48        xp = max(roots)
 49        stability.append({"u": u, "roots": roots, "desired_jacobian": jac(xp),
 50                          "locally_stable": jac(xp) < 0, "coexisting": len(roots) == 3})
 51
 52    # At u=0, x_ref=1 and separatrix=0. Thus recovery is exactly the event
 53    # 1 + eps*z > 0 (up to numerical integration tolerance).
 54    eps_values = [0.10, 0.30, 0.50, 0.70, 0.90, 1.10]
 55    basin = []
 56    for eps in eps_values:
 57        rng = np.random.default_rng(SEED + int(eps * 1000))
 58        z = rng.normal(size=4000)
 59        xref = desired_eq(0.0)
 60        final = integrate_many(xref + eps*z, 0.0)
 61        basin.append({"epsilon": eps,
 62                      "B_hat": float(np.mean(np.abs(final-xref) < .08)),
 63                      "B_predicted_gaussian": float(0.5*(1+math.erf(1/(eps*math.sqrt(2))))),
 64                      "separatrix": 0.0, "reference": xref})
 65
 66    # Finite-horizon comparison. The intervention is deliberately tested after
 67    # local stability is restored. The guard estimates basin recovery with probes
 68    # and retains damping when recovery is not demonstrated.
 69    def trial(mode, n=300):
 70        rng = np.random.default_rng(SEED + (1 if mode == 'guard' else 2))
 71        xref = desired_eq(0.0)
 72        failures = releases = 0
 73        for _ in range(n):
 74            x = xref + 0.95 * rng.normal()
 75            if mode == 'threshold':
 76                c = 1.0  # local Jacobian is safely negative, so threshold-only releases
 77                releases += 1
 78            else:
 79                probes = xref + 0.95*rng.normal(size=80)
 80                probe_final = integrate_many(probes, 0.0, T=6.0, damping=0.35)
 81                B = np.mean(np.abs(probe_final-xref) < .15)
 82                c = 1.0 if jac(xref) <= -0.5 and B >= .95 else 0.35
 83                releases += int(c == 1.0)
 84            y = integrate(x, 0.0, T=6.0, damping=c)
 85            failures += int(abs(y-xref) >= .15)
 86        return {"failure_rate": failures/n, "release_rate": releases/n}
 87
 88    result = {
 89        "seed": SEED,
 90        "model": "dx/dt = x - x^3 + u (desired positive equilibrium)",
 91        "predictions": {
 92            "saddle_node": {"analytic_uc": UC, "numeric_last_three_root_u": numeric_uc,
 93                            "absolute_error": abs(UC-numeric_uc)},
 94            "local_vs_global": "desired Jacobian remains negative throughout the three-root bistable interval; local stability does not imply global recovery",
 95            "basin_transition": "at u=0, separatrix is x=0 and desired equilibrium is x=1; recovery probability is Phi(1/epsilon)"
 96        },
 97        "stability_sweep": stability,
 98        "basin_sweep": basin,
 99        "finite_horizon_comparison": {"threshold_only": trial("threshold"), "basin_guard": trial("guard")}
100    }
101    with open("results.json", "w") as f:
102        json.dump(result, f, indent=2)
103    print(json.dumps(result, indent=2))
104
105if __name__ == '__main__':
106    main()