.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/astro/nbody/plot_04_symplectic_long_term_stability.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_04_symplectic_long_term_stability.py: Symplectic integration and long-term N-body stability ========================================================= Wisdom and Holman (1991) showed that symplectic integrators -- of which the kick-drift-kick leapfrog scheme is the simplest example -- keep a conservative system's numerically computed energy oscillating within a small bound *forever*, rather than drifting secularly away as an ordinary (non-symplectic) integrator's does. This is exactly what made simulations at the scale of the Millennium Simulation's billions of particles, evolved over cosmological timescales, numerically trustworthy at all. This example integrates the same eccentric two-body orbit with :class:`~physicskit.astro.nbody.NBodySystem`'s symplectic :func:`~physicskit.astro.nbody.leapfrog_step` for a thousand orbital periods, alongside a hand-written non-symplectic explicit-Euler stepper using the *same* force law (:func:`~physicskit.astro.nbody.gravitational_acceleration`), and compares their energy error growth directly. .. GENERATED FROM PYTHON SOURCE LINES 22-28 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.astro.nbody import NBodySystem, gravitational_acceleration from physicskit.astro.orbital_mechanics import orbital_period, state_from_orbital_elements .. GENERATED FROM PYTHON SOURCE LINES 29-31 An eccentric two-body orbit, run for 1000 periods --------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 31-53 .. code-block:: Python mu_total, m1, m2 = 1.0, 0.3, 0.7 a, e = 1.0, 0.5 r_vec, v_vec = state_from_orbital_elements(a, e, 0.0, 0.0, 0.0, 0.0, mu_total) positions0 = np.array([-(m2 / mu_total) * r_vec, (m1 / mu_total) * r_vec]) velocities0 = np.array([-(m2 / mu_total) * v_vec, (m1 / mu_total) * v_vec]) masses = np.array([m1, m2]) T_period = orbital_period(a, mu_total) steps_per_period = 200 n_periods = 1000 dt = T_period / steps_per_period n_steps = steps_per_period * n_periods def total_energy(positions, velocities, masses): kinetic = 0.5 * np.sum(masses[:, None] * velocities**2) r = np.linalg.norm(positions[0] - positions[1]) potential = -masses[0] * masses[1] / r return kinetic + potential .. GENERATED FROM PYTHON SOURCE LINES 54-56 The symplectic run ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 56-65 .. code-block:: Python symplectic = NBodySystem(positions0.copy(), velocities0.copy(), masses) E0 = symplectic.total_energy() energy_symplectic = np.empty(n_periods + 1) energy_symplectic[0] = E0 for p in range(n_periods): for _ in range(steps_per_period): symplectic.step(dt) energy_symplectic[p + 1] = symplectic.total_energy() .. GENERATED FROM PYTHON SOURCE LINES 66-72 The non-symplectic run: hand-written explicit Euler, same force law -------------------------------------------------------------------------- ``pos += v*dt; v += a*dt`` using the very same :func:`~physicskit.astro.nbody.gravitational_acceleration`, so any difference in long-term behavior is due entirely to the integration scheme, not the physics. .. GENERATED FROM PYTHON SOURCE LINES 72-82 .. code-block:: Python pos_e, vel_e = positions0.copy(), velocities0.copy() energy_euler = np.empty(n_periods + 1) energy_euler[0] = total_energy(pos_e, vel_e, masses) for p in range(n_periods): for _ in range(steps_per_period): acc = gravitational_acceleration(pos_e, masses) pos_e = pos_e + vel_e * dt vel_e = vel_e + acc * dt energy_euler[p + 1] = total_energy(pos_e, vel_e, masses) .. GENERATED FROM PYTHON SOURCE LINES 83-85 Bounded oscillation vs. secular drift -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 85-103 .. code-block:: Python periods = np.arange(n_periods + 1) rel_err_symplectic = np.abs((energy_symplectic - E0) / E0) rel_err_euler = np.abs((energy_euler - E0) / E0) print(f"after {n_periods} periods:") print(f" symplectic (leapfrog) relative energy error: {rel_err_symplectic[-1]:.2e} (max over the run: {rel_err_symplectic.max():.2e})") print(f" explicit Euler relative energy error: {rel_err_euler[-1]:.2e}") fig, ax = plt.subplots(figsize=(7, 4.5)) ax.semilogy(periods, np.maximum(rel_err_symplectic, 1e-16), color="steelblue", label="symplectic leapfrog (bounded)") ax.semilogy(periods, np.maximum(rel_err_euler, 1e-16), color="firebrick", label="explicit Euler (secular drift)") ax.set_xlabel("orbital periods elapsed") ax.set_ylabel("|relative energy error|") ax.set_title("Bounded oscillation vs. secular drift over 1000 orbits") ax.legend() fig.tight_layout() plt.show() .. image-sg:: /api/gallery/astro/nbody/images/sphx_glr_plot_04_symplectic_long_term_stability_001.png :alt: Bounded oscillation vs. secular drift over 1000 orbits :srcset: /api/gallery/astro/nbody/images/sphx_glr_plot_04_symplectic_long_term_stability_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none after 1000 periods: symplectic (leapfrog) relative energy error: 2.67e-03 (max over the run: 2.69e-03) explicit Euler relative energy error: 1.07e+00 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.142 seconds) .. _sphx_glr_download_api_gallery_astro_nbody_plot_04_symplectic_long_term_stability.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_04_symplectic_long_term_stability.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_04_symplectic_long_term_stability.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_04_symplectic_long_term_stability.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_