Relativistic orbit equation

Einstein's orbit equation u'' + u = alpha + epsilon u² of a planet, whose perihelion precesses.

• Name: relativistic_orbit
• Field: physics/celestial_mechanics
• Formalism: ordinary differential equations
• References: Einstein 1915; Dilão 2023

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('relativistic_orbit')
entry = describe('relativistic_orbit')
trajectory = create_simulation('relativistic_orbit').prepare().run().trajectory_

Factory

@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'second order', 'gravity', 'general relativity', 'precession'),
    references=('Einstein 1915', 'Dilão 2023'),
    defaults={'until': 60.0, 'dt': 0.02},
)
def relativistic_orbit(u: float = 1.5, du: float = 0.0, *, alpha: float = 1.0, epsilon: float = 0.05) -> System:
    """Return Einstein's orbit equation `u'' + u = alpha + epsilon u²` of a planet, whose perihelion precesses.

    The inverse radius `u = 1 / r` of the planet is a function of its polar angle, which the simulation time stands
    for, so the orbit in the plane is `(cos t / u, sin t / u)`. For `epsilon = 0` the orbit is a closed Newtonian
    ellipse of eccentricity `|u(0) / alpha - 1|` for `u'(0) = 0`. The term `epsilon u²` is the first relativistic
    correction, and with it the orbit is a rosette whose perihelion advances by about `2π epsilon alpha` every
    revolution. The equation has two fixed points, the roots `(1 ∓ sqrt(1 - 4 alpha epsilon)) / (2 epsilon)` of
    `epsilon u² - u + alpha`. The smaller one, close to `alpha`, is a centre, the stable circular orbit, and the
    larger one is a saddle, the unstable circular orbit close to the black hole, beyond which a planet falls in.

    Args:
        u:
            The initial inverse radius.

        du:
            The initial derivative of the inverse radius along the angle.

        alpha:
            The Newtonian inverse semi-latus rectum.

        epsilon:
            The strength of the relativistic correction.

    Returns:
        system:
            The system, with the inverse radius `planet.u` of order 2 along the angle and its derivative
            `planet.u_d1`.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import relativistic_orbit
        >>> relativistic_orbit(epsilon=0.0).elements[0].epsilon
        0.0
    """
    return System(elements=[RelativisticOrbitElement('planet', u=u, u_d1=du, alpha=alpha, epsilon=epsilon)])

Components

Elements:

• RelativisticOrbitElement, from Element

class RelativisticOrbitElement(Element):
    """Inverse radius of a planet along its orbit angle under the relativistic correction of Newtonian gravity.

    The inverse radius `u = 1 / r` follows Einstein's orbit equation `u'' + u = alpha + epsilon u²`, whose independent
    variable is the polar angle of the planet, which the simulation time stands for. For `epsilon = 0` it is the
    Newtonian conic `u = alpha (1 + e cos θ)`, and the quadratic term makes the perihelion advance every revolution.

    Args:
        label:
            The label of the planet.

        visual:
            How the planet is drawn, or `None` for `default_visual_`.

        u:
            The initial inverse radius.

        u_d1:
            The initial derivative of the inverse radius along the angle.

        alpha:
            The Newtonian inverse semi-latus rectum `G M / h²`, with `h` the angular momentum per unit mass.

        epsilon:
            The strength of the relativistic correction, `3 G M / c²`.

    Attributes:
        default_visual_:
            The visual of the planet when it sets none.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import RelativisticOrbitElement
        >>> RelativisticOrbitElement('planet').epsilon
        0.05
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        u: Annotated[float, StateSpec(order=2)] = 1.5,
        u_d1: float = 0.0,
        alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    ) -> None:
        super().__init__(label, visual=visual)
        self.u = u
        self.u_d1 = u_d1
        self.alpha = alpha
        self.epsilon = epsilon

    def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
        """Return Einstein's orbit equation.

        Args:
            namespace:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The rate of the second derivative of `u` along the angle.
        """
        s = namespace
        return [RateLaw(s.u, s.alpha + s.epsilon * s.u**2 - s.u)]

Equations

planet.u: diff(u_planet(t), t, t) == epsilon_planet*u_planet(t)^2 + alpha_planet - u_planet(t)