1import json
 2from pathlib import Path
 3import numpy as np
 4from scipy.linalg import expm, eigvals
 5
 6OUT = Path('floquet_results.json')
 7M1 = np.array([[-5.83182011, 6.67253161], [4.98376025, -2.86979823]])
 8M2 = np.array([[4.15640208, -3.67581620], [5.62849205, -1.08915620]])
 9I = np.eye(2)
10
11def edge(A):
12    return float(np.max(np.real(eigvals(A))))
13
14def gamma(M1, M2, damping, period):
15    h = period / 2.0
16    P = expm((-damping*I + M2)*h) @ expm((-damping*I + M1)*h)
17    return float(np.log(np.max(np.abs(eigvals(P)))) / period)
18
19def avg_edge(d):
20    return .5*(edge(M1-d*I) + edge(M2-d*I))
21
22def crossing(xs, ys):
23    for x0,x1,y0,y1 in zip(xs[:-1],xs[1:],ys[:-1],ys[1:]):
24        if y0*y1 <= 0:
25            return float(x0-y0*(x1-x0)/(y1-y0))
26    return None
27
28def main():
29    e1,e2=edge(M1),edge(M2)
30    comm=float(np.linalg.norm(M2@M1-M1@M2))
31    periods=[.01,.03,.1,.2,.5,1.,2.,5.,10.,20.,50.,100.]
32    # P1: scalar damping shifts the exact Floquet exponent by -d.
33    shifts=[]
34    for T in periods:
35        g0,g2=gamma(M1,M2,0,T),gamma(M1,M2,2,T)
36        shifts.append({'period':T,'measured_shift':g0-g2,'predicted_shift':2.,
37                       'abs_error':abs(g0-g2-2.)})
38    # P2: at long periods the boundary approaches the mean instantaneous edge.
39    ds=np.linspace(0,3,1201)
40    adiabatic=.5*(e1+e2)
41    boundaries=[]
42    for T in periods:
43        ys=np.array([gamma(M1,M2,d,T) for d in ds])
44        boundaries.append({'period':T,'floquet_boundary':crossing(ds,ys),
45                           'adiabatic_boundary':adiabatic})
46    # P3: noncommuting finite-period correction can stabilize even when the
47    # phasewise rightmost eigenvalues and constant-mean system are unstable.
48    mean_edge=edge((M1+M2)/2)
49    correction=[]
50    for T in periods:
51        g=gamma(M1,M2,0,T)
52        correction.append({'period':T,'frequency':2*np.pi/T,'floquet_gamma':g,
53                           'constant_mean_gamma':mean_edge,'correction':g-mean_edge})
54    d=1.4; T=1.0
55    rescue={'damping':d,'period':T,'phase_edges':[e1-d,e2-d],
56            'floquet_gamma':gamma(M1,M2,d,T),'adiabatic_estimate':avg_edge(d),
57            'constant_mean_baseline':mean_edge-d}
58    rescue['instantaneous_edges_positive']=min(rescue['phase_edges'])>0
59    rescue['floquet_stable']=rescue['floquet_gamma']<0
60    rescue['baseline_unstable']=rescue['constant_mean_baseline']>0
61    # Linearized-training proxy: norm after repeated periods, same mean B and d.
62    x0=np.array([1.,-1.]); n=20
63    h=T/2
64    P=expm((-d*I+M2)*h)@expm((-d*I+M1)*h)
65    periodic=np.linalg.matrix_power(P,n)@x0
66    baseline=expm((-d*I+(M1+M2)/2)*T*n)@x0
67    proxy={'periods':n,'initial_norm':float(np.linalg.norm(x0)),
68           'constant_mean_final_norm':float(np.linalg.norm(baseline)),
69           'periodic_final_norm':float(np.linalg.norm(periodic))}
70    result={'matrices':{'M1':M1.tolist(),'M2':M2.tolist(),'edge_M1':e1,
71                         'edge_M2':e2,'commutator_norm':comm},
72      'predictions':{'P1':'gamma(d)=gamma(0)-d exactly',
73       'P2':'slow-period boundary approaches average instantaneous edge',
74       'P3':'noncommuting finite-period forcing can stabilize beyond mean-matrix dynamics'},
75      'damping_shift':shifts,'boundaries':boundaries,'mean_matrix_edge':mean_edge,
76      'correction':correction,'rescue':rescue,'proxy_comparison':proxy}
77    OUT.write_text(json.dumps(result,indent=2)); print(json.dumps(result,indent=2))
78if __name__=='__main__': main()