I am using Verlet integration (more specifically, the last equation in the "Non-constant time differences" section in the Verlet integration Wikipedia page), with C++ and the Eigen library, to create a numerical N-body simulation.
I've tried measuring the momentum and energy conservation of the simulation. For momentum conservation, even for systems where objects get very close and fast, I've gotten about 10^(-7)% deviation from the initially calculated momentum. As for the energy, I've seen it jump all the way to 800%. It starts off low, about 0.0001%, but gets much much bigger as the simulation goes on, and reaches up to 100%+ when objects are very close.
(to clarify, by deviation, I mean (y/x) - 1, where y is the current value and x is the initial.)
It seems very, very strange that my program conserves momentum so well, yet energy is completely off.
According to what I read on Meta StackOverflow, it's not recommended to post large chunks of code. I've included here the bare essentials of what I'm doing (Verlet integration implementation, gravitational force calculation). Here is a Github repo of a minimal reproducible example if anyone wants to look into the finer details.
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);
}
I had reason to suspect the inaccuracy was due to overflow of doubles. I changed everything to long doubles, and tested on systems that had much lower energies, yet the deviation still persists. In addition, another post presents a similar problem, but the solution in their case was orderly computation of acceleration (without updating bodies in-between), which I have implemented, as can be seen in the above code.
In the potential each interaction is only contained once in the potential energy. At the moment you compute the interaction as double iteration over the full list of bodies, so that the interaction potential for each pair is computed twice. So either correct that in the computation or divide the result by two.
By the way, IMO the whole idea of variable time steps with Verlet is bunk. The quasi-energy preservation property critically depends on the time step being fixed, the ersatz energy that is nearly constant has a perturbation term that depends on the step size. With a different step size you get a different ersatz energy that is kept constant, so the first ersatz energy can change. During periodic changes from smaller to larger step sizes and back these energy discrepancies can accumulate
Momentum behaves a little bit better as symplectic method steps preserve linear and quadratic expressions in the state variable exactly.
With n-body simulations you can usually not avoid that the bodies get close to near collisions, at these times even the fixed-step Verlet breaks down, as the perturbation terms of the ersatz energy have a worse singularity than the energy function itself, usually one gets non-sensical explosions out of the system.
One gets similar divergences for two-body problems with high eccentricity, that is, very elliptical, non-circular orbits.