Задача трех тел¶
Реализуем численное решение задачи симуляции динамики трех тел, находящихся в гравитационном взаимодействии друг с другом.
Описание задачи¶
Три материальные точки с координатами $\bar x_1$, $\bar x_2$, $\bar x_3$ и скоростями $\bar v_1$, $\bar v_2$, $\bar v_3$ находятся в едином пространстве. Нужно определить их траектории во времени.
Как мы знаем, к сожалению, пока не существует аналитического решения задачи расчета координат всех объектов системы в любой момент времени t. Но мы решим эту задачу при помощи численного интегрирования.
$$\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}$ здесь обозначают радиус-векторы каждого тела, $\mathbf{\ddot r_1} \dots \mathbf{\ddot r_3}$ – их вторые производные (ускорение), $m_1 \dots m_3$ – масса каждого тела, гравитационную постоянную $G$ можно принять равной 1 без потерь для наших нужд по визуализации динамики системы. Обозначение $|{\mathbf{r_1}} - {\mathbf{r_2}}|$ подразумевает норму между радиус-векторами, проще говоря – расстоянии между телами 1 и 2.
Итак, гравитация каждого из двух тел влияет на ускорение третьего. Мы сформулировали уравнения, для которых нет аналитического решения, к тому же при $G=1$ наше решение будет абстрактным и не будет иметь правильных размерностей.
Решение системы дифференциальных уравнений¶
Для решения этой системы из трех дифференциальных уравнений второго порядка нам нужно будет также решить уравнения первого порядка. Воспользуемся библиотекой DifferentialEquations:
using DifferentialEquations
gr()
julia = Colors.JULIA_LOGO_COLORS # Пусть графики будут в цветах логотипа Julia
Решение на текущем шаге и обновление состояния системы мы создадим функцию 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
Можно заметить, что векторизация значительно укоротила бы наш код, также как и использование функции norm() из библиотеки LinearAlgebra для вычисления расстояний. Но в этом учебном примере мы пока отдадим предпочтение визуальному сходству кода и формул, которое есть в текущей реализации.
Задаем начальные параметры¶
У нас большой вектор начальных координат. Добавим несколько полей ввода чтобы было легко изменять эти параметры.
В качестве значений начальной скорости мы возьмем следующие числа (в том или ином направлении):
$$\begin{align} vx_0 = -0.148 \cdot \cos( 60^\circ ) \\ vy_0 = -0.148 \cdot \sin( 60^\circ ) \end{align}$$
Замечание:
1) Следующая ячейка автоматически не выполняется, чтобы обновить значения переменных, нажмите на кнопку ⏵, или нажмите в свободную область маски и запустите ячейку через клавиатурное сокращение Ctrl+Enter.
2) Чтобы скрыть код и оставить только элементы управления, нужно дважды кликнуть по маске
# @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"}
;
Сформулируем задачу и решим систему ОДУ¶
# Поместим все начальные значения переменных в векторы
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 );
Вспомогательный код для вывода графиков¶
Создадим функцию, которая будет извлекать отдельные траектории из решения, которое мы получили.
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
Эта функция позволит нам разложить переменные из решения sol в несколько векторов, где будут собраны координаты x и y.
x_values, y_values = get_values(sol);
Создадим еще несколько переменных, в которых будет находиться положение каждого тела в момент окончания симуляции.
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];
Выведем решение¶
Можно напрямую построить график всех координат по-отдельности и оценить, где находятся объекты, но этот график будет запутанным и займет очень много места в файле (около 40 Мб).
# plot( sol )
Вместо этого выведем несколько других графиков.
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 )
Наиболее частый исход этого эксперимента состоит в том, что одно из тел отделяется от остальных, а два других образуют бинарную систему и продолжают перемещаться в пространстве вокруг общего центра притяжения.
Чтобы лучше изучить движение тел, построим интерактивную симуляцию.
# @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 )
Заключение¶
Мы построили модель системы из трех материальных точек и заодно показали, как маскирование кодовых ячеек позволяет управлять параметрами и решением этой дифференциальной системы во времени.
Можно легко проверить, что самое незначительное изменение параметров модели (координат, скорости или массы) приводит к огромному изменению итоговых траекторий, что и является основной чертой хаотических систем.