Numerical simulations for the Toda lattices Hamiltonian system: Higher order symplectic illustrative perspective

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).


Introduction
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 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, where q a and p a 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, q a ð0Þ ¼ q a 0 ; p a ð0Þ ¼ p a0 : We define Eq (2) in abbreviated form as, 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, The above equations of motion are equivalent to Newton's second law of mechanics with conservative forces. The dynamics generated by these equations define, for each evolution time, a mapping between regions of phase space. A general feature of these mappings is that they preserve the volume of the regions being mapped (and some related properties, collectively referred to as the symplectic structure). The standard non-symplectic numerical integrators, that have been used to solve general initial value problems numerically, do not preserve this qualitative behaviour, or the constants of motion for the system. Examples of such numerical integrators are the classical Runge-Kutta methods of different order, as found in standard integration packages.
By contrast the geometrical numerical integrators have gained popularity in the scientific community, due to their geometry preserving properties, in order to find qualitatively better solutions to Hamiltonian problems. In physical systems energy preservation, symmetries, time-reversal invariance, symplecticity, angular momentum, phase-space volume and dissipation are some key and crucial components to understand geometric properties. Detailed discussions of symplectic integrators with geometric properties have been given in [1][2][3][4].
Since the symplectic solvers have been widely accepted to be superior than the conventional numerical methods for solving the Hamiltonian systems, Mushtaq et. al. [5] constructed a well behaved class of higher order symplectic integrators schemes based on the extensions of the Störmer-Verlet scheme for Hamiltonians like Eq (1). An overview of these extensions are presented in Section 3. In this paper, we apply these schemes to the Toda lattice models.
The new proposed (KiMoKi) schemes involve extensive calculations of higher order derivatives of the Hamiltonian; hence it becomes a nightmare to do correct implementations manually. A collection of Python program HOMsPy (Higher Order Methods in Python) has been developed and presented by Mushtaq and Olaussen [8] to overcome these cumbersome and error-prone calculations for higher accuracy. More details, with implementation of many applications, can be found tutorial on HOMsPY by Mushtaq [9].
The structure of the rest of this paper is as follows: In Section 2 we review the Toda lattice models. These are integrable, nonlinear systems that have a number of extra constants of motions beyond standard ones like energy and momentum. The form of these can be described in a very consise manner by Eqs (7) and (8). In Section 3 we review the construction of the KiMoKi class of symplectic numerical solvers. In Section 4 we present and discuss the numerical simulations of the Toda lattice models by use of these methods. In Section 5 we conclude the main body of paper with brief remarks. Some technical details are delegated to appendices.

Toda lattices
The periodic Toda lattice with N sites (or particles) can, in suitable dimensionless coordinates, be specified by the Hamiltonian [10], where q a and p a are phase-space coordinates of positions and momenta respectively, and the index a is interpreted modulo N (i.e., q N þ1 � q 1 ). 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) = q 2 /2 − αq 3 /3 + βq 4 /4. The original study by Fermi et. al. [11] only treated linear chains with fixed endpoints, and parameter choices for which αβ = 0. Integrability is one of the most important properties of the Toda lattices. The model describes a set of equal mass particles moving on a ring with exponentially decreasing nearest neighbour interactions. The 2N phase-space coordinates can be used to define a symmetric, periodic tridiagonal N � N (time dependent) matrix, where v a = −e (q a −q a+1 )/2 . It was shown by Flaschka [12], using theory developed by Lax [13], that the eigenvalues λ a of L remain unchanged if the evolution q a (t), p a (t) is generated by the Hamiltonian of Eq (6). This means that all quantities are constants of motions. The first two are familiar, general expressions, The third one is a special consequence of integrability, where the last term is just an uninteresting constant, since v 1 v 2 v 3 = 1 when N ¼ 3.
It may not be easy to discover a general prescription like the one above. Alternative methods to find additional conservation laws (when one suspects that such exists) are the more brute force type of searches used by Göktaş et. al. [14] and Hohler et. al. [15]. They started by deducing the general form of a possible conservation law, with unknown coefficients, and next tried to explicitly solve for the coefficients with the help of computer algebra. This is a more pedestrian and cumbersome approach, but may be more likely to succeed when applied to models with unknown properties.

The 3-particle case
For N ¼ 3 it is simple to introduce center-of-mass and relative coordinates. A common physicists choice is Jacobi coordinates with corresponding conjugate momenta (see Appendix A), This separates the Hamiltonian into center-of-mass and internal contributions, In a similar manner we may rewrite C 3 ¼ 1 A direct evaluation of dC 3? /dt, using the Hamilton equations generated by H ? , confirms that it vanishes. I.e., that C 3? indeed is a constant of motion. The Hamiltonian of Eq (10a) can be rewritten by a canonical scale transformation, followed by a change of time and mass units (see Appendix B), t ! ffi ffi ffi 3 p t, m ! 1 8 m. This transforms Eq (10a) to the expression used by Lunsford and Ford [16], With the same transformations the conserved quantity of Eq (10b) can be expressed as If the potential in Eq (11) is expanded to third order in the coordinates q 1 , q 2 , one obtains the Henon-Heiles [17] Hamiltonian This is the motivation for the form by Lunsford and Ford [16]. The KiMoKi solvers have already been implemented for the Henon-Heiles model by Mushtaq [9], and used to detect chaotic and non-chaotic regions of that model. The former occur for oscillations of larger amplitude, for which the higher order terms in the Toda lattice Hamiltonian become of importance. This explains why chaotic behaviour may occur in the Henon-Heiles model, but not in the Toda lattice models (since they are integrable).
3 An overview to construct the higher order symplectic scheme One idea for construction of a symplectic integrator for the evolution generated by a Hamiltonian, 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 ! T eff and V ! V eff , 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], V eff 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 T eff ; it must depend on both p and q. However, what is really needed is not the infinitesimal generator T eff (p, q; τ), but its corresponding, sufficiently accurate, finite (but short) time generator G(q, P; τ). The latter can be constructed in a systematic manner: such that the transformation (q, p) ! (Q, P) is defined by which preserves the symplectic structure exactly, reproduces the time evolution generated by T eff to order τ N . Here Q a is shorten for q a (t + τ), and P a shorten for p a (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.
Define K such that 2K + 2 is the order of the method, where K = 1, 2, 3 and K = 0 corresponds to the Störmer-Verlet scheme. One full time step with this modification for kick-pushmove-kick scheme is, 2. "Push": But P a is unknown yet. We have to solve a nonlinear equation to find P a . However our generating function takes the form, Hence Eq (14a) can be written as, or in a form suitable for an iterative solution, 4. "Kick": The explicit expressions for V 2k and G k were published by Mushtaq et. al. in ref [5]. For convenience, on request from a reviewer, they are included in Appendix D.
To implement these higher order methods for our Toda lattice Hamiltonian, we define the model by the code in Listing 1 of Appendix C, and use programs in the HOMsPy package to automatically generate all the KiMoKi solver code. The complete code package is available as "Supporting information" for this paper.

Numerical simulations with HOMsPy
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.
As mentioned earlier that we implemented these numerical schemes as Python routines in HOMsPy. Python is an open source programming language which has gained increasing popularity in general, including (successful) applications for scientific computing. It is fast and easy to code and use for small "prototyping" tasks, since there is no need for explicit declaration of variables or a separate compilation cycle. It also comes with a huge repository of packages covering a large area of applications. Python is equipped with other features which facilitates development and encourages documentation of large well-structured program systems. Obviously, as an interpreted language, native Python is not suitable for performing extended numerical computations. But very often the code for such computations reduces to calls to pre-compiled library routines.

Preservation of constants of motion
We have already stressed the advantage of using symplectic solvers for numerical analysis of Hamiltonian models. This is most important when simulating long time series, where conventional numerical algorithms (like the Runge-Kutta methods usually implemented in numerical packages) can lead to a continuous degradation of important qualitative properties of Hamiltonian systems, like symplecticity (preservation of phase space volume and related quantities) and constants of motion. These methods have no built-in mechanisms for preserving such properties, as is illustrated in this subsection.
As the name suggests, symplectic solvers preserve symplecticity exactly. There will, of course, always be errors caused by numerical roundoff, but such errors do not depend on the accuracy of the method, only on the numerical precision being used. Symplectic solvers do not preserve most other constants of motion exactly, but the error (deviation from the initial value) will oscillate in a narrow band around zero. The width of this band scales with the accuracy of the method (i.e., order and timestep) in the expected way. For the symplectic (KiMoKi) algorithms used in this paper, applicable to Hamiltonians of the form H = T(p) + V(q), a constant of motion is preserved exactly if it is conserved separately by T and V. In this paper, one such example is the total momentum P of Eq (9a), while H of Eq (9b) and C 3 Eq (9c) are not.
It has been proven that symplectic integrators that preserve the Hamiltonian must actually be exact solvers (modulo errors introduced by finite numerical precision). There exist special methods for integrable models, as f.i. discussed by Kuznetsov and Vanhaecke [6] and Zullo [7]. The KiMoKi integrators, aimed for a more general class of problems, are not able to provide exact solutions, at least not when the conventional coordinates are used.
In this subsection we compare the KiMoKi solvers of order 2 and 4 with the RK23 (order 2) and RK45 (order 4) Runge-Kutta methods available through the solve_ivp routine in the scipy.integrate package, for the same values of the timestep τ (for the Runge-Kutta solvers, τ is the maximum timestep).
There are additional methods available in solve_ivp, but they are-for this comparisoninferior to the Runge-Kutta ones.
As can be seen from Figs 1 and 2, for short times the Runge-Kutta accuracy may very well be better than the KiMoKi ones, but as time increases it steadily becomes worse. By comparing Figs 1 and 2, we note that the time interval in which the Runge-Kutta methods are competitive becomes larger with decreasing τ.

Poincaré section technique (surface-of-section)
The reduced 3-site Toda lattice model we investigate has 4 degrees of freedom, z � (q 1 , q 2 , p 1 , p 2 ). Even in this rather simple case it is a challenge to present and visualize how the solutions behave. One, quite popular and efficient method, is the Poincaré section technique (also known as surface-of-section), introduced by Henri Poincaré the early 20th century. Generally, for cases where energy is conserved, H(q 1 , q 2 , p 1 , p 2 ) = H 0 , each orbit is restricted to a 3-dimensional constant energy surface of 4-dimensional phase space. The points where one coordinate is kept fixed (for example q 2 = 0) define another, in general independent, 3-dimensional surface. The intersection of these two surfaces is therefore two-dimensional. It can be specified by two coordinates, for example (q 1 , p 1 ). In this example, the points (q 1 (t n ), p 1 (t n )) where the constant energy orbit crosses the q 2 = 0 surface are therefore easy to visualize in two-dimensional plots. The times t n of crossings, and the order of repeated crossings, will be lost (or visualized by other means).  Repeated crossings will generate a pattern which indicates the nature of the dynamics. In our case, where there is an additional constant of motion (C 3? ), repeated crossings will appear on two smooth curves-one for each direction in which the (q 2 = 0)-plane is crossed, determined by the initial conditions. For ergodic motion, the crossings should spread smoothly over or more regions of the plane, according to density predicted by classical statistical mechanics. According to KAM theory, perturbations of integrable models are expected to lie in-between: For a finite fraction of initial condition, the crossings will appear on a smooth curve, while the rest will appear to be spread over a region of finite area.
Cheb-Terrab and de Oliveira [18] have written a MapleV R.3 routine for visualizing Hamiltonian dynamics by the Poincaré section technique. They employ the Toda lattice model for usage demonstration. We have implemented their algorithm in python, in combination with the KiMoKi solvers. The algorithm uses linear interpolation to determine crossings; hence it is of limited accuracy and is best used with short time-steps τ. (All our code could have been implemented in Maple, but this framework is not freely available to all.)

3D camera plots of each orbit
An alternative method to visualize the solution behaviour is to make a projection to a 3-dimensional subspace, and display the orbit in a "3-dimensional" plot. This is best done in interactive sessions, since this allows one to vary the viewing direction over all possible spherical angles. Snapshots examples from such matplotlib sessions are shown in Figs 6 and 7, for phase space orbits {z(t)|0 � t � 2 000}. Each plot displays quasi-periodic motion on a two-dimensional surface determined by the initial value z 0 (or more precisely the corresponding constants of motion, H 0 and C 3? ).

Behavior of energy error
We have earlier in Section 4.1 and Figs 1 and 2 shown that the long time behaviour of the KiMoKi solvers are better than the standard Runge-Kutta solvers of the same order, with respect to preservation of constants of motion. In Fig 8 we show that this behaviour can be observed for all orders N of the KiMoKi solvers, with the accuracy increasing with N for a fixed timestep τ. As can be seen, the errors keep varying in an oscillating manner, with no noticeable increase in amplitude with time.
In Fig 9 we further show that the error scales with order N and timestep τ in the expected manner. I.e., proportional with τ N , with a N-dependent constant of proportionality. Higher order symplectic illustrative perspective

Concluding remarks
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.
A brief summary of this work is as follows: • The symplectic property is preserved provided we solve the non-linear Eq (14a) for push steps to sufficient accuracy.
• Without prior knowledge the quasi-periodic nature of the solutions can easily be detected from 3D plots of the orbits (projected to 3-dimensional subspaces). Further (but not independent) confirmation can be found by investigating the behaviour of the Poincaré section of each orbit.
• Although the KiMoKi solver do not preserve constants of motion exactly, the time oscillating error in these quantities do not systematically increase in "amplitude" with time. This amplitude can be reduced in a predictable manner by increasing the order N of the method, or decreasing the timestep τ, or both.

A Jacobi coordinates for few-body systems
Consider first a translation invariant Hamiltonian system with two non-relativistic particles of mass m 1 resp. m 2 , with position coordinates q 1 resp. q 2 . To exploit translation invariance it is Higher order symplectic illustrative perspective common to introduce center-of-mass and relative coordinates, where μ j = m j /(m 1 + m 2 ), for j = 1, 2. For common systems with conjugate momenta p 1 ¼ m 1 _ q 1 and p 2 ¼ m 2 _ q 2 , the new momenta become The linear transformation of Eq (16) is canonical (because the matrices M q in Eq (16a) and M p in Eq (16b) are related by M t p ¼ M À 1 q ) and maintains the diagonal form of the kinetic energy. In the equal-mass case, m 1 ¼ m 2 ¼ 1 2 , the matrices M q and M p 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.
The inverse of Eq (16) is The extension to three particles is obvious for the center-of-mass coordinate, and one may further maintain the previous definition of one relative coordinate. As a second relative coordinate, select the distance between the center-of mass of the first two particles, and the third one. Hence X x 1 where now μ j = m j /(m 1 + m 2 + m 3 ) for j = 1, 2, 3 andm j ¼ m j =ðm 1 þ m 2 Þ for j = 1, 2. The new conjugate momenta becomes respectively P π 1 The inverse of Eq (18) is For the case of equal masses, m j ¼ 1 3 andm j ¼ 1 2 .

B Unit transformations
Most quantities in physical expressions, like the Hamiltonian 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 t 0 , ℓ 0 , m 0 . I.e., a dimensionless potential energy V(q) = e (q a −q a + 1 ) must be interpreted to mean ðm 0 ' 2 0 =t 2 0 Þe ðq a À q aþ1 Þ=' 0 , and the factor 1 2 in the dimensionless kinetic energy must be interpreted to mean 1 2 m À 1 0 . In "units where t 0 = ℓ 0 = m 0 = 1". Consider now a change of reference units to with all physical quantities fixed. It is rather obvious that dimensionless coordinates will change as In this context, the statement t ! λ t t is shorthand for i) a change of reference units t 0 !t 0 ¼ t 0 =l t , implying ii) the transformation of Eq (21a), often followed by iii) a symbol renaming back to the original one,t ! t.
The corresponding transformations of V(q) and 1 2 p 2 are less obvious, and cannot be deduced from the dimensionless expressions without knowledge of which physical quantity they represent (energy in this case).

C Code snippets
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. The package itself can downloaded from the CPC Program Library at http://cpc.cs.qub.ac. uk/summaries/ADTV, by providing the Catalogue Id AESD v1.0. The code, with accompanying information which need not be repeated here, is found by unpacking the downloaded . tar-file. Since the package is written Python 2.7, we provide code snippets illustrating how the solvers can be accessed by Python 3 code.
One way to continue is to add the code in Listing By running this code several files will be created. The most important one is TodaLattices.py, which contains the KiMoKi solvers up to 8th order. This file should not be modified manually; is not intended to be studied in detail by humans. But (in particular) the multiprecision code should be checked against unintended conversions to floating point numbers. F.i., if the final division /24 in the above definition of V is replaced by a pre-multiplication (1/24) � , then this factor will be converted to a double precision number at an early stage, and thereby pollute all multiprecision accuracy.
The file runTodaLattices.py provide some usage examples, intended to be modified and extended. Since this file will be overwritten the next time makeTodaLattices is executed, it is recommended to work on a renamed copy.
The current version of the HOMsPY package, with all its output, is written in Python 2.7. But functions can be accessed from Python 3 through an interface like the one in Listing 2. Data exchange via pickle files may seem primitive, but this has the advantage of documenting (preserving) the arguments and data being used. Further, the last two blocks in run_todalattices.py were changed to the one of Listing 4. @@@Execute the function named by sys.argv [1], with argument sys.argv [2]. @@@ argc = len(sys.argv) func = sys.argv [1] if argc > 1 else @solve_ivp@ argfile = @%s.pkl@ % sys.argv [2] if argc > 2 else @args.pkl@ funcs[func](argfile)

D Explicit expressions
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]: Explicit (but compact) expressions for the terms of order τ 2 , τ 4 , and τ 6 were given in [5,19]. With the notation, @ a � @ @q a ; @ a � M ab @ b ; p a � M ab p b ; D � p a @ a ; � D � ð@ a VÞ@ a ; where the Einstein summation convention is employed (an index which occur twice, once in lower position and once in upper position, are implicitly summed over all available values; i.e, M ab @ b � ∑ b M ab @ b -we generally use the matrix M to rise an index from lower to upper position), they are In this last line we have introduced The kick-steps can still be integrated directly, since the V 2k 's only depend on q. However, the T 2k 's (for k � 1) in general depend on both q and p; hence the move-steps cannot be integrated directly. To overcome this problem we introduce a generating function Gðq; P; tÞ ¼ X 0�k�N G k ðq; PÞ t k ð24Þ such that the transformation (q, p) ! (Q, P) defined implicitly by preserves the symplectic structure exactly, and reproduce the time evolution generated by T eff to order τ N . Here Q a is shorthand for q a (t + τ), and P a shorthand for p a (t + τ). The explicit expressions for the coefficients G k are G 0 ¼ q a P a ; ð26aÞ The Eqs (22,26) define the kick-move-kick scheme for a general potential V. If one uses all the listed terms the local error becomes of order τ 9 , and the scheme will respect long-time conservation of energy to order τ 8 .