Simulating Three Body Orbits

Simulating Classical Mechanics

Simulating classical physics in a computer is straight forward. Every few years I write a new simulation script to investigate some interesting phenomena, typically planet orbits. The pseudocode is simple:

For a tiny change in time
   For each object
      For each other object
         sum all the forces
      compute the acceleration based on the forces
      compute the velocity change based on the acceleration
      compute the position change based on the velocity

The challenge is the error introduced by applying one force/acceleration across an entire period of time in each step. In reality this force changes continuously as the object moves, but the simulation only calculates it once for the entire "tiny time" of the step. You can use mathematical tricks to minimize this error, but it is an approximation no matter what. What's worse is that this error tends to add up across the many steps, so a small error quickly throws off the entire solution.

That being said, it is still possible to run simulations like this by making the step in time as small as possible. This minimizes how much the force should have changed during that step, and hence minimizes the error. This accuracy extends the amount of time the simulation remains accurate before the error compounds too much. Since planets often move in repeating patterns, like the Earth orbiting the sun each year, if you can calculate one period of that orbit accurately it's often enough for simple problems.

Three Bodies

When talking about planets moving due to gravity, it's possible to mathematically solve the exact orbits for two bodies only affected by each other. Given the starting position, we know the exact equations describing how they will move. However, no general solution exists when you introduce a third body. This increase in complexity makes the three body problem really interesting. Moreover, it's interesting that there are some known starting positions that CAN be described by mathematical equations, just not most of them.

I recently stumbled across a paper (https://arxiv.org/pdf/1805.07980.pdf) from 2018 that found some new starting positions that led to stable periodic orbits for three bodies in 2 dimensions. This got me itching again to write a quick simulator and see if I could reproduce the orbits they described. Moreover I wanted to see how the step size actually ends up affecting the decay of these solution.

One of my favorite orbits they found was surprisingly simple. Some of the others were quite complex, but this one didn't have too many spirals or back-and-forths. Here is a diagram of the paths the bodies follow for these starting conditions:

For the most part, the solutions in the paper all look a bit funny because the planet goes back and forth on a line rather than in a circle. The planets have zero velocity at the ends of the lines, accelerate towards each other, and then slow back down as they move apart. They all do this in a synchronized dance that puts them in the exact same position over and over.

One cool part about this type of orbit is that if you start the three planets on their "zero velocity" positions, they will orbit correctly. This is distinct from a traditional elliptical orbit, which requires an initial velocity in addition to a relative position. Even cooler, the orbits in the paper all kept the same starting position for two of the planets, so the only inputs into the entire system were the mass of each planet and the coordinates of the third planet.

The Simulation

I chose the F1 orbit I mentioned above from the paper because it looked the simplest. This makes it easy to validate visually. I ran the simulation using pypy - a python interpreter that did Just-In-Time compilation to accelerate the code. This makes the code about an order of magnitude faster on my machine. Python generally isn't the right language for this type of computation, but I was able to get it to run in under a minute for a full period to an accuracy I was happy with.

The Results

The simulation was able to reproduce the orbits for at least one period with good overlap as the bodies returned to their initial positions. My simulation produced this plot:

However, how well my simulation reproduced the orbits depended heavily on how small I made the time steps. Here are some examples with larger time steps:

If you take the total positional error at the end of the first period, you can plot out the impact the step size has on the final error.

The only factor that really keeps us from reducing the time step further is the computation time, so we can also plot how the error varies as we spend more time on the computation.

Conclusion

This orbit does in fact look stable! As we improve our simulation, the error from starting position to finish position across one period trends towards zero. This implies that the next period will behave the same way as long as our time step is small enough. I found it very cool that this paper showed some simple starting positions to make testing a simulation like this so easy!