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 (2Athough I think he makes an extremely interesting point, I note that his remark does not address the possibility that Shadwickd-1)-dimensional manifoldH= 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. Whendis 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.

Fri Dec 27 17:41:39 EST 1996

Access count (updated once a day) since 1 Jan 1997: 16501