Solving the Three Body Problem With Python

Rhett Allain
9 min readMar 10
Image: Rhett Allain. Stable 3 Body Motion

Yes, it’s time for the famous Three Body Problem. How do you model the motion of three objects interacting gravitational (those are the 3 bodies)? Well, let’s jump back for a second. How do you model just TWO objects with a gravitational interaction?

The Two Body Problem

Here are two planets (or stars — it doesn’t matter, you pick).

With the vector from planet A to B (r_AB), I can calculate the gravitational force on B as:

Just to be clear, I’m using notation that F_AB is the force ON B due to A. With that, the force on A (due to B) is just the negative of this vector. Now, let’s say that we break this planet motion into “short” time intervals (Δt) — and by “short” I mean that it’s short enough that we can approximate the forces as being constant during this time. Depending on the initial conditions this time interval could be as short as a minute or maybe as long as an hour (or more even).

Suppose planet B starts with a momentum p_B1. Since we know the force acting on it (and assuming this force is constant) we can use the momentum principle to find the momentum at the end of the time interval (p_B2).

Yes, this is a numerical calculation using the Euler method. Now, I can assume that the momentum is constant during this time interval to find the new position of planet B (using the average velocity).

Now I just need to do the same thing for planet A and that will give me the position at the end of that one time interval. If I keep repeating this process — I get the motion of the two planets. Yes, I’m going to do this in python when we get to three planets.

Actually, we don’t have to use a numerical calculation to model this motion. You can exactly solve the 2…

Rhett Allain

Physics faculty, science blogger of all things geek. Technical Consultant for CBS MacGyver and MythBusters. WIRED blogger.