The three-body problem¶
Let us realise the numerical solution of the problem of simulation of the dynamics of three bodies in gravitational interaction with each other.
Problem description¶
Three material points with coordinates $\bar x_1$, $\bar x_2$, $\bar x_3$ and velocities $\bar v_1$, $\bar v_2$, $\bar v_3$ are in a single space. We need to determine their trajectories in time.
As we know, unfortunately, there is no analytical solution to the problem of calculating the coordinates of all objects of the system at any time t. But we will solve this problem using numerical integration.
$$\begin{align} \mathbf{{\ddot{r}}_1} &= -G m_2 \frac{{\mathbf{r_1}} - {\mathbf{r_2}}} {{|{\mathbf{r_1}} - {\mathbf{r_2}}|}^3} - G m_3 \frac{{\mathbf{r_1}} - {\mathbf{r_3}}} {{|{\mathbf{r_1}} - {\mathbf{r_3}}|}^3} \\ \mathbf{{\ddot{r}}_2} &= -G m_3 \frac{{\mathbf{r_2}} - {\mathbf{r_3}}} {{|{\mathbf{r_2}} - {\mathbf{r_3}}|}^3} - G m_1 \frac{{\mathbf{r_2}} - {\mathbf{r_1}}} {{|{\mathbf{r_2}} - {\mathbf{r_1}}|}^3} \\ \mathbf{{\ddot{r}}_3} &= -G m_1 \frac{{\mathbf{r_3}} - {\mathbf{r_1}}} {{|{\mathbf{r_3}} - {\mathbf{r_1}}|}^3} - G m_2 \frac{{\mathbf{r_3}} - {\mathbf{r_1}}} {{|{\mathbf{r_3}} - {\mathbf{r_1}}|}^3} \\ \end{align}$$
$\mathbf{r_1} \dots \mathbf{r_3}$ here denote radius-vectors of each body, $\mathbf{\ddot r_1} \dots \mathbf{\ddot r_3}$ - their second derivatives (acceleration), $m_1 \dots m_3$ - mass of each body, gravitational constant $G$ can be taken equal to 1 without losses for our needs for visualisation of the system dynamics. The designation $|{\mathbf{r_1}} - {\mathbf{r_2}}|$ implies the norm between the radius-vectors, simply speaking - the distance between bodies 1 and 2.
So, the gravity of each of the two bodies affects the acceleration of the third. We have formulated equations for which there is no analytical solution, moreover, at $G=1$ our solution will be abstract and will not have the correct dimensions.
Solution of the system of differential equations¶
To solve this system of three second-order differential equations, we will also need to solve the first-order equations. Let's use the library DifferentialEquations:
Pkg.add(["DifferentialEquations"])
using DifferentialEquations
gr()
julia = Colors.JULIA_LOGO_COLORS # Пусть графики будут в цветах логотипа Julia
Solve at the current step and update the state of the system we will create the function thereebody.
function threebody( du, u, p, t )
# Положение, скорость
x1, y1, x2, y2, x3, y3,
vx1, vy1, vx2, vy2, vx3, vy3 = u
# Параметры системы (масса каждого тела)
m1, m2, m3 = p
# Приращение координат каждого тела
du[1] = vx1; du[2] = vy1; du[3] = vx2; du[4] = vy2; du[5] = vx3; du[6] = vy3;
# Расстояние от одного тела до другого
r12 = sqrt((x1 - x2) ^ 2 + (y1 - y2)^2 )
r23 = sqrt((x2 - x3) ^ 2 + (y2 - y3)^2 )
r13 = sqrt((x1 - x3) ^ 2 + (y1 - y3)^2 )
# Ускорение (при G=1)
du[7] = -m2 * ((x1-x2) / (r12^3)) - m3 * ((x1-x3) / (r13^3))
du[8] = -m2 * ((y1-y2) / (r12^3)) - m3 * ((y1-y3) / (r13^3))
du[9] = -m3 * ((x2-x3) / (r23^3)) - m1 * ((x2-x1) / (r12^3))
du[10] = -m3 * ((y2-y3) / (r23^3)) - m1 * ((y2-y1) / (r12^3))
du[11] = -m1 * ((x3-x1) / (r13^3)) - m2 * ((x3-x2) / (r23^3))
du[12] = -m1 * ((y3-y1) / (r13^3)) - m2 * ((y3-y2) / (r23^3))
end
You can see that vectorisation would shorten our code considerably, as would using the norm() function from the LinearAlgebra library to calculate distances. But in this tutorial example, we will favour the visual similarity of the code and formulas in the current implementation for the time being.
Let's set the initial parameters¶
We have a large vector of initial coordinates. Let's add some input fields to make it easy to change these parameters.
As values of initial velocity we will take the following numbers (in one or another direction):
$$\begin{align} vx_0 = -0.148 \cdot \cos( 60^\circ ) \\ vy_0 = -0.148 \cdot \sin( 60^\circ ) \end{align}$$
Note:
1) The following cell is not automatically executed, to update the values of the variables, press ⏵, or click in the free area of the mask and run the cell via keyboard shortcut Ctrl+Enter.
2) To hide the code and leave only the controls, double-click on the mask
# @markdown # Начальное положение $\bar x_1 \dots \bar x_3$
Положение_по_умолчанию = false #@param {type: "boolean"}
if Положение_по_умолчанию
x1_0 = -1.0
y1_0 = 0.0
x2_0 = 1.0
y2_0 = 0.0
x3_0 = 0.0
y3_0 = sqrt(3)
else
x1_0 = -1.0 # @param {type:"slider", min:-2, max:2, step:0.1}
y1_0 = 0.0 # @param {type:"slider", min:-2, max:2, step:0.1}
x2_0 = 1.0 # @param {type:"slider", min:-2, max:2, step:0.1}
y2_0 = 0.0 # @param {type:"slider", min:-2, max:2, step:0.1}
x3_0 = 0.0 # @param {type:"slider", min:-2, max:2, step:0.1}
y3_0 = 1.732 # @param {type:"slider", min:-2, max:2, step:0.1}
end
scatter( [x1_0, x2_0, x3_0],
[y1_0, y2_0, y3_0],
mc=[julia.red, julia.purple, julia.green],
ms=15, xlimits=(-2.5,2.5),
ylimits=(-2.5,2.5),
legend=:none, aspect_ratio=:equal )
# @markdown # Начальная скорость $\bar v_1 \dots \bar v_3$
# @markdown
Скорость_по_умолчанию = false #@param {type: "boolean"}
if Скорость_по_умолчанию
vx1_0, vy1_0, vx2_0, vy2_0, vx3_0, vy3_0 =
0.074, 0.128, -0.074, -0.128, 0.148, 0.0
else
vx1_0 = -0.074 # @param {type:"slider", min:-1.0, max:1.0, step:0.001}
vy1_0 = 0.128 # @param {type:"slider", min:-1.0, max:1.0, step:0.001}
vx2_0 = -0.074 # @param {type:"slider", min:-1.0, max:1.0, step:0.001}
vy2_0 = -0.128 # @param {type:"slider", min:-1.0, max:1.0, step:0.001}
vx3_0 = 0.148 # @param {type:"slider", min:-1.0, max:1.0, step:0.001}
vy3_0 = 0.0 # @param {type:"slider", min:-1.0, max:1.0, step:0.001}
end
scatter( [x1_0, x2_0, x3_0], [y1_0, y2_0, y3_0],
mc=[julia.red, julia.purple, julia.green], ms=15, xlimits=(-2.5,2.5), ylimits=(-2.5,2.5),
legend=:none, aspect_ratio=:equal )
plot!( [x1_0, x1_0+4*vx1_0], [y1_0, y1_0+4*vy1_0], c=julia.red, lw=2, arrow=true )
plot!( [x2_0, x2_0+4*vx2_0], [y2_0, y2_0+4*vy2_0], c=julia.purple, lw=2, arrow=true )
plot!( [x3_0, x3_0+4*vx3_0], [y3_0, y3_0+4*vy3_0], c=julia.green, lw=2, arrow=true )
# @markdown # Параметры $\bar m_1 \dots \bar m_3$
# @markdown
Масса_по_умолчанию = false #@param {type: "boolean"}
if Масса_по_умолчанию
m1, m2, m3 = 0.1, 0.1, 0.1
else
m1 = 0.1 # @param {type:"slider", min:0.01, max:2.0, step:0.01}
m2 = 0.1 # @param {type:"slider", min:0.01, max:2.0, step:0.01}
m3 = 0.1 # @param {type:"slider", min:0.01, max:2.0, step:0.01}
end
;
# @markdown # Длительность моделирования
# @markdown
t_begin = 0
t_end = 100.0 # @param {type:"number"}
;
Let's formulate the problem and solve the ODE system¶
# Поместим все начальные значения переменных в векторы
p = [m1, m2, m3]
u_begin = [ x1_0, y1_0, x2_0, y2_0, x3_0, y3_0,
vx1_0, vy1_0, vx2_0, vy2_0, vx3_0, vy3_0 ];
t_span = (t_begin, t_end)
# Сформулируем решаемую задачу
prob = ODEProblem( threebody, u_begin, t_span, p )
# Получим решение
sol = solve( prob, Vern7(), reltol = 1e-16, abstol=1e-16 );
Auxiliary code for graph output¶
Let's create a function that will extract the individual trajectories from the solution we obtained.
function get_values( sol; dt = 0.001, final = 1.0 )
final_time = sol.t[end] * final;
idx = sol.t[1]:dt:final_time
x1_values = map( value -> value[1], sol.(idx) )
y1_values = map( value -> value[2], sol.(idx) )
x2_values = map( value -> value[3], sol.(idx) )
y2_values = map( value -> value[4], sol.(idx) )
x3_values = map( value -> value[5], sol.(idx) )
y3_values = map( value -> value[6], sol.(idx) )
return(
[x1_values, x2_values, x3_values],
[y1_values, y2_values, y3_values]
)
end
This function will allow us to decompose the variables from the solution sol into several vectors, where the coordinates x and y will be collected.
x_values, y_values = get_values(sol);
Let's create some more variables that will contain the position of each body at the end of the simulation.
x1_end = x_values[1][end]; y1_end = y_values[1][end];
x2_end = x_values[2][end]; y2_end = y_values[2][end];
x3_end = x_values[3][end]; y3_end = y_values[3][end];
Conclusion the solution¶
It is possible to directly plot all coordinates separately and estimate where the objects are located, but this plot will be confusing and will take up a lot of space in the file (about 40 MB).
# plot( sol )
Instead, let's plot a few other graphs.
traj = plot( x_values, y_values,
c=[julia.red julia.purple julia.green],
legend=false,
title="Задача трех тел",
xaxis="X", yaxis="Y",
xlims=(-2.0, 2.0),
ylims=(-1.0, 2.0),
aspect_ratio = 1.1 )
scatter!( traj, [x1_end, x2_end, x3_end], [y1_end, y2_end, y3_end],
mc=[julia.red, julia.purple, julia.green], ms=7 )
The most frequent outcome of this experiment is that one of the bodies is separated from the rest, and the other two bodies form a binary system and continue to move in space around a common centre of attraction.
To better study the motion of the bodies, let us build an interactive simulation.
# @markdown # Модельное время
# @markdown
time_fraction = 0.307 # @param {type:"slider", min:0.0, max:1.0, step:0.001}
x_values, y_values = get_values(sol, dt=0.1, final=time_fraction);
x1_end = x_values[1][end]; y1_end = y_values[1][end];
x2_end = x_values[2][end]; y2_end = y_values[2][end];
x3_end = x_values[3][end]; y3_end = y_values[3][end];
traj = plot( x_values, y_values,
c=[julia.red julia.purple julia.green],
legend=false,
title="Задача трех тел",
xaxis="X", yaxis="Y",
xlims=(-2.0, 2.0),
ylims=(-1.0, 3.0),
aspect_ratio = 1.1 )
scatter!( traj, [x1_end, x2_end, x3_end], [y1_end, y2_end, y3_end],
mc=[julia.red, julia.purple, julia.green], ms=7 )
Conclusion¶
We have constructed a model of a system of three material points and at the same time we have shown how masking of code cells allows us to control the parameters and the solution of this differential system in time.
One can easily check that the smallest change of the model parameters (coordinates, velocity or mass) leads to a huge change of the final trajectories, which is the main feature of chaotic systems.