1import json, math
  2import numpy as np
  3
  4# Homogeneous cubic field: F_x=alpha*x^3, F_y=(alpha+beta)*x^2*y.
  5# At u*=e1: a=alpha, G=0, and the tangent Jacobian eigenvalue is beta.
  6def field(z, alpha, beta):
  7    x, y = z
  8    return np.array([alpha*x**3, (alpha+beta)*x*x*y], dtype=float)
  9
 10def leading_on_unit(u, alpha, beta):
 11    return field(u, alpha, beta)
 12
 13def controller(z, alpha, beta, rc=2.0, delta=.15, k=8.0, qtol=.20):
 14    r = np.linalg.norm(z)
 15    if r == 0: return field(z, alpha, beta), False
 16    u = z/r
 17    v = leading_on_unit(u, alpha, beta)
 18    a = float(u @ v)
 19    q = np.linalg.norm(v-a*u)
 20    # smooth radial damping gated by radius and the near-equilibrium test
 21    kap = k * np.logaddexp(0., (r-rc)/delta)
 22    active = (r > rc and q < qtol and a > 0)
 23    return field(z, alpha, beta) - (kap*u if active else 0.), active
 24
 25def rk4_step(z, dt, fun):
 26    k1=fun(z); k2=fun(z+dt*k1/2); k3=fun(z+dt*k2/2); k4=fun(z+dt*k3)
 27    return z + dt*(k1+2*k2+2*k3+k4)/6
 28
 29def integrate(z0, T, dt, fun, rcap=1e5):
 30    z=np.array(z0,dtype=float); maxr=np.linalg.norm(z); interventions=0; first=None
 31    n=int(math.ceil(T/dt)); actual=0.
 32    for i in range(n):
 33        h=min(dt,T-actual)
 34        if h<=0: break
 35        znew=rk4_step(z,h,fun)
 36        actual += h; z=znew
 37        r=np.linalg.norm(z); maxr=max(maxr,r)
 38        if getattr(fun,'last_active',False): interventions += 1; first=first if first is not None else actual
 39        if not np.all(np.isfinite(z)) or r>rcap:
 40            return {'z':z,'r':r,'maxr':maxr,'t':actual,'failed':True,'interventions':interventions,'first':first}
 41    return {'z':z,'r':np.linalg.norm(z),'maxr':maxr,'t':actual,'failed':False,'interventions':interventions,'first':first}
 42
 43def make_control_fun(alpha,beta):
 44    state={'active':False}
 45    def f(z):
 46        v, act=controller(z,alpha,beta)
 47        state['active']=act
 48        f.last_active=act
 49        return v
 50    f.last_active=False
 51    return f
 52
 53def math_check():
 54    # Exact infinity-atlas quantities at e1, checked by finite differences.
 55    alpha,beta=1.7,.63; u=np.array([1.,0.]); v=leading_on_unit(u,alpha,beta)
 56    a=float(u@v); eps=1e-6
 57    def G(th):
 58        uu=np.array([math.cos(th),math.sin(th)])
 59        vv=leading_on_unit(uu,alpha,beta)
 60        return vv-(uu@vv)*uu
 61    num=(G(eps)[1]-G(-eps)[1])/(2*eps)
 62    return {'a_at_u_star':a,'predicted_alpha':alpha,'tangent_eigenvalue_fd':num,'predicted_beta':beta,
 63            'relative_errors':[abs(a-alpha)/alpha,abs(num-beta)/beta]}
 64
 65def radial_scaling():
 66    # Along exact e1, blow-up time is t*=1/(2 alpha r0^2). Measure when r reaches R.
 67    r0=.8; R=20.; dt=5e-5
 68    rows=[]
 69    for alpha in [.5, .8, 1.2, 1.8]:
 70        target=(1/r0**2-1/R**2)/(2*alpha)
 71        T=target + .05
 72        out=integrate([r0,0],T,dt,lambda z:field(z,alpha,0),rcap=R)
 73        rows.append({'alpha':alpha,'predicted_time_to_R':target,'observed_time':out['t'],'relative_error':abs(out['t']-target)/target})
 74    return rows
 75
 76def angular_scaling():
 77    # For small angle, d log(theta)/d tau = beta. tau=int r^2 dt, and use a short interval.
 78    r0=.35; theta0=1e-4; T=.25; dt=1e-5
 79    rows=[]
 80    for beta in [-.8,.2,.5,1.0,1.5]:
 81        alpha=0.4
 82        z0=[r0*math.cos(theta0),r0*math.sin(theta0)]
 83        out=integrate(z0,T,dt,lambda z:field(z,alpha,beta))
 84        th0=math.atan2(z0[1],z0[0]); th=math.atan2(out['z'][1],out['z'][0])
 85        tau=(r0**2)*T # leading approximation, adequate at small T
 86        measured=math.log(abs(th/th0))/tau
 87        rows.append({'beta':beta,'predicted_slope':beta,'observed_slope':measured,'relative_error':abs(measured-beta)/max(abs(beta),1e-12)})
 88    return rows
 89
 90def controller_sweep():
 91    # Perturbed near-equilibrium direction: vanilla bursts, controller damps radial growth.
 92    rows=[]; r0=.8; T=.8; dt=2e-5; beta=.5; theta=1e-3
 93    for alpha in [.8,1.2,1.6,2.0]:
 94        z0=[r0*math.cos(theta),r0*math.sin(theta)]
 95        vanilla=integrate(z0,T,dt,lambda z:field(z,alpha,beta),rcap=1e4)
 96        cf=make_control_fun(alpha,beta)
 97        controlled=integrate(z0,T,dt,cf,rcap=1e4)
 98        rows.append({'alpha':alpha,'vanilla_max_r':vanilla['maxr'],'controlled_max_r':controlled['maxr'],
 99                     'vanilla_failed':vanilla['failed'],'controlled_failed':controlled['failed'],
100                     'controller_interventions':controlled['interventions'],'first_intervention':controlled['first']})
101    return rows
102
103def onset_boundary():
104    # Along e1, reaching rc by T requires alpha >= (1/r0^2-1/rc^2)/(2T).
105    r0=.8; rc=2.; T=.8; dt=5e-5
106    predicted=(1/r0**2-1/rc**2)/(2*T)
107    rows=[]
108    for alpha in [.6,.8,.9,1.0,1.2]:
109        out=integrate([r0,0],T,dt,lambda z:field(z,alpha,.5),rcap=1e4)
110        rows.append({'alpha':alpha,'reached_rc':bool(out['maxr']>rc),'max_r':out['maxr']})
111    return {'predicted_alpha_boundary':predicted,'sweep':rows}
112
113def main():
114    result={'math_check':math_check(),'radial_scaling':radial_scaling(),
115            'angular_scaling':angular_scaling(),'onset_boundary':onset_boundary(),'controller_sweep':controller_sweep()}
116    print(json.dumps(result,indent=2))
117
118if __name__=='__main__': main()