N 体问题 - 我应用龙格-库塔四阶方法来模拟太阳系中行星的运动是否正确？

Question

当我将时间步长设置为 dt = 3600*24（1 天，以秒为单位）并运行 365 步的模拟时，我希望地球绕太阳旋转一圈。然而，结果甚至不接近四分之一轨道。 对于我选择的 dt 值，我需要绕一个地球轨道运行的步数比我预期的要多，我可以使 dt 更大，但在我使用的单位中没有意义。

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()

也许我错了，不应该期望它是准确的，但我认为它至少应该接近。任何帮助都会很棒。谢谢。

这是我的第一个问题，如果我的问题格式错误，抱歉。

Answer 1

我不是这方面的专家，所以代码可能还可以进一步改进。但主要问题是：

1. k2r

   、

   k3r

   和

   k4r

   的计算应包含更新的速度，而不是速度更新的大小。

所以代替：

k2r = dt * (0.5 * k1v) 

k3r = dt * (0.5 * k2v) 

k4r = dt * (k3v)

应该是：

k2r = dt * (self.velocity + 0.5 * k1v)

k3r = dt * (self.velocity + 0.5 * k2v)

k4r = dt * (self.velocity + k3v)

2. 初始距离已偏离。地球与太阳的距离（即1天文单位）约为1.50e11米。当前代码中的 1e12 与此相差 6.67 倍。其他行星距离似乎有相同比例的偏差，而地月距离则相当大（应该约为 0.00266 天文单位或 3.84e8 米）。

将设置代码更正为：

au = 1.50e11

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.39*au   , 0.0, 0.0]), np.array([0.0, 47000.0, 0.0]), 3.285e23)
venus   = Planet('Venus'  , np.array([0.72*au   , 0.0, 0.0]), np.array([0.0, 35000.0, 0.0]), 4.867e24)
earth   = Planet('Earth'  , np.array([1.0*au    , 0.0, 0.0]), np.array([0.0, 29783.0, 0.0]), 5.972e24)
moon    = Planet('Moon'   , np.array([1.00266*au , 0.0, 0.0]), np.array([0.0, 30805.0, 0.0]), 7.348e22)

给出了更真实的结果（尽管我不确定月球的行为）。