三体问题：方程与周期轨道

三体问题的数学建模：从方程到周期轨道

三体问题涉及根据牛顿万有引力定律计算三个大质量天体的轨迹。初始条件指定了三维空间中各天体的位置和速度，从而产生一组 18 个非线性二阶微分方程。

在完整问题中，各天体的坐标和速度与取决于它们相互距离的加速度相关联。归一化质量（相对于最大质量的比例）简化了计算。以下是以一阶导数形式表示的方程组：

X1'  = VX1
Y1'  = VY1
Z1'  = VZ1
X2'  = VX2
Y2'  = VY2
Z2'  = VZ2
X3'  = VX3
Y3'  = VY3
Z3'  = VZ3
VX1' = -m2*(X1-X2)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5-m3*(X1-X3)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5
VY1' = -m2*(Y1-Y2)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5-m3*(Y1-Y3)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5
VZ1' = -m2*(Z1-Z2)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5-m3*(Z1-Z3)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5
VX2' = -m3*(X2-X3)/((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2)^1.5-m1*(X2-X1)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5
VY2' = -m3*(Y2-Y3)/((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2)^1.5-m1*(Y2-Y1)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5
VZ2' = -m3*(Z2-Z3)/((X2-X3)^2+(Y2-Y3)^2+(Z2-Z3)^2)^1.5-m1*(Z2-Z1)/((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2)^1.5
VX3' = -m1*(X3-X1)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5-m2*(X3-X2)/((X3-X2)^2+(Y3-Y2)^2+(Z3-Z2)^2)^1.5
VY3' = -m1*(Y3-Y1)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5-m2*(Y3-Y2)/((X3-X2)^2+(Y3-Y2)^2+(Z3-Z2)^2)^1.5
VZ3' = -m1*(Z3-Z1)/((X1-X3)^2+(Y1-Y3)^2+(Z1-Z3)^2)^1.5-m2*(Z3-Z2)/((X3-X2)^2+(Y3-Y2)^2+(Z3-Z2)^2)^1.5

简化的变体：平面问题和限制问题

平面问题（P3BP）限制在二维空间，方程减少到 12 个。在此，天体保持在初始条件定义的平面内。

平面问题的方程：

X1'  = VX1
X2'  = VX2
X3'  = VX3
Y1'  = VY1
Y2'  = VY2
Y3'  = VY3
VX1' = m2*(X2-X1)/((X1-X2)^2+(Y1-Y2)^2)^1.5+m3*(X3-X1)/((X1-X3)^2+(Y1-Y3)^2)^1.5
VX2' = m3*(X3-X2)/((X2-X3)^2+(Y2-Y3)^2)^1.5+m1*(X1-X2)/((X2-X1)^2+(Y2-Y1)^2)^1.5
VX3' = m1*(X1-X3)/((X3-X1)^2+(Y3-Y1)^2)^1.5+m2*(X2-X3)/((X3-X2)^2+(Y3-Y2)^2)^1.5
VY1' = m2*(Y2-Y1)/((X1-X2)^2+(Y1-Y2)^2)^1.5+m3*(Y3-Y1)/((X1-X3)^2+(Y1-Y3)^2)^1.5
VY2' = m3*(Y3-Y2)/((X2-X3)^2+(Y2-Y3)^2)^1.5+m1*(Y1-Y2)/((X2-X1)^2+(Y2-Y1)^2)^1.5
VY3' = m1*(Y1-Y3)/((X3-X1)^2+(Y3-Y1)^2)^1.5+m3*(Y2-Y3)/((X3-X2)^2+(Y3-Y2)^2)^1.5

限制问题（PR3BP）假设第三个天体质量很小。其运动方程：

X' = VX
Y' = VY
VX' = X+2*Y4-m1*(X+m2)/((X+m2)^2+Y^2)^1.5-m2*(X-m1)/((X-m1)^2+Y^2)^1.5
VY' = Y-2*Y3-m1*Y/((X+m2)^2+Y^2)^1.5-m2*Y/((X-m1)^2+Y^2)^1.5

圆周限制问题（CR3BP）模拟航天器在两个天体系统中，其中质量较小的天体处于圆轨道。使用参数 μ = m2 / (m1 + m2)：

X' = VX
Y' = VY
z' = VZ
VX' = 2*VY+X-(1-mu) * (X+mu)/((X+mu)^2+Y^2+Z^2)^1.5-mu * (X-1+mu)/((X-1+mu)^2+Y^2+Z^2)^1.5
VY' = -2*VX+Y-(1-mu) * Y/((X+mu)^2+Y^2+Z^2)^1.5-mu * Y/((X-1+mu)^2+Y^2+Z^2)^1.5
VZ' = -(1-mu) * Z/((X+mu)^2+Y^2+Z^2)^1.5-mu * Z/((X-1+mu)^2+Y^2+Z^2)^1.5

这些系统使用 Mathematica、Python（SciPy）或 MATLAB 中的数值方法求解。

解析解：欧拉和拉格朗日

已知的唯一解析解是欧拉的共线解（1767 年）和拉格朗日的三角解（1772 年）。在拉格朗日情形下，天体从等边三角形的顶点出发，速度匹配。轨迹是旋转 120° 的椭圆；天体保持三角形配置。

这些解是不稳定的：初始条件的微小偏差会导致混沌。实际实现需要达到 10^{-10} 的精度。

平面问题中的周期轨道

亨利·庞加莱证明了等质量情况下存在无穷多个周期解。计算机搜索发现了具有对称性的族。

• Broucke-Enò-Hadjidimitriou 族 (1975)：8 符号和 16 符号轨道。示例 R4（周期 5.4）：

```

T0 = 0, T1 = 5.4

X1(T0) = 0.8733047091, X2(T0) = -0.6254030288, X3(T0) = -0.2479016803

Y1(T0) = 0, Y2(T0) = 0, Y3(T0) = 0

VX1(T0) = 0, VX2(T0) = 0, VX3(T0) = 0

VY1(T0) = 1.010776444, VY2(T0) = -1.683353346, VY3(T0) = 0.6725769022

```

• Shuvakov-Dmitrashinovich 族 (2013)：蜻蜓（周期 22）：

```

T0 = 0, T1 = 22

X1(T0) = -1, X2(T0) = 1, X3(T0) = 0

Y1(T0) = 0, Y2(T0) = 0, Y3(T0) = 0

VX1(T0) = 0.08058, VX2(T0) = 0.08058, VX3(T0) = -0.16116

VY1(T0) = 0.58884, VY2(T0) = 0.58884, VY3(T0) = -1.17768

```

其他：蝴蝶、蛾子、护目镜、阴阳。

• Shin 族 (2016)：椭圆形、猫脸（周期 5.1）。

要重现，使用带自适应步长的 8 阶龙格-库塔法。

要点

• 完整问题需要 18 个方程；简化后减少到 4–6 个，用于实际模拟。
• 欧拉/拉格朗日解析解不稳定，仅适用于理想条件。
• 周期轨道使用数值方法发现；Broucke 和 Shuvakov 族包含数百种变体。
• 为实现周期性，初始条件必须高精度指定。
• CR3BP 对 NASA/ESA 任务至关重要：地月和日木系统中的轨道机动。

— Editorial Team