da03127Two-oscillator chaos toy: model, experiments, figures, web appJeremy Magland 1"""Experiments with the two quartic-coupled oscillators. Each writes results/<name>.npz; plots.py draws them.
3 python experiments.py sections # Poincare sections at several energies (w2 = 1 and 1.618)
4 python experiments.py lyapunov # finite-time largest Lyapunov exponent vs energy
5 python experiments.py damped # a damped 'strike' and an ensemble of nearby strikes
6 python experiments.py all # all three, in parallel processes
8beta = 1 throughout, so the energy E is the single parameter (beta E).
9"""
11from __future__ import annotations
13import os
14import sys
15from multiprocessing import Pool
17os.environ.setdefault("OMP_NUM_THREADS", "1")
19import numpy as np # noqa: E402
21from model import System, run, section_initial_conditions # noqa: E402
23HERE = os.path.dirname(os.path.abspath(__file__))
24RES = os.path.join(HERE, "results")
25W2S = (1.0, 1.618)
28def dt_for(E):
29 """Step small enough for the stiffening at energy E (the effective frequency grows like sqrt(E))."""
30 return min(0.01, 0.1 / np.sqrt(1 + 4 * E))
33def _sections_one(args):
34 w2, E = args
35 sys_ = System(1.0, w2, 1.0)
36 s0 = section_initial_conditions(sys_, E, 9)
37 r = run(sys_, s0, 3000.0, dt=dt_for(E), lyap=True, section=True)
38 pts = r["section"]
39 return w2, E, r["lyap"], np.concatenate([np.full(len(p), i) for i, p in enumerate(pts)]), \
40 np.concatenate([p for p in pts if len(p)] or [np.zeros((0, 2))])
43def sections(pool):
44 energies = (0.25, 0.5, 1.0, 1.5, 2.0, 4.0, 8.0)
45 out = {}
46 for w2, E, lam, owner, pts in pool.imap_unordered(_sections_one, [(w2, E) for w2 in W2S for E in energies]):
47 key = f"w2_{w2:g}_E_{E:g}"
48 out[key + "_lyap"], out[key + "_owner"], out[key + "_pts"] = lam, owner, pts
49 print(f" sections w2 {w2:g} E {E:g}: {len(lam)} orbits", flush=True)
50 np.savez(os.path.join(RES, "sections.npz"), energies=np.array(energies), w2s=np.array(W2S), **out)
53def _lyap_one(args):
54 w2, E = args
55 sys_ = System(1.0, w2, 1.0)
56 s0 = section_initial_conditions(sys_, E, 16)
57 r = run(sys_, s0, 2000.0, dt=dt_for(E), lyap=True)
58 return w2, E, r["lyap"]
61def lyapunov(pool):
62 energies = np.geomspace(0.1, 100, 25)
63 lam = {}
64 for w2, E, l in pool.imap_unordered(_lyap_one, [(w2, E) for w2 in W2S for E in energies]):
65 lam[(w2, E)] = l
66 print(f" lyapunov w2 {w2:g} E {E:.3g}: median {np.median(l):.3f}", flush=True)
67 np.savez(os.path.join(RES, "lyapunov.npz"), energies=energies, w2s=np.array(W2S),
68 **{f"w2_{w2:g}_E{i}": lam[(w2, E)] for w2 in W2S for i, E in enumerate(energies)})
71def _damped_one(args):
72 """A 'strike': both oscillators start at rest at the origin with momentum sqrt(2 E0) in direction a,
73 then decay with p -> p exp(-gamma dt). An ensemble perturbs a by i * 1e-10. Records the mode energies
74 and the running Lyapunov log-sum of member 0."""
75 E0, gamma, n_members = args
76 sys_ = System(1.0, 1.618, 1.0)
77 a = 0.7 + 1e-10 * np.arange(n_members)
78 p = np.sqrt(2 * E0)
79 s0 = np.array([np.zeros(n_members), np.zeros(n_members), p * np.cos(a), p * np.sin(a)])
80 r = run(sys_, s0, 1500.0, dt=dt_for(E0), gamma=gamma, lyap=True, renorm_every=10,
81 record_every=int(round(0.5 / dt_for(E0))))
82 Ex, Ey = sys_.mode_energies(np.moveaxis(r["states"], 1, 0))
83 return E0, r["t"], Ex, Ey, sys_.energy(np.moveaxis(r["states"], 1, 0)), r["logsum"]
86def damped(pool):
87 gamma = 0.005
88 out = {}
89 for E0, t, Ex, Ey, Etot, logsum in pool.imap_unordered(_damped_one, [(E0, gamma, 16) for E0 in (0.5, 3.0, 20.0)]):
90 k = f"E0_{E0:g}"
91 out[k + "_t"], out[k + "_Ex"], out[k + "_Ey"], out[k + "_E"], out[k + "_logsum"] = t, Ex, Ey, Etot, logsum
92 print(f" damped E0 {E0:g} done", flush=True)
93 np.savez(os.path.join(RES, "damped.npz"), gamma=gamma, **out)
96if __name__ == "__main__":
97 os.makedirs(RES, exist_ok=True)
98 which = sys.argv[1:] or ["all"]
99 if "all" in which:
100 which = ["sections", "lyapunov", "damped"]
101 with Pool(int(os.environ.get("WORKERS", "16"))) as pool:
102 for w in which:
103 print(w, flush=True)
104 dict(sections=sections, lyapunov=lyapunov, damped=damped)[w](pool)