基于四阶Runge-Kutta的太阳系行星运动模拟正确性排查

我用四阶Runge-Kutta（RK4）方法模拟太阳系行星运动，设置时间步长dt=3600*24（1天，单位秒），运行365步后预期地球完成一圈公转，但结果连四分之一轨道都没达到。必须增加步数才能完成公转，但增大dt不符合当前单位逻辑。以下是代码实现及调用逻辑，请问我的RK4方法应用是否正确？

行星类代码 #

class Planet:
    G = 6.67430e-11
    def __init__(self, name, position, velocity, mass):
        self.name = name
        self.position = position
        self.velocity = velocity
        self.mass = mass
        self.dr = np.array([0,0,0])
        self.dv = np.array([0,0,0])
        self.pos_data = [self.position.copy()]
        self.vel_data = [self.position.copy()]

    def acceleration(self, pos, planets):
        acceleration = np.array([0,0,0])
        for planet in planets:
            if planet != self:
                r = planet.position - pos
                r_mag = np.linalg.norm(r)
                force_magnitude = (self.G * self.mass * planet.mass) / (r_mag ** 2)
                accel_magnitude = force_magnitude / self.mass
                acceleration = acceleration + accel_magnitude * (r / r_mag)
        return acceleration

    def rk4(self, dt, planets):
        k1v = dt * self.acceleration(self.position, planets)
        k1r = dt * self.velocity        
        intermediate_position = self.position + 0.5 * k1r
        k2v = dt * self.acceleration(intermediate_position, planets)
        k2r = dt * (0.5 * k1v)        
        intermediate_position = self.position + 0.5 * k2r
        k3v = dt * self.acceleration(intermediate_position, planets)
        k3r = dt * (0.5 * k2v)        
        intermediate_position = self.position + k3r
        k4v = dt * self.acceleration(intermediate_position, planets)
        k4r = dt * (k3v)
        self.dp = (k1r + 2 * k2r + 2 * k3r + k4r) / 6
        self.dv = (k1v + 2 * k2v + 2 * k3v + k4v) / 6
    
    def update_state(self):
        self.position = self.position + self.dp
        self.velocity = self.velocity + self.dv
        self.pos_data.append(self.position)
        self.vel_data.append(self.velocity)

调用代码 #

sol     = Planet("Sun"    , np.array([0.0       , 0.0, 0.0]), np.array([0.0, 0.0    , 0.0]), 1.989e30)
mercury = Planet('Mercury', np.array([0.39e12   , 0.0, 0.0]), np.array([0.0, 47000.0, 0.0]), 3.285e23)
venus   = Planet('Venus'  , np.array([0.72e12   , 0.0, 0.0]), np.array([0.0, 35000.0, 0.0]), 4.867e24)
earth   = Planet('Earth'  , np.array([1.0e12    , 0.0, 0.0]), np.array([0.0, 29783.0, 0.0]), 5.972e24)
moon    = Planet('Moon'   , np.array([1.08e12    , 0.0, 0.0]), np.array([0.0, 30805.0, 0.0]), 7.348e22)

planets = [sol, earth, venus, moon, mercury]

dt = 3600*24 # 1 day
for _ in range(365*10):
    for planet in planets:
        planet.rk4(dt, planets)
    for planet in planets:
        planet.update_state()

plt.clf()
for planet in planets:
    x = [pos[0] for pos in planet.pos_data]
    y = [pos[1] for pos in planet.pos_data]
    z = [pos[2] for pos in planet.pos_data]
    plt.scatter(x[-1], y[-1])
    plt.plot(x,y)
plt.scatter(0,0)
plt.show()

你的核心问题是RK4方法的实现完全错误，尤其是位移增量的计算部分，导致行星运动速度被严重低估，轨道推进极慢。此外还有几处次要错误，逐一说明：

1. RK4位移增量（k2r、k3r、k4r）计算错误 #

标准RK4的位移增量计算逻辑是基于当前速度加上对应速度增量的一半/全部，而你错误地直接用了速度增量的一半来计算位移增量，完全偏离了RK4的核心公式：

错误代码： #

k2r = dt * (0.5 * k1v)        
k3r = dt * (0.5 * k2v)        
k4r = dt * (k3v)

修正后代码： #

k2r = dt * (self.velocity + 0.5 * k1v)  # 用当前速度+半份k1v计算中间位移
k3r = dt * (self.velocity + 0.5 * k2v)  # 用当前速度+半份k2v计算中间位移
k4r = dt * (self.velocity + k3v)        # 用当前速度+全份k3v计算中间位移

这是导致轨道推进过慢的直接原因——你之前的计算相当于把速度增量当成了速度本身，实际运动速度被大幅缩小。

2. 加速度计算冗余（不影响结果，但可优化） #

你计算加速度时先算了合力再除以自身质量，其实可以直接简化为：

accel_magnitude = self.G * planet.mass / (r_mag ** 2)

因为F=G*M*m/r²，a=F/m=G*M/r²，这样可以避免不必要的乘法，提升效率。

3. 初始化时vel_data存储错误 #

__init__里的self.vel_data = [self.position.copy()]应该改为self.vel_data = [self.velocity.copy()]，否则你存储的速度数据实际是位置数据，不过这不会影响轨道模拟，只是后续分析速度时会出错。

4. 地球轨道半径与实际值不符（可选修正） #

你设置的地球轨道半径是1.0e12米，而实际天文单位是1.496e12米。如果要更贴近真实太阳系，建议修正这个数值，否则地球的轨道周期会比实际更短（不过这不是你当前“365步完不成公转”的原因）。

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