太陽系の物体の動きをシミュレートするプログラムを作成しましたが、結果にさまざまな不正確さが生じています。
それはおそらく私の統合方法と関係があると思います。
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)
ご覧のとおり、結果はほとんど常に最初の 4 桁で同じです。しかし、それはまだかなり大きなミスです！数年後にシミュレートする必要があり、エラーがエスカレートする可能性があるため、心配しています。
私のシミュレーションの核心を見て、私の計算が間違っているか教えていただけませんか?
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.