平面三体问题（Julia版）
评论
# 奇怪：这个库必须放在最前面才能一次加载成功
using Plots
# gr()
pyplot()
using SciPy
using LinearAlgebra
Plots.GRBackend()
在平面上引入直角坐标系，每个星体的坐标和速度分别为：
\[ x_i \in \mathbb{R}^{2} \ \ \dot{x}_i \in \mathbb{R}^{2} \ \ i=1,2,3 \]直接根据万有引力定律和牛顿第二定律，可写出下面的三体动力学方程组（其中，引力常数可通过选择适当的单位制使得 \(G=1\) ）：
\[ \begin{aligned}\frac{d x_1}{dt} &=\dot{x}_1 \\ \frac{d x_2}{dt} &=\dot{x}_2 \\ \frac{d x_3}{dt} &=\dot{x}_3 \\ \frac{d \dot{x}_1}{dt} &=m_2 \frac{x_2-x_1}{r_3^{\ \ 3}} + m_3 \frac{x_3-x_1}{r_2^{\ \ 3}} \\ \frac{d \dot{x}_2}{dt} &=m_3 \frac{x_3-x_2}{r_2^{\ \ 3}} + m_1 \frac{x_1-x_2}{r_3^{\ \ 3}} \\ \frac{d \dot{x}_3}{dt} &=m_1 \frac{x_1-x_3}{r_2^{\ \ 3}} + m_2 \frac{x_2-x_3}{r_1^{\ \ 3}} \end{aligned} \]下面的程序中按如下次序装配成12个元素的数组u：
\[ u=[x_1,\dot{x}_1,x_2,\dot{x}_2,x_3,\dot{x}_3] \in \mathbb{R}^{12} \]function threebody(t, u, m1,m2,m3)
# 3个星体各自的坐标
x1 = u[1:2]
x2 = u[5:6]
x3 = u[9:10]
# \vec{r}/r^3
d1 = (x3-x2)/norm(x3-x2)^3
d2 = (x1-x3)/norm(x1-x3)^3
d3 = (x2-x1)/norm(x2-x1)^3
# du
du = zeros(12)
du[1:2] = u[3:4]
du[5:6] = u[7:8]
du[9:10] = u[11:12]
du[3:4] = m2*d3 - m3*d2;
du[7:8] = m3*d1 - m1*d3;
du[11:12] = m1*d2 - m2*d1
du
end;
初始静止的平面三体，其位置分别在(1,-1) (1,3) (-2,-1)。
u0 = [1.0;-1.0;0.0;0.0;1.0;3.0;0.0;0.0;-2.0;-1.0;0.0;0.0]
tspan = [0,40]
# 调用ode进行求解 odeint
sol = SciPy.integrate.solve_ivp(threebody, tspan, u0, args=(5.0, 3.0, 4.0))
t = sol["t"]
u = sol["y"];
绘制运动过程
@gif for i in eachindex(t)
plot(u[1,1:i], u[2,1:i], color = "red", label = "")
plot!(u[5,1:i], u[6,1:i], color = "green", label = "")
plot!(u[9,1:i], u[10,1:i], color = "blue", label = "")
scatter!((u[1,i], u[2,i]), color = "red", label = "")
scatter!((u[5,i], u[6,i]), color = "green", label = "")
scatter!((u[9,i], u[10,i]), color = "blue", label = "")
end every 5
┌ Info: Saved animation to
│ fn = F:\谷歌云端硬盘\Colab Notebooks\tmp.gif
└ @ Plots C:\Users\DELL\.julia\packages\Plots\5ItHH\src\animation.jl:104
如果选择三体的质心为坐标原点，那么在此坐标系下，知道三体任意两体的坐标，就知道第三个的坐标：
\[ x_3 = -\frac{1}{m_3}(m_1 x_1+m_2 x_2) \]进而只需要求解关于 \(x_1,x_2\) 的方程组，此时可重新装配待求变量的数组：
\[ u=[x_1,\dot{x}_1,x_2,\dot{x}_2] \in \mathbb{R}^{8} \]于是有：
(m1,m2,m3) = (5,3,4);
function threebody2(t, u, m1,m2,m3)
# 3个星体各自的坐标
x1 = u[1:2]
x2 = u[5:6]
x3 = -(m1*x1+m2*x2)/m3
# \vec{r}/r^3
d1 = (x3-x2)/norm(x3-x2)^3
d2 = (x1-x3)/norm(x1-x3)^3
d3 = (x2-x1)/norm(x2-x1)^3
# du (只需要其中两体的动力学方程)
du = zeros(8)
du[1:2] = u[3:4]
du[5:6] = u[7:8]
du[3:4] = m2*d3 - m3*d2;
du[7:8] = m3*d1 - m1*d3;
du
end;
# 为了方便比较，确保初值条件和前面完全一致
u0 = [1.0;-1.0;0.0;0.0;1.0;3.0;0.0;0.0;-2.0;-1.0;0.0;0.0]
o = (m1*u0[1:2]+m2*u0[5:6]+m3*u0[9:10])/(m1+m2+m3)
u0 = u0[1:8]
u0[1:2] -= o
u0[5:6] -= o
tspan = [0,40]
# 调用ode进行求解 odeint
sol = SciPy.integrate.solve_ivp(threebody2, tspan, u0, args=(m1,m2,m3))
t = sol["t"]
u = sol["y"];
x1 = u[[1,2],:];
x2 = u[[5,6],:];
x3 = -(m1*x1+m2*x2)/m3
@gif for i in eachindex(t)
plot(x1[1,1:i], x1[2,1:i], color = "red", label = "")
plot!(x2[1,1:i], x2[2,1:i], color = "green", label = "")
plot!(x3[1,1:i], x3[2,1:i], color = "blue", label = "")
scatter!((x1[1,i], x1[2,i]), color = "red", label = "")
scatter!((x2[1,i], x2[2,i]), color = "green", label = "")
scatter!((x3[1,i], x3[2,i]), color = "blue", label = "")
end every 5
┌ Info: Saved animation to
│ fn = F:\谷歌云端硬盘\Colab Notebooks\tmp.gif
└ @ Plots C:\Users\DELL\.julia\packages\Plots\5ItHH\src\animation.jl:104
我们注意到，原本给出的初值条件的质心本来就在原点：
o
2-element Array{Float64,1}:
0.0
0.0
所以按理说，前后两次的绘图应该完全一样才对，但事实上不一样。其实，这暗含三体运动的不可预测性（数值计算的微小舍入误差可能会被无限放大）。