C++中的N体模拟具有很大的动量守恒和巨大的能量偏差

Question

我正在使用 Verlet 积分（更具体地说，Verlet 积分维基百科页面中“非恒定时间差”部分中的最后一个方程）与 C++ 和 Eigen 库来创建数值 N 体模拟。

我尝试测量模拟的动量和能量守恒。对于动量守恒，即使对于物体非常接近且快速的系统，我也得到了与最初计算的动量大约 10^(-7)% 的偏差。至于能量，我见过它一路跳到800%。它一开始很低，约为 0.0001%，但随着模拟的进行，它会变得更大，当物体非常接近时，它会达到 100%+。

（为了澄清，我所说的偏差是指 (y/x) - 1，其中 y 是当前值，x 是初始值。）

这似乎非常非常奇怪，我的程序保存动量如此之好，但能量却完全消失了。

根据我在 Meta StackOverflow 上读到的内容，不建议发布大块代码。我在这里包含了我正在做的事情的基本要点（Verlet 积分实现、重力计算）。如果有人想了解更详细的细节，这里是一个最小的可重现示例的 Github 存储库。

    void verletNextStep(double dt) {

        // as per wikipedia
        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;
    };

    //acceleration exerted by body2, on 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;
    };

    //executes verlet integration on all bodies
    void doVerlet(double dt) {

        for (Body &body : bodyList)
            body.acceleration = totalGravitationalAccelerationOf(body);


        for (Body &body : bodyList)
            body.verletNextStep(dt);

    }

我有理由怀疑不准确是由于双打溢出造成的。我将所有内容都改为长双打，并在能量低得多的系统上进行测试，但偏差仍然存在。 此外，另一篇文章提出了类似的问题，但他们的解决方案是有序计算加速度（不更新中间的主体），我已经实现了，如上面的代码所示。

Answer 1

Verlet 的可变时间步长的整个想法都是废话。准能量守恒特性关键取决于固定的时间步长，几乎恒定的代用能量具有取决于步长大小的扰动项。

动量要好一些，因为辛步骤准确地保留了状态变量中的线性和二次表达式。

在 n 体模拟中，通常无法避免物体接近碰撞，在这些时候，即使是固定步长 Verlet 也会崩溃，因为代用能量的扰动项比能量函数本身具有更差的奇异性，通常，系统会发生无意义的爆炸。

对于高偏心率的二体问题，即非常椭圆的非圆形轨道，也会得到类似的散度。