三体问题
我们实现了一个数值解的问题，模拟三个体在重力相互作用的动力学。
任务说明
具有坐标的三个材料点 , , 和速度 , , 它们位于单个空间中。 有必要及时确定它们的轨迹。
正如我们所知，不幸的是，在任何给定时间计算系统中所有物体坐标的问题没有分析解决方案。 t. 但是我们将使用数值积分来解决这个问题。
这里，表示各体的半径矢量。, -它们的二阶导数（加速度）, -每个物体的质量，引力常数 它可以假设等于1而不损失我们需要可视化系统的动态。 指定；指定 它意味着半径矢量之间的范数，换句话说，体1和2之间的距离。
因此，两个物体中的每一个的重力都会影响第三个物体的加速度。 我们已经制定了方程，其中没有解析解，此外，当 我们的解决方案将是抽象的，不会有正确的维度。
求解微分方程组
为了求解这个由三个二阶微分方程组成的系统，我们还需要求解一阶方程。 让我们使用图书馆 DifferentialEquations:
In [ ]:
Pkg.add(["DifferentialEquations"])
In [ ]:
using DifferentialEquations
gr()
julia = Colors.JULIA_LOGO_COLORS # Пусть графики будут в цветах логотипа Julia
Out[0]:
解决当前步骤并更新系统状态，我们将创建一个函数 thereebody.
In [ ]:
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
Out[0]:
您可以看到矢量化会显着缩短我们的代码，以及使用该函数 norm() 从图书馆 LinearAlgebra 来计算距离。 但是在本教程中，就目前而言，我们将优先考虑代码和公式的视觉相似性，这在当前实现中。
设置初始参数
我们有一个很大的初始坐标向量。 我们将添加几个输入字段，以便于更改这些参数。
我们将使用以下数字（在一个方向或另一个方向）作为初始速度值:
备注:
1）下一个单元格不会自动执行，要更新变量的值，请按⏵按钮，或单击蒙版的空闲区域，并使用键盘快捷键Ctrl+Enter启动单元格。
2）要隐藏代码并只留下控件，请双击掩码。
In [ ]:
# @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 )
Out[0]:
In [ ]:
# @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 )
Out[0]:
In [ ]:
# @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
;
In [ ]:
# @markdown # Длительность моделирования
# @markdown
t_begin = 0
t_end = 100.0 # @param {type:"number"}
;
让我们制定的问题和解决系统的颂歌
In [ ]:
# Поместим все начальные значения переменных в векторы
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 );
显示图形的辅助代码
让我们创建一个函数，从我们收到的解决方案中提取单个轨迹。
In [ ]:
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
Out[0]:
这个函数将允许我们从解中分解变量。 sol 在几个向量，其中坐标将被收集 x 和 y.
In [ ]:
x_values, y_values = get_values(sol);
让我们创建更多的变量，这些变量将包含模拟结束时每个主体的位置。
In [ ]:
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];
让我们推导出解决方案
您可以直接单独绘制所有坐标并估计对象所在的位置，但此图会令人困惑，并且会占用文件中的大量空间（约40MB）。
In [ ]:
# plot( sol )
相反，我们将显示其他几个图表。
In [ ]:
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 )
Out[0]:
这个实验最常见的结果是其中一个身体与其他身体分开，而另外两个形成一个二元系统，并继续围绕一个共同的重心在空间中移动。
为了更好地研究身体的运动，我们将建立一个交互式模拟。
In [ ]:
# @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 )
Out[0]:
结论
我们建立了一个由三个材料点组成的系统模型，并展示了掩蔽代码单元如何使我们能够随时间控制这个差分系统的参数和解决方案。
您可以轻松检查模型参数（坐标，速度或质量）的最小变化导致最终轨迹的巨大变化，这是混沌系统的主要特征。