$N$-body gravity simulator: why does energy conservation break down when introducing an adaptive timestep?

I am playing with an N-body gravity simulator using the velocity-verlet algorithm. In the actual simulation I normalize everything so that it is nicely behaved numerically, but I'll express things in real units here. My acceleration on body $$i$$ is $$\vec{a_i}=\sum_{j\neq i}\frac{Gm_j}{||\vec{r_{ij}}||^3}\vec{r_{ij}}$$ The basic, constant timestep implementation works well, as long as things don't get too close together, but when the acceleration gets very large, i.e. when two bodies approach very close to one another, it tends to blow up.

No problem, I'll use an adaptive timestep. Because this is a second-order integrator, I expect the error to be proportional to the magnitude of the acceleration times the square of the timestep. So if I pick an acceptable error level on each update, I can just rescale the timestep based on the smallest nearest-neighbour distance in the simulation at any given time, and everything works well again. But now I have the problem that I'm spending a lot of time using tiny timesteps to resolve the motion of closely-orbiting bodies, when most of the bodies in the system are far from anything else and do not require nearly that level of precision.

I came up with an idea for using an adaptive timestep which is different for every body in the system, using a small timestep for closely orbiting bodies and a large one for far away bodies, and taking timesteps as needed in order to time-sync everything periodically with a large timestep $$\Delta t_0$$. In order to do that, I find, for each body, the timestep that satisfies the acceptable error criterion I defined above, and I round it down to the nearest real number that is of the form $$\Delta t_n=2^{-n}\Delta t_0$$ for some integer $$n$$. I then apply the verlet algorithm recursively to particles with a different exponent $$n$$. The algorithm used to advance the entire simulation by $$\Delta t_0$$ (one call to move()) looks like this:

def move(order, maxorder, bodies, dt0):
if order == maxorder:
velocity_verlet(order, bodies, dt0)
else:
move(order + 1, maxorder, bodies, dt0)
velocity_verlet(order, bodies, dt0)
move(order + 1, maxorder, bodies, dt0)

def velocity_verlet(order, bodies, dt0):
acceleration(order, bodies)
for b in bodies:
if b.order == order:
b.update_position(dt0)
b.half_update_velocity(dt0)
acceleration(order, bodies)
for b in bodies:
if b.order == order:
b.half_update_velocity(dt0)


acceleration(order) just updates the acceleration term for all particles of exponent (order) $$n$$ based on their positions after the internal update.

In a nutshell, I advance all particles of the smallest timestep by their timestep, then advance all the next-larger timestep particles, then go back and advance all the smallest-timestep particles again. Then move up a level recursively and apply this again. Of course, any physical observable can only be measured when all particles have synced up again, that is, after the top-level call to move() has been resolved.

This also works, for some simple cases over short timescales, except that energy is no longer conserved and it will eventually blow up. I am hoping that someone can help me understand why this is the case, and what, if anything, can be done to rescue it.

If this works it would be quite nice, as I get the best of both worlds: high computational effort is expended only exactly where it is needed, resolving the trajectories of particles deep in gravity wells, while leaving particles in shallower parts of the potential curve to take large steps through space without comprising accuracy. But clearly the issue of time-synchrony is causing nonphysical effects even after particles have been re-synced. Any suggestions?

if anyone is curious to play with it, you can find the code here.

• Would Computational Science be a better home for this question? (Not that I think it's off-topic here.) – Emilio Pisanty Mar 27 at 20:42
• It's definitely on the edge. Didn't even know that one existed. I have no objection to it being moved there if you think it should be. – KBriggs Mar 27 at 20:54
• I did very similar algorithm few years ago, for 10000 particles on a notebook, but using leapfrog instead of Verlet. The individual adaptive step algorithm worked very well, energy did not change much at all, apart from sudden jumps here and there which I thought were due to insufficient numerical accuracy associated with close encounters. Do you really see energy in your simulation changing much? – Ján Lalinský Mar 27 at 22:22
• @EmilioPisanty No, I do not think that Computational Science would be better for this specific question. I am sure that if there is a reasonable answer to this question it could only come from physicists. Verlet algorithm is nice, apparently simple, but its full understanding requires good knowledge of analytic mechanics, Hamiltonian systems and symplectic transformations. – GiorgioP Mar 27 at 23:13
• Interesting concept. Just a wild guess, but maybe it'd help if you interpolate the positions of the large timestep bodies when you update the small timestep ones. Another option is to use synchronized leapfrog, with Yoshida coefficients. See the artcompsci link in note #5 in the Wikipedia article for details & discussion. You can easily use recursion on the Yoshida steps to achieve higher orders of integration. IME, leapfrog + Yoshida is impressive, but I haven't tried it in your scenario with a mix of timesteps. – PM 2Ring Mar 28 at 7:01