二维N实体模拟

我遵循维基百科上n体问题的方程式,实现了一个简单的O(n)n体模拟。然而,一旦我将模拟可视化,事情就不会像预期的那样进行,即所有粒子都会远离中心,好像它们具有很高的排斥力。起初我想我可能弄错了力向量的方向,但我试着翻转它,结果几乎是一样的。

data = np.random.rand(100, 2)

velocities = np.zeros_like(data)
masses = np.ones_like(data)
dt = 60 * 60 * 24


for _ in range(10000):
    forces = np.zeros_like(data)

    for i, node1 in enumerate(data):
        for j, node2 in enumerate(data):
            d = node2 - node1
            # First term is gravitational constant, 1e-8 is a softening factor
            forces[i] += 6.67384e-11 * d / (np.sqrt(d.dot(d) + 1e-8) ** 3)

    velocities += forces * dt / masses
    data += velocities * dt

    yield data  # for visualization

我还认为它可能在2D中不起作用(虽然没有理由不起作用,所以我也在3D中尝试过,将rand dimensions设置为(100,3),但行为是一样的。

我已经查看了其他在线可用的代码,但我似乎找不到我做错了什么(或与其他人的不同),所以如果有任何帮助,我将不胜感激。

编辑1 这实际上似乎与方程式一致。我已经手工计算了[-1,1]和[1,1](忽略G)的前几步,对于p1,力分别为[0.25,0.7,81,0,0]。然而,由于第三步的速度很高,粒子p2的速度与p1相反,所以它们移动得非常快。然而,其他在线容易找到的实现并没有面临这个问题。我似乎不明白为什么。我认为这可能是初始化,但其他实现似乎没有受到影响。