引力 N 体模拟的问题

Question

背景 我正在用 Python 编写一个简单的 N 体模拟代码，核心物理求解器，例如在 Cython 中实现的加速度和位置集成。代码使用了蛙跳法来整合位置。

条件

1. 质量：10，100 用于 2 个主体； 1 表示所有机构，机构数量 > 2。
2. 位置：随机生成
3. 速度：随机生成
4. 时间步长：1e-3

错误

1. 如果物体数量 > 2：物体相互排斥而不是相互吸引。
2. 对于物体数量 = 2：它们不相互绕行，而是沿直线运动。
3. 一般错误（不考虑物体数量）：排斥力

预期行为

1. 两个物体必须相互绕行
2. 这些力量必须具有吸引力

解决问题的努力

1. 加速的负号
2. 将加速度乘以 -1
3. 添加新的临时表达式

代码 加速函数（Cython）：

def acceleration(np.ndarray[np.float64_t, ndim=2] pos, np.ndarray[np.float64_t, ndim=1] mass):
cdef int N = pos.shape[0] # Number of bodies   
# Memoryview for acceleration
cdef np.float64_t [:,:] acc = np.zeros((N,3),dtype="float64")
cdef double soft = 1e-4 # Softening length
cdef G = 6.673e-11 # Gravitational constant
# Pairwise separations
cdef double dx
cdef double dy
cdef double dz
# Total separation vectors
cdef double r
cdef double tmp
# Mass
cdef mj
# Acceleration calculation loop
for i in range(N):
    for j in range(N):
        # Remove gravitational self forces
        if i==j:
            continue
        
        # Calculate pairwise separation vectors
        dx = pos[j,0] - pos[i,0]
        dy = pos[j,1] - pos[i,1]
        dz = pos[j,2] - pos[i,2]

        # Vector magnitude of separation vector
        r = dx**2 + dy**2 + dz**2
        r = np.sqrt(r)

        # Mass
        mj = mass[j]

        tmp = G * mj * r**3

        # Calculate accelerations
        acc[i,0] += tmp * dx
        acc[i,1] += tmp * dy
        acc[i,2] += tmp * dz

return np.asarray(acc)

位置整合：

def leapfrog(np.ndarray[np.float64_t, ndim=2] pos, np.ndarray[np.float64_t, ndim=2] vel, np.ndarray[np.float64_t, ndim=2] acc, np.ndarray[np.float64_t, ndim=1] mass):
cdef double dt = 1e-3 # Timestep

# The Leapfrog integration method
# v(t+1/2) = v(t) + a(t) x dt/2
# x(t+1) = x(t) + v(t+1/2) x dt
# a(t+1) = (G * m/((dx^2 + dy^2 + dz^2)^(3/2))) * dx * x
# v(t+1) = v(t+1/2) + a(t+1) x dt/2
vel += acc * dt/2
pos += vel * dt
acc = acceleration(pos, mass)
vel += acc * dt/2

return pos, acc

主模拟循环（Python）： 对于 _ 在范围内（Nt）： # 计算位置并获得新的加速度值 pos, acc = jumpfrog(pos, vel, acc, m)

绘图（Python）：

plt.scatter(pos_arr[:,0], pos_arr[:,1])

请帮我解决这个问题。

有关该错误的更多信息，请参阅图片：

两个物体相互排斥+直线运动而不是轨道运动

错误的 3D 视图。

编辑 1：速度 = 0 这是速度 = 0 (vel = np.zeros((N,3)) 的结果，其中 N=2。 红点是位置数组的第一个元素，绿点是最后一个点。

编辑 2：tmp * dx/2： 这是通过将分离向量除以 r 获得的结果。 更新代码：

acc[i,0] += tmp * dx/r
acc[i,1] += tmp * dy/r
acc[i,2] += tmp * dz/r

编辑 3：速度 = 1e-6 这就是设置两个速度 = 1e-6 的结果。它们保持静止。

Answer 1

问题解决了：使用更好的加速算法解决了问题。该代码现在可在GitHub 获得。现在唯一的问题是动画情节。请让我知道如何在 cmets 中为其设置动画。感谢一直以来的帮助！

以下是两体模拟的结果：