根据这个答案给我的建议，我已经在我的重力模拟器中实现了一个龙格-库塔积分器。
然而，在我模拟了一年的太阳系之后，位置仍然偏离了大约11万公里，这是不可接受的。
我最初的数据是由美国宇航局的地平线系统提供的。通过它，我获得了行星、冥王星、月球、火卫一和火卫一在特定时间点的位置和速度矢量。
这些向量是三维的，然而，有些人告诉我，我可以忽略第三维度，因为行星在太阳周围排成一个板块，我就这样做了。我只是把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
算法如下:
- 更新系统中所有物体的加速度
- RK4(第一步)
- goto 1
- RK4(二)
- goto 1
- RK4(三)
- goto 1
- RK4(四)
我搞砸了我的RK4实现吗?或者我只是从损坏的数据开始(重要的主体太少，忽略了第三维度)?
如何解决这个问题?
解释我的数据等…
我所有的坐标都相对于太阳(即太阳在(0,0)处)
./my_simulator 1yr
Earth position: (-1.47589927462e+11, 18668756050.4)
HORIZONS (NASA):
Earth position: (-1.474760457316177E+11, 1.900200786726017E+10)
我用模拟器预测的地球x坐标减去NASA给出的地球x坐标，得到了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%
相对误差似乎很小，但这只是因为在我的模拟和美国宇航局的模拟中，地球离太阳都很远。距离仍然很大，使我的模拟器无用。
110 000 km绝对误差指的是什么相对误差?
我通过减去我预测的地球x得到了110,000 km的值与NASA的地球x坐标。
我不知道你在这里计算什么，或者你说的"NASA的地球x坐标"是什么意思。这是从哪个原点，在哪个坐标系，在什么时间的距离?(据我所知，地球是绕太阳运行的，所以它的x坐标w.r.t.——以太阳为中心的坐标系——一直在变化。)
无论如何，通过从"NASA的地球x坐标"中减去计算值，您计算出的绝对误差为110,000 km。你似乎觉得这是个糟糕的答案。你的期望是什么?要击中目标吗?在一米之内?一公里?什么是你可以接受的，为什么?
通过将误差差除以"NASA的地球x坐标"得到相对误差。把它想象成一个百分比。你得到了什么值?如果只有1%或更少，恭喜你自己。那太好了。
你应该知道浮点数在计算机上是不精确的。(你不能精确地用二进制表示0.1，就像你不能精确地用十进制表示1/3一样。)会有误差。作为模拟器，你的工作是了解错误并尽可能地减少它们。
你可能有步长问题。试着把你的时间步长减少一半，看看你是否做得更好。如果你这样做了，它说你的结果没有收敛。再减少一半，直到达到可接受的误差。
你的方程可能条件不好。如果这是真的，那么小的初始错误将随着时间的推移而被放大。
我建议你对你的方程进行无量纲化，并计算稳定性极限步长。你对"足够小"步长应该是多少的直觉可能会让你大吃一惊。
我还读了一些关于身体问题的书。这是微妙的。
您也可以尝试使用数值积分库来代替您的积分方案。你要把你的方程式编程，然后交给工业强度的积分器。它可以让你了解到是你的实现还是物理导致了问题。
就我个人而言，我不喜欢你的实现。如果你用数学向量来解的话会更好。相对位置的"如果"测试让我感到困惑。矢量机制会让符号自然地出现。 更新:好的，你的相对误差相当小。
当然绝对误差很重要——这取决于你的需求。如果你要让飞行器降落在一个星球上你不想偏离那么多。
因此，您需要停止假设是什么构成了过小的步长，并采取必要措施将误差驱动到可接受的水平。
你计算的所有数量都是64位IEEE浮点数吗?如果你不这样做，你永远也不会成功。
64位浮点数的精度约为16位。如果需要更多，就必须使用像Java的BigDecimal这样的无限精度对象，或者——等等——重新缩放方程，使用千米以外的长度单位。如果你用对你的问题有意义的东西来衡量你所有的距离(例如，地球的直径或地球轨道的长/短轴的长度)，你可能会做得更好。
要对太阳系进行RK4积分，你需要非常高的精度，否则你的解会很快发散。假设您已经正确地实现了所有内容，您可能会看到RK用于这种模拟的缺点。
要验证这种情况:尝试不同的积分算法。我发现，使用Verlet积分，一个太阳系的模拟将不那么混乱。Verlet也比RK4更容易实现。
对于长期预测(如完整轨道)，Verlet(及其衍生方法)通常比RK4更好的原因是它们是辛的，也就是说，RK4没有保留动量。因此，即使在它发散之后，Verlet也会给出更好的行为(一个真实的模拟，但其中有一个错误)，而RK一旦发散就会给出非物理行为。
还有:确保你尽可能好地使用浮点数。在太阳系中不要使用米为单位的距离，因为浮点数的精度在0..1区间。使用单位或其他标准化的刻度比使用米要好得多。正如另一个主题所建议的，确保使用epoch作为时间，以避免在添加时间步长时累积错误。
这样的模拟是出了名的不可靠。舍入误差累积并引入不稳定性。提高精确度并没有多大帮助;问题是你(必须)使用有限的步长，而大自然使用零步长。
您可以通过减小步长来减少问题，这样错误就需要更长的时间才能变得明显。如果你不是实时地做这些，你可以使用动态步长，当两个或更多的物体非常接近时，它会减少。
我对这类模拟所做的一件事是在每一步之后"重新正常化"，使总能量相同。整个系统的引力和动能之和应该是一个常数(能量守恒)。计算出每一步后的总能量，然后将所有物体的速度乘以一个常数，以保持总能量不变。这至少使输出看起来更可信。如果没有这种缩放，每一步都会有少量的能量被添加到系统中或从系统中移除，并且轨道往往会膨胀到无穷大或螺旋状进入太阳。非常简单的改变将会改善事情(正确使用浮点值)
- 改变单位系统，使用尽可能多的尾数位。使用仪表，你做错了……如上所述，使用AU。更棒的是，按比例缩放，这样太阳系就能放进一个1x1x1的盒子里了
- 已经在另一篇文章中说过:你的时间，计算它为time = epoch_count * time_step，而不是通过添加!这样可以避免累积错误。
- 当对几个值进行求和时，使用高精度和算法，如Kahan求和。在python中，是数学。fsum为您做。
力的分解不应该是
th = atan(b, a);
return Vector2(cos(th) * fg, sin(th) * fg)
(参见http://www.physicsclassroom.com/class/vectors/Lesson-3/Resolution-of-Forces或https://fiftyexamples.readthedocs.org/en/latest/gravity.html#solution)
顺便说一句:你取平方根来计算距离，但你实际上需要的是距离的平方…
为什么不简化
r = math.sqrt(a*a + b*b)
fg = (self.G * b1.m * b2.m) / pow(r, 2)
r2 = a*a + b*b
fg = (self.G * b1.m * b2.m) / r2
我不确定python，但在某些情况下，你可以得到更精确的计算中间结果(英特尔cpu支持80位浮点数，但当分配给变量时，它们被截断为64位):如果存储在中间"double"变量
你的行星坐标和速度在什么参照系中并不十分清楚。如果它是在日心参考系(参考系系在太阳上)，那么你有两个选择:(i)在非惯性参考系(太阳不是静止的)中进行计算，(ii)将位置和速度转换为惯性质心参考系。如果你的坐标和速度是在质心参考系中，那么你也必须有太阳的坐标和速度。