In response to @maciejmatyka on YouTube, I was attempting to create an n-segment pendulum that can conserve energy. Obviously, numerical issues would cause discrepancies but I assumed those would be small and predictable, so my goal was to make a good code to single that out, then fix that accordingly. Here is my pendulum class in JS:

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;                 // Number of pendulums (constant)
        this.g = g;                 // Gravitational acceleration (constant)
        this.dt = dt;               // Time step (short for delta-time)
        this.thetas = thetas;       // Angles of the pendulums (polar)
        this.thetaDots = thetaDots; // Angular velocities (polar)
    }
    
    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();

        // Half-step for velocities
        const halfThetaDots = thetaDots.map((dot, i) => dot + acc[i] * dt / 2);

        // Full step for positions
        this.thetas = thetas.map((theta, i) =>
            ((theta + halfThetaDots[i] * dt) + Math.PI) % (2 * Math.PI) - Math.PI
        );

        // Full step for velocities
        const newAcc = this.accelerations();
        this.thetaDots = halfThetaDots.map((dot, i) => dot + newAcc[i] * dt / 2);
    }
    
    tick() {
        // This is where I would correct velocities to match the energy
        this.leapfrogStep();
        // Check energy
        // Apply difference
    }
    
    kineticEnergy() {
        const { thetas, thetaDots } = this;
        
        // Sum of 1/2 * ((xDot_j)^2 + (yDot_j)^2) for j in n
        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;
        
        // Sum of cos(theta)+1's up to i for i in n
        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 };
        });
    }
}

I was expecting a pendulum to stay stable for at least a few hundred iterations, but it seems to be constantly losing energy, which I can't explain unless something is fundamentally wrong with my code.

Here is a simple integration that maintains conservation of energy but its not particularly rigorous so I'm still looking for a good solution:

    tick() {
        this.leapfrogStep();
        
        // Fix velocities
        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);
        }
    }