Phase-space geometry in three steps¶
Family: manifolds · Script: examples/manifolds/phase_space_geometry.py
Obtained |
two-mode model: RPO “1” in the slice, T = 3.6415119771, θ = 0.0809697054; Lorenz homoclinic orbit as a connection, r = 13.9265574075 |
Reference |
Budanur, Borrero-Echeverry & Cvitanović (arXiv:1411.3303), Table I: 3.64151221; Sparrow 1982: 13.926 |
1. Symmetry reduction: the two-mode SO(2) model of Budanur, Borrero-Echeverry & Cvitanović
(arXiv:1411.3303) in the first Fourier mode slice; its relative periodic orbit "1" becomes a
periodic orbit of the reduced flow (Table I: T = 3.64151221, θ = 0.08096967, Λ = −1.48372354).
2. Poincaré section of the slice through the relative equilibrium, located on the declared dense
output (DOP853 continuous extension, Brent).
3. Lorenz: the homoclinic orbit of the origin found as a *connection* (unstable manifold by the
parametrization method, projection condition on a sphere around the origin), r ≈ 13.926.
Run: python examples/manifolds/phase_space_geometry.py
Run it: python examples/manifolds/phase_space_geometry.py (in the conatus environment).
1"""Phase-space geometry in three steps.
2
31. Symmetry reduction: the two-mode SO(2) model of Budanur, Borrero-Echeverry & Cvitanović
4 (arXiv:1411.3303) in the first Fourier mode slice; its relative periodic orbit "1" becomes a
5 periodic orbit of the reduced flow (Table I: T = 3.64151221, θ = 0.08096967, Λ = −1.48372354).
62. Poincaré section of the slice through the relative equilibrium, located on the declared dense
7 output (DOP853 continuous extension, Brent).
83. Lorenz: the homoclinic orbit of the origin found as a *connection* (unstable manifold by the
9 parametrization method, projection condition on a sphere around the origin), r ≈ 13.926.
10
11Run: python examples/manifolds/phase_space_geometry.py
12"""
13import numpy as np
14import scipy.optimize as so
15
16from conatus.expert import (FirstFourierModeSlice, PoincareSection, find_connection,
17 find_rpo_in_slice, ir_system, ode_problem, unstable_manifold)
18from conatus.bifurcation import DynamicalSystem
19from conatus.manifolds import ODEIntegrator, Trajectory, poincare_crossings, unstable_left_basis
20from conatus.slices import RotationOps, two_mode_callables
21
22
23def two_mode_rpo():
24 f, jac = two_mode_callables() # eq. (21) with the parameters of eq. (22)
25 ds = DynamicalSystem.from_callables(f, 4, {}, jac=jac, fd_x=1e-5, fd_p=1e-6, name="two-mode")
26 ops = RotationOps([1, 2]) # (z1, z2) -> (e^{iθ} z1, e^{2iθ} z2)
27 slc = FirstFourierModeSlice(ops, border_tol=1e-8) # template (1, 0, 0, 0): Im z1 = 0, Re z1 > 0
28 # relative equilibrium = equilibrium of the in-slice flow v̂ = v − θ̇ T x̂
29 tw3 = so.fsolve(lambda y: slc.velocity(np.array([y[0], 0, y[1], y[2]]), ds.G)[0][[0, 2, 3]],
30 [0.44, -0.39, 0.07], xtol=1e-14)
31 x_tw = np.array([tw3[0], 0, tw3[1], tw3[2]])
32 print("TW in the slice:", x_tw)
33 T = np.array([ops.ddx(e) for e in np.eye(4)]).T
34 w, V = np.linalg.eig(ds.J(x_tw) - slc.velocity(x_tw, ds.G)[1] * T)
35 vu = V[:, np.argmax(w.real)]
36 integ = ODEIntegrator(ds, method="DOP853", rtol=1e-12, atol=1e-13)
37 tr = integ.integrate(x_tw + 1e-3 * np.real(vu), 250.0)
38 st = slc.reduce_trajectory(tr.t, tr.X) # unwrapped phase, border distance
39 reduced = Trajectory(tr.t, st.X_hat, "reduced DOP853 dense output",
40 lambda i, t: slc.reduce(tr._dense(i, t))[0])
41 a = np.imag(vu / np.linalg.norm(vu))
42 a[1] = 0.0
43 n = np.cross(a[[0, 2, 3]], [0, 0, 1.0]) # plane ∋ TW spanned by Im v_u and ŷ2
44 sec = PoincareSection.hyperplane([n[0], 0, n[1], n[2]], x_tw, direction=-1)
45 cr = poincare_crossings(reduced, sec, xtol=1e-13, t_min=100.0)
46 print(cr.explain())
47 # closest one-return recurrences, best first; a seed that falls back onto the relative
48 # equilibrium is refused by find_rpo_in_slice (steady_tol) and the next one is tried
49 d = [np.linalg.norm(cr.X[j + 1] - cr.X[j]) for j in range(len(cr.t) - 1)]
50 for k in np.argsort(d):
51 try:
52 rpo = find_rpo_in_slice(slc, integ, sec, cr.X[k], cr.t[k + 1] - cr.t[k], tol=1e-11,
53 fd_step=1e-6, max_iterations=10, steady_tol=1e-6)
54 except Exception as e: # noqa: BLE001 -- a refused seed (any reason) is printed; try the next return
55 print("seed refused:", str(e).splitlines()[0])
56 continue
57 print(rpo.explain())
58 return rpo
59
60
61def lorenz_homoclinic():
62 pb = ode_problem(["x", "y", "z"], {"sigma": 10.0, "r": 13.9, "b": 8 / 3},
63 lambda X, P: [P["sigma"] * (X["y"] - X["x"]), P["r"] * X["x"] - X["y"] - X["x"] * X["z"],
64 X["x"] * X["y"] - P["b"] * X["z"]], "Lorenz")
65 ds = ir_system(pb, ["r"])
66 origin = np.zeros(3)
67
68 def setup(q): # q = [r]
69 par = {"r": float(q[0])}
70 man = unstable_manifold(ds, origin, order=10, degree=2, eigvec_scale=1.0,
71 unstable_tol=1e-9, par=par)
72 return dict(x0=man([0.5]), x_plus=origin, M=ds.M, end_radius=0.1,
73 integrator=ODEIntegrator(ds, method="DOP853", rtol=1e-11, atol=1e-13, par=par),
74 Psi=unstable_left_basis(ds, origin, unstable_tol=1e-9, par=par),
75 start_defect=man.defect([0.5]))
76 con = find_connection(setup, 8.0, bracket=(13.92, 13.93), xtol=1e-11)
77 print(con.explain())
78 return con
79
80
81if __name__ == "__main__":
82 rpo = two_mode_rpo()
83 con = lorenz_homoclinic()
84 # Budanur, Borrero-Echeverry & Cvitanovic, Table I: T = 3.64151221 (8 digits printed);
85 # Sparrow (1982): homoclinic orbit at r = 13.926 (3 decimals printed)
86 assert rpo is not None and abs(rpo.T - 3.64151221) < 5e-7 and rpo.residual < 1e-10, rpo
87 assert abs(con.unknowns[0] - 13.926) < 1e-3, con.unknowns