solving the gravitational equations

A brief description of our approach to solving the gravity equations describing the motion of the spheres

Once we finally decided to move the spheres in our application around based on gravity (and therefore call it Newton’s Playground!), it felt strange to return to trying to solve differential equations after all these years. More than 30 years ago I did a PhD on star formation and the methods used for the calculation were based on these sort of N-Body calculations. For the next 15 years or so I made my living solving equations related to the physics of optics, lasers and communication. However for the last 15 years any calculations have been in optimisation (simulated annealing and those sort of techniques) and shortest/diverse-path routing calculations - which are far removed from differential equations. There was also a period not too long ago doing commercial stuff but in the end the pull of doing technical work again was too strong. Believe it, or not optimisation and routing algorithms are very interesting (oh well..). Anyway enough …

Er… so lets solve some differential equations then. We don’t really need to do this too accurately as we are only moving balls around to make them interesting. I know, Runge-Kutta techniques are sort of standard for differential equations so lets use this. Actually even with a relatively simple scheme like this it is surprisingly easy to screw up the code, as the forces need to be calculated based on the sphere position at each time step. At least I found it surprisingly easy to screw up…

The only way to make sure you don’t screw up is to calculate independently the total energy in the system (Kinetic + (negative) Potential) energy. as the bodies move, energy swaps between Kinetic and Potential energy but their total remains constant (or at least should remain constant). So what you do is to watch the total energy changing rapidly as the calculations progress - then go back to the code and take out typos. When it stops changing rapidly and stays pretty constant then you are getting somewhere. It is still quite easy to screw up more subtly though. You can easily make a small error on the Runge-Kutta formula. Then it still sort of work but you lose the nice 4th order accuracy that it is supposed to have. What you have to do is to do some experiments by halving the step size h. The error (e.g. in the total Energy) should then go down by a factor of 16 (h4). So you go back to the code and take out more typos until that is what happens.

Once you do this it works reasonably well, though you have to change the step size after each interval based on the evolving forces. However, the gravitational n-Body problem is special and when it goes wrong it does so in style with bodies getting trapped in tight binaries and then eventually flying off into infinity with any small inaccuracy. So even while it may not matter that the calculations are too accurate, the problem doesn’t just nicely be a “bit inaccurate” but goes crazy. Hence once started we need to do this properly.