Three-Body Problem: Equations and Periodic Orbits

# Mathematical Modeling of the Three-Body Problem: From Equations to Periodic Orbits

The three-body problem involves calculating the trajectories of three massive objects under Newton's law of gravitational attraction. Initial conditions specify the positions and velocities of the bodies in 3D space, leading to a system of 18 nonlinear second-order differential equations.

In the full problem, the coordinates and velocities of each body are linked to accelerations that depend on their mutual distances. Normalized masses (ratio to the maximum) simplify the calculations. Here is the system in the form of first derivatives:

X1'  = VX1
Y1'  = VY1
Z1'  = VZ1
X2'  = VX2
Y2'  = VY2
Z2'  = VZ2
X3'  = VX3
Y3'  = VY3
Z3'  = VZ3
VX1' = -m2*(X1-X2)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5-m3*(X1-X3)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5
VY1' = -m2*(Y1-Y2)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5-m3*(Y1-Y3)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5
VZ1' = -m2*(Z1-Z2)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5-m3*(Z1-Z3)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5
VX2' = -m3*(X2-X3)/((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2)^1.5-m1*(X2-X1)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5
VY2' = -m3*(Y2-Y3)/((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2)^1.5-m1*(Y2-Y1)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5
VZ2' = -m3*(Z2-Z3)/((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2)^1.5-m1*(Z2-Z1)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5
VX3' = -m1*(X3-X1)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5-m2*(X3-X2)/((X3-X2)^2+(Y3-Y2)^2+(Z3-Z2)^2)^1.5
VY3' = -m1*(Y3-Y1)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5-m2*(Y3-Y2)/((X3-X2)^2+(Y3-Y2)^2+(Z3-Z2)^2)^1.5
VZ3' = -m1*(Z3-Z1)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5-m2*(Z3-Z2)/((X3-X2)^2+(Y3-Y2)^2+(Z3-Z2)^2)^1.5

Simplified Variants: Planar and Restricted Problems

The planar problem (P3BP) is restricted to 2D space, reducing the equations to 12. Here, the bodies remain in the plane defined by the initial conditions.

For the planar problem:

X1'  = VX1
X2'  = VX2
X3'  = VX3
Y1'  = VY1
Y2'  = VY2
Y3'  = VY3
VX1' = m2*(X2-X1)/((X1-X2)^2+(Y1-Y2)^2)^1.5+m3*(X3-X1)/((X1-X3)^2+(Y1-Y3)^2)^1.5
VX2' = m3*(X3-X2)/((X2-X3)^2+(Y2-Y3)^2)^1.5+m1*(X1-X2)/((X2-X1)^2+(Y2-Y1)^2)^1.5
VX3' = m1*(X1-X3)/((X3-X1)^2+(Y3-Y1)^2)^1.5+m2*(X2-X3)/((X3-X2)^2+(Y3-Y2)^2)^1.5
VY1' = m2*(Y2-Y1)/((X1-X2)^2+(Y1-Y2)^2)^1.5+m3*(Y3-Y1)/((X1-X3)^2+(Y1-Y3)^2)^1.5
VY2' = m3*(Y3-Y2)/((X2-X3)^2+(Y2-Y3)^2)^1.5+m1*(Y1-Y2)/((X2-X1)^2+(Y2-Y1)^2)^1.5
VY3' = m1*(Y1-Y3)/((X3-X1)^2+(Y3-Y1)^2)^1.5+m3*(Y2-Y3)/((X3-X2)^2+(Y3-Y2)^2)^1.5

The restricted problem (PR3BP) assumes a small mass for the third body. Equations for its motion:

X' = VX
Y' = VY
VX' = X+2*Y4-m1*(X+m2)/((X+m2)^2+Y^2)^1.5-m2*(X-m1)/((X-m1)^2+Y^2)^1.5
VY' = Y-2*Y3-m1*Y/((X+m2)^2+Y^2)^1.5-m2*Y/((X-m1)^2+Y^2)^1.5

The circular restricted problem (CR3BP) models a spacecraft in a system of two bodies with the smaller one in a circular orbit. The parameter μ = m2 / (m1 + m2) is used:

X' = VX
Y' = VY
z' = VZ
VX' = 2*VY+X-(1-mu) * (X+mu)/((X+mu)^2+Y^2+Z^2)^1.5-mu * (X-1+mu)/((X-1+mu)^2+Y^2+Z^2)^1.5
VY' = -2*VX+Y-(1-mu) * Y/((X+mu)^2+Y^2+Z^2)^1.5-mu * Y/((X-1+mu)^2+Y^2+Z^2)^1.5
VZ' = -(1-mu) * Z/((X+mu)^2+Y^2+Z^2)^1.5-mu * Z/((X-1+mu)^2+Y^2+Z^2)^1.5

These systems are solved using numerical methods in Mathematica, Python (SciPy), or MATLAB.

Analytical Solutions: Euler and Lagrange

The only known analytical solutions are Euler's collinear solutions (1767) and Lagrange's triangular solutions (1772). In the Lagrange case, the bodies start at the vertices of an equilateral triangle with matching velocities. The trajectories are ellipses rotated by 120°; the bodies maintain the triangular configuration.

These solutions are unstable: small deviations in initial conditions lead to chaos. Practical implementations require precision to 10^{-10}.

Periodic Orbits in the Planar Problem

Henri Poincaré proved the existence of infinitely many periodic solutions for equal masses. Computer searches have revealed families with symmetry.

• Broucke-Enò-Hadjidimitriou family (1975): 8- and 16-symbol orbits. Example R4 (period 5.4):

```

T0 = 0, T1 = 5.4

X1(T0) = 0.8733047091, X2(T0) = -0.6254030288, X3(T0) = -0.2479016803

Y1(T0) = 0, Y2(T0) = 0, Y3(T0) = 0

VX1(T0) = 0, VX2(T0) = 0, VX3(T0) = 0

VY1(T0) = 1.010776444, VY2(T0) = -1.683353346, VY3(T0) = 0.6725769022

```

• Shuvakov-Dmitrashinovich family (2013): Dragonfly (period 22):

```

T0 = 0, T1 = 22

X1(T0) = -1, X2(T0) = 1, X3(T0) = 0

Y1(T0) = 0, Y2(T0) = 0, Y3(T0) = 0

VX1(T0) = 0.08058, VX2(T0) = 0.08058, VX3(T0) = -0.16116

VY1(T0) = 0.58884, VY2(T0) = 0.58884, VY3(T0) = -1.17768

```

Others: Butterfly, Moth, Goggles, Yin-Yang.

• Shin family (2016): Oval, cat's face (period 5.1).

To reproduce, use 8th-order Runge-Kutta with adaptive step size.

Key Points

• The full problem requires 18 equations; simplifications reduce it to 4–6 for practical simulations.
• Euler/Lagrange analytical solutions are unstable and only apply in ideal conditions.
• Periodic orbits are found using numerical methods; Broucke and Shuvakov families contain hundreds of variants.
• Initial conditions must be specified with high precision for periodicity.
• CR3BP is critical for NASA/ESA missions: maneuvers in Earth-Moon and Sun-Jupiter systems.

— Editorial Team