蛙跳积分法是否适用于三体问题？附实现代码求验证

我一直在尝试用蛙跳积分法记录三体问题（3BP）的哈密顿量随时间的变化，但始终没搞懂标准半步法的实现方式，于是自己写了一个变体，不确定它是否正确。以下是我用到的变量和函数，其中p、v、m都是numpy数组：

p = [array([x,y]), array([x,y]), array([x,y])]
v = [array([x,y]), array([x,y]), array([x,y])]
m = [1, 1, 1]

完整代码如下：

from numpy import sum
from numpy.linalg import norm
from copy import deepcopy

def H(p, v, m): #hamiltonian function
    #sum of kinetic energy for all bodies
    T = sum([m[i]*norm(v[i])**2/2 for i in range(3)])
    #sum of potential energy between all bodies 
    V = -sum([m[0-i]*m[1-i]/norm(p[0-i]-p[1-i]) for i in range(3)])
    return T + V

def a(p, n, m): #sum of the acceleration arrays for body n and the two other bodies
    return m[n-1]*(p[n-1]-p[n])/(norm(p[n-1]-p[n])**3) + m[n-2]*(p[n-2]-p[n])/(norm(p[n-2]-p[n])**3)

def collision(p): #checks for collisions
    for i in range(3):
        if norm(p[0-i] - p[1-i]) < 0.1:
            return True

#leapfrog
def Leapfrog(P, V, dt, steps, m):
    p,v=deepcopy(P),deepcopy(V)
    H_L = [H(p, v, m)]

    for t in range(steps):
        atemp = [a(p, i, m) for i in range(3)] #acceleration at time step i

        #calc new values
        for i in range(3):
            p[i] = p[i] + v[i]*dt + 0.5*atemp[i]*dt**2

        for i in range(3):
            v[i] = v[i] + 0.5*(atemp[i] + a(p, i, m))*dt #acceleration at timestep i+1


        if collision(p):
            return H_L
        H_L.append(H(p, v, m))
    return H_L

代码问题分析 #

1. 势能计算的写法隐患
   你写的势能计算虽然结果正确，但通过负索引循环的方式可读性差，容易出错。更清晰的写法是直接枚举所有两两不同的天体组合，避免依赖索引的负数取值。

2. 加速度函数的注释错误
   a函数的注释描述有误，它实际计算的是其他两个天体对第n个天体的加速度之和，而非“第n个天体和另外两个天体的加速度总和”。

3. 你的变体并非标准蛙跳
   标准蛙跳的核心是位置和速度错开半步更新，步骤是：

   • 先把速度推进半步：v(t+dt/2) = v(t) + a(t)*dt/2
   • 用半步速度更新位置：p(t+dt) = p(t) + v(t+dt/2)*dt
   • 计算新位置的加速度a(t+dt)，再把速度推进到全步：v(t+dt) = v(t+dt/2) + a(t+dt)*dt/2

   而你的实现是用泰勒展开更新位置，再用梯形法更新速度，本质是Verlet积分的一种形式，并非标准蛙跳。这种变体的哈密顿量守恒性不如标准蛙跳——标准蛙跳是辛积分，能在长时间模拟中更好地保持能量守恒，而你的写法可能会出现能量随时间漂移的情况。

标准蛙跳积分实现 #

下面是修正后的标准蛙跳代码，同时优化了部分逻辑的可读性：

from numpy import sum
from numpy.linalg import norm
from copy import deepcopy

def H(p, v, m):
    # 计算总动能
    kinetic = sum([m[i] * norm(v[i])**2 / 2 for i in range(3)])
    # 计算总势能（枚举所有两两组合，避免重复计算）
    potential = - (
        m[0]*m[1]/norm(p[0]-p[1]) 
        + m[0]*m[2]/norm(p[0]-p[2]) 
        + m[1]*m[2]/norm(p[1]-p[2])
    )
    return kinetic + potential

def get_acceleration(p, n, m):
    # 计算第n个天体受到的合加速度（引力常数G=1）
    acc = 0.0
    for i in range(3):
        if i != n:
            r_vec = p[i] - p[n]
            r_norm = norm(r_vec)
            acc += m[i] * r_vec / (r_norm ** 3)
    return acc

def check_collision(p):
    # 检查任意两个天体是否发生碰撞
    for i in range(3):
        for j in range(i+1, 3):
            if norm(p[i] - p[j]) < 0.1:
                return True
    return False

def standard_leapfrog(init_p, init_v, dt, steps, m):
    p, v = deepcopy(init_p), deepcopy(init_v)
    hamiltonian_history = [H(p, v, m)]
    
    # 初始化加速度，先把速度推进半步
    current_acc = [get_acceleration(p, i, m) for i in range(3)]
    for i in range(3):
        v[i] += current_acc[i] * dt / 2
    
    for _ in range(steps):
        # 用半步速度更新位置
        for i in range(3):
            p[i] += v[i] * dt
        
        # 计算新位置的加速度
        new_acc = [get_acceleration(p, i, m) for i in range(3)]
        
        # 把速度推进到全步
        for i in range(3):
            v[i] += new_acc[i] * dt / 2
        
        if check_collision(p):
            return hamiltonian_history
        
        # 记录当前时刻的哈密顿量
        hamiltonian_history.append(H(p, v, m))
    
    return hamiltonian_history

验证建议 #

• 对比你的变体和标准蛙跳的哈密顿量曲线，标准蛙跳的能量波动应该更小，且不会随时间出现明显漂移
• 用已知的稳定三体构型（比如等边三角形解、拉格朗日点解）测试，观察模拟结果是否能长期保持稳定

内容的提问来源于stack exchange，提问作者Nnx_0C