.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/astro/nbody/plot_01_three_body_sensitivity.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_astro_nbody_plot_01_three_body_sensitivity.py: Poincare's discovery: sensitive dependence in the three-body problem ========================================================================= Studying the three-body problem for King Oscar II's 1889 prize, Poincare found that no general closed-form solution exists at all: two trajectories starting arbitrarily close together can separate at an exponential rate, the first rigorous demonstration of what would later be called deterministic chaos. Lagrange (1772) found one very special three-body configuration that *is* exactly solvable -- three masses at the vertices of an equilateral triangle, rigidly rotating about their common center of mass -- but Routh (1875) showed its linear stability depends on the mass ratios (stable only when one mass strongly dominates, as for the Sun-Jupiter-Trojan system). For masses like 3, 4, and 5 (in :math:`G=1` units), it is deep in the unstable regime, and provides a numerically clean, close-encounter-free setting to watch Poincare's sensitive dependence directly with :class:`~physicskit.astro.nbody.NBodySystem`, without needing to resolve the extreme close encounters of a genuinely colliding configuration. .. GENERATED FROM PYTHON SOURCE LINES 23-29 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.astro.nbody import NBodySystem from physicskit.astro.visualizers import plot_nbody_trajectories .. GENERATED FROM PYTHON SOURCE LINES 30-38 The (unstable) Lagrange equilateral-triangle equilibrium --------------------------------------------------------------- Three masses at the vertices of an equilateral triangle of side ``d``, each moving on a circle about the (mass-weighted) center of mass with the same angular velocity :math:`\Omega=\sqrt{GM_{\rm tot}/d^3}` -- an exact rotating equilibrium for *any* choice of the three masses, a remarkable property special to the equilateral configuration. .. GENERATED FROM PYTHON SOURCE LINES 38-55 .. code-block:: Python masses = np.array([3.0, 4.0, 5.0]) M_tot = masses.sum() d = 1.0 R_geom = d / np.sqrt(3) angles = np.deg2rad([90.0, 210.0, 330.0]) vertices = np.stack([R_geom * np.cos(angles), R_geom * np.sin(angles), np.zeros(3)], axis=1) center_of_mass = (masses[:, None] * vertices).sum(axis=0) / M_tot positions0 = vertices - center_of_mass Omega = np.sqrt(M_tot / d**3) z_hat = np.array([0.0, 0.0, 1.0]) velocities0 = Omega * np.cross(z_hat, positions0) side_lengths = [np.linalg.norm(positions0[i] - positions0[j]) for i, j in [(0, 1), (1, 2), (2, 0)]] print(f"initial side lengths (should all equal d={d}): {np.round(side_lengths, 6)}") .. rst-class:: sphx-glr-script-out .. code-block:: none initial side lengths (should all equal d=1.0): [1. 1. 1.] .. GENERATED FROM PYTHON SOURCE LINES 56-60 Two runs, perturbed by one part in :math:`10^8` ----------------------------------------------------- Identical initial conditions and identical integration, except that the perturbed run's first body starts displaced by :math:`10^{-8}`. .. GENERATED FROM PYTHON SOURCE LINES 60-78 .. code-block:: Python dt, n_steps = 0.001, 16_000 baseline = NBodySystem(positions0, velocities0, masses, softening=1e-6) E0 = baseline.total_energy() history_baseline = baseline.simulate(dt, n_steps) E1 = baseline.total_energy() print(f"energy conservation check: relative drift = {abs((E1 - E0) / E0):.2e} (leapfrog is trustworthy here)") perturbation = 1e-8 positions0_perturbed = positions0.copy() positions0_perturbed[0, 0] += perturbation perturbed = NBodySystem(positions0_perturbed, velocities0.copy(), masses, softening=1e-6) history_perturbed = perturbed.simulate(dt, n_steps) fig1, ax1 = plot_nbody_trajectories(history_baseline) ax1.set_title("Unstable Lagrange triangle: the triangle itself deforms") fig1.tight_layout() .. image-sg:: /api/gallery/astro/nbody/images/sphx_glr_plot_01_three_body_sensitivity_001.png :alt: Unstable Lagrange triangle: the triangle itself deforms :srcset: /api/gallery/astro/nbody/images/sphx_glr_plot_01_three_body_sensitivity_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none energy conservation check: relative drift = 1.18e-06 (leapfrog is trustworthy here) .. GENERATED FROM PYTHON SOURCE LINES 79-88 Exponential separation, then saturation --------------------------------------------- A single part-in-:math:`10^8` initial nudge grows, on a semi-log plot, as a straight line -- exponential divergence -- for as long as the separation stays much smaller than the system's own size, then saturates once it does not: exactly the two-stage signature of deterministic chaos, and exactly what makes long-term three-body prediction impossible in practice, regardless of how precisely the initial conditions are known. .. GENERATED FROM PYTHON SOURCE LINES 88-104 .. code-block:: Python separation = np.linalg.norm(history_baseline[:, 0, :] - history_perturbed[:, 0, :], axis=1) t = np.arange(n_steps + 1) * dt fig2, ax2 = plt.subplots(figsize=(6.5, 4.5)) ax2.semilogy(t, np.maximum(separation, 1e-300)) ax2.set_xlabel("t") ax2.set_ylabel(f"separation of body 1 (initial perturbation = {perturbation:.0e})") ax2.set_title("Sensitive dependence on initial conditions") fig2.tight_layout() growth_mask = (t > 1.0) & (t < 7.0) & (separation > 0) slope = np.polyfit(t[growth_mask], np.log(separation[growth_mask]), 1)[0] print(f"\nseparation grows from {separation[growth_mask][0]:.2e} to {separation[growth_mask][-1]:.2e} over t in [1, 7]") print(f"fitted exponential growth rate (Lyapunov-like exponent): {slope:.3f} per unit time") plt.show() .. image-sg:: /api/gallery/astro/nbody/images/sphx_glr_plot_01_three_body_sensitivity_002.png :alt: Sensitive dependence on initial conditions :srcset: /api/gallery/astro/nbody/images/sphx_glr_plot_01_three_body_sensitivity_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none separation grows from 3.85e-08 to 5.52e-02 over t in [1, 7] fitted exponential growth rate (Lyapunov-like exponent): 2.414 per unit time .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.576 seconds) .. _sphx_glr_download_api_gallery_astro_nbody_plot_01_three_body_sensitivity.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_three_body_sensitivity.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_three_body_sensitivity.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_three_body_sensitivity.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_