.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/classical/newtonian/plot_02_kepler_precession.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_classical_newtonian_plot_02_kepler_precession.py: Kepler orbits and perihelion precession =========================================== :class:`~physicskit.classical.systems.newtonian.KeplerSystem` is the planar, reduced 2-body problem in the center-of-mass frame, a separable Hamiltonian in Cartesian coordinates :math:`q=(x,y)`, :math:`p=(p_x,p_y)`, .. math:: H(q, p) = \frac{|p|^2}{2\mu} + V(r), \qquad r = |q|, \qquad V(r) = -\frac{k}{r} + \frac{c_\mathrm{pn}}{r^3} , with :math:`k` the gravitational coupling, :math:`\mu` the reduced mass, and :math:`c_\mathrm{pn}` an effective post-Newtonian correction that reproduces the qualitative apsidal (perihelion) precession General Relativity predicts for Mercury's orbit. Both terms are purely radial, so the exactly conserved Laplace-Runge-Lenz vector of the pure Kepler problem, .. math:: \mathbf{A} = \mathbf{p} \times L - \mu k \hat{\mathbf{r}}, \qquad L = x p_y - y p_x , slowly precesses once :math:`c_\mathrm{pn} \neq 0`, even though the system remains conservative and separable. Shows a pure 1/r orbit's closed ellipse next to the precessing rosette produced by the post-Newtonian correction, then verifies the key physical claim: total energy stays conserved (to machine precision, via Yoshida4) even as the Laplace-Runge-Lenz vector's direction slowly rotates -- the signature of apsidal precession, e.g. Mercury's orbit. .. GENERATED FROM PYTHON SOURCE LINES 37-43 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.classical.systems.newtonian import KeplerSystem from physicskit.classical.utils.conservation import relative_energy_drift .. GENERATED FROM PYTHON SOURCE LINES 44-46 Pure ellipse vs. precessing rosette --------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 46-65 .. code-block:: Python pure = KeplerSystem.from_orbital_elements(a=1.0, e=0.5) perturbed = KeplerSystem.from_orbital_elements(a=1.0, e=0.5, c_pn=0.008) res_pure = pure.integrate((0, 40), dt=1e-3, method="yoshida4") res_pert = perturbed.integrate((0, 40), dt=1e-3, method="yoshida4") fig1, axes = plt.subplots(1, 2, figsize=(9, 4.2)) axes[0].plot(res_pure.q[:, 0], res_pure.q[:, 1], color="steelblue", lw=0.8) axes[0].plot(0, 0, "o", color="gold") axes[0].set_title("Pure 1/r: closed ellipse") axes[0].set_aspect("equal") axes[1].plot(res_pert.q[:, 0], res_pert.q[:, 1], color="firebrick", lw=0.6) axes[1].plot(0, 0, "o", color="gold") axes[1].set_title("With c_pn: precessing rosette") axes[1].set_aspect("equal") fig1.tight_layout() .. image-sg:: /api/gallery/classical/newtonian/images/sphx_glr_plot_02_kepler_precession_001.png :alt: Pure 1/r: closed ellipse, With c_pn: precessing rosette :srcset: /api/gallery/classical/newtonian/images/sphx_glr_plot_02_kepler_precession_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 66-72 Energy stays flat while the LRL vector precesses ----------------------------------------------------- Both terms in ``V(r)`` are still purely radial (velocity-independent), so the system remains conservative and separable -- Yoshida4 keeps the energy to machine precision even though the LRL vector itself is no longer a constant of motion. .. GENERATED FROM PYTHON SOURCE LINES 72-95 .. code-block:: Python system = KeplerSystem.from_orbital_elements(a=1.0, e=0.5, c_pn=0.008) lrl0 = system.lrl_vector() result = system.integrate((0, 200), dt=1e-3, method="yoshida4") lrl_vecs = np.array([system.lrl_vector(q, p) for q, p in zip(result.q, result.p)]) angle = np.arctan2(lrl_vecs[:, 1], lrl_vecs[:, 0]) - np.arctan2(lrl0[1], lrl0[0]) drift = relative_energy_drift(result.energy) print(f"total perihelion precession over t=200: {np.unwrap(angle)[-1]:.4f} rad") print(f"max relative energy drift: {np.max(drift):.3e}") fig2, axes2 = plt.subplots(1, 2, figsize=(9, 3.8)) axes2[0].plot(result.t, np.unwrap(angle), color="firebrick") axes2[0].set_xlabel("t") axes2[0].set_ylabel("LRL precession angle (rad)") axes2[0].set_title("Perihelion advance") axes2[1].semilogy(result.t, drift, color="steelblue") axes2[1].set_xlabel("t") axes2[1].set_ylabel("|H(t) - H(0)| / |H(0)|") axes2[1].set_title("Energy still conserved") fig2.tight_layout() .. image-sg:: /api/gallery/classical/newtonian/images/sphx_glr_plot_02_kepler_precession_002.png :alt: Perihelion advance, Energy still conserved :srcset: /api/gallery/classical/newtonian/images/sphx_glr_plot_02_kepler_precession_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none total perihelion precession over t=200: -6.1370 rad max relative energy drift: 1.823e-11 .. GENERATED FROM PYTHON SOURCE LINES 96-99 This is the same physics as :doc:`plot_01_cannonball` (a conservative central force integrated symplectically), generalized from uniform gravity to an inverse-square field with a relativistic correction. .. GENERATED FROM PYTHON SOURCE LINES 99-101 .. code-block:: Python plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.415 seconds) .. _sphx_glr_download_api_gallery_classical_newtonian_plot_02_kepler_precession.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_02_kepler_precession.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_02_kepler_precession.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_02_kepler_precession.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_