MASTER SEED: 20260914 ============================================================================== VALIDATION 1: TWO-BODY PROBLEM AGAINST THE ANALYTIC KEPLER ORBIT ============================================================================== Analytic Kepler period: T = 4.4428829382 Eccentricity: e = 0.5 Semi-major axis (relative): a = 1.0 Integrator step: dt = 2.221441e-04 (20000 steps/period, fixed) Sampled 25 points across one full orbit and compared the club integrator's relative position and velocity against the closed-form Kepler-equation solution (Newton's method on E - e sin E = M) at each one. Max relative position error over the orbit: 5.058e-12 Max relative velocity error over the orbit: 3.871e-12 Tolerance: 1.0e-08 PASS -- two-body integrator matches the analytic solution. Return-to-start residual after exactly one period: 6.942e-12 (relative to orbit size 5.000e-01: 1.388e-11) Energy at t=0: -0.500000000000 Energy after 1 period: -0.500000000000 Relative energy drift over one period: 7.505e-14 ============================================================================== VALIDATION 2: THE FIGURE-EIGHT THREE-BODY SOLUTION ============================================================================== Literature period (Simo 2001), commonly quoted: T_lit = 6.3259138848 Residual |state(T_lit) - state(0)| integrating for exactly that literature period: 1.517e-03 (relative to orbit scale 1.0000: 1.517e-03), too large to call a clean return. Self-measured return time (parabolic fit around the closest numerical return to the initial state, using this exact code): T = 6.3265551531 Difference from the literature value: +0.000641 time units (+0.0101%), consistent with the widely circulated initial condition being quoted to only 8-9 significant figures. Re-integrating cleanly at the self-measured period, 60000 steps: Residual |state(T) - state(0)|: 1.401e-04 Relative to the orbit's own length scale (1.0000): 1.401e-04 Tolerance: 5.0e-04 (set by the ~8-9 significant figures published for this initial condition, not by the integrator -- its own energy conservation, below, holds to 1e-13). PASS -- figure-eight returns to its initial state within tolerance after one period. Initial energy E0 = -1.2870468378, initial angular momentum L0 = 0.0000000000 Over one period: max |relative energy drift| = 4.624e-14 Over one period: max |linear momentum drift| = 7.289e-14 (started at 0.00e+00) Over one period: max |angular momentum drift| = 5.235e-14 (L0 = 0.000e+00: this choreography has zero net angular momentum by construction, so drift is judged against a scale of 1 when L0 itself is numerically zero: 5.235e-14) Over 20 periods (126.53 time units): max |relative energy drift| = 8.540e-14 (bounded, not growing secularly). ============================================================================== LYAPUNOV EXPONENT MEASUREMENT (Benettin two-trajectory method) ============================================================================== --- 3a. Figure-eight orbit: is the exact periodic solution chaotic? --- 10 independent random perturbation directions, each integrated for 25 periods (158.2 time units). Mean Lyapunov exponent: 0.03736 (SD 0.00444, SEM 0.00140) Lyapunov TIME (1/lambda): 26.77 time units, or 4.23 figure-eight periods. (elapsed 4.3 s) --- 3b. The Pythagorean three-body problem (Burrau 1913 initial data) --- Triangle side lengths (should be 3, 4, 5): 3.000000, 4.000000, 5.000000 Initial energy (all velocities zero): E0 = -12.816667 Fixed-step diagnostic (this is the failure described in the docstring): fixed dt=0.00100, softening 0.001, to t=30: relative energy drift = 3.503e+02 fixed dt=0.00010, softening 0.001, to t=30: relative energy drift = 8.875e+02 Adaptive-step result (eta=0.005, used for every quantitative number below): 10 independent perturbation directions, each run to t=65.0. 10 of 10 completed without hitting the numerical safety cap. Mean Lyapunov exponent: 0.47159 (SD 0.02087, SEM 0.00660) Lyapunov time (1/lambda): 2.120 time units. Closest triple-encounter distance reached, across trials: 0.000413 Reference (unperturbed) trajectory to t=65.0, adaptive eta=0.005: 9141 steps, closest approach 0.000413, relative energy drift = 8.752e-04 (elapsed 3.7 s) --- 3c. Random generic bound triples: the headline Lyapunov exponent --- 40 of 40 random bound triangles (alpha=0.5, i.e. half the virial speed) produced a finite exponent over t=35.0. Mean E of this batch: -2.8601 (SD 0.8995) Mean Lyapunov exponent: 2.54475 (SD 3.02082, SEM 0.47763) Median Lyapunov exponent: 1.32596 Range across the 40 trials: 0.1977 to 14.4680 14 of 40 trials exceed twice the median: a heavy right tail driven by the rare unusually close triple encounter, not a well-behaved bell curve. The mean is pulled well above the median by exactly those few trials; we use the MEDIAN as the headline 'typical' exponent below, and report the mean alongside it rather than hiding the skew. Lyapunov time from the median: 0.754 time units. (elapsed 82.7 s) ============================================================================== PREDICTION HORIZON: FROM A LYAPUNOV EXPONENT TO A USABLE TIME LIMIT ============================================================================== Horizon formula: t_horizon = (1/lambda) * ln(L / delta), where L is the system's own size scale and delta is the initial uncertainty in the same units. Using the random-triple MEDIAN exponent lambda = 1.32596 (Lyapunov time 0.754 time units) and L = 1 (the sampling scale): delta (fraction of L) | t_horizon (time units) | t_horizon / Lyapunov time 1e-02 | 3.47 | 4.61 1e-04 | 6.95 | 9.21 1e-06 | 10.42 | 13.82 1e-08 | 13.89 | 18.42 1e-10 | 17.37 | 23.03 1e-12 | 20.84 | 27.63 Reading it: even an initial uncertainty of one part in a trillion (delta=1e-12) buys only 27.6 Lyapunov times, because the horizon grows only as the LOGARITHM of your measurement precision. Ten more digits of precision buys a fixed, small number of extra Lyapunov times, not ten times longer. ============================================================================== MONTE CARLO: EJECTION FRACTION AS A FUNCTION OF INITIAL ENERGY ============================================================================== alpha=0.3 mean E=-3.1070 n_ok=170/170 ejected=45 fraction=0.2647 +/- 0.0338 mean |energy drift|=1.32e-02 alpha=0.5 mean E=-2.8360 n_ok=170/170 ejected=26 fraction=0.1529 +/- 0.0276 mean |energy drift|=7.32e-03 alpha=0.7 mean E=-2.2112 n_ok=170/170 ejected=25 fraction=0.1471 +/- 0.0272 mean |energy drift|=1.55e-02 alpha=0.9 mean E=-1.6113 n_ok=170/170 ejected=37 fraction=0.2176 +/- 0.0316 mean |energy drift|=4.09e-03 alpha=1.1 mean E=-0.7408 n_ok=170/170 ejected=106 fraction=0.6235 +/- 0.0372 mean |energy drift|=8.37e-03 alpha=1.3 mean E=+0.4675 n_ok=170/170 ejected=138 fraction=0.8118 +/- 0.0300 mean |energy drift|=1.51e-03 (Monte Carlo elapsed 586.0 s, 1020 trials total) Convergence of the ejection-fraction estimate as trials accumulate, alpha=0.9: after 10 trials: fraction = 0.2000 +/- 0.1265 after 25 trials: fraction = 0.1600 +/- 0.0733 after 50 trials: fraction = 0.2000 +/- 0.0566 after 100 trials: fraction = 0.1900 +/- 0.0392 after 150 trials: fraction = 0.1867 +/- 0.0318 after 170 trials: fraction = 0.2176 +/- 0.0316 ============================================================================== SENSITIVITY: WHAT CHANGING ETA AND THE SAFETY LENGTH DO ============================================================================== Default: eta=0.005, safety length=1e-06. Re-measuring the alpha=0.9 ejection fraction at a coarser and a finer eta: eta=0.02 (4x coarser): fraction=0.344 +/- 0.050 (n_ok=90/90, mean steps/trial=1383) eta=0.00125 (4x finer): fraction=0.278 +/- 0.047 (n_ok=90/90, mean steps/trial=23363) Re-measuring with the safety length changed by 100x each way (90 trials at alpha=0.9, default eta): safety length=1e-04 (100x larger): fraction=0.256 +/- 0.046 (n_ok=90/90) safety length=1e-08 (100x smaller): fraction=0.278 +/- 0.047 (n_ok=90/90) (sensitivity section elapsed 259.5 s) ============================================================================== FIGURE DATA ============================================================================== FIGURE DATA A: figure-eight orbit trajectory, one period, for the shape figure. (x1,y1,x2,y2,x3,y3) sampled every 1200 steps: t=0.0000 b1=(+0.97004,-0.24309) b2=(-0.97004,+0.24309) b3=(+0.00000,+0.00000) t=0.1265 b1=(+1.01987,-0.18654) b2=(-0.90082,+0.29455) b3=(-0.11906,-0.10801) t=0.2531 b1=(+1.05343,-0.12802) b2=(-0.80941,+0.33529) b3=(-0.24403,-0.20726) t=0.3796 b1=(+1.07324,-0.06895) b2=(-0.69541,+0.35680) b3=(-0.37783,-0.28784) t=0.5061 b1=(+1.08091,-0.00983) b2=(-0.56355,+0.34971) b3=(-0.51736,-0.33988) t=0.6327 b1=(+1.07709,+0.04929) b2=(-0.42400,+0.30874) b3=(-0.65309,-0.35803) t=0.7592 b1=(+1.06146,+0.10840) b2=(-0.28759,+0.23668) b3=(-0.77387,-0.34508) t=0.8857 b1=(+1.03272,+0.16720) b2=(-0.15977,+0.14249) b3=(-0.87296,-0.30969) t=1.0122 b1=(+0.98864,+0.22465) b2=(-0.03932,+0.03638) b3=(-0.94933,-0.26103) t=1.1388 b1=(+0.92623,+0.27829) b2=(+0.07903,-0.07256) b3=(-1.00526,-0.20573) t=1.2653 b1=(+0.84243,+0.32341) b2=(+0.20149,-0.17576) b3=(-1.04392,-0.14765) t=1.3918 b1=(+0.73574,+0.35239) b2=(+0.33235,-0.26372) b3=(-1.06809,-0.08868) t=1.5184 b1=(+0.60892,+0.35574) b2=(+0.47076,-0.32618) b3=(-1.07969,-0.02956) t=1.6449 b1=(+0.47062,+0.32616) b2=(+0.60907,-0.35572) b3=(-1.07969,+0.02956) t=1.7714 b1=(+0.33222,+0.26366) b2=(+0.73588,-0.35234) b3=(-1.06810,+0.08868) t=1.8980 b1=(+0.20139,+0.17568) b2=(+0.84255,-0.32333) b3=(-1.04394,+0.14765) t=2.0245 b1=(+0.07894,+0.07247) b2=(+0.92634,-0.27820) b3=(-1.00528,+0.20573) t=2.1510 b1=(-0.03938,-0.03649) b2=(+0.98874,-0.22455) b3=(-0.94935,+0.26103) t=2.2776 b1=(-0.15982,-0.14260) b2=(+1.03281,-0.16710) b3=(-0.87299,+0.30969) t=2.4041 b1=(-0.28763,-0.23678) b2=(+1.06154,-0.10829) b3=(-0.77391,+0.34508) t=2.5306 b1=(-0.42402,-0.30884) b2=(+1.07716,-0.04918) b3=(-0.65315,+0.35803) t=2.6572 b1=(-0.56355,-0.34982) b2=(+1.08098,+0.00994) b3=(-0.51743,+0.33989) t=2.7837 b1=(-0.69540,-0.35693) b2=(+1.07331,+0.06906) b3=(-0.37791,+0.28787) t=2.9102 b1=(-0.80938,-0.33544) b2=(+1.05350,+0.12812) b3=(-0.24412,+0.20732) t=3.0367 b1=(-0.90080,-0.29473) b2=(+1.01995,+0.18664) b3=(-0.11915,+0.10809) t=3.1633 b1=(-0.97003,-0.24328) b2=(+0.97012,+0.24318) b3=(-0.00009,+0.00010) t=3.2898 b1=(-1.01987,-0.18674) b2=(+0.90091,+0.29464) b3=(+0.11897,-0.10790) t=3.4163 b1=(-1.05344,-0.12823) b2=(+0.80951,+0.33538) b3=(+0.24393,-0.20715) t=3.5429 b1=(-1.07326,-0.06917) b2=(+0.69554,+0.35690) b3=(+0.37771,-0.28773) t=3.6694 b1=(-1.08094,-0.01004) b2=(+0.56370,+0.34984) b3=(+0.51724,-0.33979) t=3.7959 b1=(-1.07713,+0.04908) b2=(+0.42415,+0.30889) b3=(+0.65297,-0.35797) t=3.9225 b1=(-1.06151,+0.10819) b2=(+0.28775,+0.23686) b3=(+0.77376,-0.34504) t=4.0490 b1=(-1.03278,+0.16699) b2=(+0.15992,+0.14268) b3=(+0.87287,-0.30967) t=4.1755 b1=(-0.98872,+0.22445) b2=(+0.03947,+0.03657) b3=(+0.94925,-0.26102) t=4.3021 b1=(-0.92633,+0.27811) b2=(-0.07887,-0.07239) b3=(+1.00520,-0.20572) t=4.4286 b1=(-0.84255,+0.32325) b2=(-0.20132,-0.17562) b3=(+1.04387,-0.14763) t=4.5551 b1=(-0.73587,+0.35227) b2=(-0.33216,-0.26361) b3=(+1.06803,-0.08866) t=4.6817 b1=(-0.60905,+0.35566) b2=(-0.47058,-0.32612) b3=(+1.07963,-0.02954) t=4.8082 b1=(-0.47073,+0.32612) b2=(-0.60890,-0.35570) b3=(+1.07963,+0.02958) t=4.9347 b1=(-0.33231,+0.26365) b2=(-0.73573,-0.35235) b3=(+1.06804,+0.08870) t=5.0612 b1=(-0.20145,+0.17569) b2=(-0.84243,-0.32336) b3=(+1.04387,+0.14767) t=5.1878 b1=(-0.07898,+0.07248) b2=(-0.92623,-0.27823) b3=(+1.00521,+0.20575) t=5.3143 b1=(+0.03937,-0.03647) b2=(-0.98864,-0.22459) b3=(+0.94928,+0.26105) t=5.4408 b1=(+0.15983,-0.14257) b2=(-1.03273,-0.16714) b3=(+0.87290,+0.30971) t=5.5674 b1=(+0.28766,-0.23675) b2=(-1.06147,-0.10834) b3=(+0.77381,+0.34509) t=5.6939 b1=(+0.42407,-0.30880) b2=(-1.07710,-0.04922) b3=(+0.65303,+0.35802) t=5.8204 b1=(+0.56362,-0.34976) b2=(-1.08092,+0.00990) b3=(+0.51731,+0.33986) t=5.9470 b1=(+0.69548,-0.35683) b2=(-1.07326,+0.06902) b3=(+0.37778,+0.28781) t=6.0735 b1=(+0.80946,-0.33531) b2=(-1.05346,+0.12809) b3=(+0.24399,+0.20723) t=6.2000 b1=(+0.90087,-0.29458) b2=(-1.01991,+0.18660) b3=(+0.11904,+0.10798) t=6.3266 b1=(+0.97010,-0.24312) b2=(-0.97008,+0.24315) b3=(-0.00001,-0.00003) FIGURE DATA B: running Lyapunov-exponent estimate vs time for one representative random chaotic trial (alpha=0.5): t= 0.50 running lambda estimate=3.11104 t= 1.00 running lambda estimate=1.27532 t= 1.50 running lambda estimate=1.09326 t= 2.00 running lambda estimate=4.64169 t= 2.50 running lambda estimate=4.09286 t= 3.00 running lambda estimate=3.56250 t= 3.50 running lambda estimate=3.16947 t= 4.00 running lambda estimate=2.90944 t= 4.50 running lambda estimate=2.79468 t= 5.00 running lambda estimate=2.91266 t= 5.50 running lambda estimate=2.21666 t= 6.00 running lambda estimate=1.85913 t= 6.50 running lambda estimate=1.86460 t= 7.00 running lambda estimate=1.78855 t= 7.50 running lambda estimate=1.71178 t= 8.00 running lambda estimate=1.65487 t= 8.50 running lambda estimate=1.63626 t= 9.00 running lambda estimate=1.74431 t= 9.50 running lambda estimate=1.66086 t= 10.00 running lambda estimate=1.50208 t= 10.50 running lambda estimate=1.41379 t= 11.00 running lambda estimate=1.34998 t= 11.50 running lambda estimate=1.29864 t= 12.00 running lambda estimate=1.25386 t= 12.50 running lambda estimate=1.22411 t= 13.00 running lambda estimate=1.20966 t= 13.50 running lambda estimate=1.24974 t= 14.00 running lambda estimate=0.97872 t= 14.50 running lambda estimate=0.98232 t= 15.00 running lambda estimate=0.98426 t= 15.50 running lambda estimate=0.97451 t= 16.00 running lambda estimate=0.96430 t= 16.50 running lambda estimate=0.95875 t= 17.00 running lambda estimate=0.95062 t= 17.50 running lambda estimate=0.93860 t= 18.00 running lambda estimate=0.84029 t= 18.50 running lambda estimate=0.85674 t= 19.00 running lambda estimate=0.81871 t= 19.50 running lambda estimate=0.77128 t= 20.00 running lambda estimate=0.70761 t= 20.50 running lambda estimate=0.67470 t= 21.00 running lambda estimate=0.73402 t= 21.50 running lambda estimate=0.77041 t= 22.00 running lambda estimate=0.76488 t= 22.50 running lambda estimate=0.78751 t= 23.00 running lambda estimate=0.75835 t= 23.50 running lambda estimate=0.74139 t= 24.00 running lambda estimate=0.72339 t= 24.50 running lambda estimate=0.69721 t= 25.00 running lambda estimate=0.62428 t= 25.50 running lambda estimate=0.66968 t= 26.00 running lambda estimate=0.72905 t= 26.50 running lambda estimate=0.78915 t= 27.00 running lambda estimate=0.67335 t= 27.50 running lambda estimate=0.70316 t= 28.00 running lambda estimate=0.69777 t= 28.50 running lambda estimate=0.69110 t= 29.00 running lambda estimate=0.69343 t= 29.50 running lambda estimate=0.69126 t= 30.00 running lambda estimate=0.68996 Final running estimate at t=30.0: lambda=0.68996 (single-trial demo; the headline number above is the 40-trial median). FIGURE DATA C: ejection fraction vs mean energy, with standard errors (x=mean E of the bin, y=fraction ejected, yerr=SE): alpha=0.3 E=-3.1070 fraction=0.2647 se=0.0338 n=170 alpha=0.5 E=-2.8360 fraction=0.1529 se=0.0276 n=170 alpha=0.7 E=-2.2112 fraction=0.1471 se=0.0272 n=170 alpha=0.9 E=-1.6113 fraction=0.2176 se=0.0316 n=170 alpha=1.1 E=-0.7408 fraction=0.6235 se=0.0372 n=170 alpha=1.3 E=+0.4675 fraction=0.8118 se=0.0300 n=170 FIGURE DATA D: convergence of the alpha=0.9 ejection-fraction estimate as trials accumulate (x=n trials, y=running fraction, yerr=running SE): n= 1 fraction=0.0000 se=0.0000 n= 11 fraction=0.1818 se=0.1163 n= 21 fraction=0.1905 se=0.0857 n= 31 fraction=0.2581 se=0.0786 n= 41 fraction=0.2195 se=0.0646 n= 51 fraction=0.1961 se=0.0556 n= 61 fraction=0.1967 se=0.0509 n= 71 fraction=0.1972 se=0.0472 n= 81 fraction=0.1975 se=0.0442 n= 91 fraction=0.1868 se=0.0409 n= 101 fraction=0.1881 se=0.0389 n= 111 fraction=0.1712 se=0.0358 n= 121 fraction=0.1736 se=0.0344 n= 131 fraction=0.1756 se=0.0332 n= 141 fraction=0.1702 se=0.0316 n= 151 fraction=0.1921 se=0.0321 n= 161 fraction=0.2050 se=0.0318 n= 170 fraction=0.2176 se=0.0316 FIGURE DATA E: Lyapunov exponent distribution across the 40 random trials (alpha=0.5), sorted: rank 1/40: lambda=0.19765 rank 2/40: lambda=0.24105 rank 3/40: lambda=0.36447 rank 4/40: lambda=0.43692 rank 5/40: lambda=0.45959 rank 6/40: lambda=0.47549 rank 7/40: lambda=0.51297 rank 8/40: lambda=0.51834 rank 9/40: lambda=0.63740 rank 10/40: lambda=0.75572 rank 11/40: lambda=0.78262 rank 12/40: lambda=0.80480 rank 13/40: lambda=0.91240 rank 14/40: lambda=0.95088 rank 15/40: lambda=1.02841 rank 16/40: lambda=1.03127 rank 17/40: lambda=1.09221 rank 18/40: lambda=1.14397 rank 19/40: lambda=1.19722 rank 20/40: lambda=1.32244 rank 21/40: lambda=1.32948 rank 22/40: lambda=1.59425 rank 23/40: lambda=1.71657 rank 24/40: lambda=1.91160 rank 25/40: lambda=2.43023 rank 26/40: lambda=2.54833 rank 27/40: lambda=2.66932 rank 28/40: lambda=3.09117 rank 29/40: lambda=3.38874 rank 30/40: lambda=3.57752 rank 31/40: lambda=3.70572 rank 32/40: lambda=3.75719 rank 33/40: lambda=3.76035 rank 34/40: lambda=4.31564 rank 35/40: lambda=4.34365 rank 36/40: lambda=4.44537 rank 37/40: lambda=5.37806 rank 38/40: lambda=5.43935 rank 39/40: lambda=13.05385 rank 40/40: lambda=14.46797 ============================================================================== SUMMARY ============================================================================== Two-body validation: max relative error 5.06e-12 (tolerance 1e-08) -- PASS Figure-eight validation: relative return residual 1.40e-04 (tolerance 5e-04) -- PASS Energy conservation (fig-8, 20 periods, fixed-step symplectic): max relative drift 8.54e-14 Energy conservation (Pythagorean, adaptive step, to t=65.0): relative drift 8.75e-04 Figure-eight Lyapunov exponent: 0.03736 +/- 0.00140 (SEM, n=10) Pythagorean-problem exponent: 0.47159 +/- 0.00660 (SEM, n=10) Random-triple mean exponent: 2.54475 +/- 0.47763 (SEM, n=40) Random-triple MEDIAN exponent (headline): 1.32596 Lyapunov time (from median): 0.754 time units Prediction horizon at delta=1e-6 relative uncertainty: 10.42 time units TOTAL RUNTIME: 964.5 seconds