Symplectic Integrators

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.

Runge-Kutta (4th Order Nonsymplectic)