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

问：

我正在使用 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);

    }

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