Each point on this sphere is one orbit shape of the 2D isotropic oscillator. A small perturbation εV₁ makes the shapes drift along the contours of its average ⟨V₁⟩, and the shapes that stop drifting are the periodic orbits that survive. However many survive, elliptic minus hyperbolic is χ(S²) = 2. Below, the same count is followed past the point where averaging fails, onto a torus, where it becomes 0, and onto a sphere with a cone point, where one orbit has to be left out.
Drag to rotate. Click anywhere on the sphere to start that orbit in the full equations (ε = 0.01) and trace its shape for about 950 periods. Color runs dark to bright with ⟨V₁⟩. Poles: x and y are straight-line motion along each axis, / and \ the diagonals, ccw and cw the circular motions.
Averaging stops predicting multiplier rates by ε ≈ 0.05 and stops predicting which orbits exist by ε ≈ 0.4. These sweeps continue every orbit to ε = 3 on the energy surface H = ½, with a blind search for new ones at each step.
Every orbit sits on one Floquet axis. Below 0 it is hyperbolic, between 0 and 1 elliptic, above 1 inverse-hyperbolic. Reaching 0 is a fold, where two orbits meet and vanish; crossing 1 is a period doubling.
As written, and with 0.1·x⁴ added, every term of V₁ is even in (q, p), so each orbit's monodromy is the square of a half-period map and can touch 1 but never cross it. Where the half-period map doubles instead (random 13, at ε ≈ 1.87 as written and 2.04 with x⁴), a mirror pair of new orbits appears.
With an odd term added, the squaring is gone and orbits that only touched 1 cross it, turning inverse-hyperbolic with no pair spawned. How far they cross follows how strongly the symmetry is broken. With 0.5·x³ the crossings are deep, down to k = −10. With 0.1·x⁵ they are shallow: in random 13 one orbit is inverse-hyperbolic only between ε = 1.46 and 1.79, and at most 6 × 10⁻⁴ past −4, still well above the numerical noise of 10⁻⁷. In all four systems E − H drifts while E + I − H stays at 2.
Every orbit of a charge in a uniform field is a cyclotron circle with the same period, so orbits are labeled by their centers, and the centers fill a torus. Averaging a weak potential over a circle of radius ρ multiplies each Fourier mode by J₀(|k|ρ). Near a zero of J₀ a mode switches off and returns with the opposite sign, so critical points are created and destroyed. The torus has χ = 0, so E + I − H should stay 0 however many there are.
Fig. 4 ⟨V⟩ over the torus of guiding centres (both edges wrap). Filled marks are the critical points averaging predicts; rings are the orbits found in the full equations at ε = 0.01. Circles are elliptic, diamonds hyperbolic.
Where averaging's list of orbits fails. At 9 of the 55 radii the orbits found do not match ⟨V⟩'s critical points one to one, so that prediction failed as written. Every failure sits next to a degeneracy of ⟨V⟩. At ρ ≈ 0.35 and 1.7 a min/saddle pair is about to merge, and in the full system it has already merged. Around ρ ≈ 2.4, the zero of J₀ for the two strongest modes, the landscape is so flat that orbits sit up to 0.23 away from the predicted points. Halving ε halves those offsets and brings the missing pairs back, so these are O(ε) corrections to averaging, not missed orbits. E + I − H is 0 at all 55 radii either way.
With frequencies 1 and 2, every orbit has period 2π except the y-mode, which closes in half the time and sits at a cone point of the orbit space. Counted at period 2π it appears as a double cover. The rule from Frauenfelder, Koh and Moreno leaves a double cover out when the orbit itself is negative-hyperbolic. The coupling κ·x²y pushes the y-mode into its parametric-resonance tongue near κ* = 1.125, which is exactly that case.
Fig. 7 ⟨V₁⟩ along the resonance axis of the orbit space, from the cone point through the opposite pole and back to the cone. Every critical point away from the cone lies on this axis. The cone is the y-mode.
This is the SFT-Euler characteristic of Frauenfelder, Koh and Moreno, which is −(E + I − H) in the notation here. Elliptic and negative-hyperbolic orbits have odd Conley–Zehnder index and positive-hyperbolic ones even, which is why E groups with I against H. It cannot change through a bifurcation, by invariance of local Floer homology. At small ε its value is set by the averaged function: the Euler characteristic of the orbit space.
They run that bookkeeping on systems it had not been run on: the pair born at a symmetric doubling, inverse-hyperbolic parents once the symmetry is broken, zero on a torus as J₀ reshapes the landscape, and the one bad orbit at a 1:2 resonance. Every index count held. The one prediction that failed was averaging's list of orbits on the torus, next to its own degeneracies. That checks the numerics against the theorem; it does not prove anything new.
A rigorous count at small ε (averaging with explicit bounds), the invariance theorem, and a proof that no periodic orbit enters or leaves the period window across the sweep. Each listed orbit's existence and type can be proven with interval arithmetic: a Krawczyk test on the shooting map and an enclosure of tr M − 4 away from 0 and −4.
Weinstein and Moser for the small-ε count on S²; Ginzburg for at least four orbits of a charge on a flat torus at low energy; Frauenfelder, Koh and Moreno (2022) for χSFT, good and bad orbits, and the generic bifurcations it survives.
Every prediction below was written into the script before its run, together with the rule for when a point counts: no orbit may sit within numerical noise of a type change.
One fix after a failed check. The first run of the torus and 1:2 sweeps reported the 1:2 count failing at one point, κ = 1.1276, just past the y-mode's doubling. The new orbit was there, but its period and mean position matched the y-mode's double cover to 2 × 10⁻⁶, and my duplicate filter merged the two. Adding each orbit's width to the filter separates them by 10⁻⁴. All four sweeps were rerun with the fixed filter; the table shows the rerun.
Python script check_hopf.py. Periodic orbits are solved on the perturbed energy surface H = H₀ + εV₁ = ½. A case passes only if every check holds:
Checks on this page's own code. The in-browser solver for the sphere's critical points agreed with the exact enumeration on 5,000 random slider settings and 2,200 settings near saddle-node bifurcations, disagreeing only where two points were 3 × 10⁻⁶ rad apart, which it labels degenerate. The RK4 probe differs from a run at 4× finer step by at most 0.5° over 955 periods, H drifts by at most 6 × 10⁻⁵, and the probe's ⟨V₁⟩ spread is unchanged at the finer step.
"""Perturbed 1:1 oscillator on the Hopf base S^2.
Claim: periodic orbits near period 2*pi on one energy level satisfy
#elliptic - #hyperbolic = chi(S^2) = 2.
H = (px^2+py^2+x^2+y^2)/2 + EPS*V1,
V1 = x^4 + y^4 + c x^2 y^2 + d (x^2-y^2) + e x y + r (x py - y px).
"""
import json, time
import numpy as np
from scipy.integrate import solve_ivp
EPS, S0, H_LEVEL = 0.01, 1.0, 0.5
rng = np.random.default_rng(0)
class Sys:
def __init__(s, c, d, e, r):
s.p = dict(c=c, d=d, e=e, r=r)
s.c, s.d, s.e, s.r = c, d, e, r
def V1(s, X):
x, y, px, py = X
return x**4 + y**4 + s.c*x*x*y*y + s.d*(x*x - y*y) + s.e*x*y + s.r*(x*py - y*px)
def H(s, X):
x, y, px, py = X
return 0.5*(px*px + py*py + x*x + y*y) + EPS*s.V1(X)
def rhs(s, X):
x, y, px, py = X
c, d, e, r = s.c, s.d, s.e, s.r
return np.array([px - EPS*r*y,
py + EPS*r*x,
-x - EPS*(4*x**3 + 2*c*x*y*y + 2*d*x + e*y + r*py),
-y - EPS*(4*y**3 + 2*c*x*x*y - 2*d*y + e*x - r*px)])
def jac(s, X):
x, y, px, py = X
c, d, e, r = s.c, s.d, s.e, s.r
J = np.zeros((4, 4))
J[0, 1], J[0, 2] = -EPS*r, 1
J[1, 0], J[1, 3] = EPS*r, 1
J[2, 0] = -1 - EPS*(12*x*x + 2*c*y*y + 2*d)
J[2, 1] = -EPS*(4*c*x*y + e)
J[2, 3] = -EPS*r
J[3, 0] = -EPS*(4*c*x*y + e)
J[3, 1] = -1 - EPS*(12*y*y + 2*c*x*x - 2*d)
J[3, 2] = EPS*r
return J
def gradH(s, X):
f = s.rhs(X)
return np.array([-f[2], -f[3], f[0], f[1]])
# averaged V1 on the Hopf base: f(u) = u.A.u + b.u + const, u = unit Stokes vector
def Ab(s):
A = S0**2*np.diag([(3 - s.c/2)/16, s.c/16, 0.0])
b = S0*np.array([s.d/2, s.e/4, s.r/2])
const = S0**2*(3/16 + s.c/32)
return A, b, const
def lift(u):
u = u/np.linalg.norm(u)
a = np.sqrt(S0*(1 + u[0])/2)
b = np.sqrt(S0*max(1 - u[0], 0)/2)*np.exp(1j*np.arctan2(u[2], u[1]))
return np.array([a.real, b.real, a.imag, b.imag]) # x, y, px, py
def hopf(X):
x, y, px, py = X
a, b = x + 1j*px, y + 1j*py
ab = np.conj(a)*b
S = np.array([abs(a)**2 - abs(b)**2, 2*ab.real, 2*ab.imag])
return S/np.linalg.norm(S, axis=0)
def fib(n):
i = np.arange(n) + 0.5
z = 1 - 2*i/n
t = np.pi*(1 + 5**0.5)*i
return np.stack([np.sqrt(1 - z*z)*np.cos(t), np.sqrt(1 - z*z)*np.sin(t), z], 1)
# ---- S1: analytic average vs numerical average over Hopf circles
def check_average(sy, n=50):
A, b, k = sy.Ab()
th = np.linspace(0, 2*np.pi, 64, endpoint=False)
err = 0.0
for u in rng.normal(size=(n, 3)):
u /= np.linalg.norm(u)
X = lift(u)
za, zb = (X[0] + 1j*X[2])*np.exp(-1j*th), (X[1] + 1j*X[3])*np.exp(-1j*th)
num = sy.V1(np.array([za.real, zb.real, za.imag, zb.imag])).mean()
err = max(err, abs(num - (u@A@u + b@u + k)))
return err
# ---- S2: critical points of the averaged function on S^2
def critical_points(sy):
A, b, _ = sy.Ab()
out = []
for u in fib(400):
lam = u@(2*A@u + b)
for _ in range(60):
F = np.concatenate([2*A@u + b - lam*u, [u@u - 1]])
if np.linalg.norm(F) < 1e-14:
break
J = np.zeros((4, 4))
J[:3, :3] = 2*A - lam*np.eye(3)
J[:3, 3] = -u
J[3, :3] = 2*u
try:
st = np.linalg.solve(J, -F)
except np.linalg.LinAlgError:
break
u, lam = u + st[:3], lam + st[3]
if np.linalg.norm(F) > 1e-12 or any(np.linalg.norm(u - v["u"]) < 1e-6 for v in out):
continue
t1 = np.cross(u, [1, 0, 0] if abs(u[0]) < 0.9 else [0, 1, 0]); t1 /= np.linalg.norm(t1)
t2 = np.cross(u, t1)
T = np.stack([t1, t2], 1)
ev = np.linalg.eigvalsh(T.T@(2*A - lam*np.eye(3))@T)
kind = "min" if ev.min() > 0 else "max" if ev.max() < 0 else "saddle"
out.append(dict(u=u, kind=kind, hess=ev))
return out
# ---- S2b: exact enumeration, independent of Newton seeding.
# (p_k - lam) u_k = -b_k with p = 2 diag(A), so |u| = 1 becomes g(lam) = sum b_k^2/(p_k - lam)^2 - 1 = 0.
# g is monotone outside the poles (one root each side) and convex between adjacent poles (0 or 2 roots),
# so bisection on each interval finds every root. b_k = 0 adds solutions at lam = p_k.
# Returns (points, gap); gap -> 0 at a saddle-node of <V1>. Returns None for a circle of critical points.
def algebraic_critical(sy):
A, b, _ = sy.Ab()
p = 2*np.diag(A)
act = [k for k in range(3) if b[k] != 0]
g = lambda l: sum(b[k]**2/(p[k] - l)**2 for k in act) - 1
dg = lambda l: sum(2*b[k]**2/(p[k] - l)**3 for k in act)
def bis(f, lo, hi): # needs f < 0 near lo, f > 0 near hi; never evaluates the endpoints
for _ in range(400):
m = 0.5*(lo + hi)
if m <= lo or m >= hi:
break
lo, hi = (m, hi) if f(m) < 0 else (lo, m)
return 0.5*(lo + hi)
lams, gap = [], np.inf
q = sorted(set(p[k] for k in act))
if q:
R = 2*np.linalg.norm(b) + 1
lams += [bis(g, q[0] - R, q[0]), bis(lambda l: -g(l), q[-1], q[-1] + R)]
for lo, hi in zip(q, q[1:]):
ls = bis(dg, lo, hi)
gm = g(ls)
gap = min(gap, abs(gm))
if gm < 0:
lams += [bis(lambda l: -g(l), lo, ls), bis(g, ls, hi)]
pts = [np.array([-b[k]/(p[k] - l) if b[k] != 0 else 0.0 for k in range(3)]) for l in lams]
for k in range(3):
if b[k] != 0:
continue
l = p[k]
if any(j != k and p[j] == l for j in range(3)):
if any(j != k and p[j] == l and b[j] == 0 for j in range(3)):
return None
continue
u = np.array([-b[j]/(p[j] - l) if j != k and b[j] != 0 else 0.0 for j in range(3)])
rem = 1 - u@u
gap = min(gap, abs(rem))
if rem > 0:
for s in (1, -1):
v = u.copy(); v[k] = s*np.sqrt(rem); pts.append(v)
return pts, gap
# ---- S3: periodic orbits of the full system by Newton shooting
def flow(sy, X0, T, stm=True, rtol=1e-11):
if not stm:
return solve_ivp(lambda t, X: sy.rhs(X), (0, T), X0, method="DOP853",
rtol=1e-11, atol=1e-12, dense_output=True)
def f(t, Y):
return np.concatenate([sy.rhs(Y[:4]), (sy.jac(Y[:4])@Y[4:].reshape(4, 4)).ravel()])
Y = solve_ivp(f, (0, T), np.concatenate([X0, np.eye(4).ravel()]),
method="DOP853", rtol=rtol, atol=rtol/10).y[:, -1]
return Y[:4], Y[4:].reshape(4, 4)
def shoot(sy, X, T=2*np.pi, maxit=40, maxstep=0.1):
for _ in range(maxit):
XT, M = flow(sy, X, T)
F = np.concatenate([XT - X, [sy.H(X) - H_LEVEL], [0.0]])
if np.linalg.norm(F) < 1e-10:
return dict(X=X, T=T, M=M, res=np.linalg.norm(F))
J = np.zeros((6, 5))
J[:4, :4] = M - np.eye(4)
J[:4, 4] = sy.rhs(XT)
J[4, :4] = sy.gradH(X)
J[5, :4] = sy.rhs(X)
dz = np.linalg.lstsq(J, -F, rcond=None)[0]
n = np.linalg.norm(dz)
if n > maxstep:
dz *= maxstep/n
X, T = X + dz[:4], T + dz[4]
return None
def orbit_center(sy, orb):
sol = flow(sy, orb["X"], orb["T"], stm=False)
S = hopf(sol.sol(np.linspace(0, orb["T"], 64, endpoint=False))).mean(1)
return S/np.linalg.norm(S)
def floquet(orb):
k = np.trace(orb["M"]) - 4 # = lambda + 1/lambda - 2 for the nontrivial multiplier pair
return ("elliptic" if -4 < k < 0 else "hyperbolic" if k > 0 else "inverse-hyperbolic"), k
def blind_search(sy, nseed=400): # max basins can be ~0.1 rad wide; 400 seeds ≈ 0.1 rad apart
found = []
for u in fib(nseed):
orb = shoot(sy, lift(u))
if orb is None:
continue
cen = orbit_center(sy, orb)
if all(np.linalg.norm(cen - f["center"]) > 0.05 for f in found):
found.append(dict(center=cen, type=floquet(orb)[0]))
return found
def run_case(name, params, blind=False):
sy = Sys(**params)
cps = critical_points(sy)
n_ext = sum(c["kind"] != "saddle" for c in cps)
n_sad = sum(c["kind"] == "saddle" for c in cps)
rows, ok = [], True
for cp in cps:
orb = shoot(sy, lift(cp["u"]))
if orb is None:
ok = False
rows.append(dict(kind=cp["kind"], converged=False))
continue
cen = orbit_center(sy, orb)
typ, k = floquet(orb)
match = (typ == "elliptic") == (cp["kind"] != "saddle") and typ != "inverse-hyperbolic"
ok &= match
XT2, M2 = flow(sy, orb["X"], orb["T"], rtol=1e-13) # same orbit, 100x tighter tolerance
k2 = floquet(dict(M=M2))[1]
# averaging predicts multiplier angle (elliptic) or exponent (hyperbolic) = 4*eps*sqrt|h1 h2| per unit time
ang = np.arccos(1 + k/2) if typ == "elliptic" else np.arccosh(1 + k/2) if typ == "hyperbolic" else np.nan
kappa = ang/(orb["T"]*EPS*np.sqrt(abs(np.prod(cp["hess"]))))
rows.append(dict(kind=cp["kind"], u=cp["u"].round(4).tolist(), floquet=typ,
trM_minus_4=float(k), trM_minus_4_tight=float(k2), T=float(orb["T"]), res=float(orb["res"]),
res_tight=float(np.linalg.norm(XT2 - orb["X"])), dH=float(abs(sy.H(orb["X"]) - H_LEVEL)),
kappa=float(kappa), shift=float(np.linalg.norm(cen - cp["u"])), match=bool(match),
center=cen))
n_ell = sum(r.get("floquet") == "elliptic" for r in rows)
n_hyp = sum(r.get("floquet") == "hyperbolic" for r in rows)
centers = [r["center"] for r in rows if "center" in r]
distinct = all(np.linalg.norm(a - b) > 0.05 for i, a in enumerate(centers) for b in centers[i+1:])
out = dict(name=name, params=params, n_ext=n_ext, n_saddle=n_sad,
n_elliptic=n_ell, n_hyperbolic=n_hyp, distinct=bool(distinct),
max_res=max((r.get("res", 1) for r in rows), default=1),
max_shift=max((r.get("shift", 9) for r in rows), default=9),
types_match=bool(ok))
alg = algebraic_critical(sy)
out["alg_count"] = None if alg is None else len(alg[0])
out["alg_match"] = bool(alg is not None and len(alg[0]) == len(cps) and all(
any(np.linalg.norm(a - c["u"]) < 1e-8 for a in alg[0]) for c in cps))
out["max_dH"] = max((r.get("dH", 1) for r in rows), default=1)
out["max_res_tight"] = max((r.get("res_tight", 1) for r in rows), default=1)
out["min_floquet_snr"] = min((abs(r["trM_minus_4"])/max(abs(r["trM_minus_4"] - r["trM_minus_4_tight"]), 1e-300)
for r in rows if "trM_minus_4" in r), default=0)
ks = [r["kappa"] for r in rows if "kappa" in r]
out["kappa_range"] = [min(ks), max(ks)] if ks else None
if blind:
bf = blind_search(sy)
same = len(bf) == len(centers) and all(
any(np.linalg.norm(f["center"] - c) < 0.05 for c in centers) for f in bf)
out["blind_found"] = len(bf)
out["blind_matches"] = bool(same)
for r in rows:
r.pop("center", None)
out["orbits"] = rows
return out
def well_separated(params):
cps = critical_points(Sys(**params))
if any(abs(c["hess"]).min() < 2e-3 for c in cps):
return False
us = [c["u"] for c in cps]
return all(np.linalg.norm(a - b) > 0.2 for i, a in enumerate(us) for b in us[i+1:])
if __name__ == "__main__":
t0 = time.time()
presets = [
("quartic c=1", dict(c=1.0, d=0.0, e=0.0, r=0.0)),
("quartic c=8", dict(c=8.0, d=0.0, e=0.0, r=0.0)),
("anisotropy d=0.5", dict(c=1.0, d=0.5, e=0.0, r=0.0)),
("anisotropy d=1", dict(c=1.0, d=1.0, e=0.0, r=0.0)),
("mixed", dict(c=2.0, d=0.15, e=0.3, r=0.05)),
]
avg_err = max(check_average(Sys(**p)) for _, p in presets)
print(f"S1 max |analytic - numerical average| = {avg_err:.2e}")
results = [run_case(n, p, blind=True) for n, p in presets]
rejected, k = 0, 0
while k < 16:
p = dict(c=float(rng.uniform(-3, 10)), d=float(rng.uniform(-.4, .4)),
e=float(rng.uniform(-.4, .4)), r=float(rng.uniform(-.4, .4)))
if not well_separated(p):
rejected += 1
continue
k += 1
results.append(run_case(f"random {k}", {q: round(v, 4) for q, v in p.items()}))
print(f"{'case':18s} ext sad alg | ell hyp diff | res res@tight |H-1/2| floq.snr kappa match distinct blind")
allok = avg_err < 1e-12
for R in results:
diff = R["n_elliptic"] - R["n_hyperbolic"]
b = f'{R["blind_found"]} {"ok" if R["blind_matches"] else "MISMATCH"}' if "blind_found" in R else "-"
kr = R["kappa_range"]
print(f'{R["name"]:18s} {R["n_ext"]:3d} {R["n_saddle"]:3d} {"ok " if R["alg_match"] else "BAD"} |'
f' {R["n_elliptic"]:3d} {R["n_hyperbolic"]:3d} {diff:4d} |'
f' {R["max_res"]:.1e} {R["max_res_tight"]:.1e} {R["max_dH"]:.1e} {R["min_floquet_snr"]:8.1e} {kr[0]:.2f}-{kr[1]:.2f}'
f' {R["types_match"]} {R["distinct"]} {b}')
allok &= (R["n_ext"] - R["n_saddle"] == 2 and diff == 2 and R["types_match"]
and R["distinct"] and R["max_res"] < 1e-10 and R.get("blind_matches", True)
and R["alg_match"] and R["max_dH"] < 1e-10 and R["max_res_tight"] < 1e-9
and R["min_floquet_snr"] > 100)
print(f"rejected near-degenerate random draws: {rejected}")
# kappa should approach the averaging prediction 4 as eps -> 0, with an O(eps) gap
kappa_eps = []
for e in (0.01, 0.005, 0.0025):
EPS = e
kr = run_case("anisotropy d=0.5", presets[2][1])["kappa_range"]
kappa_eps.append([e, kr])
print(f"eps={e:<7} kappa range {kr[0]:.4f}-{kr[1]:.4f} max |kappa-4| = {max(abs(k-4) for k in kr):.4f}")
EPS = 0.01
print("ALL CHECKS PASS" if allok else "SOME CHECK FAILED")
print(f"elapsed {time.time()-t0:.0f}s")
json.dump(dict(eps=EPS, avg_err=avg_err, rejected=rejected, results=results, kappa_eps=kappa_eps),
open("hopf_results.json", "w"), indent=1)
"""eps sweep for the perturbed 1:1 oscillator of check_hopf.py, on the energy surface H = 1/2.
Every term of V1 is even, so P(q, p) = -(q, p) is a symmetry. An orbit is symmetric if flow(T/2)(X) = -X.
For those, G = P o flow(T/2) fixes X and the monodromy is M = (DG)^2, so M cannot have a negative real multiplier.
Classification (k = tr - 4 of the nontrivial multiplier pair; E elliptic, H hyperbolic, I inverse-hyperbolic):
full k = tr M - 4 every orbit of period ~T
quotient k' = tr DG - 4 = -tr flow'(T/2) - 4 symmetric orbits only
Predictions, fixed before the run:
P1 quotient: E'-H' = 2 until the first I' (mu = -1 of G); E'+I'-H' = 2 at every eps
P2 full: no symmetric orbit is ever I; E+I-H = 2 over all orbits of period ~T, nonsymmetric pairs included
"""
import json, sys, time
import numpy as np
from multiprocessing import Pool
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import check_hopf as h
CASES = {"mixed": dict(c=2.0, d=0.15, e=0.3, r=0.05),
"random 5": dict(c=1.8834, d=0.3526, e=-0.1185, r=-0.0551),
"random 13": dict(c=8.2995, d=-0.0494, e=0.3065, r=0.0604)}
BASE = np.unique(np.concatenate([np.geomspace(0.01, 2.0, 31), np.linspace(2.0, 3.0, 6)]))
def stm(sy, X, T, rtol):
def f(t, Y):
return np.concatenate([sy.rhs(Y[:4]), (sy.jac(Y[:4]) @ Y[4:].reshape(4, 4)).ravel()])
s = solve_ivp(f, (0, T), np.concatenate([X, np.eye(4).ravel()]), method="DOP853",
rtol=rtol, atol=rtol/10, t_eval=[T/2, T])
return [(s.y[:4, i], s.y[4:, i].reshape(4, 4)) for i in range(2)]
def kind(k, dk):
if min(abs(k), abs(k + 4)) < max(1e3*dk, 1e-8):
return "unsure"
return "E" if -4 < k < 0 else "H" if k > 0 else "I"
def signature(sy, X, T):
pts = h.flow(sy, X, T, stm=False).sol(np.linspace(0, T, 128, endpoint=False))
c = h.hopf(pts).mean(1)
return dict(T=float(T), center=c/np.linalg.norm(c), mean=pts.mean(1))
def same(a, b):
return (abs(a["T"] - b["T"]) < 1e-6 and np.linalg.norm(a["center"] - b["center"]) < 1e-5
and np.linalg.norm(a["mean"] - b["mean"]) < 1e-5)
def describe(sy, X, T, sig=None):
(Xh, Mh), (XT, M) = stm(sy, X, T, 1e-11)
(_, Mh2), (_, M2) = stm(sy, X, T, 1e-13)
k, k2 = np.trace(M) - 4, np.trace(M2) - 4
d = dict(sig or signature(sy, X, T), X=X, res=float(np.linalg.norm(XT - X)),
dH=float(abs(sy.H(X) - h.H_LEVEL)), sym=float(np.linalg.norm(Xh + X)),
k=float(k), full=kind(k, abs(k - k2)))
if d["sym"] < 1e-7:
q, q2 = -np.trace(Mh) - 4, -np.trace(Mh2) - 4
d.update(kq=float(q), quot=kind(q, abs(q - q2)), square_err=float(abs((q + 2)**2 - 4 - k)))
w, v = np.linalg.eig(-Mh) # DG; eigenvector for the multiplier nearest -1
d["v_pd"] = np.real(v[:, np.argmin(abs(w + 1))])
return d
# ---- symmetric shooting: solve flow(T/2)(X) = -X. Its Jacobian flow'(T/2) + I = I - DG stays nonsingular
# through mu = -1 of G, where full-period shooting (M - I, M = DG^2) is singular and can jump to the new pair.
def shoot_sym(sy, X, T, maxit=40, maxstep=0.1):
for _ in range(maxit):
Xh, Mh = h.flow(sy, X, T/2)
F = np.concatenate([Xh + X, [sy.H(X) - h.H_LEVEL], [0.0]])
if np.linalg.norm(F) < 1e-10:
return dict(X=X, T=T)
J = np.zeros((6, 5))
J[:4, :4], J[:4, 4] = Mh + np.eye(4), 0.5*sy.rhs(Xh)
J[4, :4], J[5, :4] = sy.gradH(X), sy.rhs(X)
dz = np.linalg.lstsq(J, -F, rcond=None)[0]
dz *= min(1, maxstep/np.linalg.norm(dz))
X, T = X + dz[:4], T + dz[4]
return None
# ---- continuation of each averaged critical point's orbit in eps (secant predictor, step halving)
def advance(sy, path, e_target):
while path[-1][0] < e_target:
e0, X0, T0 = path[-1]
step = e_target - e0
while True:
e1 = e_target if step >= e_target - e0 else e0 + step
if len(path) >= 2:
ea, Xa, Ta = path[-2]
w = (e1 - e0)/(e0 - ea)
Xg, Tg = X0 + w*(X0 - Xa), T0 + w*(T0 - Ta)
else:
Xg, Tg = X0, T0
h.EPS = e1
o = shoot_sym(sy, Xg, Tg)
if o is not None and np.linalg.norm(o["X"] - Xg) < 0.05 and abs(o["T"] - Tg) < 0.05:
path.append((e1, o["X"], o["T"]))
del path[:-2]
break
step /= 2
if step < 1e-7:
return False
return True
def continue_case(name, grid):
sy = h.Sys(**CASES[name])
out = []
for cp in h.critical_points(sy):
h.EPS = grid[0]
o = h.shoot(sy, h.lift(cp["u"]))
path, rows, end = [(grid[0], o["X"], o["T"])], [], None
for e in grid:
if not advance(sy, path, e):
end = float(path[-1][0])
break
h.EPS = e
d = describe(sy, path[-1][1], path[-1][2])
if "kq" in d: # averaging: quotient angle = 2 eps sqrt|h1h2| T
phi = np.arccos(np.clip(1 + d["kq"]/2, -1, 1)) if d["quot"] == "E" else np.nan
d["kappa"] = float(2*phi/(d["T"]*e*np.sqrt(abs(np.prod(cp["hess"])))))
rows.append(dict(eps=float(e), **d))
out.append(dict(kind=cp["kind"], u=cp["u"].tolist(), rows=rows, end=end))
return name, out
def events(name):
"""First pass on BASE: locate period-doublings of G (k' crosses -4) by bisection, and fold ends."""
_, brs = continue_case(name, BASE)
sy = h.Sys(**CASES[name])
ev = []
for b in brs:
for r0, r1 in zip(b["rows"], b["rows"][1:]):
if "kq" in r0 and "kq" in r1 and (r0["kq"] + 4)*(r1["kq"] + 4) < 0:
lo, hi, X, T = r0["eps"], r1["eps"], r0["X"], r0["T"]
for _ in range(30):
m, path = 0.5*(lo + hi), [(lo, X, T)]
if not advance(sy, path, m):
break
h.EPS = m
d = describe(sy, path[-1][1], path[-1][2])
if (d["kq"] + 4)*(r0["kq"] + 4) > 0:
lo, X, T = m, path[-1][1], path[-1][2]
else:
hi = m
ev.append(("period-doubling of G", float(lo)))
if b["end"] is not None:
ev.append(("fold", b["end"]))
return name, ev
# ---- batched Newton shooting: all seeds integrated together as one ODE in normalized time tau = t/T
def jac_v(sy, X):
x, y, px, py = X
E, c, d, e, r = h.EPS, sy.c, sy.d, sy.e, sy.r
J = np.zeros((X.shape[1], 4, 4))
J[:, 0, 1], J[:, 0, 2], J[:, 1, 0], J[:, 1, 3] = -E*r, 1, E*r, 1
J[:, 2, 0] = -1 - E*(12*x*x + 2*c*y*y + 2*d)
J[:, 2, 1] = J[:, 3, 0] = -E*(4*c*x*y + e)
J[:, 2, 3], J[:, 3, 2] = -E*r, E*r
J[:, 3, 1] = -1 - E*(12*y*y + 2*c*x*x - 2*d)
return J
def batch_flow(sy, X, T, dense=0):
n = len(T)
def f(tau, Y):
Y = Y.reshape(n, 20)
Xs = Y[:, :4].T
dX = sy.rhs(Xs).T*T[:, None]
dP = T[:, None, None]*(jac_v(sy, Xs) @ Y[:, 4:].reshape(n, 4, 4))
return np.concatenate([dX, dP.reshape(n, 16)], 1).ravel()
Y0 = np.concatenate([X, np.tile(np.eye(4).ravel(), (n, 1))], 1).ravel()
tev = np.linspace(0, 1, dense, endpoint=False) if dense else [1.0]
s = solve_ivp(f, (0, 1), Y0, method="DOP853", rtol=1e-11, atol=1e-12, t_eval=tev).y.reshape(n, 20, -1)
return s if dense else (s[:, :4, -1], s[:, 4:, -1].reshape(n, 4, 4))
def batch_shoot(sy, X, T, maxit=40, maxstep=0.1):
X, T, T0 = np.array(X), np.array(T, float), np.array(T, float)
act, done, best = np.ones(len(T), bool), np.zeros(len(T), bool), np.full(len(T), np.inf)
for it in range(maxit):
idx = np.where(act)[0]
if not len(idx):
break
XT, M = batch_flow(sy, X[idx], T[idx])
for j, i in enumerate(idx):
F = np.concatenate([XT[j] - X[i], [sy.H(X[i]) - h.H_LEVEL], [0.0]])
nf = np.linalg.norm(F)
if nf < 1e-9:
act[i], done[i] = False, True
continue
best[i] = min(best[i], nf)
if not np.isfinite(nf) or (it >= 15 and nf > 1e-3) or np.abs(X[i]).max() > 5 or abs(T[i] - T0[i]) > 0.5*T0[i]:
act[i] = False
continue
Jm = np.zeros((6, 5))
Jm[:4, :4], Jm[:4, 4] = M[j] - np.eye(4), sy.rhs(XT[j])
Jm[4, :4], Jm[5, :4] = sy.gradH(X[i]), sy.rhs(X[i])
dz = np.linalg.lstsq(Jm, -F, rcond=None)[0]
dz *= min(1, maxstep/np.linalg.norm(dz))
X[i], T[i] = X[i] + dz[:4], T[i] + dz[4]
return X[done], T[done]
# ---- searches at one eps: 400 on-shell seeds, plus seeds pushed off each symmetric orbit along DG's mu~-1 eigenvector
def search(args):
name, e, rows = args
h.EPS = e
sy = h.Sys(**CASES[name])
seeds = [] # (state, period guess)
for u in h.fib(400):
X = h.lift(u)
Tg = min(rows, key=lambda r: np.linalg.norm(r["center"] - u))["T"] # period of the nearest branch
try:
seeds.append((brentq(lambda s: sy.H(s*X) - h.H_LEVEL, 0.1, 3.0)*X, Tg))
except ValueError:
pass
for r in rows:
if "v_pd" in r:
for dlt in (0.003, 0.01, 0.03, 0.1):
seeds += [(r["X"] + dlt*r["v_pd"], r["T"]), (r["X"] - dlt*r["v_pd"], r["T"])]
Tlo, Thi = 0.75*min(r["T"] for r in rows), 1.25*max(r["T"] for r in rows)
Xc, Tc = batch_shoot(sy, [s[0] for s in seeds], [s[1] for s in seeds])
keep = (Tc > Tlo) & (Tc < Thi)
Xc, Tc = Xc[keep], Tc[keep]
found = []
if len(Tc): # dedupe on a coarse signature, then polish one per orbit
pts = batch_flow(sy, Xc, Tc, dense=128)[:, :4, :]
for X, T, P in zip(Xc, Tc, pts):
c = h.hopf(P).mean(1)
s = dict(T=float(T), center=c/np.linalg.norm(c), mean=P.mean(1))
if not any(abs(s["T"] - f["T"]) < 1e-5 and np.linalg.norm(s["center"] - f["center"]) < 1e-4
and np.linalg.norm(s["mean"] - f["mean"]) < 1e-4 for f in found):
found.append(dict(s, X=X))
out = []
for f in found:
o = h.shoot(sy, f["X"], T=f["T"])
if o is not None:
d = describe(sy, o["X"], o["T"])
if not any(same(d, g) for g in out):
out.append(d)
return name, float(e), out
def tojson(o):
if isinstance(o, dict):
return {k: tojson(v) for k, v in o.items()}
if isinstance(o, (list, tuple)):
return [tojson(v) for v in o]
if isinstance(o, np.ndarray):
return o.tolist()
if isinstance(o, np.floating):
return float(o)
return o
if __name__ == "__main__":
t0 = time.time()
names = sys.argv[2:] or list(CASES)
with Pool(2) as pool:
ev = dict(pool.map(events, names))
print(f"events {ev} {time.time()-t0:.0f}s", flush=True)
grids = {}
for n in names:
extra = [x for k, e in ev[n] for x in
((e - 0.002, e + 0.002, e + 0.01, e + 0.03) if k != "fold" else (e - 0.002, e + 0.002))]
grids[n] = np.unique(np.concatenate([BASE, [x for x in extra if BASE[0] < x < BASE[-1]]]))
cont = dict(pool.starmap(continue_case, [(n, grids[n]) for n in names]))
print(f"continuation done {time.time()-t0:.0f}s", flush=True)
jobs = []
for n in names:
for i, e in enumerate(grids[n]):
jobs.append((n, float(e), [b["rows"][i] for b in cont[n] if i < len(b["rows"])]))
found = pool.map(search, jobs, chunksize=1)
json.dump(tojson(dict(grids=grids, events=ev, cases=CASES, continuation=cont, search=found)),
open(sys.argv[1], "w"))
print(f"done {time.time()-t0:.0f}s", flush=True)
"""Evaluate P1/P2 on sweep_results.json. Rules are fixed here, before looking at the results:
an eps counts only if no orbit in the set is 'unsure' (multiplier within numerical noise of a transition);
P1 quotient (symmetric orbits, half-period map G): E'+I'-H' = 2 at every counted eps;
E'-H' = 2 at every counted eps before the first I', and != 2 wherever I' > 0
P2 full system, all distinct orbits of period ~T: no symmetric orbit is ever I; E+I-H = 2 at every counted eps;
every nonsymmetric orbit has a partner with the same period and opposite mean position
"""
import json, sys
import numpy as np
D = json.load(open(sys.argv[1]))
def same(a, b):
return (abs(a["T"] - b["T"]) < 1e-6 and np.linalg.norm(np.subtract(a["center"], b["center"])) < 1e-5
and np.linalg.norm(np.subtract(a["mean"], b["mean"])) < 1e-5)
search = {(n, round(e, 12)): f for n, e, f in D["search"]}
ok_all = True
for name in D["cases"]:
grid, brs = D["grids"][name], D["continuation"][name]
print(f"\n== {name} events: " + ", ".join(f"{k} at eps={e:.4f}" for k, e in D["events"][name]))
print(" eps | orbits sym+non | full E H I E-H E+I-H | quot E' H' I' E'-H' E'+I'-H' | search new/missed | notes")
first_I, avg_break, p1, p2, pairs_ok, symI, counted = None, None, True, True, True, False, 0
for i, e in enumerate(grid):
F = [b["rows"][i] for b in brs if i < len(b["rows"])]
nb = len(F)
found = search[(name, round(e, 12))]
missed = sum(not any(same(r, f) for f in found) for r in F)
new = 0
for f in found:
if not any(same(f, g) for g in F):
F.append(f); new += 1
Q = [r for r in F if "quot" in r]
nonsym = [r for r in F if "quot" not in r]
unsure = [r for r in F if r["full"] == "unsure" or r.get("quot") == "unsure"]
cnt = lambda S, key, t: sum(r.get(key) == t for r in S)
E, H_, I = (cnt(F, "full", t) for t in "EHI")
Eq, Hq, Iq = (cnt(Q, "quot", t) for t in "EHI")
symI |= any(r["full"] == "I" for r in Q)
for r in nonsym: # partner: same T, opposite mean
pairs_ok &= any(abs(r["T"] - s["T"]) < 1e-6 and np.linalg.norm(np.add(r["mean"], s["mean"])) < 1e-5
for s in nonsym if s is not r)
if first_I is None and Iq:
first_I = e
# averaging's prediction: the eps=0.01 orbit set, extrema E/E and saddles H/H, nothing added or lost
pred_ok = nb == len(brs) and not nonsym and all(
(r["full"], r["quot"]) == (("H", "H") if b["kind"] == "saddle" else ("E", "E"))
for b, r in zip(brs, F[:nb]))
if avg_break is None and not pred_ok:
avg_break = e
note = []
if unsure:
note.append(f"{len(unsure)} unsure (not counted)")
else:
counted += 1
p1 &= (Eq + Iq - Hq == 2) and ((Eq - Hq == 2) if Iq == 0 else (Eq - Hq != 2))
p2 &= (E + I - H_ == 2)
if Eq + Iq - Hq != 2 or E + I - H_ != 2:
note.append("COUNT FAILS")
print(f" {e:7.4f} | {len(Q):3d} + {len(nonsym):<3d} | {E:2d} {H_:2d} {I:2d} {E-H_:3d} {E+I-H_:4d} |"
f" {Eq:2d} {Hq:2d} {Iq:2d} {Eq-Hq:4d} {Eq+Iq-Hq:6d} | {new:3d} / {missed:<3d} | {'; '.join(note)}")
ks = [(r["eps"], r["kappa"]) for b in brs for r in b["rows"] if r.get("kappa") == r.get("kappa") and "kappa" in r]
k10 = min((e for e, k in ks if abs(k - 4) > 0.4), default=None)
print(f" first I' (G period-doubles): {first_I} averaging's orbit set first wrong at eps = {avg_break}"
f" kappa first off by >10%: eps = {k10}")
print(f" P1 {'PASS' if p1 else 'FAIL'} P2 {'PASS' if p2 and pairs_ok and not symI else 'FAIL'}"
f" (symmetric orbit ever I: {symI}; nonsymmetric orbits paired: {pairs_ok}) counted eps: {counted}/{len(grid)}")
ok_all &= p1 and p2 and pairs_ok and not symI
print("\nALL PREDICTIONS HOLD" if ok_all else "\nA PREDICTION FAILED")
"""Symmetry-broken eps sweep: the model and code of sweep_hopf.py with G*x^N added to V1 (default G = 0.5, N = 3).
The odd term breaks P(q, p) = -(q, p), so no orbit is symmetric and M is no longer forced to be (DG)^2.
G = 0.5 was fixed before any orbit was computed: the largest of {0.5, 0.3, 0.2} for which the H = 1/2
surface stays star-shaped for all three cases up to eps = 3.
Changes from sweep_hopf.py, made here by patching that module so the shared code is unchanged:
- the added term in V1, its force, and both Jacobians;
- continuation uses full-period shooting (symmetric shooting has no solutions now);
- events: k crosses -4 (full-period doubling) or 0, located by bisection, and folds;
- targeted seeds push off each orbit along both directions transverse to the flow in the energy surface.
Controls (env SWEEP_POWER, SWEEP_G), predictions fixed before running:
N = 4, G = 0.1 even, so P survives: symmetric shooting, half-period classification and sweep_hopf's events and
search are kept, and the analytic average gains G*<x^4>. Judged by sweep_check.py's rules:
E'+I'-H' = 2; no symmetric orbit ever I; E+I-H = 2 with mirror pairs counted.
N = 5, G = 0.1 odd, like x^3. Judged by check_odd.py's rules: E+I-H = 2; at each full-period doubling the
parent itself turns I and no same-period orbit is spawned.
"""
import json, os, sys, time
import numpy as np
from multiprocessing import Pool
import check_hopf as h
import sweep_hopf as s
N = int(os.environ.get("SWEEP_POWER", 3))
G = float(os.environ.get("SWEEP_G", 0.5))
assert N % 2 == 1 or N == 4
class SysOdd(h.Sys):
def V1(self, X):
return super().V1(X) + G*X[0]**N
def rhs(self, X):
f = super().rhs(X)
f[2] = f[2] - h.EPS*N*G*X[0]**(N - 1)
return f
def jac(self, X):
J = super().jac(X)
J[2, 0] -= h.EPS*N*(N - 1)*G*X[0]**(N - 2)
return J
def Ab(self): # odd powers average to zero; <x^4> = (3/32)(S0 + S1)^2
A, b, const = super().Ab()
if N == 4:
A = A.copy(); b = b.copy()
A[0, 0] += 3*G/32*h.S0**2
b[0] += 3*G/16*h.S0**2
const += 3*G/32*h.S0**2
return A, b, const
_jac_v = s.jac_v
def jac_v(sy, X):
J = _jac_v(sy, X)
J[:, 2, 0] -= h.EPS*N*(N - 1)*G*X[0]**(N - 2)
return J
h.Sys = SysOdd # every Sys built in check_hopf / sweep_hopf gets the term
s.jac_v = jac_v
if N % 2:
s.shoot_sym = lambda sy, X, T: h.shoot(sy, X, T=T)
def events(name):
_, brs = s.continue_case(name, s.BASE)
sy = h.Sys(**s.CASES[name])
ev = []
for b in brs:
for thr, label in ((-4, "full-period doubling"), (0, "k crosses 0")):
for r0, r1 in zip(b["rows"], b["rows"][1:]):
if (r0["k"] - thr)*(r1["k"] - thr) < 0:
lo, hi, X, T = r0["eps"], r1["eps"], r0["X"], r0["T"]
for _ in range(30):
m, path = 0.5*(lo + hi), [(lo, X, T)]
if not s.advance(sy, path, m):
break
h.EPS = m
d = s.describe(sy, path[-1][1], path[-1][2])
if (d["k"] - thr)*(r0["k"] - thr) > 0:
lo, X, T = m, path[-1][1], path[-1][2]
else:
hi = m
ev.append((label, float(lo)))
if b["end"] is not None:
ev.append(("fold", b["end"]))
return name, ev
def search(args):
name, e, rows = args
h.EPS = e
sy = h.Sys(**s.CASES[name])
aug = []
for r in rows:
Q = np.linalg.qr(np.column_stack([sy.gradH(r["X"]), sy.rhs(r["X"]), np.eye(4)]))[0]
aug += [dict(r, v_pd=Q[:, 2]), dict(r, v_pd=Q[:, 3])]
return s.search((name, e, aug))
if __name__ == "__main__":
t0 = time.time()
names = list(s.CASES)
if N % 2 == 0: # symmetry kept: sweep_hopf's half-period events and seeds
events, search = s.events, s.search
with Pool(2) as pool:
ev = dict(pool.map(events, names))
print(f"events {ev} {time.time()-t0:.0f}s", flush=True)
grids = {}
for n in names:
extra = [x for k, e in ev[n] for x in
((e - 0.002, e + 0.002, e + 0.01, e + 0.03) if k != "fold" else (e - 0.002, e + 0.002))]
grids[n] = np.unique(np.concatenate([s.BASE, [x for x in extra if s.BASE[0] < x < s.BASE[-1]]]))
cont = dict(pool.starmap(s.continue_case, [(n, grids[n]) for n in names]))
print(f"continuation done {time.time()-t0:.0f}s", flush=True)
jobs = [(n, float(e), [b["rows"][i] for b in cont[n] if i < len(b["rows"])])
for n in names for i, e in enumerate(grids[n])]
found = pool.map(search, jobs, chunksize=1)
json.dump(s.tojson(dict(g=G, power=N, grids=grids, events=ev, cases=s.CASES, continuation=cont, search=found)),
open(sys.argv[1], "w"))
print(f"done {time.time()-t0:.0f}s", flush=True)
"""Evaluate the symmetry-broken sweep. Rules fixed before the run:
family orbits found by continuation or search with T < 1.6 x the shortest branch period at that eps
(a period-doubled orbit has T ~ 2T and falls outside; a warning prints if a branch itself would)
counted an eps counts only if no family orbit is 'unsure'
P1 E + I - H = 2 at every counted eps
P2 at every full-period doubling (a branch's k crosses -4), checked at eps* -/+ 0.002:
that branch is the one whose type changes between E and I, and the family has the same
number of orbits on both sides (no same-period orbit spawned)
falsified any counted eps with E + I - H != 2, or any doubling where P2's conditions fail
"""
import json, sys
import numpy as np
D = json.load(open(sys.argv[1]))
search = {(n, round(e, 12)): f for n, e, f in D["search"]}
def same(a, b):
return (abs(a["T"] - b["T"]) < 1e-6 and np.linalg.norm(np.subtract(a["center"], b["center"])) < 1e-5
and np.linalg.norm(np.subtract(a["mean"], b["mean"])) < 1e-5)
def family(name, i):
grid, brs = D["grids"][name], D["continuation"][name]
e = grid[i]
B = [b["rows"][i] for b in brs if i < len(b["rows"])]
F = list(B)
for f in search[(name, round(e, 12))]:
if not any(same(f, g) for g in F):
F.append(f)
tmin = min(r["T"] for r in B)
if max(r["T"] for r in B) >= 1.6*tmin:
print(f" WARNING {name} eps={e}: branch periods span more than 1.6x")
fam = [r for r in F if r["T"] < 1.6*tmin]
missed = sum(not any(same(r, f) for f in search[(name, round(e, 12))]) for r in B)
return B, fam, len(F) - len(fam), missed
p1_all = p2_all = True
for name in D["cases"]:
grid, brs = D["grids"][name], D["continuation"][name]
print(f"\n== {name} events: " + (", ".join(f"{k} at eps={e:.4f}" for k, e in D["events"][name]) or "none"))
print(" eps | branches +new | E H I | E+I-H | outside window | search missed | notes")
p1, counted = True, 0
for i, e in enumerate(grid):
B, fam, outside, missed = family(name, i)
unsure = [r for r in fam if r["full"] == "unsure"]
E, H_, I = (sum(r["full"] == t for r in fam) for t in "EHI")
note = ""
if unsure:
note = f"{len(unsure)} unsure, not counted"
else:
counted += 1
if E + I - H_ != 2:
p1, note = False, "COUNT FAILS"
print(f" {e:7.4f} | {len(B):3d} +{len(fam)-len(B):<3d} | {E:2d} {H_:2d} {I:2d} | {E+I-H_:4d} | {outside:6d} | {missed:6d} | {note}")
print(f" P1 {'PASS' if p1 else 'FAIL'} (counted {counted}/{len(grid)})")
p1_all &= p1
for kind, es in D["events"][name]:
if kind != "full-period doubling":
continue
i0 = int(np.argmin(np.abs(np.array(grid) - (es - 0.002))))
i1 = int(np.argmin(np.abs(np.array(grid) - (es + 0.002))))
parents = [b for b in brs if i1 < len(b["rows"]) and (b["rows"][i0]["k"] + 4)*(b["rows"][i1]["k"] + 4) < 0]
_, f0, _, _ = family(name, i0)
_, f1, _, _ = family(name, i1)
ok = len(parents) == 1 and {parents[0]["rows"][i0]["full"], parents[0]["rows"][i1]["full"]} == {"E", "I"} \
and len(f0) == len(f1)
if ok:
pr = parents[0]["rows"]
print(f" P2 doubling at eps*={es:.5f}: parent ({parents[0]['kind']} branch) {pr[i0]['full']} -> {pr[i1]['full']},"
f" k {pr[i0]['k']:+.4f} -> {pr[i1]['k']:+.4f}; family size {len(f0)} -> {len(f1)} PASS")
else:
print(f" P2 doubling at eps*={es:.5f}: parents {len(parents)}, family {len(f0)} -> {len(f1)} FAIL")
p2_all &= ok
n_dbl = sum(k == "full-period doubling" for n in D["cases"] for k, _ in D["events"][n])
print(f"\nP1 {'HOLDS' if p1_all else 'FALSIFIED'}; P2 " +
("UNTESTED (no full-period doubling occurred)" if n_dbl == 0 else ("HOLDS" if p2_all else "FALSIFIED")) +
f" over {n_dbl} doubling(s)")
"""Generic periodic-orbit toolkit for two-degree-of-freedom Hamiltonian flows (state X in R^4).
The methods of check_hopf.py / sweep_hopf.py, made system-agnostic for torus_exp.py and res12_exp.py:
Newton shooting with energy and phase conditions; Floquet type from k = tr M - 4 at two tolerances;
detection of double covers (orbit closes at T/2) with the prime orbit's type; batched shooting for blind
searches; secant continuation with step halving; bisection of events; a parallel sweep driver.
An experiment module provides: make(p) -> system, E(p) -> energy, seeds(p, sy) -> [(X, T guess)], TWIN = (lo, hi).
A system provides: rhs(X), jacv(X[4, n]) -> [n, 4, 4], H(X), gradH(X), wrap(dX), canon(mean), ok(X).
"""
import importlib
import numpy as np
from multiprocessing import Pool
from scipy.integrate import solve_ivp
def stm(sy, X, T, rtol=1e-11):
def f(t, Y):
return np.concatenate([sy.rhs(Y[:4]), (sy.jacv(Y[:4, None])[0] @ Y[4:].reshape(4, 4)).ravel()])
s = solve_ivp(f, (0, T), np.concatenate([X, np.eye(4).ravel()]), method="DOP853",
rtol=rtol, atol=rtol/10, t_eval=[T/2, T])
return [(s.y[:4, i], s.y[4:, i].reshape(4, 4)) for i in (0, 1)]
def shoot(sy, X, T, E, maxit=40, maxstep=0.1):
"""Newton on X(T) - X = 0 with H(X) = E and a phase condition; returns (X, T) or None."""
X = np.array(X, float)
for _ in range(maxit):
XT, M = stm(sy, X, T)[1]
F = np.concatenate([sy.wrap(XT - X), [sy.H(X) - E], [0.0]])
if np.linalg.norm(F) < 1e-10:
return X, T
J = np.zeros((6, 5))
J[:4, :4], J[:4, 4] = M - np.eye(4), sy.rhs(XT)
J[4, :4], J[5, :4] = sy.gradH(X), sy.rhs(X)
dz = np.linalg.lstsq(J, -F, rcond=None)[0]
dz *= min(1, maxstep/max(np.linalg.norm(dz), 1e-300))
X, T = X + dz[:4], T + dz[4]
return None
def kind(k, dk):
if min(abs(k), abs(k + 4)) < max(1e3*dk, 1e-8):
return "unsure"
return "E" if -4 < k < 0 else "H" if k > 0 else "I"
def signature(sy, X, T):
"""Phase-independent orbit identity: period, mean state, and width (centred second moment). The width is what
separates a newborn period-doubled orbit from the double cover it was born from: their periods and means agree
to ~1e-6 near the bifurcation, their widths do not."""
P = batch_traj(sy, X[None], np.array([T]), 128)[0]
return dict(T=float(T), mean=sy.canon(P.mean(1)), var=((P - P.mean(1, keepdims=True))**2).mean(1))
def same(sy, a, b):
return (abs(a["T"] - b["T"]) < 1e-6 and np.linalg.norm(sy.wrap(np.subtract(a["mean"], b["mean"]))) < 1e-5
and np.linalg.norm(np.subtract(a["var"], b["var"])) < 1e-7)
def describe(sy, X, T, E):
(Xh, Mh), (XT, M) = stm(sy, X, T)
(_, Mh2), (_, M2) = stm(sy, X, T, rtol=1e-13)
k, k2 = np.trace(M) - 4, np.trace(M2) - 4
d = dict(signature(sy, X, T), X=X, res=float(np.linalg.norm(sy.wrap(XT - X))), dH=float(abs(sy.H(X) - E)),
k=float(k), type=kind(k, abs(k - k2)), half_res=float(np.linalg.norm(sy.wrap(Xh - X))))
if d["half_res"] < 1e-7: # double cover of a period-T/2 orbit
kh, kh2 = np.trace(Mh) - 4, np.trace(Mh2) - 4
d.update(kh=float(kh), prime=kind(kh, abs(kh - kh2)))
d["good"] = d["prime"] != "I" # even covers of negative-hyperbolic orbits are bad
return d
# ---- batched integration in normalized time tau = t/T
def batch_traj(sy, X, T, m):
n = len(T)
f = lambda tau, Y: (sy.rhs(Y.reshape(n, 4).T).T*T[:, None]).ravel()
s = solve_ivp(f, (0, 1), np.asarray(X, float).ravel(), method="DOP853", rtol=1e-11, atol=1e-12,
t_eval=np.linspace(0, 1, m, endpoint=False))
return s.y.reshape(n, 4, m)
def batch_flow(sy, X, T):
n = len(T)
def f(tau, Y):
Y = Y.reshape(n, 20)
Xs = Y[:, :4].T
dX = sy.rhs(Xs).T*T[:, None]
dP = T[:, None, None]*(sy.jacv(Xs) @ Y[:, 4:].reshape(n, 4, 4))
return np.concatenate([dX, dP.reshape(n, 16)], 1).ravel()
Y0 = np.concatenate([X, np.tile(np.eye(4).ravel(), (n, 1))], 1).ravel()
Y = solve_ivp(f, (0, 1), Y0, method="DOP853", rtol=1e-11, atol=1e-12).y[:, -1].reshape(n, 20)
return Y[:, :4], Y[:, 4:].reshape(n, 4, 4)
def batch_shoot(sy, E, X, T, maxit=40, maxstep=0.1):
X, T, T0 = np.array(X, float), np.array(T, float), np.array(T, float)
act, done = np.ones(len(T), bool), np.zeros(len(T), bool)
for it in range(maxit):
idx = np.where(act)[0]
if not len(idx):
break
XT, M = batch_flow(sy, X[idx], T[idx])
for j, i in enumerate(idx):
F = np.concatenate([sy.wrap(XT[j] - X[i]), [sy.H(X[i]) - E], [0.0]])
nf = np.linalg.norm(F)
if nf < 1e-9:
act[i], done[i] = False, True
continue
if not np.isfinite(nf) or (it >= 15 and nf > 1e-3) or not sy.ok(X[i]) or abs(T[i] - T0[i]) > 0.5*T0[i]:
act[i] = False
continue
J = np.zeros((6, 5))
J[:4, :4], J[:4, 4] = M[j] - np.eye(4), sy.rhs(XT[j])
J[4, :4], J[5, :4] = sy.gradH(X[i]), sy.rhs(X[i])
dz = np.linalg.lstsq(J, -F, rcond=None)[0]
dz *= min(1, maxstep/max(np.linalg.norm(dz), 1e-300))
X[i], T[i] = X[i] + dz[:4], T[i] + dz[4]
return X[done], T[done]
def search(exp, p, rows):
"""Blind seeds from the experiment, plus seeds pushed off each known orbit along both transverse directions."""
sy, E = exp.make(p), exp.E(p)
seeds = list(exp.seeds(p, sy))
for r in rows:
Q = np.linalg.qr(np.column_stack([sy.gradH(r["X"]), sy.rhs(r["X"]), np.eye(4)]))[0]
for v in (Q[:, 2], Q[:, 3]):
for dl in (0.003, 0.01, 0.03, 0.1):
seeds += [(r["X"] + dl*v, r["T"]), (r["X"] - dl*v, r["T"])]
Xc, Tc = batch_shoot(sy, E, [s[0] for s in seeds], [s[1] for s in seeds])
keep = (Tc > exp.TWIN[0]) & (Tc < exp.TWIN[1])
Xc, Tc = Xc[keep], Tc[keep]
found = []
if len(Tc):
P = batch_traj(sy, Xc, Tc, 128)
for X, T, Pi in zip(Xc, Tc, P):
s = dict(T=float(T), mean=sy.canon(Pi.mean(1)), var=((Pi - Pi.mean(1, keepdims=True))**2).mean(1), X=X)
if not any(abs(s["T"] - f["T"]) < 1e-6 and np.linalg.norm(sy.wrap(s["mean"] - f["mean"])) < 1e-6
and np.linalg.norm(s["var"] - f["var"]) < 1e-7 for f in found):
found.append(s)
out = []
for f in found:
o = shoot(sy, f["X"], f["T"], E)
if o is not None:
d = describe(sy, o[0], o[1], E)
if exp.TWIN[0] < d["T"] < exp.TWIN[1] and not any(same(sy, d, g) for g in out):
out.append(d)
return out
# ---- continuation in the experiment parameter p
def advance(exp, path, target):
while path[-1][0] != target:
p0, X0, T0 = path[-1]
step = target - p0
while True:
p1 = target if abs(step) >= abs(target - p0) else p0 + step
if len(path) >= 2:
pa, Xa, Ta = path[-2]
w = (p1 - p0)/(p0 - pa)
Xg, Tg = X0 + w*(X0 - Xa), T0 + w*(T0 - Ta)
else:
Xg, Tg = X0, T0
o = shoot(exp.make(p1), Xg, Tg, exp.E(p1), maxit=15) # from a good predictor Newton needs < 10
if o is not None and np.linalg.norm(o[0] - Xg) < 0.05 and abs(o[1] - Tg) < 0.05:
path.append((p1, o[0], o[1]))
del path[:-2]
break
step /= 2
if abs(step) < 1e-5:
return False
return True
def branch_job(modname, grid, X, T, find_events):
exp = importlib.import_module(modname)
path, rows, end = [(grid[0], X, T)], [], None
for p in grid:
if not advance(exp, path, p):
end = float(path[-1][0])
break
rows.append(dict(p=float(p), **describe(exp.make(p), path[-1][1], path[-1][2], exp.E(p))))
ev = []
if find_events:
for key, thr, label in (("k", -4, "full-period doubling"), ("k", 0, "k crosses 0"),
("kh", -4, "prime orbit doubles")):
for r0, r1 in zip(rows, rows[1:]):
if key in r0 and key in r1 and (r0[key] - thr)*(r1[key] - thr) < 0:
lo, hi, Xl, Tl = r0["p"], r1["p"], r0["X"], r0["T"]
for _ in range(30):
m, pth = 0.5*(lo + hi), [(lo, Xl, Tl)]
if not advance(exp, pth, m):
break
d = describe(exp.make(m), pth[-1][1], pth[-1][2], exp.E(m))
if key not in d:
break
if (d[key] - thr)*(r0[key] - thr) > 0:
lo, Xl, Tl = m, pth[-1][1], pth[-1][2]
else:
hi = m
ev.append((label, float(lo)))
if end is not None:
ev.append(("fold", end))
return rows, end, ev
def search_job(modname, p, rows):
return search(importlib.import_module(modname), p, rows)
def tojson(o):
if isinstance(o, dict):
return {k: tojson(v) for k, v in o.items()}
if isinstance(o, (list, tuple)):
return [tojson(v) for v in o]
if isinstance(o, np.ndarray):
return o.tolist()
if isinstance(o, (np.floating, np.integer, np.bool_)):
return o.item()
return o
def run(modname, base, continuation=True):
"""Find orbits at base[0], continue them across base, locate events, refine the grid, re-continue, search everywhere.
With continuation=False every point is searched independently (seeded by the experiment's predictions)."""
exp = importlib.import_module(modname)
if not continuation:
with Pool(2) as pool:
found = pool.starmap(search_job, [(modname, float(p), []) for p in base], chunksize=1)
return tojson(dict(grid=base, events=[], branches=[], search=found))
init = search(exp, base[0], [])
with Pool(2) as pool:
first = pool.starmap(branch_job, [(modname, base, d["X"], d["T"], True) for d in init])
ev = sorted({e for _, _, es in first for e in es}, key=lambda x: x[1])
lo, hi = min(base[0], base[-1]), max(base[0], base[-1])
extra = [x for _, e in ev for x in (e - 0.002, e + 0.002, e + 0.01) if lo < x < hi]
grid = np.unique(np.concatenate([base, extra]))
if base[0] > base[-1]:
grid = grid[::-1]
cont = pool.starmap(branch_job, [(modname, grid, d["X"], d["T"], False) for d in init])
jobs = [(modname, float(p), [b[0][i] for b in cont if i < len(b[0])]) for i, p in enumerate(grid)]
found = pool.starmap(search_job, jobs, chunksize=1)
return tojson(dict(grid=grid, events=ev, branches=[dict(rows=r, end=e) for r, e, _ in cont], search=found))
"""Charged particle on the flat torus [0, 2pi)^2 in a uniform magnetic field B = 1, perturbed by eps*V(q).
dq/dt = v, dv/dt = (v_y, -v_x) - eps*grad V(q), energy |v|^2/2 + eps*V(q) = rho^2/2.
Unperturbed, every orbit is a cyclotron circle of radius rho and period 2pi, q = c + rho(sin t, cos t), so orbits
correspond to guiding centres c = (x + v_y, y - v_x) in T^2. Averaging V over the circle multiplies each Fourier mode
by J0(|k| rho): <V>(c) = sum A J0(|k| rho) cos(k.c + phase). Centres drift along its contours, so extrema of <V>
give elliptic orbits and saddles hyperbolic ones; chi(T^2) = 0.
Predictions, fixed before any orbit was computed:
T1 eps = 0.01, every rho: the full-system orbits match the critical points of <V> one to one, with types;
E - H = 0; at least 4 orbits
T2 rho sweep 0.3 -> 3.0 at eps = 0.01: the number of orbits changes, E + I - H = 0 at every counted rho
T3 eps sweep at rho = 1.6, up to 0.9 of the eps where the speed could reach zero: E + I - H = 0 at every counted eps
"""
import numpy as np
from scipy.special import j0
MODES = [(1, 0, 1.0, 0.0), (0, 1, 0.7, 0.5), (1, 1, 0.5, 1.3), (2, -1, 0.3, 0.4), (1, -2, 0.25, 2.1)]
K = np.array([(m[0], m[1]) for m in MODES], float)
A = np.array([m[2] for m in MODES])
PH = np.array([m[3] for m in MODES])
TWO_PI = 2*np.pi
TWIN = (0.7*TWO_PI, 1.3*TWO_PI)
SWEEP = "rho" # "rho": p = rho at EPS0; "eps": p = eps at RHO_E
EPS0, RHO_E = 0.01, 1.6
def V(x, y):
return sum(a*np.cos(kx*x + ky*y + ph) for (kx, ky), a, ph in zip(K, A, PH))
def dV(x, y):
s = [-a*np.sin(kx*x + ky*y + ph) for (kx, ky), a, ph in zip(K, A, PH)]
return sum(si*kx for si, (kx, _) in zip(s, K)), sum(si*ky for si, (_, ky) in zip(s, K))
def ddV(x, y):
c = [-a*np.cos(kx*x + ky*y + ph) for (kx, ky), a, ph in zip(K, A, PH)]
return (sum(ci*kx*kx for ci, (kx, ky) in zip(c, K)), sum(ci*kx*ky for ci, (kx, ky) in zip(c, K)),
sum(ci*ky*ky for ci, (kx, ky) in zip(c, K)))
VMAX = float(np.max(V(*np.meshgrid(np.linspace(0, TWO_PI, 801), np.linspace(0, TWO_PI, 801)))))
EPS_MAX = 0.9*(RHO_E**2/2)/VMAX
class Torus:
def __init__(self, eps):
self.eps = eps
def rhs(self, X):
x, y, vx, vy = X
Vx, Vy = dV(x, y)
return np.array([vx, vy, vy - self.eps*Vx, -vx - self.eps*Vy])
def jacv(self, X):
x, y = X[0], X[1]
Vxx, Vxy, Vyy = ddV(x, y)
J = np.zeros((X.shape[1], 4, 4))
J[:, 0, 2] = J[:, 1, 3] = J[:, 2, 3] = 1
J[:, 3, 2] = -1
J[:, 2, 0], J[:, 2, 1] = -self.eps*Vxx, -self.eps*Vxy
J[:, 3, 0], J[:, 3, 1] = -self.eps*Vxy, -self.eps*Vyy
return J
def H(self, X):
return 0.5*(X[2]**2 + X[3]**2) + self.eps*V(X[0], X[1])
def gradH(self, X):
Vx, Vy = dV(X[0], X[1])
return np.array([self.eps*Vx, self.eps*Vy, X[2], X[3]])
def wrap(self, dX):
dX = np.array(dX, float)
dX[..., :2] = (dX[..., :2] + np.pi) % TWO_PI - np.pi
return dX
def canon(self, m):
m = np.array(m, float)
m[:2] %= TWO_PI
return m
def ok(self, X):
return np.abs(X[2:]).max() < 10
def rho_eps(p):
return (p, EPS0) if SWEEP == "rho" else (RHO_E, p)
def make(p):
return Torus(rho_eps(p)[1])
def E(p):
return rho_eps(p)[0]**2/2
# ---- averaged potential on T^2 and its critical points
def avg(c, rho, d=0):
"""<V>, its gradient (d=1) or Hessian (d=2) at guiding centre c = (cx, cy)."""
w = A*j0(np.linalg.norm(K, axis=1)*rho)
ph = np.asarray(c) @ K.T + PH
if d == 0:
return np.sum(w*np.cos(ph), axis=-1)
if d == 1:
s = -w*np.sin(ph)
return np.stack([np.sum(s*K[:, 0], -1), np.sum(s*K[:, 1], -1)])
co = -w*np.cos(ph)
return np.array([[np.sum(co*K[:, 0]*K[:, 0]), np.sum(co*K[:, 0]*K[:, 1])],
[np.sum(co*K[:, 0]*K[:, 1]), np.sum(co*K[:, 1]*K[:, 1])]])
def avg_critical(rho, n=240):
"""Critical points of <V> by Newton from every grid cell whose gradient winding is nonzero and from every
local minimum of |grad|^2; checked against the windings (a degree count over all cells)."""
g = np.linspace(0, TWO_PI, n, endpoint=False)
cx, cy = np.meshgrid(g, g, indexing="ij")
C = np.stack([cx, cy], -1)
G = avg(C, rho, 1) # [2, n, n]
ang = np.arctan2(G[1], G[0])
def dang(a, b):
return (b - a + np.pi) % TWO_PI - np.pi
a00, a10, a11, a01 = ang, np.roll(ang, -1, 0), np.roll(np.roll(ang, -1, 0), -1, 1), np.roll(ang, -1, 1)
wind = np.rint((dang(a00, a10) + dang(a10, a11) + dang(a11, a01) + dang(a01, a00))/TWO_PI).astype(int)
g2 = G[0]**2 + G[1]**2
locmin = np.ones_like(g2, bool)
for sx in (-1, 0, 1):
for sy in (-1, 0, 1):
if sx or sy:
locmin &= g2 <= np.roll(np.roll(g2, sx, 0), sy, 1)
h = TWO_PI/n
starts = [C[i, j] + h/2 for i, j in zip(*np.nonzero(wind))] + [C[i, j] for i, j in zip(*np.nonzero(locmin))]
pts = []
for c in starts:
c = np.array(c, float)
for _ in range(50):
gr = avg(c, rho, 1)
if np.linalg.norm(gr) < 1e-13:
break
try:
c = c - np.linalg.solve(avg(c, rho, 2), gr)
except np.linalg.LinAlgError:
break
if np.linalg.norm(avg(c, rho, 1)) > 1e-11:
continue
c %= TWO_PI
if any(np.linalg.norm((c - q["c"] + np.pi) % TWO_PI - np.pi) < 1e-7 for q in pts):
continue
Hs = avg(c, rho, 2)
dt, tr = np.linalg.det(Hs), np.trace(Hs)
pts.append(dict(c=c, kind="saddle" if dt < 0 else ("min" if tr > 0 else "max"), det=float(dt)))
# degree check: every cell with nonzero winding holds points whose indices sum to that winding
idx = np.zeros_like(wind)
for q in pts:
i, j = (np.floor(q["c"]/h).astype(int)) % n
idx[i, j] += -1 if q["kind"] == "saddle" else 1
return pts, bool(np.array_equal(idx, wind)), int(wind.sum())
def seeds(p, sy):
rho, eps = rho_eps(p)
cs = [q["c"] for q in avg_critical(rho)[0]] + \
[np.array([i, j])*TWO_PI/20 for i in range(20) for j in range(20)]
out = []
for c in cs:
x, y = c[0], c[1] + rho
sp2 = 2*(E(p) - eps*V(x, y))
if sp2 > 0:
out.append((np.array([x, y, np.sqrt(sp2), 0.0]), TWO_PI))
return out
"""1:2 resonance: H = (px^2 + x^2)/2 + (py^2 + 4 y^2)/2 + eps*V1,
V1 = kappa x^2 y + d (x^2 - y^2) + x^4 + y^4 + c x^2 y^2 + g x^3, d = 0.2, c = 1, g = 0.5, energy E = 1/2.
Unperturbed, every orbit has period 2pi except the y-mode (x = px = 0), which has period pi. The orbit space is a
sphere with one cone point, the y-mode. x^2 y is the resonant (Fermi) term; g x^3 averages to zero at first order
but breaks x -> -x, which would otherwise make every orbit's monodromy a square, as P did for the 1:1 system.
Count rule (Frauenfelder-Koh-Moreno): count good orbits of period ~2pi with E + I - H. The y-mode appears there as
its double cover, which is bad (left out) when the y-mode itself is negative-hyperbolic.
Averaging near the cone point (w = x-amplitude * e^{i alpha}, the y-mode at amplitude b0, b0^2 = E/2):
<V1> ~ const + |w|^2 (D + R cos 2 alpha), D = 5d/8 - 3 b0^2/16 + c b0^2/4, R = kappa b0/4
-> cone point an extremum (y-mode elliptic) for |R| < |D|, a saddle (y-mode negative-hyperbolic: its
parametric-resonance tongue) for |R| > |D|; kappa* = 4|D|/b0 = 1.125.
Predictions, fixed before any orbit was computed:
R1 kappa sweep 0.05 -> 2.5 at eps = 0.01, and eps sweep at kappa = 2 (up to 0.9 of the star-shape limit, max 1):
good-orbit count E + I - H = 2 at every counted point
R2 near kappa* the y-mode turns elliptic -> negative-hyperbolic and one new orbit of period ~2pi appears;
the naive count (bad double cover counted as hyperbolic) drops from 2 to 1 there
"""
import numpy as np
from scipy.optimize import brentq
D_, C_, G_, EN = 0.2, 1.0, 0.5, 0.5
TWIN = (0.7*2*np.pi, 1.3*2*np.pi)
SWEEP = "kappa" # "kappa": p = kappa at EPS0; "eps": p = eps at KAPPA_E
EPS0, KAPPA_E = 0.01, 2.0
B0 = np.sqrt(EN/2)
DCOEF = 5*D_/8 - 3*B0**2/16 + C_*B0**2/4
KAPPA_STAR = 4*abs(DCOEF)/B0
class Res12:
def __init__(self, eps, kappa):
self.eps, self.kappa = eps, kappa
def V1(self, x, y):
return self.kappa*x*x*y + D_*(x*x - y*y) + x**4 + y**4 + C_*x*x*y*y + G_*x**3
def rhs(self, X):
x, y, px, py = X
e, k = self.eps, self.kappa
return np.array([px, py,
-x - e*(2*k*x*y + 2*D_*x + 4*x**3 + 2*C_*x*y*y + 3*G_*x*x),
-4*y - e*(k*x*x - 2*D_*y + 4*y**3 + 2*C_*x*x*y)])
def jacv(self, X):
x, y = X[0], X[1]
e, k = self.eps, self.kappa
J = np.zeros((X.shape[1], 4, 4))
J[:, 0, 2] = J[:, 1, 3] = 1
J[:, 2, 0] = -1 - e*(2*k*y + 2*D_ + 12*x*x + 2*C_*y*y + 6*G_*x)
J[:, 2, 1] = J[:, 3, 0] = -e*(2*k*x + 4*C_*x*y)
J[:, 3, 1] = -4 - e*(-2*D_ + 12*y*y + 2*C_*x*x)
return J
def H(self, X):
x, y, px, py = X
return 0.5*(px*px + x*x) + 0.5*(py*py + 4*y*y) + self.eps*self.V1(x, y)
def gradH(self, X):
f = self.rhs(X)
return np.array([-f[2], -f[3], f[0], f[1]])
def wrap(self, dX):
return np.asarray(dX, float)
def canon(self, m):
return np.asarray(m, float)
def ok(self, X):
return np.abs(X).max() < 5
def eps_kappa(p):
return (EPS0, p) if SWEEP == "kappa" else (p, KAPPA_E)
def make(p):
return Res12(*eps_kappa(p))
def E(p):
return EN
def star_limit(kappa, eps_grid):
"""Largest eps in eps_grid for which every ray from the origin crosses H = E once (energy surface star-shaped)."""
rng = np.random.default_rng(0)
n = rng.normal(size=(200000, 4)); n /= np.linalg.norm(n, axis=1, keepdims=True)
nx, ny, npp = n[:, 0], n[:, 1], (n[:, 2:]**2).sum(1)
last = None
for e in eps_grid:
A2 = npp/2 + (nx**2 + 4*ny**2)/2 + e*D_*(nx**2 - ny**2)
B3 = e*(kappa*nx**2*ny + G_*nx**3)
C4 = e*(nx**4 + ny**4 + C_*nx**2*ny**2)
disc = 9*B3**2 - 32*A2*C4
two = (A2 > 0) & (disc > 0) & (B3 < 0)
sq = np.sqrt(np.where(two, disc, 0))
t1, t2 = (-3*B3 - sq)/(8*C4), (-3*B3 + sq)/(8*C4)
f = lambda t: A2*t**2 + B3*t**3 + C4*t**4
if (two & (f(t1) > EN) & (f(t2) < EN)).any() or (A2 <= 0).any():
break
last = e
return last
def avg_f(kappa, P, Q):
"""<V1> in the pole chart (P, Q) = b (cos psi, sin psi), x = sqrt(s) cos t, y = b cos(2t - psi), s = 2E - 4b^2."""
b2 = P*P + Q*Q
s = 2*EN - 4*b2
return kappa*s*P/4 + D_*(s - b2)/2 + 3*(s*s + b2*b2)/8 + C_*s*b2/4
def avg_critical(kappa):
"""Regular critical points of <V1>. For kappa != 0, df/dQ = 0 forces Q = 0, so they are the roots of
d/dP f(P, 0) on (-b0, b0); type from the 2x2 Hessian. The cone point (y-mode) is always critical:
extremum if |R| < |D|, saddle if |R| > |D|."""
h = 1e-6
g = lambda P: (avg_f(kappa, P + h, 0) - avg_f(kappa, P - h, 0))/(2*h)
Ps = np.linspace(-B0, B0, 4001)[1:-1]
gv = np.array([g(P) for P in Ps])
pts = []
for i in np.nonzero(np.sign(gv[:-1]) != np.sign(gv[1:]))[0]:
P = brentq(g, Ps[i], Ps[i + 1], xtol=1e-14)
fpp = (avg_f(kappa, P + 1e-4, 0) - 2*avg_f(kappa, P, 0) + avg_f(kappa, P - 1e-4, 0))/1e-8
fqq = (avg_f(kappa, P, 1e-4) - 2*avg_f(kappa, P, 0) + avg_f(kappa, P, -1e-4))/1e-8
kind = "saddle" if fpp*fqq < 0 else ("min" if fpp > 0 else "max")
pts.append(dict(s=float(2*EN - 4*P*P), psi=0.0 if P >= 0 else float(np.pi), kind=kind))
R = kappa*B0/4
cone = "saddle" if abs(R) > abs(DCOEF) else ("min" if DCOEF > 0 else "max")
return pts, cone
def seeds(p, sy):
out = []
pred = avg_critical(eps_kappa(p)[1])[0]
for q in pred: # predicted orbits first
a, b = np.sqrt(q["s"]), np.sqrt(max(EN - q["s"]/2, 0)/2)
X = np.array([a, b*np.cos(-q["psi"]), 0.0, 2*b*np.sin(q["psi"])])
out.append((brentq(lambda k: sy.H(k*X) - EN, 0.3, 2.0)*X, 2*np.pi))
ends = [2*EN*f for k in (1, 2, 3) for f in (10.0**-k, 1 - 10.0**-k)] # rings near the y-mode and the pole
for s in np.unique(np.concatenate([np.linspace(0, 2*EN, 21), ends])): # s = x-amplitude^2; 0 = y-mode, 2E = pole
for psi in np.linspace(0, 2*np.pi, 20, endpoint=False) if 0 < s < 2*EN else [0.0]:
a, b = np.sqrt(s), np.sqrt(max(EN - s/2, 0)/2)
X = np.array([a, b*np.cos(-psi), 0.0, 2*b*np.sin(psi)])
try:
X = brentq(lambda k: sy.H(k*X) - EN, 0.3, 2.0)*X
except ValueError:
continue
out.append((X, 2*np.pi))
return out
"""Run the four sweeps (torus rho, torus eps, 1:2 kappa, 1:2 eps); one JSON file per sweep."""
import json, time
import numpy as np
import orbits, torus_exp, res12_exp
t0 = time.time()
plan = [("torus_rho", torus_exp, "rho", np.linspace(0.3, 3.0, 55), False), # predicted seeds at every rho; folds
("torus_eps", torus_exp, "eps", np.geomspace(0.01, torus_exp.EPS_MAX, 30), True), # make continuation useless
("res12_kappa", res12_exp, "kappa", np.linspace(0.05, 2.5, 50), True)]
lim = res12_exp.star_limit(res12_exp.KAPPA_E, np.linspace(0.01, 1.0, 100))
plan.append(("res12_eps", res12_exp, "eps", np.geomspace(0.01, min(1.0, 0.9*lim), 30), True))
for name, mod, sweep, base, cont in plan:
mod.SWEEP = sweep
res = orbits.run(mod.__name__, base, continuation=cont)
json.dump(res, open(f"exp_{name}.json", "w"))
print(f"{name}: {len(res['grid'])} points, events {res['events']} ({time.time()-t0:.0f}s)", flush=True)
print("done", flush=True)
"""Evaluate the torus and 1:2 sweeps. Rules fixed before looking at the results:
family continued orbits plus search finds at each point (deduplicated), period inside TWIN
counted a point counts only if no family orbit (or prime orbit of a double cover) is 'unsure'
T1 torus rho sweep (eps = 0.01): full orbits match <V>'s critical points one to one (guiding centre within 0.05,
extremum <-> E, saddle <-> H); E - H = 0; at least 4 orbits
T2 torus rho sweep: E + I - H = 0 at every counted rho; the number of orbits changes along the sweep
T3 torus eps sweep: E + I - H = 0 at every counted eps
R1 1:2 sweeps: good-orbit E + I - H = 2 at every counted point
R2 1:2 kappa sweep: the y-mode turns E -> I near kappa* = 1.125 (prime k crosses -4), the family gains one
orbit across that crossing, and the naive count goes 2 -> 1
"""
import json, sys
import numpy as np
import orbits, torus_exp, res12_exp
TW = torus_exp.Torus(0.0)
RW = res12_exp.Res12(0.0, 0.0)
def families(R, sy):
out = []
for i, p in enumerate(R["grid"]):
F = [b["rows"][i] for b in R["branches"] if i < len(b["rows"])]
for f in R["search"][i]:
if not any(orbits.same(sy, f, g) for g in F):
F.append(f)
out.append((p, F))
return out
def counts(F, good_only=False):
S = [f for f in F if f.get("good", True)] if good_only else F
E, H, I = (sum(f["type"] == t for f in S) for t in "EHI")
return E, H, I, E + I - H
def unsure(F):
return any(f["type"] == "unsure" or f.get("prime") == "unsure" for f in F)
summary, verdict = {}, {}
# ---------------- torus
for name, target in (("torus_rho", 0), ("torus_eps", 0)):
R = json.load(open(f"exp_{name}.json"))
fams, ok, counted, rows = families(R, TW), True, 0, []
t1 = True
for p, F in fams:
E, H, I, c = counts(F)
u = unsure(F)
rec = dict(p=p, n=len(F), E=E, H=H, I=I, count=c, unsure=u)
if not u:
counted += 1
ok &= c == target
if name == "torus_rho":
torus_exp.SWEEP = "rho"
pts, deg_ok, _ = torus_exp.avg_critical(p)
used, match = set(), len(pts) == len(F)
for f in F:
d = [np.linalg.norm(TW.wrap(np.r_[f["mean"][:2] - q["c"], 0, 0]))[()] for q in pts]
j = int(np.argmin(d)) if d else -1
if j < 0 or d[j] > 0.05 or j in used or ((f["type"] == "E") != (pts[j]["kind"] != "saddle")):
match = False
used.add(j)
rec.update(n_avg=len(pts), match=bool(match), degree_ok=deg_ok,
avg=[dict(c=q["c"].tolist(), kind=q["kind"]) for q in pts])
t1 &= match and len(F) >= 4 and (E - H == 0)
rec["orbits"] = [dict(T=f["T"], type=f["type"], k=f["k"], c=f["mean"][:2]) for f in F]
rows.append(rec)
ns = [r["n"] for r in rows]
print(f"\n== {name}: events {R['events'] or 'none'}")
print(f" points {len(rows)}, counted {counted}, orbits per point {min(ns)}-{max(ns)}, "
f"E+I-H values seen {sorted({r['count'] for r in rows if not r['unsure']})}")
if name == "torus_rho":
bad = [r["p"] for r in rows if not r["match"]]
print(f" T1 {'PASS' if t1 else 'FAIL'} (mismatches at rho = {[round(b, 4) for b in bad] or 'none'})")
print(f" T2 {'PASS' if ok and len(set(ns)) > 1 else 'FAIL'}")
verdict["T1"], verdict["T2"] = t1, ok and len(set(ns)) > 1
else:
print(f" T3 {'PASS' if ok else 'FAIL'}")
verdict["T3"] = ok
summary[name] = dict(events=R["events"], rows=rows,
branches=[[dict(p=r["p"], k=r["k"], type=r["type"]) for r in b["rows"]] for b in R["branches"]])
# ---------------- 1:2 resonance
r1 = True
for name in ("res12_kappa", "res12_eps"):
R = json.load(open(f"exp_{name}.json"))
fams, ok, counted, rows = families(R, RW), True, 0, []
for p, F in fams:
Eg, Hg, Ig, cg = counts(F, True)
En, Hn, In, cn = counts(F)
ym = [f for f in F if "prime" in f]
u = unsure(F)
if not u:
counted += 1
ok &= cg == 2
rows.append(dict(p=p, n=len(F), good=cg, naive=cn, E=Eg, H=Hg, I=Ig, unsure=u,
ymode=(ym[0]["prime"] if ym else None), ymode_kh=(ym[0]["kh"] if ym else None),
orbits=[dict(T=f["T"], type=f["type"], k=f["k"], prime=f.get("prime"), good=f.get("good", True))
for f in F]))
r1 &= ok
print(f"\n== {name}: events {R['events'] or 'none'}")
print(f" points {len(rows)}, counted {counted}, good counts seen {sorted({r['good'] for r in rows if not r['unsure']})},"
f" naive counts seen {sorted({r['naive'] for r in rows})}")
if name == "res12_kappa":
cross = [(a, b) for a, b in zip(rows, rows[1:]) if a["ymode"] == "E" and b["ymode"] == "I"]
if len(cross) == 1:
a, b = cross[0]
ev = [e for k, e in R["events"] if k == "prime orbit doubles"]
r2 = b["n"] == a["n"] + 1 and a["naive"] == 2 and b["naive"] == 1
print(f" y-mode E -> I between kappa {a['p']:.4f} and {b['p']:.4f} (bisected: {ev}); averaging said {res12_exp.KAPPA_STAR:.4f}")
print(f" family {a['n']} -> {b['n']}, naive {a['naive']} -> {b['naive']}, good {a['good']} -> {b['good']}")
else:
r2 = False
print(f" y-mode E -> I crossings found: {len(cross)}")
print(f" R2 {'PASS' if r2 else 'FAIL'}")
verdict["R2"] = r2
summary[name] = dict(events=R["events"], rows=rows,
branches=[[dict(p=r["p"], k=r["k"], type=r["type"], kh=r.get("kh")) for r in b["rows"]]
for b in R["branches"]])
verdict["R1"] = r1
print(f"\n R1 {'PASS' if r1 else 'FAIL'}")
print("\n" + " ".join(f"{k} {'PASS' if v else 'FAIL'}" for k, v in sorted(verdict.items())))
json.dump(orbits.tojson(dict(verdict=verdict, **summary)), open("experiments_summary.json", "w"))
"""Compact the result files into the data the page draws: page_data.json.
Floquet coordinate y (one axis for every orbit type): elliptic y = theta/pi in (0, 1) where the multiplier is e^{i theta};
hyperbolic y = -arccosh(1 + k/2)/pi < 0; inverse-hyperbolic y = 1 + arccosh(-1 - k/2)/pi > 1.
y = 0 is a multiplier at +1 (folds), y = 1 a multiplier at -1 (period doubling)."""
import json
import numpy as np
def fy(k):
if k > 0:
return -np.arccosh(1 + k/2)/np.pi
if k < -4:
return 1 + np.arccosh(-1 - k/2)/np.pi
return np.arccos(np.clip(1 + k/2, -1, 1))/np.pi
def r4(v):
return round(float(v), 4)
def same(a, b):
return (abs(a["T"] - b["T"]) < 1e-6 and np.linalg.norm(np.subtract(a["center"], b["center"])) < 1e-5
and np.linalg.norm(np.subtract(a["mean"], b["mean"])) < 1e-5)
def one_to_one(path):
D, out = json.load(open(path)), {}
search = {(n, round(e, 12)): f for n, e, f in D["search"]}
for name in D["cases"]:
grid, brs = D["grids"][name], D["continuation"][name]
tracks = [[r4(fy(r["k"])) for r in b["rows"]] for b in brs]
ttype = ["".join(r["full"][0] if r["full"] != "unsure" else "u" for r in b["rows"]) for b in brs]
extra, cnt, counted = [], {"E": [], "H": [], "I": []}, 0
for i, e in enumerate(grid):
F = [b["rows"][i] for b in brs if i < len(b["rows"])]
for f in search[(name, round(e, 12))]:
if not any(same(f, g) for g in F):
F.append(f)
extra.append([r4(e), r4(fy(f["k"])), f["full"][0]])
for t in "EHI":
cnt[t].append(sum(r["full"] == t for r in F))
counted += not any(r["full"] == "unsure" or r.get("quot") == "unsure" for r in F)
out[name] = dict(x=[r4(e) for e in grid], tracks=tracks, ttype=ttype, extra=extra, count=cnt, counted=counted,
events=[[k, r4(e)] for k, e in D["events"][name]],
kinds=[b["kind"] for b in brs])
return out
S = json.load(open("experiments_summary.json"))
def tracks_of(sw, key="k"):
return [[r4(fy(r[key])) for r in b if r.get(key) is not None] for b in sw["branches"]]
tr, te, rk, re_ = S["torus_rho"], S["torus_eps"], S["res12_kappa"], S["res12_eps"]
page = dict(
sym=one_to_one("sweep_results.json"),
odd=one_to_one("odd_results.json"),
torus_rho=dict(x=[r4(r["p"]) for r in tr["rows"]], n=[r["n"] for r in tr["rows"]],
E=[r["E"] for r in tr["rows"]], H=[r["H"] for r in tr["rows"]], I=[r["I"] for r in tr["rows"]],
match=[r["match"] for r in tr["rows"]],
orbits=[[[r4(o["c"][0]), r4(o["c"][1]), o["type"][0]] for o in r["orbits"]] for r in tr["rows"]],
avg=[[[r4(q["c"][0]), r4(q["c"][1]), q["kind"][0]] for q in r["avg"]] for r in tr["rows"]]),
torus_eps=dict(x=[r4(r["p"]) for r in te["rows"]], E=[r["E"] for r in te["rows"]], H=[r["H"] for r in te["rows"]],
I=[r["I"] for r in te["rows"]], tracks=tracks_of(te),
tx=[[r4(r["p"]) for r in b] for b in te["branches"]],
ttype=["".join(r["type"][0] for r in b) for b in te["branches"]], events=te["events"]),
res12_kappa=dict(x=[r4(r["p"]) for r in rk["rows"]], n=[r["n"] for r in rk["rows"]],
good=[r["good"] for r in rk["rows"]], naive=[r["naive"] for r in rk["rows"]],
ymode=[r["ymode"] for r in rk["rows"]],
ymode_y=[r4(fy(r["ymode_kh"])) if r["ymode_kh"] is not None else None for r in rk["rows"]],
events=rk["events"]),
res12_eps=dict(x=[r4(r["p"]) for r in re_["rows"]], n=[r["n"] for r in re_["rows"]],
good=[r["good"] for r in re_["rows"]], naive=[r["naive"] for r in re_["rows"]],
ymode=[r["ymode"] for r in re_["rows"]],
ymode_y=[r4(fy(r["ymode_kh"])) if r["ymode_kh"] is not None else None for r in re_["rows"]],
tracks=tracks_of(re_), tx=[[r4(r["p"]) for r in b] for b in re_["branches"]],
ttype=["".join(r["type"][0] for r in b) for b in re_["branches"]], events=re_["events"]),
verdict=S["verdict"])
import subprocess, re, torus_exp, res12_exp
page["x4"], page["x5"] = one_to_one("x4_results.json"), one_to_one("x5_results.json")
out4 = subprocess.run(["python3", "sweep_check.py", "x4_results.json"], capture_output=True, text=True).stdout
out5 = subprocess.run(["python3", "check_odd.py", "x5_results.json"], capture_output=True, text=True).stdout
m5 = re.search(r"P1 (\w+);\s+P2 (.+?) over (\d+) doubling", out5)
page["controls"] = dict(x4=dict(held="ALL PREDICTIONS HOLD" in out4),
x5=dict(held=bool(m5) and m5.group(1) == "HOLDS" and m5.group(2).startswith("HOLDS"),
p1=bool(m5) and m5.group(1) == "HOLDS",
p2_untested=bool(m5) and m5.group(2).startswith("UNTESTED"), doublings=int(m5.group(3)) if m5 else 0))
page["torus_modes"] = [list(m) for m in torus_exp.MODES]
page["res12_const"] = dict(E=res12_exp.EN, d=res12_exp.D_, c=res12_exp.C_, kappa_star=r4(res12_exp.KAPPA_STAR))
page["counted"] = {k: sum(v["counted"] for v in page[k].values()) for k in ("sym", "odd", "x4", "x5")}
page["counted"].update({k + "_total": sum(len(v["x"]) for v in page[k].values()) for k in ("sym", "odd", "x4", "x5")})
for key in ("torus_rho", "torus_eps", "res12_kappa", "res12_eps"):
page["counted"][key] = sum(not r["unsure"] for r in S[key]["rows"])
page["counted"][key + "_total"] = len(S[key]["rows"])
json.dump(page, open("page_data.json", "w"), separators=(",", ":"))
print({k: len(json.dumps(v)) for k, v in page.items()})