The Three-Body Problem
Chaos and stability seem like concepts that can’t coexist within the same system. That’s only partly true. A system can be entirely deterministic and still be unpredictable if it is sensitive enough to its initial conditions. Small perturbations, given enough time, can send such a system into results no one could have foreseen. It’s a strange idea: a system governed by fixed physical laws, yet resistant to prediction.
However, chaos doesn’t necessarily mean instability. Given specific initial conditions, a system that is chaotic can remain stable over a relatively long period. This piece explores the three-body problem using fundamental physical laws, and simulates its behaviour through numerical integration.
In my previous articles, I investigated coupled systems governed by a restoring force, like the coupled harmonic oscillator, where a spring pulls each mass back toward equilibrium. Here, the coupling is different: a central force, like gravity, where no such equilibrium exists and nothing pulls the system back to where it started.
The classical three-body problem
Newton’s law of universal gravitation and his laws of motion provide a comprehensive framework for approaching the physics of any classical coupled system. These laws yield equations of motion that fully describe how the system evolves. But having the equations isn’t the same as solving them. Further mathematical methods are needed to extract information about the state as closed-form analytic solutions. For relatively simpler systems such as the two-body problem or the coupled harmonic oscillator, there exist analytic solutions for which the states of the systems are defined. The three-body problem, however, has no general solution in the familiar closed-form sense. Or does it?
The three-body problem is often presented as an unsolvable mystery that has eluded scientists since Newton formulated his gravitational law for the behaviour of any two masses. For two bodies, the resulting system of differential equations can be solved exactly. The motion is always confined to a plane, since angular momentum is conserved, so the state is described by eight phase-space dimensions (the position and velocity components of the two bodies). The two-body system reduces to an equivalent one-body problem, which is solvable.
However, when a third body is added, each body tugs on the other two in proportion to their masses. Nothing forces this system to remain in a plane. To keep things visual, I will explore this problem in a plane. This is not an approximation, since planar initial position and velocity components produce motion that is planar forever. In the planar three-body problem, we obtain differential equations with twelve phase-space dimensions. This system of equations is much harder to reduce to a solvable form. Why?
The differential equations describing motion in a three-body system are:
Here, rᵢ is the position vector of body i, mᵢ is its mass, and G is the gravitational constant. The double dot denotes the second derivative with respect to time (acceleration), and the index runs over the three bodies (i = 1, 2, 3).
Each of these is a vector equation. Written in components, the three equations become nine scalar second-order equations in three dimensions, or six in the plane. Since each second-order equation needs a position and a velocity to specify its state, this gives the twelve phase-space dimensions of the planar problem.
The right-hand side of the equations above shows that the acceleration of any body depends on the positions of the other two, which are also changing with time. It is therefore impossible to solve any one equation in isolation, so the equations must be solved simultaneously. This is the coupling: every body’s motion is bound up with the other two.
Coupling alone isn’t the obstacle, since coupled oscillators and the two-body system both have exact solutions. The real limitation is a shortage of conserved quantities. The planar problem has six: energy, angular momentum, and the center-of-mass motion (two components of linear momentum and two of position). For two bodies, these six, roughly speaking, reduce the eight-dimensional phase space to two dimensions, which can be integrated analytically. For three bodies, the same six leave six of the twelve dimensions unresolved. Adding a body adds four dimensions but no new conservation laws. At the end of the nineteenth century, Poincaré showed that no further independent conserved quantities of the usual kind exist, so the system cannot be reduced by the classical methods.
Get ayilara’s stories in your inbox
Join Medium for free to get updates from this writer.
If the problem cannot be reduced analytically, the only path left is to compute the motion step by step.
Numerical solutions of the three-body problem
Numerical integration doesn’t solve the three-body problem in general. It computes one trajectory at a time, for a chosen set of initial conditions, by advancing the equations of motion in small time steps. I use the fourth-order Runge-Kutta method (RK4), which estimates the acceleration at four points within each step and combines them into a weighted average. The error shrinks rapidly as the step size h is reduced (the accumulated error scales as h⁴), so RK4 is accurate without being complicated.
I set the gravitational constant G and all three masses to 1, use a step size of h = 0.01, and integrate for 30 time units. The resulting trajectories are shown in the animation below, starting from the initial conditions in the table. The full code is in my GitHub repository. You can interact with it and test out other initial condition values.
It is not perfect, though. RK4 doesn’t preserve the geometric structure of gravitational motion, so the total energy slowly drifts even though it should stay constant. That drift is useful. Plotting the relative energy drift against time gives a direct check on the simulation: if it stays small to within a small tolerance, the numerics are obeying the physics. Close encounters are where it fails first, because the force grows as 1/r² when two bodies pass near each other, and a fixed step size can no longer resolve the motion.
Figure 2 shows the relative energy drift on a logarithmic scale. It starts at machine precision and grows to about 10⁻⁷ by t ≈ 2, then stays between 10⁻⁷ and 10⁻⁶ until t ≈ 10. At t ≈ 10, the drift jumps by six to seven orders of magnitude, to about 0.5. This is the signature of a close encounter: two bodies pass near each other, the 1/r² force spikes, and a step of h = 0.01 is too coarse to follow it. After that point, the accumulated numerical error becomes large enough that the computed trajectory can no longer be considered a reliable approximation to the physical solution.
Numerical methods let us compute the motion, but they don’t tell us whether a formula exists.
The search for analytic solutions
Numerical methods have grown far more sophisticated, but they still compute one trajectory at a time. Is there a general formula? Poincaré’s result seemed to close the question, though what he showed was only the assertion that no further conserved quantities of the classical kind exist to reduce the problem.
In 1912, the Finnish mathematician Karl Sundman proved that the general three-body problem does have an analytic solution, in the form of a convergent infinite series known as the Puiseux series. It holds for all time and for every initial condition with non-zero total angular momentum. The only catch with this solution is convergence: the series needs an astronomically large number of terms (on the order of 10⁸⁰⁰⁰⁰⁰⁰) to reach useful precision. It is a solution in principle and useless in practice, because the amount of terms needed cannot be handled by any physical machine.
What I think
What the three-body problem ultimately demonstrates is that determinism and predictability are not the same thing. The equations completely determine the future once the initial state is specified, yet the system can remain extraordinarily sensitive to tiny changes in that state. Numerical error adds another layer of uncertainty: even a tiny error introduced by the integration method can behave like a perturbation to the initial conditions and grow over time.
That is what makes the problem so interesting. The equations themselves are deterministic. The difficulty is not that nature has stopped following a rule; it is that the consequences of those rules can become extraordinarily difficult to predict.
