Symplectic integrators revisited

Symplectic Integrators Revisited

T.J. Stuchi

Instituto de Física, Universidade Federal do Rio de Janeiro,

Caixa Postal 6528, 21945-970, Rio de Janeiro, RJ, Brazil

Received on 10 October, 2001. Revised version received on 10 April, 2002.

This paper is a survey on Symplectic Integrator Algorithms (SIA): numerical integrators designed for Hamiltonian systems. As it is well known, n degrees of freedom Hamiltonian systems have an important property: their ows preserve not only the total volume of the phase space, which is only one of the Poincaré invariants, but also the volume of sub-spaces less then 2n. These invariants are inherited from the conservation of the symplectic area. It is usually demanded of integrators that they should preserve energy. In this survey the main point is to convince the readers that the preservation of the symplectic area or canonicity of the Hamiltonian ow can be equally important, mainly when the concern is not one particular trajectory but the behavior of the phase space as a whole for long intervals of time. The KAM theorem asserts that for any integrable Hamiltonian perturbed by a small Hamiltonian term, such as that caused by the construction of the SIA, the perturbed dynamics preserves most of the incommensurate, nondegenerate, invariant tori. Unstable objects and their invariant manifolds are structurally stable and will be well represented by symplectic integrators.

I Introduction

A great number of phenomena are modelled by ordinary differential equations. When solved, analytically or numerically, they describe the time evolution of the quantities used to model the phenomena. Among these systems there are those called conservative or Hamiltonian; there also exists systems that are dissipative but have a Hamiltonian core, for example, the movement of an electron in an electron ring. The applications of these latter systems cover many areas such as engineering, plasmas, celestial mechanics, particle accelerators, and many problems originating from hydrodynamics, elasticity, relativity and quantum mechanics.

Analytical solutions for Hamilton's equations are exception rather than the rule. This means that there are very few integrable Hamiltonian systems, and we have to use numerical codes to solve the equations. Numerical codes reduce the differential equations to finite difference equations through algorithms which are now standard, for instance the Runge-Kutta family, linear multistep codes, extrapolation methods, Taylor series, to mention some of the main methods. Most of these algorithms are well known, not only practically but also theoretically. Their stability is a research area on its own, including techniques originated from the well established theory of Dynamical Systems (Grifiths et al., Sanz-Serna and Vadillo, 1987). However, when we have to find numerical solutions of a Hamiltonian set of differential equations, the standard integrators, no matter how well established, do not have any information regarding the natural geometry of the Hamiltonian system. Usually, only the conservation of energy is monitored, and there are techniques to force the numerical solution to remain on a selected energy surface.

At the end of the eighties, the first ideas on symplectic integrators were launched, based on the work of De Vogleare, who pioneered the subject in 1953. A pioneering Symplectic Integrator Algorithm (SIA), due to J. Wisdom, appeared in the early eighties, and it was applied to Celestial Mechanics with good results. He obtained the symplectic scheme by means of periodic Dirac functions. The periodic deltas were used by Wisdom to replace the high frequency time dependent perturbation, but as a matter of fact he obtained a first order SIA. The boom of these schemes reached its maximum at the beginning of the nineties, and culminated with a meeting held at the Fields Institute, Canada, at the end of 1993. The SIA's can be obtained by ''symplectification'' of some existing algorithms, generating functions of canonical transformations (a natural symplectic integrator), and some researchers started even at a more basic level, the variational principle itself (Wu 1990, Maclahan and Scovel 1993).

The results showed that the SIA's have generally a better performance in qualitative long term investigations, although the energy is preserved only on an average. There are interesting examples that point towards the necessity of preserving the symplectic area or angular momentum rather than the energy which will be discussed later. The main weakness is their poor performance with adaptative step schemes. Nevertheless, their performance on preserving adiabatic invariants has not been excelled by any other kind of algorithm [Benettin et al, l994, Shimada and Yoshida, l996]. Benetin studying models of classical gases uses explicitly the interpretation that when one numerically integrates the system of Hamilton's equation, using a symplectic algorithm, one replaces the "true" Hamiltonian H by a different Hamiltonian K, -close to it, and then iterates "exactly" (within the round-off error of the machine) the time-one mapping of K. Some builders of accelerators (Forest and Ruth, 1990) had used SIA's as a tool long before the subject became more popular. The electron ring is a nice example of a system that, in spite of losses by radiation, has a Hamiltonian core, as already mentioned.

I do not mean to cover all the richness of numerical integrators, their accuracy and stability, but only mention the main aspects bearing on the discussion of the SIA's. I give some examples on the main features using a Yoshida-Ruth SIA of 4th order and a symplectic midpoint RK2. Included is a fairly complete set of references with some comments. To make the text self-contained, I start with a note on Hamiltonian systems and their invariants. The group aspect of the dynamics will also be briefly discussed since it is used to derive higher order SIA's from lower order ones.

There is an alternative to symplectic integrators, which are the exact energy-momentum conserving algorithms, designed to preserve these constants of motion. There is a theorem by Zhong Ge and Marsden (l988) which says that, unfortunately, there cannot exist integration schemes which are both symplectic and energy conserving for non-integrable systems: the only one to preserve both is the flow of the system itself. Therefore, we always have to pay a price for non-integrable systems: a kind of numerical uncertainty principle.

II Hamiltonian Systems: invariants and canonical transformations

A Hamiltonian system is a set of 2n equations given by:

where H(p,q) is the Hamiltonian function and (p,q) = (qi₁,¼,qi_n,pi₁,¼,pi_n) denotes the n canonically conjugated pair of generalized positions and momenta. The flow of a Hamiltonian system preserves volume according to Liouville's theorem. However, there are systems that preserve phase space volume and are not Hamiltonian. This is due to the fact that volume preservation is only one of the many possible invariants that can be derived from the conservation of the symplectic area which gives to the Hamiltonian phase space its typical geometric structure. More precisely, the flow of a Hamiltonian system with n degrees of freedom preserves the sum of the projected areas on the n (p_i,q_i), i = 1¼n, planes. This structure, called symplectic structure, defined by the non-degenerate closed 2-form:

also known as the wedge product. Note that the 2-form is the differential of the reduced action dS = pdq. The wedge product of two arbitrary vectors u, v of the phase space is defined by the operation:

where

0 and the unity matrix I are n×n blocks. J defines the standard canonical structure and is known as the standard symplectic matrix.

If one integrates the k¢th order invariant

^2k, over an arbitrary 2k-dimensional region of the 2n-dimensional phase space, 1 £ i_k £ n, one obtains the invariant

on the 2k-dimensional subspace of the phase space. Making k = n this becomes Liouville's Theorem on the preservation of phase space volume. The various powers of ² are known as the Poincaré invariants. An even more valuable invariant can be written if time and minus the Hamiltonian itself are taken as the (n+1)th pair, (q_n+1,p_n+1); taking k = 1 and applying Stoke's theorem:

this equality for all circuits taken on a tube of trajectories is known as the Poincaré-Cartan theorem. If dt = 0, the circuits are at constant t and the left hand side reproduces (5) in the particular case of k = 1. The form pdq-H dt is known as the integral invariant of Poincaré-Cartan . From the right hand side of equation (6) one can derive all the generating functions S(p, q,t) and rules for obtaining all canonical transformations, and in particular the ones found in most textbooks. The interested reader can find more details in Arnold (l978), Abraham and Marsden (l978), Marsden and Ratiu (l994) and Goldstein (l980).

In terms of the standard symplectic matrix, a canonical transformation of variables : (p,q)® (P,Q) is such that

(Ñ_(p,q)) being the Jacobian matrix of and T stands for the transposing operation.

In Hamilton-Jacobi theory it is necessary to find a canonical time-dependent transformation which takes a set of variables (p_o,q_o) at t = 0, to another set (p,q) at time t. In fact, this is one of the ways of getting a symplectic integrator: we can go forward from t_o to t, and, more, we do so canonically. This map acting for a time interval h = t-t₀, the path of the integrator, is usually called ''the step'' of the integrator (symplectic or not) which will be denoted M(t), then:

The integrating step M(t) is then naturally a symplectic mapping, that is

in accordance with (7). The above symplectic condition should be verified to the same order as the precision of the symplectic operator M(t). In this way, it can be guaranteed that the numerical solution given after k steps of M(t), i.e., for t_k = kt, has the same qualitative behavior of the exact solution: they should agree with each other on the Poincaré invariants up to the same order of the integrator.

Recall that for any pair of functions on the phase space, one can define a bilinear operation {,}, the Poisson bracket, by

Then define z = (p,q) so that the Hamiltonian system can be written in the compact form

It is easy to solve the above equation up to first order, taking the right hand side evaluated at the initial condition

This first order approximation is known as the Euler step for a first order numerical integrator, and if it is substituted back in the original equation and integrated again, one obtains

Continuing this procedure, the formal solution is

called the Hamiltonian flow g which, for an arbitrary H, is a one parameter Abelian Lie group, or . In a more convenient notation (Dragt (l993, 1996), Forest and Ruth (l990)) to be used later in this survey we write z = exp : -t H: z_o, and the Hamilton's equations read

where :-H: z = {-H(z),z}.

I recall also that any function of the phase space that is invariant under the action of this flow is called a first integral or a constant of movement, i.e.,

which is equivalent to {f,H} = 0. On the other hand, any function f(z) which is not constant obeys the differential equation:

Suppose now that there is a symplectic mapping M(t) such that f(z(t)) = M f(z_o), then we have the following sequence of equalities :

This last equation is the invariance of the Poisson bracket under symplectic transformation. However, since nothing is supposed about f(z_o), the equation is true for any arbitrary function of z_o:

which means that M(t) acts on functions of z_o and transforms their functional form in accordance with equation (15). This last equation integrated in the case of an autonomous Hamiltonian function is:

which is the required relation between the symplectic mapping and the Hamiltonian flow. So, symplectic integration means expressing M(t) by a product of mappings which approximates it up to a given order, and in such a way that each factor can be explicitly evaluated on functions of z₀.

The Baker-Campbell-Hausdorff formula (BCH) is a rule which is useful to express this product in terms of the factors. This formula (Dragt l993, l996) organizes the exponents of two elements of the Lie group given in the form e

^A and e

^B where the exponents are some linear operator. The product e^C is given by e^C = e

^Ae

^B, where C is

[A,B] = AÑB-BÑA, the commutator of two vector fields which can be proven to be {H_A,H_B} if A and B stem from Hamiltonians H_A and H_B.

The above formulas do not depend on the symplectic structure and are valid for any Lie Group. More precisely, the BCH formula depends on the structure of the Lie group of the symplectic group, that is, the possibility of writing mappings near the identity as the exponent of some Poisson bracket operator. For example, if we are given two symplectic mappings, M₁(t) = exp(-:tH₁:) and M₂(t) = exp(-:tH₂:), then

i.e., g+1+H₂ = exp (-:tH:) = exp(-: t(H₁+H₂):) up to first order only.

To end this brief review of Hamiltonian systems, it is important to add that the class of integrators related to the conservation laws given by equations (16) are known as reduced algorithms and the main concern is to preserve a chosen first integral. Usually, the classical integration algorithms have some built-in routine constructed to take care of the reduction H = constant at most. Otherwise, the most popular numerical algorithms to solve ordinary differential equations do not ''know'' that they have to preserve either the Poincaré invariants or first integrals other than the energy. Marsden and Ge (1990) have reported an interesting example performed by J. Simó (1990). In this example Simó made simulations of the free vibrations and rotations of a rod using an energy preserving algorithm. The numerical experiment shows that rotations cease after a few periods, violating the angular momentum conservation. So, energy conservation is not always the invariant that must be preserved at any cost. Simó and Wong (1989) have a reduced SIA that also preserves the angular moment first integral. For these I give some references but do not discuss them in this survey. However, it was proved by Marsden and Ge (1988) that the only ''integrator'' which preserves all invariants is the true solution itself. The review of Ge (l996) in the book Numerical Methods in Mechanics has many interesting references and is a nice short introduction.

Note again that the class of invariants to be preserved in a particular situation remains a matter of choice depending on what is being investigated.

III SIA's from generating functions

I follow here Channel and Scovel (1990) and refer to the works of Feng (1986), Feng and Meng-zhao (l988,l991), Feng et al. (1989), Feng and Wang (1994) and Zhong Ge (1996) for a more general presentation, which includes generalizations of the standard symplectic form.

A time dependent generating function of any type (see Arnold (1978), Goldstein (1980)), for example, of the first kind, F(q_o,q,t), generates a canonical transformation from the initial conditions (p_o, q_o) at t = 0 to (p, q) at any given time t, the flow g for some H(q,p). The corresponding transformation equations are

The Hamiltonian flow and F(q_o,q,t) must satisfy the Hamilton-Jacobi equation

The interesting aspect to be noticed here is that equations (22) are a natural integrator, and, better than this, it is canonical or symplectic. This is one instance of the step M(t) defined above. Of course, it is semi-implicit, since one has to solve the first set of equations for q_o, and then the other half is explicit. This means that it is implicit for the q¢s and explicit for the p¢s. However, implicit codes are not unusual. The most well known are the implicit Runge-Kuttas which were developed for stiff systems of ordinary equations.

In general, any mapping derived from a generating function is symplectic. It remains to be found a way of writing the generating function from the Hamiltonian or Hamiltonian's equations so that F(q_o,q,t) is a computable formula as we have seen in the previous section.

The formal series approach supposes the generating function expressed as a series on the small parameter dt = t-t₀:

where G_o = -p_oq which generates the identity transformation. This choice of type of generating function for G is not mandatory. They are comfortable because the identity transformation is the limit of (q,p) ® (q_o,p_o) as t ® t₀. Although generating functions of the first kind do not have this property, they have good symmetrical properties which are useful to construct symplectic integrators that also preserve angular momentum.

In any case, there are many techniques to find out the unknown G_m: for example, with the aid of Hamilton's equation and the formal expansion for (22) up to an arbitrary order m. I quote a few terms taken from Channel and Scovel (1990):

Note that at each order k one has to calculate up to the (k-1)th order derivative of the Hamiltonian. Channel gives a routine that can be easily adapted to Maple. One can also build a specialist code to evaluate G with numerical coefficients. They can be implemented in any language and are faster then any software which is designed for several tasks. See for example Stuchi (2002). It is worth noting that the latest Maple-64 has the advantage of not being limited by its memory. The implementation of this type of SIA is costly anyway: an implicit system has to be solved by Newton's method (one more derivative of the Hamiltonian, totalling k per step of a kth order SIA) which requires the evaluation of matrix inverses. In the next section I discuss explicit first order methods which are easily implemented and more elegant. Moreover, they are the basic mapping, or step, from which higher order explicit SIA's can be built.

IV First order SIA's

I start discussing the easiest case and will use just one variable to avoid cumbersome expressions, so that the main idea can be grasped before more involved discussions. A first order explicit SIA can be generated if the Hamiltonian is of the type H(p,q) = T(p)+V(q). The following first order SIA

where V_q = and T_p = , has as its generating function:

which is the first term in the expansion of Channel and Scovel (1990) given by equation ( 27) of the previous section. Due to the special form imposed on the Hamiltonian, the above SIA1 is explicit in both canonical variables. Note that if one takes the ratios dp/t and dq/t, at each step, and let t go to zero, one recovers Hamilton's equations. It is important to maintain the order in which the step for q_i+1 and p_i+1 are performed, i.e., first the dp increment during time t is calculated at constant q_i, followed by the position change dq at constant velocity, for the same time interval t, as if the potential were turned on only at the position q_i. If a third kind of generating function is chosen, G(p_i,q_i+1), the order of the steps is inverted.

Recall that the Euler integrator for Hamilton's equation derived from H(p,q)=T(p)+V(q), given by

differs from (26) because Dq is evaluated at p_i+1 given by the first stage of the integrator. In this sense it is a leapfrog type method. The Jacobian of the Euler transformation (q_i,p_i) ® (q_i+1,p_i+1) applied for example to the Hamiltonian equations of the simple harmonic oscillator, = p and = -q is

which shows clearly that the mapping is not area preserving, and the value of the Hamiltonian has a secular growth. On the other hand the SIA1 applied to the same equations give:

which is obviously area preserving and a substitution in the Hamiltonian function shows that the energy has no secular growth.

From the above, it is perhaps worthwhile to recall that any integration algorithm is derived by matching the Taylor series development of the true solution. The order k of any integrator means that it agrees with the Taylor series up to order k. There are integrations which are performed directly with the Taylor series written as derivatives of the right hand side of the system, usually called force function in the literature of integrators. This procedure, although explicit, is in general expensive because like the Channel and Scovel (1990) implicit SIA's one has to calculate a large number of derivatives.

As a preliminary illustration, the explicit first order SIA (SIA1) given by equations (26) will be applied to the simple pendulum and to the Hénon-Heiles Hamiltonian which describes the motion of a star in an axial symmetrical galaxy.

IV.1 Some preliminary numerical examples

The first example is the pendulum Hamiltonian

Since the region near the origin is very well approximated by the first order approximation, near the two equilibrium points q_e = ±2kp and q_i = ±(2k+1)p, it can be used to check the behavior near these points which are integrable stable and unstable equilibria, respectively.

The left column of Fig. 1 shows, from top to bottom, the phase space, DE/E and the symplectic area evolution. They have been generated with a forth order Runge-Kutta (RK4) with fixed path h = 0.1, running 120000 steps. The right column shows the same items of the left one, obtained with a symmetric SIA2, a second order symplectic integrator composed with three SIA1 steps to be discussed later. It has been run 350000 steps with path h = 0.3. In spite of the larger path and number of steps, note that the SIA2 represents the unstable fixed point ±p and its invariant manifolds accurately, while the RK4 produces a stochastic layer typical of non-integrable systems. Note also the secular growth of the energy relative error and the decrease of the symplectic area given by the RK4. In contrast, the second figure of the left column shows an oscillation for the energy relative error. This oscillation is a typical feature of symplectic integrators: no matter how long the integration is run, the energy evolution has no secular component. The third figure of each column shows the evolution of the symplectic area. The area dissipation and the energy increase of the trajectories generated by the RK4 depend on the initial condition.

Figure 2 (a) shows the evolution of the symplectic area for the RK4 with four different integrating steps when run in single precision. Note that a smaller path does not improve the symplectic area conservation of the RK4, although the secular growth of the energy error is somewhat improved.

It is even more amazing to compare the performance of a SIA1 with a RK4. In Fig. 3(a) the phase space of the pendulum Hamiltonian has been run with a RK4 with path h = 0.2, while Fig. 3(b) shows the same with a SIA1. As before, the RK4 spreads the phase trajectories as if a dissipation were added, while the SIA1 preserves all the invariant manifolds. Note the interesting break of the symmetry: the trajectories of Fig. 3 (b) are distorted clockwise; this distortion increases with path of integration and appears from a certain path onwards.

The Hénon-Heiles Hamiltonian (Hénon and Heiles, l964)

is a well known system motivated by its astronomical applications. The potential represents the gravitational attraction seen by a star in a galaxy with cylindrical symmetry. Since it has two degrees of freedom the complete three dimensional constant energy trajectories cannot be seen as in the case of the pendulum. Therefore, it is convenient to make a Poincaré section ( Hénon and Heiles (l964), Ozório de Almeida (l987)) through q₁ = 0, ₁ > 0, and project into the (q₂,p₂) plane. As it is well known, this section takes an invariant surface with the topology of a 2-torus of integrable regions, which inhabits the constant energy surface, to invariant closed curves where the iterates of the section lies.

Figure 4(a) shows this section for E = 0.08333333 and path h = 0.1 calculated a with a RK4 running 100000 steps. Fig. 4(b) shows the same section with a SIA1. Note that the SIA1 can compete with the RK4 algorithm for qualitative studies. It can do even better since the SIA1 catches a trajectory with is on a torus very near the separatrix, while the RK4 provokes a separatrix crossing due to its dissipation of area. Compare these with the effect of the RK4 integration in the pendulum case. Fig. 5(a) and (b) show the energy loss and oscillations typical of the RK4 and SIA's, respectively. In Fig. 6 (a) and 6(b) the behavior of the symplectic area is shown, respectively, for the RK4 and the SIA1.

Before explaining the various ways of getting higher order SIA's, I present another interesting way of obtaining a first order one which was devised by Chirikov (1979).

IV.2 Wisdom-Chirikov Delta function technique

Chirikov (l979) proposed a way to solve the problem of a perturbed pendulum Hamiltonian to study what he called the overlap of resonances. The pendulum is in fact a model for one resonance cell or islands that can be seen in the Poincaré sections of two degrees of freedom systems.

A non-perturbed system in action-angle variables is:

whose solution is: q =

t+q₀, I = constant and H(I,q) = H₀(I). If a perturbing potential V(I,q,n t) is introduced, where n is an external or internal perturbing frequency, the system takes the form:

T = is the period of the perturbation, || < 1; since V is periodic in q and nt the perturbing potential can be expanded, for example, in a cosine Fourier series :

If mn+n(I) = 0, for any pair of integers (m,n) we have a resonance, or a so called long period term. Applying perturbation theory, for example Von Zeipel's method, we keep the resonant terms in the ''averaged'' Hamiltonian and the non-resonant ones goes to the generating function of infinitesimal transformations (Bocalletati and Pucacco, 1999; Chirikov, 1979). The averaged or resonant Hamiltonian is then the old one without the non-resonant (short period) terms:

where J = I-I_r, y = y_mn = m nt + nq (the resonant argument) , I_r is the value of I at the exact resonance and ¢ = I = I_r. Thus, we note that hamiltonian (35) represents a pendulum with mass (n¢)^–1 and g = n v_mn(I_r). The following further simplified version of H_R multiplied by a Delta periodic function is the mapping Hamiltonian:

and the periodic Delta function means the addition of a complete series of high frequency terms. In this way, the following time dependent Hamiltonian is obtained, using Poisson transformation (Ozório de Almeida (l987)):

This perturbation is in fact a sequence of Delta functions as follows:

i.e., an expansion in Fourier series of a sequence of delta pulses, which is a periodic delta function of period T_M = 2p/W_M. The equations for H_M(t) are easy to integrate since the change dJ occurs when t = nT_M only, else the system rotates as a free rotor without gravity. The result is the well known standard mapping:

Note that this mapping is similar to the first order one we have discussed in the last two sections. it is easy to see that the determinant of its Jacobian matrix is one, therefore, area-preserving.

Usually the Hamiltonians of resonant problems in Celestial Mechanics are obtained from the classical expansion, which is a Fourier series. This expansion can be separated by the time scale of the arguments of the sines and cosines,

where H_orbital is in general averaged out since it does not affect the long term behavior. Wisdom's method is to write the mapping Hamiltonian as follows:

where d_2p is a sequence of delta functions of period 2p whose Fourier Series is as before with W_M = n being the frequency of the delta's (in the asteroidal case, n is the mean frequency of Jupiter). In general, the time dependence in Hamiltonians (39) is not explicit as in the potential (34).

Some results of this technique applied to the asteroidal resonances 2/1, 3/2 and 4/3, the so called (p+1)/p, can be seen in Murray (l986), Stuchi and Sessin (l989). Stuchi (l991) shows that a mapping derived with this technique and a first order one derived from a generating function give the same result when applied to a sample of real asteroids from the above mentioned resonances. In 1990 Wisdom developed a higher order SIA for the n-body problem and compared with direct integrations performed in his digital orrery: a specially designed parallel computing device to simulate the Solar system (Sussman and Wisdom, l992). For applications in Celestial Mechanics see also Yoshida and Kinoshita (1991).

V The various higher order compositions

Having discussed the Channel and Scovel (l990) integrator up to first order, I now discuss how they can be used to obtain higher order SIA's. I note again that higher order SIA's obtained directly with this procedure are very complicated, and demand hard algebraic manipulations that generate quite cumbersome expressions for the integrator. Fortunately, Forest and Ruth (1983, 1990) and Yoshida (1990) found that higher order SIA's could be generated by performing several intermediate steps during the fundamental step t, in a fashion similar to the one used in generating high order Runge-Kutta integrators. The work of Candy and Rozmus (1991) contains very interesting examples. Unless strictly necessary, only one degree of freedom will be used to keep the notation simple and outline the reasoning.

Take the first order SIA as the map M₁(t) given by

i = 0¼N, where N is the number of intermediate stages necessary to complete one step h. The c_i and d_i are to be determined (less then one by definition). This step can be seen as a sequence of canonical transformations, and through its inverse sequence:

the final Hamiltonian function H_o(q_o,p_o;c₁ ¼ c_N,d₁ ¼ d_N) = 0, up to some order k+1, because the equations of motion are = 0 and = 0. This allow us to determine the values of the c_i and d_i. In fact, we are proceeding as in the resolution of Hamilton-Jacobi equation. For the simple type of Hamiltonian we are considering, H=T(p)+V(q), a second order SIA is obtained for the following set of values:

The idea is simple, but again the complexity of the system of equations for the N coefficients c_i and d_i grows dramatically with the order, and the aid of algebraic manipulator is necessarily already at order four. Fortunately, there is a way out through the use of the Lie group view of de Hamiltonian dynamics which also conveys some elegance to the basic idea. The first researcher who used Lie groups was Neri (1987) and was shortly followed by Forest (1987), Forest and Ruth (1990) and Yoshida (1990), who noticed that the SIA4 previously developed by Forest was a composition of two SIA2.

We have seen how to compose a first order maps for simple Hamiltonians of the type H = T+V. This can also be extended to a more general Hamiltonian of the form H = H₁+H₂. However, in both cases the first order maps, composed of two stages, one for each parcel of the Hamiltonian, do not take into account the fact that these two parcels do not commute. Recently, Laskar and Robutel (2001) have developed a SIA for this type of Hamitonians but with a parameter , that is, H = H₁+ H₂. Their SIA is actually being used for sophisticated investigations of the solar system.

Recall that the Hamiltonian flow, g^t, is a one parameter group of operators acting on phase space. For a specific t, let us call M(t) = g^t the one parameter operator exp (t :-H:) and that ''symplectic integration is the substitution of M(t) by a composition of maps which approximates M(t) up to a given arbitrary order, and this composition may be applied to z_o in an exact and explicit way''.

Note that M(t) = exp:T+V: is M(t) = N₁ (t) N₂(t) = exp t:T:exp t : V: up to first order only and that it can be corrected to higher orders by means of the BCH formula given by equation (21) in the section II. The SIA, N(t), defined by the coefficients (43)

is a symmetrical composition of N₁, N₂ given by equations (41). Applying the BCH formula to (44), it can be shown that

which means that a second order SIA (SIA2) is generated. It can be shown that such symmetric compositions generate SIA's without odd powers of the step t, and these integrators are always of even order. The important consequence is that they are time reversible, i.e.,

I is the identity matrix.

Using this fact Yoshida (1990) developed the theory of SIA's of order 2m+2 by composing maps of order 2m symmetrically. More precisely, given the map N_2m (t) its symmetric composition yields N_2m+2(t), i.e.,

Supposing N_2m known, Yoshida showed with the BCH-formula that:

For this to have order 2m+2 it is imposed that

which when solved for z_o and z₁ gives:

Forest (1991) noticed that the form of the operator R(t = 0) is not important, the important consequence being that we are not limited to the class of Hamiltonians of the form H= T(p) + V(q), or more generally, H = H₁+H₂, as dealt with in Forest and Ruth (l990). Therefore, second order standard integrators which are proven symplectic can be used as a seed for higher order symmetric compositions, for example the Euler step evaluated at the midpoint in phase-space

which is a second order reversible SIA as it will be seen in Section 6. In fact, this approximation was known already to Poincaré, and according to Feng it goes back to Von Zeipel. It has been rediscovered recently in many studies on symplectic integrators. Before closing this section, I again recommend the work of Feng et al (1991) for another view of the midpoint and symplectic algorithms.

V.1 Examples with a symmetric SIA4

First, an example studied by Candy and Rozmus (l991),

which describes the motion of a particle of mass m, charge –e, in the field of a standing wave. In units such that = k = m = 1, the first order Hamiltonian equations of motion can be derived from the Hamiltonian function:

where is the strength of the time dependent perturbation.

In this case the time can be considered as the second coordinate, q₂, and the work is done in the extended phase space so that the potential becomes a function of the q_i's alone. Then, one can make a second order midpoint RK2 to generate the symmetric second order "seed" of the higher order SIA's. Figs. 7(a) and 7(b), = 0.73/2p, show the Poincaré section q = 2kp, which can be obtained easily by taking the step as h = 2p/N and integrating N steps. For h = 2p/25 the SIA4 Fig. 7(b) shows the Poincaré section without any dissipation, while the RK4 with the same step shows chaos and separatrix crossing. Fig. 7(c) shows the energy evolution typical of fixed path RK4 codes as discussed in the example of the pendulum. Fig. 7(d) shows the usual energy behavior, that is, it oscillates around a kind of average value without secular growth, no matter how long the system is integrated.

The Poincaré section shown in Fig. 8(a) was calculated with the RK4, with path h = p/25. We doubled the path used in Fig. 7 and this seems sufficient to keep the trajectory in a invariant torus as in Fig. 7(b)(obtained with the SIA4). However, if we look at evolution of the symplectic area, Figs. 8(b) and 8(c), we see a secular decrease and a sudden series of oscillations and it finally diverges. We do not show the time interval where this divergence occurs. The energy evolution is shown in Figure 8(b) and it shows the usual secular behavior.

In Fig. 9(a) the Poincaré map obtained with the SIA4 reveals an invariant KAM torus which persists for the 30000 Poincaré iterations. This is an evidence of the good stability properties of the symplectic algorithm. The same initial conditions were run with the RK4 and the numerical dissipation due to the non-conservation of the symplectic area makes the iterates migrate to a neighboring invariant curve. The initial conditions are q₁ = 0.0, p = 0.3125, q₂ = t = 0.0, = 1/2p and the path h = 2p/30. In Figs. 10(a) and 10(b) a chaotic trajectory near the separatrix is shown, and in this case the violation of the symplectic area is not so easily identified since the orbit is in a region of widespread chaos due to the migration of iterates from the neighborhood of one set of islands to another, that is, due to the overlap of resonances.

The next example is the classical Hénon-Heiles Hamiltonian (24), and the section used is again the classical one, viz., {x = 0, > 0}. (I have used Hénon's trick to make the section: when the flow cuts the section the integration variable is changed from t to x and a single path is integrated taking h = -x either before or after the section.) Fig. 11(a) shows a regular section obtained with a path h = 0.2, initial conditions (0,0,p₁,0) for E = 0.159 and 100000 steps integrated with the SIA4. The invariant tori are not destroyed even with such large a path. When the energy is slightly increased to E = 0.1597, and the same initial condition (0,0,p₁,0), for the same number of steps, a fuzzy region starts as it can be see in Fig. 11(b). However, when the number of steps is increased to 500000, the region of remanent invariant tori is not invaded by the chaotic iterates as seen in Fig. 11(c). Fig. 11(e) shows the poor, in fact inacceptable, performance of the RK4 for the same conditions of Fig. 11(a). In another trial, I have run 1200000 steps with h = 1/6 and the separate regions of chaos and integrability are still preserved. Fig. 11-d shows the evolution with time of DE; note the typical oscillations of symplectic integrators.

Figure 12(a) shows the Poincaré section obtained during 1200000 steps of h = 1/6 for E = 0.029952 and the initial condition (0,0,p₁,0) using a RK4. As now expected we see the spread of area dissipation. In Fig. 12(b) the poor energy performance can be seen once more. Fig. 12(d) shows the preservation of the symplectic area in the case 12(c) where a sharp Poincaré section can be seen for the same conditions and energy of the RK4 run.

Figure 13 shows two initial conditions near (0,0,p₁,0) for E = 0.159, h = 0.05 for 200000 steps integrated with the SIA4. Note that the tiny islands shown in Fig. 13(a) are finally depicted in Fig. 13(b) and we can see four satellites, in the region near (.001,0,p₁,0), made of even more minute islands. The RK4 run of the same conditions shows that the minute islands are not so sharply captured.

VI Symplectification of standard integrators

An alternative approach to the two main lines presented so far, is the one which tries to impose the preservation of the symplectic area on existing integrators. In this way the extensive set of results and experiments gained on these integrators can be incorporated in their symplectic version, questions like stability, stiffness, etc.

As it is well known the integrating algorithms for ordinary differential equations are divided into two main groups: single path and linear multipath. I do not know of any attempt made to adapt symplectification to extrapolation methods.

VI.1 Symplectic single path integrators

The most used class of single path integrators is the Runge-Kutta family. The leading group in the symplectification of these integrators is that of Sanz-Serna (1988, 1991). There are also the works of Lasagni (1988) and Suris (1989). The three authors independently found a necessary and sufficient condition to make a Runge-Kutta integrator be symplectic. Sanz-Serna and Valdillo already in 1987 started by trying to study the stability of standard integrators using KAM theory. To achieve this, they duplicated the phase space of a non-canonical set of ordinary differential equations and observed that the iterates in this doubled phase space behave like a Hamiltonian system. This procedure leads naturally to the symplectification of standard algorithms. As it is well known, a Runge-Kutta of m stages is given by:

where c_i is the partition of the step h, A is a square matrix of dimension m and b_m are the weight of the various intermediate solutions in the final solution. This is represented schematically by:

These parameters should satisfy the consistency condition c_i = = 1A_ij and = 1 b_i = 1. If A is left triangular, the algorithm is explicit, otherwise it is implicit. The latter is costly and it is used for stiff systems only.

The necessary and sufficient condition for a symplectic Runge-Kutta found by Sanz-Serna, Lazagni and Suris is that the m-dimensional matrix M whose elements are

be zero. In fact, this matrix was established in studies concerning algebraic stability of Runge-Kuttas' algorithms. It can be proved that all RK which are of the type Gauss-Legendre (Decker and Verver, 1984) satisfy condition ( 55). However, no explicit RK satisfy this condition, while the simply-implicit methods can be made to obey condition ( 55) (Cooper, 1987). Sanz-Serna (1991) is entirely dedicated to this question, using tree theory.

The simplest RK which is symplectic is the midpoint, one stage, simply implicit Runge-Kutta:

where c₁ = , A₁₁ = and b₁ = 1.

Figure 14 shows the performance of the midpoint RK2 compared to a RK4 used so far in all experiments in this review. The same conditions of Fig. 9 and some nearby trajectories were integrated with the midpoint RK2 and are shown in Fig. 14 (a). The neatness of the midpoint compared with the RK4, Fig, 14 (b) is quite impressive. Note that the thin stochastic layer near the separatrix of the big island does not invade invariant curves as it does with the RK4.

Lazagni found a generating function for the canonical transformation defined by the symplectic RK's, and it is given by the following expression:

where Y_i (P,Q)^T are the stages of the RK, while H_p, H_q are column vectors with the partial derivatives of H and should be seen as functions of (p_n,q_n+1) and h.

Sanz-Serna (1988) found that RK's preserve quadratic first integrals, or more generally that bilinear invariants are preserved for any system of ordinary differential equations. His interesting result is that imposing energy conservation for general Hamiltonian functions leads to degradation of the performance of the algorithm. However, symplectic conservation does not do any harm to the algorithm. This result is once more in accordance to the theorem by Marsden and Zhong Ge (1988), stating that for non-integrable Hamiltonians only the true solution is able to preserve all first integrals besides the symplectic form. As it can be seen from the derivation of the Yoshida-Ruth algorithm, the energy is preserved only up to the order of the discrete symplectic map S(t). Therefore, one can not in principle take care of all invariants at the same time, and for numerical integration one has always to pay a price.

Suris (1989) has results that are analogous to these of Sanz-Serna and Lazagni. He finds symplectic conditions for RK-Nystron which are algorithms which work directly with second order differential equations. They are convenient for cases where obtaining Hamilton's equations is not easy. These integrators have the following form:

and they apply to systems of the form

Note that this RK differs from the usual ones by the term in h² in the evaluation of the intermediate stages, but the final formula is of the same type. The symplectic condition is very similar to condition (55):

VI.2 Multi-step linear algorithms

The multi-step algorithms are the least developed case of symplectic integration. A linear multi-step formula requires the k previous paths for the evaluation of y_n+k and its general form is:

where f_r = f(t_r,y_r). If the last coefficient b_k = 0 then the formula is explicit since the left hand side does not involve the value of the point y_n+k which is to be evaluated; otherwise the multi-step is called implicit. The method has two characteristic polynomials

and the condition s(t) = 1 is used to normalize (61).

The simplest example of this class of formulas is the trapezoidal method given by

which is an implicit multi-step algorithm which depends on two points; the values of its parameters according to equation (61) are: b_o = b₁ = 1/2 and |a_i| = 1. The trapezoidal method involves two evaluations of the force function f(y,t) at the beginning and at the end of the time step h.

There is a class of multi-step algorithm which are known as backward differentiation (BD) and they have b_i = 0 for i ¹ k. These methods need only one evaluation of the force function. There are also the one-leg multi-path algorithm which are similar to the BD's because they also use only one evaluation of the force function, but in the form of an average of all steps, i.e.,

The easiest example is the ''center midpoint'', or midpoint, a one-leg counterpart of the trapezoidal method which has the following form

In fact, this method can also be seen as the one stage implicit Runge-Kutta given by equation (56).

This midpoint formula has been shown to be symplectic by many authors . For example, Aziu (1985) started with the general definition of a multi-step formula and using the conditions for consistency and stability concluded that they can only be of two steps in order to be symplectic. Again the midpoint algorithm! However, Suris in his paper on symplectic Runge-Kuttas comments that the Aizu result is applicable to formulas which can be reduced to a single step formula, i.e.,

which is in fact a two stages RK given by the matrix

and satisfies condition ( 55).

Eirola and Sanz-Serna (1992) show that linear one-leg multi-step formulas, when symmetrical and irreducible, preserve not only any quadratic invariant but also the symplectic structure. A linear multi-step formula is symmetric if it satisfies the following relations for its parameters

and is said to be irreducible if the polynomials r(s) and s(s), which are characteristic of the method, do not have any common root. Once more the omnipresent midpoint satisfies these conditions.

Therefore, the linear multi-path formulas are not in general symplectic. Unfortunately, the only instance where they can be shown to be symplectic, the symmetrical cases, have a very poor performance as regards stability for more than one leg. This is the issue which can impair efficient implementation of multi-step symplectic formulas since a one-stage formula is too poor. Moreover, if the Hamiltonian system is stiff, linear multi-steps are cheaper than implicit Runge-Kuttas.

VII Final comments, stability and variable step size

The idea of symplectic integrators is theoretically attractive and the first results were most encouraging. The ones shown in this survey were collected in the literature and verified. However, one would like to see them implemented in variable step size like the standard codes that are largely used, but the results were quite disappointing: they perform exactly as the standard codes as far the preservation of the symplectic structure is concerned. This means that preservation of the symplectic structure seems to go with fixed step only.

Gear (1992), a traditional researcher in numerical analysis, constructed an argument to justify the bad performance with variable step size. He finds that the step h_n depends on the previous step, h_n-1, and iterate y_n, and not on the derivatives of the Hamiltonian function. In fact, his results are in accordance with the ones by Calvo and Sanz-Serna (1992). However, according to Gear, their results prove only that a SIA designed for a fixed step performs badly when operating at variable step. Therefore, there is still hope that a SIA can be specially designed to operate at variable step. The only work so far in this direction is due to W. McEvoy (1992), adapting a SIA to operate in stages like the Runge-Kuttas.

Stability is a concept which should not be mixed up with precision. Sometimes, trying to get a better precision, one can be working outside the stability region of a given method. This word of caution is due not only for symplectic algorithms but for any integrating algorithm as well. The general references at the end can supply the reader with some nice examples. Since we have mentioned precision, McLachlan e Atlea (l992) made a very extensive work on precision of SIA's including a highly precise fourth order method for a special class of Hamiltonians.

To conclude these notes, it should be said that the subject is still alive due to its good performance in problems where adiabatic invariants are of interest (see Yoshida, l998). Yoshida carried some numerical experiments with a one-dimensional harmonic oscillator with a slowly varying frequency and a one-particle system with a slowly varying isochronic potential. He found that the adiabatic invariant of these systems when integrated with SIA's, do not show secular growth of the error, unlike the Runge-Kutas. He also found that the best order depends on the error tolerance required.

Since the midpoint scheme is useful to compose higher order symplectic integrators the paper by Acher and Reich (l999) presents a study on advantages and pitfalls of this SIA2. A recent paper by M.Guzzo (2001) improves the midpoint scheme and applies it in a perturbed Kepler problem. Although the paper is focused on Celestial Mechanics, the main ideas can be used in other problems. He takes into account the perturbation introduced by the symplectic integrator and is able to express it in the original Hamiltonian in a separate form. There is another recent development of symplectic integrators by Laskar and Robutel (2001 , 2002) and it is actually being used in their current research on the Solar system.

A revival of the "kick" technique of Chirikov is taken by Abdullaev (l998) treating it in the framework of canonical transformations to obtain rigorous symmetric mappings, thus improving the standard mapping as an integration device for Hamilton's equation. This is also done by Wisdom et al. (1993) and it is instructive to compare the two improvements of Chirikov's trick.

In Dattoli et al. (1997) symplectic integrators methods are carried to unitary operators in quantum mechanics with many interesting examples in optics. Also, the idea of preserving the symplectic geometry inspired people working in non-holonomic mechanics (Barth et al (1999) and Cortéz et al. (2001)). In this context perhaps it should be worth while to mention McLahlan and Scovel (1993) for an introductory paper.

To close, I find appropriate a quote taken from Channel and Neri (l996) due to Hamming (l986):" an algorithm which transforms properly with respect to a class of transformations is more basic than one that does not. In a sense the invariant algorithm attacks the problem and not the particular representation used".

Acknowledgments: I thank Dr. Jair Koiller for calling this subject to my attention and Dr. Carlos Eduardo Aguiar for suggesting the numerical calculation of the symplectic area. I also thank Dr. A.M. Ozório de Almeida for his careful reading of the manuscript and several suggestions that greatly improved this review. Fundação José Bonifácio is acknowledged for providing computational support.

VIII Bibliography

VIII.1 Numerical integrators

VIII.2 Hamiltonian Systems

VIII.3 Symplectic integrators obtained via generating functions and composition of first order SIA's

The best list of references are in the paper of Scovel and in the review article by Yoshida. The papers of the Feng group are more mathematical and general. The Physica D, 60 of l992 has many interesting results so far. The special issue of the conference held at the Fields Institute: Marsden, J.E., Patrick, G.W. and Shadwick, W.F. (eds), (l996) Integration Algorithms and Classical Mechanics, American Mathematical Society, Providence, RI "Numerical Methods in Mechanics" is also a worth reference book.

VIII.4 Symplectification of standard algorithms

VIII.5 Momentum preserving algorithms

• Marsden, J.E., O'Reilly, O.M., Wicking, F.J. and Zombro, B.W., (1990), Symmetry, stability, geometric phases and mechanical integrators, Non-linearity - pre-print.

VIII.6 Variational Approach

VIII.7 Variable step for Symplectic Integrators

• MacEvoy, W., (1992), An adaptative symplectic integrator, Los Alamos, personal communication.

VIII.8 Numerical Stability of Integrators

Besides the general texts mentioned above, its worthwhile to mention the following ones.

VIII.9 Numerical experiments and symplectic integrators accuracy

VIII.10 Recent papers: miscellaneous

• Laskar, J. and Robutel P.(2002), personal communication.

•  Baker, D.T.H., Phillips C., eds., (1981) The numerical solution of non-stiff problems, Oxford Science Publications, Clarendon Press, Oxford.
•  Gear, D.W., (1971) Initial value problems in ordinary differential equations, Academic Press.
•  Gradwell I., Sayers, D.K., eds., (1981) Solving ordinary differential equations I, Nonstiff Problems, Spring-Verlag. (it has a section with systems which have strange attractors).
•  Hamming, R.W., Numerical Methods for Scientists and Engineers, Second edition, p.73, Dover, New York.
•  Arnold, V.I. : (1978) Mathematical Methods of Classical Mechanics, Springer-Verlag, New York, l978.
•  Abrahm, R. and Marsden, J.E. (1978) Foundations of Mechanics, Addison-Wesly, Reading, Mass.
•  Marsden, J.E. and Ratiu, T.S, (1994) Introduction to Mechanics and Symmetry, Springer-Verlag, New York.
•  Ozório de Almeida, A.M., (1987), Sistemas Hamiltonianos: caos e quantizaçăo, Editora UNICAMP, Campinas, SP.
•  Goldstein, H.: (1980) Classical Mechanics, Addison-Wesley, Reading, Mass.
• Candy, J. and Rozmus, W., (1991), J. Computational Physics, 92, 230-256.
• de Vogleaere, (1956) Methods of integration which preserve the contact transformation property of Hamiltonian systems, Tech. Report No 4, Dept. Mathem., Unev. of Notre Dame, Notre Dame, Ind.
• Dragt, A. J. (1993), Lectures on Nonlinear Dynamics and Lie Methods with Applications to Accelerator Physics, University of Maryland Physics Department.
• Dragt, A.J., Abell, D.T., (1996) Symplectic Maps and Computation of Orbits in Particle Accelerators, Fields Institute Communications, Volume 10.
• Chirikov, B.V., (1979), Phys. Rep., 52, 263.
• Channel, P.J. and Scovel, J.C., (1990) Symplectic integration of Hamiltonian systems, Nonlinearity 3, 231-29. (at the end of this article there is a algebraic code which can be translated into Maple and can be used for experiments with your choice Hamiltonian.)
• Forest E., Bergtsson, J., Reusch, M.F., (1990) Application of the Yoshida-Ruth techniques to implicit integration and multi-map explicit integration, Physics Letters A, 158, 99-101.
• Feng, K., (1986), Difference Schemes for Hamiltonian Formalism and symplectic geometry, J. Comp. Math 4, 27-28.
• Feng, K.and Meng-zhao, Q., (1988), The symplectic methods for the computation of Hamiltonian equations, in Numerical Methods for partial differential equations, Lecture Notes in Math. 1297, Springer-Verlag.
• Feng, K., Hua-mo Wu, and Meng-zhao, Q., and Dao-liu Wang, (1989), Construction of canonical difference schemes for Hamiltonian formalism via generating functions, J.Comput.Math. 7, 371-380.
• Feng, K. and Dao-liu, W., (1991), Symplectic difference schemes for Hamiltonian systems in general symplectic structure, J. Comp. Math 9(1), 86-96.
• Feng, K. and Meng-zhao, Q., (1991) Hamiltonian algorithms for Hamiltonian systems and comparative numerical study, Comp. Phys. Communications, 65, 173-187.
• Forest, E., Ruth, R.D., (1990), Fourth-order symplectic structure, Physica D, 443, 105-117.
• Hénon, M. and Heiles, (1964) The applicability of the third integral of motion: some numerical experiments, Astron. J. 69, 73/79.
• Kinoshita, H., Yoshida H., and Nakai, H., (1990) Symplectic integrators and application to dynamical astronomy, Celestial Mech., 50, 59/71.
• Murray, C.D., (1986) Structure of the 2/1 and 3/2 Jovian resonances, Icarus, 65, 70.
• Ruth, R.D, A canonical integration technique, IEEE Trans. Nucl. Sci., 30, 2669-2671.
• Scovel, J.C., (1989), Numerical integration of Hamiltonian systems, Proc. Workshop on Geometry of Hamiltonian Systems, Ratiu, T, ed., Springer, 463-496.
• Stuchi, T.J. and Sessin, W. (1989), The reduced Hamiltonian system as an alternative mapping for the first order resonance, in Orbital Dynamics of Natural and Artificial Objects, R.Vieira Martins, D. Lazzaro and W. Sessin editors, Observatório Nacional, Rio de Janeiro, Brasil.
• Stuchi, T.J., (1991) Doctor Thesis, Instituto Tecnológico de Aeronautica (ITA), Săo José dos Campos, Săo Paulo, Brasil.
• Stuchi, T.J., (2000) Invariant KAM tori in the center manifold of the Hill-3D problem, Advances in Space Mechanics, Proceeding of the X Colóquio Brasileiro e Dinâmica Orbital, Nazaré Paulista, S. Paulo, to appear.
• Sussman, G.J. and Wisdom, J., (1992) Chaotic evolution of the solar system, Science 257, no.5066, 56-62.
• Yoshida, H., (1990), Construction of higher order symplectic integrators, Phys. Lett. A, 463-268.
• Yoshida, H., (1990), Conserved quantities of symplectic integrators for Hamiltonian systems, Physica D.
• Yoshida, H., (1993), Recent progress in the theory and application of symplectic integrators, Celestial Mechanics and Dynamical Astronomy, 56, 27-43.
• Wisdom J., (1982), Astron. J. 87, 577.
• Wisdom J., (1983), Icarus, 56, 51.
• Wisdom J., (1987), Urey Prize Lecture: Chaotic Dynamics in the Solar System, Icarus 72, 241.
• Wisdom J., Holman M., (1991) Symplectic maps for the n-body problem, Astron. J., 102(4), 1528-1538.
• Aizu, K., (1985), Canonical transformation invariance and linear multi-step formula for integration of Hamiltonian systems. J. Comp. Phys. 58, 270-274.
• Eirola, T. and Sanz-Serna, J.M., (1992) Conservation of integrals and symplectic structure in the integration of differential equations by multisteps methods, Numer. Math. 61, 281-290.
• Zhong Ge and Marsden, J.E., (1988), Lie-Poisson Hamilton-Jacobi Theory and Lie-Poisson integrators, Phys. Letters A, 133 (3), 134-139.
• Lazagni, F., (1988) Canonical Runge-Kutta methods, ZAMP 39, 952-953.
• Okunbor, D. and Stell, R.D., (1992), An explicit Runge-Kutta-Nystrom method is canonical if and only if it is adjoint and explicit, SIAM J. Number, Anal., 29 (2), 521-527.
• Sanz-Serna, J.M., (1988), Runge-Kutta schemes for Hamiltonian systems, BIT 28, 877-883.
• Sanz-Serna, J.M. and Abia, L., (1991), Order conditions for canonical Runge-Kutta schemes, SIAM J. Numer. Anal. 28, 1081, 1096.
• Sanz-Serna, J.M., (1992), Symplectic integrators for Hamiltonian problems: an overview, Acta Numerica, 1992, Cambridge Univ.Press, Cambridge, pp. 243-286.
• Sanz-Serna, J.M. and Calvo, M. (1984), Numerical Hamiltonian Problems, Chapman and Hall, London.
• Suris, Y.B., (1989), The canonicity of mappings generated by Runge-Kutta type methods when integrating the systems x˘˘ = -śU/śx, USSR Comput. Maths. Math. Phys., 29 (1), 138-144.
• Zhong Ge and Marsden, J.E., (1988), Lie-Poisson Hamilton-Jacobi Theory and Lie-Poisson integrators, Phys. Letters A, 133 (3), 134-139.
• Marsden, J.E., Lectures on Mechanics, chapter 9 on Mechanical Integrators, London Math Society, Lecture Notes Series 174, 1992, Cambridge Press.
• Simó, J.C. and Wong, K.K., (1989), Unconditionally stable algorithms for the orthogonal group that exactly preserves energy and momentum, Int. J. Num. Engng, 31, 19-52.
• Simó, J.C., N.Tarnow and Wong, K.K., (1992) Exact Energy-Momentum Conserving Algorithms and Symplectic Schemes for Nonlinear Dynamics, Computer Methods in Applied Mechanics and Engineering, 1, 63-116.
• Wu, Y. (1990), The discrete variational approach to the Euler-Lagrange equations, Computers Math Applic., 20 (8), 61-75.
• Mclachlan, R. I. and Scovel, C. (l993), Equivariant Constrained Symplectic Integration, Technical Report, Program in Applied Mathematics, University of Colorado at Boulder.
• Calvo, M.P. and Sanz-Serna, J.M. (1991), Variable steps for symplectic integrators, Applied Mathematics and Computational Reports, report 1991/3.
• Gear, C.W., (1992), Invariants and Numerical Methods for ODE's Physica D, 60, 303-310.
• Skeel, R.D. and Gear, C.W., (1992), Does variable step size ruin a symplectic integrator?, Physica D, 60, 311-313.
• Cartwright, J.H.E. and Piro, O. (1992), The dynamics of Runge-Kutta methods, International Journal of Bifurcation and Chaos, 2, 3, 427-449.
• Cooper, G.J., (1987), Stability of Runge-Kutta methods for trajectory problems, IMA J. Numer. Anal. 7, 1-13.
• Corless, R.M., (1992), Defect-controlled numerical methods and shadowing for chaotic differential equations, Physica D60, 323-344.
• Dekker K. and Verwer, J.G., Stability of Runge-Kutta methods for stiff nonlinear differential equations, North-Holland, 1984.
• Gear, C.W., (1971) Numerical initial value problems in ordinary differential equations, Prentice-Hall, Inc. Englewood Cliffs, N.J.
• Griffiths, D.F., Sweby, P.K. and Yee, H.C., (1992), On spurious asymptotic numerical solutions of explicit Runge-Kutta methods, IMA J. Numer. Anal. 12, 319-338. (this paper uses concepts from Dynamical Systems).
• Sanz-Serna, J.M., (1985), Studies in numerical nonlinear instability I. Why do leapfrog schemes go unstable?, SIAM J. Sci. Stat. Comput., 6 (4), 924-938.
• Sanz-Serna, J.M. and Vadillo, F., (1987), SIAM J. Appl. Math 47 (1), 92-108.
• Wu, Yuhua, (1991), A note on the stability of two-level symplectic schemes, Appl. Math. Lett. 4(5), 1-5.
• Z. Mei-ging and Q. Meng-zhao, (1991) A note on convergence of symplectic schemes for Hamiltonian systems, J. Comp. Math, 9 (1), 1-4.
• Feng, K. and Meng-zhao, Q., (1991), Hamiltonian algorithms for Hamiltonian systems and comparative numerical study, Comp. Phys. Communications, 65, 173-187.
• Mclachlan, R.I. and Atela, P., (1992), The accuracy of symplectic integrators, Nonlinearity 5, 541-562.
• Okunbor, D., Canonical methods for Hamiltonian systems: numerical experiments, Physica D 60, 314-322.
• Skeel, R.D, Zhang, G. and Schlick, (1998), A family of symplectic integrators: stability accuracy, and molecular dynamics applications, SIAM J. Sci.Compot. 18, no. 1, 203-222, Dedicated to C. William Gear on the occasion of his 60th birthday.
• Abdullaev, S.S., (l999), A new integration method of Hamiltonian systems by a symplectic mapping, J.Phys. A, Math.Gen-1999- V 32,2745-66.
• Ascher, U.M. and Reich, S., (l999), The midpoint scheme and variants for Hamiltonian systems: advantages and pitfalls, SIAM J. Sci. COMPUT, Vol. 21, No.3, 1045-1065.
• Chin, S.A., (l997), Phys. Lett. A, 226, 334.
• Chin, S.A., Kidwell, D.W., (2000) Higher-order force gradient symplectic algorithms, PRE, 62, 6, 8746-8752.
• Dattoli, G., Ottaviani, P.L., Torre, A. and Vásquez, A. , Evolution operator equation: integration with algebraic and finite -difference methods. Application to optical problems in classical mechanics and quantum field theory, La Rivista del Nuovo Cimento, l997.
• Benettin, G., Adiabatic invariants and time scales for energy sharing in models of classical gases, in Hamiltonian Mechanics, J.Seimenis editor, Plenum Press, New York, l994.
• Guzzo, M, (2001), "Improved leap-frog symplectic integrator for orbits of small eccentricity in the perturbed Kepler problem", Celestial Mechanics and Dynamical Astronomy 80,63-80.
• Holman, M. and Murray, N. (1999), The origin of chaos in the outer solar system, Science 283, no.5409, l877-1881.
• Laskar, J. and Robutel P.(2001), "High order symplectic integrators for perturbed Hamiltonian systems", Celestial Mechanics and Dynamical Astronomy 80, 39-62.
• Marsden, J.E., Patrick, G.W. and Shadwick, W.F. (eds), (1996) Integration Algorithms and classical mechanics, American Mathematical Society, Providence, RI.
• Marsden, J.E. and Patrick, G.W., (1998) Multisymplectic geometry, variational integrators and nonlinear PDEs, Comm.Math.Phys. 199, no.2, 351-395.
• Channel, P.J. and Neri, F.R., An Introduction to Symplectic Integrators, Fields Institute Communications, Volume 10, 1996.
• Shimada, S. and Yoshida H., Long-ter conservation of adiabatic invariants by using symplectic integrators,Publ. Astron. Soc. Japan 48, 147-155, (1996).
• Skeel, R.D, Zhang, G. and Schlick, (1997), A family of symplectic integrators: stability accuracy, and molecular dynamics applications, SIAM J. Sci.Compot. 18, no. 1, 203-222, Dedicated to C. William Gear on the occasion of his 60th birthday.
• Schilier, Ch. and Seiter, A., (1998) Symplectic integration of classical trajectories: A case study, J.Phys.Chem. A 102, 9399-9404.
• Wisdom, J., Holman, M., and Touma, J. (1996) Symplectic correctors, Integration algorithms and classical mechanics, J.E.Marsden, G.W.Patrick and W.F.Shadwick editors, Fields Institute Communications, 10, 1996.

Publication Dates

• Publication in this collection
  11 Feb 2003
• Date of issue
  Dec 2002

History

• Received
  10 Apr 2002