我的重力模拟哪出问题了？
根据我在这个回答中得到的建议，我在我的引力模拟器中实现了一个龙格-库塔积分器。
但是，当我模拟了一年的太阳系后，位置仍然偏差大约110,000公里，这个误差是不可接受的。
我最初的数据是通过NASA的HORIZONS系统提供的。通过这个系统，我获得了行星、冥王星、月球、火卫一和火卫二在某个特定时刻的位置和速度向量。
这些向量是三维的，不过，有些人告诉我可以忽略第三维，因为行星围绕太阳排列成一个平面，所以我就这么做了。我只是把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（第一步）
- 返回第1步
- RK4（第二步）
- 返回第1步
- RK4（第三步）
- 返回第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)
我通过从NASA给出的地球x坐标中减去我模拟器预测的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%
相对误差看起来很小，但这只是因为地球离太阳真的很远，无论是在我的模拟中还是在NASA的数据中。这个距离依然很大，使得我的模拟器变得无用。
6 个回答
这种模拟通常不太可靠。因为在计算过程中会出现四舍五入的错误，这些错误会不断累积，导致结果不稳定。即使你提高计算的精度，也帮助不大；问题在于你必须使用一个有限的步长，而自然界实际上是用零步长在运作。
你可以通过减小步长来缓解这个问题，这样错误显现出来的时间就会变长。如果你不是在实时进行模拟，可以使用动态步长，当两个或多个物体非常接近时，步长会自动减小。
我在进行这类模拟时，通常会在每一步之后进行“重新归一化”，也就是调整一下总能量，使其保持不变。整个系统的重力能和动能的总和应该是一个常数（能量守恒）。在每一步之后计算总能量，然后按一个固定的比例调整所有物体的速度，以保持总能量不变。这样至少可以让结果看起来更合理。如果不进行这种调整，每一步后系统会多出或少掉一点能量，导致轨道要么无限膨胀，要么最终掉进太阳里。
要用RK4方法来模拟太阳系，你需要非常高的精度，否则你的结果会很快出错。如果你已经正确实现了所有内容，但还是遇到问题，那可能是RK方法在这种模拟中有缺陷。
要验证是不是这样，可以试试其他的积分算法。我发现使用Verlet积分法来模拟太阳系会更稳定，混乱程度会小很多。而且，Verlet的实现也比RK4简单。
Verlet（以及它的衍生方法）在长期预测（比如完整的轨道）时通常比RK4更好，原因在于它是“辛”的，也就是说它能保持动量，而RK4做不到。因此，即使Verlet在某些情况下出现偏差，它的表现仍然比较合理（虽然有误差），而RK一旦偏差就会出现不符合实际的情况。
另外，确保你尽量使用浮点数。不要在太阳系中用米作为距离单位，因为浮点数在0到1的范围内精度更高。使用天文单位（AU）或其他标准化的尺度比用米要好得多。还有，像其他话题提到的那样，确保使用一个时间基准，避免在添加时间步时累积错误。
110 000 km 的绝对误差意味着什么相对误差？
我通过将我预测的地球x坐标减去NASA的地球x坐标得到了110000公里这个值。
我不太明白你在计算什么，或者你所说的“NASA的地球x坐标”是什么意思。这是相对于哪个起点的距离？使用的是什么坐标系统？在什么时间？（据我所知，地球在绕太阳公转，所以它相对于以太阳为中心的坐标系统的x坐标是一直在变化的。）
无论如何，你通过将你计算的值减去“NASA的地球x坐标”得到了110,000公里的绝对误差。你似乎觉得这个结果不好。你期望的是什么？完全准确？在一米之内？一公里？对你来说什么是可以接受的，为什么？
你可以通过将你的误差差值除以“NASA的地球x坐标”来得到相对误差。可以把它想象成一个百分比。你得到的值是多少？如果是1%或更少，恭喜你，那已经很不错了。
你应该知道，计算机上的浮点数并不精确。（你无法用二进制精确表示0.1，就像你无法用十进制精确表示1/3一样。）会有误差。作为一个模拟器，你的任务是理解这些误差，并尽量将它们减少到最小。
你可能存在步长问题。试着将你的时间步长减半，看看结果是否更好。如果更好，说明你的结果没有收敛。再减半，直到你达到可接受的误差。
你的方程可能条件不良。如果是这样，小的初始误差会随着时间的推移而放大。
我建议你对你的方程进行无量纲化，并计算稳定性极限步长。你对“足够小”的步长的直觉可能会让你感到惊讶。
我还建议你多了解一下多体问题。这个问题比较微妙。
你也可以尝试使用数值积分库，而不是你自己的积分方案。你可以编写方程，然后交给一个工业级的积分器。这可能会让你更清楚是你的实现问题，还是物理问题导致了错误。
就我个人而言，我不太喜欢你的实现。如果你考虑数学向量，解决方案会更好。相对位置的“如果”测试让我觉得不太舒服。使用向量力学会让符号自然出来。
更新：
好的，你的相对误差相当小。
当然，绝对误差是重要的——这取决于你的需求。如果你是在一个星球上着陆，你可不想偏差那么多。
所以你需要停止对什么是过小的步长做假设，做你必须做的事情，把误差控制在可接受的范围内。
你计算中的所有量都是64位IEEE浮点数吗？如果不是，你永远也达不到目标。
一个64位浮点数大约有16位的精度。如果你需要更多，就得使用无限精度的对象，比如Java的BigDecimal，或者——等一下——重新调整你的方程，使用其他单位而不是公里。如果你将所有的距离按对你问题有意义的东西（例如地球的直径或地球轨道的长短轴）进行缩放，可能会得到更好的结果。