1import json, math, os
 2import numpy as np
 3
 4# Bifurcation-aware adaptive RK4 toy experiment.
 5# The scalar normal form is dz/dt = mu + z^2, with passage measured from z=-1 to z=+1.
 6
 7def rhs(z, mu):
 8    return mu + z*z
 9
10def rk4_step(z, h, mu):
11    k1 = rhs(z, mu)
12    k2 = rhs(z + .5*h*k1, mu)
13    k3 = rhs(z + .5*h*k2, mu)
14    k4 = rhs(z + h*k3, mu)
15    return z + h*(k1 + 2*k2 + 2*k3 + k4)/6
16
17def exact_residence(mu):
18    # Integral dz/(mu+z^2), for positive mu, between -1 and +1.
19    q = math.sqrt(mu)
20    return 2.0*math.atan(1.0/q)/q
21
22def controller_h(mu, hmin=5e-4, hmax=5e-2, mu0=1e-2, delta=1e-12):
23    # Formula in the proposal, with mu_hat clipped to its positive ghost regime.
24    return min(hmax, max(hmin, hmax*math.sqrt((max(mu, 0.0)+delta)/mu0)))
25
26def passage(mu, mode, fixed_h=2e-3):
27    z, t, n = -1.0, 0.0, 0
28    while z < 1.0 and n < 100000000:
29        h = fixed_h if mode == 'fixed' else controller_h(mu)
30        # Do not step far beyond the event; this also makes errors comparable.
31        h = min(h, (1.0-z)/max(rhs(z, mu), 1e-30))
32        z = rk4_step(z, h, mu)
33        t += h; n += 1
34    return t, n, z
35
36def log_slope(x, y):
37    return float(np.polyfit(np.log(x), np.log(y), 1)[0])
38
39def slow_passage(eps):
40    # mu=eps*t. Start at mu=-0.1 on the attracting quasi-static branch,
41    # then measure time from the fold (t=0) until z crosses zero.
42    t = -0.1/eps
43    z = -math.sqrt(.1)
44    # A modest step relative to the universal eps^(1/3) time scale.
45    h = min(0.02, 0.03*eps**(-1/3))
46    n = 0
47    while z < 0.0 and n < 20000000:
48        # RK4 with time-varying mu.
49        def f(tt, zz): return eps*tt + zz*zz
50        k1=f(t,z); k2=f(t+h/2,z+h*k1/2); k3=f(t+h/2,z+h*k2/2); k4=f(t+h,z+h*k3)
51        zn=z+h*(k1+2*k2+2*k3+k4)/6
52        if zn >= 0.0:
53            # linear interpolation of the crossing, sufficient for scaling.
54            frac=(0.0-z)/(zn-z)
55            t=t+frac*h
56            break
57        z=zn; t += h; n += 1
58    return t, eps*t, n
59
60def main():
61    np.random.seed(0)
62    mus=np.logspace(-4, -0.5, 8)
63    residence=[]; fixed=[]; adaptive=[]
64    for mu in mus:
65        exact=exact_residence(mu)
66        tf,nf,_=passage(mu,'fixed')
67        ta,na,_=passage(mu,'adaptive')
68        residence.append(exact)
69        fixed.append({'time':tf,'n':nf,'relerr':abs(tf-exact)/exact})
70        adaptive.append({'time':ta,'n':na,'relerr':abs(ta-exact)/exact,'h':controller_h(mu)})
71    # Prediction 1: ghost residence has slope -1/2 at small positive mu.
72    small=mus[:5]
73    residence_slope=log_slope(small, np.array(residence[:5]))
74    # Prediction 2: controller h ~ sqrt(mu), hence controller allocations over the
75    # same physical ghost crossing scale approximately mu^-1 (T~mu^-1/2 and 1/h~mu^-1/2).
76    alloc_slope=log_slope(mus[:5], np.array([x['n'] for x in adaptive[:5]],float))
77    # Prediction 3: slow saddle-node passage time scales eps^-1/3.
78    eps=np.logspace(-4,-1,7)
79    slow=[slow_passage(e) for e in eps]
80    slow_times=np.array([x[0] for x in slow])
81    slow_slope=log_slope(eps, slow_times)
82    out={
83      'predictions':{
84        'ghost_residence_slope_pred':-0.5,'ghost_residence_slope_obs':residence_slope,
85        'controller_allocation_slope_pred':-1.0,'controller_allocation_slope_obs':alloc_slope,
86        'slow_delay_slope_pred':-1/3,'slow_delay_slope_obs':slow_slope},
87      'mu_sweep':[{'mu':float(m),'exact_T':float(residence[i]),'fixed':fixed[i], 'adaptive':adaptive[i],
88                   'fixed_over_adaptive_n':fixed[i]['n']/max(adaptive[i]['n'],1)} for i,m in enumerate(mus)],
89      'slow_sweep':[{'epsilon':float(e),'delay_time':float(slow[i][0]),'mu_at_crossing':float(slow[i][1]),'steps':slow[i][2]} for i,e in enumerate(eps)],
90      'summary':{
91        'mean_fixed_relerr':float(np.mean([x['relerr'] for x in fixed])),
92        'mean_adaptive_relerr':float(np.mean([x['relerr'] for x in adaptive])),
93        'far_mu_eval_saving':float(1-adaptive[-1]['n']/fixed[-1]['n']),
94        'near_mu_eval_saving':float(1-adaptive[0]['n']/fixed[0]['n'])}}
95    with open('results.json','w') as f: json.dump(out,f,indent=2)
96    print(json.dumps(out,indent=2))
97
98if __name__=='__main__': main()