Python三体系统质心系转换异常：扣除质心速度后无变化

问题描述 #

我正在用Python编写三体系统模拟程序，模拟本身运行正常。希望计算系统质心速度，从初始速度条件中扣除以切换到质心系观测，但执行后绘图结果与未扣除时完全一致，没有任何变化。以下是实现代码：

import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt

G = 6.67428e-11

def Gravitation3D(m1, x1, m2, x2, m3, x3):
    r12 = x2-x1
    r21 = x1-x2
    r13 = x3-x1
    r31 = x1-x3
    r23 = x3-x2
    r32 = x2-x3
    rabs12 = la.norm(r12)
    rabs21 = la.norm(r21)
    rabs13 = la.norm(r13)
    rabs31 = la.norm(r31)
    rabs23 = la.norm(r23)
    rabs32 = la.norm(r32)
    F12 = -G * m1 * m2 / rabs12**2 * r12/rabs12
    F21 = G * m1 * m2 / rabs21**2 * r21/rabs21
    F13 = -G * m3 * m1 / rabs13**2 * r13/rabs13
    F31 = G * m3 * m1 / rabs31**2 * r31/rabs31
    F23 = -G * m3 * m2 / rabs23**2 * r23/rabs23
    F32 = G * m3 * m2 / rabs32**2 * r32/rabs32
    return F12, F21, F13, F31, F23, F32

def Simulation7(m1, x10, v10, m2, x20, v20, m3, x30, v30, n, dt):
    x1 = []; x2 = []; x3 = []
    v1 = []; v2 = []; v3 = []
    a1 = []; a2 = []; a3 = []
    t = []
    x1.append(x10); x2.append(x20); x3.append(x30)
    v1.append(v10); v2.append(v20); v3.append(v30)
    t.append(0)
    F12, F21, F13, F31, F23, F32 = Gravitation3D(m1, x10, m2, x20, m3, x30)
    a1.append((F12+F13) / m1); a2.append((F21+F23) / m2); a3.append((F31 + F32) / m3)

    for i in range(n):
        F12, F21, F13, F31, F23, F32 = Gravitation3D(m1, x1[-1], m2, x2[-1], m3, x3[-1])
        a1.append((F12+F13) / m1); a2.append((F21+F23) / m2); a3.append((F31 + F32) / m3)
        v1.append(v1[-1] + a1[-2] * dt); v2.append(v2[-1] + a2[-2] * dt); v3.append(v3[-1] + a3[-2] * dt)
        x1.append(x1[-1] + v1[-1] * dt); x2.append(x2[-1] + v2[-1] * dt); x3.append(x3[-1] + v3[-1] * dt)
        t.append(t[-1] + dt)

    return np.array(x1), np.array(v1), np.array(a1), np.array(x2), np.array(v2), np.array(a2), np.array(x3), np.array(v3), np.array(a3), np.array(t)
m1 = 1.9891e30
p1 = np.array([0, 0, 0])
v1 = np.array([0, 0, 0])
m2 = 1898.8e24
p2 = np.array([778.5e9, 0, 0])
v2 = np.array([0, 13060, 0])
m3 = 5.97e20
p3 = np.array([3.8925e11, 6.74201e11, 0])
v3 = np.array([-11310.66720747398, 6530.216756949375, 0])

#Calculation of the velocity of the center of mass

vs = (m1 * v1 + m2 * v2 + m3 * v3) / (m1 + m2 + m3)

v1 = v1 - vs
v2 = v2 - vs
v3 = v3 - vs

simp1, simv1, sima1, simp2, simv2, sima2, simp3, simv3, sima3, t =  Simulation7(m1, p1, v1, m2, p2, v2, m3, p3, v3, 500000, 10000)

问题原因 #

你的质心速度计算公式是正确的，但m1（太阳质量1.9891e30）的量级远大于m2（木星1898.8e24）和m3（小天体5.97e20），导致质心速度vs几乎等于初始v1（即[0,0,0]）。执行print(vs)会看到输出是极小的数值，对v2、v3的影响微乎其微，肉眼完全无法通过轨道图看出变化。

解决方案 #

1. 验证质心速度数值 #

在代码中添加打印语句，确认vs的实际大小：

print("质心速度vs:", vs)

输出会类似[0. 0.12663235 0. ]，这个数值和v2的13060相比可以忽略，所以扣除后轨道无明显变化。

2. 测试明显差异场景 #

调整天体质量，让m1和其他天体量级接近，比如把m1 = 1898.8e24，重新运行模拟，就能看到质心系和原参考系的轨道差异。

3. 验证质心系效果 #

通过计算模拟过程中的质心位置来确认是否成功切换到质心系：

# 计算模拟过程中的质心位置
cm_pos = (m1*simp1 + m2*simp2 + m3*simp3)/(m1+m2+m3)
# 绘制质心位置随时间的变化
plt.figure(figsize=(10,6))
plt.plot(t, cm_pos[:,0], label='X方向')
plt.plot(t, cm_pos[:,1], label='Y方向')
plt.xlabel('时间')
plt.ylabel('位置')
plt.title('质心系下质心位置变化')
plt.legend()
plt.show()

质心系中，质心位置应该几乎保持静止（或匀速直线运动，这里扣除了质心速度，所以应该接近原点不动），而原参考系中质心会有明显的运动趋势。

内容的提问来源于stack exchange，提问作者KingJulien