Leap Frog算法求解两体问题程序输出异常，求排查原因

你遇到的问题很典型——相同初始条件下其他算法能跑对，唯独Leap Frog出问题，肯定是算法的实现细节没处理对。我帮你梳理下代码里的几个关键错误：

一、代码里的明显语法/命名问题 #

首先看你的加速度函数定义：

def a(j,i,2):

Python函数的参数不能是字面量（比如这里的2），得用变量名，比如n_particles；另外函数内部又用了a作为变量名，这会和函数名冲突，导致后续调用出错，得改成别的名字比如acc。修正后的函数应该是：

def a(j, i, n_particles):
    # j是目标粒子索引，i是时间步索引，n_particles是粒子总数
    acc = 0.0
    for k in range(n_particles):
        if k != j:
            delta_r = r[j,i] - r[k,i]
            norm_r = np.linalg.norm(delta_r)
            acc -= G * m[k] * delta_r / (norm_r ** 3)
    return acc

二、Leap Frog算法的核心实现错误 #

你现在的代码是逐个粒子更新位置后立刻计算加速度，这是致命错误！因为两体问题中每个粒子的加速度依赖于另一个粒子的当前位置——如果你先更新了粒子0的位置，粒子1的位置还停留在上一时间步，此时计算粒子0的加速度时用的是新旧混合的位置，结果自然不对。

标准的速度-Verlet（属于Leap Frog算法家族）的正确步骤应该是：

1. 先计算所有粒子的初始加速度
2. 对每个时间步：
   • 先批量更新所有粒子的位置
   • 再批量计算所有粒子的新加速度（此时所有粒子的位置都是当前时间步的）
   • 最后批量更新所有粒子的速度

三、修正后的完整代码 #

import numpy as np
import math
from matplotlib import pyplot as plt

# ~ ~ ~ ~ ~ ~ FUNCTIONS ~ ~ ~ ~ ~ ~
def a(j, i, n_particles):
    # j是目标粒子索引，i是时间步索引，n_particles是粒子总数
    acc = 0.0
    for k in range(n_particles):
        if k != j:
            delta_r = r[j,i] - r[k,i]
            norm_r = np.linalg.norm(delta_r)
            acc -= G * m[k] * delta_r / (norm_r ** 3)
    return acc

# ~ ~ ~ ~ ~ ~ MAIN BODY ~ ~ ~ ~ ~ ~
Dt = 300
h = 0.01
N = int(Dt/h)
G = 1  # 单位简化为1
m = np.ones(2, float)  # 两个粒子质量都是1
t = np.zeros(N, float)
r = np.zeros((2, N, 2), float)  # 位置数组：[粒子索引, 时间步, x/y维度]
u = np.zeros((2, N, 2), float)  # 速度数组

# 初始条件
r[0,0] = [1.0, 1.0]
r[1,0] = [-1.0,-1.0]
u[0,0] = [-0.5, 0.0]
u[1,0] = [0.5, 0.0]

# 先计算初始加速度
a_prev = np.zeros((2, 2), float)
for j in range(2):
    a_prev[j] = a(j, 0, 2)

for i in range(1, N):
    t[i] = t[i-1] + h
    # 1. 批量更新所有粒子的位置
    for j in range(2):
        r[j,i] = r[j,i-1] + h * u[j,i-1] + 0.5 * h**2 * a_prev[j]
    # 2. 批量计算所有粒子的新加速度
    a_curr = np.zeros((2, 2), float)
    for j in range(2):
        a_curr[j] = a(j, i, 2)
    # 3. 批量更新所有粒子的速度
    for j in range(2):
        u[j,i] = u[j,i-1] + 0.5 * h * (a_prev[j] + a_curr[j])
    # 更新加速度，用于下一个时间步
    a_prev = a_curr.copy()

# ~ ~ ~ ~ ~ ~ ~ PLOTS ~ ~ ~ ~ ~ ~ ~
plt.plot(r[0,:,0], r[0,:,1], 'black', label='Particle 0')
plt.plot(r[1,:,0], r[1,:,1], 'red', label='Particle 1')
plt.legend()
plt.axis('equal')  # 保证x/y轴比例一致，轨迹显示更准确
plt.show()

四、额外优化建议 #

• 加上plt.axis('equal')：两体问题的轨迹是对称的，用等比例轴能更直观看到正确的轨道形状
• 给曲线加上图例，方便区分两个粒子

内容的提问来源于stack exchange，提问作者Δημήτρης Μ