为何我的n段摆无法保持能量守恒？代码排查求助

我正在为YouTube用户@maciejmatyka实现一个能维持能量守恒的n段摆，原本认为数值误差应该微小且可预测，于是编写了如下JavaScript Pendulum类。但运行过程中摆的能量持续流失，我无法确定代码是否存在根本性问题。

我尝试过一种简易的能量修正方法，虽然能勉强维持能量守恒，但逻辑不够严谨，现寻求更合理的解决方案。

Pendulum类实现代码 #

class Pendulum {
    constructor({
        n = 20,
        g = -10,
        dt = 0.002,
        thetas = Array(n).fill(0.5 * Math.PI),
        thetaDots = Array(n).fill(0)
    } = {}) {
        this.n = n;                 // 摆的段数（常量）
        this.g = g;                 // 重力加速度（常量）
        this.dt = dt;               // 时间步长（delta-time）
        this.thetas = thetas;       // 各段摆的角度（极坐标）
        this.thetaDots = thetaDots; // 各段摆的角速度（极坐标）
    }
    
    matrixA() {
        const { n, thetas } = this;
        return Array.from({ length: n }, (_, i) =>
            Array.from({ length: n }, (_, j) =>
                (n - Math.max(i, j)) * Math.cos(thetas[i] - thetas[j])
            )
        );
    }
    
    vectorB() {
        const { n, thetas, thetaDots, g } = this;
        return thetas.map((theta, i) => {
            let b_i = 0;
            for (let j = 0; j < n; j++) {
                const weight = n - Math.max(i, j);
                const delta = thetas[i] - thetas[j];
                b_i -= weight * Math.sin(delta) * thetaDots[j] ** 2;
            }
            b_i -= g * (n - i) * Math.sin(theta);
            return b_i;
        });
    }
    
    accelerations() {
        const A = this.matrixA();
        const b = this.vectorB();
        return math.lusolve(A, b).flat();
    }
    
    leapfrogStep() {
        const { thetas, thetaDots, dt } = this;
        const acc = this.accelerations();

        // 速度半步更新
        const halfThetaDots = thetaDots.map((dot, i) => dot + acc[i] * dt / 2);

        // 位置全步更新
        this.thetas = thetas.map((theta, i) =>
            ((theta + halfThetaDots[i] * dt) + Math.PI) % (2 * Math.PI) - Math.PI
        );

        // 速度全步更新
        const newAcc = this.accelerations();
        this.thetaDots = halfThetaDots.map((dot, i) => dot + newAcc[i] * dt / 2);
    }
    
    tick() {
        // 原本计划在这里修正速度以匹配能量
        this.leapfrogStep();
        // 检查能量
        // 应用修正
    }
    
    kineticEnergy() {
        const { thetas, thetaDots } = this;
        
        // 计算所有摆段的动能之和：Σ(1/2 * ((xDot_j)^2 + (yDot_j)^2))
        return thetas.reduce((T, _, i) => {
            let xDot = 0, yDot = 0;
            for (let j = 0; j <= i; j++) {
                xDot += thetaDots[j] * Math.cos(thetas[j]);
                yDot += thetaDots[j] * Math.sin(thetas[j]);
            }
            return T + 0.5 * (xDot*xDot + yDot*yDot);
        }, 0);
    }
    
    potentialEnergy() {
        const { thetas, n, g } = this;
        
        // 计算所有摆段的势能之和
        return thetas.reduce(
            (V, theta, i) => V - g * (n - i) * (Math.cos(theta)+1),
            0
        );
    }
    
    totalEnergy() {
        return this.kineticEnergy() + this.potentialEnergy();
    }

    get coordinates() {
        let x = 0, y = 0;
        return this.thetas.map(theta => {
            x += Math.sin(theta);
            y += Math.cos(theta);
            return { x, y };
        });
    }
}

临时能量修正代码 #

tick() {
        this.leapfrogStep();
        
        // 修正速度
        const curEnergy = this.totalEnergy();
        const energyDifference = curEnergy - this.initEnergy;
        
        const tolerance = 1e-5;
        if (Math.abs(energyDifference) > tolerance) {
            const scalingFactor = Math.sqrt(this.initEnergy / curEnergy);
            this.thetaDots = this.thetaDots.map(dot => dot * scalingFactor);
        }
    }

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