三体问题¶
让我们来实现模拟三个物体在引力作用下的动力学问题的数值解决方案。
问题描述¶
坐标分别为$\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 库：
Pkg.add(["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
您可以看到，矢量化将大大缩短我们的代码，使用LinearAlgebra 库中的norm() 函数计算距离也是如此。但在本教程示例中，我们将暂时倾向于当前实现中代码和公式的视觉相似性。
设置初始参数¶
我们有一个很大的初始坐标向量。让我们添加一些输入字段，方便更改这些参数。
作为初始速度值，我们将使用下列数字（在一个或另一个方向上）：
$$\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"}
;
让我们提出问题并求解 ODE 系统¶
# Поместим все начальные значения переменных в векторы
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 MB）。
# 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 )
结论¶
我们已经构建了一个由三个材料点组成的系统模型，同时我们还展示了如何通过屏蔽代码单元来控制这个微分系统的参数和及时求解。
我们可以很容易地发现，模型参数（坐标、速度或质量）的最小变化都会导致最终轨迹的巨大变化，而这正是混沌系统的主要特征。