Three-body problem

step 0
Loading the engineShared by all simulations, so the next one opens instantly.
The physics

Two bodies under gravity move on ellipses forever. Add a third and the motion is in general chaotic: bodies pair up, swap partners, and often one is eventually thrown out of the system. Henri Poincaré’s prize memoir on the problem, published in 1890, contained the first description of what is now called chaos.

What you see

Colored discs are the bodies, and the lines behind them are their trails. By default the leapfrog integrator (velocity Verlet) advances them; it is symplectic, so the energy error stays bounded instead of drifting. Total energy and angular momentum are plotted alongside.

Sun-Earth-Moon: Earth starts on a circular orbit 200 units from the Sun and the Moon on a circular orbit 30 units from Earth. The proportions are scaled for the screen: the Sun is about 330 times heavier than Earth (the real factor is about 333,000) and Earth 75 times heavier than the Moon. At this distance the Moon starts outside Earth’s Hill sphere (radius about 20 units), the region where Earth’s pull dominates the Sun’s tide, so it does not keep a regular monthly orbit.

Figure-8 orbit: three equal masses start from the positions of the Chenciner and Montgomery figure-eight, the middle body moving and the outer two moving at half its speed in the opposite direction, so the total momentum is zero. The bodies then chase each other along the figure eight.

Lagrange point L4: a small body starts 60° ahead of Earth on Earth’s orbit, moving at Earth’s orbital speed. L4 is stable when the larger mass outweighs the smaller by a factor of about 25 or more (here about 330), so the body librates around L4 instead of drifting off, like the Trojan asteroids that share Jupiter’s orbit around its L4 and L5 points.

Random triple: three equal masses at random positions 65 to 130 units from the center, with random velocities that add up to zero momentum. Every reset draws a new configuration.

Governing equation

mi𝐫¨i=G∑j≠imimj|𝐫j−𝐫i|2𝐫ˆij

Here 𝐫ˆij is the unit vector from body i to body j. A softening length ε is added inside the denominator to remove the 1/r2 singularity during close approaches.