我的重力模拟出了什么问题？

Question

根据此答案中给我的建议，我在重力模拟器中实现了龙格-库塔积分器。

然而，在我模拟太阳系一年后，位置仍然偏离了约110 000公里，这是不可接受的。

我的初始数据是由 NASA 的 HORIZONS 系统提供的。通过它，我获得了行星、冥王星、月球、火卫二和火卫一在特定时间点的位置和速度矢量。

这些矢量是三维的，然而，有些人告诉我，当行星在围绕太阳的板块中排列时，我可以忽略三维，所以我这样做了。我只是将 xy 坐标复制到我的文件中。

这是我改进的更新方法的代码：

"""
Measurement units:

[time] = s
[distance] = m
[mass] = kg
[velocity] = ms^-1
[acceleration] = ms^-2
"""

class Uni:
    def Fg(self, b1, b2):
        """Returns the gravitational force acting between two bodies as a Vector2."""

        a = abs(b1.position.x - b2.position.x) #Distance on the x axis
        b = abs(b1.position.y - b2.position.y) #Distance on the y axis

        r = math.sqrt(a*a + b*b)

        fg = (self.G * b1.m * b2.m) / pow(r, 2)

        return Vector2(a/r * fg, b/r * fg)

    #After this is ran, all bodies have the correct accelerations:
    def updateAccel(self):
        #For every combination of two bodies (b1 and b2) out of all bodies:
        for b1, b2 in combinations(self.bodies.values(), 2):
            fg = self.Fg(b1, b2) #Calculate the gravitational force between them

            #Add this force to the current force vector of the body:
            if b1.position.x > b2.position.x:
                b1.force.x -= fg.x
                b2.force.x += fg.x
            else:
                b1.force.x += fg.x
                b2.force.x -= fg.x


            if b1.position.y > b2.position.y:
                b1.force.y -= fg.y
                b2.force.y += fg.y
            else:
                b1.force.y += fg.y
                b2.force.y -= fg.y

        #For body (b) in all bodies (self.bodies.itervalues()):
        for b in self.bodies.itervalues():
            b.acceleration.x = b.force.x/b.m
            b.acceleration.y = b.force.y/b.m
            b.force.null() #Reset the force as it's not needed anymore.

    def RK4(self, dt, stage):
        #For body (b) in all bodies (self.bodies.itervalues()):
        for b in self.bodies.itervalues():
            rd = b.rk4data #rk4data is an object where the integrator stores its intermediate data

            if stage == 1:
                rd.px[0] = b.position.x
                rd.py[0] = b.position.y
                rd.vx[0] = b.velocity.x
                rd.vy[0] = b.velocity.y
                rd.ax[0] = b.acceleration.x
                rd.ay[0] = b.acceleration.y

            if stage == 2:
                rd.px[1] = rd.px[0] + 0.5*rd.vx[0]*dt
                rd.py[1] = rd.py[0] + 0.5*rd.vy[0]*dt
                rd.vx[1] = rd.vx[0] + 0.5*rd.ax[0]*dt
                rd.vy[1] = rd.vy[0] + 0.5*rd.ay[0]*dt
                rd.ax[1] = b.acceleration.x
                rd.ay[1] = b.acceleration.y

            if stage == 3:
                rd.px[2] = rd.px[0] + 0.5*rd.vx[1]*dt
                rd.py[2] = rd.py[0] + 0.5*rd.vy[1]*dt
                rd.vx[2] = rd.vx[0] + 0.5*rd.ax[1]*dt
                rd.vy[2] = rd.vy[0] + 0.5*rd.ay[1]*dt
                rd.ax[2] = b.acceleration.x
                rd.ay[2] = b.acceleration.y

            if stage == 4:
                rd.px[3] = rd.px[0] + rd.vx[2]*dt
                rd.py[3] = rd.py[0] + rd.vy[2]*dt
                rd.vx[3] = rd.vx[0] + rd.ax[2]*dt
                rd.vy[3] = rd.vy[0] + rd.ay[2]*dt
                rd.ax[3] = b.acceleration.x
                rd.ay[3] = b.acceleration.y

            b.position.x = rd.px[stage-1]
            b.position.y = rd.py[stage-1]

    def update (self, dt):
        """Pushes the uni 'dt' seconds forward in time."""

        #Repeat four times:
        for i in range(1, 5, 1):
            self.updateAccel() #Calculate the current acceleration of all bodies
            self.RK4(dt, i) #ith Runge-Kutta step

        #Set the results of the Runge-Kutta algorithm to the bodies:
        for b in self.bodies.itervalues():
            rd = b.rk4data
            b.position.x = b.rk4data.px[0] + (dt/6.0)*(rd.vx[0] + 2*rd.vx[1] + 2*rd.vx[2] + rd.vx[3]) #original_x + delta_x
            b.position.y = b.rk4data.py[0] + (dt/6.0)*(rd.vy[0] + 2*rd.vy[1] + 2*rd.vy[2] + rd.vy[3])

            b.velocity.x = b.rk4data.vx[0] + (dt/6.0)*(rd.ax[0] + 2*rd.ax[1] + 2*rd.ax[2] + rd.ax[3])
            b.velocity.y = b.rk4data.vy[0] + (dt/6.0)*(rd.ay[0] + 2*rd.ay[1] + 2*rd.ay[2] + rd.ay[3])

        self.time += dt #Internal time variable

算法如下：

1. 更新系统中所有物体的加速度
2. RK4（第一步）
3. 转到1
4. RK4（第二）
5. 转到1
6. RK4（第三）
7. 转到1
8. RK4（第四）

我的 RK4 实施是否搞砸了？或者我只是从损坏的数据开始（重要的实体太少并且忽略了第三维）？

如何解决这个问题？

---

我的数据等的解释...

我的所有坐标都是相对于太阳的（即太阳位于 (0, 0)）。

./my_simulator 1yr

Earth position: (-1.47589927462e+11, 18668756050.4)

HORIZONS (NASA):

Earth position: (-1.474760457316177E+11,  1.900200786726017E+10)

我110 000 km通过从模拟器预测的坐标中减去 NASA 给出的地球 x 坐标得到了错误。

relative error = (my_x_coordinate - nasa_x_coordinate) / nasa_x_coordinate * 100
               = (-1.47589927462e+11 + 1.474760457316177E+11) / -1.474760457316177E+11 * 100
               = 0.077%

相对误差似乎很小，但这只是因为无论是在我的模拟中还是在美国宇航局的模拟中，地球都离太阳很远。距离仍然很大，使我的模拟器毫无用处。

Answer 1

要对太阳系进行 RK4 积分，您需要非常好的精度，否则您的解决方案将很快发散。假设您已正确实现所有内容，您可能会看到 RK 对于此类模拟的缺点。

要验证情况是否如此：尝试不同的积分算法。我发现使用Verlet 积分，太阳系模拟将不会那么混乱。Verlet 的实现也比 RK4 简单得多。

Verlet（和派生方法）在长期预测（如全轨道）方面通常优于 RK4 的原因是它们是辛的，也就是说，动量守恒，而 RK4 则不然。因此，即使在发散之后，Verlet 也会给出更好的行为（真实的模拟，但存在错误），而 RK 一旦发散，就会给出非物理行为。

另外：确保您尽可能地使用浮点。不要使用太阳系中以米为单位的距离，因为浮点数的精度在 0..1 区间内要好得多。使用 AU 或其他标准化比例比使用米好得多。正如其他主题所建议的，请确保使用纪元作为时间，以避免在添加时间步长时累积错误。