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

Question

我制作了一个模拟太阳系中物体运动的程序,然而,我的结果却出现了各种不准确之处.

我相信它可能与我的集成方法有关.

---

tl; dr我的模拟和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 2006地球数据,得到了以下结果:

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

整合方法非常重要.您正在使用Euler显式方法,该方法具有低阶精度,对于正确的物理模拟而言太低.现在,你可以做出选择

• 一般行为最重要:Verlet方法,或Beeman方法(Verlet 方法更精确),它具有非常好的节能但位置和速度精度较低.
• 精确位置最重要:Runge-Kutta订单4或更多.能量将无法保存,因此您的模拟系统将表现得像能量增加.

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