Time-symmetric integration in astrophysics

Time-symmetric integration in astrophysics

David M. Hernandez and Edmund Bertschinger 11footnotemark: 1
Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Ave.,
Cambridge, Massachusetts 02139, USA
Email: dmhernan@mit.edu (DMH); edbert@mit.edu (EB)

Calculating the long term solution of ordinary differential equations, such as those of the -body problem, is central to understanding a wide range of dynamics in astrophysics, from galaxy formation to planetary chaos. Because generally no analytic solution exists to these equations, researchers rely on numerical methods which are prone to various errors. In an effort to mitigate these errors, powerful symplectic integrators have been employed. But symplectic integrators can be severely limited because they are not compatible with adaptive stepping and thus they have difficulty accommodating changing time and length scales. A promising alternative is time-reversible integration, which can handle adaptive time stepping, but the errors due to time-reversible integration in astrophysics are less understood. The goal of this work is to study analytically and numerically the errors caused by time-reversible integration, with and without adaptive stepping. We derive the modified differential equations of these integrators to perform the error analysis. As an example, we consider the trapezoidal rule, a reversible non-symplectic integrator, and show it gives secular energy error increase for a pendulum problem and for a Hénon–Heiles orbit. We conclude that using reversible integration does not guarantee good energy conservation and that, when possible, use of symplectic integrators is favored. We also show that time-symmetry and time-reversibility are properties that are distinct for an integrator.

methods: numerical—celestial mechanics—globular clusters: general— planets and satellites: dynamical evolution and stability— galaxies

d [1] #1 [1] #1 [1] \s#1 ∂ [1] e^#1 [2] φ_#1^#2 - [1] ~#1 12 ϵ 12 V_pl V_inter ¯V ¯V_n-1 T_pl ~T_pl T_ch H_pl [2]\diff#1\diff#2 [2]\p#1\p#2 [3]\p2#1\p#2\p#3 [2] ()

1 Introduction

Obtaining solutions to initial value problems of ordinary differential equations (ODEs) over long time periods is central to dynamical calculations in astrophysics. These ODEs might represent problems such as the -body problem, point particles interacting through pairwise forces, or the problem of particle orbits in a time-independent galactic potential. The ODEs are frequently described by a time-dependent or time-independent Hamiltonian.

Obtaining a solution to the -body problem is essential for many purposes, from calculating the evolution of dark matter in the Universe to understanding stability and chaos of orbits in planetary systems. Different techniques, relying on different assumptions, have been developed to obtain approximate -body solutions. The -body problem is generally chaotic and non-integrable, so we rely on these approximations to obtain its solutions. But, in general -body cases, it is unknown how reliable such approximations are. In fact, the approximations themselves give rise to chaos, separate from the physical chaos of the problem itself. If the numerical method itself can be responsible for chaos, then the error from the original trajectory can grow exponentially, leading to call into question the validity of the calculated solution.

Galactic potentials usually have only a few degrees of freedom, but can still be chaotic and non-integrable, and suffer from the same problems described above. In fact, much of the study of chaos began with the study of the Hénon–Heiles problem, which was motivated by the study of galactic potentials.

It would appear numerical approximation to chaotic ODEs should be suspect, but fortunately, geometric numerical integration, integration aimed at respecting the geometry of the underlying ODEs, has been developed and helps restore confidence in these numerical solutions (Channell & Scovel, 1990). Depending on the equations, geometric properties include the Hamiltonian flow, time-reversibility, and quadratic and linear invariants in the phase space. In the last 30 years, astrophysics researchers have made geometric integration a standard in various fields of dynamics, including planets (Wisdom & Holman, 1991; Chambers, 1999; Duncan et al., 1998; Hernandez, 2016), stellar clusters (Kokubo et al., 1998; Hut et al., 1995; Hernandez & Bertschinger, 2015; Dehnen & Hernandez, 2017), or galaxy formation (Springel, 2005).

Geometric integrators that respect Hamiltonian flow are also called symplectic integrators, and they conserve generalizations of volumes in phase space, also known as Poincaré invariants. The theory of symplectic integration is well developed (Hairer et al., 2006). The citations above (Wisdom & Holman, 1991; Chambers, 1999; Duncan et al., 1998; Hernandez, 2016; Kokubo et al., 1998; Hut et al., 1995; Hernandez & Bertschinger, 2015; Dehnen & Hernandez, 2017; Springel, 2005), are all concerned with time-independent Hamiltonian problems, so a symplectic integrator is ideal. However, symplectic integrators applied to the above problems have a limitation; if the step sizes are chosen as a function of the phase space, the evolution of the trajectory is no longer Hamiltonian. This limitation is severe for the -body problem because gravity has no length scale: two-body relaxation is affected by close and far encounters. Thus, the range of time and length scales is large, posing a severe challenge for fixed time step integration.

Thus, some researchers (Pelupessy et al., 2012; Hut et al., 1995; Kokubo et al., 1998) have abandoned the requirement of a symplectic integrator and instead focused on integrators that preserve time-reversibility, if the underlying equations have this symmetry. It is important to note there exist irreversible conservative differential equations— see Section 2.2. However, many important problems such as the -body problem are conservative and reversible. Time-reversible integration appears to be less studied than symplectic integration, but it has been observed that time-reversible integrators generally can reduce errors for integrations in astrophysics. An explanation for such behavior is sometimes not provided. It is possible to adapt time steps while still preserving time-reversibility (Makino et al., 2006; Funato et al., 1996), so clearly we would like to abandon symplecticity if possible and if time-reversible integration is good enough. But a clear error analysis is needed with these integrators in order for one to have confidence in their use.

The goal of this paper is to provide that error analysis, and to use it to show that the behavior of a time-reversible integrator in astrophysics can be worse than a symplectic integrator. This suggests that researchers in astrophysics should use symplectic integrators when possible. To perform the error analysis, we derive the modified differential equations (MDEs) obeyed by these methods using adaptive time steps, and we use these equations to calculate how well the methods conserve energy. We study a simple pendulum problem and Hénon–Heiles orbits. We show how various reversible integrators do not conserve energy to all orders. It was already noted by Faou et al. (2004) that some fixed-step Runge–Kutta reversible methods do not conserve an energy.

Section 2.1 shows the tools necessary for deriving the MDE. Section 3 derives the MDE for the trapezoidal rule, a non-symplectic but symmetric second order Runge–Kutta method. We also derive the MDE for the trapezoidal rule with adaptive steps. We derive various properties of the trapezoidal rule. In Section 4, we apply our numerical analysis to understand the error in energy of the modified pendulum problem and Hénon–Heiles orbits and to find that time-symmetric integration can give energy drift. While the analysis in this section is limited to the trapezoidal rule, there is no reason to believe other time symmetric integrators would not suffer from energy drift. In the Appendices, we show that for Runge–Kutta methods, time-symmetry, reversibility, and symplecticity are independent concepts. We conclude in Section 5.

2 The modified differential equation

2.1 Time-symmetric integration

A system of autonomous ordinary differential equations can be written


where and are both vectors of length . We are concerned with the case where is a vector of positions and velocities or the phase space defined by canonical coordinates and momenta: in equations, or . The problem is to find given . We assume that the system is autonomous, so that depends only on and not on .

A numerical one-step method estimates the solution at , , where


for some that is related to . The method may be iterated to estimate the solution at etc. For now, we assume that is a constant independent of and .

A goal of numerical analysis is to find a that is inexpensive to evaluate so that is smaller than a specified tolerance. The study of the errors is known as forward error analysis. One can turn around the problem. Given , find a modified differential equation whose exact solution is . That modified differential equation is written


The goal then becomes to minimize . This is done by determining from . Studying the errors is known as backward error analysis.

A symmetric one-step integrator is one for which a forward step followed by a backward step restores the initial conditions. The requirement is


The associated modified differential equation is even:

-reversibility (Hairer et al., 2006, Section V.1) means that if we change the sign of velocities, while keeping the position coordinates constant, the solution trajectory must stay the same— only the direction of motion is inverted. Let be an invertible linear transformation that changes the signs of velocities: . All autonomous Newtonian physics problems are described by positions and velocities and can be written as a system of first order ODE’s: and . They are not all reversible; if they are, then,


While many problems satisfy this requirement, not all do. For example, the system of differential equations for a charged particle moving in a magnetic field are


where is the charge of the particle, is the mass of the particle, and is the external magnetic field. These equations do not satisfy (5); the solution trajectory is different when we switch the sign of the velocities (the resolution of this apparent irreversibility is that the sign of changes if we reverse the currents causing it).

If we use a one-step method to solve a -reversible set of differential equations, the symmetric integrator is -reversible. The -reversibility condition for an integrator is connected to (5):


which implies . This only holds if the integrator is time-symmetric, or . Thus, in what follows, until Section 3.4, ‘symmetric’ one-step methods will be equivalent to‘time-reversible’ one-step methods because we are only concerned with -reversible differential equations. However, it is important to bear in mind that a symmetric method is not necessarily the same as a time-reversible method; one way to break the equivalency is by letting the step vary as a function of phase space.

2.2 Derivation of modified differential equation

Our goal is to understand time-symmetric integrators. Unlike symplectic methods, symmetric integrators generally have no surrogate Hamiltonian (Hairer et al., 2006, Section IX.8) which informs us of the dynamics, so we instead derive the differential equations the integrator obeys. We call this the modified differential equations (MDEs), and its study has been referred to as backward error analysis (Hairer et al., 2006, Chapter IX).

Proceed as follows: first write the formally exact solution of equation (3) with initial condition ,


This is just the usual Taylor expansion solution. Next expand and in power series in :


is the from (1). Use (9) to expand the derivative operator


Combining equations (2) and (8)–(10) gives, for ,


Our goal is to obtain from . One way is to solve equations (2.2) recursively, starting with substituting into , and so on. This is useful for determining for small .

As an example, consider the explicit Euler method


for which , . This method is first order because . For an th order method,


Recursive solution is impractical to extend to high order. An alternative approach (which may also be difficult, but is conceptually appealing) is to sum the series for in equation (8) by defining the differential operator




The logarithm of an operator is defined by its series expansion. Applying the operators to gives and so that


The relation defines .

3 A study of the trapezoidal rule: a time-symmetric but non-symplectic integrator

3.1 Relating the trapezoidal and midpoint rule

We introduce several one-step integrators. Let , and indicate the trapezoidal, midpoint, explicit Euler, and implicit Euler one-step integration methods, respectively. The midpoint rule is symplectic, while the trapezoidal rule is not, but they have a close connection. The two integrators are defined by




The explicit and implicit Euler methods are first-order, not time-symmetric, and non-symplectic. They are


We see that


so that . Thus, the trapezoidal and midpoint rules are said to be conjugate (Hairer et al., 2006, Section VI.8) to each other. To get a trapezoidal orbit, we need only apply a correction at the beginning and ending of a midpoint rule integration. This means the trapezoidal rule solution should have similar error properties to a symplectic method like the midpoint rule; we will show this more carefully in Section 3.3.

3.2 Runge–Kutta methods

The numerical algorithm (2) is a mapping of the vector space onto itself. A broad class of integrators, that encompasses various common algorithms including the ones of Section 3.1, defines the mapping using only and derivative operators that are scalars under coordinate transformations of . They are called Runge–Kutta (RK) methods. An RK method of stages is defined by constants and for when there is no explicit time-dependence in the governing ODE’s:


which is explicit, and thus less computationally expensive, if and only if for (a strictly triangular matrix). For this method, depends on and on differential operators like that are scalars under general coordinate transformations . The popular leapfrog method is not an RK method—it is known as a partitioned Runge–Kutta method— because it uses a different rule for updating the positions and momenta. In the partitioned RK case, the differential operators defining the are no longer covariant under general linear transformations of the full space.

The methods (17) and (18) are RK methods. We can check that for the former, , , , and . For the latter, , , and . Both are implicit and thus will need to be solved through iteration, whether fixed-point or Newton-Rhapson (Press et al., 2002).

It is easy to see both the implicit midpoint (as opposed to explicit midpoint, a different RK method) and trapezoidal rule are time-symmetric (and reversible, cf. Section 2.1 ) if used with fixed time step. Take a step forwards from to obtain and a step backwards to obtain . The rules require .

Next, we investigate whether the methods are symplectic. It has been shown (Hairer et al., 2006, Section VI.4), that if and only if an RK method conserves quadratic invariants in the phase space variables of the underlying differential equations, it is symplectic. The reason is related to the fact that the functions of , , defined by equations (60), which must be invariant for symplecticity to hold, are first integrals of the variational equations. One such typical quadratic invariant in some problems is the angular momentum. Any quadratic invariant can be written , with a symmetric matrix. Write the implicit midpoint rule as


Multiply from the left by . The left hand side gives


which follows from the fact that the transpose of a scalar is the scalar. The right hand is zero because . Thus, we are left with , which means the implicit midpoint rule conserves quadratic invariants and is thus symplectic. Any numerical experiment with a symplectic integrator that is an RK method will show conservation of all quadratic invariants; an example is the angular momentum, for differential equations that have this symmetry, such as the Kepler problem. We show a more direct proof of the symplecticity of the midpoint rule in Appendix A.

Write the trapezoidal rule as


If we multiply on the left by , we find that and is generally not conserved, meaning quadratic invariants are not conserved, and the trapezoidal rule is not symplectic. Numerical experiments indeed show the trapezoidal rule does not conserve quadratic invariants such as the angular momentum.

We will largely focus on the trapezoidal rule for the remainder of the paper, because we are interested in a time-symmetric, but non-symplectic integrator, and this method is one of the simplest examples of this. Some researchers have used leapfrog, which is symplectic when using fixed time step, with reversible steps. Once the steps are adapted, however, the symplectic property is lost, so there is no advantage from this standpoint to use leapfrog. On the other hand, even when used with adaptive steps, leapfrog conserves angular momentum exactly, while the trapezoidal rule does not. However, the tests we consider in what follows have no angular momentum invariant. Also, the trapezoidal rule has a related invariant for every quadratic, in the phase space, invariant in the underlying equations; see Section 3.3. Both leapfrog and trapezoidal rule conserve linear invariants, such as the total linear momentum, exactly (all RK methods do). A disadvantage of the trapezoidal rule is that it requires solving implicit equations, unlike leapfrog. But any time-symmetric Runge–Kutta method is implicit. We focus our efforts on Runge–Kutta methods because their properties have been well established, and they treat all phase space components with the same functional update rule, which will simplify our analysis of their symplecticity and energy conservation properties in Appendices C and D. An advantage of the trapezoidal rule, as shown in Section 3.1, is its connection to a symplectic method.

3.3 A conserved quantity for the trapezoidal rule

Consider the broad class of separable Hamiltonians,


where is the phase space dimension, and define


Additional derivatives of with respect to the are denoted by more subscript indices. In Section 3, let and refer to from the trapezoidal and implicit midpoint rule, respectively. Let other functions be functions of . As an example,

In Appendix B, we derive the modified differential equations for the trapezoidal and implicit midpoint rule, and the Hamiltonian for the implicit midpoint rule. Using the results from Appendix B, we find


where is the midpoint Hamiltonian given by (72). Note we do not follow Einstein summation convention, but we could restore the convention, for example, by substituting for . Using this information, we can compute that along the trapezoidal trajectory,


This equation describes the energy drift of trapezoidal rule. The term can be integrated with respect to time, and we find that




the trapezoidal rule has a conserved energy at least to second order. We will check this numerically in Section 4.1. This means a time-symmetric, non-symplectic method can also have a conserved energy at some order, but this fact may not in of itself be useful. The trapezoidal and midpoint rule, and DKD and KDK leapfrog all have a conserved energy to second order, which, for Hamiltonian (25) has form


and their coefficients and are shown in Table 1. and differ from each other for the leapfrog methods because they are partitioned RK methods, as mentioned in Section 3.2.

KDK Leapfrog
DKD Leapfrog
Table 1: The midpoint and trapezoidal rules, and KDK and DKD leapfrogs, have a conserved energy to second order described by (31). They only differ in the values of the coefficients of and , whose absolute value is either or , and we list them here.

We can do better and show that trapezoidal rule has a conserved energy to at least fourth order. Substituting its MDE (B) into equations (C), reveals,


For a conventional Hamiltonian (25), this becomes

In fact, we are able to show that the trapezoidal rule conserves an energy function to all orders in , and we can write it down. Rewrite the trapezoidal rule as a sequence of three RK steps:


The first step is a backwards explicit Euler step, the second is a symplectic midpoint method, and the third is an implicit Euler step. Because the implicit midpoint rule has a conserved Hamiltonian (assuming convergence of the series), it is natural to assume that the trapezoidal rule respects an energy function with the same functional form, but with shifted initial conditions. Indeed, let


Then, we can check , which implies that the trapezoidal rule conserves the energy function . If the underlying equations have a quadratic invariant , we also see the trapezoidal rule has a related invariant,


To fourth order, (34) agrees with equation (32), but it is exact to all orders. We will derive in Appendix C that there exist time-symmetric methods which are not energy conserving to all orders. These results are summarized in Table LABEL:table:RK. This means we can find energy drift with a symmetric integrator with fixed time step— this result has already been discussed by Faou et al. (2004) and others.

3.3.1 An example: the simple harmonic oscillator

We derive the conserved energy of the trapezoidal rule for the simple harmonic oscillator (SHO). The Hamiltonian for the SHO is


For this Hamiltonian, the trapezoidal rule becomes explicit, since the coordinate derivatives are linear in coordinates. Also, in this case, the implicit midpoint rule gives an identical rule. The rules say


When solved for and , they say




(3.3.1) is also the exact trajectory after time for a Hamiltonian


so long as




for , so the numerical value of the modified Hamiltonian is smaller than . When , an satisfying (41) does not exist, so the governing equations are no longer Hamiltonian. Thus, the trapezoidal and implicit midpoint rules’ MDEs are governed by (40). This implies they exactly conserve the energy of the SHO, as one can verify numerically.

For a general Hamiltonian (e.g. the Hénon–Heiles problem), these simple exact results no longer hold. However, for a time-independent Hamiltonian, symplectic methods always have a conserved energy, and so do conjugate methods like the trapezoidal rule.

3.4 Modified differential equation with adaptive time steps

In previous sections and the Appendix, we discuss integrators with fixed step-sizes, but for fixed step-sizes, there already exist excellent symplectic integrators in astrophysics, starting with leapfrog. Time-symmetric integrators are popular in astrophysics due to their ability to accommodate adaptive stepping. An exactly time-symmetric integrator was proposed by Hut et al. (1995), and approximately time-symmetric integrators have been developed by Pelupessy et al. (2012) and Kokubo et al. (1998). We focus on the proposal by Hut et al. (1995), because it is exactly time-symmetric, under certain conditions we describe. In conjunction with leapfrog, they propose to write the time step as an implicit equation,


is a function that we can specify using a priori knowledge about the solution trajectory (e.g., the relevant timescales) or even without this knowledge (Stoffer, 1995). An implicit step criterion can be used with an implicit one-step method, like the trapezoidal rule, not necessarily resulting in more iterations when solving the update equations. (43) can be written as an explicit infinite series in ,


The direction of time is now determined by the sign of . Of course, the depend on the method. Letting , for trapezoidal rule,


where is the phase space dimension. is the th component of in eq. (1) and . Symmetry requires


For criterion (43), this requirement is automatically satisfied. -reversibility would require . For steps (43) this means (Hairer et al., 2006, Section VIII.3)


This condition is easy to satisfy, but it is not always satisfied. Hut et al. (1995) were interested in the -body problem and proposed a that is the minimum of the close encounter and free fall times. This satisfies (47) if we are taking the absolute values of relative velocities. We will explore what happens when eq. (47) is not obeyed in Section 4. The equations (43), (46), and (47) apply whether the underlying method is an RK method, like trapezoidal rule, or a partitioned Runge–Kutta method, like leapfrog. But the underlying method must be time-symmetric for either (46) or (47) to hold.

The stepping rule (43) is implicit, which is more cumbersome to analyze than an explicit rule. However, we use an implicit criterion for the following reasons:

  • The trapezoidal rule is already implicit, so choice (43) does not necessarily add more computational work.

  • (43) has already been used by Hut et al. (1995) and others.

  • Dehnen (2017) studies explicit stepping criteria with step sizes that can only take certain values. The discreteness of the step sizes breaks the time-symmetry and reversibility symmetries. We want to construct exactly symmetric and reversible methods for our tests to be able to conclude that errors are not due to breaks in these symmetries.

  • The explicit stepping criteria in (Dehnen, 2017, Sections 4 and 5), even in the continuous, non-discrete case, risk becoming unsynchronized with , leading to stepping of questionable efficiency and accuracy. This stepping can be regarded as a multistep method.

We can construct the MDE with adaptive time steps, following the procedure of Section 2.2 and using the form (44), so that the series are now written in terms of . Now, instead of eq. (27), we have


As in eq. (28), we can calculate the energy drift along the trapezoidal orbit as


If we let and , this expression is just (28). This shows that the energy drift is a function of the problem and the choice of . In general, this term cannot be integrated in terms of elementary functions. There exist reversible (cf. eq. 47) which lead to secular drift and irreversible which lead to no energy drift, as we will show in Section 4. (49) holds for any , not just (43), so long as to lowest order in , . For example, (49) applies to the geometric mean time step,


4 Numerical Demonstration

In this section we apply the error analysis of Section 3 to see that energy conservation is violated in a number of situations, even when a method is symmetric and reversible. For the tests, we consider the pendulum solution and Hénon–Heiles orbits.

4.1 The modified pendulum

In the following, we consider the trapezoidal rule (17) along with the adaptive step criteria (43). Consider the simple pendulum, with a modified potential:


This modified pendulum was considered by Faou et al. (2004). The reason for choosing this potential with will become apparent below. This Hamiltonian is -reversible: , and .

The potential is not symmetric in over the periodic range of , as seen in Fig. 1. The minimum is and occurs at .

Figure 1: Potential as a function of the periodic range of for the modified pendulum Hamiltonian (51). There is no symmetry in the potential and it has a minimum of .

First, we choose , so that the step is constant. Because this satisfies (47), the integrator is reversible. We choose : there are roughly 100 steps per period. We let . As initial conditions, we choose and so that . So for the exact solution (and in all our numerical solutions), the sign of the momentum does not change (the orbit is circulating). In Fig. 2, we show that the change in various phase space quantities in time mimics the behavior of a symplectic integrator.

Figure 2: The evolution of some phase space quantities as a function of time when we integrate the modified pendulum with the symmetric non-symplectic trapezoidal rule. We use a constant and initial conditions and . The top panel gives the change in energy as given by (49). The middle panel gives the energy error, and the bottom panel gives the error of the conserved second order energy. No energy drift is observed in the middle panel, in agreement with the top panel. The second order energy is conserved better than the energy, as expected.

is given by the second order term of eq. (49):


which oscillates symmetrically around 0. We also show the energy and error, from (31). and are the initial energy and values, respectively. is conserved better than , supporting the finding that a conserved energy exists. We checked that, for a fixed integration time, and .

We next choose an adaptive step strategy. We choose , so that . This choice is both reversible and time-symmetric, according to the discussion above (47). We checked if we integrate forwards, change the sign of , and integrate the same number of steps backwards, we recover the initial conditions up to roundoff error. We choose so that the average time step is approximately still the same as the previous test. The initial phase space coordinates remain the same in this test: and . We initialize the integration with guess for the initial step, , and thereafter we use the previous time step as the initial guess. We integrate for the same total time as Fig. 2, and show the results in Fig. 3.

Figure 3: Same as Figure 2, but now with a reversible, time-symmetric step size strategy, . . We observe a linear energy drift, which roughly agrees with the prediction from the top panel, given by the red curve in the lower panel. This integration is fully symmetric and reversible, yet shows a linear energy error drift.


now is not symmetric around 0. This leads to a linear drift in energy error as seen on the bottom subplot. If is the value of at time step number , and is the value of the time step, define


which we expect to be close to , the energy at step . We see in Fig. 3 this is the case at small time, but the approximation breaks down for larger time. As we decrease , the difference between the two curves becomes undetectable on the same type of plot. We checked the slope of the curve scales as , confirming the energy error is . All other we tested, reversible and irreversible, gave a linear drift in energy for this problem. For the geometric mean time step (50), the errors are similar, as expected. The final energy error at changes by less than .

We integrate (51) with one-step symplectic methods of different properties. We use stepsize and initial conditions and (as in Section 4.1). The methods are described in Table LABEL:tab:modpend, where we write the method’s name, order, Runge–Kutta classification, and whether the method is time-symmetric. None of the methods yield energy drift.

Method Symmetric? Order Classification
Leapfrog DKD and KDK Yes 2 Partitioned Runge–Kutta
Gauss-Legendre Yes 4 Runge–Kutta
Symplectic Euler No 1 Runge–Kutta
Table 2: Description of symplectic integrators used in tests of the modified pendulum Hamiltonian.

By contrast, it has been reported that Lobatto IIIA and IIIB, two fourth order symmetric and reversible, but non-symplectic methods show energy drift when used with fixed step size. The tests were performed in Faou et al. (2004), using and . Faou et al. (2004) differs from our work in that it considered only fixed time-steps. We confirmed the leading order error for Lobatto IIIB is while the leading order error for Lobatto IIIA is and has no error drift (at leading order), except for roundoff error contributions. Lobatto IIIA and IIIB have the same symmetries as the trapezoidal rule with adaptive symmetric and reversible steps: time-symmetry and time-reversibility.

The simplified Takahashi–Imada method is symmetric and volume preserving. This means it conserves one Poincaré invariant, which does not guarantee symplecticity. It was reported in Hairer et al. (2009) that this method gives secular drift in a function close to the energy for Hamiltonian (51). They used and . In our own tests, we found the method gives energy drift. Volume preservation is equivalent to symplecticity for one-degree-of-freedom problems such as Hamiltonian (51), so this would appear to be an example of a symplectic integrator giving energy drift. But this integrator is generally non-symplectic.

The only example we found of a reversible, symmetric, and non-symplectic method conserving energy for this problem is studied in Fig. 2. We have not tested symplectic integrators with adaptive steps because symplecticity is not conserved and this advantage of the method is lost.

When the orbits of these initial conditions are computed with symplectic integrators, we have observed that different steps along the orbit produce increases or decreases in energy. The net increase is zero. For the orbit of Fig. 3, the net increase is negative over an orbit, leading to energy drift.

For the unmodified pendulum,


consider the same initial conditions given above. The sign of is still invariant in all tests. Both reversible and irreversible , such as and with and constants, give no drift in energy error. However, we again get drift if we let be an asymmetric function of , such as the we used above, even though this is reversible. These experiments show that time reversibility and energy conservation are independent concepts. To summarize, for circulating pendulum orbits, all time-symmetric methods, except the fixed time trapezoidal rule, gave undesirable error behavior. This exception may be related to the fact it is a conjugate symplectic method (see Section 3.1). All tested symplectic methods except the simplified Takahashi-Imada method yielded desirable energy conservation. The simplified Takahashi-Imada method is generally non-symplectic.

4.2 Hénon–Heiles orbits

The Hénon–Heiles problem Hamiltonian (Henon & Heiles, 1964) is a two-degree-of-freedom problem—a simplified model of a galactic potential. The Hamiltonian is


with . This Hamiltonian allows both chaotic and regular trajectories and is -reversible. We consider a regular orbit with initial conditions , , , , and (although we only need to keep three significant figures to get the same qualitative results). Using , and , we show the trajectory in Fig. 4 with . Also plotted is the bounding equipotential curve, .

Figure 4: A regular box Hénon–Heiles orbit. The initial conditions are , , , and . We also plot the bounding equipotential. The integration is run until with a constant step-size trapezoidal method.

The trajectory does not span the entire allowed area, which tells us it is a regular orbit. We can verify this in a surface of section plot. In Fig. 5, we plot a point in the plane every time is crossed with , up to time . We use a fifth and sixth order pair of Runge–Kutta methods, in what’s known as Verner’s embedded Runge–Kutta method (Verner, 1978), for this plot. This is an adaptive step method: two methods allow an estimate of the local truncation error which is then used to determine a step size. The final energy error is . This orbit is a box orbit: the sign (and magnitude) of the angular momentum oscillates. If we change the sign of the initial momenta, the trajectory is confined to the same bounding curve, which means the second isolating integral besides the energy (for these initial conditions) does not depend on the sign of the momenta. This is a consequence of the -reversibility of the equations due to (55). We checked this by running the trajectory with a sign change in the initial momenta and checking that the minimum and maximum of the trajectory is the same to 15 significant figures.

Figure 5: Surface of section plot for the orbit of Fig. 4. A point is plotted everytime is crossed with . The symmetry in indicates a box orbit, and the closed curve indicates a regular orbit far from resonance. The surface of section is computed up to time with a high accuracy Runge–Kutta method.


For this Hamiltonian, for eq. (49), we have