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).