Integration of differential equations is maybe the most important tool that
numericists have in their toolkit. For astrophysicists, these are almost
always the equations of motion for some collection of objects (gas, stars,
planets, dark matter particles, etc.). These problems have well-defined
Hamiltonians, and a number of conserved quantitities (energy, angular/linear
momentum, etc.). A kind of numerical integrator that automatically
conserves these quantities is called a Symplectic
Integrator. Since a Hamiltonian system has time-reversal symmetry, we should
also be able to run the integrator in reverse and reproduce the initial conditions to
within the errors from floating point machine precision.
Even a low-order symplectic integrator will conserve
the phase-space volume of a problem, and as such, never introduce or remove
energy systematically. Below is a demo of how a second order symplectic
integrator can outperform a fourth order non-symplectic integrator.
Beware, of course, that if your Hamiltonian is not formulated to conserve a
quantity, a symplectic integrator will not guarantee that quantity is
conserved (in this case, the phase angle of our orbit).
A Simple Orbital Problem
A particle orbiting another in a central inverse-square potential has the following Hamiltonian,
for the position vector $p$ and momentum vector $q$.
$$H = -\frac{1}{q}+\frac{p^2}{2}$$
Solving Hamilton's equations gives these equations of motion:
$$\dot p = {|q|}^{-2}$$
$$\dot q = p$$
These are the equations of motion for a planet orbiting a star, for example. Let's take a look at what
happens when we throw an elliptical orbit set of initial conditions at two different integrators, and see
how they conserve energy and angular momentum.