C++基于Verlet积分的N体模拟：动量守恒佳但能量偏差极大

我用C++结合Eigen库，基于维基百科Verlet积分页面「Non-constant time differences」小节的最终公式实现了数值N体模拟。测试守恒量时发现：

• 动量守恒表现极佳：即使天体近距离高速运动，动量仅偏离初始值约10^(-7)%
• 能量偏差随模拟急剧恶化：初始偏差仅0.0001%，后期飙升至800%以上，天体近距离交汇时偏差甚至超过100%
  偏差定义为(当前值/初始值)-1

核心代码 #

Verlet积分步进实现 #

void verletNextStep(double dt) {
    // 参考维基百科公式实现
    if (!firstStepDone) {
        lastPosition = position;
        position = lastPosition
                   + velocity * dt
                   + 0.5 * acceleration * dt * dt;
        lastDt = dt;
        firstStepDone = true;
        return;
    }

    Vector2<long double> nextPosition = position
                            + (position - lastPosition) * dt/lastDt
                            + acceleration * dt * (dt + lastDt)/2;
    lastDt = dt;

    lastPosition = position;
    position = nextPosition;

    velocity = (position - lastPosition)/dt;
};

引力加速度计算与积分执行 #

Vector2<long double> totalGravitationalAccelerationOf(Body &body) {
    Vector2<long double> acc(0, 0);
    for (Body &otherBody : bodyList)
        if (body.vectorTo(otherBody).norm() != 0)
            acc += gravitationalAccelerationOf(body, otherBody);
    return acc;
};

// body2对body1的引力加速度
Vector2<long double> gravitationalAccelerationOf(Body &body1, Body &body2) {
    Vector2<long double> rHat = body1.vectorTo(body2).normalized();
    long double rSquared = body1.vectorTo(body2).squaredNorm();
    return rHat * G * body2.mass/rSquared;
};

// 对所有天体执行Verlet积分
void doVerlet(double dt) {
    // 先批量计算所有天体的加速度，避免中途更新状态影响计算
    for (Body &body : bodyList)
        body.acceleration = totalGravitationalAccelerationOf(body);

    // 再批量更新位置
    for (Body &body : bodyList)
        body.verletNextStep(dt);
}

已尝试的排查措施 #

• 将所有数值类型从double改为long double，并在低能量系统测试，能量偏差问题仍存在
• 遵循同类问题解决方案：先计算所有天体的加速度，再统一更新位置，避免中途状态更新的耦合干扰，问题未解决

原因分析与解决思路 #

1. 变步长Verlet的辛性破坏：标准Verlet积分是辛积分，能保证长期能量守恒有界，但维基百科中的变步长公式并非辛格式，会导致能量无界漂移。若使用固定步长Verlet，能量偏差应会控制在稳定范围内（短期有小波动，长期不发散）。
2. 近距离引力计算的数值误差：当天体距离极近时，rSquared趋近于0，除法运算会放大浮点误差，导致加速度计算严重失准，进而引入能量误差。可尝试：
   • 对近距离天体引入软核修正（rSquared += epsilon，epsilon为极小值），避免分母过小
   • 使用更高精度的浮点运算（如尝试GMP等任意精度库）
3. 速度计算的近似误差：当前速度通过(position - lastPosition)/dt计算，这是Verlet积分中速度的近似值（并非直接积分得到），变步长下该近似的误差会累积。可尝试在计算速度时加入加速度修正，比如用velocity = (position - lastPosition)/dt + 0.5*acceleration*dt（类似速度Verlet的速度更新方式）。
4. 步长选择不合理：天体高速运动时，固定步长无法捕捉引力场的快速变化，导致积分误差激增。可尝试自适应步长策略：根据天体当前加速度、速度调整步长，保证每步的局部误差在阈值内。

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