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

Question

我制作了一个模拟太阳系中身体运动的程序，但是，我的结果中出现了各种不准确的地方。

我认为这可能与我的集成方法有关。如果你可以看看我的代码，并告诉我我的数学是否是错误的，那么我的模拟和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 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 method，或Beeman method（verlet的更精确），它具有很好的节能，但精度较低的位置和速度。
• 精确位置最重要：Runge-Kutta订购4个或更多。能量不会被节约，所以你的模拟系统会像能量增加一样起作用。

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