da03127Two-oscillator chaos toy: model, experiments, figures, web appJeremy Magland 1"""Two oscillators coupled by a quartic potential, vectorized over many orbits.
3 H = (px^2 + py^2) / 2 + (w1^2 x^2 + w2^2 y^2) / 2 + beta x^2 y^2
5With w1 = w2 this is the system studied by Pullen and Edmonds (1981). Scaling x, y, px, py by sqrt(E)
6shows that for fixed w1, w2 the dynamics depends only on beta E, the analogue of the strike strength
7of the plate. The forces are
9 Fx = -w1^2 x - 2 beta x y^2, Fy = -w2^2 y - 2 beta x^2 y.
11Integration: 4th-order symplectic composition (Yoshida 1990) of the kick-drift-kick leapfrog, which
12conserves energy to high accuracy over long times. Optional linear damping p -> p exp(-gamma dt) after
13each step (dp/dt = -gamma p: the amplitude decays at rate gamma / 2 and the energy at rate gamma,
14the analogue of the plate's sigma0). The tangent
15dynamics (for Lyapunov exponents) is integrated with the same splitting, linearized.
16"""
18from __future__ import annotations
20from dataclasses import dataclass
22import numpy as np
24_C = 2 ** (1 / 3)
25YOSHIDA = (1 / (2 - _C), -_C / (2 - _C), 1 / (2 - _C))
29class System:
30 w1: float = 1.0
31 w2: float = 1.0
32 beta: float = 1.0
34 def force(self, x, y):
35 return -self.w1**2 * x - 2 * self.beta * x * y * y, -self.w2**2 * y - 2 * self.beta * x * x * y
37 def tangent_force(self, x, y, dx, dy):
38 """Linearized force applied to a tangent vector (dx, dy) at (x, y)."""
39 b = self.beta
40 fxx = -self.w1**2 - 2 * b * y * y
41 fyy = -self.w2**2 - 2 * b * x * x
42 fxy = -4 * b * x * y
43 return fxx * dx + fxy * dy, fxy * dx + fyy * dy
45 def energy(self, s):
46 x, y, px, py = s
47 return 0.5 * (px**2 + py**2) + 0.5 * (self.w1**2 * x**2 + self.w2**2 * y**2) + self.beta * x**2 * y**2
49 def mode_energies(self, s):
50 """Linear (uncoupled) energies of the two oscillators."""
51 x, y, px, py = s
52 return 0.5 * (px**2 + self.w1**2 * x**2), 0.5 * (py**2 + self.w2**2 * y**2)
55def _leapfrog(sys, s, d, h):
56 """One kick-drift-kick step of size h on state s (4, M) and tangent d (4, M) or None."""
57 x, y, px, py = s
58 fx, fy = sys.force(x, y)
59 px = px + 0.5 * h * fx
60 py = py + 0.5 * h * fy
61 if d is not None:
62 dx, dy, dpx, dpy = d
63 gx, gy = sys.tangent_force(x, y, dx, dy)
64 dpx = dpx + 0.5 * h * gx
65 dpy = dpy + 0.5 * h * gy
66 dx = dx + h * dpx
67 dy = dy + h * dpy
68 x = x + h * px
69 y = y + h * py
70 fx, fy = sys.force(x, y)
71 px = px + 0.5 * h * fx
72 py = py + 0.5 * h * fy
73 if d is not None:
74 gx, gy = sys.tangent_force(x, y, dx, dy)
75 dpx = dpx + 0.5 * h * gx
76 dpy = dpy + 0.5 * h * gy
77 d = np.array([dx, dy, dpx, dpy])
78 return np.array([x, y, px, py]), d
81def step(sys, s, d, dt, gamma=0.0):
82 for c in YOSHIDA:
83 s, d = _leapfrog(sys, s, d, c * dt)
84 if gamma:
85 f = np.exp(-gamma * dt)
86 s[2:] *= f
87 if d is not None:
88 d[2:] *= f
89 return s, d
92def run(sys, s0, t_end, dt=0.01, gamma=0.0, lyap=False, renorm_every=10, record_every=None,
93 section=False, d0=None):
94 """Integrate orbits s0 (4, M) to t_end.
96 lyap: integrate a tangent vector and return the finite-time largest Lyapunov exponent of each orbit
97 (Benettin: renormalize every renorm_every steps and average the log growth), plus the running
98 log-growth sum at each record time if record_every is given.
99 section: record crossings of y = 0 with py > 0 (Poincare section), as (x, px) per orbit.
100 record_every: record the state every that many steps.
101 """
102 s = np.array(s0, dtype=float)
103 M = s.shape[1]
104 n = int(round(t_end / dt))
105 d = None
106 logsum = np.zeros(M)
107 if lyap:
108 if d0 is None:
109 rng = np.random.default_rng(0)
110 d = rng.standard_normal((4, M))
111 else:
112 d = np.array(d0, dtype=float)
113 d /= np.linalg.norm(d, axis=0)
114 rec_t, rec_s, rec_log = [], [], []
115 sec = [[] for _ in range(M)] if section else None
116 for i in range(n):
117 y_old = s[1].copy()
118 s_old = s.copy()
119 s, d = step(sys, s, d, dt, gamma)
120 if section:
121 cross = np.nonzero((y_old < 0) & (s[1] >= 0) & (s[3] > 0))[0]
122 for m in cross:
123 a = -y_old[m] / (s[1, m] - y_old[m])
124 sec[m].append(((1 - a) * s_old[0, m] + a * s[0, m], (1 - a) * s_old[2, m] + a * s[2, m]))
125 if lyap and (i + 1) % renorm_every == 0:
126 nrm = np.linalg.norm(d, axis=0)
127 logsum += np.log(nrm)
128 d /= nrm
129 if record_every and (i + 1) % record_every == 0:
130 rec_t.append((i + 1) * dt)
131 rec_s.append(s.copy())
132 if lyap:
133 rec_log.append(logsum.copy())
134 out = dict(s=s)
135 if lyap:
136 out["lyap"] = logsum / (n * dt)
137 if record_every:
138 out["t"] = np.array(rec_t)
139 out["states"] = np.array(rec_s) # (T, 4, M)
140 if lyap:
141 out["logsum"] = np.array(rec_log) # (T, M)
142 if section:
143 out["section"] = [np.array(p) for p in sec]
144 return out
147def section_initial_conditions(sys, E, n_grid):
148 """Initial conditions on the section y = 0, py > 0 at energy E: a grid over the allowed (x, px) ellipse."""
149 xm = np.sqrt(2 * E) / sys.w1
150 pm = np.sqrt(2 * E)
151 xs = np.linspace(-xm, xm, n_grid + 2)[1:-1]
152 ps = np.linspace(-pm, pm, n_grid + 2)[1:-1]
153 X, P = np.meshgrid(xs, ps)
154 X, P = X.ravel(), P.ravel()
155 rest = E - 0.5 * P**2 - 0.5 * sys.w1**2 * X**2
156 ok = rest > 1e-9 * E
157 X, P, rest = X[ok], P[ok], rest[ok]
158 return np.array([X, np.zeros_like(X), P, np.sqrt(2 * rest)])