爲什麼我的天文模擬不準確？

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，其中時期是一個整數，可以使時間變量的精度與步數無關。