Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2881582

Formats

Article sections

- Abstract
- Introduction
- 1. Simplified Takahashi-Imada method
- 2. Review on energy conservation with symmetric methods
- 3. Modified differential equation
- 4. Energy conservation of the simplified Takahashi-Imada integrator
- 5. Numerical illustration
- References

Authors

Related links

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/2009019PMCID: PMC2881582

NIHMSID: NIHMS136293

Ernst Hairer: hc.eginu@reriaH.tsnrE; Robert I. McLachlan: zn.ca.yessam@nalhcaLcM.R; Robert D. Skeel: ude.eudrup.sc@leeks

See other articles in PMC that cite the published article.

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.

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.

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

$$H(p,q)={\scriptstyle \frac{1}{2}}{p}^{\top}p+U(q)$$

(1.1)

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

$$\begin{array}{l}\stackrel{.}{p}=f(q)\hfill \\ \stackrel{.}{q}=p\hfill \end{array}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{with}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}f(q)=-\nabla U(q).$$

(1.2)

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

For a differential equation = *f* (*q*), the replacement of the second derivative by a finite difference yields the most natural discretization *q _{n}*

$$\begin{array}{l}{p}_{n+1/2}={p}_{n}+{\scriptstyle \frac{h}{2}}f({q}_{n})\\ {q}_{n+1}={q}_{n}+h{p}_{n+1/2}\\ {p}_{n+1}={p}_{n+1/2}+{\scriptstyle \frac{h}{2}}f({q}_{n+1}).\end{array}$$

(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 (*h*^{3}), so that the global error in the solution as well as in the energy is (*h*^{2}).

Another way to formulate the Störmer-Verlet method as a one-step method is to introduce *q _{n}*

The accuracy of the Störmer-Verlet discretization is greatly improved by adding the expression *αh*^{2}*f*′(*q*)*f*(*q*) with
$\alpha ={\scriptstyle \frac{1}{12}}$ to every force evaluation in (1.3). This yields the method (*cf.* [10,14])

$$\begin{array}{l}{p}_{n+1/2}={p}_{n}+{\scriptstyle \frac{h}{2}}(I+{\alpha h}^{2}{f}^{\prime}({q}_{n}))f({q}_{n})\\ {q}_{n+1}={q}_{n}+h{p}_{n+1/2}\\ {p}_{n+1}={p}_{n+1/2}+{\scriptstyle \frac{h}{2}}(I+{\alpha h}^{2}{f}^{\prime}({q}_{n+1}))f({q}_{n+1}).\end{array}$$

(1.4)

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

$$V(q)=U(q)-{\scriptstyle \frac{\alpha}{2}}{h}^{2}{\left|\right|\nabla U(q)\left|\right|}^{2}.$$

Although it is still of classical order two, it is an effectively fourth order integrator for
$\alpha ={\scriptstyle \frac{1}{12}}$. Denoting the method (1.4) by Φ* _{h}*: (

$$p=\widehat{p}+c{h}^{2}{f}^{\prime}(q)\widehat{p},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\widehat{q}=q+c{h}^{2}f(q),$$

(1.5)

a straightforward computation by Taylor expansion shows that
${\chi}_{h}\circ {\mathrm{\Phi}}_{h}\circ {\chi}_{h}^{-1}$ is a method of order three if and only if
$c=\alpha ={\scriptstyle \frac{1}{12}}$. Since this composed method is symmetric (which follows from *χ*_{−}* _{h}* =

Notice that in a constant step size implementation ${({\chi}_{h}\circ {\mathrm{\Phi}}_{h}\circ {\chi}_{h}^{-1})}^{n}={\chi}_{h}\circ {\mathrm{\Phi}}_{h}^{n}\circ {\chi}_{h}^{-1}$ 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]).

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

$$\begin{array}{l}{p}_{n+1/2}={p}_{n}+{\scriptstyle \frac{h}{2}}f({q}_{n}+{\alpha h}^{2}f({q}_{n}))\\ {q}_{n+1}={q}_{n}+h{p}_{n+1/2}\\ {p}_{n+1}={p}_{n+1/2}+{\scriptstyle \frac{h}{2}}f({q}_{n+1}+{\alpha h}^{2}f({q}_{n+1}))\end{array}$$

(1.6)

which is a (*h*^{5}) 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*) (*p* + *a*(*q*)*, q*) and (*p, q*) (*p, q* + *hp*)), and shears always preserve the volume.

On the other hand, a shear (*p, q*) (*p* + *a*(*q*)*, q*) is symplectic only if *a*′(*q*) is symmetric. Thus the first step of (1.6) is symplectic only if *f*′(*q _{n}* +

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.

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)={\scriptstyle \frac{1}{2}}{\omega}^{2}{q}^{2}$. For the Störmer-Verlet method we get the recursion *q _{n}*

$${p}_{n}^{2}+\left(1-{\scriptstyle \frac{{(\omega h)}^{2}}{4}}\right){\omega}^{2}{q}_{n}^{2}=\text{Const}.$$

that is exactly preserved. The maximal relative error in the Hamiltonian equals (*ωh*)^{2}/(4 − (*ωh*)^{2}) and it is attained when *p _{n}* = 0. This error behaves like (

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

$${q}_{n+1}-(2-{(\omega h)}^{2}+\alpha {(\omega h)}^{4}){q}_{n}+{q}_{n-1}=0.$$

For $\alpha ={\scriptstyle \frac{1}{12}}$ this recursion is stable if and only if $\omega h<2\sqrt{3}$. It also has a discrete energy which is

$${p}_{n}^{2}+\left(1-\beta {\scriptstyle \frac{{(\omega h)}^{2}}{4}}\right)\beta {\omega}^{2}{q}_{n}^{2}=\text{Const}.\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{with}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\beta =1-{\scriptstyle \frac{{(\omega h)}^{2}}{12}}.$$

As for the Störmer-Verlet method the maximal relative error in the Hamiltonian is of size (*h*^{2}) for *ωh* → 0, and tends to infinity when
$\omega h\to 2\sqrt{3}$. However, when we consider the transformed variables of (1.5), * _{n}* =

$${\widehat{p}}_{n}^{2}+\left(1-\beta {\scriptstyle \frac{{(\omega h)}^{2}}{4}}\right){\beta}^{-3}{\omega}^{2}{\widehat{q}}_{n}^{2}=\text{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 (*h*^{6}) for small step sizes and remains bounded (converges to one) when
$\omega h\to 2\sqrt{3}$.

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 *ωh*_{eff}, where the effective step size *h*_{eff} is *h* for the Takahashi-Imada integrator and 2*h* 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

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.

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

$$\begin{array}{l}\stackrel{.}{u}=f(u,v)\hfill \\ \stackrel{.}{v}=g(u,v)\hfill \end{array}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{where}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\begin{array}{l}f(u,-v)=-f(u,v)\hfill \\ g(u,-v)=g(u,v).\hfill \end{array}$$

(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*) (*u,* −*v*) it satisfies
$\rho \circ {\varphi}_{t}\circ \rho ={\varphi}_{t}^{-1}$. Many symmetric methods, and in particular also the simplified Takahashi-Imada method, are reversible and thus satisfy
$\rho \circ {\mathrm{\Phi}}_{h}\circ \rho ={\mathrm{\Phi}}_{h}^{-1}$.

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

$$\stackrel{.}{a}=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\stackrel{.}{\theta}=\omega (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}* (

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 (*th ^{r}*) when

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

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 = −*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 *r*th order method the energy error is (*h ^{r}*) on time intervals of length

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

$${q}_{n+1}-2{q}_{n}+{q}_{n-1}=\frac{{h}^{2}}{12}\left(f({q}_{n+1})+10f({q}_{n})+f({q}_{n-1})\right).$$

With suitable velocity approximations *p _{n}* it is conjugate to a symplectic integrator and thus has an excellent long-time behavior (see [11]).

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
${\chi}_{h}\circ {\mathrm{\Phi}}_{h}\circ {\chi}_{h}^{-1}-{\mathrm{\Psi}}_{h}=\mathcal{O}({h}^{s+1})$. If the method is of order

$$H({p}_{n},{q}_{n})=\text{Const}.+\mathcal{O}({h}^{r})+\mathcal{O}({th}^{s}),$$

(2.2)

so that the energy error is bounded by (*h ^{r}*) on intervals of length

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.

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 (*h ^{r}*

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

$$\left(\begin{array}{c}\stackrel{.}{p}\\ \stackrel{.}{q}\end{array}\right)=\left(\begin{array}{c}f(q)\\ p\end{array}\right)+h\left(\begin{array}{c}{f}_{2}(p,q)\\ {g}_{2}(p,q)\end{array}\right)+{h}^{2}\left(\begin{array}{c}{f}_{3}(p,q)\\ {g}_{3}(p,q)\end{array}\right)+\dots $$

(3.1)

Denoting by (*t*), (*t*) the solution of this system, we formally^{4} have for the numerical solution that *p _{n}* = (

- 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

$$\stackrel{\sim}{H}(p,q)={\scriptstyle \frac{1}{2}}{p}^{\top}p+U(q)+h{H}_{2}(p,q)+{h}^{2}{H}_{3}(p,q)+{h}^{3}{H}_{4}(p,q)+\dots $$

such that the modified differential equation is given by

$$\stackrel{.}{p}=-{\nabla}_{q}\stackrel{\sim}{H}(p,q),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\stackrel{.}{q}={\nabla}_{p}\stackrel{\sim}{H}(p,q).$$

Notice that the functions *f _{j}* (

The Störmer-Verlet method is symplectic, so that its modified differential equation is Hamiltonian. Using the notation
$T(p)={\scriptstyle \frac{1}{2}}{p}^{\top}p$ for the kinetic energy, and {*A, B*} =* _{q}A^{T}_{p}B* −

$$\stackrel{\sim}{H}=T+U+{h}^{2}\left({\scriptstyle \frac{1}{12}}\{T,\{T,U\}\}-{\scriptstyle \frac{1}{24}}\{U,\{U,T\}\}\right)+{h}^{4}\left(-{\scriptstyle \frac{1}{720}}\{T,\{T,\{T,\{T,U\}\}\}\}+{\scriptstyle \frac{1}{360}}\{U,\{T,\{T,\{T,U\}\}\}\}-{\scriptstyle \frac{1}{480}}\{U,\{U,\{T,\{T,U\}\}\}\}+{\scriptstyle \frac{1}{120}}\{T,\{T,\{U,\{U,T\}\}\}\}\right)+\mathcal{O}({h}^{6}).$$

(3.2)

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

Substituting *U* in (3.2) with
$U-{\scriptstyle \frac{\alpha}{2}}{h}^{2}{\left|\right|\nabla U\left|\right|}^{2}=U-{\scriptstyle \frac{\alpha}{2}}{h}^{2}\{U,\{U,T\}\}$ (*cf.* Sect. 1.2) yields the modified Hamiltonian of the symplectic Takahashi-Imada integrator (1.4):

$$\widehat{H}=T+U+{h}^{2}\left({\scriptstyle \frac{1}{12}}\{T,\{T,U\}\}-\left({\scriptstyle \frac{1}{24}}+{\scriptstyle \frac{\alpha}{2}}\right)\{U,\{U,T\}\}\right)+{h}^{4}\left(-{\scriptstyle \frac{1}{720}}\{T,\{T,\{T,\{T,U\}\}\}\}+{\scriptstyle \frac{1}{360}}\{U,\{T,\{T,\{T,U\}\}\}\}-{\scriptstyle \frac{1}{480}}\{U,\{U,\{T,\{T,U\}\}\}\}+\left({\scriptstyle \frac{1}{120}}-{\scriptstyle \frac{\alpha}{24}}\right)\{T,\{T,\{U,\{U,T\}\}\}\}-{\scriptstyle \frac{\alpha}{24}}\{U,\{T,\{U,\{U,T\}\}\}\}\right)+\mathcal{O}({h}^{6}).$$

(3.3)

Let (*p _{n}, q_{n}*) denote the numerical solution of (1.4). With the values
${\widehat{p}}_{n}={p}_{n}+{\scriptstyle \frac{1}{12}}{h}^{2}{\nabla}^{2}U({q}_{n}){p}_{n}$ and
${\widehat{q}}_{n}={q}_{n}-\left({\scriptstyle \frac{1}{12}}+{\scriptstyle \frac{\alpha}{2}}\right){h}^{2}\nabla U({q}_{n})$ we have

The choice
$\alpha =-{\scriptstyle \frac{1}{12}}$ gives * _{n}* =

To find the modified differential equation of (1.6) we notice that it is a (*h*^{5}) perturbation of the Takahashi-Imada method. If we denote by _{n}_{+1}_{/}_{2}, _{n}_{+1}, _{n}_{+1} the values obtained from (1.4) subject to initial values *q _{n}, p_{n}*, then the numerical approximation obtained from (1.6) satisfies with

$$\begin{array}{l}{p}_{n+1}={\widehat{p}}_{n+1}+{\scriptstyle \frac{1}{2}}{h}^{5}{\alpha}^{2}({f}^{\u2033}(f,f))({q}_{n+1/2})+\mathcal{O}({h}^{7})\\ {q}_{n+1}={\widehat{q}}_{n+1}+{\scriptstyle \frac{1}{4}}{h}^{6}{\alpha}^{2}({f}^{\u2033}(f,f))({q}_{n+1/2})+\mathcal{O}({h}^{8}).\end{array}$$

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

$$\begin{array}{l}\stackrel{.}{p}=-{\nabla}_{q}\widehat{H}(p,q)+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}({f}^{\u2033}(f,f))(q)+\mathcal{O}({h}^{6})\\ \stackrel{.}{q}={\nabla}_{p}\widehat{H}(p,q)+\mathcal{O}({h}^{6})\end{array}$$

(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*.

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

$${\scriptstyle \frac{\text{d}}{\text{d}t}}\widehat{H}(p,q)={\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}{p}^{\top}{f}^{\u2033}(f,f)+\mathcal{O}({h}^{6})=-{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}{U}_{\mathit{qqq}}(p,{U}_{q},{U}_{q})+\mathcal{O}({h}^{6}),$$

(4.1)

where we have used *f*(*q*) = −*U* (*q*) = −*U _{q}* (

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

$$\begin{array}{c}{\scriptstyle \frac{\text{d}}{\text{d}t}}\left(2{U}_{qq}({U}_{q},{U}_{q})-4{U}_{\mathit{qqq}}(p,p,{U}_{q})-{U}_{\mathit{qqqq}}(p,p,p,p)+2{({U}_{qq}p)}^{\top}({U}_{qq}p)\right)\\ =10{U}_{\mathit{qqq}}(p,{U}_{q},{U}_{q})-{U}_{\mathit{qqqqq}}(p,p,p,p,p)\end{array}$$

which holds along solutions of the Hamiltonian system (1.2). A (*h*^{2}) 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 *U _{qqq}* (

$${U}_{\mathit{qqq}}(p,{U}_{q},{U}_{q})={\scriptstyle \frac{\text{d}}{\text{d}t}}K(p,q)+{\scriptstyle \frac{1}{10}}{U}_{\mathit{qqqqq}}(p,p,p,p,p)+\mathcal{O}({h}^{2}),$$

(4.2)

where

$$K(p,q)={\scriptstyle \frac{1}{5}}{U}_{qq}({U}_{q},{U}_{q})-{\scriptstyle \frac{2}{5}}{U}_{\mathit{qqq}}(p,p,{U}_{q})-{\scriptstyle \frac{1}{10}}{U}_{\mathit{qqqq}}(p,p,p,p)+{\scriptstyle \frac{1}{5}}{({U}_{qq}p)}^{\top}({U}_{qq}p).$$

It therefore follows from (4.1) that

$${\scriptstyle \frac{\text{d}}{\text{d}t}}\left(\widehat{H}(p,q)+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K(p,q)\right)=-{\scriptstyle \frac{1}{20}}{h}^{4}{\alpha}^{2}{U}_{\mathit{qqqqq}}(p,p,p,p,p)+\mathcal{O}({h}^{6})$$

(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 *t _{n}* =

*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*h*^{4}-term of the right-hand side of (4.3) vanishes, so that the Hamiltonian $H(p,q)={\scriptstyle \frac{1}{2}}{p}^{\top}p+U(q)$ satisfies*H*(*p*) =_{n}, q_{n}*H*_{0}+ (*h*^{2}) + (*th*^{6}) 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 (*h*^{2}) with lifetime*t*=_{n}*nh*≤ (*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 (*th*^{4}), so that the lifetime of energy conservation (with error of size (*h*^{2})) of the simplified Takahashi-Imada integrator is at least*t*=_{n}*nh*≤ (*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 havewhere$$\underset{t\to \infty}{lim}\frac{1}{t}{\int}_{0}^{t}F(p(s),q(s))\text{d}s={\int}_{A}F(x)\mu (\text{d}x),$$(4.4)*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 (*th*^{4}) 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 $\mathcal{O}(\sqrt{t}{h}^{4})$. This implies that the lifetime with an (*h*^{2}) error of the simplified Takahashi-Imada integrator is at least*t*=_{n}*nh*≤ (*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.

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.

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 (*h*^{2}) there is a linear drift of the form (*th*^{4}). Figure 2 shows the truncated modified energy
${\scriptstyle \frac{1}{2}}{p}^{\top}p+U(q)+{h}^{2}{H}_{3}(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.

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 *h*^{4} 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.

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

$$U({q}_{1},{q}_{2})={\scriptstyle \frac{1}{2}}\left({q}_{1}^{2}+{q}_{2}^{2}\right)+{q}_{1}^{2}{q}_{2}-\frac{1}{k}{q}_{2}^{k}$$

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.

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 *q*_{1}(0) = 0, *q*_{2}(0) = 0.2, *p*_{2}(0) = 0.3, and a positive *p*_{1}(0) is chosen such that the Hamiltonian is *H*_{0} = 0.125. Figure 3 shows the error in the modified energy
$\widehat{H}({p}_{n},{q}_{n})+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K({p}_{n},{q}_{n})$ (with *α* = 1/12, *cf.* Eq. (4.3)). It behaves like (*h*^{6}) 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.

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 *q*_{1}(0) = 0, *q*_{2}(0) = 0.2, *p*_{2}(0) = 0.3, and *p*_{1}(0) such that the value of the Hamiltonian is *H*_{0} = 1/8 (implying a chaotic solution). Figure 4 shows the modified energy
$\widehat{H}({p}_{n},{q}_{n})+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K({p}_{n},{q}_{n})$ 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 × 10^{6}) and the standard deviation (*σ* = 0.54 × 10^{−6} at *t* = 2.5 × 10^{6}) 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.

Hénon-Heiles problem with quintic potential: error in the modified energy
$\widehat{H}({p}_{n},{q}_{n})+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K({p}_{n},{q}_{n})$ 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 *q*_{1} = 0 in the positive direction *p*_{1} > 0, we consider the point (*q*_{2}*, p*_{2}). 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 *p*_{2} < 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.

Poincaré cut for the numerical solution of the simplified Takahashi-Imada method corresponding to the interval [1.0 × 10^{6}*,* 1.1 × 10^{6}] 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.

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

$$U({q}_{1},\dots ,{q}_{N})=\sum _{1\le j<i\le N}V({r}_{ij}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{r}_{ij}=\left|\right|{q}_{i}-{q}_{j}\left|\right|.$$

In our experiment we assume that the positions *q _{i}* of the particles are in

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
$\widehat{H}({p}_{n},{q}_{n})+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K({p}_{n},{q}_{n})$ 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.

Error in the modified energy
$\widehat{H}({p}_{n},{q}_{n})+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K({p}_{n},{q}_{n})$ 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 (*α*^{2}*h*^{4}*t*^{1/2}).

We emphasize that in Figure 6 we have plotted the errors in the modified energy
$\widehat{H}({p}_{n},{q}_{n})+{\scriptstyle \frac{1}{2}}{h}^{4}{\alpha}^{2}K({p}_{n},{q}_{n})$. Looking at the error in the Hamiltonian *H*(*p _{n}, q_{n}*), which for this experiment is of size 10

^{*}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.

^{4}In 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.

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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |