我的重力模拟有什么问题？

Question

根据this answer中给我的建议，我在重力模拟器中实现了一个龙格库塔积分器。在

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

我最初的数据是由美国宇航局的地平线系统提供的。通过它，我获得了行星，冥王星，月亮，戴莫斯和火卫一在特定时间点的位置和速度矢量。在

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

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

"""
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））。在

^{pr2}$

我得到了110 000 km错误，从我的模拟器预测的x坐标中减去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

你需要很快对太阳系进行精确的积分。假设您已经正确地实现了所有的东西，您可能会看到这种模拟的RK的缺点。在

若要验证是否存在这种情况，请尝试不同的积分算法。我发现使用Verlet integration模拟太阳系将不会太混乱。Verlet的实现也比RK4简单得多。在

Verlet（和导出的方法）在长期预测（如全轨道）方面通常比RK4好的原因是它们是辛的，也就是说，保持了RK4没有的动量。因此，Verlet将给出一个更好的行为，即使在它发散（一个真实的模拟，但有一个错误），而RK将给出一个非物理行为一旦它发散。在

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

Answer 2

110 000 km绝对误差是什么相对误差？在

I got the 110 000 km value by subtracting my predicted Earth's x coordinate with NASA's Earth x coordinate.

我不知道你在计算什么或者你说的“美国宇航局的地球x坐标”是什么意思。这是从哪个原点，在什么坐标系，什么时间的距离？（据我所知，地球是绕太阳运行的，所以它的x坐标w.r.t.是以太阳为中心的坐标系一直在变化。）

在任何情况下，你计算出的绝对误差是11万公里，从“美国宇航局的地球x坐标”中减去你的计算值。你似乎认为这是个不好的回答。你的期望是什么？直接击中目标？在一米之内？一公里？你能接受什么？为什么？在

用误差差除以“美国宇航局的地球x坐标”得到相对误差。把它想象成一个百分比。你有什么价值？如果是1%或更少，恭喜你自己。那就太好了。在

你应该知道floating point numbers aren't exact on computers。（你不能精确地用二进制来表示0.1，正如你不能精确地用十进制表示1/3一样。）会有错误。作为一个模拟器，你的工作就是理解错误并尽可能减少它们。在

你可能会有一个台阶大小的问题。试着把你的时间步长减少一半，看看你是否做得更好。如果你这样做了，它说明你的结果还没有收敛。再减少一半，直到达到可接受的误差。在

你的方程式可能条件不好。如果这是真的，小的初始误差会随着时间的推移而放大。在

我建议你把方程无量纲化，计算稳定极限步长。你对“足够小”步长的直觉可能会让你大吃一惊。在

我还读了更多关于many body problem的文章。很微妙。在

您也可以尝试使用数值积分库来代替积分方案。你将编程你的方程，并把它们交给一个工业强度积分器。它可以让你了解到底是你的实现还是物理造成了问题。在

我个人不喜欢你的实现。如果你考虑到数学向量的话，这会是一个更好的解决方案。相对位置的“如果”测试让我很冷淡。矢量力学可以让符号自然地显现出来。在

更新：

好吧，你的相对误差很小。在

当然，绝对误差是很重要的-取决于您的要求。如果你要在一个星球上着陆，你不想离开那么远。在

因此，您需要停止对什么构成了太小的步骤大小的假设，并做您必须做的，以推动错误到可接受的水平。在

计算中的所有数量是64位IEEE浮点数吗？否则，你永远也到不了那里。在

一个64位浮点数大约有16 digits of accuracy。如果你需要更多，你必须使用一个无限精度的对象，比如Java的BigDecimal，或者——等等——重新缩放你的方程，使其使用千米以外的长度单位。如果你用对你的问题有意义的东西来衡量你所有的距离（例如，地球的直径或地球轨道长轴/短轴的长度），你可能会做得更好。在

Answer 3

这种模拟是出了名的不可靠。舍入误差累积并引入不稳定性。提高精度没有多大帮助；问题是你（必须）使用有限的步长，而大自然使用的是零步长。在

您可以通过减小步长来减少问题，因此需要更长的时间才能使错误变得明显。如果不是实时执行此操作，则可以使用动态步长，如果两个或多个实体非常接近，则可以减小该步长。在

我对这些模拟所做的一件事就是在每一步之后“重新标准化”，使总能量相同。整个系统的重力和动能之和应该是一个常数（能量守恒）。计算出每一步后的总能量，然后将所有物体的速度按一定的比例缩放以保持总能量不变。这至少让输出看起来更合理。如果没有这种缩放，在每一步之后都会有少量的能量被添加到系统中，或者从系统中移除，而轨道往往会膨胀到无穷大，或者盘旋进入太阳。在