【发布时间】：2018-08-08 20:33:31
【问题描述】：
我遵循了维基百科上的 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 的相反，它们移动得非常快。但是，在网上很容易找到的其他实现不会遇到这个问题。我似乎无法弄清楚为什么。我认为这可能是初始化，但其他实现似乎没有受到此影响。
【问题讨论】：
-
你应该做一些调试。单步执行一个非常简单的示例（例如两个主体）的代码，以查看行为与预期不同的地方。
-
我已经做了相当多的调试，实际上这是我的调试代码。对于 2 个物体的最简单情况，它会发散，它们显然被排斥（对于力矢量的任一方向）。我的猜测是我一定遗漏了更新规则的一些关键部分，但我似乎找不到我遗漏的内容或写错了。
-
我的建议是，对于两个物体，你应该能够逐帧推断行为（你甚至可以手动计算）——一切都移动了吗例如，第一帧中的正确方向？
-
我现在已经手动完成了几个步骤，似乎对于最简单的情况 [-1, 1], [1, 1]，p1 上的力 = [0.25, 0.7 , 81, 0, 0]。所以这实际上是预期的行为？然而，速度仍然非常高，所以它一直在移动。另一个粒子则相反，因此它们彼此远离。这就引出了一个问题，为什么我们可以在网上找到的其他实现没有这个问题？我认为没有原则性的方法可以解决这个问题，但其他人似乎很容易做到。
-
也许你的时间步长太大了？对于从零初始条件开始的两个物体，我相信预期的行为是它们无限期地振荡（基本上是钟摆）。如果你的时间步长大于振荡周期，事情就会出错。
标签： python algorithm simulation physics