为什么我的天文模拟不准确？ [英] Why is my astronomy simulation inaccurate?
问题描述
我相信可能与我的整合方法有关。
>
tl; dr在我的模拟之间，地球的位置和速度有一点点差异和美国航空航天局的数据，如果你可以请看看我的代码，并告诉我，我的数学是否是错误的。
测试我跑了10天（864000秒）的模拟，从 Thu Mar 13 18:30:59 2006 开始，结束于星期四三月二十三日18:30:59 2006 。
模拟结束后，该计划报告了地球的以下统计数据：
地球的位置是：（-1.48934630382e + 11，-7437423391.22）
地球速度：（990.996767368，-29867.6967867）
测量单位当然是米和米第二。
我已经使用HORIZONS系统获取大多数大型机构的起始位置和速度矢量，位于 c> Thu Mar 23 18:30:59 2006 Thu Mar 13 18:30：在测试之后，我再次询问HORIZONS是为了$ 在太阳系中，并将它们放在模拟中。 Earth数据，得到如下结果：
地球的位置是：（-1.489348720130393E + 11，-7437325664.023257）
地球速度：（990.4160633376971，-2986.736541327986）
正如您所看到的，结果在前四位数字中几乎总是相同的。但是，这仍然是一个相当大的错过！我很担心，因为我将不得不模拟几年的时间，错误可能会升级。
请您看看我的模拟的核心并告诉我是否我的数学不正确？
$ b $ pre $ def update（self，dt）：
推uni'dt ''
self.time + = dt
为b1，b2组合（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 self.bodies.itervalues（）：
ax = b.force .x / bm
ay = b.force.y / bm
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
我有这个方法的另一个版本，应该执行得更好，但是性能要差得多：
$ b $ pre $ def update（self，dt）：
$ .bodies.valu es（），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 self.bodies.itervalues（）：
#Acceleration at（t ）：
ax = b.force.x / bm
ay = b.force.y / bm
#Velocity at（t）：
ovx = b.velocity.x
ovy = b.velocity.y
#Velocity at（t + dt）：
nvx = ovx + ax * dt
nvy = ovy + ay * dt
＃位置在（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（）＃重置部队。
>重要。您正在使用欧拉显式方法，该方法具有低阶精度，对于正确的物理模拟来说过低。现在，您可以选择
另外，time = time + dt将会大量步骤的精度损失增加。考虑时间= epoch * dt，其中epoch是一个整数，可以使时间变量的精度与步数无关。
I've made a program that simulates movement of bodies in the solar system, however, I'm getting various inaccuracies in my results.
I believe that it probably has something to do with my integration method.
tl;dr there's a slight difference between the position and velocity of Earth between my simulation and NASA's data, if you could please look at my code below and tell me whether my math is wrong.
The test I've ran is a 10 day long (864000 seconds) simulation, that starts at Thu Mar 13 18:30:59 2006 and ends at Thu Mar 23 18:30:59 2006.
After the simulation the program reported the following statistics for the planet Earth:
Earth position: (-1.48934630382e+11, -7437423391.22)
Earth velocity: (990.996767368, -29867.6967867)
The measurement units are, of course meters and meters per second.
I've used the HORIZONS system to get the starting position and velocity vectors of most large bodies at Thu Mar 13 18:30:59 2006 in the solar system and put them in the simulation.
After the test, I've queried HORIZONS again for Thu Mar 23 18:30:59 2006 Earth data, and got the following results for it:
Earth position: (-1.489348720130393E+11, -7437325664.023257)
Earth velocity: (990.4160633376971, -2986.736541327986)
As you can see, the results are almost always the same in their first four digits. However, that is still a pretty big miss! I'm worried because I'll have to simulate a few years of time and the error could escalate.
Could you please look at the core of my simulation and tell me whether my maths are incorrect?
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
I have another version of this method, that should perform better, but it performs much worse:
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.
The integration method is very important. You are using Euler explicit method, which is of low order precision, too low for proper physics simulation. Now, you get choices
- General behaviour matters most : Verlet method, or Beeman method (Verlet with more precision), which have very good energy conservation but lesser precision for position and speed.
- Precise positions matters most : Runge-Kutta order 4 or more. Energy won't be conserved, so your simulated system will act as if energy increases.
Furthermore, time = time + dt will have increasing loss of precisions for large number of steps. Consider time = epoch * dt where epoch is an integer, would make the precision of the time variable independent from the number of steps.
这篇关于为什么我的天文模拟不准确？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
