.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/chaos/advanced/plot_symplectic_integrators.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_chaos_advanced_plot_symplectic_integrators.py: Symplectic vs. Non-Symplectic Integrators: Long-Horizon Energy Drift ======================================================================== For conservative (Hamiltonian) systems, the choice of integrator matters far more over long integration horizons than local per-step accuracy suggests. This example integrates a simple, exactly-solvable Kepler orbit with three integrators from :mod:`physicskit.chaos.core.integrators` -- RK4, the 2nd-order symplectic leapfrog, and the 4th-order symplectic Yoshida scheme -- and compares how their energy-conservation error behaves as the integration horizon grows. RK4's error grows *secularly* (roughly proportional to integration time), while the symplectic integrators' error stays *bounded* regardless of how long you integrate, because they exactly preserve phase-space volume at every step: over a handful of orbits RK4 actually has the smallest error of the three, but by a few dozen orbits it has already grown past Yoshida4's, and keeps growing indefinitely from there. A second figure makes the abstract "energy drift" number concrete: since a Kepler orbit's specific energy fixes its semi-major axis, tracking the aphelion distance orbit by orbit shows RK4's secular energy error directly shrinking the orbit over time, while Yoshida4's bounded error keeps it flat. .. GENERATED FROM PYTHON SOURCE LINES 22-30 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from numba import njit from physicskit.chaos.core.integrators import leapfrog_integrate, rk4_integrate, yoshida4_integrate .. GENERATED FROM PYTHON SOURCE LINES 31-45 A Kepler orbit -------------- A planet of negligible mass orbiting a unit point mass under Newtonian gravity, in units where G*M = 1, so the position :math:`\mathbf{r}=(x,y)` obeys .. math:: \ddot{\mathbf{r}} = -\frac{\mathbf{r}}{|\mathbf{r}|^3}, with conserved specific orbital energy :math:`E = \tfrac{1}{2}|\dot{ \mathbf{r}}|^2 - 1/|\mathbf{r}|`. Gravity here depends only on position, so this is a *separable* Hamiltonian -- exactly the kind of system the symplectic integrators are built for. .. GENERATED FROM PYTHON SOURCE LINES 45-75 .. code-block:: Python @njit(cache=True) def _kepler_force(pos, t, params): r = np.sqrt(pos[0] ** 2 + pos[1] ** 2) return -pos / r**3 @njit(cache=True) def _kepler_rhs(state, t, params): x, y, vx, vy = state[0], state[1], state[2], state[3] r = np.sqrt(x * x + y * y) out = np.empty(4) out[0], out[1] = vx, vy out[2], out[3] = -x / r**3, -y / r**3 return out def _kepler_energy(pos, vel): r = np.hypot(pos[:, 0], pos[:, 1]) v2 = vel[:, 0] ** 2 + vel[:, 1] ** 2 return 0.5 * v2 - 1.0 / r # An eccentric orbit (e = 0.6, semi-major axis a = 1, period 2*pi), starting # at perihelion (r_min = a*(1-e) = 0.4) with the vis-viva speed # v_peri = sqrt(2/r_min - 1/a) = 2.0. pos0 = np.array([0.4, 0.0]) vel0 = np.array([0.0, 2.0]) params = np.array([0.0]) dt = 0.01 .. GENERATED FROM PYTHON SOURCE LINES 76-80 Integrate over successively longer horizons ---------------------------------------------- The orbital period here is ``2*pi ~= 6.283``; these horizons span roughly 10 to 10,000 orbits. .. GENERATED FROM PYTHON SOURCE LINES 80-97 .. code-block:: Python period = 2.0 * np.pi horizons = [int(n_orbits * period / dt) for n_orbits in (10, 100, 1000, 10000)] drift_rk4, drift_lf, drift_y4 = [], [], [] for n_steps in horizons: _, states = rk4_integrate(_kepler_rhs, np.concatenate([pos0, vel0]), 0.0, dt, n_steps, params) e_rk4 = 0.5 * (states[:, 2] ** 2 + states[:, 3] ** 2) - 1.0 / np.hypot(states[:, 0], states[:, 1]) drift_rk4.append(np.max(np.abs(e_rk4 - e_rk4[0]))) _, pos_lf, vel_lf = leapfrog_integrate(_kepler_force, pos0, vel0, 0.0, dt, n_steps, params) e_lf = _kepler_energy(pos_lf, vel_lf) drift_lf.append(np.max(np.abs(e_lf - e_lf[0]))) _, pos_y4, vel_y4 = yoshida4_integrate(_kepler_force, pos0, vel0, 0.0, dt, n_steps, params) e_y4 = _kepler_energy(pos_y4, vel_y4) drift_y4.append(np.max(np.abs(e_y4 - e_y4[0]))) .. GENERATED FROM PYTHON SOURCE LINES 98-100 Plot ---- .. GENERATED FROM PYTHON SOURCE LINES 100-110 .. code-block:: Python t_max = np.array(horizons) * dt fig, ax = plt.subplots(figsize=(7, 5)) ax.loglog(t_max, drift_rk4, "o-", label="RK4 (not symplectic)") ax.loglog(t_max, drift_lf, "s-", label="Leapfrog (2nd-order symplectic)") ax.loglog(t_max, drift_y4, "^-", label="Yoshida (4th-order symplectic)") ax.set_xlabel("integration horizon (t_max)") ax.set_ylabel("max |energy drift|") ax.set_title("RK4 starts more accurate, but grows past Yoshida4 within ~50 orbits") ax.legend() .. image-sg:: /api/gallery/chaos/advanced/images/sphx_glr_plot_symplectic_integrators_001.png :alt: RK4 starts more accurate, but grows past Yoshida4 within ~50 orbits :srcset: /api/gallery/chaos/advanced/images/sphx_glr_plot_symplectic_integrators_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 111-122 What the drift actually does to the orbit: a shrinking aphelion ------------------------------------------------------------------------ The log-log plot above only shows a scalar (the worst energy error) at a handful of horizons; it doesn't show what that error actually does to the orbit. A negative specific energy fixes the semi-major axis via :math:`E=-1/(2a)`, so any energy drift directly rescales the orbit's size; tracking the aphelion distance (the local maxima of :math:`r(t)`, one per orbit) makes that rescaling visible one orbit at a time. Over 1000 orbits at the same ``dt`` used above, RK4's aphelion visibly creeps inward as its secular energy error accumulates, while Yoshida4's stays flat -- bounded noise at the :math:`10^{-6}` level, with no trend at all. .. GENERATED FROM PYTHON SOURCE LINES 122-154 .. code-block:: Python n_orbits_shape = 1000 n_steps_shape = int(n_orbits_shape * period / dt) _, states_rk4_shape = rk4_integrate(_kepler_rhs, np.concatenate([pos0, vel0]), 0.0, dt, n_steps_shape, params) _, pos_y4_shape, vel_y4_shape = yoshida4_integrate(_kepler_force, pos0, vel0, 0.0, dt, n_steps_shape, params) def _aphelion_per_orbit(pos): """Distance from the origin at each local maximum of r(t) (one per orbit).""" r = np.hypot(pos[:, 0], pos[:, 1]) is_max = (r[1:-1] > r[:-2]) & (r[1:-1] > r[2:]) return r[1:-1][is_max] aphelion_rk4 = _aphelion_per_orbit(states_rk4_shape[:, :2]) aphelion_y4 = _aphelion_per_orbit(pos_y4_shape) print( f"Aphelion distance (true value = 1.6): RK4 drifts from {aphelion_rk4[0]:.6f} to " f"{aphelion_rk4[-1]:.6f}; Yoshida4 stays within {np.ptp(aphelion_y4):.1e} of {aphelion_y4[0]:.6f}" ) fig2, ax2 = plt.subplots(figsize=(8, 5)) ax2.plot(np.arange(1, len(aphelion_rk4) + 1), aphelion_rk4, color="crimson", label="RK4 (not symplectic)") ax2.plot(np.arange(1, len(aphelion_y4) + 1), aphelion_y4, color="steelblue", label="Yoshida4 (4th-order symplectic)") ax2.axhline(1.6, color="gray", lw=0.8, ls="--", label="exact aphelion = 1.6") ax2.set_xlabel("orbit number") ax2.set_ylabel("aphelion distance") ax2.set_title(f"Same Kepler orbit, same dt: RK4's aphelion creeps inward over {n_orbits_shape} orbits") ax2.legend() fig2.tight_layout() plt.show() .. image-sg:: /api/gallery/chaos/advanced/images/sphx_glr_plot_symplectic_integrators_002.png :alt: Same Kepler orbit, same dt: RK4's aphelion creeps inward over 1000 orbits :srcset: /api/gallery/chaos/advanced/images/sphx_glr_plot_symplectic_integrators_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Aphelion distance (true value = 1.6): RK4 drifts from 1.600000 to 1.599956; Yoshida4 stays within 2.9e-06 of 1.599998 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.596 seconds) .. _sphx_glr_download_api_gallery_chaos_advanced_plot_symplectic_integrators.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_symplectic_integrators.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_symplectic_integrators.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_symplectic_integrators.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_