1import json, math
  2from pathlib import Path
  3import numpy as np
  4from scipy.linalg import eigvals, expm
  5
  6SEED = 1418
  7np.random.seed(SEED)
  8
  9
 10def spectral_abscissa(M):
 11    return float(np.max(np.real(eigvals(M))))
 12
 13
 14def linear_sweep():
 15    # x' = A x, y = [1,0]x; z' = A z + G(y-Cz)
 16    # The kernel-projected H is always [a22], independent of G.
 17    A = np.array([[0., 1.], [1., .5]])
 18    C = np.array([[1., 0.]])
 19    g1 = 2.0
 20    gs = np.linspace(0.0, 4.0, 81)
 21    rows = []
 22    for g2 in gs:
 23        G = np.array([[g1], [g2]])
 24        M = A - G @ C
 25        # U_K = [0,1]^T, hence H = .5.
 26        H = np.array([[0.5]])
 27        rows.append((float(g2), spectral_abscissa(M), float(H[0, 0])))
 28    stable = [g for g, ab, _ in rows if ab < 0]
 29    observed_boundary = min(stable) if stable else None
 30    # exact characteristic polynomial gives det(M)=g2-2 and trace=-1.5.
 31    predicted_boundary = 2.0
 32
 33    # Check exponential error decay at a stable and unstable gain with exact flow.
 34    x0 = np.array([1., -1.])
 35    times = np.linspace(0, 8, 161)
 36    checks = {}
 37    for g2 in (1.5, 3.0):
 38        M = A - np.array([[g1], [g2]]) @ C
 39        norms = np.array([np.linalg.norm(expm(M*t) @ x0) for t in times])
 40        # fit after transient, avoiding zero/roundoff
 41        slope = float(np.polyfit(times[20:], np.log(norms[20:] + 1e-30), 1)[0])
 42        checks[str(g2)] = {"spectral_abscissa": spectral_abscissa(M), "log_norm_slope": slope,
 43                           "norm_t0": float(norms[0]), "norm_t8": float(norms[-1])}
 44    return {"predicted_stability_boundary_g2": predicted_boundary,
 45            "observed_first_grid_stable_g2": observed_boundary,
 46            "grid_step": float(gs[1]-gs[0]),
 47            "projected_H": 0.5,
 48            "checks": checks}
 49
 50
 51def gain_sweep():
 52    A = np.array([[0., 1.], [1., .5]])
 53    C = np.array([[1., 0.]])
 54    U = np.array([[0.], [1.]])
 55    rows = []
 56    for g2 in np.linspace(0., 4., 17):
 57        G = np.array([[2.], [g2]])
 58        M = A - G @ C
 59        H = U.T @ M @ U
 60        rows.append({"g2": float(g2), "full_abscissa": spectral_abscissa(M),
 61                     "projected_H": float(H[0, 0])})
 62    return rows
 63
 64
 65def euler_boundary_sweep():
 66    # Prediction: explicit Euler is stable iff max |1+dt*lambda_i|<1.
 67    A = np.array([[0., 1.], [1., .5]])
 68    C = np.array([[1., 0.]])
 69    g1 = 2.0
 70    g2 = 3.0
 71    M = A - np.array([[g1], [g2]]) @ C
 72    dts = np.linspace(.001, 3.5, 3500)
 73    stable = []
 74    for dt in dts:
 75        rho = max(abs(eigvals(np.eye(2) + dt*M)))
 76        if rho < 1: stable.append(dt)
 77    # solve numerically at fine grid; report observed and direct formula check.
 78    observed = max(stable)
 79    eig = eigvals(M)
 80    # Exact disk condition for each (possibly complex) eigenvalue:
 81    # |1 + dt*lambda| < 1 => dt < -2 Re(lambda)/|lambda|^2.
 82    predicted = min(float(-2*l.real/(abs(l)**2)) for l in eig)
 83    return {"g2": g2, "eigenvalues": [[float(l.real), float(l.imag)] for l in eig],
 84            "predicted_dt_max": predicted, "observed_dt_max_grid": observed,
 85            "dt_grid": float(dts[1]-dts[0])}
 86
 87
 88def lorenz_rhs(x):
 89    sig, rho, beta = 10., 28., 8./3.
 90    return np.array([sig*(x[1]-x[0]), x[0]*(rho-x[2])-x[1], x[0]*x[1]-beta*x[2]])
 91
 92
 93def rk4_step(fun, x, dt):
 94    k1=fun(x); k2=fun(x+dt*k1/2); k3=fun(x+dt*k2/2); k4=fun(x+dt*k3)
 95    return x + dt*(k1+2*k2+2*k3+k4)/6
 96
 97
 98def lorenz_observer():
 99    # Known vector field, partial continuous observations x,y. This is a practical
100    # nudging check, not a trained neural ODE benchmark.
101    dt=.005; T=10.; n=int(T/dt)
102    x=np.array([1.,1.,1.]); z0=np.array([-8.,7.,18.])
103    C=np.array([[1.,0.,0.],[0.,1.,0.]])
104    gains={"open_loop":np.zeros((3,2)), "nudging":np.array([[8.,0.],[0.,8.],[12.,-4.]])}
105    out={}
106    for name,G in gains.items():
107        truth=x.copy(); z=z0.copy(); errors=[]
108        for _ in range(n):
109            y=C@truth
110            def obs(q): return lorenz_rhs(q)+G@(y-C@q)
111            truth=rk4_step(lorenz_rhs,truth,dt)
112            z=rk4_step(obs,z,dt)
113            errors.append(np.linalg.norm(truth-z))
114        out[name]={"initial_error":float(np.linalg.norm(x-z0)),
115                   "final_error":float(errors[-1]),
116                   "median_last_20pct":float(np.median(errors[int(.8*n):]))}
117    return out
118
119
120def main():
121    result={"seed":SEED, "linear_contraction_sweep":linear_sweep(),
122            "gain_sweep":gain_sweep(),
123            "euler_stability_sweep":euler_boundary_sweep(),
124            "lorenz_partial_observation":lorenz_observer(),
125            "interpretation":{
126              "prediction_1":"Full 2D observer stability begins at g2=2 from det(A-GC)>0; observed grid boundary should agree.",
127              "prediction_2":"At stable g2=3, Euler stability has a finite dt boundary predicted by -2 Re(lambda)/|lambda|^2; observed sweep should agree.",
128              "prediction_gain":"Projected H remains +0.5 across all gains, while full abscissa crosses zero at g2=2; this falsifies H-only sufficiency here.",
129              "prediction_3":"The proposed H on ker(C) remains +0.5 for every gain, despite full observer becoming stable; therefore the stated projected certificate is not sufficient in this coupled example."}}
130    Path('results.json').write_text(json.dumps(result,indent=2))
131    print(json.dumps(result,indent=2))
132
133if __name__=='__main__': main()