""" three-body-chaos.py Science Journaling Club, Volume 3 Issue 1, Fall 2026. QUESTION -------- The three-body problem is the standard example handed out whenever someone wants a one-line illustration of chaos: three point masses, one force, no closed-form solution, sensitive dependence on initial conditions. That claim is usually just asserted. This program measures it. How fast do two nearby starting conditions actually pull apart, expressed as a number (a Lyapunov exponent) rather than a gesture? How far into the future can such a system be predicted, in concrete units of time, before a given initial measurement error has grown to the size of the system itself? And when three bodies are dropped at random and left to interact, what fraction of the time does one of them get thrown out permanently, and how does that depend on how tightly bound the system started out? WHAT THIS PROGRAM IS --------------------- The club has no telescope and no laboratory. This is not an observation of real orbiting bodies; nothing here was tracked through a telescope or measured from real astronomical data. It is a numerical experiment: a Newtonian gravity simulator that this file writes, runs, and validates against known solutions, entirely on a laptop. Where the text below says "measured", it means "computed by this code and printed below", the same way you would measure the frequency of heads by flipping a coin many times rather than reasoning about coins. MODEL ----- Three point masses under mutual Newtonian gravity, G = 1, restricted to a plane (a 2D three-body problem; the genuinely 3D problem is qualitatively similar but plane motion keeps every figure drawable and does not remove the chaos, which is already present with 2 degrees of freedom per body). Equations of motion, body i: x_i'' = sum_{j != i} G m_j (x_j - x_i) / r_ij^3. TWO INTEGRATION REGIMES, CHOSEN FOR HONESTY, NOT CONVENIENCE -------------------------------------------------------------- This file uses a fixed-step 4th-order symplectic integrator (Yoshida 1990: three velocity-Verlet leapfrog sub-steps composed with the "triple jump" coefficients w1 = 1/(2-2^(1/3)), w0 = 1-2w1) for SECTION 1 and SECTION 2 (the two-body Kepler orbit and the figure-eight periodic solution), where it is superb: energy conserved to 1e-13, because neither of those two orbits ever passes through a close encounter. It is NOT adequate, at any fixed step size we could afford, for a system that passes through a close triple encounter -- and generic three-body motion does. We found this out directly, the way this whole project insists things get found out: by checking. A first version of the Monte Carlo sections below used a fixed step with a small softening length, and its own conserved-quantity monitor caught energy drifting by factors of 10 to 1000 during close encounters -- a fixed step simply cannot resolve how fast the force changes when two bodies swing close together. The Pythagorean three-body problem (SECTION 3b) shows this worst: at a fixed step of 0.001 time units it reaches a relative energy error of hundreds; even at 0.0001 it is still large. This is a genuine, well known numerical-methods fact about this configuration, not a bug we shipped and caught quietly, and both fixed step diagnostics are printed below exactly as they came out. The fix used everywhere a close encounter is possible (SECTION 3 onward) is an ADAPTIVE step size tied to the local two-body free-fall time at the closest pair: at every step, dt = eta * 2*pi*sqrt(r_min^3 / (G*(m_i+m_j))) (capped at dt_max) where r_min is the current smallest pairwise separation and m_i, m_j are that pair's masses -- essentially a small fraction (eta = 0.005 by default) of the local two-body orbital period for the tightest pair at that instant. The same three-substep Yoshida composition still advances each step, but because the step size now depends on the state, the method is no longer exactly symplectic in the strict sense (that requires a time transformation this file does not implement, e.g. Mikkola's algorithmic regularization). We do not claim symplecticity for these sections. We claim, and print, the thing that actually matters: measured conservation of energy, linear momentum and angular momentum, checked on every run. On the same Pythagorean problem, this scheme reaches a relative energy error under 1e-3, using far fewer steps than the failed fixed-step attempt. ASSUMPTIONS ----------- * Point masses, no collisions, no relativity, no tides, no external field. * Motion confined to a plane (a plane already gives a fully chaotic, 4-degree-of-freedom-per-body-pair system; a real 3D triple would differ quantitatively but not, we expect, in kind). * A tiny safety length (1e-6) is added in quadrature to every squared separation purely to guard against division by exactly zero; SS9 checks directly that it has no measurable effect on the Monte Carlo answer. * Escape/ejection (SS5) is judged by a simplified two-body criterion: a body is called ejected once it sits beyond a fixed radius from the center of mass of the other two, still moving away, with positive specific orbital energy relative to that pair treated as one point mass. Standard in three-body scattering experiments, an approximation to the true three-body energy partition. * "Ejected within the run" is a finite-time measurement (every Monte Carlo trial runs at most T_MAX time units). Classical results say that almost every bound, non-periodic triple system eventually disrupts given unlimited time; a system that survives to T_MAX here is reported as "not yet ejected by T_MAX", never as "stable forever". * The Lyapunov exponent of a randomly drawn chaotic triple is not one number with a tight scatter: SECTION 3c finds a heavy right tail, driven by the rare trial that happens to pass unusually close during a triple encounter. We report the mean AND the median, and use the median as the more representative "typical" value, precisely because the mean of a heavy-tailed quantity is pulled around by its largest few draws. LIMITATIONS ----------- * Two dimensions, not three. * The adaptive step size responds to the single closest pair; a simultaneous close approach of a different pair elsewhere in a wide, hierarchical configuration is not separately protected (not an issue for the compact, comparable-mass systems studied here, but worth naming). * The escape criterion is a simplification; a handful of borderline cases at the margin could be mis-classified. * A Lyapunov exponent measured over a finite window is itself an estimate with sampling noise; every value below carries the trial count it used. * eta = 0.005 is one choice; SS9 measures what a coarser and a finer eta do to the headline Monte Carlo answer. SEED ---- Master seed 20260914 (today's date), hard-coded below. Every random stream used anywhere in this file (initial triangles, initial velocity directions, perturbation directions) is spawned from it through numpy's SeedSequence, so the entire output is deterministic and reproducible from this file alone. RUN --- python three-body-chaos.py > three-body-chaos-output.txt Needs Python 3.12 and numpy (numpy is used only for random-number generation and for assembling printed summary statistics; every physics step -- forces, the integrator, the adaptive step-size rule, the Lyapunov renormalisation, the escape check -- is plain Python, hand-unrolled for 3 bodies in 2D, because for an object this small numpy's own per-call array overhead is slower than direct arithmetic on Python floats). Expect on the order of 5 to 9 minutes on a laptop. """ import math import time import numpy as np SEED = 20260914 print(f"MASTER SEED: {SEED}") _ss = np.random.SeedSequence(SEED) _children = _ss.spawn(8) RNG_TRIANGLES = np.random.default_rng(_children[0]) RNG_PERTURB = np.random.default_rng(_children[2]) RNG_MC = np.random.default_rng(_children[3]) T_START = time.perf_counter() def section(title): print() print("=" * 78) print(title) print("=" * 78) # --------------------------------------------------------------------------- # CORE PHYSICS: forces (hand-unrolled for exactly 3 bodies in 2D) and the # 4th-order symplectic (Yoshida) step. Plain Python throughout -- for a # 3-body, 2D system the arrays are so small that numpy's per-call overhead # is slower than direct arithmetic. Every number in this study passes # through accelerations() and yoshida4_step() below. # --------------------------------------------------------------------------- SAFETY_EPS2 = (1.0e-6) ** 2 # guards only against exact division by zero def accelerations(pos, masses, eps2=0.0): """pos: [[x1,y1],[x2,y2],[x3,y3]]; masses: [m1,m2,m3]. G = 1.""" x1, y1 = pos[0] x2, y2 = pos[1] x3, y3 = pos[2] m1, m2, m3 = masses dx12 = x2 - x1; dy12 = y2 - y1 r2_12 = dx12 * dx12 + dy12 * dy12 + eps2 inv12 = r2_12 ** -1.5 dx13 = x3 - x1; dy13 = y3 - y1 r2_13 = dx13 * dx13 + dy13 * dy13 + eps2 inv13 = r2_13 ** -1.5 dx23 = x3 - x2; dy23 = y3 - y2 r2_23 = dx23 * dx23 + dy23 * dy23 + eps2 inv23 = r2_23 ** -1.5 ax1 = m2 * dx12 * inv12 + m3 * dx13 * inv13 ay1 = m2 * dy12 * inv12 + m3 * dy13 * inv13 ax2 = -m1 * dx12 * inv12 + m3 * dx23 * inv23 ay2 = -m1 * dy12 * inv12 + m3 * dy23 * inv23 ax3 = -m1 * dx13 * inv13 - m2 * dx23 * inv23 ay3 = -m1 * dy13 * inv13 - m2 * dy23 * inv23 return [[ax1, ay1], [ax2, ay2], [ax3, ay3]] def vv_substep(pos, vel, acc, masses, dt, eps2): """One velocity-Verlet (2nd order, symplectic) sub-step. acc must be the acceleration AT pos on entry; returns new pos, vel, acc.""" p0, p1, p2 = pos v0, v1, v2 = vel a0, a1, a2 = acc half_dt2 = 0.5 * dt * dt newpos = [ [p0[0] + v0[0] * dt + a0[0] * half_dt2, p0[1] + v0[1] * dt + a0[1] * half_dt2], [p1[0] + v1[0] * dt + a1[0] * half_dt2, p1[1] + v1[1] * dt + a1[1] * half_dt2], [p2[0] + v2[0] * dt + a2[0] * half_dt2, p2[1] + v2[1] * dt + a2[1] * half_dt2], ] newacc = accelerations(newpos, masses, eps2) na0, na1, na2 = newacc half_dt = 0.5 * dt newvel = [ [v0[0] + (a0[0] + na0[0]) * half_dt, v0[1] + (a0[1] + na0[1]) * half_dt], [v1[0] + (a1[0] + na1[0]) * half_dt, v1[1] + (a1[1] + na1[1]) * half_dt], [v2[0] + (a2[0] + na2[0]) * half_dt, v2[1] + (a2[1] + na2[1]) * half_dt], ] return newpos, newvel, newacc _W1 = 1.0 / (2.0 - 2.0 ** (1.0 / 3.0)) _W0 = 1.0 - 2.0 * _W1 def yoshida4_step(pos, vel, acc, masses, dt, eps2=0.0): """One 4th-order composition step, three velocity-Verlet sub-steps.""" pos, vel, acc = vv_substep(pos, vel, acc, masses, _W1 * dt, eps2) pos, vel, acc = vv_substep(pos, vel, acc, masses, _W0 * dt, eps2) pos, vel, acc = vv_substep(pos, vel, acc, masses, _W1 * dt, eps2) return pos, vel, acc def integrate_fixed(pos0, vel0, masses, dt, nsteps, eps2=0.0, record_every=0): """Fixed-step Yoshida-4 integration. Used only where no close encounter occurs (SECTION 1, SECTION 2).""" pos = [list(p) for p in pos0] vel = [list(v) for v in vel0] acc = accelerations(pos, masses, eps2) history = [] if record_every > 0: history.append((0.0, [list(p) for p in pos], [list(v) for v in vel])) for step in range(1, nsteps + 1): pos, vel, acc = yoshida4_step(pos, vel, acc, masses, dt, eps2) if record_every > 0 and (step % record_every == 0 or step == nsteps): history.append((step * dt, [list(p) for p in pos], [list(v) for v in vel])) if record_every > 0: return pos, vel, history return pos, vel def min_pair_dist(pos, masses): """Smallest pairwise separation and the summed mass of that pair.""" best_r2 = 1e18 best_M = 1.0 for i in range(3): for j in range(i + 1, 3): dx = pos[i][0] - pos[j][0] dy = pos[i][1] - pos[j][1] r2 = dx * dx + dy * dy if r2 < best_r2: best_r2 = r2 best_M = masses[i] + masses[j] return math.sqrt(best_r2), best_M def adaptive_dt(pos, masses, eta, dt_max): """A fraction eta of the local two-body free-fall/orbital time at the closest pair, capped at dt_max. Returns (dt, r_min).""" r_min, M_pair = min_pair_dist(pos, masses) r_eff = math.sqrt(r_min * r_min + SAFETY_EPS2) t_dyn = 2 * math.pi * math.sqrt(r_eff ** 3 / M_pair) dt = eta * t_dyn if dt > dt_max: dt = dt_max return dt, r_min def energy(pos, vel, masses, eps2=0.0): ke = 0.0 for i in range(3): ke += 0.5 * masses[i] * (vel[i][0] ** 2 + vel[i][1] ** 2) pe = 0.0 for i in range(3): for j in range(i + 1, 3): dx = pos[j][0] - pos[i][0] dy = pos[j][1] - pos[i][1] r = math.sqrt(dx * dx + dy * dy + eps2) pe -= masses[i] * masses[j] / r return ke + pe def linear_momentum(vel, masses): px = sum(masses[i] * vel[i][0] for i in range(3)) py = sum(masses[i] * vel[i][1] for i in range(3)) return px, py def angular_momentum(pos, vel, masses): L = 0.0 for i in range(3): L += masses[i] * (pos[i][0] * vel[i][1] - pos[i][1] * vel[i][0]) return L def state_vector(pos, vel): v = [] for i in range(len(pos)): v.append(pos[i][0]); v.append(pos[i][1]) for i in range(len(vel)): v.append(vel[i][0]); v.append(vel[i][1]) return v def state_diff_norm(pos_a, vel_a, pos_b, vel_b): sa = state_vector(pos_a, vel_a) sb = state_vector(pos_b, vel_b) return math.sqrt(sum((a - b) ** 2 for a, b in zip(sa, sb))) def all_finite(pos, vel): for b in pos: if not (math.isfinite(b[0]) and math.isfinite(b[1])): return False for b in vel: if not (math.isfinite(b[0]) and math.isfinite(b[1])): return False return True # --------------------------------------------------------------------------- # SECTION 1: VALIDATE ON THE TWO-BODY PROBLEM (fixed step, bare 1/r^2) # Same accelerations()/yoshida4_step() functions, third body given zero # mass and parked far away so it neither feels nor exerts gravity over the # run: a genuine test of the exact code the rest of the study builds on. # --------------------------------------------------------------------------- section("VALIDATION 1: TWO-BODY PROBLEM AGAINST THE ANALYTIC KEPLER ORBIT") M1, M2 = 1.0, 1.0 GM = M1 + M2 ECC = 0.5 A_REL = 1.0 r_peri = A_REL * (1 - ECC) v_peri = math.sqrt(GM * (2.0 / r_peri - 1.0 / A_REL)) T_KEPLER = 2 * math.pi * math.sqrt(A_REL ** 3 / GM) rel_pos0 = (r_peri, 0.0) rel_vel0 = (0.0, v_peri) pos0_2b = [ [-(M2 / GM) * rel_pos0[0], -(M2 / GM) * rel_pos0[1]], [(M1 / GM) * rel_pos0[0], (M1 / GM) * rel_pos0[1]], [1.0e6, 1.0e6], ] vel0_2b = [ [-(M2 / GM) * rel_vel0[0], -(M2 / GM) * rel_vel0[1]], [(M1 / GM) * rel_vel0[0], (M1 / GM) * rel_vel0[1]], [0.0, 0.0], ] masses_2b = [M1, M2, 0.0] DT_2B = T_KEPLER / 20000.0 N_SAMPLES_2B = 25 def kepler_solve(M, e, tol=1e-13, maxiter=100): E = M if e < 0.8 else math.pi for _ in range(maxiter): f = E - e * math.sin(E) - M fp = 1 - e * math.cos(E) dE = -f / fp E += dE if abs(dE) < tol: break return E def kepler_relative_state(t, a, e, GM_): n = math.sqrt(GM_ / a ** 3) M = math.fmod(n * t, 2 * math.pi) E = kepler_solve(M, e) cosE, sinE = math.cos(E), math.sin(E) x = a * (cosE - e) y = a * math.sqrt(1 - e * e) * sinE Edot = n / (1 - e * cosE) vx = -a * sinE * Edot vy = a * math.sqrt(1 - e * e) * cosE * Edot return (x, y), (vx, vy) pos_2b, vel_2b, hist_2b = integrate_fixed(pos0_2b, vel0_2b, masses_2b, DT_2B, 20000, eps2=0.0, record_every=20000 // N_SAMPLES_2B) max_rel_pos_err = 0.0 max_rel_vel_err = 0.0 for (t, p, v) in hist_2b: rel_p_num = (p[1][0] - p[0][0], p[1][1] - p[0][1]) rel_v_num = (v[1][0] - v[0][0], v[1][1] - v[0][1]) rel_p_an, rel_v_an = kepler_relative_state(t, A_REL, ECC, GM) dp = math.hypot(rel_p_num[0] - rel_p_an[0], rel_p_num[1] - rel_p_an[1]) dv = math.hypot(rel_v_num[0] - rel_v_an[0], rel_v_num[1] - rel_v_an[1]) scale_p = math.hypot(*rel_p_an) + 1e-300 scale_v = math.hypot(*rel_v_an) + 1e-300 max_rel_pos_err = max(max_rel_pos_err, dp / scale_p) max_rel_vel_err = max(max_rel_vel_err, dv / scale_v) print(f"Analytic Kepler period: T = {T_KEPLER:.10f}") print(f"Eccentricity: e = {ECC}") print(f"Semi-major axis (relative): a = {A_REL}") print(f"Integrator step: dt = {DT_2B:.6e} ({T_KEPLER/DT_2B:.0f} steps/period, fixed)") print("Sampled 25 points across one full orbit and compared the club integrator's") print("relative position and velocity against the closed-form Kepler-equation") print("solution (Newton's method on E - e sin E = M) at each one.") print(f"Max relative position error over the orbit: {max_rel_pos_err:.3e}") print(f"Max relative velocity error over the orbit: {max_rel_vel_err:.3e}") TWO_BODY_TOL = 1e-8 print(f"Tolerance: {TWO_BODY_TOL:.1e}") TWO_BODY_PASS = max(max_rel_pos_err, max_rel_vel_err) < TWO_BODY_TOL print("PASS" if TWO_BODY_PASS else "FAIL", "-- two-body integrator matches the analytic solution.") pos_ret, vel_ret = integrate_fixed(pos0_2b, vel0_2b, masses_2b, DT_2B, 20000, eps2=0.0) return_err = state_diff_norm(pos0_2b[:2], vel0_2b[:2], pos_ret[:2], vel_ret[:2]) scale = math.hypot(*rel_pos0) print(f"Return-to-start residual after exactly one period: {return_err:.3e}" f" (relative to orbit size {scale:.3e}: {return_err/scale:.3e})") E0_2b = energy(pos0_2b, vel0_2b, masses_2b, eps2=0.0) E1_2b = energy(pos_ret, vel_ret, masses_2b, eps2=0.0) print(f"Energy at t=0: {E0_2b:.12f}") print(f"Energy after 1 period: {E1_2b:.12f}") print(f"Relative energy drift over one period: {abs((E1_2b-E0_2b)/E0_2b):.3e}") # --------------------------------------------------------------------------- # SECTION 2: VALIDATE ON THE FIGURE-EIGHT PERIODIC SOLUTION (fixed step) # --------------------------------------------------------------------------- section("VALIDATION 2: THE FIGURE-EIGHT THREE-BODY SOLUTION") pos0_fig8 = [ [0.9700436, -0.24308753], [-0.9700436, 0.24308753], [0.0, 0.0], ] vel0_fig8 = [ [0.466203685, 0.43236573], [0.466203685, 0.43236573], [-0.93240737, -0.86473146], ] masses_fig8 = [1.0, 1.0, 1.0] T_FIG8_LIT = 6.3259138848 # commonly quoted literature period (Simo 2001) DT_FIG8_SEARCH = T_FIG8_LIT / 60000.0 # The published initial condition is quoted to only 8-9 significant figures, # so it does not close the loop exactly at the literature period -- it closes # as well as those figures allow. Rather than assume the literature period is # exact, we measure our own best return time: integrate a little past the # nominal period, record the distance back to the start at every step near # it, and locate the true minimum with a parabolic fit through the three # closest samples. pos_s = [list(p) for p in pos0_fig8]; vel_s = [list(v) for v in vel0_fig8] acc_s = accelerations(pos_s, masses_fig8, 0.0) search_records = [] for step in range(1, 62001): pos_s, vel_s, acc_s = yoshida4_step(pos_s, vel_s, acc_s, masses_fig8, DT_FIG8_SEARCH, 0.0) if step > 58000: d = state_diff_norm(pos_s, vel_s, pos0_fig8, vel0_fig8) search_records.append((step * DT_FIG8_SEARCH, d)) best_idx = min(range(len(search_records)), key=lambda i: search_records[i][1]) t1, d1 = search_records[best_idx - 1] t2, d2 = search_records[best_idx] t3, d3 = search_records[best_idx + 1] y1, y2, y3 = d1 * d1, d2 * d2, d3 * d3 denom = (t1 - t2) * (t1 - t3) * (t2 - t3) A_fit = (t3 * (y2 - y1) + t2 * (y1 - y3) + t1 * (y3 - y2)) / denom B_fit = (t3 * t3 * (y1 - y2) + t2 * t2 * (y3 - y1) + t1 * t1 * (y2 - y3)) / denom T_FIG8 = -B_fit / (2 * A_fit) best_residual = min(d1, d2, d3) # Step 60000 of this same search integration lands at exactly t = T_FIG8_LIT # (DT_FIG8_SEARCH was built as T_FIG8_LIT / 60000), so the article's claim # about integrating for "exactly that literature period" is this one entry # of search_records, read directly rather than recomputed -- printed here so # that specific number is traceable to this output file instead of requiring # a rerun to reproduce. _lit_idx = 60000 - 58001 _t_at_lit, _d_at_lit = search_records[_lit_idx] assert abs(_t_at_lit - T_FIG8_LIT) < 1e-9, "search_records index for T_FIG8_LIT drifted" _scale_f8_lit = max(math.hypot(*p) for p in pos0_fig8) print(f"Literature period (Simo 2001), commonly quoted: T_lit = {T_FIG8_LIT}") print(f"Residual |state(T_lit) - state(0)| integrating for exactly that literature " f"period: {_d_at_lit:.3e} (relative to orbit scale {_scale_f8_lit:.4f}: " f"{_d_at_lit/_scale_f8_lit:.3e}), too large to call a clean return.") print(f"Self-measured return time (parabolic fit around the closest numerical") print(f"return to the initial state, using this exact code): T = {T_FIG8:.10f}") print(f"Difference from the literature value: {T_FIG8-T_FIG8_LIT:+.6f} time units " f"({(T_FIG8-T_FIG8_LIT)/T_FIG8_LIT*100:+.4f}%), consistent with the widely " f"circulated initial condition being quoted to only 8-9 significant figures.") DT_FIG8 = T_FIG8 / 60000.0 NSTEPS_FIG8 = 60000 pos_f8, vel_f8, hist_f8 = integrate_fixed(pos0_fig8, vel0_fig8, masses_fig8, DT_FIG8, NSTEPS_FIG8, eps2=0.0, record_every=200) return_err_f8 = state_diff_norm(pos0_fig8, vel0_fig8, pos_f8, vel_f8) scale_f8 = max(math.hypot(*p) for p in pos0_fig8) print(f"\nRe-integrating cleanly at the self-measured period, {NSTEPS_FIG8} steps:") print(f"Residual |state(T) - state(0)|: {return_err_f8:.3e}") print(f"Relative to the orbit's own length scale ({scale_f8:.4f}): {return_err_f8/scale_f8:.3e}") FIG8_TOL = 5e-4 print(f"Tolerance: {FIG8_TOL:.1e} (set by the ~8-9 significant figures published for this") print("initial condition, not by the integrator -- its own energy conservation, below,") print("holds to 1e-13).") FIG8_PASS = return_err_f8 / scale_f8 < FIG8_TOL print("PASS" if FIG8_PASS else "FAIL", "-- figure-eight returns to its initial state within tolerance after one period.") E0_f8 = energy(pos0_fig8, vel0_fig8, masses_fig8, eps2=0.0) Px0_f8, Py0_f8 = linear_momentum(vel0_fig8, masses_fig8) L0_f8 = angular_momentum(pos0_fig8, vel0_fig8, masses_fig8) max_dE = max_dP = max_dL = 0.0 for (t, p, v) in hist_f8: E = energy(p, v, masses_fig8, eps2=0.0) Px, Py = linear_momentum(v, masses_fig8) L = angular_momentum(p, v, masses_fig8) max_dE = max(max_dE, abs((E - E0_f8) / E0_f8)) max_dP = max(max_dP, math.hypot(Px - Px0_f8, Py - Py0_f8)) max_dL = max(max_dL, abs(L - L0_f8)) print(f"Initial energy E0 = {E0_f8:.10f}, initial angular momentum L0 = {L0_f8:.10f}") print(f"Over one period: max |relative energy drift| = {max_dE:.3e}") print(f"Over one period: max |linear momentum drift| = {max_dP:.3e} (started at {math.hypot(Px0_f8,Py0_f8):.2e})") L0_scale = abs(L0_f8) if abs(L0_f8) > 1e-9 else 1.0 print(f"Over one period: max |angular momentum drift| = {max_dL:.3e} " f"(L0 = {L0_f8:.3e}: this choreography has zero net angular momentum by " f"construction, so drift is judged against a scale of 1 when L0 itself is " f"numerically zero: {max_dL/L0_scale:.3e})") pos_f8_20, vel_f8_20, hist_f8_20 = integrate_fixed(pos0_fig8, vel0_fig8, masses_fig8, DT_FIG8, NSTEPS_FIG8 * 20, eps2=0.0, record_every=NSTEPS_FIG8) max_dE_20 = 0.0 for (t, p, v) in hist_f8_20: E = energy(p, v, masses_fig8, eps2=0.0) max_dE_20 = max(max_dE_20, abs((E - E0_f8) / E0_f8)) print(f"Over 20 periods ({20*T_FIG8:.2f} time units): max |relative energy drift| " f"= {max_dE_20:.3e} (bounded, not growing secularly).") # --------------------------------------------------------------------------- # SECTION 3: LYAPUNOV EXPONENT -- Benettin's two-trajectory method with # periodic renormalisation, on the ADAPTIVE-STEP integrator (see docstring). # --------------------------------------------------------------------------- section("LYAPUNOV EXPONENT MEASUREMENT (Benettin two-trajectory method)") DELTA0 = 1.0e-8 RENORM_EVERY = 0.5 ETA_DEFAULT = 0.005 DT_MAX_DEFAULT = 0.02 MAX_SUBSTEPS_PER_RENORM = 200000 def lyapunov_run(pos0, vel0, masses, t_total, rng, delta0=DELTA0, renorm_every=RENORM_EVERY, eta=ETA_DEFAULT, dt_max=DT_MAX_DEFAULT): pos_ref = [list(p) for p in pos0] vel_ref = [list(v) for v in vel0] acc_ref = accelerations(pos_ref, masses, SAFETY_EPS2) direction = rng.normal(size=12) direction = direction / np.linalg.norm(direction) pos_per = [list(p) for p in pos0] vel_per = [list(v) for v in vel0] for i in range(3): pos_per[i][0] += delta0 * direction[4 * i] pos_per[i][1] += delta0 * direction[4 * i + 1] vel_per[i][0] += delta0 * direction[4 * i + 2] vel_per[i][1] += delta0 * direction[4 * i + 3] acc_per = accelerations(pos_per, masses, SAFETY_EPS2) n_renorms = max(1, round(t_total / renorm_every)) log_sum = 0.0 running = [] t = 0.0 min_dist_overall = float("inf") for r in range(n_renorms): t_target = (r + 1) * renorm_every substeps = 0 while t < t_target - 1e-13: dt, rmin = adaptive_dt(pos_ref, masses, eta, dt_max) if rmin < min_dist_overall: min_dist_overall = rmin if t + dt > t_target: dt = t_target - t pos_ref, vel_ref, acc_ref = yoshida4_step(pos_ref, vel_ref, acc_ref, masses, dt, SAFETY_EPS2) pos_per, vel_per, acc_per = yoshida4_step(pos_per, vel_per, acc_per, masses, dt, SAFETY_EPS2) t += dt substeps += 1 if substeps > MAX_SUBSTEPS_PER_RENORM or not all_finite(pos_ref, vel_ref): return float("nan"), running, min_dist_overall D = state_diff_norm(pos_ref, vel_ref, pos_per, vel_per) if D <= 0 or not math.isfinite(D): break log_sum += math.log(D / delta0) running.append((t, log_sum / t)) scale = delta0 / D sref = state_vector(pos_ref, vel_ref) sper = state_vector(pos_per, vel_per) newper = [sref[k] + (sper[k] - sref[k]) * scale for k in range(12)] pos_per = [[newper[0], newper[1]], [newper[2], newper[3]], [newper[4], newper[5]]] vel_per = [[newper[6], newper[7]], [newper[8], newper[9]], [newper[10], newper[11]]] acc_per = accelerations(pos_per, masses, SAFETY_EPS2) lam = log_sum / t if t > 0 else float("nan") return lam, running, min_dist_overall # --- 3a: the figure-eight orbit itself ------------------------------------- print("\n--- 3a. Figure-eight orbit: is the exact periodic solution chaotic? ---") t0 = time.perf_counter() N_FIG8_LYAP_TRIALS = 10 fig8_lams = [] for k in range(N_FIG8_LYAP_TRIALS): lam, running, mind = lyapunov_run(pos0_fig8, vel0_fig8, masses_fig8, t_total=25 * T_FIG8, rng=RNG_PERTURB) fig8_lams.append(lam) fig8_lams = np.array(fig8_lams) print(f"{N_FIG8_LYAP_TRIALS} independent random perturbation directions, " f"each integrated for 25 periods ({25*T_FIG8:.1f} time units).") print(f"Mean Lyapunov exponent: {fig8_lams.mean():.5f} " f"(SD {fig8_lams.std(ddof=1):.5f}, SEM {fig8_lams.std(ddof=1)/math.sqrt(len(fig8_lams)):.5f})") print(f"Lyapunov TIME (1/lambda): {1.0/fig8_lams.mean():.2f} time units, " f"or {1.0/fig8_lams.mean()/T_FIG8:.2f} figure-eight periods.") print(f"(elapsed {time.perf_counter()-t0:.1f} s)") # --- 3b: the Pythagorean three-body problem, a famous standard config ------ print("\n--- 3b. The Pythagorean three-body problem (Burrau 1913 initial data) ---") masses_pyth = [3.0, 4.0, 5.0] pos0_pyth = [[1.0, 3.0], [-2.0, -1.0], [1.0, -1.0]] vel0_pyth = [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]] d12 = math.hypot(pos0_pyth[0][0]-pos0_pyth[1][0], pos0_pyth[0][1]-pos0_pyth[1][1]) d13 = math.hypot(pos0_pyth[0][0]-pos0_pyth[2][0], pos0_pyth[0][1]-pos0_pyth[2][1]) d23 = math.hypot(pos0_pyth[1][0]-pos0_pyth[2][0], pos0_pyth[1][1]-pos0_pyth[2][1]) print(f"Triangle side lengths (should be 3, 4, 5): {d23:.6f}, {d13:.6f}, {d12:.6f}") E0_pyth = energy(pos0_pyth, vel0_pyth, masses_pyth, eps2=0.0) print(f"Initial energy (all velocities zero): E0 = {E0_pyth:.6f}") # Demonstrate, honestly, why a fixed step fails here before using the # adaptive integrator for the real measurement. print("\nFixed-step diagnostic (this is the failure described in the docstring):") for dt_fixed in [0.001, 0.0001]: pos_d = [list(p) for p in pos0_pyth]; vel_d = [list(v) for v in vel0_pyth] fixed_eps2 = (0.001) ** 2 acc_d = accelerations(pos_d, masses_pyth, fixed_eps2) nsteps_d = int(30.0 / dt_fixed) for _ in range(nsteps_d): pos_d, vel_d, acc_d = yoshida4_step(pos_d, vel_d, acc_d, masses_pyth, dt_fixed, fixed_eps2) Ef_d = energy(pos_d, vel_d, masses_pyth, eps2=fixed_eps2) E0_d = energy(pos0_pyth, vel0_pyth, masses_pyth, eps2=fixed_eps2) print(f" fixed dt={dt_fixed:.5f}, softening 0.001, to t=30: " f"relative energy drift = {abs((Ef_d-E0_d)/E0_d):.3e}") print("\nAdaptive-step result (eta=0.005, used for every quantitative number below):") t0 = time.perf_counter() N_PYTH_LYAP_TRIALS = 10 T_PYTH = 65.0 pyth_lams = [] pyth_min_dists = [] for k in range(N_PYTH_LYAP_TRIALS): lam, running, mind = lyapunov_run(pos0_pyth, vel0_pyth, masses_pyth, t_total=T_PYTH, rng=RNG_PERTURB) pyth_lams.append(lam) pyth_min_dists.append(mind) pyth_lams = np.array(pyth_lams) print(f"{N_PYTH_LYAP_TRIALS} independent perturbation directions, each run to t={T_PYTH}.") n_finite_pyth = int(np.sum(np.isfinite(pyth_lams))) pyth_lams_f = pyth_lams[np.isfinite(pyth_lams)] print(f"{n_finite_pyth} of {N_PYTH_LYAP_TRIALS} completed without hitting the numerical safety cap.") print(f"Mean Lyapunov exponent: {pyth_lams_f.mean():.5f} " f"(SD {pyth_lams_f.std(ddof=1):.5f}, SEM {pyth_lams_f.std(ddof=1)/math.sqrt(len(pyth_lams_f)):.5f})") print(f"Lyapunov time (1/lambda): {1.0/pyth_lams_f.mean():.3f} time units.") print(f"Closest triple-encounter distance reached, across trials: {min(pyth_min_dists):.6f}") # Reference trajectory conservation check at the production eta. pos_ref_p = [list(p) for p in pos0_pyth]; vel_ref_p = [list(v) for v in vel0_pyth] acc_ref_p = accelerations(pos_ref_p, masses_pyth, SAFETY_EPS2) t_p = 0.0 min_d_pyth_ref = float("inf") steps_p = 0 while t_p < T_PYTH: dt_p, rmin_p = adaptive_dt(pos_ref_p, masses_pyth, ETA_DEFAULT, DT_MAX_DEFAULT) if rmin_p < min_d_pyth_ref: min_d_pyth_ref = rmin_p if t_p + dt_p > T_PYTH: dt_p = T_PYTH - t_p pos_ref_p, vel_ref_p, acc_ref_p = yoshida4_step(pos_ref_p, vel_ref_p, acc_ref_p, masses_pyth, dt_p, SAFETY_EPS2) t_p += dt_p steps_p += 1 E_end_pyth = energy(pos_ref_p, vel_ref_p, masses_pyth, eps2=0.0) print(f"Reference (unperturbed) trajectory to t={T_PYTH}, adaptive eta={ETA_DEFAULT}: " f"{steps_p} steps, closest approach {min_d_pyth_ref:.6f}, " f"relative energy drift = {abs((E_end_pyth-E0_pyth)/E0_pyth):.3e}") print(f"(elapsed {time.perf_counter()-t0:.1f} s)") # --- 3c: random generic configurations, for the headline exponent ---------- print("\n--- 3c. Random generic bound triples: the headline Lyapunov exponent ---") t0 = time.perf_counter() def sample_triangle(rng, scale=1.0, min_sep=0.3): pts = None for _ in range(1000): pts = rng.uniform(-scale, scale, size=(3, 2)) d01 = math.hypot(*(pts[0]-pts[1])) d02 = math.hypot(*(pts[0]-pts[2])) d12v = math.hypot(*(pts[1]-pts[2])) if min(d01, d02, d12v) > min_sep: break return [list(pts[0]), list(pts[1]), list(pts[2])] def potential_energy_bare(pos, masses): pe = 0.0 for i in range(3): for j in range(i+1, 3): dx = pos[j][0]-pos[i][0]; dy = pos[j][1]-pos[i][1] r = math.hypot(dx, dy) pe -= masses[i]*masses[j]/r return pe def make_random_ic(rng, masses, alpha, scale=1.0): """Random triangle, random velocity directions scaled by alpha times the virial speed, net momentum removed. Returns pos, vel, actual E.""" pos = sample_triangle(rng, scale=scale) PE = potential_energy_bare(pos, masses) M_tot = sum(masses) v0 = math.sqrt(2.0 * (-PE) / M_tot) thetas = rng.uniform(0, 2*math.pi, size=3) vel = [[alpha*v0*math.cos(th), alpha*v0*math.sin(th)] for th in thetas] px = sum(masses[i]*vel[i][0] for i in range(3)) / M_tot py = sum(masses[i]*vel[i][1] for i in range(3)) / M_tot for i in range(3): vel[i][0] -= px vel[i][1] -= py E = energy(pos, vel, masses, eps2=0.0) return pos, vel, E masses_std = [1.0, 1.0, 1.0] N_RAND_LYAP = 40 T_RAND_LYAP = 35.0 rand_lams = [] rand_Es = [] for k in range(N_RAND_LYAP): pos, vel, E = make_random_ic(RNG_TRIANGLES, masses_std, alpha=0.5) lam, running, mind = lyapunov_run(pos, vel, masses_std, t_total=T_RAND_LYAP, rng=RNG_PERTURB) if math.isfinite(lam): rand_lams.append(lam) rand_Es.append(E) rand_lams = np.array(rand_lams) rand_Es = np.array(rand_Es) rand_lams_median = float(np.median(rand_lams)) print(f"{len(rand_lams)} of {N_RAND_LYAP} random bound triangles (alpha=0.5, " f"i.e. half the virial speed) produced a finite exponent over t={T_RAND_LYAP}.") print(f"Mean E of this batch: {rand_Es.mean():.4f} (SD {rand_Es.std(ddof=1):.4f})") print(f"Mean Lyapunov exponent: {rand_lams.mean():.5f} " f"(SD {rand_lams.std(ddof=1):.5f}, SEM {rand_lams.std(ddof=1)/math.sqrt(len(rand_lams)):.5f})") print(f"Median Lyapunov exponent: {rand_lams_median:.5f}") print(f"Range across the {len(rand_lams)} trials: {rand_lams.min():.4f} to {rand_lams.max():.4f}") n_above_2x_median = int(np.sum(rand_lams > 2 * rand_lams_median)) print(f"{n_above_2x_median} of {len(rand_lams)} trials exceed twice the median: a heavy right") print("tail driven by the rare unusually close triple encounter, not a well-behaved") print("bell curve. The mean is pulled well above the median by exactly those few") print("trials; we use the MEDIAN as the headline 'typical' exponent below, and") print("report the mean alongside it rather than hiding the skew.") print(f"Lyapunov time from the median: {1.0/rand_lams_median:.3f} time units.") print(f"(elapsed {time.perf_counter()-t0:.1f} s)") LYAP_STANDARD_MEAN = rand_lams.mean() LYAP_STANDARD_SEM = rand_lams.std(ddof=1) / math.sqrt(len(rand_lams)) LYAP_STANDARD_MEDIAN = rand_lams_median LYAP_PYTH_MEAN = pyth_lams_f.mean() LYAP_PYTH_SEM = pyth_lams_f.std(ddof=1) / math.sqrt(len(pyth_lams_f)) LYAP_FIG8_MEAN = fig8_lams.mean() LYAP_FIG8_SEM = fig8_lams.std(ddof=1) / math.sqrt(len(fig8_lams)) # --------------------------------------------------------------------------- # SECTION 4: PREDICTION HORIZON # --------------------------------------------------------------------------- section("PREDICTION HORIZON: FROM A LYAPUNOV EXPONENT TO A USABLE TIME LIMIT") print("Horizon formula: t_horizon = (1/lambda) * ln(L / delta),") print("where L is the system's own size scale and delta is the initial") print("uncertainty in the same units. Using the random-triple MEDIAN") print(f"exponent lambda = {LYAP_STANDARD_MEDIAN:.5f} (Lyapunov time " f"{1.0/LYAP_STANDARD_MEDIAN:.3f} time units) and L = 1 (the sampling scale):") print() print(f"{'delta (fraction of L)':>26} | {'t_horizon (time units)':>24} | {'t_horizon / Lyapunov time':>26}") horizon_table = [] for delta in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12]: th = (1.0/LYAP_STANDARD_MEDIAN) * math.log(1.0/delta) horizon_table.append((delta, th)) print(f"{delta:>26.0e} | {th:>24.2f} | {th*LYAP_STANDARD_MEDIAN:>26.2f}") print() print("Reading it: even an initial uncertainty of one part in a trillion") print(f"(delta=1e-12) buys only {horizon_table[-1][1]*LYAP_STANDARD_MEDIAN:.1f} " "Lyapunov times, because the horizon grows only as the LOGARITHM of your") print("measurement precision. Ten more digits of precision buys a fixed, small") print("number of extra Lyapunov times, not ten times longer.") # --------------------------------------------------------------------------- # SECTION 5: MONTE CARLO -- EJECTION FRACTION VS INITIAL ENERGY (adaptive step) # --------------------------------------------------------------------------- section("MONTE CARLO: EJECTION FRACTION AS A FUNCTION OF INITIAL ENERGY") R_ESC = 20.0 T_MAX_MC = 20.0 MAX_STEPS_MC = 400000 def check_escape(pos, vel, masses, r_esc): for k in range(3): others = [j for j in range(3) if j != k] i, j = others mi, mj = masses[i], masses[j] M = mi + mj comx = (mi*pos[i][0] + mj*pos[j][0]) / M comy = (mi*pos[i][1] + mj*pos[j][1]) / M dx = pos[k][0]-comx; dy = pos[k][1]-comy dist = math.hypot(dx, dy) if dist > r_esc: comvx = (mi*vel[i][0] + mj*vel[j][0]) / M comvy = (mi*vel[i][1] + mj*vel[j][1]) / M rvx = vel[k][0]-comvx; rvy = vel[k][1]-comvy speed2 = rvx*rvx + rvy*rvy E_rel = 0.5*speed2 - M/dist radial = dx*rvx + dy*rvy if E_rel > 0 and radial > 0: return k return -1 def run_trial(rng, masses, alpha, eps2, t_max=T_MAX_MC, r_esc=R_ESC, eta=ETA_DEFAULT, dt_max=DT_MAX_DEFAULT): pos, vel, E0 = make_random_ic(rng, masses, alpha) acc = accelerations(pos, masses, eps2) t = 0.0 ejected_at = None min_dist = float("inf") ok = True steps = 0 while t < t_max: dt, rmin = adaptive_dt(pos, masses, eta, dt_max) if rmin < min_dist: min_dist = rmin if t + dt > t_max: dt = t_max - t pos, vel, acc = yoshida4_step(pos, vel, acc, masses, dt, eps2) t += dt steps += 1 if not all_finite(pos, vel): ok = False break k = check_escape(pos, vel, masses, r_esc) if k >= 0: ejected_at = t break if steps > MAX_STEPS_MC: ok = False break energy_final = energy(pos, vel, masses, eps2=0.0) if ok else float("nan") drift = abs((energy_final-E0)/E0) if ok and E0 != 0 else float("nan") return { "E": E0, "ejected": ejected_at is not None, "t_eject": ejected_at, "min_dist": min_dist, "energy_drift": drift, "ok": ok, "steps": steps, } ALPHAS = [0.3, 0.5, 0.7, 0.9, 1.1, 1.3] N_TRIALS_PER_ALPHA = 170 t0 = time.perf_counter() mc_results = {} convergence_track = None for ai, alpha in enumerate(ALPHAS): trials = [run_trial(RNG_MC, masses_std, alpha, SAFETY_EPS2) for _ in range(N_TRIALS_PER_ALPHA)] mc_results[alpha] = trials n_ok = sum(1 for t in trials if t["ok"]) n_ejected = sum(1 for t in trials if t["ok"] and t["ejected"]) frac = n_ejected / n_ok if n_ok else float("nan") se = math.sqrt(frac*(1-frac)/n_ok) if n_ok else float("nan") mean_E = np.mean([t["E"] for t in trials]) mean_drift = np.mean([t["energy_drift"] for t in trials if t["ok"] and math.isfinite(t["energy_drift"])]) print(f"alpha={alpha:.1f} mean E={mean_E:+.4f} n_ok={n_ok}/{N_TRIALS_PER_ALPHA} " f"ejected={n_ejected} fraction={frac:.4f} +/- {se:.4f} " f"mean |energy drift|={mean_drift:.2e}") if alpha == 0.9: convergence_track = trials print(f"(Monte Carlo elapsed {time.perf_counter()-t0:.1f} s, " f"{len(ALPHAS)*N_TRIALS_PER_ALPHA} trials total)") running_frac = [] n_ejected_running = 0 n_ok_running = 0 for i, t in enumerate(convergence_track, start=1): if t["ok"]: n_ok_running += 1 if t["ejected"]: n_ejected_running += 1 frac = n_ejected_running/n_ok_running if n_ok_running else float("nan") se = math.sqrt(frac*(1-frac)/n_ok_running) if n_ok_running else float("nan") running_frac.append((i, frac, se)) print("\nConvergence of the ejection-fraction estimate as trials accumulate, alpha=0.9:") for i in [10, 25, 50, 100, 150, N_TRIALS_PER_ALPHA]: if i <= len(running_frac): _, f, se = running_frac[i-1] print(f" after {i:>4} trials: fraction = {f:.4f} +/- {se:.4f}") # --------------------------------------------------------------------------- # SECTION 6: SENSITIVITY TO MODELLING CHOICES # --------------------------------------------------------------------------- section("SENSITIVITY: WHAT CHANGING ETA AND THE SAFETY LENGTH DO") t0 = time.perf_counter() print(f"\nDefault: eta={ETA_DEFAULT}, safety length={math.sqrt(SAFETY_EPS2):.0e}.") print("Re-measuring the alpha=0.9 ejection fraction at a coarser and a finer eta:") N_SENS = 90 for eta_test, label in [(0.02, "4x coarser"), (0.00125, "4x finer")]: trials = [run_trial(RNG_MC, masses_std, 0.9, SAFETY_EPS2, eta=eta_test) for _ in range(N_SENS)] n_ok = sum(1 for t in trials if t["ok"]) n_ej = sum(1 for t in trials if t["ok"] and t["ejected"]) f = n_ej/n_ok se = math.sqrt(f*(1-f)/n_ok) mean_steps = np.mean([t["steps"] for t in trials if t["ok"]]) print(f" eta={eta_test:<8} ({label}): fraction={f:.3f} +/- {se:.3f} " f"(n_ok={n_ok}/{N_SENS}, mean steps/trial={mean_steps:.0f})") print(f"\nRe-measuring with the safety length changed by 100x each way " f"({N_SENS} trials at alpha=0.9, default eta):") for eps_test, label in [((1e-4) ** 2, "100x larger"), ((1e-8) ** 2, "100x smaller")]: trials = [run_trial(RNG_MC, masses_std, 0.9, eps_test) for _ in range(N_SENS)] n_ok = sum(1 for t in trials if t["ok"]) n_ej = sum(1 for t in trials if t["ok"] and t["ejected"]) f = n_ej/n_ok se = math.sqrt(f*(1-f)/n_ok) print(f" safety length={math.sqrt(eps_test):.0e} ({label}): fraction={f:.3f} +/- {se:.3f} (n_ok={n_ok}/{N_SENS})") print(f"(sensitivity section elapsed {time.perf_counter()-t0:.1f} s)") # --------------------------------------------------------------------------- # SECTION 7: FIGURE DATA # --------------------------------------------------------------------------- section("FIGURE DATA") print("\nFIGURE DATA A: figure-eight orbit trajectory, one period, for the") print("shape figure. (x1,y1,x2,y2,x3,y3) sampled every 1200 steps:") for (t, p, v) in hist_f8[::6]: print(f" t={t:.4f} b1=({p[0][0]:+.5f},{p[0][1]:+.5f}) " f"b2=({p[1][0]:+.5f},{p[1][1]:+.5f}) b3=({p[2][0]:+.5f},{p[2][1]:+.5f})") print("\nFIGURE DATA B: running Lyapunov-exponent estimate vs time for one") print("representative random chaotic trial (alpha=0.5):") pos_demo, vel_demo, E_demo = make_random_ic(RNG_PERTURB, masses_std, alpha=0.5) lam_demo, running_demo, mind_demo = lyapunov_run(pos_demo, vel_demo, masses_std, t_total=30.0, rng=RNG_PERTURB) for (t, lam_run) in running_demo: print(f" t={t:6.2f} running lambda estimate={lam_run:.5f}") print(f"Final running estimate at t={running_demo[-1][0]:.1f}: " f"lambda={running_demo[-1][1]:.5f} (single-trial demo; the headline " f"number above is the {N_RAND_LYAP}-trial median).") print("\nFIGURE DATA C: ejection fraction vs mean energy, with standard errors") print("(x=mean E of the bin, y=fraction ejected, yerr=SE):") for alpha in ALPHAS: trials = mc_results[alpha] n_ok = sum(1 for t in trials if t["ok"]) n_ej = sum(1 for t in trials if t["ok"] and t["ejected"]) frac = n_ej/n_ok se = math.sqrt(frac*(1-frac)/n_ok) mean_E = np.mean([t["E"] for t in trials]) print(f" alpha={alpha:.1f} E={mean_E:+.4f} fraction={frac:.4f} se={se:.4f} n={n_ok}") print("\nFIGURE DATA D: convergence of the alpha=0.9 ejection-fraction estimate") print("as trials accumulate (x=n trials, y=running fraction, yerr=running SE):") for i in range(0, len(running_frac), 10): n, f, se = running_frac[i] print(f" n={n:>4} fraction={f:.4f} se={se:.4f}") n, f, se = running_frac[-1] print(f" n={n:>4} fraction={f:.4f} se={se:.4f}") print("\nFIGURE DATA E: Lyapunov exponent distribution across the " f"{len(rand_lams)} random trials (alpha=0.5), sorted:") for i, lam in enumerate(sorted(rand_lams)): print(f" rank {i+1:>2}/{len(rand_lams)}: lambda={lam:.5f}") # --------------------------------------------------------------------------- # SECTION 8: SUMMARY # --------------------------------------------------------------------------- section("SUMMARY") print(f"Two-body validation: max relative error {max(max_rel_pos_err, max_rel_vel_err):.2e} " f"(tolerance {TWO_BODY_TOL:.0e}) -- {'PASS' if TWO_BODY_PASS else 'FAIL'}") print(f"Figure-eight validation: relative return residual {return_err_f8/scale_f8:.2e} " f"(tolerance {FIG8_TOL:.0e}) -- {'PASS' if FIG8_PASS else 'FAIL'}") print(f"Energy conservation (fig-8, 20 periods, fixed-step symplectic): max relative drift {max_dE_20:.2e}") print(f"Energy conservation (Pythagorean, adaptive step, to t={T_PYTH}): " f"relative drift {abs((E_end_pyth-E0_pyth)/E0_pyth):.2e}") print(f"Figure-eight Lyapunov exponent: {LYAP_FIG8_MEAN:.5f} +/- {LYAP_FIG8_SEM:.5f} (SEM, n={N_FIG8_LYAP_TRIALS})") print(f"Pythagorean-problem exponent: {LYAP_PYTH_MEAN:.5f} +/- {LYAP_PYTH_SEM:.5f} (SEM, n={len(pyth_lams_f)})") print(f"Random-triple mean exponent: {LYAP_STANDARD_MEAN:.5f} +/- {LYAP_STANDARD_SEM:.5f} (SEM, n={len(rand_lams)})") print(f"Random-triple MEDIAN exponent (headline): {LYAP_STANDARD_MEDIAN:.5f}") print(f"Lyapunov time (from median): {1.0/LYAP_STANDARD_MEDIAN:.3f} time units") print(f"Prediction horizon at delta=1e-6 relative uncertainty: " f"{(1.0/LYAP_STANDARD_MEDIAN)*math.log(1e6):.2f} time units") total_elapsed = time.perf_counter() - T_START print(f"\nTOTAL RUNTIME: {total_elapsed:.1f} seconds")