在物理学中，三体问题是经典力学中一个臭名昭著的问题，它涉及预测三个质量在彼此引力影响下的运动。

三体问题之所以如此臭名昭著，是因为它的混沌性质使其难以解决。虽然可以找到对应的二体解决方案的通用封闭形式解决方案，但不幸的是，当我们引入第三体时，我们无法做到同样的事情。

然而，尽管三体问题表现出这种几乎无法预测的混沌性质，但实际上使用数值方法相对容易模拟。在这篇文章中，我想介绍如何在 Python 上直观地模拟三体问题，以生成如下所示的模拟：

数值求解 ODE

在讨论三体问题之前，我将首先回顾一下数值求解常微分方程 (ODE) 的基础知识。

当涉及到数值求解 ODE 时，拥有一种表达任何类型 ODE 的通用方法会变得很有用。这是通过将一组任意阶的 ODE 转换为 N 个同时的一阶 ODE 来完成的。

例如，考虑牛顿第二定律（二阶常微分方程）：

为了将其转化为一阶 ODE，我们将做出两个定义：

我们还可以用向量表示法表示如下：

（请注意，我们从 0 而不是 1 开始，以便与 Python 基于 0 的索引保持一致）

使用这两个新定义，我们可以创建以下两个一阶 ODE：

然后我们可以将这两个方程的右侧定义为另一个向量函数 f 的分量，如下所示：

这有时被称为“力函数”。

因此，求解原始二阶 ODE 相当于求解以下一组 ODE：

使用这种紧凑形式很方便，因为现在我们有一种通用的方法来表示我们的 ODE 求解器始终可以理解的原始 ODE。因此，只要我们有一个可以求解一组一阶 ODE 的算法，我们就不需要太担心算法本身需要如何针对不同情况进行更改。

在以后的文章中，我可能会深入探讨特定的 ODE 算法，但为了防止这篇文章变得太长，我将假设我们已经有一个可以使用的工作 ODE 求解器。

二维中的两个物体

如果我们不确定自己在做什么，那么直接跳到三维空间中的三个物体的情况可能会出错，所以让我们从一个更简单的情况开始，只有两个物体的质量为二维质量 m1 和 m2。

尽管在 3D 中扩展到三个实体非常简单，但这有望使事情变得更容易理解。

设置 ODE

首先，让我们写下我们希望使用牛顿万有引力定律求解的常微分方程。为此，请考虑下图：

假设我们想要找到由于物体 2 作用在物体 1 上的力的 x 分量。因为我们知道要找到力的 x 分量，我们需要将 cosθ 乘以总力，我们可以写成 F = Gm1m2cosθ/r^2 其中 r 是两个物体之间的相对距离。然而，通过对上述三角形使用三角函数，我们还知道 cosθ = Δx/r，从而得出以下结论：

对主体 2 和 y 分量执行相同的操作，我们将得到以下四个需要求解的 ODE：

接下来，我们将设置 y 向量，如下所示：

（在某种程度上，y向量可以被认为是我们系统的状态向量，因为它拥有物体的位置和速度的所有信息）

使用 y 的分量，我们可以将四个二阶 ODE 重写为以下八个一阶 ODE：

如果我们将这八个 ODE 的右侧组合成一个向量函数 f，则我们的 ODE 的形式为 dy/dt = f。

用Python实现

通过在 2D 情况下为两个物体导出 y 和 f 向量，我们现在可以尝试在 Python 中实现这一点。这可以通过四个主要步骤来完成：

1. 设置初始条件，包括 y 和 f 函数

其操作如下所示：

from numpy import *

t = 0.; tf = 500. # Start and end time
h = 0.1 # Time increments
n = 8 # Number of components for y and force function
r1 = array([10., 0.]);   v1 = array([0., -0.1]); m1 = 1. # Initial conditions for body 1
r2 = array([-10., 0.]);  v2 = array([0., 0.1]); m2 = 1. # Initial conditions for body 2
y = concatenate((r1,r2,v1,v2))

def f(t, y):
    r = ((y[0]-y[2])**2+(y[1]-y[3])**2)**(1/2) # Relative distance
    return concatenate((y[4:8],
                        m2*(y[2:4]-y[0:2])/r**3,
                        m1*(y[0:2]-y[2:4])/r**3))

（请注意，我让 G = 1，因为使用 G 的实际值会使引力明显减弱）

虽然我现在已经输入了初始条件，但用户可以根据自己的需要进行调整。开始/结束时间和时间增量也是如此。

2. 设置视觉屏幕，我们将在其中看到两个物体的运动

为此，我将使用 VPython，但也可以在 matplotlib 中完成：

from vpython import *

scene = canvas(x=0,y=0,width = 500, height = 500) # Setting up the visual canvas
path1 = curve(color = color.blue, radius = 0.13) # Path of body 1
path2 = curve(color = color.red, radius = 0.13) # Path of body 2 
body1 = sphere(pos = vec(r1[0],r1[1],0), color = color.blue, radius = 0.5*m1**(1/3)) # Body 1
body2 = sphere(pos = vec(r2[0],r2[1],0), color = color.red, radius = 0.5*m2**(1/3)) # Body 2

3. 导入或复制并粘贴适当的 ODE 求解器

为了求解 ODE dy/dt = f，我们需要使用适当的 ODE 算法。

虽然许多 Python 库（例如 SciPy）已经内置了 ODE 求解器，但我将仅使用标准 Runge-Kutta (rk4) 算法作为我的 ODE 求解器。对于那些不了解rk4的人，只需将其视为欧拉方法的改进版本即可。

def rk4(t, h):
    ydumb = zeros(n, float); k1 = zeros(n, float); k2 = zeros(n, float)
    k3 = zeros(n, float);    k4 = zeros(n, float)
    for i in range(n):
        k1[i] = h*f(t,y)[i]
    for i in range(n):
        ydumb[i] = y[i] + k1[i]/2.
    k2 = h*f(t+h/2., ydumb)
    for i in range(n):
        ydumb[i] = y[i] + k2[i]/2.
    k3 = h*f(t+h/2., ydumb)
    for i in range(n):
        ydumb[i] = y[i] + k3[i]
    k4 = h*f(t+h, ydumb)
    for i in range(n):
        y[i] = y[i] + (k1[i] + 2.*(k2[i] + k3[i]) + k4[i])/6.
    return y

通过在当前时间 t 和增量 h 调用此函数，它将根据算法将状态向量 y 更新到下一个时间增量 t + h。

4. 创建一个 while 循环，为每个时间增量运行 ODE 求解器

最后，一切设置完毕后，我们现在可以运行一个 while 循环来重复求解 ODE，以找到两个物体在下一时刻的位置和速度。然后，这些结果将通过 VPython 在屏幕上更新：

while t < tf:
    rate(200) # Loops up to 200 times per seconds
    y = rk4(t,h) # Updates y according to the ODE algorithm
    t += h
    body1.pos = vec(y[0], y[1], 0) # Update the current positions of the bodies
    body2.pos = vec(y[2], y[3], 0)
    path1.append(pos = vec(y[0], y[1], 0)) # Update the paths the bodies took
    path2.append(pos = vec(y[2], y[3], 0))

扩展到三维空间中的三个物体

现在我们知道了需要采取的具体步骤，将其扩展到原始的三体问题应该很简单。

我们唯一需要注意的是，每个物体现在同时受到两种力，并且我们现在需要跟踪三个相对距离，而不仅仅是一个。

由于这两个事实，我们的力函数 f 现在如下所示：

def f(t, y):
    r12 = ((y[0]-y[3])**2+(y[1]-y[4])**2+(y[2]-y[5])**2)**(1/2) # Relative distance between 1 and 2
    r23 = ((y[3]-y[6])**2+(y[4]-y[7])**2+(y[5]-y[8])**2)**(1/2) # Relative distance between 2 and 3
    r31 = ((y[6]-y[0])**2+(y[7]-y[1])**2+(y[8]-y[2])**2)**(1/2) # Relative distance between 3 and 1
    return concatenate((y[9:18],
                m2*(y[3:6]-y[0:3])/r12**3+m3*(y[6:9]-y[0:3])/r31**3,
                m1*(y[0:3]-y[3:6])/r12**3+m3*(y[6:9]-y[3:6])/r23**3,
                m1*(y[0:3]-y[6:9])/r31**3+m2*(y[3:6]-y[6:9])/r23**3))

请注意，我们现在拥有三个相对距离，并且 fnow 中的最后 9 个元素有两项而不是只有一项。

如果我们然后调整其他所有内容以适应具有三个坐标的三个实体，我们将得到以下代码：

from vpython import *
from numpy import *

t0 = 0.; tf = 1000. # Start and end time
t = t0; h = 0.025
n = 18 # Number of components for y and force function
r1 = array([10., 0., 0.]);  v1 = array([-0.1, 0., 0.1]);   m1 = 1.
r2 = array([-10., 0., 0.]); v2 = array([0., 0.1, -0.1]);   m2 = 1.
r3 = array([0., 0., 5.]);   v3 = array([0.2, -0.1, 0.]);   m3 = 1.
y = concatenate((r1,r2,r3,v1,v2,v3))

scene = canvas(x=0,y=0,width = 700, height = 700)
path1 = curve(color = color.blue, radius = 0.13)
path2 = curve(color = color.red, radius = 0.13)
path3 = curve(color = color.green, radius = 0.13)
body1 = sphere(pos = vec(r1[0],r1[1],r1[2]), color = color.blue, radius = 0.5*m1**(1/3))
body2 = sphere(pos = vec(r2[0],r2[1],r2[2]), color = color.red, radius = 0.5*m2**(1/3))
body3 = sphere(pos = vec(r3[0],r3[1],r3[2]), color = color.green, radius = 0.5*m3**(1/3))

def f(t, y): # Force function
    r12 = ((y[0]-y[3])**2+(y[1]-y[4])**2+(y[2]-y[5])**2)**(1/2)
    r23 = ((y[3]-y[6])**2+(y[4]-y[7])**2+(y[5]-y[8])**2)**(1/2)
    r31 = ((y[6]-y[0])**2+(y[7]-y[1])**2+(y[8]-y[2])**2)**(1/2)
    return concatenate((y[9:18],
                m2*(y[3:6]-y[0:3])/r12**3+m3*(y[6:9]-y[0:3])/r31**3,
                m1*(y[0:3]-y[3:6])/r12**3+m3*(y[6:9]-y[3:6])/r23**3,
                m1*(y[0:3]-y[6:9])/r31**3+m2*(y[3:6]-y[6:9])/r23**3))

def rk4(t, h): # Runge-Kutta algorithm
    ydumb = zeros(n, float); k1 = zeros(n, float); k2 = zeros(n, float)
    k3 = zeros(n, float);    k4 = zeros(n, float)
    for i in range(n):
        k1[i] = h*f(t,y)[i]
    for i in range(n):
        ydumb[i] = y[i] + k1[i]/2.
    k2 = h*f(t+h/2., ydumb)
    for i in range(n):
        ydumb[i] = y[i] + k2[i]/2.
    k3 = h*f(t+h/2., ydumb)
    for i in range(n):
        ydumb[i] = y[i] + k3[i]
    k4 = h*f(t+h, ydumb)
    for i in range(n):
        y[i] = y[i] + (k1[i] + 2.*(k2[i] + k3[i]) + k4[i])/6.
    return y

while t < tf:
    rate(500)
    y = rk4(t,h)
    t += h
    body1.pos = vec(y[0], y[1], y[2])
    body2.pos = vec(y[3], y[4], y[5])
    body3.pos = vec(y[6], y[7], y[8])
    path1.append(pos = vec(y[0], y[1], y[2]))
    path2.append(pos = vec(y[3], y[4], y[5]))
    path3.append(pos = vec(y[6], y[7], y[8]))

运行此代码会产生以下模拟：

正如您所看到的，这三个物体在重力作用下明显相互吸引，但它们的路径似乎极其复杂，几乎无法预测。这正是我们想要模拟的混沌本质。

当然，如果您愿意，您可以尝试使用初始条件并查看系统演化的各种方式（如果您注意到物体突然飞走，这是因为它们彼此靠得太近，导致程序分裂接近 0 的数字。这可以通过使用更精确的 ODE 求解器或减小时间增量来解决