2D N体模拟

Question

我按照维基百科上发现的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（虽然没有理由它根本不应该，所以我通过将rand尺寸设置为（100,3）来尝试3D，但行为是相同的。

我查看了其他在线可用的代码，但我似乎无法找到我做错了（或与其他人不同），所以任何帮助都将不胜感激。

---

编辑1这实际上似乎与方程式一致。我已经为[-1,1]和[1,1]（忽略G）手动设计了前几个步骤，对于p1，力分别为[0.25,0.7,81,0,0]。然而，由于从第三步开始速度很快，并且粒子p2与p1相反，它们移动的速度非常快。但是，在线轻松找到的其他实现不会遇到此问题。我似乎无法弄清楚为什么。我认为它可能是初始化，但其他实现似乎并没有受此影响。

Answer 1

我的dt太大了。将dt设置为较小的值，例如0.05做到了。