Unless otherwise stated, the term ``system'' in this section will mean an autonomous Hamiltonian system. The discussion of symplectic integration for Hamiltonian systems in this section draws heavily from Sanz-Serna [66], which is a very readable introduction to the subject. Another seminal paper is Channel and Scovel [13].
The question of what kind of integrator to use to integrate N-body systems is one that has had plenty of attention. No discussion of stellar system N-body integrators would be complete without mention of Aarseth's classic integration scheme for collisional systems [1], which includes regularization. Makino [55] derived formulae for optimal order and timestep criteria for Aarseth-type integrators. For collisionless N-body systems, and indeed many non-astronomical Hamiltonian systems, symplectic integrators have enjoyed widespread popularity. Recent years have also seen a growing interest in exactly conserving integrators.
The popularity of symplectic integrators
stems from the fact that the time-t flow of any Hamiltonian system is
a symplectic or canonical mapping. One consequence of this
is that the mapping is volume preserving in phase space. This is a
consequence of Liouville's Theorem (see, for example,
[56]), which states that the phase-space density of
trajectories inside an ensemble remains constant with time. A
symplectic integrator also preserves phase space volume to within the
machine precision. As a result, symplectic integrators tend to preserve
qualitative properties of phase space trajectories: trajectories
do not cross, and
although energy is not exactly conserved, energy fluctuations are
bounded. In fact, it can be shown that the phase space of a linear
Hamiltonian problem simply undergoes slight stretches and/or contractions
along various axes, so for example a circular trajectory in 2 dimensions
becomes an ellipse.
Furthermore, if the Hamiltonian function is sufficiently
smooth, then the discrete solution produced by a symplectic integrator
lies on the exact solution of a problem that has a slightly perturbed
Hamiltonian. In general, even for non-smooth Hamiltonians,
the discrete solution lies exponentially close in 
to an exact
solution to a nearby Hamiltonian problem, where is the local
error [66].
Non-symplectic integrators generally turn a conservative system into a dissipative one, while a symplectic integrator, on average over time, preserves the conservative nature of the system. Earn and Tremaine [23] point out that roundoff error actually destroys the exact symplecticness of a time- and/or space-continuous system. They demonstrate how one can build a precisely symplectic discrete map over a lattice discretization of space using integer arithmetic, so that if the modeled problem is Hamiltonian, then the discrete model is also exactly Hamiltonian, but with a slightly different Hamiltonian function. This builds confidence in the qualitative nature of the solution being similar to that of the continuous system.
Since symplectic integrators preserve the topological structure of trajectories in phase space, they are probably more reliable for long-term integrations than non-symplectic integrators, since a non-symplectic integrator causes the system to become dissipative, which in turn can badly distort the topological structure of the phase flow over long periods. However, symplectic integrators are often implicit, and require more function evaluations and smaller timesteps to produce the same computational accuracy as compared to standard non-symplectic integrators [66]. Thus, for short-term integrations where accuracy is paramount, a standard non-symplectic integrator may be more efficient.
Another drawback of symplectic integration is that any integration scheme
which is symplectic for constant timestep does not remain symplectic
if the timestep is changed based solely and explicitly on the state vector
[71].
Recall that the time- 
map produced by a symplectic integrator
lies exactly (or exponentially close to) the true time- flow
of a slightly perturbed Hamiltonian problem. However, if
changes, then the problem we are close to also changes. Thus, if we
vary the timestep, it appears that there is no single Hamiltonian
problem our symplecticly-generated solution remains close to.
Nonetheless, several researchers have begun to experiment with variable timestep symplectic integrators. Hut, Makino and McMillan [40] introduce an implicit method of symmetrization which is intended to preserve approximate time symmetry. For the Leapfrog integrator, their method reduces to
where h(y) is a function used to compute the timestep, based solely
on the state vector.
They tested their symmetrized leapfrog on a Kepler problem
with eccentricity e=0.9, using 1 implicit iteration per timestep and
timesteps per revolution. The errors in energy, semi-major
axis, eccentricity, time of perihelion passage and longitude of
pericenter all have moderately sized periodic errors, but secular
drift is nil or linear with a small slope,
compared to the non-symmetrized error of
linear or quadratic, respectively.
They also demonstrate how a 4th-order generalization
of Leapfrog can use their time-symmetrization technique,
and use it to perform an impressive integration of an orbit with eccentricity
. With 3 iterations per
timestep of the implicit time-stepping symmetrization equation, they
were able to integrate this trajectory using
timesteps per orbit, with a drift in the above quantities
of only one part in 
, over revolutions. This
seems very impressive. They do not mention any theoretical results
on the symplecticness of their method, however.
Skeel and Biesiadecki [72] introduce another method for using different timesteps while preserving the symplectic nature of the integrator. They show that, if it is possible to split the Hamiltonian into independent additive terms, then each term may be integrated with a different, but constant, timestep, and the integration remains exactly symplectic. Each timestep must be an integer multiple of the smallest timestep. For example, if there are only two step sizes and the large timestep is a multiple M of the smaller, then the force at timestep n is computed as
They show that
an artificial splitting of the Hamiltonian into ``quickly''
and ``slowly'' varying components can also utilize their method.
They test their method on the Kepler problem,
setting the ``quickly'' varying component to be the half of the
orbit that is closest to the central mass, and the slowly varying component
to be the other half of the orbit. Their results are similar
to the second-order results of Hut et al. [40].
It would be interesting to see if a 4th order version of their
method would also give comparable results to Hut et al.
[40].
The error of this method is analyzed more fully in Littell, Skeel, and
Zhang [54], where they come to the unsurprising
conclusion that the smallest stepsize times the highest frequency of
the ``quickly'' changing components must be small. In other words,
their method still needs an adequate sample of the highest frequency
changes in order to integrate them. They prove the also unsurprising
result that, if 
is
the maximum frequency in the system, integrated with timestep h,
and 
is the maximum frequency of one of the ``slower''
components, then the most efficient timstep for that component is
.
In addition to symplectic integration, there is an active and quickly blooming body of literature on exactly conserving algorithms. In general, there does not exist an integrator that is both symplectic and exactly energy conserving for Hamiltonian systems [69, 71]. However, since energy is not the only quantity conserved in most systems, there is still the possibility of building a symplectic integrator that conserves other quantities, although there are limits to what the integrator can conserve in comparison to the real system [25]. Furthermore, if the system is not Hamiltonian, then a symplectic integrator provides no added benefit, whereas a conserving integrator may still be of great value. One must be careful, however, how one enforces energy conservation. It is plausible that a method which naïvely enforces global energy conservation may introduce spurious mechanisms of energy transport within the system [69].
Shadwick, Bowman, and Morrison [69] introduce a non-traditional integration scheme which is explicit, and exactly conserving of all quantities one is willing to program it to conserve. The drawback is that the method is not general: one needs to derive a complicated set of explicit formulae for each new integrator-problem pair. One applies the integration method to one's problem on paper, and then derives exactly how the conserved quantities of interest are not conserved. Then one simply adds terms to the integration routine that compensate for the non-conservation. If the standard version of the integrator is explicit, then so is the modified version. Although it is a very tedious process, it seems to have impressive results. The authors present 3 examples: the ``three-wave'' problem, the Lotka-Volterra predator-prey problem, and the Kepler problem, each demonstrated with both Euler's method and a 2nd order predictor-corrector. In the Kepler problem they force conservation of energy, angular momentum, and a quantity known as the Runge-Lenz vector. They note that non-conservation of phase-space volume in Hamiltonian systems can be used as a measure of error, since one can no longer use conservation of other quantities as a measure of error. The authors, however, do not clearly state the limitations of their method. For example, it is not clear precisely why one couldn't use their method to force conservation of phase space volume as well as energy, which we know to be impossible [25, 71].
We will close our discussion of integrators with a remark by Sanz-Serna [66, p.278,] (d is the dimension of configuration space):
Conservation of energy restricts the dynamics of the numerical solution by forcing the computed points to be on the correct (2d-1)-dimensional manifold H= constant, but otherwise poses no restriction to the dynamics: within the manifold the points are free to move anywhere and only motions orthogonal to the manifold are forbidden. When d is large this is clearly a rather weak restriction. On the other hand, symplecticness restricts the dynamics in a more global way: all directions in phase space are taken into account.Athough I think he makes an extremely interesting point, I note that his remark does not address the possibility that Shadwick et al. [69] have demonstrated: that energy is not the only quantity that can be conserved. Thus, it is this author's opinion that much interesting work remains in the area of integration of conservative systems.