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 have gone through a couple of iterations of formulas. For some reason, potential energy looks too high, and kinetic energy gets even higher at the lowest part of its swing. I’m not the most veteran coder, and I’m not even sure if this is the right site, but I’m getting really cross-eyed by this, and I don’t know what’s causing it to fail so dramatically.
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.
I thought to try adjusting the velocity dependent on the difference in energy from an initial constant, but that created some weird effects. Speeding up the pendulum far more than is realistic led me to assume my energy calculations were written wrong.