.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/relativity/schwarzschild/plot_orbit_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_relativity_schwarzschild_plot_orbit_precession.py: Perihelion precession: the anomaly that first confirmed General Relativity ================================================================================ Mercury's orbit precesses by about 43 arcseconds per century more than Newtonian gravity (including the perturbations of the other planets) can account for -- a puzzle known since Le Verrier's 1859 analysis. Einstein's 1915 calculation of exactly this excess, from General Relativity alone with no free parameters, was his first observational triumph and a major reason the theory was taken seriously. Unlike a Newtonian ellipse, a bound Schwarzschild orbit of semi-major axis :math:`a` and eccentricity :math:`e` advances its periapsis by, at leading post-Newtonian order, .. math:: \Delta\phi = \frac{6\pi M}{a (1 - e^2)} each orbit. This example integrates an eccentric Schwarzschild geodesic exactly (no approximation -- solving the full timelike geodesic equation of motion) and measures its perihelion precession directly from the trajectory, comparing against this leading-order post-Newtonian formula. To make the (otherwise tiny) effect visible, the orbit used here is deliberately far more relativistic than Mercury's, starting much closer to the black hole. .. GENERATED FROM PYTHON SOURCE LINES 26-32 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.relativity.chapters.schwarzschild import SchwarzschildBlackHole .. GENERATED FROM PYTHON SOURCE LINES 33-35 Integrate an eccentric orbit and watch it precess ------------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 35-51 .. code-block:: Python bh = SchwarzschildBlackHole(M=1.0) y0 = bh.eccentric_orbit_initial_state(r0=20.0, eccentricity_boost=0.15) traj = bh.integrate_geodesic(y0, dtau=0.02, n_steps=200000) fig, ax = plt.subplots(figsize=(6, 6)) x = traj["r"] * np.cos(traj["phi"]) y = traj["r"] * np.sin(traj["phi"]) ax.plot(x, y, linewidth=0.7, color="steelblue") circle = plt.Circle((0, 0), bh.horizon_radius, color="black") ax.add_patch(circle) ax.set_aspect("equal") ax.set_xlabel("x [M]") ax.set_ylabel("y [M]") ax.set_title("A precessing (rosette) Schwarzschild orbit") plt.tight_layout() .. image-sg:: /api/gallery/relativity/schwarzschild/images/sphx_glr_plot_orbit_precession_001.png :alt: A precessing (rosette) Schwarzschild orbit :srcset: /api/gallery/relativity/schwarzschild/images/sphx_glr_plot_orbit_precession_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 52-54 Measure the precession per orbit and compare to the weak-field formula ---------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 54-66 .. code-block:: Python precession = bh.perihelion_precession(traj) if len(precession) == 0: print("No full periapsis passage completed -- increase n_steps to see precession.") else: e_est = (traj["r"].max() - traj["r"].min()) / (traj["r"].max() + traj["r"].min()) a_semi = (traj["r"].max() + traj["r"].min()) / 2.0 weak_field = bh.weak_field_precession_per_orbit(a_semi, e_est) print(f"Measured precession per orbit: {np.degrees(np.mean(precession)):.3f} deg") print(f"Leading-order PN estimate: {np.degrees(weak_field):.3f} deg") print("(Mercury's actual GR excess is only ~43 arcsec/century -- this orbit is") print(" deliberately much closer to a strong-field black hole to make the effect visible.)") .. rst-class:: sphx-glr-script-out .. code-block:: none Measured precession per orbit: 139.731 deg Leading-order PN estimate: 86.128 deg (Mercury's actual GR excess is only ~43 arcsec/century -- this orbit is deliberately much closer to a strong-field black hole to make the effect visible.) .. GENERATED FROM PYTHON SOURCE LINES 67-69 Energy and angular momentum conservation (integration accuracy check) ---------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 69-81 .. code-block:: Python fig, axes = plt.subplots(1, 2, figsize=(10, 3.5)) axes[0].plot(traj["t"], (traj["energy"] - traj["energy"][0]) / traj["energy"][0]) axes[0].set_xlabel("coordinate time t") axes[0].set_ylabel("fractional energy drift") axes[1].plot( traj["t"], (traj["angular_momentum"] - traj["angular_momentum"][0]) / traj["angular_momentum"][0], ) axes[1].set_xlabel("coordinate time t") axes[1].set_ylabel("fractional angular momentum drift") plt.tight_layout() plt.show() .. image-sg:: /api/gallery/relativity/schwarzschild/images/sphx_glr_plot_orbit_precession_002.png :alt: plot orbit precession :srcset: /api/gallery/relativity/schwarzschild/images/sphx_glr_plot_orbit_precession_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.162 seconds) .. _sphx_glr_download_api_gallery_relativity_schwarzschild_plot_orbit_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_orbit_precession.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_orbit_precession.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_orbit_precession.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_