Согласно совету, данному мне в this answer, я применил интегратор Runge-Kutta в моем симуляторе силы тяжести.Что не так с моим симуляцией гравитации?
Однако после того, как я смоделировал один год солнечной системы, позиции по-прежнему отключены примерно на 110 000 километров, что неприемлемо.
Мои исходные данные были предоставлены системой ГОРИЗОНТОВ НАСА. Через него я получил положения и скорости векторов планет, Плутона, Луны, Деймоса и Фобоса в определенный момент времени.
Эти векторы были 3D, однако некоторые люди сказали мне, что я мог бы игнорировать третье измерение, поскольку планеты выровнялись в пластине вокруг Солнца, и поэтому я это сделал. Я просто скопировал координаты x-y в свои файлы.
Это код моего усовершенствованного способа обновления:
"""
Measurement units:
[time] = s
[distance] = m
[mass] = kg
[velocity] = ms^-1
[acceleration] = ms^-2
"""
class Uni:
def Fg(self, b1, b2):
"""Returns the gravitational force acting between two bodies as a Vector2."""
a = abs(b1.position.x - b2.position.x) #Distance on the x axis
b = abs(b1.position.y - b2.position.y) #Distance on the y axis
r = math.sqrt(a*a + b*b)
fg = (self.G * b1.m * b2.m)/pow(r, 2)
return Vector2(a/r * fg, b/r * fg)
#After this is ran, all bodies have the correct accelerations:
def updateAccel(self):
#For every combination of two bodies (b1 and b2) out of all bodies:
for b1, b2 in combinations(self.bodies.values(), 2):
fg = self.Fg(b1, b2) #Calculate the gravitational force between them
#Add this force to the current force vector of the body:
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 body (b) in all bodies (self.bodies.itervalues()):
for b in self.bodies.itervalues():
b.acceleration.x = b.force.x/b.m
b.acceleration.y = b.force.y/b.m
b.force.null() #Reset the force as it's not needed anymore.
def RK4(self, dt, stage):
#For body (b) in all bodies (self.bodies.itervalues()):
for b in self.bodies.itervalues():
rd = b.rk4data #rk4data is an object where the integrator stores its intermediate data
if stage == 1:
rd.px[0] = b.position.x
rd.py[0] = b.position.y
rd.vx[0] = b.velocity.x
rd.vy[0] = b.velocity.y
rd.ax[0] = b.acceleration.x
rd.ay[0] = b.acceleration.y
if stage == 2:
rd.px[1] = rd.px[0] + 0.5*rd.vx[0]*dt
rd.py[1] = rd.py[0] + 0.5*rd.vy[0]*dt
rd.vx[1] = rd.vx[0] + 0.5*rd.ax[0]*dt
rd.vy[1] = rd.vy[0] + 0.5*rd.ay[0]*dt
rd.ax[1] = b.acceleration.x
rd.ay[1] = b.acceleration.y
if stage == 3:
rd.px[2] = rd.px[0] + 0.5*rd.vx[1]*dt
rd.py[2] = rd.py[0] + 0.5*rd.vy[1]*dt
rd.vx[2] = rd.vx[0] + 0.5*rd.ax[1]*dt
rd.vy[2] = rd.vy[0] + 0.5*rd.ay[1]*dt
rd.ax[2] = b.acceleration.x
rd.ay[2] = b.acceleration.y
if stage == 4:
rd.px[3] = rd.px[0] + rd.vx[2]*dt
rd.py[3] = rd.py[0] + rd.vy[2]*dt
rd.vx[3] = rd.vx[0] + rd.ax[2]*dt
rd.vy[3] = rd.vy[0] + rd.ay[2]*dt
rd.ax[3] = b.acceleration.x
rd.ay[3] = b.acceleration.y
b.position.x = rd.px[stage-1]
b.position.y = rd.py[stage-1]
def update (self, dt):
"""Pushes the uni 'dt' seconds forward in time."""
#Repeat four times:
for i in range(1, 5, 1):
self.updateAccel() #Calculate the current acceleration of all bodies
self.RK4(dt, i) #ith Runge-Kutta step
#Set the results of the Runge-Kutta algorithm to the bodies:
for b in self.bodies.itervalues():
rd = b.rk4data
b.position.x = b.rk4data.px[0] + (dt/6.0)*(rd.vx[0] + 2*rd.vx[1] + 2*rd.vx[2] + rd.vx[3]) #original_x + delta_x
b.position.y = b.rk4data.py[0] + (dt/6.0)*(rd.vy[0] + 2*rd.vy[1] + 2*rd.vy[2] + rd.vy[3])
b.velocity.x = b.rk4data.vx[0] + (dt/6.0)*(rd.ax[0] + 2*rd.ax[1] + 2*rd.ax[2] + rd.ax[3])
b.velocity.y = b.rk4data.vy[0] + (dt/6.0)*(rd.ay[0] + 2*rd.ay[1] + 2*rd.ay[2] + rd.ay[3])
self.time += dt #Internal time variable
алгоритм выглядит следующим образом:
- Update ускорения всех тел в системе
- RK4 (первый этап)
- goto 1
- RK4 (второй)
- Гото 1
- RK4 (третий)
- Гото 1
- RK4 (четвёртое)
ли я что-то запутались с моей реализации RK4? Или я только начал с испорченных данных (слишком мало важных органов и игнорировало 3-е измерение)?
Как это можно исправить?
Объяснение моих данных и т.д. ...
Все мои координаты относительно Солнца (то есть Солнце находится в точке (0, 0)).
./my_simulator 1yr
Earth position: (-1.47589927462e+11, 18668756050.4)
HORIZONS (NASA):
Earth position: (-1.474760457316177E+11, 1.900200786726017E+10)
я получил ошибку 110 000 km путем вычитания х Земли координат задается НАСА из одного предсказанного моего тренажере.
relative error = (my_x_coordinate - nasa_x_coordinate)/nasa_x_coordinate * 100
= (-1.47589927462e+11 + 1.474760457316177E+11)/-1.474760457316177E+11 * 100
= 0.077%
Относительная ошибка кажется мизерной, но это просто потому, что Земля находится очень далеко от Солнца, как в моем моделировании и в НАСА. Расстояние все еще огромный и делает мой симулятор бесполезным.
Это не вопрос программирования, это проблема математики (от некоторых ва riety). –
Хорошая мысль, но я считаю, что она хорошо подходит для SO, потому что речь идет о конкретной реализации алгоритма - это то, чего вы не делаете в математике, но в основном в программировании. Я также задал очень похожий вопрос на этот (но принципиально иной), и он был очень хорошо принят. – corazza
Хорошо, но вы по сути просто просите людей обнаружить ошибки в вашем коде путем проверки, что не является большим вопросом. –