1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5
  6def rot(theta):
  7    return np.array([[math.cos(theta), -math.sin(theta)],
  8                     [math.sin(theta),  math.cos(theta)]], dtype=float)
  9
 10
 11def exact(v, R):
 12    n = np.linalg.norm(v)
 13    return v if n <= R else (R / n) * v
 14
 15
 16def smooth(v, R):
 17    # Differentiable radial surrogate from the idea.
 18    return v / math.sqrt(1.0 + (np.linalg.norm(v) / R) ** 2)
 19
 20
 21def map_x(x, gain, R, K, theta, limiter):
 22    # Centered affine map: x0 is an exact fixed point, so perturbation
 23    # dynamics are governed by the displayed Jacobian at x0.
 24    Q = rot(theta) @ np.diag([K, 1.0])
 25    x0 = np.array([1.2, 0.0])
 26    return x0 + gain * Q @ (limiter(x, R) - limiter(x0, R))
 27
 28
 29def fd_jacobian(fun, x, h=1e-6):
 30    J = np.zeros((2, 2))
 31    for j in range(2):
 32        xp, xm = x.copy(), x.copy()
 33        xp[j] += h
 34        xm[j] -= h
 35        J[:, j] = (fun(xp) - fun(xm)) / (2.0 * h)
 36    return J
 37
 38
 39def analytic_jacobian(R, r, K, theta, kind):
 40    t = r / R
 41    Q = rot(theta) @ np.diag([K, 1.0])
 42    if kind == "exact":
 43        # radial/tangential derivatives of projection outside the ball
 44        L = np.diag([0.0, 1.0 / t])
 45    else:
 46        # S(v)=v/sqrt(1+(||v||/R)^2)
 47        radial = (1.0 + t * t) ** (-1.5)
 48        tangential = (1.0 + t * t) ** (-0.5)
 49        L = np.diag([radial, tangential])
 50    return Q @ L
 51
 52
 53def rho(J):
 54    return float(np.max(np.abs(np.linalg.eigvals(J))))
 55
 56
 57def complex_pair(J):
 58    ev = np.linalg.eigvals(J)
 59    return bool(abs(ev[0].imag) > 1e-7 and abs(ev[1].imag) > 1e-7)
 60
 61
 62def predicted_exact_boundary(R, r, theta):
 63    # For Q=Rot(theta)diag(K,1), exact projection is rank one and
 64    # rho(g Q DP)=g*|cos(theta)|/(r/R), independent of K.
 65    return (r / R) / abs(math.cos(theta))
 66
 67
 68def predicted_smooth_boundary(R, r, K):
 69    t = r / R
 70    a = (1.0 + t * t) ** (-1.5)
 71    s = (1.0 + t * t) ** (-0.5)
 72    # In the complex-eigenvalue regime, rho(g Q DS)=g*sqrt(K*a*s).
 73    return 1.0 / math.sqrt(K * a * s)
 74
 75
 76def trajectory(gain, R, K, theta, kind, n=80, perturb=1e-5):
 77    limiter = exact if kind == "exact" else smooth
 78    x = np.array([1.2 + perturb, 0.0])
 79    norms = []
 80    for _ in range(n):
 81        norms.append(float(np.linalg.norm(x)))
 82        x = map_x(x, gain, R, K, theta, limiter)
 83        if not np.all(np.isfinite(x)) or np.linalg.norm(x) > 1e12:
 84            norms.extend([float("inf")] * (n - len(norms)))
 85            break
 86    return norms
 87
 88
 89def main():
 90    np.random.seed(7)
 91    R, r, theta = 1.0, 1.2, 1.2
 92    x0 = np.array([r, 0.0])
 93    Ks = [1.0, 2.0, 4.0, 8.0]
 94    rows = []
 95
 96    # Core math sanity check: analytic Jacobians versus finite differences.
 97    fd_errors = {}
 98    for kind, limiter in [("exact", exact), ("smooth", smooth)]:
 99        for K in [2.0, 4.0]:
100            f = lambda x: map_x(x, 1.0, R, K, theta, limiter)
101            Ja = analytic_jacobian(R, r, K, theta, kind)
102            Jfd = fd_jacobian(f, x0)
103            fd_errors[f"{kind}_K{int(K)}"] = float(np.max(np.abs(Ja - Jfd)))
104
105    # Prediction 1: exact threshold independent of anisotropy K.
106    exact_pred = predicted_exact_boundary(R, r, theta)
107    exact_obs = []
108    for K in Ks:
109        J = analytic_jacobian(R, r, K, theta, "exact")
110        exact_obs.append({"K": K, "predicted_gain": exact_pred,
111                          "observed_gain_from_rho": 1.0 / rho(J)})
112
113    # Prediction 2: smooth crossing scales as K^(-1/2), and is complex.
114    smooth_rows = []
115    for K in Ks:
116        pred = predicted_smooth_boundary(R, r, K)
117        Junit = analytic_jacobian(R, r, K, theta, "smooth")
118        observed = 1.0 / rho(Junit)
119        ev = np.linalg.eigvals(pred * Junit)
120        smooth_rows.append({"K": K, "predicted_gain": pred,
121                            "observed_gain_from_rho": observed,
122                            "complex_at_crossing": complex_pair(pred * Junit),
123                            "eigenvalues_at_crossing": [[float(z.real), float(z.imag)] for z in ev]})
124
125    # Prediction 3: increasing radius moves both boundaries linearly in r/R
126    # for exact projection, while the smooth boundary follows the analytic
127    # nonlinear formula; compare a radius sweep directly.
128    radius_rows = []
129    for rr in [0.8, 1.0, 1.2, 1.6, 2.0]:
130        ep = predicted_exact_boundary(R, rr, theta)
131        sp = predicted_smooth_boundary(R, rr, 4.0)
132        radius_rows.append({"r_over_R": rr / R, "exact_predicted": ep,
133                            "smooth_predicted_K4": sp})
134
135    # A direct continuation-style grid demonstrates the guard interval.
136    K = 4.0
137    ep, sp = exact_pred, predicted_smooth_boundary(R, r, K)
138    grid = np.linspace(0.8 * sp, 1.05 * ep, 25)
139    for g in grid:
140        Jp = g * analytic_jacobian(R, r, K, theta, "exact")
141        Js = g * analytic_jacobian(R, r, K, theta, "smooth")
142        rows.append({"gain": float(g), "rho_exact": rho(Jp),
143                     "rho_smooth": rho(Js), "smooth_complex": complex_pair(Js)})
144
145    # Perturbed trajectories: choose a gain strictly between the smooth and
146    # exact crossings, where the predicted artificial instability should occur.
147    middle = 0.5 * (sp + ep)
148    traj = {kind: trajectory(middle, R, K, theta, kind) for kind in ["exact", "smooth"]}
149    result = {
150        "setup": {"R": R, "r": r, "theta": theta, "Ks": Ks},
151        "finite_difference_max_errors": fd_errors,
152        "prediction_1_exact_boundary": exact_obs,
153        "prediction_2_smooth_complex_crossing": smooth_rows,
154        "prediction_3_radius_sweep": radius_rows,
155        "continuation_grid_K4": rows,
156        "trajectory_gain_between_boundaries": middle,
157        "trajectory_final_norms": {k: v[-1] for k, v in traj.items()},
158        "trajectory_max_norms": {k: float(np.max(v)) for k, v in traj.items()},
159        "trajectory_norms": traj,
160        "guard_interval": {"smooth_crossing": sp, "exact_crossing": ep,
161                           "width": ep - sp}
162    }
163    Path("results.json").write_text(json.dumps(result, indent=2))
164    print(json.dumps({"fd_errors": fd_errors, "exact_boundary": exact_pred,
165                      "smooth_K4_boundary": sp, "smooth_K4_complex": smooth_rows[2]["complex_at_crossing"],
166                      "guard_width": ep-sp, "trajectory_final_norms": result["trajectory_final_norms"]}, indent=2))
167
168
169if __name__ == "__main__":
170    main()