为什么我的天文学模拟不准确？

Question

我已经编写了一个程序来模拟太阳系中物体的运动，但是，我的结果出现了各种不准确之处。

我认为这可能与我的集成方法有关。

---

tl;博士在我的模拟和 NASA 的数据之间，地球的位置和速度之间存在细微差别，如果您可以查看下面的代码并告诉我我的数学是否错误。

---

我运行的测试是一个为期 10 天（864000 秒）的模拟，从 Thu Mar 13 18:30:59 2006 开始，到 Thu Mar 23 18:30:59 2006 结束。

模拟后，程序报告了地球的以下统计数据：

Earth position: (-1.48934630382e+11, -7437423391.22)
Earth velocity: (990.996767368, -29867.6967867)

测量单位当然是米和米每秒。

我已经使用 HORIZONS 系统获得了太阳系中Thu Mar 13 18:30:59 2006 处大多数大型天体的起始位置和速度矢量，并将它们放入模拟中。

测试后，我再次向HORIZONS查询Thu Mar 23 18:30:59 2006Earth数据，得到如下结果：

Earth position: (-1.489348720130393E+11, -7437325664.023257)
Earth velocity: (990.4160633376971, -2986.736541327986)

如您所见，结果的前四位几乎总是相同的。但是，这仍然是一个很大的失误！我很担心，因为我必须模拟几年的时间，而且错误可能会升级。

能否请您看看我的模拟核心并告诉我我的数学是否不正确？

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

    self.time += dt

    for b1, b2 in combinations(self.bodies.values(), 2):
        fg = self.Fg(b1, b2)

        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 b in self.bodies.itervalues():
        ax = b.force.x/b.m
        ay = b.force.y/b.m

        b.position.x += b.velocity.x*dt
        b.position.y += b.velocity.y*dt

        nvx = ax*dt
        nvy = ay*dt

        b.position.x += 0.5*nvx*dt
        b.position.y += 0.5*nvy*dt

        b.velocity.x += nvx
        b.velocity.y += nvy

        b.force.x = 0
        b.force.y = 0

我有这个方法的另一个版本，它应该表现更好，但它表现得更差：

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

    self.time += dt

    for b1, b2 in combinations(self.bodies.values(), 2):
        fg = self.Fg(b1, b2)

        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 b in self.bodies.itervalues():
        #Acceleration at (t):
        ax  = b.force.x/b.m
        ay  = b.force.y/b.m
        #Velocity at (t):
        ovx = b.velocity.x
        ovy = b.velocity.y
        #Velocity at (t+dt):
        nvx = ovx + ax*dt
        nvy = ovy + ay*dt
        #Position at (t+dt):
        b.position.x = b.position.x + dt*(ovx+nvx)/2
        b.position.y = b.position.y + dt*(ovy+nvy)/2


        b.force.null() #Reset the forces.

Answer 1

集成方法非常很重要。您正在使用欧拉显式方法，该方法的阶精度低，对于适当的物理模拟来说太低了。现在，您可以选择

• 一般行为最重要：Verlet method 或 Beeman method（精度更高的 Verlet），它们具有非常好的能量守恒，但位置和速度的精度较低。
• 精确的位置最重要：Runge-Kutta 订单 4 或更多。能量不会守恒，因此您的模拟系统会表现得好像能量增加了一样。

此外，对于大量步骤，时间 = 时间 + dt 会增加精度损失。考虑 time = epoch * dt 其中 epoch 是一个整数，会使时间变量的精度与步数无关。