Python 重力模拟器表现奇怪

Question

我正在用 Python 进行重力模拟(准确地说是用 VPython 进行 3D 模拟)，我确信代码没有任何问题，但当两个物体彼此靠近时，它的行为会很奇怪。

我的灵感是http://testtubegames.com/gravity.html 。注意如何放置两颗没有速度的行星，它们相互移动、超越、减速和返回。在我的程序中，它们超车并减速，但仅与距离成正比，因此从技术上讲，它无论如何都不应该回头。

我意识到，如果 r(距离)太接近 0，则 f=G*(m1*m2)/r**2 定律将不起作用，因此我添加了最大输出，因此如果它小于 1，则设置为 1(顺便说一句，单位不是像素)，但它仍然不起作用。

简单的逻辑还表明对象不应该以这种方式使用react，所以接下来的事情是我一定错过了一些东西。

以下是代码摘录:

from visual import *
a = sphere(x=-10,mass=10, vel=vector())
b = sphere(x=10, mass=10, vel=vector())

while 1:
    rate(20)

    #distance between the two objects, a and b, where a.r.mag would be the magnitude of the vector
    a.r = b.pos - a.pos
    b.r = a.pos - b.pos

    a.force = a.r
    if a.r.mag > 1:
        a.force.mag = (a.mass * b.mass) / a.r.mag**2
    else:
        a.force.mag = (a.mass * b.mass) / 1
    a.vel = a.vel + a.force / a.mass


    b.force = b.r
    if b.r.mag > 1:
        b.force.mag = (a.mass * b.mass) / b.r.mag**2
    else:
        b.force.mag = (a.mass * b.mass) / 1
    b.vel = b.vel + b.force / b.mass

    a.pos = a.pos + a.vel
    b.pos = b.pos + b.vel

编辑:为响应 shockburner 重新编写的代码:

from visual import *
import sys

limit2 = sys.float_info.min
limit = limit2**0.5
timestep = 0.0005

a = sphere(x=-5,mass=10, vel=vector())
b = sphere(x=5, mass=10, vel=vector())

def force(ob1, ob2):
    ob1.r = ob2.pos - ob1.pos
    ob1.force = ob1.r + vector()
    if ob1.r.mag > limit:
        ob1.force.mag = (ob1.mass * ob2.mass) / ob1.r.mag2
    else:
        ob1.force.mag = (ob1.mass * ob2.mass) / limit2
    return ob1.force

while 1:
    rt = int(1/timestep)
    rate(rt)

    a.acc = force(a, b) / a.mass
    b.acc = force(b, a) / b.mass

    a.pos = a.pos + timestep * (a.vel + timestep * a.acc / 2)
    b.pos = b.pos + timestep * (b.vel + timestep * b.acc / 2)

    a.acc1 = force(a,b) / a.mass
    b.acc1 = force(b,a) / b.mass

    a.vel = a.vel + timestep * (a.acc + a.acc1) / 2
    b.vel = b.vel + timestep * (b.acc + b.acc1) / 2

任何正确方向的帮助或指示将不胜感激，如果答案非常简单(对我来说通常是这种情况)，请记住我无论如何都是个白痴。

Answer 1

我的猜测是您的问题源于积分方法中的数值错误。看来您使用的是欧拉方法，该方法容易出现较大的数值误差，因为它是一阶积分方法。我会推荐velocity verlet用于轨道数值积分，因为它是一种二阶方法，还可以保留总能量(动能 + 引力势)以保证机器精度。这种能量守恒通常使速度 Verlet 比四阶更稳定 Runge–Kutta ，因为束缚轨道保持束缚状态。

此外，您可能需要考虑使用动态时间步长而不是静态时间步长。当粒子靠近在一起时，速度和位置变化得更快。因此，为了减少数值误差，您需要采取较小的时间步长。

最后，我会让你的限制器(if a.r.mag > 1:)尽可能小/实用。我会尝试以下操作:

import sys
limit2 = sys.float_info.min
limit = limit2**.5
...
if a.r.mag > limit:
    a.force.mag = (a.mass * b.mass) / a.r.mag**2
else:
    a.force.mag = (a.mass * b.mass) / limit2
...