Newton's Formulation

In 1687, Isaac Newton published his magnum opus — Philosophiæ Naturalis Principia Mathematica (Mathematical Principles of Natural Philosophy). In this work, he formulated the law of universal gravitation and the laws of motion, laying the foundations of classical mechanics. It is here that the three-body problem was first posed: Newton noticed that the gravitational interaction of three celestial bodies defies precise analytical description.

Law of Universal Gravitation

At the core of Newton's formulation lies the law of universal gravitation. Two bodies with masses m₁ and m₂, separated by distance r, attract each other with a force:

F=Gm1m2r2F = G \frac{m_1 m_2}{r^2}

where G is the gravitational constant. The force acts along the line connecting the centers of mass of the two bodies.

Equations of Motion

For a system of N bodies, Newton's second law yields a system of differential equations:

mid2ridt2=∑j≠iGmimj∣rj−ri∣3(rj−ri)m_i \frac{d^2 \mathbf{r}_i}{dt^2} = \sum_{j \neq i} G \frac{m_i m_j}{|\mathbf{r}_j - \mathbf{r}_i|^3} (\mathbf{r}_j - \mathbf{r}_i)

For three bodies (i = 1,2,3) we obtain three second-order vector equations, equivalent to eighteen first-order scalar differential equations. This is the mathematical formulation of the three-body problem.

Why Is It Difficult?

Unlike the two-body problem, which has an elegant analytic solution (Kepler's laws), the three-body problem has no general closed-form solution. Newton understood this difficulty but lacked the analytical tools to overcome it. It would take over two centuries before Henri Poincaré proved that the three-body problem is generally non-integrable and exhibits chaotic behavior. This simulation uses the Euler method — the simplest numerical integration scheme, closest to Newton's original concept of motion under forces: acceleration is computed from force, then velocity is updated, followed by position.

Numerical Example

For the two-body problem, Newton derived an exact analytical solution — Kepler's laws. Orbits are conic sections, and positions can be computed exactly at any time. For the three-body problem, no analytical solution exists, so the simulation uses 2 bodies with calculations following Kepler's laws. Let's examine a binary system of two equal-mass bodies on circular orbits: body 1 at (40, 0), m = 200; body 2 at (−40, 0), m = 200. G = 80.

Gravitational force between the bodies:

F=Gm1m2r2=80⋅200⋅200802=500F = G \frac{m_1 m_2}{r^2} = 80 \cdot \frac{200 \cdot 200}{80^2} = 500

Kepler's third law gives the orbital period (a = 80 is the semi-major axis of the relative orbit):

T2=4π2a3G(m1+m2)=4π2⋅80380⋅400=64π2T=8π≈25.13\begin{aligned} T^2 &= \frac{4\pi^2 a^3}{G(m_1 + m_2)} = \frac{4\pi^2 \cdot 80^3}{80 \cdot 400} = 64\pi^2 \\ T &= 8\pi \approx 25.13 \end{aligned}

Orbital velocity of each body (a₁ = 40 is the distance to the center of mass):

v=2πa1T=2π⋅408π=10v = \frac{2\pi a_1}{T} = \frac{2\pi \cdot 40}{8\pi} = 10

Simulation