**Date Published:** April 18, 2019

**Publisher:** Public Library of Science

**Author(s):** Asif Mushtaq, Amna Noreen, Muhammad Asif Farooq, Joseph Páez Chávez.

http://doi.org/10.1371/journal.pone.0215054

**Abstract**

**In this paper we apply some higher order symplectic numerical methods to analyze the dynamics of 3-site Toda lattices (reduced to relative coordinates). We present benchmark numerical simulations that has been generated from the HOMsPY (Higher Order Methods in Python) library. These results provide detailed information of the underlying Hamiltonian system. These numerical simulations reinforce the claim that the symplectic numerical methods are highly accurate qualitatively and quantitatively when applied not only to Hamiltonian of the Toda lattices, but also to other physical models. Excepting exactly integrable models, these symplectic numerical schemes are superior, efficient, energy preserving and suitable for a long time integrations, unlike standard non-symplectic numerical methods which lacks preservation of energy (and other constants of motion, when such exist).**

**Partial Text**

Hamiltonian equations of motion belong to a class of ordinary differential equations (ODEs) which in general are difficult or mostly impossible to solve analytically. Consider a separable Hamiltonian written in the form

H(q,p)=12pTMp+V(q)=T(p)+V(q),(1)

where T(p) is the non-relativistic kinetic energy, V(q) is the potential energy and M is the inverse mass matrix. The autonomous Hamiltonian equations of motion constitute a system of first order ordinary differential equations,

q˙a=∂H∂pa,p˙a=−∂H∂qa,a=1,…N(2)

where qa and pa are generalized coordinates of positions and momenta, respectively. The ” ˙ ” denotes differentiation with respect to time t, and H = H(q, p). The initial conditions at t = 0 can be written as,

qa(0)=q0a,pa(0)=pa0.

We define Eq (2) in abbreviated form as,

z=(qp),∇H=(∂H/∂qa∂H/∂pa),(3)J=(0I−I0),(4)

where J is a skew-symmetric matrix. Further I and 0 represent the (N×N) unit and zero matrices, respectively. The compact conservative Hamiltonian system in differential form is,

z˙=J−1∇H(z).(5)

The periodic Toda lattice with N sites (or particles) can, in suitable dimensionless coordinates, be specified by the Hamiltonian [10],

H(q,p)=12∑a=1Npa2+∑a=1N(exp(qa−qa+1)−1),(6)

where qa and pa are phase-space coordinates of positions and momenta respectively, and the index a is interpreted modulo N (i.e., qN+1≡q1). This mode belongs to a more general class of lattice models where the nearest-neighbour potential, exp(q) − 1, is replaced by an arbitrary function V(q). Another famous member of this class is the Fermi-Pasta-Ulam-Tsingou problem, with V(q) = q2/2 − αq3/3 + βq4/4. The original study by Fermi et. al. [11] only treated linear chains with fixed endpoints, and parameter choices for which αβ = 0.

One idea for construction of a symplectic integrator for the evolution generated by a Hamiltonian,

H=T(p)+V(q),

is to replace it with an iterated sequence of short-time evolutions generated by respectively T(p) (moves, which changes q without changing p) and V(q) (kicks, which changes p without changing q), since each of these are exactly integrable. This is the Störmer-Verlet scheme, which in its symmetrized form has a global error scaling like the timestep squared, τ2. One way to achieve higher accuracy is by replacing T and V by effective quantities, T → Teff and V → Veff, in a systematic manner. The effective quantities will depend on the timestep τ, and the wanted order N of accuracy, τN. In the kick-push-move-kick scheme proposed by Mushtaq et. al. [5], Veff is still a function of q only (in addition to τ); hence it can still be treated a potential, only slightly changed. Then this is no longer possible for Teff; it must depend on both p and q. However, what is really needed is not the infinitesimal generator Teff(p, q; τ), but its corresponding, sufficiently accurate, finite (but short) time generator G(q, P; τ). The latter can be constructed in a systematic manner:

G(q,P)=∑k=0NGk(q,P)τk,(13)

such that the transformation (q, p) → (Q, P) is defined by

pa=∂G∂qa,(14a)Qa=∂G∂Pa,(14b)

which preserves the symplectic structure exactly, reproduces the time evolution generated by Teff to order τN. Here Qa is shorten for qa(t + τ), and Pa shorten for pa(t + τ). The change in momentum p (of order τ3—i.e. a gentle push) is then defined through an implicit equation (but one which has turned out to be unproblematic to solve by iteration for all cases tried), while the change in position q continues to be explicit. Hence, the evolution step generated by G consists of a move, accompanied with a gentle push.

As has been mentioned before, numerical simulations on several Hamiltonian systems with the algorithm outlined by Eq (15) has compared favourable to conventional non-symplectic methods. We here present the results of additional simulations, of the model defined by Eq (11), which strengthen this evidence further.

In this paper, the KiMoKi algorithms for numerical solutions of the Hamilton equations for a Toda lattice model have been discussed and tested. These methods preserve the symplectic structure exactly (within the accuracy given by the employed numerical precision); For order N = 2 the method is equal to the Störmer-Verlet scheme, with long-time accuracy of order τ2; it has been extended to methods of order τ4, τ6 and τ8. As demonstrated, the method works as expected (sometimes even better than expected) for the reduced 3-site Toda lattice model.

Consider first a translation invariant Hamiltonian system with two non-relativistic particles of mass m1 resp. m2, with position coordinates q1 resp. q2. To exploit translation invariance it is common to introduce center-of-mass and relative coordinates,

(Xx)=(μ1μ21−1)(q1q2),(16a)

where μj = mj /(m1 + m2), for j = 1, 2. For common systems with conjugate momenta p1=m1q˙1 and p2=m2q˙2, the new momenta become

(Pπ)=(11μ2−μ1)(p1p2).(16b)

The linear transformation of Eq (16) is canonical (because the matrices Mq in Eq (16a) and Mp in Eq (16b) are related by Mpt=Mq−1) and maintains the diagonal form of the kinetic energy. In the equal-mass case, μ1=μ2=12, the matrices Mq and Mp do not become orthogonal. The latter, which can be obtained by scale transformations of X and x, may look simpler and more natural from a mathematical point of view. However, this would obscure physical interpretation of the coordinates.

Most quantities in physical expressions, like the Hamiltonian

H=12mp2+V(q),

are dimensionful. I.e., they carry units of time, length, and mass. When expressed in dimensionless form like in Eq (6) or Eq (10a), this means that the dimensionless time t, length ℓ, and mass m actually are expressed in terms of some reference quantities t0, ℓ0, m0. I.e., a dimensionless potential energy V(q) = e(qa−qa + 1) must be interpreted to mean (m0ℓ02/t02)e(qa−qa+1)/ℓ0, and the factor 12 in the dimensionless kinetic energy must be interpreted to mean 12m0−1. In “units where t0 = ℓ0 = m0 = 1”. Consider now a change of reference units to

(t˜0,ℓ˜0,m˜0)=(t0/λt,ℓ0/λℓ,m0/λm),(20)

with all physical quantities fixed. It is rather obvious that dimensionless coordinates will change as

t→t˜=λtt,(21a)q→q˜=λqq,(21b)p→p˜=(λmλℓ/λt)p.(21c)

In this context, the statement t → λtt is shorthand for i) a change of reference units t0→t˜0=t0/λt, implying ii) the transformation of Eq (21a), often followed by iii) a symbol renaming back to the original one, t˜→t.

In this section we provide some information of how the routines in the HOMsPY package can be used to create the symplectic solvers for the Toda lattice Hamiltonian, and how these solvers can be used to solve an initial value problem from provided initial data.

On request from a reviewer we here for convenience include the explicit expressions used in the algorithms of Eq (15). The rest of this section is an essentially unedited copy of a section with the same name, previously published by Mushtaq and Olaussen [8]:

Source:

http://doi.org/10.1371/journal.pone.0215054