by Paulo Marques, 2014/03/13

Updated 2016/02/15 for supporting Julia 0.4. Thanks mbeltagy for the contribution.

This notebook implements a numeric simulation of the three-body problem. This is done using the Julia Language and the Sundials library.

The problem corresponds to determine the trajectories of three bodies in space under the influence of their mutual gravity forces (e.g., the motion of the Sun, Moon and Earth, taken together). Isaac Newton was the first person to try to solve this problem. In 1887, Bruns and Poincaré showed that there is no general analytical closed solution for this problem. In fact, the motion of the three bodies is generally non-repeating! That's pretty amazing.

Newton taught us that two bodies with masses $m_0$ and $m_1$ attract themselves with a force proportional to their masses and inversionally proportional to the square of their distance $d$. This is called the Universal Law of Gravitation. Mathematically this is expressed by:

$$ F = G \frac{m_0 \, m_1}{d^2} $$where $G$ is the gravitational constant. But, if we want to consider the movement of three bodies in space, we'll need to formulate this in terms of vector calculus. Let's suppose that mass $m_0$ is at position $\vec{r_0}$ and $m_1$ at position $\vec{r_1}$. This is shown in the next diagram (which also shows a mass $m_2$ and its position -- for now, assume it doesn't exist):

In this case, the law of gravitation can be formulated as:

$$ \vec{F_{01}} = - G \frac{m_0 m_1}{{\vert \vec{r_0} - \vec{r_1} \vert}^2} \hat{r_{01}} $$where $\vec{F_{01}}$ is the force exerted on $m_1$ due to $m_0$ and $\hat{r_{01}}$ is the unit vector from $m_0$ to $m_1$:

$$ \hat{r_{01}} = \frac{\vec{r_1} - \vec{r_0}}{\vert \vec{r_1} - \vec{r_0} \vert} $$By simple manipulation we can write:

$$ \vec{F_{01}} = - G \, m_0 \, m_1 \frac{\vec{r_1} - \vec{r_0}}{{\vert \vec{r_1} - \vec{r_0} \vert}^3} $$Finally, there's an important thing that we also know: force equals mass times accelaration (Newton's second law) and that accelaration is the second derivative of position:

$$ \vec{F} = m \frac{d^2\vec{r}}{dt^2} $$Thus we can write:

$$ m_1 \frac{d^2\vec{r_1}}{dt^2} = - G \, m_0 \, m_1 \frac{\vec{r_1} - \vec{r_0}}{{\vert \vec{r_1} - \vec{r_0} \vert}^3} $$If we only had two bodies this would be great. But, with three bodies, we need to consider all the interactions and forces. Since we are calculating the accelaration of $m_1$ we just need to add the force corresponding to $\vec{F_{21}}$:

$$ m_1 \frac{d^2\vec{r_1}}{dt^2} = - G \, m_0 \, m_1 \frac{\vec{r_1} - \vec{r_0}}{{\vert \vec{r_1} - \vec{r_0} \vert}^3} \ - G \, m_2 \, m_1 \frac{\vec{r_1} - \vec{r_2}}{{\vert \vec{r_1} - \vec{r_2} \vert}^3} $$From the last equation above and by symmetry we get the follwing three equations:

$$ \frac{d^2\vec{r_0}}{dt^2} = G \, m_1 \frac{\vec{r_1} - \vec{r_0}}{{\vert \vec{r_1} - \vec{r_0} \vert}^3} + \ G \, m_2 \frac{\vec{r_2} - \vec{r_0}}{{\vert \vec{r_2} - \vec{r_0} \vert}^3} $$$$ \frac{d^2\vec{r_1}}{dt^2} = G \, m_0 \frac{\vec{r_0} - \vec{r_1}}{{\vert \vec{r_0} - \vec{r_1} \vert}^3} + \ G \, m_2 \frac{\vec{r_2} - \vec{r_1}}{{\vert \vec{r_2} - \vec{r_1} \vert}^3} $$$$ \frac{d^2\vec{r_2}}{dt^2} = G \, m_0 \frac{\vec{r_0} - \vec{r_2}}{{\vert \vec{r_0} - \vec{r_2} \vert}^3} + \ G \, m_1 \frac{\vec{r_1} - \vec{r_2}}{{\vert \vec{r_1} - \vec{r_2} \vert}^3} $$These can be integrated numerically. Let's do so.

Let's start by importing a ploting library and a differential equation solver.

In [1]:

```
using PyPlot
using Sundials
```

Now we need to define a function $\dot{g} = f(t, g)$ where $g$ is a function of time $t$ and $\dot{g}$ represents its derivative. Using the Sundials library this is done by defining a function `f(t, g, gdot)`

. The function should update the return parameter `gdot`

with the derivative of `g(t)`

at time `t`

. Note that this design is not particularly great -- it would be much better if we could just return `gdot`

. (Note that typically $\dot{g} = f(t, g)$ is written as $\dot{y} = f(t, y)$. We are not doing so because we'll be using $x$ and $y$ to denote coordinates.)

Looking at the equations above it's easy to see that both `y`

and `ydot`

are going to be 12-entry vectors. We are going to work in two dimentions: X and Y. Thus, each vector $\vec{r_k(t)}$ is going to correspond to two coordinates: $\vec{r_k(t)} = (x_k(t), y_k(t))$. Since we have 3 of these position vectors this corresponds to 6 entries. Since the differencial equations above are of second order, we will need the derivatives of each one of $\vec{r_k(t)}$, which corresponds to the bodies velocities $\vec{v_k(t)}$. Since we have 3 bodies, each with an X and Y velocity, we get 6 other entries.

So, let's write this function using the 3 differential equations above. For simplicity, we'll assume that $G=1$.

In [2]:

```
function f(t, g, gdot)
# Extract the position and velocity vectors from the g array
r0, v0 = g[1:2], g[3:4]
r1, v1 = g[5:6], g[7:8]
r2, v2 = g[9:10], g[11:12]
# The derivatives of the position are simply the velocities
dr0 = v0
dr1 = v1
dr2 = v2
# Now calculate the the derivatives of the velocities, which are the accelarations
# Start by calculating the distance vectors between the bodies (assumes m0, m1 and m2 are global variables)
# Slightly rewriten the expressions dv0, dv1 and dv2 comprared to the normal equations so we can reuse d0, d1 and d2
d0 = (r2 - r1) / ( norm(r2 - r1)^3.0 )
d1 = (r0 - r2) / ( norm(r0 - r2)^3.0 )
d2 = (r1 - r0) / ( norm(r1 - r0)^3.0 )
dv0 = m1*d2 - m2*d1
dv1 = m2*d0 - m0*d2
dv2 = m0*d1 - m1*d0
# Reconstruct the derivative vector
gdot[:] = [dr0; dv0; dr1; dv1; dr2; dv2]
end;
```

In [3]:

```
# Masses of the bodies
m0 = 5.0
m1 = 4.0
m2 = 3.0
# Initial positions and velocities of each body (x0, y0, vx0, vy0)
gi0 = [ 1.0; -1.0; 0.0; 0.0]
gi1 = [ 1.0; 3.0; 0.0; 0.0]
gi2 = [-2.0; -1.0; 0.0; 0.0]
# Simulation from time t=0 to t=30
tf = 30.0
dt = 500
t = collect(linspace(0.0, tf, round(Int,tf*dt)))
g0 = [gi0; gi1; gi2]
res = Sundials.cvode(f, g0, t, reltol=1e-10)
# Extract the position and velocity vectors from the solution
r0, v0, r1, v1, r2, v2 = res[:,1:2], res[:,3:4], res[:,5:6], res[:,7:8], res[:,9:10], res[:,11:12]
# Calculate the center of mass
cx = [(r0[i,1]*m0 + r1[i,1]*m1 + r2[i,1]*m2) / (m0 + m1 + m2) for i=1:length(t)]
cy = [(r0[i,2]*m0 + r1[i,2]*m1 + r2[i,2]*m2) / (m0 + m1 + m2) for i=1:length(t)]
# Write the results to file if we want to do something with them later on (e.g., create some videos)
# writecsv("results.csv", res);
```

`plot_trajectory`

and will pass as parameter the start time and stop time to plot.

In [4]:

```
function plot_trajectory(t1, t2)
t1i = round(Int,(length(t)-1) * t1/tf) + 1
t2i = round(Int,(length(t)-1) * t2/tf) + 1
# Plot the initial and final positions
# In these vectors, the first coordinate will be X and the second Y
X = 1
Y = 2
figure(figsize=(6,6))
plot(r0[t1i,X], r0[t1i,Y], "ro")
plot(r0[t2i,X], r0[t2i,Y], "rs")
plot(r1[t1i,X], r1[t1i,Y], "go")
plot(r1[t2i,X], r1[t2i,Y], "gs")
plot(r2[t1i,X], r2[t1i,Y], "bo")
plot(r2[t2i,X], r2[t2i,Y], "bs")
# Plot the trajectories
plot(r0[t1i:t2i,X], r0[t1i:t2i,Y], "r-")
plot(r1[t1i:t2i,X], r1[t1i:t2i,Y], "g-")
plot(r2[t1i:t2i,X], r2[t1i:t2i,Y], "b-")
# Plot cente of mass
plot(cx[t1i:t2i], cy[t1i:t2i], "kx")
# Setup the axis and titles
xmin = minimum([r0[t1i:t2i,X]; r1[t1i:t2i,X]; r2[t1i:t2i,X]]) * 1.10
xmax = maximum([r0[t1i:t2i,X]; r1[t1i:t2i,X]; r2[t1i:t2i,X]]) * 1.10
ymin = minimum([r0[t1i:t2i,Y]; r1[t1i:t2i,Y]; r2[t1i:t2i,Y]]) * 1.10
ymax = maximum([r0[t1i:t2i,Y]; r1[t1i:t2i,Y]; r2[t1i:t2i,Y]]) * 1.10
axis([xmin, xmax, ymin, ymax])
title(@sprintf "3-body simulation for t=[%.1f .. %.1f]" t1 t2)
end;
```

Let's plot the trajectory from 0 to 10 time units.

In [5]:

```
plot_trajectory(0, 10);
```

That's pretty crazy! The trajectory of each body is drawn in a different color. The start position is marked by a small circle and the end position by a square. The center of gravity is marked with a black X. As you can see the behavior is pretty insane. "Almost collisions" and unstable aperiodic trajectories.

And it goes on and on...

In [6]:

```
plot_trajectory(10, 20);
```

It almost looks like they are dancing...

In [7]:

```
plot_trajectory(20, 30);
```

I've also decided to make three small videos containing the simulation. [The videos weren't generated here because (i)Julia's support for doing this is not great...]

Look at all the close encounters!

In [3]:

```
# Function for creating an embedded video given a filename
function html_video(filename)
base64_video = base64encode(open(readbytes, filename))
"""<video controls src="data:video/x-m4v;base64,$base64_video">"""
end
videos = [html_video(v) for v=[
"imgs/3d_body_0_30_small.mp4";
"imgs/3d_body_0_10_small.mp4";
"imgs/3d_body_10_20_small.mp4";
"imgs/3d_body_20_30_small.mp4"]]
table = """
<table>
<tr> <td>$(videos[1])</td> <td>$(videos[2])</td> </tr>
<tr> <td>$(videos[3])</td> <td>$(videos[4])</td> </tr>
</table>"""
display("text/html", table)
```

Copyright (C) 2014 Paulo Marques (pjp.marques@gmail.com)

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice sha