引力 N 体模拟的问题

Question

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

条件

质量：10，2 个物体 100； 1 表示所有物体，物体数量 > 2。
职位：随机生成
速度：随机生成
时间步长：1e-3

错误

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

预期行为

两个天体必须相互环绕
势力一定要有吸引力

努力解决bug

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

代码加速函数（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 = leapfrog(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 获得。现在唯一的问题是动画情节。请在评论中告诉我如何制作动画。感谢大家一直以来的帮助！

以下是二体模拟的结果：

引力 N 体模拟的问题

Problems with gravitational N-body simulation

python

cython

physics