1

I'm making a gravity simulation in Python (in 3D with VPython, to be exact) and I'm sure there's nothing wrong with the code, but it behaves strangely when two objects get close to each other.

My inspiration is http://testtubegames.com/gravity.html. Note how you can place two planets with no velocity, they move towards each other, overtake, decelerate and turn back. In my program, they overtake, and decelerate, but only proportionately to the distance, so technically it should never turn back anyway.

I realize that the law f=G*(m1*m2)/r**2 wouldn't work if r (the distance) is or gets too close to 0, so I have included a max-out, so if it is less than 1 it is set to 1 (units not in pixels, by the way), but it still does not work.

Simple logic also suggests that the objects should not react in this way, so the next thing that follows is that I must be missing something.

Here is an extract of the code:

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

EDIT: Code re-written in response to 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

Any help or pointer in the right direction would be greatly appreciated, and if the answer turns out to be idiotically simple (which with me is usually the case) remember I am quite an idiot anyway.

4

2 回答 2

2

我的猜测是您的问题源于您的积分方法中的数值错误。您似乎正在使用欧拉方法，因为它是一阶积分方法，因此容易产生较大的数值误差。我建议对轨道进行数值积分的速度 verlet，因为它是一种二阶方法，也可以将总能量（动能 + 引力势）保留到机器精度。这种能量守恒通常使速度 verlet 比 4 阶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
...
于 2013-09-05T23:20:28.207 回答
0

我以前也遇到过这个问题。如果你直接去龙格库塔，一切都会自行解决。此 pdf 将解释如何合并该方法：http ://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf 。祝你好运！

于 2013-09-25T14:52:57.840 回答