Skip to content

N-body and solar system

Bodies in the plane that attract each other pairwise by Newtonian gravity.

  • Name: n_body
  • Field: physics/celestial_mechanics
  • Formalism: ordinary differential equations
  • References: Argun et al. 2021; Dilão 2023; Chenciner and Montgomery 2000

Factory

@register_model(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'gravity', 'chaos'),
    references=('Argun et al. 2021', 'Dilão 2023', 'Chenciner and Montgomery 2000'),
    defaults={'until': 6.33, 'dt': 0.01},
)
def n_body(
    preset: str = 'figure_eight',
    *,
    masses: Sequence[float] | None = None,
    positions: Sequence[Sequence[float]] | None = None,
    velocities: Sequence[Sequence[float]] | None = None,
    G: float | None = None,
) -> System:
    """Return bodies in the plane that attract each other pairwise by Newtonian gravity.

    The preset `'figure_eight'` is the periodic orbit of three equal masses found by Chenciner and Montgomery, with
    period 6.3259 for `G = 1`. The preset `'inner_solar_system'` is the Sun with Mercury, Venus, the Earth, and Mars on
    circular orbits, in astronomical units and years, with `G = 4 π²`. Masses, positions, velocities, or `G` given
    explicitly replace those of the preset.

    Args:
        preset:
            `'figure_eight'` or `'inner_solar_system'`.

        masses:
            The masses of the bodies, or `None` for those of the preset.

        positions:
            The initial positions of the bodies, or `None` for those of the preset.

        velocities:
            The initial velocities of the bodies, or `None` for those of the preset.

        G:
            The gravitational constant, or `None` for that of the preset.

    Returns:
        system:
            The model, with the particles `body_<i>` in two dimensions.

    Raises:
        ValidationError:
            The preset is unknown.

    Examples:
        >>> from skcomplex.models.physics.celestial_mechanics import n_body
        >>> len(n_body('inner_solar_system').interactions)
        10
    """
    preset_G, preset_masses, preset_positions, preset_velocities = _preset(preset)
    masses = preset_masses if masses is None else masses
    positions = preset_positions if positions is None else positions
    velocities = preset_velocities if velocities is None else velocities
    G = preset_G if G is None else G
    bodies = [
        MassPointElement(f'body_{index}', x=list(position), x_d1=list(velocity), m=mass)
        for index, (mass, position, velocity) in enumerate(zip(masses, positions, velocities, strict=True))
    ]
    attractions = [
        GravityInteraction(
            f'gravity_{first}_{second}',
            first=f'body_{first}',
            second=f'body_{second}',
            G=G,
        )
        for first, second in itertools.combinations(range(len(bodies)), 2)
    ]
    return System(elements=bodies, interactions=attractions, space=ContinuousSpace(n_dim=2))

Components

class MassPointElement(PointElement):
    """Point mass, a point whose position `x` has order 2, with its velocity `x_d1` and a mass `m`.

    The forces on the point divide by the mass through their divisor `'m'`.

    Args:
        label:
            The label of the point.

        x:
            The initial position. A number applies to every component.

        x_d1:
            The initial velocity. A number applies to every component.

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

        m:
            The mass.

    Attributes:
        default_visual_:
            Points at the positions.

    Examples:
        >>> from skcomplex.models.physics.common.particles import MassPointElement
        >>> MassPointElement('ball', x=[0.0, 1.0], m=2.0).m
        2.0
    """

    def __init__(
        self,
        label: str,
        x: npt.ArrayLike = 0.0,
        x_d1: npt.ArrayLike = 0.0,
        visual: VisualSpec | None = None,
        *,
        m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    ) -> None:
        super().__init__(label, x=x, x_d1=x_d1, visual=visual)
        self.m = m

class GravityInteraction(PotentialEnergyInteraction):
    """Newtonian gravity between two particles.

    The energy is `-G * m_first * m_second / |x_first - x_second|`, and its gradient moves both particles.

    Args:
        label:
            The label of the interaction.

        first:
            The path of the first particle.

        second:
            The path of the second particle.

        G:
            The gravitational constant.

        states:
            The name of the vector state that holds the positions of the particles.

        divisor:
            The name of the state or parameter that holds the masses of the particles.

    Examples:
        >>> from skcomplex.models.physics.common.particles import GravityInteraction
        >>> GravityInteraction('gravity', first='sun', second='earth', G=1.0).G
        1.0
    """

    def __init__(
        self,
        label: str,
        *,
        first: Annotated[str, RoleSpec()],
        second: Annotated[str, RoleSpec()],
        G: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        states: str = 'x',
        divisor: str = 'm',
    ) -> None:
        super().__init__(label, target=first, states=states, divisor=divisor)
        self.first = first
        self.second = second
        self.G = G

    def _get_potential(self, s: SimpleNamespace) -> Expression:
        """Return `-G * m_first * m_second / |x_first - x_second|`.

        Args:
            s:
                The namespace of the symbols the energy can use.

        Returns:
            energy:
                The potential energy as a symbolic expression.
        """
        gap = getattr(s.first, self._get_names()[0]) - getattr(s.second, self._get_names()[0])
        masses = getattr(s.first, str(self.divisor)) * getattr(s.second, str(self.divisor))
        return -s.G * masses / sqrt(gap.dot_product(gap))

Equations

body_0.x[0]: diff(x_body_0_0(t), t, t) == -G_gravity_0_1*m_body_1*(x_body_0_0(t) - x_body_1_0(t))/((x_body_0_0(t) - x_body_1_0(t))^2 + (x_body_0_1(t) - x_body_1_1(t))^2)^(3/2) - G_gravity_0_2*m_body_2*(x_body_0_0(t) - x_body_2_0(t))/((x_body_0_0(t) - x_body_2_0(t))^2 + (x_body_0_1(t) - x_body_2_1(t))^2)^(3/2)
body_0.x[1]: diff(x_body_0_1(t), t, t) == -G_gravity_0_1*m_body_1*(x_body_0_1(t) - x_body_1_1(t))/((x_body_0_0(t) - x_body_1_0(t))^2 + (x_body_0_1(t) - x_body_1_1(t))^2)^(3/2) - G_gravity_0_2*m_body_2*(x_body_0_1(t) - x_body_2_1(t))/((x_body_0_0(t) - x_body_2_0(t))^2 + (x_body_0_1(t) - x_body_2_1(t))^2)^(3/2)
body_1.x[0]: diff(x_body_1_0(t), t, t) == G_gravity_0_1*m_body_0*(x_body_0_0(t) - x_body_1_0(t))/((x_body_0_0(t) - x_body_1_0(t))^2 + (x_body_0_1(t) - x_body_1_1(t))^2)^(3/2) - G_gravity_1_2*m_body_2*(x_body_1_0(t) - x_body_2_0(t))/((x_body_1_0(t) - x_body_2_0(t))^2 + (x_body_1_1(t) - x_body_2_1(t))^2)^(3/2)
body_1.x[1]: diff(x_body_1_1(t), t, t) == G_gravity_0_1*m_body_0*(x_body_0_1(t) - x_body_1_1(t))/((x_body_0_0(t) - x_body_1_0(t))^2 + (x_body_0_1(t) - x_body_1_1(t))^2)^(3/2) - G_gravity_1_2*m_body_2*(x_body_1_1(t) - x_body_2_1(t))/((x_body_1_0(t) - x_body_2_0(t))^2 + (x_body_1_1(t) - x_body_2_1(t))^2)^(3/2)
body_2.x[0]: diff(x_body_2_0(t), t, t) == G_gravity_0_2*m_body_0*(x_body_0_0(t) - x_body_2_0(t))/((x_body_0_0(t) - x_body_2_0(t))^2 + (x_body_0_1(t) - x_body_2_1(t))^2)^(3/2) + G_gravity_1_2*m_body_1*(x_body_1_0(t) - x_body_2_0(t))/((x_body_1_0(t) - x_body_2_0(t))^2 + (x_body_1_1(t) - x_body_2_1(t))^2)^(3/2)
body_2.x[1]: diff(x_body_2_1(t), t, t) == G_gravity_0_2*m_body_0*(x_body_0_1(t) - x_body_2_1(t))/((x_body_0_0(t) - x_body_2_0(t))^2 + (x_body_0_1(t) - x_body_2_1(t))^2)^(3/2) + G_gravity_1_2*m_body_1*(x_body_1_1(t) - x_body_2_1(t))/((x_body_1_0(t) - x_body_2_0(t))^2 + (x_body_1_1(t) - x_body_2_1(t))^2)^(3/2)