C++模拟库仑力核碰撞 最小间距处终止循环问题修复

原有代码核心问题 #

• 循环判断条件d >= d为恒真逻辑（浮点数为NaN时除外），既无法实现“间距持续减小则迭代、间距开始增大则终止”的需求，还大概率导致死循环
• 时间步长选择使用12层嵌套if，冗余度极高，可维护性差
• 全局时间变量t未初始化，初始值为内存随机垃圾值，计算结果无意义
• 幂运算统一用pow处理平方项，精度和效率都不如直接乘法

具体修正点 #

1. 循环逻辑重构：新增变量存储上一轮迭代的间距值prev_d，每完成一次位置更新，就对比新间距和上一轮间距：只要新间距更小，说明两核仍在靠近，继续迭代；一旦新间距大于等于上一轮间距，说明已经越过最近距离点，直接终止循环，此时上一轮存储的间距值就是所求的最小间距。
2. 步长计算简化：去掉多层嵌套分支，通过当前间距的数量级直接计算对应时间步长，和原分支逻辑完全等价但代码量大幅缩减。
3. 细节补全：初始化未赋值的全局变量，平方计算替换为直接乘法，增加数值误差容错阈值，避免浮点截断误差导致的循环异常。

修正后完整代码 #

#include <iostream>
#include <cmath>
using namespace std;

void At_col();

// 物理常数
const long double m_He = 6.64e-27;
const long double m_Li = 1.16e-26;
const long double K = 2.31e-23;
const int dk = 13;
const long double err_threshold = 1e-16;

// 运动状态变量
long double a_He, a_Li;
long double v_He = 4.39e6;
long double v_Li = 0;
long double x_He = 49;
long double x_Li = 50;
long double d = x_Li - x_He;
long double dt;
long double t = 0; // 补全初始化

// 根据当前间距计算适配的时间步长，替代原多层嵌套if
long double calc_dt(long double curr_d) {
    if (curr_d > 0.1) return powl(10, -dk);
    if (curr_d > 1e-2) return powl(10, -(dk+1));
    if (curr_d > 1e-3) return powl(10, -(dk+2));
    if (curr_d > 1e-4) return powl(10, -(dk+3));
    if (curr_d > 1e-5) return powl(10, -(dk+4));
    if (curr_d > 1e-6) return powl(10, -(dk+5));
    if (curr_d > 1e-7) return powl(10, -(dk+6));
    if (curr_d > 1e-8) return powl(10, -(dk+7));
    if (curr_d > 1e-9) return powl(10, -(dk+8));
    if (curr_d > 1e-10) return powl(10, -(dk+9));
    if (curr_d > 1e-11) return powl(10, -(dk+10));
    if (curr_d > 1e-12) return powl(10, -(dk+11));
    if (curr_d > 1e-13) return powl(10, -(dk+12));
    return powl(10, -(dk+13));
}

int main() {
    long double prev_d;
    // 执行第一次迭代拿到初始比较值
    dt = calc_dt(d);
    At_col();

    // 核心循环：间距持续减小就继续迭代，留误差阈值避免浮点震荡
    while (d < prev_d - err_threshold) {
        prev_d = d;
        dt = calc_dt(d);
        At_col();
    }

    cout << "\n=== 计算完成，两核最小间距为：" << prev_d << " m ===" << endl;
    return 0;
}

void At_col() {
    // 平方项直接用d计算，替代pow提升精度和效率
    long double d_sq = d * d;
    a_He = -K / m_He / d_sq;
    a_Li = K / m_Li / d_sq;

    // 半隐式欧拉积分更新状态
    v_Li += a_Li * dt;
    x_Li += v_Li * dt;
    v_He += a_He * dt;
    x_He += v_He * dt;

    t += dt;
    d = x_Li - x_He;

    cout << "t = " << t << "|d = " << d << "|x_He = " << x_He << "|x_Li = " << x_Li
         << "|v_Li = " << v_Li << "|a_Li = " << a_Li << "|v_He= " << v_He << endl;
}

补充说明 #

如果需要更高的计算精度，可以把半隐式欧拉积分替换为Verlet积分，能大幅减小长时迭代的数值漂移；如果需要找更精确的最近点，也可以在循环终止后，在最后两个迭代点之间做一次线性插值，进一步缩小最小间距的误差。

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