PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Esaim Math Model Numer Anal. Author manuscript; available in PMC 2010 July 1.
Published in final edited form as:
Esaim Math Model Numer Anal. 2009 July 1; 43(4): 631–644.
doi:  10.1051/m2an/2009019
PMCID: PMC2881582
NIHMSID: NIHMS136293

ON ENERGY CONSERVATION OF THE SIMPLIFIED TAKAHASHI-IMADA METHOD*, **

Abstract

In long-time numerical integration of Hamiltonian systems, and especially in molecular dynamics simulation, it is important that the energy is well conserved. For symplectic integrators applied with sufficiently small step size, this is guaranteed by the existence of a modified Hamiltonian that is exactly conserved up to exponentially small terms. This article is concerned with the simplified Takahashi-Imada method, which is a modification of the Störmer-Verlet method that is as easy to implement but has improved accuracy. This integrator is symmetric and volume-preserving, but no longer symplectic. We study its long-time energy conservation and give theoretical arguments, supported by numerical experiments, which show the possibility of a drift in the energy (linear or like a random walk). With respect to energy conservation, this article provides empirical and theoretical data concerning the importance of using a symplectic integrator.

Keywords and phrases: Symmetric and symplectic integrators, geometric numerical integration, modified differential equation, energy conservation, Hénon-Heiles problem, N -body problem in molecular dynamics

Introduction

For accurate simulations in molecular dynamics it is generally accepted that it is essential to have an integrator that nearly conserves energy on long time intervals and nearly preserves phase space volume (“what makes molecular dynamics work”, see [12,15]). Being symplectic implies volume-preservation and ensures energy conservation for suficiently small time steps, but is it necessary to use a symplectic integrator? Another interesting question concerns the existence of a modified conserved energy for integrators that are merely volume preserving. The aim of this article is to present an example of a potentially important time integrator suggesting that the answer is yes for the first question. It also gives evidence that volume-preserving dynamics does not always possess a modified energy that is conserved for exponentially long time.

Since the beginning of simulations in molecular dynamics, the Störmer-Verlet integrator has been a method of choice and much in use. It is explicit for systems with separable Hamiltonian (and thus easy to implement); it is symmetric, time-reversible, symplectic, and volume-preserving. All essential geometric properties of the exact flow are reflected in the numerical discretization. Its only disadvantage is its low order two, which results in lack of accuracy. (This is not always a concern though, because of much greater inaccuracy from statistical error.)

In Section 1 we introduce the Takahashi-Imada integrator and its simplified version. This simplified Takahashi-Imada method is symmetric and volume-preserving, but not symplectic. We therefore review in Section 2 what is known about the near energy conservation of symmetric time-reversible integration methods. To study the long-time behavior we make use of backward error analysis and compute in Section 3 the modified differential equation of the simplified Takahashi-Imada method. Since the method is not symplectic, the modified differential equation is not Hamiltonian, and we show that it does not have a modified conserved energy. In Section 4 we discuss different behaviors that may arise: (i) the error of the energy can remain bounded and small; (ii) there is the possibility of a linear drift and; (iii) most important for molecular dynamics simulations, for chaotic systems it typically behaves like a random walk, so that the typical error increases like square root of time. Such a random walk behavior has also been observed by [8] for a non-holonomic mechanical system which is reversible and energy-preserving but not Hamiltonian, and by [4] for a Hamiltonian system with non-smooth potential. Section 5 illustrates the results of this article for an unsymmetric pendulum equation, for a modification of the Hénon-Heiles problem, and for an N -body problem with bodies interacting via the Lennard-Jones potential.

1. Simplified Takahashi-Imada method

Consider a Hamiltonian system of the form q = −[nabla]U (q) with smooth potential U (q). The total energy (kinetic plus potential)

H(p,q)=12pp+U(q)
(1.1)

is exactly conserved along its flow. Written as a first order differential equation, it becomes

p.=f(q)q.=pwithf(q)=U(q).
(1.2)

All statements of this article can be extended straightforwardly to problems Mq = −[nabla]U (q) where M is a constant mass matrix. For notational convenience we assume M = I, I being the identity matrix.

1.1. Störmer-Verlet method

For a differential equation q = f (q), the replacement of the second derivative by a finite difference yields the most natural discretization qn+1 − 2qn + qn−1 = h2f (qn). With pn = (qn+1qn−1)/2h as approximation to the derivative p = q, this integrator can be written in the form of a one-step method as

pn+1/2=pn+h2f(qn)qn+1=qn+hpn+1/2pn+1=pn+1/2+h2f(qn+1).
(1.3)

It is usually called the Störmer-Verlet, velocity-Verlet or leapfrog method. It has all nice properties of a geometric integrator. It is symplectic and volume-preserving, it is symmetric and reversible, but it has only order two. This means that the local discretization error (contribution of one step) is of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h3), so that the global error in the solution as well as in the energy is An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2).

Another way to formulate the Störmer-Verlet method as a one-step method is to introduce qn+1/2 = (qn + qn+1)/2 and to consider the mapping (pn−1/2, qn−1/2) [mapsto] (pn+1/2, qn+1/2). This integrator (also called position-Verlet) has the same properties and all statements of the present article carry over in a straightforward way to this situation.

1.2. Takahashi-Imada integrator

The accuracy of the Störmer-Verlet discretization is greatly improved by adding the expression αh2f′(q)f(q) with α=112 to every force evaluation in (1.3). This yields the method (cf. [10,14])

pn+1/2=pn+h2(I+αh2f(qn))f(qn)qn+1=qn+hpn+1/2pn+1=pn+1/2+h2(I+αh2f(qn+1))f(qn+1).
(1.4)

It is symmetric and reversible, and it is symplectic and volume-preserving, because for f(q) = −[nabla]U(q) it can be interpreted as applying the Störmer-Verlet method to the Hamiltonian system with potential

V(q)=U(q)α2h2||U(q)||2.

Although it is still of classical order two, it is an effectively fourth order integrator for α=112. Denoting the method (1.4) by Φh: (pn, qn) [mapsto] (pn+1, qn+1), this means that there exists a change of coordinates χh such that χhΦhχh1 approximates the exact flow to order four. In other words, the integrator (1.4) is conjugate to a method of order four. Indeed, with the mapping χh: (p, q) [mapsto] ([p with hat], q) defined by

p=p^+ch2f(q)p^,q^=q+ch2f(q),
(1.5)

a straightforward computation by Taylor expansion shows that χhΦhχh1 is a method of order three if and only if c=α=112. Since this composed method is symmetric (which follows from χh = χh), it is automatically of order four. For f(q) = −[nabla]U(q) the transformation (1.5) is symplectic.

Notice that in a constant step size implementation (χhΦhχh1)n=χhΦhnχh1 holds, so that it is sufficient to apply the change of coordinates in the beginning and the end of the integration. The concept of effective order has been introduced in [2]. It is systematically used in the design of splitting methods for Hamiltonian equations (see for example [1]).

1.3. Simplified Takahashi-Imada integrator

To avoid the derivative evaluation of the vector field f (q), which corresponds to a Hessian evaluation in the case of f (q) = −[nabla]U (q), we replace (I + αh2f ′(q))f (q) with f (q + αh2f (q)) and thus consider (cf. [16])

pn+1/2=pn+h2f(qn+αh2f(qn))qn+1=qn+hpn+1/2pn+1=pn+1/2+h2f(qn+1+αh2f(qn+1))
(1.6)

which is a An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h5) perturbation of the method (1.4). It is still symmetric, reversible, and of effective order four, and it is volume preserving and thus symplectic for problems with one degree of freedom [16]. To see this, note that (1.6) is a composition of shears (mappings of the form (p, q) [mapsto] (p + a(q), q) and (p, q) [mapsto] (p, q + hp)), and shears always preserve the volume.

On the other hand, a shear (p, q) [mapsto] (p + a(q), q) is symplectic only if a′(q) is symmetric. Thus the first step of (1.6) is symplectic only if f′(qn + αh2f (qn) f′(qn) is symmetric. This is the product of two (in general, different) symmetric matrices and is not symmetric for all f (q) = −[nabla]U (q). Consequently, the method is not symplectic for more than one degree of freedom. It is even not conjugate to a symplectic integrator. Otherwise, as a method that can be expanded into a P-series, it would have a conserved quantity that can be written as a series in terms of (globally defined) elementary Hamiltonians (see [7], Sect. IX.10). That this is not the case is evident from the numerical experiments of Section 5.

Despite the lack of symplecticity, the simplified Takahashi-Imada method is an interesting integrator because it is markedly cheaper than the known symplectic integrators of (effective) order 4, requiring only 2 force evaluations and no evaluations of the derivative of the force per time step. The aim of this article is to study the long-time energy conservation of this integrator.

1.4. Linear stability analysis

To get some insight into the maximal time step of practical value, we consider the harmonic oscillator with frequency ω, which is Hamiltonian with quadratic potential U(q)=12ω2q2. For the Störmer-Verlet method we get the recursion qn+1 − (2 − (ωh)2)qn + qn−1 = 0, which is stable only for ωh ≤2. It has a discrete energy

pn2+(1(ωh)24)ω2qn2=Const.

that is exactly preserved. The maximal relative error in the Hamiltonian equals (ωh)2/(4 − (ωh)2) and it is attained when pn = 0. This error behaves like An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) for small step size, but tends to infinity when ωh → 2.

Applied to the harmonic oscillator, both Takahashi-Imada integrators are identical and yield the recursion

qn+1(2(ωh)2+α(ωh)4)qn+qn1=0.

For α=112 this recursion is stable if and only if ωh<23. It also has a discrete energy which is

pn2+(1β(ωh)24)βω2qn2=Const.withβ=1(ωh)212.

As for the Störmer-Verlet method the maximal relative error in the Hamiltonian is of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) for ωh → 0, and tends to infinity when ωh23. However, when we consider the transformed variables of (1.5), [p with hat]n = β1pn and qn = βqn, this discrete energy becomes

p^n2+(1β(ωh)24)β3ω2q^n2=Const.

The maximal relative error in the original Hamiltonian is (4−β(ωh)2−4β3)/(4−β(ωh)2). Due to an unexpected cancellation of the fourth order terms, it behaves like An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h6) for small step sizes and remains bounded (converges to one) when ωh23.

When comparing the Takahashi-Imada integrator with the Störmer-Verlet method, one has to take into account that the simplified Takahashi-Imada method is twice as expensive. Therefore, the effective step size of the Störmer-Verlet method is twice as large. The right picture of Figure 1 shows the maximal relative error in the Hamiltonian as a function of ωheff, where the effective step size heff is h for the Takahashi-Imada integrator and 2h for the Störmer-Verlet method. For all step sizes that give rise to bounded numerical solutions, the Takahashi-Imada integrator is more accurate for the same amount of work. A similar picture is obtained when the global error is plotted as function of the effective step size. If εp and εq denote the errors in the p and q component, respectively, we consider the scaled norm (p|2 + 2q|2)1/2, so that the error only depends on the product ωh. The global error, which mainly consists of a phase error, increases linearly with time. In the left picture of Figure 1 we show the maximal global error on the interval 50/ω. Notice that for the Takahashi-Imada method this error behaves like An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h4). The higher order An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h6) is observed only in the Hamiltonian.

Figure 1
Harmonic oscillator: error in the solution (left picture) and error in the Hamiltonian (right picture) as a function of the effective step size for the Störmer-Verlet method (dashed line) as well as for the solution ([p with hat]n, q ...

2. Review on energy conservation with symmetric methods

One of the main motivations of using symplectic time integrators is their good energy conservation over long times. The simplified Takahashi-Imada integrator is no longer symplectic. What can be said about its long-time energy conservation? This section reviews energy conservation by symmetric (non-symplectic) methods.

2.1. Completely integrable differential equations

For completely integrable Hamiltonian systems there exists a perturbation theory which permits one to prove that, for symplectic methods applied with small step size, the numerical solution nearly conserves all action variables over long times and in particular also the Hamiltonian. There is an analogue for reversible systems solved with symmetric integrators.

Consider a differential equation (not necessarily Hamiltonian) with variables split into (u, v) such that

u.=f(u,v)v.=g(u,v)wheref(u,v)=f(u,v)g(u,v)=g(u,v).
(2.1)

The problem (1.2) is of this form with u = q and v = p, but often there are other splittings that yield this form, too. The exact flow of this differential equation is reversible, which means that with the transformation ρ: (u, v) [mapsto] (u,v) it satisfies ρϕtρ=ϕt1. Many symmetric methods, and in particular also the simplified Takahashi-Imada method, are reversible and thus satisfy ρΦhρ=Φh1.

The system (2.1) is called an integrable reversible system, if there exist a change of coordinates a = μ(u, v), θ = ν(u, v) preserving reversibility (i.e., μ(u,v) = μ(u, v) and ν(u,v) = −ν(u, v)) and a frequency vector ω(a) such that in the new variables the system is of the form

a.=0,θ.=ω(a).

The action variables a are constants of motion, and the angle variables θ increase linearly with time. Based on a reversible KAM theory it is proved in [7], Chapter XI, that, under a non-resonance condition on the frequencies ωj (a) and for sufficiently small step sizes, reversible numerical methods applied to integrable reversible systems have the same long-time behavior as symplectic methods applied to integrable Hamiltonian systems. In particular, they nearly preserve all action variables (including the total energy in the case of a Hamiltonian system) on exponentially long time intervals and the global error increases only linearly with time. More precisely, as long as the solution remains in a compact set, the error in the action variables is bounded by An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hr) + An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(teγ/h), where r is the order of the method.

2.2. General (chaotic) Hamiltonian systems

Problems in molecular dynamics are in general chaotic, which means that they are very sensitive to perturbations in initial data. Assuming that the numerical solution remains in a compact set and that the step size is sufficiently small, the local error contributions to the Hamiltonian just sum up. A worst case analysis therefore shows that the energy error at t = nh is bounded by An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(thr) when r is the order of the method.

2.2.1. Conjugate-symplectic methods

Since the symplecticity requirement is rather strong, Stoffer [13] suggests considering methods that are conjugate to a symplectic method. Denoting the discrete flow of the method by Φh, this means that there exists a change of coordinates χh (close to the identity) such that Ψh:=χhΦhχh1 is a symplectic transformation. The numerical solution obtained with Φh then has the same long-time behavior as that obtained with the symplectic integrator Ψh, and the energy error will be of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hr) + An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(teγ/h), so that the lifetime of energy conservation is exponential. A well-known example is the symmetric (non-symplectic) trapezoidal rule which is conjugate to the symplectic midpoint rule.

The underlying one-step method of symmetric linear multistep methods is another example that is conjugate to a symplectic method (see [3,6] and [7], Sect. XV.4). The lifetime of energy conservation would be exponential if there were not disturbing parasitic solutions of the difference equation. Only for symmetric methods for second order differential equations q = −[nabla]U (q), where the ρ-polynomial does not have double roots other than 1, these parasitic solutions are under control. It is proved in [6] (see also [7], Sect. XV.6) that for an rth order method the energy error is An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hr) on time intervals of length t = nh = An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hr−2) or even longer.

A special case of such a method for second order equations q = f (q) = −[nabla]U (q) is the implicit 4th order method (known as Cowell’s method or Numerov method)

qn+12qn+qn1=h212(f(qn+1)+10f(qn)+f(qn1)).

With suitable velocity approximations pn it is conjugate to a symplectic integrator and thus has an excellent long-time behavior (see [11]).

2.2.2. Methods that are conjugate to a symplectic method up a certain order

Relaxing the property of conjugate-symplecticity, we can ask for methods that are conjugate to a symplectic method up to order s, i.e., there exists a change of coordinates and a symplectic integrator Ψh such that χhΦhχh1Ψh=O(hs+1). If the method is of order r and the trajectory remains in a compact set, the numerical Hamiltonian will behave like

H(pn,qn)=Const.+O(hr)+O(ths),
(2.2)

so that the energy error is bounded by An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hr) on intervals of length t = nh = An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hrs), and an energy drift can appear only beyond such an interval (we always have sr). For integrable systems not only the energy behaves like (2.2), but all first integrals that can be expressed in terms of action variables. For general systems, quadratic first integrals like the angular momentum in N -body dynamics are also well conserved and behave like (2.2).

The question of conjugate-symplecticity up to a certain order is addressed in the thesis of Pierre Leone, and the most important results are collected in [7], Section VI.8. Among B-series methods (which include Runge-Kutta methods) all symmetric second order methods are conjugate-symplectic up to order 4. The fourth order Lobatto IIIB method is not conjugate-symplectic up to an order higher than four, but the fourth order Lobatto IIIA method is conjugate-symplectic up to order six.

2.2.3. Methods that preserve a modified energy up a certain order

We can still further relax the requirement and consider methods that conserve a modified energy up to order s. For B-series methods this property is discussed in [5]. In this situation we have (2.2) for the numerical Hamiltonian, and the energy is well conserved on intervals of length An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(hrs). However, one has to be aware of the fact that further properties of symplectic integrators are lost for computations on intervals longer than An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(1). For example, there is no guarantee that the angular momentum is well conserved.

3. Modified differential equation

To get insight into the long-time behavior of numerical integrators (including their energy conservation), backward error analysis is the most powerful tool. For a detailed and complete presentation of this theory we refer to [7], Chapter IX. The idea is to “interpolate” the numerical solution by the exact flow of a modified differential equation which is given by a formal series in powers of the step size

(p.q.)=(f(q)p)+h(f2(p,q)g2(p,q))+h2(f3(p,q)g3(p,q))+
(3.1)

Denoting by [p with tilde](t), q(t) the solution of this system, we formally4 have for the numerical solution that pn = [p with tilde](tn), qn = q(tn) with tn = t0 + nh. This permits to get information for the numerical solution by studying the exact flow of the modified differential equation (3.1). The reason why backward error analysis is so successful is that geometric properties of the numerical integrator are reflected in the modified differential equation. For example,

  • for symmetric methods, the modified differential equation is in even powers of h;
  • for reversible methods applied to a reversible differential equation, the modified equation is reversible;
  • for volume-preserving methods applied to a divergence-free system, the modified differential equation is divergence-free;
  • for symplectic methods applied to a Hamiltonian system, the modified equation is (locally) Hamiltonian.

In particular, for a symplectic method applied to (1.2) there exists a modified Hamiltonian

H(p,q)=12pp+U(q)+hH2(p,q)+h2H3(p,q)+h3H4(p,q)+

such that the modified differential equation is given by

p.=qH(p,q),q.=pH(p,q).

Notice that the functions fj (p, q) and gj (p, q) of the modified differential equation depend on p and q although the original problem (1.2) is a second order differential equation. Similarly, the modified Hamiltonian is in general not separable independent of whether the Hamiltonian of the original problem is separable or not.

3.1. Störmer-Verlet method

The Störmer-Verlet method is symplectic, so that its modified differential equation is Hamiltonian. Using the notation T(p)=12pp for the kinetic energy, and {A, B} =[nabla]qAT[nabla]pB[nabla]pAT[nabla]qB for the Poisson bracket of two functions depending on p and q, its modified Hamiltonian is

H=T+U+h2(112{T,{T,U}}124{U,{U,T}})+h4(1720{T,{T,{T,{T,U}}}}+1360{U,{T,{T,{T,U}}}}1480{U,{U,{T,{T,U}}}}+1120{T,{T,{U,{U,T}}}})+O(h6).
(3.2)

This formula is obtained from the symmetric Baker-Campbell-Hausdor3 formula (cf. [7], Sect. III.4.2).

3.2. Takahashi-Imada method

Substituting U in (3.2) with Uα2h2||U||2=Uα2h2{U,{U,T}} (cf. Sect. 1.2) yields the modified Hamiltonian of the symplectic Takahashi-Imada integrator (1.4):

H^=T+U+h2(112{T,{T,U}}(124+α2){U,{U,T}})+h4(1720{T,{T,{T,{T,U}}}}+1360{U,{T,{T,{T,U}}}}1480{U,{U,{T,{T,U}}}}+(1120α24){T,{T,{U,{U,T}}}}α24{U,{T,{U,{U,T}}}})+O(h6).
(3.3)

Remark 3.1

Let (pn, qn) denote the numerical solution of (1.4). With the values p^n=pn+112h22U(qn)pn and q^n=qn(112+α2)h2U(qn) we have H ([p with hat]n, qn) = Ĥ(pn, qn) + An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h4) which implies an improved energy conservation H([p with hat]n, qn) = Const. + An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h4) for any value of α. Notice, however, that the numerical solution (pn, qn) is of effective order 4 only if α=112 (cf. Sect. 1.2).

The choice α=112 gives qn = qn, so that the positions need not be modified. Approximating [nabla]2U (qn)pn in the formula for [p with hat]n by a difference of force evaluations, this choice yields the advanced SV algorithm of [9].

3.3. Simplified Takahashi-Imada method

To find the modified differential equation of (1.6) we notice that it is a An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h5) perturbation of the Takahashi-Imada method. If we denote by [p with hat]n+1/2, qn+1, [p with hat]n+1 the values obtained from (1.4) subject to initial values qn, pn, then the numerical approximation obtained from (1.6) satisfies with qn+1/2 = (qn + qn+1)/2

pn+1=p^n+1+12h5α2(f(f,f))(qn+1/2)+O(h7)qn+1=q^n+1+14h6α2(f(f,f))(qn+1/2)+O(h8).

The modified differential equation of (1.6) is therefore given by

p.=qH^(p,q)+12h4α2(f(f,f))(q)+O(h6)q.=pH^(p,q)+O(h6)
(3.4)

where Ĥ(p, q) is the modified Hamiltonian (3.3) of the Takahashi-Imada method. Due to the symmetry of the modified Takahashi-Imada method (1.6) the expansion (3.4) is in even powers of h.

4. Energy conservation of the simplified Takahashi-Imada integrator

To study energy conservation of the simplified Takahashi-Imada integrator (1.6) we notice that along the (formal) exact solution of its modified differential equation (3.4) we have

ddtH^(p,q)=12h4α2pf(f,f)+O(h6)=12h4α2Uqqq(p,Uq,Uq)+O(h6),
(4.1)

where we have used f(q) = −[nabla]U (q) = −Uq (q). In general, the expression on the right-hand side of (4.1) cannot be written as a total differential. However, there are several possibilities to extract parts that can be written as a total differential which then gives further information on the long-time energy conservation of the method.

To better explain the numerical results in Section 5.2.1, we suggest to use the identity

ddt(2Uqq(Uq,Uq)4Uqqq(p,p,Uq)Uqqqq(p,p,p,p)+2(Uqqp)(Uqqp))=10Uqqq(p,Uq,Uq)Uqqqqq(p,p,p,p,p)

which holds along solutions of the Hamiltonian system (1.2). A An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) error term has to be added, if the formula is used along the solution of the modified differential equation (3.4). Relations of this type can be derived with the use of elementary Hamiltonians (cf. [7], Sect. IX.10.2). Extracting the term Uqqq (p, Uq, Uq) from this formula, we obtain

Uqqq(p,Uq,Uq)=ddtK(p,q)+110Uqqqqq(p,p,p,p,p)+O(h2),
(4.2)

where

K(p,q)=15Uqq(Uq,Uq)25Uqqq(p,p,Uq)110Uqqqq(p,p,p,p)+15(Uqqp)(Uqqp).

It therefore follows from (4.1) that

ddt(H^(p,q)+12h4α2K(p,q))=120h4α2Uqqqqq(p,p,p,p,p)+O(h6)
(4.3)

along solutions of the modified differential equation (3.4). Since the numerical solution of the simplified Takahashi-Imada method is (formally) identical to the exact solution (at tn = nh) of (3.4), this formula yields much insight into the long-time energy conservation of (1.6) and can explain the different behavior depending on the structure of the problem.

  • Bounded error in the energy without any drift. If the potential U (q) is a polynomial of degree at most 4 (e.g., the Hénon-Heiles potential), the h4-term of the right-hand side of (4.3) vanishes, so that the Hamiltonian H(p,q)=12pp+U(q) satisfies H(pn, qn) = H0 + An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) + An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(th6) with t = nh along the numerical solution of (1.6), provided it remains in a compact set. A worst case analysis shows that the energy error is of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) with lifetime tn = nhAn external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h−4). As in [5] we expect that also higher order terms can be treated as above, so that no energy drift at all can be observed in this situation.
  • Linear energy drift. In general, without any particular assumption on the potential, the numerical energy will have a drift of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(th4), so that the lifetime of energy conservation (with error of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2)) of the simplified Takahashi-Imada integrator is at least tn = nhAn external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h−2).
  • Random walk behavior – drift like square root of time. Under the assumption that the solution of the modified differential equation is ergodic on an invariant set A with respect to an invariant measure μ, we have
    limt1t0tF(p(s),q(s))ds=AF(x)μ(dx),
    (4.4)
    where x = (p, q) and the function F (x) is the right-hand side of (4.3). Again there will be a linear drift of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(th4) in general. However, in the presence of symmetries, the integral of the right-hand side in (4.4) is likely to vanish and the numerical Hamiltonian will look like a random walk (a behavior observed by [4,8] in different situations). We therefore can expect that in addition to the bounded error terms in Ĥ(p, q) and K(p, q), there will be a drift of size O(th4). This implies that the lifetime with an An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) error of the simplified Takahashi-Imada integrator is at least tn = nhAn external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h−4). Notice that it is very difficult (or impossible) to verify the ergodicity assumption of the modified differential equation, but numerical experiments show that this behavior is a typical situation.

5. Numerical illustration

The aim of this section is to illustrate that the different behaviors concerning energy conservation, predicted by the backward error analysis of the previous sections, actually arise in practical computations. In all experiments the step size is sufficiently small so that the truncated modified equation is a reliable model for the numerical solution. For the case where the product of step size and the highest frequency in the system is not small (say, of size one or even closer to the linear stability bound of Sect. 1.4) a different analysis is required. We do not consider this situation in the present article.

5.1. Unsymmetric pendulum

Let us start with a generic example that does not have any symmetries in the solution. We consider the second order differential equation (1.2) with one degree of freedom and potential

U(q)=cosq+0.2sin(2q).
(5.1)

To the potential for the pendulum a perturbation is added so that the symmetry in q ↔ −q is destroyed. As initial values we take q(0) = 0 and p(0) = 2.5. The velocity is sufficiently large so that the pendulum turns round and p(t) remains positive for all time. Consequently, the symmetry p ↔ −p cannot have any influence onto a numerical solution. It turns out that the integral of the expression (4.3) over one period of the solution is non-zero. For the simplified Takahashi-Imada method we therefore expect that in addition to a bounded error term of size An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) there is a linear drift of the form An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(th4). Figure 2 shows the truncated modified energy 12pp+U(q)+h2H3(p,q) (shifted in the sense that its value at (p(0), q(0)) is subtracted) obtained with step size h = 0.2. This experiment confirms the theoretical investigations and shows that without any symmetry in the solution one cannot expect a better long-time behavior.

Figure 2
Unsymmetric pendulum: error in the shifted modified energy H(pn, qn)+h2 H3(p, q) along the numerical solution of the simplified Takahashi-Imada method (1.6).

There seems to be a contradiction to the fact that the simplified Takahashi-Imada method is symplectic for problems with one degree of freedom. To explain this apparent contradiction we recall that symplecticity alone is not sufficient for a good energy conservation. The numerical solution has to remain in a compact set and the modified Hamiltonian (in the sense of backward error analysis) has to be globally defined.

In this particular situation the phase space is either a cylinder, because the angle q is defined modulo 2π, or the Euclidean plane, if we do not take care of the periodicity of the vector field. On a cylinder, the coefficient function of h4 in the modified Hamiltonian is not globally defined, because the integral of the function (f″(f, f))(q) (see (3.4)) over an interval of length 2π is not zero. If we consider the Euclidean plane as the phase space of this problem, then the angle q grows linearly with time and the solution does not remain is a compact set.

5.2. Hénon-Heiles-type problems

Our next example is a Hamiltonian system (1.2) with two degrees of freedom and potential function

U(q1,q2)=12(q12+q22)+q12q21kq2k

which, for k = 3, reduces to the well-known Hénon-Heiles potential. If the value of the Hamiltonian is large enough, the solution is chaotic.

5.2.1. Cubic potential

We consider the case k = 3 and recall that in this case the dominant error term in formula (4.3) vanishes. We apply the simplified Takahashi-Imada method (1.6) with step size h = 0.2 and initial values q1(0) = 0, q2(0) = 0.2, p2(0) = 0.3, and a positive p1(0) is chosen such that the Hamiltonian is H0 = 0.125. Figure 3 shows the error in the modified energy H^(pn,qn)+12h4α2K(pn,qn) (with α = 1/12, cf. Eq. (4.3)). It behaves like An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h6) and no drift is visible, even not in a higher power of h. For this problem the simplified Takahashi-Imada method behaves exactly like a symplectic integrator. Since the problem is not integrable, a theoretical justification of this behavior would require a thorough investigation of higher order terms in the modified differential equation.

Figure 3
Hénon-Heiles problem with cubic potential: error in the modified energy H^(pn,qn)+12h4α2K(pn,qn) along the numerical solution of the simplified Takahashi-Imada method (1.6).

5.2.2. Quintic potential

To ensure that the dominating term in (4.3) does not vanish identically, we now put k = 5. The modified energy is no longer as well conserved as for k = 3 and the long-time error looks rather erratic. We therefore consider 100 different initial values randomly chosen close to the values q1(0) = 0, q2(0) = 0.2, p2(0) = 0.3, and p1(0) such that the value of the Hamiltonian is H0 = 1/8 (implying a chaotic solution). Figure 4 shows the modified energy H^(pn,qn)+12h4α2K(pn,qn) subtracted by its value at t = 0 along the numerical solution of the simplified Takahashi-Imada method applied with step size h = 0.2. Most of the trajectories have an error that behaves like a random walk throughout the considered interval. The lower picture of Figure 4 shows in addition to the errors also the average (μ = −0.11 × 10−6 at t = 2.5 × 106) and the standard deviation (σ = 0.54 × 10−6 at t = 2.5 × 106) as a function of time. One nicely observes that the average error is close to zero and that the standard deviation increases like the square root of time.

Figure 4
Hénon-Heiles problem with quintic potential: error in the modified energy H^(pn,qn)+12h4α2K(pn,qn) along the numerical solution of the simplified Takahashi-Imada method (1.6). The upper picture corresponds to 2 trajectories for which ...

However, for a few initial values (those plotted in the upper picture of Fig. 4), the error in the modified energy increases linearly with time on relatively large subintervals. We can see that for a particular initial value this happens on the interval [1 225 782, 1 415 265] which corresponds to 947 415 integration steps. To better understand this behavior, we study the Poincaré section of this particular solution. Whenever the solution (considered as a path in the phase space) crosses the surface q1 = 0 in the positive direction p1 > 0, we consider the point (q2, p2). These points are plotted in Figure 5. On the subinterval where the linear growth of the modified Hamiltonian is observed, these Poincaré cuts happen to lie in the region p2 < 0 (right picture), and one cannot expect that the integral over the expression in (4.3) vanishes. In general (left picture of Fig. 5) these cut points are symmetrically distributed in the phase space leading to a cancellation in the integral.

Figure 5
Poincaré cut for the numerical solution of the simplified Takahashi-Imada method corresponding to the interval [1.0 × 106, 1.1 × 106] with 12 827 cut points (left picture) and to the interval [1 225 782, 1 415 265] with 25 220 ...

This example illustrates that it is very difficult (even impossible) to verify rigorously any ergodicity assumption for the modified differential equation. Nevertheless, such an assumption permits one to explain the numerically observed long-time behavior very well.

5.3. A model problem for molecular dynamics

The N -body problem is a Hamiltonian system (1.2) with potential function

U(q1,,qN)=1j<iNV(rij),rij=||qiqj||.

In our experiment we assume that the positions qi of the particles are in R2, and that the interaction between them is via the Lennard-Jones potential V (r) = 0.4(r−12 − 2r−6). The initial position of the N = n2 particles (we use n = 3 for our experiment) is on the rectangular grid (i, j) for i, j [set membership] {1,, n}, and we assume zero initial momenta so that there is no evaporation. No escape of any particles is observed.

To study the energy conservation for this N -body problem, we apply the simplified Takahashi-Imada method (1.6) with step size h = 0.2 and many different initial values. The positions of the particles are chosen close to the rectangular grid with a random perturbation of maximal value 0.01 in each component. The initial momenta are always exactly zero. The result is shown in Figure 6. We have plotted the shifted modified energy H^(pn,qn)+12h4α2K(pn,qn) along the numerical solution for 50 trajectories. As bold curves we have included the average of the energy error over 400 trajectories together with its standard deviation.

Figure 6
Error in the modified energy H^(pn,qn)+12h4α2K(pn,qn) of the N -body problem along the numerical solution of the simplified Takahashi-Imada method (50 trajectories). Included is the average over 400 trajectories as a function of time (μ ...

In contrast to the previous experiment (modified Hénon-Heiles problem) we do not observe the existence of time intervals with linear growth in the energy error (there is much more mixing for problems with a larger number of degrees of freedom). These errors perfectly behave like a random walk, so that the error in the modified energy typically grows like An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(α2h4t1/2).

We emphasize that in Figure 6 we have plotted the errors in the modified energy H^(pn,qn)+12h4α2K(pn,qn). Looking at the error in the Hamiltonian H(pn, qn), which for this experiment is of size 10−3, no drift and no random walk behavior can be observed. However, the Takahashi-Imada method has been designed as to eliminate the An external file that holds a picture, illustration, etc.
Object name is nihms136293ig1.jpg(h2) term in the modified Hamiltonian (3.3) with help of the pre- and post-processing procedure χh of (1.5). Using the more accurate variables ([p with hat]n, qn) = χh(pn, qn), the random walk behavior will be present already in the original energy H([p with hat]n, qn).

Footnotes

*This material is based upon work supported by the Fonds National Suisse, project No. 200020-121561, by the Marsden Fund of the Royal Society of New Zealand, and by grant R01GM083605 from the National Institute of General Medical Sciences.

**This work emerged from discussions during the final workshop of the HOP programme on “highly oscillatory problems”, which was hosted by the Isaac Newton Institute in Cambridge.

4In general, the series (3.1) does not converge. To make backward error analysis rigorous, one has to truncate the series and one has to estimate the resulting errors. To emphasize the main ideas, we suppress these technicalities in the present work.

References

1. Blanes S, Casas F, Murua A. On the numerical integration of ordinary differential equations by processed methods. SIAM J Numer Anal. 2004;42:531–552.
2. Butcher JC. The effective order of Runge-Kutta methods. In: Morris JL, editor. Lect Notes Math; Proceedings of Conference on the Numerical Solution of Differential Equations; 1969. pp. 133–139.
3. Chartier P, Faou E, Murua A. An algebraic approach to invariant preserving integrators: the case of quadratic and Hamiltonian invariants. Numer Math. 2006;103:575–590.
4. Cottrell D, Tupper PF. Energy drift in molecular dynamics simulations. BIT. 2007;47:507–523.
5. Faou E, Hairer E, Pham TL. Energy conservation with non-symplectic methods: examples and counter-examples. BIT. 2004;44:699–709.
6. Hairer E, Lubich C. Symmetric multistep methods over long times. Numer Math. 2004;97:699–723.
7. Hairer E, Lubich C, Wanner G. Springer Series in Computational Mathematics. 2. Vol. 31. Springer-Verlag; Berlin: 2006. Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations.
8. McLachlan RI, Perlmutter M. Energy drift in reversible time integration. J Phys A. 2004;37:L593–L598.
9. Omelyan IP. Extrapolated gradientlike algorithms for molecular dynamics and celestial mechanics simulations. Phys Rev E. 2006;74:036703. [PubMed]
10. Rowlands G. A numerical algorithm for Hamiltonian systems. J Comput Phys. 1991;97:235–239.
11. Skeel RD, Zhang G, Schlick T. A family of symplectic integrators: stability, accuracy, and molecular dynamics applications. SIAM J Sci Comput. 1997;18:203–222.
12. Skeel RD. What makes molecular dynamics work? SIAM J Sci Comput. 2009;31:1363–1378. [PMC free article] [PubMed]
13. Stoffer D. Technical Report SAM-Report No. 88-05. ETH-Zürich; Switzerland: 1988. On reversible and canonical integration methods.
14. Takahashi M, Imada M. Monte Carlo calculation of quantum systems. II. Higher order correction. J Phys Soc Jpn. 1984;53:3765–3769.
15. Tupper PF. Ergodicity and the numerical simulation of Hamiltonian systems. SIAM J Appl Dyn Syst. 2005;4:563–587.
16. Wisdom J, Holman M, Touma J. Symplectic correctors. In: Marsden JE, Patrick GW, Shadwick WF, editors. Integration Algorithms and Classical Mechanics. Amer. Math. Soc.; Providence R.I.: 1996. pp. 217–244.