
Citation: Moshe Rosenberg. Journal summary 2018 from Editor in Chief[J]. AIMS Agriculture and Food, 2019, 4(1): 163-164. doi: 10.3934/agrfood.2019.1.163
[1] | Biao Tang, Weike Zhou, Yanni Xiao, Jianhong Wu . Implication of sexual transmission of Zika on dengue and Zika outbreaks. Mathematical Biosciences and Engineering, 2019, 16(5): 5092-5113. doi: 10.3934/mbe.2019256 |
[2] | Georgi Kapitanov, Christina Alvey, Katia Vogt-Geisse, Zhilan Feng . An age-structured model for the coupled dynamics of HIV and HSV-2. Mathematical Biosciences and Engineering, 2015, 12(4): 803-840. doi: 10.3934/mbe.2015.12.803 |
[3] | Christopher M. Kribs-Zaleta . Alternative transmission modes for Trypanosoma cruzi . Mathematical Biosciences and Engineering, 2010, 7(3): 657-673. doi: 10.3934/mbe.2010.7.657 |
[4] | Moatlhodi Kgosimore, Edward M. Lungu . The Effects of Vertical Transmission on the Spread of HIV/AIDS in the Presence of Treatment. Mathematical Biosciences and Engineering, 2006, 3(2): 297-312. doi: 10.3934/mbe.2006.3.297 |
[5] | Sara Y. Del Valle, J. M. Hyman, Nakul Chitnis . Mathematical models of contact patterns between age groups for predicting the spread of infectious diseases. Mathematical Biosciences and Engineering, 2013, 10(5&6): 1475-1497. doi: 10.3934/mbe.2013.10.1475 |
[6] | Muntaser Safan . Mathematical analysis of an SIR respiratory infection model with sex and gender disparity: special reference to influenza A. Mathematical Biosciences and Engineering, 2019, 16(4): 2613-2649. doi: 10.3934/mbe.2019131 |
[7] | Kelu Li, Junyuan Yang, Xuezhi Li . Effects of co-infection on vaccination behavior and disease propagation. Mathematical Biosciences and Engineering, 2022, 19(10): 10022-10036. doi: 10.3934/mbe.2022468 |
[8] | Brandy Rapatski, James Yorke . Modeling HIV outbreaks: The male to female prevalence ratio in the core population. Mathematical Biosciences and Engineering, 2009, 6(1): 135-143. doi: 10.3934/mbe.2009.6.135 |
[9] | Fabio Sanchez, Luis A. Barboza, Paola Vásquez . Parameter estimates of the 2016-2017 Zika outbreak in Costa Rica: An Approximate Bayesian Computation (ABC) approach. Mathematical Biosciences and Engineering, 2019, 16(4): 2738-2755. doi: 10.3934/mbe.2019136 |
[10] | Romulus Breban, Ian McGowan, Chad Topaz, Elissa J. Schwartz, Peter Anton, Sally Blower . Modeling the potential impact of rectal microbicides to reduce HIV transmission in bathhouses. Mathematical Biosciences and Engineering, 2006, 3(3): 459-466. doi: 10.3934/mbe.2006.3.459 |
A huge number of important problems in physics, astronomy, chemistry, etc. are modeled via sets of ordinary differential equations (ODEs) governed by the Hamiltonian formalism. Due to non-integrability, the investigation of the time evolution of these problems, and generally their properties, often rely solely on numerical techniques. As modern research requires numerical simulations to be pushed to their very limit (e.g., large integration times, macroscopic limits), a methodical assessment of the properties of different numerical methods becomes a fundamental issue. Such studies allow to highlight the most suitable scheme for both general purposes and specific tasks, according to criteria of stability, accuracy, simplicity and efficiency.
The beginning of the era of computational physics is considered to be the computer experiment performed in the 1950s by Fermi, Pasta, Ulam and Tsingou (FPUT) [1,2,3] to observe energy equipartition due to weak-anharmonic coupling in a set of harmonic oscillators. Indeed, the breaking of integrability in Hamiltonian systems is often performed with the introduction of nonlinear terms in the equations of motion. These additional nonlinear terms describe new physical processes, and led to important questions and significant advancements in condensed matter physics. For instance, the FPUT problem has been used to answer questions related to ergodicity, dynamical thermalization and chaos occurrence (see e.g., [4,5,6] and references therein) and led to the observation of solitons [7,8], and progress in Hamiltonian chaos [9]. The interested reader can find a concise review of numerical results concerning the FPUT system in [10].
Another important example concerns disordered media. In 1958, Anderson [11] theoretically showed the appearance of energy localization in one-dimensional lattices with uncorrelated disorder (a phenomenon which is now referred to as Anderson localization). This phenomenon was later investigated also for two-dimensional lattices [12]. An important question which has attracted wide attention in theory, numerical simulations and experiments is what happens when nonlinearity (which appears naturally in real world problems) is introduced in a disordered system [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30].
Two basic Hamiltonian models are at the center of these investigations: the disordered Klein-Gordon (DKG) chain and the disordered, discrete nonlinear Schrödinger equation (DDNLS). In Refs. [22,23,24,25,26,28,29,30,31,32,33,34] it was shown that Anderson localization is destroyed by the presence of nonlinearity, resulting to the subdiffusive spreading of energy due to deterministic chaos. Such results rely on the accurate, long time integration of the systems' equations of motion and variational equations. We note that the variational equations are a set of ODEs describing the time evolution of small deviations from a given trajectory, something which is needed for the computation of chaos indicators like the maximum Lyapunov exponent (mLE) [35,36,37]. The numerical integration of these sets of ODEs can be performed by any general purpose integrator. For example Runge-Kutta family schemes are still used [38,39]. Another category of integrators is the so-called symplectic integrators (SIs), which were particularly designed for Hamiltonian systems (see e.g., [40,41,42,43,44] and references therein). The main characteristic of SIs is that they conserve the symplectic structure of Hamiltonian systems, so that the numerical approximation they provide corresponds to the exact solution of a system which can be considered as a perturbation of the original one. Consequently the error in the computed value of the Hamiltonian function (usually refereed to as the system's energy) is kept almost constant over long integration times. This almost constant error can be used as an indicator of the accuracy of the used integration scheme.
SIs have been successfully implemented for the long time integration of multidimensional Hamiltonian systems like the α- and β-FPUT systems and FPUT-like models, the DKG, the DDNLS and systems of coupled rotors (see e.g., [22,23,24,25,26,28,29,30,32,33,45,46,47,48,49,50,51,52,53,54]). In these studies SIs belonging to the so-called SABA family [55] were mostly used, since the FPUT, the DKG and systems of coupled rotors can be split into two integrable parts (the system's kinetic and the potential energy), while in some more recent studies [33,56] the ABA864 [57] SI was implemented as it displayed an even better performance. As the DDNLS model is not explicitly decomposed into two integrable parts, the implementation of SABA schemes requires a special treatment involving the application of fast Fourier transforms which are computationally time consuming [29]. In [49,50] it was shown that for the DDNLS model, the split of the Hamiltonian function into three integrable parts is possible, and this approach proved to be the best choice for the model's integration with respect to both speed and accuracy. It is worth noting here that the numerical integration of the variational equations by SIs is done by the so-called Tangent Map Method [58,59,60].
The intention of this work is to present several efficient (symplectic and non-symplectic) integration techniques for the time evolution of both the equations of motion and the variational equations of multidimensional Hamiltonian systems. In particular, in our study we use these techniques to investigate the chaotic dynamics of the one-dimensional (1D) FPUT system, as well as one- and two-dimensional (2D) DDNLS models. We carry out this task by considering particular initial conditions for both the systems and for the deviation vectors, which we evolve in time requiring a predefined level of energy accuracy. Then, we record the needed CPU time to perform these numerical simulations and highlight those methods that ensure the smallest CPU times for obtaining accurate results. A similar analysis for the 1D and 2D DKG systems was performed in [56].
The paper is organized as follows. In Section 2 we introduce the models considered in our study. In Section 3 we introduce in detail the symplectic and non-symplectic schemes implemented in our investigations. In Section 4 we present our numerical results and report on the efficiency of each studied scheme. In Section 5 we review our findings and discuss the outcomes of our study. Finally, in Appendix B we provide the explicit forms of the used tangent map operators. Moreover, we present there the explicit expressions of the tangent map operators for some additional important many-body systems, the so-called β-FPUT chain, the KG system and the classical XY model (a Josephson junction chain-JJC) in order to facilitate the potential application of SIs for these systems from the interested reader, although we do not present numerical results for these models.
In this work we focus our analysis on the 1D FPUT system and the 1D and 2D DDNLS models. In what follows we briefly present the Hamiltonian functions of these systems.
The 1D FPUT model [1,2,3] was first introduced in 1955 to study the road toward equilibrium of a chain of harmonic oscillators in the presence of weak anharmonic interactions. Since then, this system has been widely used as a toy model for investigating energy equipartition and chaos in nonlinear lattices. In the literature there exist two types of the FPUT model the so-called α- and β-FPUT systems. In our study we consider the α-FPUT model whose Hamiltonian function reads
H1F=N∑i=0[p2i2+12(qi+1−qi)2+α3(qi+1−qi)3]. | (2.1) |
In Eq.(2.1), qi and pi are respectively the generalized coordinates and momenta of the i lattice site and α is a real positive parameter. In our study we consider fixed boundary conditions q0=p0=pN+1=qN+1=0. We also note that this model conserves the value of the total energy H1F.
In contrast to the α-FPUT system, the β-FPUT model is characterized by a quartic nonlinear term [see Eq. (A.26) in Appendix B.4]. The fact that the value of the Hamiltonian function of the β-FPUT model is bounded from below leads to differences between the two models. For example, phenomena of chopping time which occur in the α- model do not appear in the β-FPUT system [61]. Further discussion of the differences between the α- and the β- models can be found in [62,63,64].
The DDNLS model describes anharmonic interactions between oscillators in disordered media and has been used to study the propagation of light in optical media or Bose-Einstein condensates through the famous Gross-Pitaevskii equation [52], as well as investigate, at a first approximation, several physical processes (e.g., electron tight binding systems [39]). The Hamiltonian function of the 1D DDNLS model reads
H1D=N∑i=1[ϵi2(q2i+p2i)+β8(q2i+p2i)2−pi+1pi−qi+1qi], | (2.2) |
where qi and pi are respectively the generalized coordinates and momenta of the i lattice site and the onsite energy terms ϵi are uncorrelated numbers uniformly distributed in the interval [−W2,W2]. The real, positive numbers W and β denote the disorder and the nonlinearity strength respectively. We also consider here fixed boundary conditions i.e., q0=p0=pN+1=qN+1=0.
The two-dimensional version of the DDNLS model was considered in [30,65]. Its Hamiltonian function is
H2D=N∑i=1M∑j=1{ϵi,j2[q2i,j+p2i,j]+β8[q2i,j+p2i,j]2−[qi,j+1qi,j+qi+1,jqi,j+pi,j+1pi,j+pi+1,jpi,j]}, | (2.3) |
where qi,j and pi,j are respectively the generalized positions and momenta at site (i,j) and ϵi,j are the disorder parameters uniformly chosen in the interval [−W2,W2]. We again consider fixed boundary conditions i.e., q0,j=p0,j=qN+1,j=pN+1,j=0 for 1≤j≤M and qi,0=pi,0=qi,M+1=pi,M+1=0 for 1≤i≤N.
Additionally to the energies H1D [Eq. (2.2)] and H2D [Eq. (2.3)] both systems conserve their respective norms S1D and S2D:
S1D=12N∑i=1(q2i+p2i) ;S2D=12N∑i=1M∑j=1(q2i,j+p2i,j) . | (2.4) |
The Hamilton equations of motion
dqdt=∂H∂p ,dpdt=−∂H∂q, | (3.1) |
of the N degree of freedom (dof) Hamiltonian H=H(q,p), with q=(q1,q2,…,qN) and p=(p1,p2,…,pN) being respectively the system's generalized positions and momenta, can be expressed in the general setting of first order ODEs as
dxdt=˙x=J2N⋅DH(x(t)), | (3.2) |
where x=(q,p)=(x1,x2,…,xN,xN+1,…,x2N)=(q1,q2,…,qN,p1,p2,…,pN) is a vector representing the position of the system in its phase space and (˙) denotes differentiation with respect to time t. In Eq. (3.2)
J2N=[ONIN−INON], | (3.3) |
is the symplectic matrix with IN and ON being the N×N identity and the null matrices respectively, and
DH=[∂H∂q1,…,∂H∂qN,∂H∂p1,…,∂H∂pN]T, | (3.4) |
with (T) denoting the transpose matrix.
The variational equations (see for example [37,58]) govern the time evolution of a small perturbation w(t) to the trajectory x(t) with w(t)=(δq(t),δp(t))=(δq1(t),δq2(t),…,δqN(t),δp1(t),δp2(t),…,δpN(t)) (which can also be written as w(t)=δx(t)=(δx1(t),…,δxN(t),δxN+1(t),…,δx2N(t))) and have the following form
˙w(t)=[J2N⋅D2H(x(t))]⋅w(t), | (3.5) |
where
[D2H(x(t))]i,j=∂2H∂xi∂xj|x(t),i,j=1,2,…,N, | (3.6) |
are the 2N×2N elements of the Hessian matrix D2H(x(t)) of the Hamiltonian function H computed on the phase space trajectory x(t). Eq. (3.5) is linear in w(t), with coefficients depending on the system's trajectory x(t). Therefore, one has to integrate the variational equations (3.5) along with the equations of motion (3.2), which means to evolve in time the general vector X(t)=(x(t),δx(t)) by solving the system
˙X=(˙x(t),˙δx(t))=f(X)=[J2N⋅DH(x(t))[J2N⋅D2H(x(t))]⋅δx(t)]. | (3.7) |
In what follows we will briefly describe several numerical schemes for integrating the set of equations (3.7).
Here we present the non-symplectic schemes we will use in this work: the Taylor series method and the Runge-Kutta discretization scheme. These methods are referred to be non-symplectic because they do not preserve the geometrical properties of Hamiltonian equations of motion (e.g., their symplectic structure). The immediate consequence of that is that they do not conserve constants of motion (e.g., the system's energy).
The Taylor series method consists in expanding the solution X(t) at time t0+τ in a Taylor series of X(t) at t=t0
X(t0+τ)=X(t0)+τ11!dX(t0)dt+τ22!d2X(t0)dt2+…+τnn!dnX(t0)dtn+O(τn+1(n+1)!dn+1X(t0)dtn+1). | (3.8) |
Once the solution X at time t0+τ is approximated, X(t0+τ) is considered as the new initial condition and the procedure is repeated to approximate X(t0+2τ) and so on so forth. Further information on this integrator can be found in [66,Sec. Ⅰ.8]. If we consider in Eq. (3.8) only the first n+1 terms we then account the resulting scheme to be of order n. In addition, in order to explicitly express this numerical scheme, one has to perform n−1 successive differentiations of the field vector f. This task becomes very elaborate for complex structured vector fields. One can therefore rely on successive differentiation to automate the whole process (see e.g., [67]). For the simulations reported in this paper we used the software package called TIDES [67,68] which is freely available [69]. We particularly focused on the implementation of the TIDES package as it has been already used in studies of lattice dynamics [60]. The TIDES package comes as a Mathematica notebook in which the user provides the Hamiltonian function, the potential energy or the set of ODEs themselves. It then produces FORTRAN (or C) codes which can be compiled directly by any available compiler producing the appropriate executable programs. In addition, the TIDES package allows us to choose both the integration time step τ and the desired 'one-step' precision of the integrator δ. In practice, the integration time step is accepted if the estimated local truncation error is smaller than δ.
It is worth noting that there exists an equivalent numerical scheme to the Taylor series method, derived from Lie series [70] (for more details see e.g., the appendix of [71]). Indeed, Eq. (3.7) can be expressed as
dXdt=LHZX, | (3.9) |
where LHZ is the Lie operator [72] defined as
LHZ=2N∑i=1(dxidt∂∂xi+dδxidt∂∂δxi). | (3.10) |
The formal solution of Eq. (3.9) reads
X(t0+τ)=eτLHZX(t0), | (3.11) |
and can be expanded as
X(t0+τ)=L0HZX(t0)+τ11!L1HZX(t0)+τ22!L2HZX(t0)+…+τnn!LnHZX(t0)+O(τn+1(n+1)!Ln+1HZX(t0)). | (3.12) |
This corresponds to a Lie series integrator of order n. Similarly to the Taylor series method, one has to find the analytical expression of the successive action of the operator LHZ onto the vector X(t0). For further details concerning Lie series one can refer to [72,73]. The equivalence between the Lie and Taylor series approaches can be seen in the following way: for each element xi and δxi of the phase space vector X we compute
L0HZxi=Id⋅xi=xi(t0),1≤i≤2N,L1HZxi=2N∑j=1dxjdt∂xi∂xj=dxidt=fi,1≤i≤2N,L2HZxi=2N∑j=1dxjdt∂fi∂xj=dfidt=d2xidt2,1≤i≤2N,⋮L0HZδxi=Id⋅δxi=δxi(t0),1≤i≤2N,L1HZδxi=2N∑j=1dδxjdt∂δxi∂δxj=dδxidt=f2N+i,1≤i≤2N,L2HZδxi=2N∑j=1dxjdt∂f2N+i∂xj+dδxjdt∂f2N+i∂δxj=df2N+idt=d2δxidt2,1≤i≤2N,⋮ |
Therefore L0HZX=X, L1HZX=dXdt, L2HZX=d2Xdt2 etc.
A huge hurdle concerning the applicability of the Taylor and Lie series methods can be the explicit derivation of the differential operators (see [73] and references therein). Over the years other methods have been developed in order to overcome such issues and to efficiently and accurately approximate Eq. (3.8) up to a certain order in n. One way to perform this task was through the use of the well-known 'Runge-Kutta family' of algorithms (see e.g., [40,73]) which nowadays is one of the most popular general purpose schemes for numerical computations. This is a sufficient motivation for us to introduce a s−stage Runge-Kutta method of the form
X(t0+τ)=X(t0)+τs∑i=1biki,withki=f(t0+ciτ,X(t0)+τi−1∑j=1ai,jkj)andci=i−1∑j=1ai,j, | (3.13) |
where the real coefficients ai,j, bi with i,j=1,…,s are appropriately chosen in order to obtain the desired accuracy (see e.g., [66,Sec. Ⅱ.1]). Eq. (3.13) can be understood in the following way: in order to propagate the phase space vector X from time t=t0 to t=t0+τ, we compute s intermediate values X1,X2,…,Xs and s intermediate derivatives k1,k2,…,ks, such that ki=f(Xi), at times ti=t0+τ∑i−1j=1aij. Then each Xi is found through a linear combination of the known ki values, which are added to the initial condition X(t0).
In this work we use a 12-stage explicit Runge-Kutta algorithm of order 8, called DOP853 [66,Sec. Ⅱ.5]. This method is the most precise scheme among the Runge-Kutta algorithms presented in [66] (see Fig. 4.2 of [66]). A free access Fortran77 implementation of DOP853 is available in [74] (see also the appendix of [66]). As in the case of the TIDES package, apart from the integration time step τ, DOP853 also admits a 'one-step' precision δ based on embedded formulas of orders 5 and 3 (see [66,Sec. Ⅱ.10], or [60] and references therein for more details).
Hamiltonian systems are characterized by a symplectic structure (see e.g., [40,Chap. VI] and references therein) and they may also possess integrals of motion, like for example the energy, the angular momentum, etc. SIs are appropriate for the numerical integration of Hamiltonian systems as they keep the computed energy (i.e., the value of the Hamiltonian) of the system almost constant over the integration in time. Let us remark that, in general, SIs do not preserve any additional conserved quantities of the system (for a specific exception see [75]). SIs have already been extensively used in fields such as celestial mechanics [76], molecular dynamics [77], accelerator physics [43], condensed matter physics [22,29,33], etc.
Several methods can be used to build SIs [40]. For instance, it has been proved that under certain circumstances the Runge-Kutta algorithm can preserve the symplectic structure (see e.g., [78,79]). However, the most common way to construct SIs is through the split of the Hamiltonian function into the sum of two or more integrable parts. For example, many Hamiltonian functions are expressed as the sum of the kinetic and potential energy terms, with each one of them corresponding to an integrable system (see Appendices B.1 and B.4). Let us remark that, in general, even if each component of the total Hamiltonian is integrable the corresponding analytical solution might be unknown. Nevertheless, we will not consider such cases in this work. In our study we also want to track the evolution of the deviation vector w(t) by solving Eq. (3.7). Indeed this can be done by using SIs since upon splitting the Hamiltonian into integrable parts we know analytically for each part the exact mapping x(t)→x(t+τ), along with the mapping δx(t)→δx(t+τ).
Let us outline the whole process by considering a general autonomous Hamiltonian function H(q,p), which can be written as sum of I integrable intermediate Hamiltonian functions Ai, i.e., H=∑Ii=1Ai. This decomposition implies that the operator eτLAiZ in the formal solution in Eq. (3.11) of each intermediate Hamiltonian function Ai is known. A symplectic integration scheme to integrate the Hamilton equations of motion from t0 to t0+τ consists in approximating the action of eτLHZ in Eq. (3.11) by a product of the operators eγiτLAiZ for a set of properly chosen coefficients γi. In our analysis we will call as number of steps of a particular SI the total number of successive application of each individual operator eτLAiZ. Further details about this class of integrators can be found in [80,81,82] and references therein. In what follows, we consider the most common cases where the Hamiltonian function can be split into two (I=2) or three (I=3) integrable parts. Let us remark here that, in general, the split H=∑Ii=1Ai is not necessarily unique (see Appendix B.1). Studying the efficiency and the stability of different SIs upon different choices of splitting the Hamiltonian is an interesting topic by itself, which, nevertheless, is beyond the scope of our work.
Let us consider that the Hamiltonian H(q,p) can be separated into two integrable parts, namely H=A(q,p)+B(q,p). Then we can approximate the action of the operator eτLHZ=eτ(LAZ+LBZ) by the successive actions of products of the operators eτLAZ and eτLBZ [80,81,82,83]
eτLHZ=p∏j=1ecjτLAZedjτLBZ+O(τn+1), | (3.14) |
for appropriate choices of the real coefficients cj,dj with j=1,…,p. Different choices of p and coefficients cj, dj lead to schemes of different accuracy. In Eq. (3.14) the integer n is called the order of a symplectic integrator.
The Hamiltonian function of the 1D α-FPUT system [Eq. (2.1)] can be split into two integrable parts H1F=A(p)+B(q), with each part possessing N cyclic coordinates. The kinetic part
A(p)=N∑i=0p2i2, | (3.15) |
depends only on the generalized momenta, whilst the potential part
B(q)=N∑i=012(qi+1−qi)2+α3(qi+1−qi)3, | (3.16) |
depends only on the generalized positions. This type of split is the most commonly used in the literature, therefore a large number of SIs have been developed for such splits. Below we briefly review the SIs used in our analysis, based mainly on results presented in [56].
Symplectic integrators of order two. These integrators constitute the most basic schemes we can develop from Eq. (3.14)
LF: The simplest example of Eq. (3.14) is the so-called Störmer-Verlet or leap-frog scheme (e.g., see [40,Sect. I.3.1] and [83]) having 3 individual steps
LF(τ)=ea1τLAZeb1τLBZea1τLAZ, | (3.17) |
where a1=12 and b1=1.
SABA2/SBAB2: We consider the SABA2 and the SBAB2 SIs with 5 individual steps
SABA2(τ)=ea1τLAZeb1τLBZea2τLAZeb1τLBZea1τLAZ, | (3.18) |
where a1=12−12√3, a2=1√3 and b1=12, and
SBAB2(τ)=eb1τLBZea1τLAZeb2τLBZea1τLAZeb1τLBZ, | (3.19) |
with a1=12, b1=16 and b2=23. These schemes were presented in [84], where they were named the (4, 2) methods, and also used in [55]. We note that the SABA2 and SBAB2 SIs (as well as other two part split SI schemes) have been introduced for Hamiltonian systems of the form H=A+εB, with ε being a small parameter. Both the SABA2 and SBAB2 integrators have only positive time steps and are characterized by an accuracy of order O(τ4ε+τ2ε2) [55]. Although these integrators are particularly efficient for small perturbations (ε≪1), they have also shown a very good performance in cases of ε=1 (see e.g., [24]).
ABA82: In addition, we use in our analysis the SI
ABA82(τ)=ea1τLAZeb1τLBZea2τLAZeb2τLBZea3τLAZeb2τLBZea2τLAZeb1τLBZea1τLAZ, | (3.20) |
with 9 individual steps [84,85], where the constants ai,bi with i=1,2,3 can be found in Table 2 of [85]. We note that the ABA82 method is called SABA4 in [55].
Symplectic integrators of order four. The order of symmetric SIs can be increased by using a composition technique presented in [80]. According to that approach starting from a symmetric SI S2n(τ) of order 2n(n≥1), we can construct a SI S2n+2(τ) of order 2n+2 as *
*We adopt the notation of Iserles and Quispel [86].
S2n+2(τ)=S2n((1+d)τ)S2n(−(1+2d)τ)S2n((1+d)τ),whered=21/(2n+1)−12−21/(2n+1). | (3.21) |
FR4: Using the composition given in Eq. (3.21) for the LF SI of Eq. (3.17) we construct a SI which we name FR4 [80,87] having 7 individual steps
FR4(τ)=ea1τLAZeb1τLBZea2τLAZeb2τLBZea2τLAZeb1τLBZea1τLAZ, | (3.22) |
with coefficients a1=12(2−21/3), a2=1−21/32(2−21/3), b1=12−21/3 and b2=−21/32−21/3.
SABA2Y4/SBAB2Y4: Applying the composition given in Eq. (3.21) to the SABA2 [Eq. (3.18)] and the SBAB2 [Eq. (3.19)] integrators we obtain the fourth order SIs SABA2Y4 and SBAB2Y4 having 13 individual steps. In particular, we get
SABA2Y4(τ)=ed1a1τLAZed1b1τLBZed1a2τLAZed1b1τLBZea0τLAZed0b1τLBZed0a2τLAZed0b1τLBZ×ea0τLAZed1b1τLBZed1a2τLAZed1b1τLBZed1a1τLAZ, | (3.23) |
with coefficients d0=−21/32−21/3, d1=12−21/3, a1=12−12√3, a2=1√3, b1=12, and a0=d1a1+d0a1, and
SBAB2Y4(τ)=ed1b1τLBZed1a1τLAZed1b2τLBZed1a1τLAZeb0τLBZed0a1τLAZed0b2τLBZed0a1τLAZ×eb0τLBZed1a1τLAZed1b2τLBZed1a1τLAZed1b1τLBZ, | (3.24) |
with coefficients d0=−21/32−21/3, d1=12−21/3, a1=12, b1=16, b2=23 and b0=d1b1+d0b1.
ABA82Y4: Using the composition given in Eq. (3.21) for the second order scheme ABA82 of Eq. (3.20) we obtain an integrator with 25 individual steps having the form
ABA82Y4(τ)=ed1a1τLAZed1b1τLBZed1a2τLAZed1b2τLBZed1a3τLAZed1b2τLBZed1a2τLAZed1b1τLBZea0τLAZ×ed0b1τLBZed0a2τLAZed0b2τLBZed0a3τLAZed0b2τLBZed0a2τLAZed0b1τLBZea0τLAZ×ed1b1τLBZed1a2τLAZed1b2τLBZed1a3τLAZed1b2τLBZed1a2τLAZed1b1τLBZed1a1τLAZ, | (3.25) |
where d0=−21/32−21/3, d1=12−21/3, while ai,bi with i=1,2,3 can be found in Table 2 of [85]. Here a0=d1a1+d0a1.
SABA2K/SBAB2K: As was explained in [55] the accuracy of the SABAn (and the SBABn) class of SIs can be improved by a corrector term K={B,{B,A}}, defined by two successive applications of Poisson brackets ({⋅,⋅}), if K corresponds to a solvable Hamiltonian function. In that case, the second order integration schemes can be improved by the addition of two extra operators with negative time steps in the following way
SABAnK(τ)≡e−gτ22LKZSABAn(τ)e−gτ22LKZ, | (3.26) |
with analogous result holding for the SBABn scheme. By following this approach for the SABA2 and SBAB2 SIs [which are of the order O(τ4ε+τ4ε2)] we produce the fourth order SIs SABA2K and SBAB2K, with g=(2−√3)/24 and g=1/72 respectively. These new integration schemes are of the order O(τ4ε+τ4ε2) [55].
ABA864/ABAH864: The fourth order SIs ABA864 and ABAH864 were proposed in [57,85]. They have respectively 15 and 17 individual steps and have the form
ABA864(τ)=ea1τLAZeb1τLBZea2τLAZeb2τLBZea3τLAZeb3τLBZea4τLAZeb4τLBZea4τLAZeb3τLBZ×ea3τLAZeb2τLBZea2τLAZeb1τLBZea1τLAZ, | (3.27) |
with coefficients ai,bii=1,2,3,4 taken from Table 3 of [57], and
ABAH864(τ)=ea1τLAZeb1τLBZea2τLAZeb2τLBZea3τLAZeb3τLBZea4τLAZeb4τLBZea5τLAZeb4τLBZea4τLAZeb3τLBZ×ea3τLAZeb2τLBZea2τLAZeb1τLBZea1τLAZ, | (3.28) |
with coefficients ai,bi,i=1,2,…,5 found in Table 4 of [57]. We note that both schemes were designed for near-integrable systems of the form H=A+εB, with ε being a small parameter, but the construction of ABAH864 was based on the assumption that the integration of the B part cannot be done explicitly, but can be approximated by the action of some second order SI, since B is expressed as the sum of two explicitly integrable parts, i.e., B=B1+B2. The ABA864 and ABAH864 SIs are of order four, but their construction satisfy several other conditions at higher orders, improving in this way their performance.
Symplectic integrators of order six. Applying the composition technique of Eq. (3.21) to the fourth order SIs FR4 [Eq. (3.22)], SABA2Y4 [Eq. (3.23)], SBAB2Y4 [Eq. (3.24)], ABA82Y4 [Eq. (3.25)], SABA2K and ABA864 [Eq. (3.27)], we construct the sixth order SIs FR4Y6, SABA2Y4Y6, SBAB2Y4Y6, ABA82Y4Y6, SABA2KY6 and ABA864Y6 with 19, 37, 37, 73, 19 and 43 individual steps.
In [80] a composition technique using fewer individual steps than the one obtained by the repeated application of Eq. (3.21) to SIs of order two was proposed, having the form
S6(τ)=S2(w3τ)S2(w2τ)S2(w1τ)S2(w0τ)S2(w1τ)S2(w2τ)S2(w3τ), | (3.29) |
whose coefficients wi,i=0,1,2,3 are given in Table 1 in [80] for the case of the so-called 'solution A' of that table. Here S2 and S6 respectively represent a second and a sixth order symmetric SI. Note that Eq. (3.29) corresponds to the composition scheme s6odr6 of [88]. Applying the composition given in Eq. (3.29) to the SABA2 [Eq. (3.18)], the SBAB2 [Eq. (3.19)] and the ABA82 [Eq. (3.20)] SIs we generate the order six schemes SABA2Y6, SBAB2Y6 and ABA82Y6 having 29, 29 and 57 individual steps respectively.
We also consider in our study the composition scheme s9odr6b of [88] which is based on 9 successive applications of S2
s9odr6b(τ)=S2(δ1τ)S2(δ2τ)S2(δ3τ)S2(δ4τ)S2(δ5τ)S2(δ4τ)S2(δ3τ)S2(δ2τ)S2(δ1τ). | (3.30) |
The values of δi,i=1,2,…,5 in Eq. (3.30) can be found in the Appendix of [88]. Furthermore, we also implement the composition method
s11odr6(τ)=S2(γ1τ)S2(γ2τ)S2(γ3τ)S2(γ4τ)S2(γ5τ)S2(γ6τ)S2(γ5τ)S2(γ4τ)S2(γ3τ)S2(γ2τ)S2(γ1τ) | (3.31) |
of [89], which involves 11 applications of a second order SI S2, whose coefficients γi, i=1,2,…,6 are reported in Section 4.2 of [89]. Using the SABA2 of Eq. (3.18) as S2 in Eqs. (3.30) and (3.31), we respectively build two SIs of order six, namely the s9SABA26 and the s11SABA26 SIs with 37 and 45 individual steps. In addition, using the ABA82 integrator of Eq. (3.20) as S2 in Eqs. (3.30) and (3.31) we construct two other order six SIs namely the s9ABA82_6 and s11ABA82_6 schemes with 73 and 89 individual steps respectively.
Runge-Kutta-Nyström methods: In addition, we consider in our analysis two SIs of order six, belonging in the category of the so-called Runge-Kutta-Nyström (RKN) methods (see e.g., [40,90,44,71] and references therein), which respectively have 21 and 29 individual steps
SRKNb11(τ)=eb1τLBZea1τLAZeb2τLBZea2τLAZeb3τLBZea3τLAZeb4τLBZea4τLAZeb5τLBZea5τLAZeb6τLBZea5τLBZ×eb5τLAZea4τLBZeb4τLAZea3τLBZeb3τLAZea2τLAZeb2τLAZea1τLAZeb1τLAZ, | (3.32) |
and
SRKNa14(τ)=ea1τLAZeb1τLAZea2τLAZeb2τLAZ×…×ea7τLAZeb7τLAZea8τLAZeb7τLAZea7τLAZ×…×eb2τLAZea2τLAZeb1τLAZea1τLAZ. | (3.33) |
The values of the coefficients appearing in Eqs. (3.32) and (3.33) can be found in Table 3 of [90]. This class of integrators has, for example, been successfully implemented in a recent investigation of the chaotic behavior of the DNA molecule [91].
Symplectic integrators of order eight. Following [80] we can construct an eighth order SI S8, starting from a second order one S2, by using the composition
S8(τ)=S2(w7τ)S2(w6τ)S2(w5τ)S2(w4τ)S2(w3τ)S2(w2τ)S2(w1τ)S2(w0τ)S2(w1τ)S2(w2τ)×S2(w3τ)S4(w4τ)S2(w5τ)S2(w6τ)S2(w7τ). | (3.34) |
In our study we consider two sets of coefficients wi, i=1,…,7, and in particular the ones corresponding to the so-called 'solution A' and 'solution D' in Table 2 of [80]. Using in Eq. (3.34) the SABA2 [Eq. (3.18)] SI as S2 we construct the eighth order SIs SABA2Y8_A (corresponding to 'solution A') and SABA2Y8_D (corresponding to 'solution D') with 61 individual steps each. In a similar way the use of ABA82 [Eq. (3.20)] in Eq. (3.34) generates the SIs ABA82Y8_A and ABA82Y8_D with 121 individual steps each.
In addition, considering the composition scheme s15odr8 of [88], having 15 applications of S2, we construct the eighth order SI s15SABA28 [s15ABA82_8] having 61 [121] individual steps when SABA2 of Eq. (3.18) [ABA82 of Eq. (3.20)] is used in the place of S2.
Furthermore, implementing the composition technique s19odr8 presented in [89], which uses 19 applications of a second order SI S2, we construct the SI s19SABA28 [s19ABA82_8] with 77 [153] individual steps, when SABA2 of Eq. (3.18) [ABA82 of Eq. (3.20)] is used in the place of S2.
The SIs of this section will be implemented to numerically integrate the α-FPUT model [Eq. (2.1)] since this Hamiltonian system can be split into two integrable parts A(p) and B(q). In Section B.1 of the Appendix the explicit forms of the operators eτLAZ and eτLBZ are given, along with the operator eτLKZ of the corrector term used by the SABA2K and SBAB2K SIs. In addition, in Section B.4 of the Appendix the explicit forms of the tangent map method operators of some commonly used lattice systems, whose Hamiltonians can be split in two integrable parts, are also reported.
Let us now consider the case of a Hamiltonian function H(q,p) which can be separated into three integrable parts, namely H(q,p)=A(q,p)+B(q,p)+C(q,p). This for example could happen because the Hamiltonian function may not be split into two integrable parts, or to simplify the solution of one of the two components, A or B, of the two part split schemes discussed in Section 3.2.1. In such cases we approximate the action of the operator eτLHZ of Eq. (3.11) by the successive application of operators eτLAZ, eτLBZ and eτLCZ i.e.,
eτLHZ=p∏j=1ecjτLAZedjτLBZeejτLCZ+O(τn+1), | (3.35) |
for appropriate choices of the real coefficients cj, dj and ej with j=1,…,p. As in Eq. (3.14), in Eq. (3.35) the integer n is the order of a symplectic integrator.
As examples of Hamiltonians which can be split in three integrable parts we mention the Hamiltonian function of a free rigid body [92] and the Hamiltonian functions of the 1D [Eq. (2.2)] and 2D [Eq. (2.3)] DDNLS models we consider in this work. For example the 1D DDNLS Hamiltonian of Eq. (2.2) can be split in the following three integrable Hamiltonians: a system of N independent oscillators
A1=N∑i=1ϵiJi+β2J2i, | (3.36) |
where Ji=(q2i+p2i)/2, i=1,…,N are N constants of motion, and the Hamiltonian functions of the q− and p−hoppings
B1=−N∑i=1pipi+1,andC1=−N∑i=1qiqi+1, | (3.37) |
with each one of them having N cyclic coordinates. The three part split of the 2D DDNLS of Eq. (2.3) can be found in Section B.3 of the Appendix [Eq. (A.17)].
We note that a rather thorough survey on three part split SIs can be found in [49,50]. We decided to include in our study a smaller number of schemes than the one presented in these works, focusing on the more efficient SIs. We briefly present these integrators below.
Symplectic integrators of order two. We first present the basic three part split scheme obtained by the application of the Störmer-Verlet/leap-frog method to three-part separable Hamiltonians. This scheme has 5 individual steps and we call it ABC2 [93]
ABC2(τ)=ea1τLAZeb1τLBZec1τLCZeb1τLBZea1τLAZ, | (3.38) |
where a1=12, b1=12 and c1=1.
Symplectic integrators of order four. In order to built higher order three part split SIs we apply some composition techniques on the basic ABC2 SI of Eq. (3.38).
ABCY4: Using the composition given in Eq. (3.21) for n=1, we construct
ABCY4(τ)=ABC2(d1τ)ABC2(d0τ)ABC2(d1τ), | (3.39) |
with d0=−21/32−21/3 and d1=12−21/3. This integrator has 13 individual steps, it has been explicitly introduced in [93] and implemented in [50], where it was called ABC4[Y].
ABCS4: Implementing a composition scheme which was introduced in [81] and studied in [88] (where it was named s5odr4) we obtain the SI
ABCS4(τ)=ABC2(p2τ)ABC2(p2τ)ABC2((1−4p2)τ)ABC2(p2τ)ABC2(p2τ), | (3.40) |
where p2=14−41/3 and 1−4p2=−41/34−41/3, having 21 individual steps. This integrator was denoted as ABC4[S] in [50].
SS864S: Using the ABAH864 integrator of Eq. (3.28) where B is considered to be the sum of functions B1 and C1 of Eq. (3.37), i.e. B=B1+C1, and its solution is approximated by the second order SABA2 SI of Eq. (3.18), we construct a SI with 49 steps, which we call SS864S. This integrator has been implemented for the integration of the equations of motion of the 1D DDNLS system [Eq. (2.2)] in [49], where it was called SS4864.
Symplectic integrators of order six.
ABCY4Y6/ABCS4Y6: Applying the composition technique of Eq. (3.21) to the fourth order SIs ABCY4 [Eq. (3.39)] and ABCS4 [Eq. (3.40)], we respectively construct the schemes ABCY4Y6 and ABCS4Y6 with 37 and 49 individual steps.
ABCY6_A: Using the composition given in Eq. (3.29) we build a sixth order SI with 29 individual steps, considering the integrator ABC2 in the place of S2
ABCY6_A(τ)=ABC2(w3τ)ABC2(w2τ)ABC2(w1τ)ABC2(w0τ)ABC2(w1τ)ABC2(w2τ)ABC2(w3τ). | (3.41) |
In particular, we consider in this construction the coefficients wi, i=0,1,2,3, corresponding to the 'solution A' of Table 1 in [80]. Note that this SI has already been implemented in [49,50], where it was denoted as ABC6[Y].
s9ABC6: Implementing the composition given in Eq. (3.30), with ABC2 in the place of S2, we get
s9ABC6(τ)=ABC2(δ1τ)ABC2(δ2τ)ABC2(δ1τ)ABC2(δ4τ)×ABC2(δ5τ)ABC2(δ4τ)ABC2(δ3τ)ABC2(δ2τ)ABC2(δ1τ), | (3.42) |
which was referred to as ABC6[KL] in [49,50].
s11ABC6: From the composition given in Eq. (3.31) we create the SI scheme
s11ABC6(τ)=ABC2(γ1τ)ABC2(γ2τ)×⋯×ABC2(γ5τ)ABC2(γ6τ)ABC2(γ5τ)×⋯×ABC2(γ2τ)ABC2(γ1τ), | (3.43) |
which has 45 individual steps. This integrator was referred to as ABC6[SS] in [49,50].
Symplectic integrators of order eight.
ABCY8_A/ABCY8_D: Based on the composition given in Eq. (3.34) we construct the SI
ABCY8(τ)=ABC2(w7τ)ABC2(w6τ)×⋯×ABC2(w1τ)ABC2(w0τ)ABC2(w1τ)×⋯×ABC2(w6τ)ABC2(w7τ), | (3.44) |
setting ABC2 in the place of S2. This SI has 61 individual steps. Considering the 'solution A' of Table 2 in [80] for the coefficients wi, 0≤i≤7, we obtain the ABCY8_A SI, while the use of 'solution D' of the same table leads to the construction of the SI ABCY8_D.
s17ABC8: Consider the composition method s17odr8b of [88] we build the SI (referred to as ABC8[KL] in [49,50]) s17ABC8 having 69 individual steps
s17ABC8(τ)=ABC2(δ1τ)ABC2(δ2τ)×⋯×ABC2(δ8τ)ABC2(δ9τ)ABC2(δ8τ)×⋯×ABC2(δ2τ)ABC2(δ1τ). | (3.45) |
s19ABC8: Finally, we also implement the composition s19odr8b reported in [89,Eq. (13)] and construct the SI
s19ABC8(τ)=ABC2(γ1τ)ABC2(γ2τ)×⋯×ABC2(γ8τ)ABC2(γ9τ)ABC2(γ8τ)×⋯×ABC2(γ2τ)ABC2(γ1τ), | (3.46) |
which has 72 individual steps. We note that this scheme corresponds to the SI ABC8[SS] considered in [49,50].
The explicit forms of the operators related to the three part split of the DDNLS Hamiltonians [Eqs. (2.2) and (2.3)] are given in Sections. B.2 and B.3 of the Appendix.
We test the efficiency of the integrators presented in Section 3 by using them to follow the dynamical evolution of the α-FPUT model [Eq. (2.1)], the 1D DDNLS system [Eq. (2.2)] and the 2D DDNLS Hamiltonian [Eq. (2.3)]. For each model, given an initial condition X(t0)=(x(t0),δx(t0)) we compute the trajectory {x(tn)}n∈N with x(t)=(q1(t),q2(t),…,qN(t),p1(t),p2(t),…,pN(t)) and check the integrators' efficiency through their ability to correctly reproduce certain observables of the dynamics. We also follow the evolution of a small initial perturbation to that trajectory w(t0)=δx(t0)=(δq1(t0),δq2(t0),…,δqN(t0),δp1(t0),δp2(t0),…,δpN(t0)) and use it to compute the time evolution of the finite time mLE [35,36,37]
X1(t)=1tln[‖w(t0+τ)‖‖w(t0)‖], | (4.1) |
in order to characterize the regular or chaotic nature of the trajectory through the estimation of the most commonly used chaos indicator, the mLE χ, which is defined as χ=limt→+∞X1(t). In Eq. (4.1) ‖⋅‖ is the usual Euclidian norm, while w(t0) and w(t0+τ) are respectively the deviation vectors at t=t0 and t0+τ>t0. In the case of regular trajectories X1(t) tends to zero following the power law [37,94]
X1(t)∝t−1, | (4.2) |
whilst it takes positive values for chaotic ones.
We present here results on the computational efficiency of the symplectic and non-symplectic schemes of Section 3 for the case of the α-FPUT chain [Eq. (2.1)]. As this system can be split into two integrable parts we will use for its study the two part split SIs of Section 3.2.1. In our investigation we consider a lattice of N=210 sites with α=0.25 and integrate up to the final time tf=106 two sets of initial conditions:
● Case IF: We excite all lattice sites by attributing to their position and momentum coordinates a randomly chosen value from a uniform distribution in the interval [−1,1]. These values are rescaled to achieve a particular energy density, namely H1F/N=0.1.
● Case IIF: Same as in case IF, but for H1F/N=0.05.
We consider these two initial conditions in an attempt to investigate the potential dependence of the performance of the tested integrators on initial conditions [73,Sec. 8.3]. Since we have chosen non-localized initial conditions, we also use an initial normalized deviation vector w(t) whose components are randomly selected from a uniform distribution in the interval [−1,1].
To evaluate the performance of each integrator we investigate how accurately it follows the considered trajectories by checking the numerical constancy of the energy integral of motion, i.e., the value of H1F [Eq. (2.1)]. This is done by registering the time evolution of the relative energy error
Er(t)=|H1F(t)−H1F(0)H1F(0)|, | (4.3) |
at each time step. In our analysis we consider two energy error thresholds Er≈10−5 and Er≈10−9. The former, Er≈10−5, is typically considered to be a good accuracy in many studies in the field of lattice dynamics, like for example in investigations of the DKG and the DDNLS models, as well as in systems of coupled rotors (see for example [24,25,26,28,29,30]). In some cases, e.g., for very small values of conserved quantities, one may desire more accurate computations. Then Er≈10−9 is a more appropriate accuracy level. In addition, in order to check whether the variational equations are properly evolved, we compute the finite time mLE X1(t) [Eq. (4.1)].
In Figure 1 we show the time evolution of the relative energy error Er(t) [panels (a) and (d)], the finite time mLE X1(t) [panels (b) and (e)], and the required CPU time TC [panels (c) and (f)], for cases IF and IIF respectively, when the following four integrators were used: the fourth order SI ABA864 (blue curves), the sixth order SI SABA2Y6 (red curves), the DOP853 scheme (green curves) and the TIDES package (brown curves). These results are indicative of our analysis as in our study we considered in total 37 different integrators (see Tables 1 and 2). In Figure 1 the integration time steps τ of the SIs (reported in Tables 1 and 2) were appropriately chosen in order to achieve Er≈10−9, while for the DOP853 algorithm and the TIDES package Er(t) eventually grows in time as a power law [Figure 1(a),(d)]. Nevertheless, all schemes succeed in capturing correctly the chaotic nature of the dynamics as they do not present any noticeable difference in the computation of the finite time mLE X1 in Figure 1(b), (e). For both sets of initial conditions X1 eventually saturates to a constant positive value indicating that both trajectories are chaotic. The CPU time TC needed for the integration of the equations of motion and the variational equations are reported in Figure 1(c), (f). From these plots we see that the SIs need less computational time to perform the simulations than the DOP853 and TIDES schemes.
In Table 1 (Table 2) we present information on the performance of all considered integration schemes for the initial condition of case IF (case IIF). From the results of these tables we see that the performance and ranking (according to TC) of the integrators do not practically depend on the considered initial condition. It is worth noting that although the non-symplectic schemes manage to achieve better accuracies than the symplectic ones, as their Er values are smaller [Figure 1(a),(d)], their implementation is not recommended for the long time evolution of the Hamiltonian system, because they require more CPU time and eventually their Er values will increase above the bounded Er values obtained by the symplectic schemes.
From the results of Tables 1 and 2 we see that the best performing integrators are the fourth order SIs ABA864 and ABAH864 for Er≈10−5, and the sixth order SIs SRKNa14 and SRKNb11 for Er≈10−9. We note that the best SI for Er≈10−5, the ABA864 scheme, shows a quite good behavior also for Er≈10−9, making this integrator a valuable numerical tool for dynamical studies of multidimensional Hamiltonian systems. We remark that the eighth order SIs we implemented to achieve the moderate accuracy level Er≈10−5 exhibited an unstable behavior failing to keep their Er values bounded. A similar behavior was also observed for the two RKN schemes SRKNa14, SRKNb11 when they were used to obtain Er≈10−5. Thus, the higher order SIs are best suited for more accurate computations. It is also worth mentioning here that the ranking presented in Tables 1 and 2 is indicative of the performance of the various SIs in the sense that small changes in the implementation (e.g., a change in the last digit of the used τ value) of integrators with similar behaviors (i.e., similar TC values) could interchange their ranking positions without any noticeable difference in the produced results.
Er≈10−5 | Er≈10−9 | ||||||||
Integrator | n | Steps | τ | TC | Integrator | n | Steps | τ | TC |
ABA864 | 4 | 15 | 0.6 | 88 | SRKNa14 | 6 | 29 | 0.45 | 160 |
ABAH864 | 4 | 17 | 0.55 | 115 | SRKNb11 | 6 | 23 | 0.35 | 177 |
SABA2Y6 | 6 | 29 | 0.575 | 167 | s11SABA26 | 6 | 45 | 0.3 | 536 |
ABA864Y6 | 6 | 43 | 0.625 | 202 | s19SABA28 | 8 | 77 | 0.45 | 594 |
s9SABA26 | 6 | 37 | 0.575 | 205 | s9ABA82_6 | 6 | 73 | 0.35 | 607 |
FR4 | 4 | 7 | 0.14 | 228 | s15SABA28 | 8 | 61 | 0.35 | 611 |
SBAB2Y6 | 6 | 29 | 0.5 | 233 | SABA2Y6 | 6 | 29 | 0.14 | 683 |
SABA2K | 4 | 9 | 0.3 | 234 | ABA864 | 4 | 15 | 0.08 | 717 |
ABA82Y4 | 4 | 25 | 0.375 | 240 | s19ABA82_8 | 8 | 153 | 0.65 | 773 |
s11SABA26 | 6 | 45 | 0.65 | 247 | s9SABA26 | 6 | 37 | 0.16 | 779 |
SABA2Y4 | 4 | 13 | 0.18 | 265 | ABA864Y6 | 6 | 43 | 0.16 | 791 |
ABA82 | 2 | 5 | 0.125 | 278 | s15ABA82_8 | 8 | 121 | 0.475 | 838 |
ABA82Y6 | 6 | 57 | 0.675 | 283 | ABA82Y6 | 6 | 57 | 0.2 | 841 |
s15SABA28 | 8 | 61 | 0.65 | 339 | s11ABA82_6 | 6 | 89 | 0.275 | 941 |
SABA2 | 2 | 5 | 0.07 | 347 | SBAB2Y6 | 6 | 29 | 0.12 | 965 |
s19SABA28 | 8 | 77 | 0.775 | 356 | ABAH864 | 4 | 17 | 0.055 | 1013 |
SBAB2Y4 | 4 | 13 | 0.18 | 358 | SABA2Y8_D | 8 | 61 | 0.175 | 1223 |
FR4Y6 | 6 | 19 | 0.21 | 366 | ABA82Y8_D | 8 | 121 | 0.25 | 1575 |
s9ABA82_6 | 6 | 73 | 0.575 | 369 | SABA2Y4Y6 | 6 | 37 | 0.07 | 1701 |
s11ABA82_6 | 6 | 89 | 0.675 | 382 | FR4Y6 | 6 | 19 | 0.45 | 1787 |
SBAB2 | 2 | 5 | 0.07 | 387 | ABA82Y4Y6 | 6 | 73 | 0.125 | 1932 |
SABA2Y4Y6 | 6 | 37 | 0.3 | 394 | ABA82Y4 | 4 | 25 | 0.0375 | 2156 |
ABA82Y4Y6 | 6 | 73 | 0.525 | 405 | SBAB2Y4Y6 | 6 | 37 | 0.065 | 2239 |
SABA2Y8_D | 8 | 61 | 0.525 | 408 | SABA2K | 4 | 9 | 0.03 | 2344 |
SBAB2K | 4 | 9 | 0.2 | 416 | SABA2KY6 | 6 | 19 | 0.09 | 2465 |
s19ABA82_8 | 8 | 153 | 1.15 | 439 | FR4 | 4 | 7 | 0.01 | 2597 |
SABA2KY6 | 6 | 19 | 0.4 | 535 | SABA2Y4 | 4 | 13 | 0.018 | 2654 |
SBAB2Y4Y6 | 6 | 37 | 0.275 | 553 | SBAB2Y4 | 4 | 13 | 0.018 | 3156 |
s15ABA82_8 | 8 | 121 | 0.775 | 618 | SABA2Y8_A | 8 | 61 | 0.06 | 3570 |
ABA82Y8_D | 8 | 121 | 0.6 | 656 | SBAB2K | 4 | 9 | 0.02 | 4167 |
SABA2Y8_A | 8 | 61 | 0.225 | 1090 | ABA82Y8_A | 8 | 121 | 0.07 | 5624 |
LF | 2 | 3 | 0.018 | 1198 | ABA82 | 2 | 5 | 0.00125 | 27796 |
ABA82Y8_A | 8 | 121 | 0.225 | 1749 | DOP853 | 8 | δ=10−16 | 0.05 | 31409 |
SABA2 | 2 | 5 | 0.0007 | 34595 | |||||
SBAB2 | 2 | 5 | 0.0007 | 39004 | |||||
LF | 2 | 3 | 0.0002 | 95096 | |||||
TIDES | - | δ=10−16 | 0.05 | 232785 |
Er≈10−5 | Er≈10−9 | ||||||||
Integrator | n | Steps | τ | TC | Integrator | n | Steps | τ | TC |
ABA864 | 4 | 15 | 0.6 | 77 | SRKNa14 | 6 | 29 | 0.475 | 156 |
ABAH864 | 4 | 17 | 0.55 | 101 | SRKNb11 | 6 | 23 | 0.35 | 179 |
SABA2Y6 | 6 | 29 | 0.575 | 183 | s11SABA26 | 6 | 45 | 0.3 | 480 |
ABA864Y6 | 6 | 43 | 0.625 | 195 | s19SABA28 | 8 | 77 | 0.45 | 596 |
s9SABA26 | 6 | 37 | 0.575 | 214 | s9ABA82_6 | 6 | 73 | 0.35 | 611 |
ABA82Y4 | 4 | 25 | 0.375 | 224 | ABA864 | 4 | 15 | 0.08 | 613 |
s11SABA26 | 6 | 45 | 0.65 | 227 | s15SABA28 | 8 | 61 | 0.35 | 618 |
FR4 | 4 | 7 | 0.14 | 231 | SABA2Y6 | 6 | 29 | 0.14 | 676 |
SABA2K | 4 | 9 | 0.3 | 241 | s19ABA82_8 | 8 | 153 | 0.65 | 760 |
SBAB2Y6 | 6 | 29 | 0.5 | 255 | ABA864Y6 | 6 | 43 | 0.16 | 811 |
SABA2Y4 | 4 | 13 | 0.18 | 270 | s15ABA82_8 | 8 | 121 | 0.475 | 828 |
ABA82 | 2 | 5 | 0.125 | 280 | s9SABA26 | 6 | 37 | 0.16 | 838 |
ABA82Y6 | 6 | 57 | 0.675 | 285 | ABA82Y6 | 6 | 57 | 0.2 | 937 |
SBAB2Y4 | 4 | 13 | 0.18 | 316 | SBAB2Y6 | 6 | 29 | 0.12 | 964 |
s15SABA28 | 8 | 61 | 0.65 | 329 | ABAH864 | 4 | 17 | 0.055 | 1023 |
s19SABA28 | 8 | 77 | 0.775 | 336 | s11ABA82_6 | 6 | 89 | 0.275 | 1062 |
FR4Y6 | 6 | 19 | 0.21 | 337 | SABA2Y8_D | 8 | 61 | 0.175 | 1230 |
SBAB2K | 4 | 9 | 0.2 | 366 | FR4Y6 | 6 | 19 | 0.45 | 1577 |
s9ABA82_6 | 6 | 73 | 0.575 | 373 | ABA82Y8_D | 8 | 121 | 0.25 | 1613 |
SABA2 | 2 | 5 | 0.07 | 392 | ABA82Y4Y6 | 6 | 73 | 0.125 | 1702 |
SBAB2 | 2 | 5 | 0.07 | 393 | ABA82Y4 | 4 | 25 | 0.0375 | 2113 |
ABA82Y4Y6 | 6 | 73 | 0.525 | 398 | SABA2Y4Y6 | 6 | 37 | 0.07 | 2159 |
SABA2Y8_D | 8 | 61 | 0.525 | 407 | SBAB2Y4Y6 | 6 | 37 | 0.065 | 2310 |
s11ABA82_6 | 6 | 89 | 0.675 | 415 | SABA2KY6 | 6 | 19 | 0.09 | 2380 |
s19ABA82_8 | 8 | 153 | 1.15 | 431 | FR4 | 4 | 7 | 0.01 | 2615 |
SABA2Y4Y6 | 6 | 37 | 0.3 | 444 | SABA2K | 4 | 9 | 0.03 | 2651 |
SBAB2Y4Y6 | 6 | 37 | 0.275 | 533 | SABA2Y4 | 4 | 13 | 0.018 | 3016 |
SABA2KY6 | 6 | 19 | 0.4 | 540 | SBAB2Y4 | 4 | 13 | 0.018 | 3585 |
s15ABA82_8 | 8 | 121 | 0.775 | 598 | SABA2Y8_A | 8 | 61 | 0.06 | 3647 |
ABA82Y8_D | 8 | 121 | 0.6 | 629 | SBAB2K | 4 | 9 | 0.02 | 3663 |
SABA2Y8_A | 8 | 61 | 0.225 | 952 | ABA82Y8_A | 8 | 121 | 0.07 | 5691 |
LF | 2 | 3 | 0.018 | 1059 | ABA82 | 2 | 5 | 0.00125 | 31576 |
ABA82Y8_A | 8 | 121 | 0.225 | 1787 | DOP853 | 8 | δ=10−16 | 0.05 | 31709 |
SABA2 | 2 | 5 | 0.0007 | 39333 | |||||
SBAB2 | 2 | 5 | 0.0007 | 45690 | |||||
LF | 2 | 3 | 0.0002 | 106653 | |||||
TIDES | - | δ=10−16 | 0.05 | 225565 |
We investigate the performance of various integrators of Section 3 for the 1D DDNLS system [{Eq. (2.2)] by considering a lattice of N=210 sites and integrating two sets of initial conditions (for the same reason we did that for the α-FPUT system) up to the final time tf=106. We note that, as was already mentioned in Section 3.2.2, this model can be split into three integrable parts, so we will implement the SIs presented in that section. In particular, we consider the following two cases of initial conditions:
● Case I1D: We initially excite 21 central sites by attributing to each one of them the same constant norm sj=(q2j+p2j)/2=1, 1≤i≤N, for W=3.5 and β=0.62. This choice sets the total norm S1D=21. The random disorder parameters ϵi, 1≤i≤N, are chosen so that the total energy is H1D≈0.0212.
● Case II1D: Similar set of initial conditions as in case I1D but for W=3, β=0.03. The random disorder parameters ϵi, 1≤i≤N, are chosen such that H1D≈3.4444.
We note that cases I1D and II1D have been studied in [33] and respectively correspond to the so-called 'strong chaos' and 'weak chaos' dynamical regimes of this model. As initial normalized deviation vector we use a vector having non-zero coordinates only at the central site of the lattice, while its remaining elements are set to zero.
To evaluate the performance of each implemented integrator we check if the obtained trajectory correctly captures the statistical behavior of the normalized norm density distribution ζj=sj/S1D, 1≤j≤N, by computing the distribution's second moment
m2=N∑j=1(j−¯j)2ζj, | (4.4) |
where ¯j=∑Nj=1jζj is the position of the center of the distribution [22,24,26,29,33,49,50]. We also check how accurately the values of the system's two conserved quantities, i.e., its total energy H1D [Eq. (2.2)] and norm S1D [Eq. (2.4)], are kept constant throughout the integration by evaluating the relative energy error Er(t) [similarly to Eq. (4.3)] and the relative norm error
Sr(t)=|S1D(t)−S1D(0)S1D(0)|. | (4.5) |
In addition, we compute the finite time mLE X1(t) [Eq. (4.1)] in order to characterize the system's chaoticity and check the proper integration of the variational equations.
We consider several three part split SIs, which we divide into two groups: (i) those integrators with order n≤6, which we implement in order to achieve an accuracy of Er≈10−5, and (ii) SIs of order eight used for Er≈10−9. In addition, the two best performing integrators of the first group are also included in the second group. We do not use higher order SIs for obtaining the accuracy level of Er≈10−5 because, as was shown in [49] and also discussed in Section 4.1, usually this task requires large integration time steps, which typically make the integrators unstable. Moreover, increasing the order n of SIs beyond eight does not improve significantly the performance of the symplectic schemes for high precision (Er≈10−9) simulations [49]. Therefore we do not consider such integrators in our study.
In Figure 2 we show the time evolution of the relative energy error Er(t) [panel (a)], the relative norm error Sr(t) [panel (b)], the second moment m2(t) [panel (d)], as well as the norm density distribution ζj at time tf≈106 [panel (c)] for case I1D (we note that analogous results were also obtained for case II1D, although we do not report them here). These results are obtained by the implementation of the the second order SI ABC2 (red curves), the fourth order SI ABCY4 (blue curves), and the non-symplectic schemes DOP853 (green curves) and TIDES (brown curves). The results of Figure 2 are indicative of the results obtained by the integrators listed in Table 3. The integration time step τ of the SIs was selected so that the relative energy error is kept at Er≈10−5 [Figure 2(a)]. From the results of Figure 2(b) we see that the SIs do not keep Sr constant. Nevertheless, our results show that we lose no more than two orders of precision (in the worst case of the ABC2 scheme) during the whole integration. On the other hand, both the relative energy [Er(t)] and norm [Sr(t)] errors of the TIDES and DOP853 integrators increase in time, with the TIDES scheme behaving better than the DOP853 one. Figure 2(c), (d) show that all integrators correctly reproduce the dynamics of the system, as all of them practically produce the same norm density distribution at tf=106 [Figure 2(c)] and the same evolution of the m2(t) [Figure 2(d)]. We note that m2 increases by following a power law m2∝tα with α=1/2, as was expected for the strong chaos dynamical regime (see for example [33] and references therein).
In Figure 3 we show the evolution of the finite time mLE X1 [panels (a) and (c)] and the required CPU time TC [panels (b) and (d)] for the integration of the Hamilton equations of motion and the variational equations for cases I1D [panels (a) and (b)] and II1D [panels (c) and (d)] obtained by using the same integrators of Figure 2. Again the results obtained by these integrators are practically the same for both sets of initial conditions, reproducing the tendency of the finite time mLE to asymptotically decrease according to the power law X1(t)∝tαL with αL≈−0.3 (case I1D) and αL≈−0.25 (case II1D), in accordance to the results of [32,33].
We now check the efficiency of the used symplectic and non-symplectic methods by comparing the CPU time TC they require to carry out the simulations. These results are reported in Table 3 for case I1D and in Table 4 for case II1D. These tables show that the comparative performance of the integrators does not depend on the chosen initial condition, as the ranking of the schemes is practically the same in both tables. As in the case of the α-FPUT model, the DOP853 and TIDES integrators required, in general, more CPU time than the SIs, although they produced more accurate results (smaller Er and Sr values) at least up to tf=106, with TIDES being more precise. The integrators exhibiting the best performance for Er≈10−5 are the sixth order SIs s11ABC6 and s9ABC6, while for Er≈10−9 we have the eighth order SIs s19ABC8 and s17ABC8, with the s11ABC6 scheme performing quite well also in this case.
Er≈10−5 | Er≈10−9 | ||||||||
Integrator | n | Steps | τ | TC | Integrator | n | Steps | τ | TC |
s11ABC6 | 6 | 45 | 0.115 | 3395 | s19ABC8 | 8 | 77 | 0.09 | 7242 |
s9ABC6 | 6 | 37 | 0.095 | 3425 | s17ABC8 | 8 | 69 | 0.08 | 7301 |
ABCY6_A | 6 | 29 | 0.07 | 3720 | s11ABC6 | 6 | 45 | 0.025 | 15692 |
SS864S | 4 | 17 | 0.05 | 6432 | s9ABC6 | 6 | 37 | 0.02 | 16098 |
ABCY4 | 4 | 13 | 0.0125 | 10317 | DOP853 | 8 | δ=10−16 | 0.05 | 18408 |
ABCS4Y6 | 6 | 49 | 0.015 | 35417 | ABCY8_D | 8 | 61 | 0.002 | 258891 |
ABCY4Y6 | 6 | 37 | 0.008 | 40109 | TIDES | − | δ=10−16 | 0.05 | 419958 |
ABCS4 | 4 | 21 | 0.00085 | 267911 | |||||
ABC2 | 2 | 5 | 0.0002 | 320581 |
Er≈10−5 | Er≈10−9 | ||||||||
Integrator | n | Steps | τ | TC | Integrator | n | Steps | τ | TC |
s11ABC6 | 6 | 45 | 0.4 | 1132 | s19ABC8 | 8 | 77 | 0.3 | 2184 |
s9ABC6 | 6 | 37 | 0.285 | 1147 | s17ABC8 | 8 | 69 | 0.225 | 2632 |
ABCY6_A | 6 | 29 | 0.2 | 1308 | s11ABC6 | 6 | 45 | 0.1 | 4137 |
SS864S | 4 | 17 | 0.265 | 1365 | s9ABC6 | 6 | 37 | 0.075 | 4462 |
ABCY4 | 4 | 13 | 0.055 | 2354 | ABCY8_D | 8 | 61 | 0.065 | 8528 |
ABCS4Y6 | 6 | 49 | 0.105 | 4965 | DOP853 | 8 | δ=10−16 | 0.05 | 14998 |
ABCY4Y6 | 6 | 37 | 0.04 | 8091 | TIDES | − | δ=10−16 | 0.05 | 420050 |
ABCS4 | 4 | 21 | 0.02 | 9774 | |||||
ABC2 | 2 | 5 | 0.0055 | 11700 |
We now investigate the performance of the integrators used in Section 4.2 for the computationally much more difficult case of the 2D DDNLS lattice of Eq. (2.3), as its Hamiltonian function can also be split into three integrable parts. In order to test the performance of the various schemes we consider a lattice with N×M=200×200 sites, resulting to a system of 4×40000=160000 ODEs (equations of motion and variational equations). The numerical integration of this huge number of ODEs is a very demanding computational task. For this reason we integrate this model only up to a final time tf=105, instead of the tf=106 used for the α-FPUT and the 1D DDNLS systems. It is worth noting that due to the computational difficulty of the problem very few numerical results for the 2D DDNLS system exist in the literature (e.g., [39,65]). We consider again two sets of initial conditions:
● Case I2D: We initially excite 7×7 central sites attributing to each one of them the same norm si,j=(q2i,j+p2i,j)/2=1/6 so that the total norm is S2D=49/6, for W=15 and β=6. The disorder parameters ϵi,j, 1≤i≤N, 1≤j≤M, are chosen so that the initial total energy is H2D≈1.96.
● Case II2D: We initially excite a single central site of the lattice with a total norm S2D=1, i.e., ζ100,100=1, for W=16, β=1.25 and H2D=0.625.
The initial normalized deviation vector considered in our simulations has random non-zero values only at the 7×7 initially excited sites for case I2D, and only at site i=100, j=100 for case II2D. In both cases, all others elements of the vectors are initially set to zero. Both considered cases belong to a Gibbsian regime where the thermalization processes are well defined by Gibbs ensembles [52,97]. Therefore, we expect a subdiffusive spreading of the initial excitations to take place for both cases, although their initial conditions are significantly different.
As was done in the case of the 1D DDNLS system (Section 4.2), in order to evaluate the performance of the used integrators we follow the time evolution of the normalized norm density distribution ζi,j=si,j/S2D, 1≤i≤N, 1≤j≤M and compute the related second moment m2 and participation number P
m2=N∑i=1M∑j=1‖(i,j)T−(¯i,¯j)T‖2ζi,j,P=1∑Ni=1∑Mj=1ζ2i,j, | (4.6) |
where (¯i,¯j)T=∑i,j(i,j)Tζi,j is the mean position of the norm density distribution. We also evaluate the relative energy [Er(t)] and norm [Sr(t)] errors and compute the finite time mLE X1(t).
In Figure 4 we present results obtained for case I2D by the four best performing SIs (see Table 5), the s11ABC6 (red curves), s9ABC6 (blue curves), ABCY6_A (green curves) and ABCY4 (brown curves) schemes, along with the DOP853 integrator (grey curves). The integration time steps τ of the SIs were adjusted in order to obtain an accuracy of Er≈10−5 [Figure 4(a)], while results for the conservation of the second integral of motion, i.e., the system's total norm, are shown in Figure 4(b). We see that for all SIs the Sr values increase slowly, remaining always below Sr≈10−4, which indicates a good conservation of the system's norm. As in the case of the 1D DDNLS system, the Er and Sr values obtained by the DOP853 integrator increase, although the choice of δ=10−16 again ensures high precision computations.
Er≈10−5 | Er≈10−9 | ||||||||
Integrator | n | Steps | τ | TC | Integrator | n | Steps | τ | TC |
s9ABC6 | 6 | 45 | 0.105 | 13914 | s17ABC8 | 8 | 77 | 0.075 | 36528 |
s11ABC6 | 6 | 37 | 0.125 | 14000 | s19ABC8 | 8 | 69 | 0.08 | 38270 |
ABCY6_A | 6 | 29 | 0.08 | 15344 | s9ABC6 | 6 | 45 | 0.0235 | 65287 |
ABCY4 | 4 | 13 | 0.025 | 23030 | s11ABC6 | 6 | 37 | 0.0275 | 67314 |
SS864S | 4 | 17 | 0.085 | 23887 | ABCY8_D | 8 | 61 | 0.008 | 140506 |
ABCS4Y6 | 6 | 49 | 0.03 | 77424 | DOP853 | 8 | δ=10−16 | 0.05 | 218704 |
ABCY4Y6 | 6 | 37 | 0.0165 | 87902 | |||||
ABCS4 | 4 | 21 | 0.0065 | 132713 | |||||
ABC2 | 2 | 5 | 0.005 | 157694 |
The norm density distributions at the final integration time tf=105 along the axis i=100 [Figure 4(c)] and j=100 [Figure 4(d)] obtained by the various integrators practically overlap indicating the ability of all numerical schemes to correctly capture the system's dynamics, as well as the fact that the initial excitations expand along all directions of the 2D lattice. From Figure 4(e) [Figure 4(f)] we see that m2(t) [P(t)] is increasing according to the power law m2=t1/3 [P=t1/3] as expected from the analysis presented in [30], indicating that the 2D lattice is being thermalized. The results of Figure 4(e), (f) provide additional numerical evidences that all numerical methods reproduce correctly the dynamics. This is also seen by the similar behavior of the finite time mLE curves in Figure 4(g). From the results of this figure we see that X1 exhibits a tendency to decrease following a completely different decay from the X1∝t−1 power law observed for regular motion. This behavior was also observed for the 2D DKG model [56], as well as for the 1D DKG and DDNLS systems in [32,33] where a power law X1(t)∝tαL with αL≈−0.25 and αL≈−0.3 for, respectively, the weak and strong chaos dynamical regimes was established. Further investigations of the behavior of the finite mLE in 2D disordered systems are required in order to determine a potentially global behavior of X1, since here and in [56] only some isolated cases were discussed. Such studies will require the statistical analysis of results obtained for many different disorder realizations, parameter sets and initial conditions. Thus, the utilization of efficient and accurate numerical integrators, like the ones presented in this study, will be of utmost importance for the realization of this goal.
From Tables 5 and 6, where the CPU times TC required by the tested integrators are reported, we see that, as in the case of the 1D DDNLS model, the SIs s11ABC6 and s9ABC6 have the best performance for Er≈10−5 and the SIs s19ABC8 and s17ABC8 for Er≈10−9.
Er≈10−5 | Er≈10−9 | ||||||||
Integrator | n | Steps | τ | TC | Integrator | n | Steps | τ | TC |
s11ABC6 | 6 | 45 | 0.1515 | 11443 | s19ABC8 | 8 | 77 | 0.135 | 20952 |
s9ABC6 | 6 | 37 | 0.11 | 13408 | s17ABC8 | 8 | 69 | 0.0875 | 28966 |
ABCY6_A | 6 | 29 | 0.0775 | 14607 | s11ABC6 | 6 | 45 | 0.0335 | 50301 |
SS864S | 4 | 17 | 0.0915 | 15564 | s9ABC6 | 6 | 37 | 0.024 | 58187 |
ABCY4 | 4 | 13 | 0.0215 | 25898 | ABCY8_D | 8 | 61 | 0.009 | 150045 |
ABCS4Y6 | 6 | 49 | 0.035 | 64423 | DOP853 | 8 | δ=10−16 | 0.05 | 166145 |
ABCY4Y6 | 6 | 37 | 0.01375 | 102580 | |||||
ABCS4 | 4 | 21 | 0.005 | 185615 | |||||
ABC2 | 2 | 5 | 0.00155 | 198534 |
In this work we carried out a methodical and detailed analysis of the performance of several symplectic and non-symplectic integrators, which were used to integrate the equations of motion and the variational equations of some important many-body classical Hamiltonian systems in one and two spatial dimensions: the α-FPUT chain, as well as the 1D and 2D DDNLS models. In the case of the α-FPUT system we used two part split SIs, while for the integration of the DDNLS models we implemented several three part split SIs. In order to evaluate the efficiency of all these integrators we evolved in time different sets of initial conditions and evaluated quantities related to (a) the dynamical evolution of the studied systems (e.g., the second moment of norm density distributions for the DDNLS models), (b) the quantification of the systems' chaotic behavior (i.e., the finite time mLE), and (c) the accurate computation of the systems' integrals of motion (relative energy and norm errors), along with the CPU times needed to perform the simulations.
For the α-FPUT system several two part split SIs showed very good performances, among which we mention the ABA864 and ABAH864 SIs of order four to be the best schemes for moderate energy accuracies (Er≈10−5), while the SRKNa14 and SRKNb11 SIs of order six were the best integration schemes for higher accuracies (Er≈10−9). In particular, the ABA864 scheme appears to be an efficient, general choice as it showed a quite good behavior also for Er≈10−9. Concerning the 1D and the 2D DDNLS models our simulations showed that the SIs s9ABC6 and s11ABC6 (order six), along with the SIs s17ABC8 and s19ABC8 (order eight) are the best integrators for moderate (Er≈10−5) and high (Er≈10−9) accuracy levels respectively.
The DOP853 and TIDES non-symplectic integrators required, in general, much longer CPU times to carry out the simulations, although they produced more accurate results (i.e., smaller Er and Sr values) than the symplectic schemes. Apart from the drawback of the high CPU times, the fact that Er (and Sr) values exhibit a constant increase in time signifies that such schemes should not be preferred over SIs when very long time simulations are needed.
It is worth noting that two part split SIs of order six and higher often do now not produce reliable results for relative low energy accuracies like Er≈10−5 for the α-FPUT system (similar behaviors were reported in [56] for the DKG model). This happens because the required integration time step τ needed to keep the relative energy error at Er≈10−5 is typically large, resulting to an unstable behavior of the integrator i.e., the produced Er values do not remain bounded. Thus, SIs of order n≥6 are more suitable for calculations that require higher accuracies (e.g., Er≈10−9 or lower).
We note that we presented here a detailed comparison of the performance of several two and three part split SIs for the integration of the variational equations through the tangent map method and consequently for the computation of a chaos indicator (the mLE), generalizing, and completing in some sense, some sporadic previous investigation of the subject [56,58,59,60], which were only focused on two part split SIs.
We hope that the clear description of the construction of several two and three part split SIs in Section 3, along with the explicit presentation in the Appendix of the related differential operators for many commonly used classical, many-body Hamiltonians will be useful to researchers working on lattice dynamics. The numerical techniques presented here can be used for the computation of several chaos indicators, apart from the mLE (e.g., the SALI and the GALI methods [98]) and for the dynamical study of various lattice models, like for example of arrays of Josephson junctions in regimes of weak non-integrability [53], granular chains [99] and DNA models [91], to name a few.
We present here the exact expressions of the operators needed by the various SIs we implemented in our study to simultaneously solve the Hamilton equations of motion and the variational equations, or in order words to solve the system of Eq. (3.7)
˙X=(˙x(t),˙δx(t))=f(X)=[J2N⋅DH(x(t))[J2N⋅D2H(x(t))]⋅δx(t)]. | (A.1) |
The Hamiltonian of the α-FPUT chain [Eq. (2.1)] can be split into two integrable parts as
A(p)=N∑i=0 p2i2,B(q)=N∑i=0 12(qi+1−qi)2+α3(qi+1−qi)3. | (A.2) |
As we have already stated, the split into two integrable parts is not necessarily unique. In this particular case another possible choice of integrable splits for the α-FPUT chain is to group together the quadratic terms of the Hamiltonian [i.e., A(p,q)=∑Ni=0 p2i2+12(qi+1−qi)2] and keep separately the nonlinear terms [i.e., B(q)=∑Ni=0 α3(qi+1−qi)3]. The set of equations of motion and variational equations for the Hamiltonian function A(p) is
dXdt=LAZX:{˙qi=pi˙pi=0˙δqi=δpi˙δpi=0,for1≤i≤N, | (A.3) |
and the corresponding operator eτLAZ, which propagates the values of qi, pi, δqi and δpi for τ time units in the future, obtaining q′i, p′i, δq′i and δp′i, takes the form
eτLAZ:{q′i=qi+τpip′i=piδq′i=δqi+τδpiδp′i=δpi,for1≤i≤N. | (A.4) |
In a similar way for the B(q) Hamiltonian of Eq. (A.2) we get
dXdt=LBZX:{˙qi=0˙pi=(qi+1+qi−1−2qi)+α[(qi+1−qi)2−(qi−qi−1)2]˙δqi=0˙δpi=[2α(qi−1−qi+1)−2]δqi+[1+2α(qi+1−qi)]δqi+1+[1+2α(qi−qi−1)]δqi−1, | (A.5) |
and
eτLBZ:{q′i=qip′i=pi+τ{(qi+1+qi−1−2qi)+α[(qi+1−qi)2−(qi−qi−1)2]}δq′i=δqiδp′i=δpi+τ{[2α(qi−1−qi+1)−2]δqi+[1+2α(qi+1−qi)]δqi+1+[1+2α(qi−qi−1)]δqi−1}. | (A.6) |
According to Eq. (3.26) the accuracy of the SABAn and SBABn integrators can be improved by using a corrector Hamiltonian K [55]. In the case of a separable Hamiltonian H(q,p)=A(p)+B(q) with A(p)=∑Ni=1p2i/2, the corrector K becomes
K(q)={B{B,A}}=N∑i=1 (∂B∂qi)2. | (A.7) |
For the α-FPUT chain, the corrector Hamiltonian K is
K(q)=N∑i=1[(2qi−qi+1−qi−1)(1+α(qi+1−qi−1))]2 . | (A.8) |
As the equations of motion and variational equations associated to the corrector Hamiltonian K are cumbersome, we report here only the form of the operator eτLKZ
eτLKZ:{q′i=qip′i=pi+2τ{2(qi+1+qi−1−2qi)[1+α(qi+1−qi−1)]2−(qi+2+qi−2qi+1)[1+α(qi+2−qi)][1−2α(qi−qi+1)]−(qi−2+qi−2qi−1)[1+α(qi−qi−2)][1−2α(qi−1−qi)]}δq′i=δqiδp′i=δpi+τ{γiδqi+γi+1δqi+1+γi+2δqi+2+γi−1δqi−1+γi−2δqi−2}, | (A.9) |
where
γi=−2{4[1+α(qi+1−qi−1)]2+[1+α(qi+2−qi)][1−2α(qi−qi+1)]+[1+α(qi−qi−2)][1−2α(qi−1−qi)]+α(2qi+1−qi+2−qi)[3−4αqi+2αqi+1+2αqi+2]−α(2qi−1−qi−2−qi)[3+4αqi−2αqi−1−2αqi−2]}γi+1=4{[1+α(qi+1−qi−1)][1−α(4qi−3qi+1−qi−1)]+[1+α(qi+2−qi)][1+α(4qi+1−3qi−qi+2)]}γi−1=4{[1+α(qi+1−qi−1)][1+α(4qi−3qi−1−qi+1)]+[1+α(qi−qi−2)][1−α(4qi−1−3qi−qi−2)]}γi+2=2[1−2α(qi−qi+1)][2α(qi+1−qi+2)−1]γi−2=2[1−2α(qi−1−qi)][2α(qi−2−qi−1)−1]. | (A.10) |
We did not specify the range of index i in Eqs. (A.5), (A.6) and (A.9) intentionally, because it depends on the type of the used boundary conditions. In particular, the expression of Eqs. (A.5), (A.6) and (A.9) are accurate for the case of periodic boundary conditions, i.e., q0=qN, p0=pN, δq0=δqN, δp0=δpN, qN+1=q1, pN+1=p1, δqN+1=δq1, δpN+1=δp1. In the case of fixed boundary conditions we considered in our numerical simulations, some adjustments have to be done for the i=1 and i=N equations, like the ones presented in the Appendix of [56] where the operators for the 1D and 2D DKG models were reported (with the exception of the corrector term).
For completeness sake in Section B.4 we provide the explicit expression of the eτLBZ and eτLKZ operators for some commonly used Hamiltonians, which can be split into two integrable parts, one of which is the usual kinetic energy A(p)=∑Ni=1p2i/2.
Here we focus on the 1D DDNLS system, whose Hamiltonian [Eq. (2.2)] can be split into three integrable parts as
A1=N∑i=1ϵi2(q2i+p2i)+β8(q2i+p2i)2,B1=−N∑i=1pi+1pi ,C1=−N∑i=1qi+1qi. | (A.11) |
The set of equations of motion and variational equations associated with the Hamiltonian function A1 is
dXdt=LA1ZX:{˙qi=piθi˙pi=−qiθi,˙δqi=[θi+βp2i]δpi+βqipiδqi˙δpi=−[θi+βq2i]δqi−βqipiδpi,for1≤i≤N | (A.12) |
with θi=ϵi+β(q2i+p2i)/2 for i=1,2,…,N being constants of the motion. The corresponding operator eτLA1Z takes the form
eτLA1Z:{q′i=qicos(ταi)+pisin(ταi)p′i=picos(ταi)−qisin(ταi)δq′i=qicos(ταi)+pisin(ταi)2JiδJi+(picos(ταi)−qisin(ταi))(βδJiτ+δθi)δp′i=picos(ταi)−qisin(ταi)2JiδJi−(qicos(ταi)+pisin(ταi))(βδJiτ+δθi),for1≤i≤N | (A.13) |
with Ji≠0 and
Ji=12(q2i+p2i) ,αi=ϵi+βJi ,δJi=qiδqi+piδpi ,δθi=pi2Jiδqi−qi2Jiδpi. | (A.14) |
We note that in the special case of Ji=0 we have qi=pi=0. Then the system of Eq. (A.12) takes the simple form ˙qi=0, ˙pi=0, ˙δqi=ϵiδpi, ˙δpi=−ϵiδqi, leading to q′i=qi, p′i=pi, δq′i=δqicos(ϵiτ)+δpisin(ϵiτ), δp′i=δpicos(ϵiτ)−δqisin(ϵiτ).
The set of equations of motion and variational equations associated to the intermediate Hamiltonian functions B1 and C1 are respectively
dXdt=LB1ZX:{˙qi=−pi−1−pi+1˙pi=0˙δqi=−δpi−1−δpi+1˙δpi=0,anddXdt=LC1ZX:{˙qi=0˙pi=qi−1+qi+1˙δqi=0˙δpi=δqi−1+δqi+1. | (A.15) |
These yield to the operators eLB1Z and eLC1Z given by
eτLB1Z:{q′i=qi−τ(pi−1+pi+1)p′i=piδq′i=δqi−τ(δpi−1+δpi+1)δp′i=δpieτLC1Z:{q′i=qip′i=pi+τ(qi−1+qi+1)δq′i=δqiδp′i=δpi+τ(δqi−1+δqi+1). | (A.16) |
As in the case of the α-FPUT model, Eqs. (A.15) and (A.16) correspond to the case of periodic boundary conditions and adjustments similar to the ones presented in the Appendix of [56] for the 1D DKG model, should be implemented when fixed boundary conditions are imposed.
The Hamiltonian H2D of Eq. (2.3) can be split into three integrable parts A2, B2 and C2 as
A2=N∑i=1M∑j=1ϵi,j2[q2i,j+p2i,j]+β8[q2i,j+p2i,j]2,B2=N∑i=1M∑j=1−pi,j+1pi,j−pi+1,jpi,j,C2=N∑i=1M∑j=1−qi,j+1qi,j−qi+1,jqi,j. | (A.17) |
The equations of motion and the variational equations associated with the A2 Hamiltonian are
dXdt=LA2ZX:{˙qi,j=pi,jθi,j˙pi,j=−qi,jθi,j,˙δqi,j=[θi,j+βp2i,j]δpi,j+βqi,jpi,jδqi,j˙δpi,j=−[θi,j+βq2i,j]δqi,j−βqi,jpi,jδpi,j, | (A.18) |
for 1 \leq i \leq N , 1 \leq j \leq M , with \theta_{i, j} = \epsilon _{i, j} + \beta (q_{i, j} ^ 2 + p_{i, j} ^2)/2 being constants of motion for the Hamiltonian \mathcal{A}_2 . Then, the operator e^{\tau L_{\mathcal{A}_2Z}} is
\begin{equation} \begin{split} e^{\tau L_{\mathcal{A}_2Z}} :\left \{ \begin{array}{rl} q_{i, j}' & = q_{i, j} \cos (\tau \alpha _{i, j}) + p_{i, j} \sin (\tau \alpha _{i, j}) \\ p_{i, j}' & = p_{i, j} \cos (\tau \alpha _{i, j}) - q_{i, j} \sin (\tau \alpha _{i, j}) \\ \delta q_{i, j}' & = \frac{q_{i, j} \cos (\tau \alpha _{i, j}) + p_{i, j} \sin (\tau \alpha _{i, j} )}{ 2 J_{i, j}} \delta J_{i, j} + \left(p_{i, j} \cos (\tau \alpha _{i, j} ) - q_{i, j} \sin (\tau \alpha _{i, j} )\right) \left(\beta \delta J_{i, j} \tau + \delta \theta _{i, j} \right) \\ \delta p_{i, j}' & = \frac{p_{i, j} \cos (\tau \alpha _{i, j} ) - q_{i, j} \sin (\tau \alpha _{i, j} )}{2J_{i, j} } \delta J_{i, j} - \left( q_{i, j} \cos ( \tau \alpha _{i, j} ) + p_{i, j} \sin (\tau \alpha _{i, j} ) \right) \left(\beta \delta J_{i, j} \tau + \delta \theta _{i, j} \right) \end{array} \right. \end{split} \end{equation} | (A.19) |
with J_{i, j}\neq 0 and
\begin{equation} \begin{split} &\qquad \quad J_{i, j} = \frac{1}{2} (q_{i, j} ^2 + p_{i, j} ^2)\ , \qquad \alpha _{i, j} = \epsilon _{i, j} + \beta J_{i, j} \ , \\ &\delta J_{i, j} = q_{i, j} \delta q_{i, j} + p_{i, j} \delta p_{i, j} \ , \qquad \delta \theta _{i, j} = \frac{p_{i, j} }{2J_{i, j} }\delta q_{i, j} - \frac{q_{i, j} }{2J_{i, j} } \delta p_{i, j} \end{split}. \end{equation} | (A.20) |
Again, in the special case of J_{i, j} = 0 , the system of Eq. (A.18) takes the form \dot{q_{i, j}} = 0 , \dot{p_{i, j}} = 0 , \dot{\delta q}_{i, j} = \epsilon_{i, j} \delta p_{i, j} , \dot{\delta p}_{i, j} = - \epsilon_{i, j} \delta q_{i, j} , leading to q_{i, j}' = q_{i, j} , p_{i, j}' = p_{i, j} , \delta q_{i, j}' = \delta q_{i, j} \cos (\epsilon_{i, j}\tau) + \delta p_{i, j}\sin (\epsilon_{i, j}\tau) , \delta p_{i, j}' = \delta p_{i, j}\cos (\epsilon_{i, j}\tau) - \delta q_{i, j}\sin (\epsilon_{i, j}\tau) .
The equations of motion and the variational equations for Hamiltonians \mathcal{B}_2 and \mathcal{C}_2 are
\begin{equation} \begin{split} \frac{d\boldsymbol{X}}{dt} = L_{\mathcal{B}_2Z}\boldsymbol{X} &: \left\{ \begin{array}{rl} \dot{q}_{i, j} & = - p_{i-1, j} - p_{i, j-1} - p_{i, j+1} - p_{i+1, j} \\ \dot{p}_{i, j} & = 0 \\ \dot{\delta q}_{i, j} & = - \delta p_{i-1, j} - \delta p_{i, j-1} - \delta p_{i, j+1} - \delta p_{i+1, j} \\ \dot{\delta p}_{i, j} & = 0 \end{array} \right. \end{split}, \end{equation} | (A.21) |
and
\begin{equation} \begin{split} \frac{d\boldsymbol{X}}{dt} = L_{\mathcal{C}_2Z}\boldsymbol{X} & : \left\{ \begin{array}{rl} \dot{p}_{i, j} & = q_{i-1, j} +q_{i, j-1}+q_{i, j+1} +q_{i+1, j} \\ \dot{q}_{i, j} & = 0 \\ \dot{\delta p}_{i, j} & = \delta q_{i-1, j} +\delta q_{i, j-1} +\delta q_{i, j+1} +\delta q_{i+1, j} \\ \dot{\delta q}_{i, j} & = 0 \\ \end{array} \right. \end{split}, \end{equation} | (A.22) |
while the corresponding operators e^{L_{\mathcal{B}_2Z}} and e^{L_{\mathcal{C}_2Z}} are respectively
\begin{equation} \begin{split} e^{\tau L_{\mathcal{B}_2Z}} &:\left \{ \begin{array}{rl} q_{i, j}' & = q_{i, j} - \tau \left( p_{i-1, j} + p_{i, j-1} + p_{i, j+1} + p_{i+1, j} \right) \\ p_{i, j}' & = p_{i, j} \\ \delta q_{i, j}' & = \delta q_{i, j} - \tau \left( \delta p_{i-1, j} + \delta p_{i, j-1} + \delta p_{i, j+1} + \delta p_{i+1, j} \right) \\ \delta p_{i, j}' & = \delta p_{i, j} \\ \end{array} \right. \\ \end{split}, \end{equation} | (A.23) |
and
\begin{equation} \begin{split} e^{\tau L_{\mathcal{C}_2Z}} & : \left \{ \begin{array}{rl} q_{i, j}' & = q_{i, j} \\ p_{i, j}' & = p_{i, j} + \tau \left( q_{i-1, j} +q_{i, j-1}+q_{i, j+1} +q_{i+1, j} \right) \\ \delta q_{i, j}' & = \delta q_{i, j} \\ \delta p_{i, j}' & = \delta p_{i, j} + \tau \left( \delta q_{i-1, j} +\delta q_{i, j-1}+\delta q_{i, j+1} +\delta q_{i+1, j} \right) \end{array} \right.. \end{split} \end{equation} | (A.24) |
Here again Eqs. (A.21)–(A.24) correspond to the case of periodic boundary conditions. For fixed boundary conditions adjustments similar to the ones report in the Appendix of [56] for the 2D DKG lattice should be performed.
We present here the exact expressions of the tangent map operators needed in symplectic integration schemes which can be used to numerically integrate some important models in the field of classical many-body systems: the \beta -FPUT chain, the KG model, and the classical XY model (a JJC system) [53,100,101,102]. Similarly to the \alpha -FPUT chain of Eq. (2.1), the Hamiltonians H(\boldsymbol{q}, \boldsymbol{p}) of each of these systems can be split as
\begin{equation} H(\boldsymbol{q}, \boldsymbol{p}) = A(\boldsymbol{p})+B(\boldsymbol{q}) = \sum\limits_{i = 1}^N \frac{p_i^2}{2} + B\left(\boldsymbol{ q}\right) \end{equation} | (A.25) |
with appropriately defined potential terms B\left(\boldsymbol{ q}\right) :
\begin{align} \text{$\beta$-FPUT:}\ \qquad B_\beta\left(\boldsymbol{ q}\right) & = \sum\limits_{i = 0}^{N} \frac{1}{2}(q_{i+1} - q_i)^2 + \frac{\beta}{4}(q_{i+1} - q_i)^4, \end{align} | (A.26) |
\begin{align} \text{KG:}\ \qquad B_{K}\left(\boldsymbol{ q}\right) & = \sum\limits_{i = 1}^{N} \frac{q_i^2}{2} + \frac{q_i^4 }{4} + \frac{k}{2}(q_{i+1} - q_i)^2 , \end{align} | (A.27) |
\begin{align} \text{JJC:}\ \qquad B_R\left(\boldsymbol{ q}\right) & = \sum\limits_{i = 1}^{N} E_J \left[1-\cos(q_{i+1}-q_i) \right], \end{align} | (A.28) |
where \beta , k and E_J are real coefficients. Obviously for all these systems the operator e^{\tau L_{AZ}} of the kinetic energy part is the same as for the \alpha -FPUT chain in Eq. (A.4). Thus, we report below only the expressions of the operators e^{\tau L_{BZ}} and e^{\tau L_{KZ}} when periodic boundary conditions are imposed.
1. The \beta -Fermi-Pasta-Ulam-Tsingou chain
The operator e^{\tau L_{BZ}} of the \beta -FPUT chain of Eq. (A.26) is
\begin{equation} \begin{split} e^{ \tau L_{BZ}}: \left\{ \begin{array}{rl} q_i' & = q_i \\ p_i' & = \big\{ q_{i-1} + q_{i+1} - 2q_i + \beta \big[ (q_{i+1} - q_i)^3 - (q_i - q_{i-1})^3\big] \big\} \tau + p_i\\ \delta q_i' & = \delta q_i \\ \delta p_i' & = \big\{ \big[- 3 \beta \big[ (q_i - q_{i-1})^2 + (q_{i+1} - q_i)^2\big]- 2 \big]\delta q_i \\ &\qquad \quad + \big[ 1 + 3\beta (q_{i+1} - q_i)^2 \big]\delta q_{i+1} + \big[1 + 3\beta (q_i - q_{i-1})^2\big]\delta q_{i-1} \big\}\tau + \delta p_i \end{array} \right. \end{split}. \end{equation} | (A.29) |
The corrector Hamiltonian K of Eq. (A.7) becomes
\begin{equation} \begin{split} K(\boldsymbol{ q}) & = \sum\limits_{i = 1}^{N}\Big\{ 2q_i - q_{i-1} - q_{i+1} + \beta \big[ (q_i - q_{i-1})^3 - (q_{i+1} - q_i)^3\big] \Big\}^2, \ \\ \end{split} \end{equation} | (A.30) |
while the corresponding operator is
\begin{equation} \begin{split} e^{ \tau L_{KZ}}: \left\{ \begin{array}{rl} q_i' & = q_i \\ p_i' & = 2\Big\{ \Big[ q_{i-1} + q_{i+1} - 2q_i + \beta \big[ (q_{i+1} - q_i)^3 - (q_i - q_{i-1})^3\big] \Big] \Big[2 + 3\beta \big[ (q_i - q_{i-1})^2 + (q_{i+1} - q_i)^2 \big] \Big] \\ &- \Big[q_{i} + q_{i+2} - 2q_{i+1} + \beta \big[(q_{i+2} - q_{i+1})^3 - (q_{i+1} - q_{i})^3\big] \Big] \Big[ 1 + 3 \beta ( q_{i+1} - q_{i})^2 \Big] \\ &- \Big[ q_{i} + q_{i-2} - 2q_{i-1} + \beta \big[ (q_{i} - q_{i-1})^3 - (q_{i-1} - q_{i-2})^3 \big] \Big] \Big[ 1 + 3 \beta ( q_{i} - q_{i-1})^2 \Big] \Big\}\tau + p_i \\ \delta q_i & = \delta q_i \\ \delta p_i & = \big\{ \gamma_i \delta q_{i} + \gamma_{i+1} \delta q_{i+1} + \gamma_{i+2} \delta q_{i+2} + \gamma_{i-1} \delta q_{i-1} + \gamma_{i-2} \delta q_{i-2} \big\}\tau + \delta p_i \end{array} \right. \end{split}, \end{equation} | (A.31) |
with
\begin{equation} \begin{split} \gamma_{i}& = -2 \Big\{ \Big[2 + 3 \beta \big[ (q_i - q_{i-1})^2 +(q_{i+1} - q_i)^2 \big] \Big]^2 + \Big[ 1 + 3 \beta ( q_{i+1} - q_{i})^2 \Big]^2 + \Big[ 1 + 3 \beta ( q_{i} - q_{i-1})^2 \Big]^2 \\ &\qquad\qquad \quad +6\beta \big(2q_i - q_{i-1} - q_{i+1} \big) \Big[2q_i - q_{i-1} - q_{i+1} + \beta \big[ (q_i - q_{i-1})^3 - (q_{i+1} - q_i)^3\big] \Big] \\ &\qquad\qquad \quad- 6\beta ( q_{i} - q_{i+1}) \Big[2q_{i+1} - q_{i} - q_{i+2} + \beta \big[ (q_{i+1} - q_{i})^3 - (q_{i+2} - q_{i+1})^3\big] \Big] \\ &\qquad\qquad \quad-6\beta ( q_{i} - q_{i-1}) \Big[2q_{i-1} - q_{i} - q_{i-2} + \beta \big[ (q_{i-1} - q_{i-2})^3 - (q_{i} - q_{i-1})^3\big] \Big] \Big\} \\ \gamma_{i+1}& = 2\Big\{ \Big[ 1 + 3 \beta (q_{i+1} - q_i)^2 \Big] \Big[4 +3 \beta \big[ (q_i - q_{i-1})^2 + 2(q_{i+1} - q_i)^2 + (q_{i+2} - q_{i+1})^2 \big] \Big] \\ &\qquad\qquad \quad- 6 \beta(q_{i+1} - q_i) \Big[3q_i - q_{i-1} - 3q_{i+1} + q_{i+2} + \beta \big[(q_i - q_{i-1})^3 - 2(q_{i+1} - q_i)^3 + (q_{i+2} - q_{i+1})^3\big] \Big] \Big\} \\ \gamma_{i-1}& = 2 \Big\{ \Big[ 1 + 3 \beta (q_i - q_{i-1})^2 \Big] \Big[4 +3 \beta \big[ (q_{i+1} - q_i)^2 + 2 (q_i - q_{i-1})^2 + (q_{i-1} - q_{i-2})^2 \big] \\ &\qquad\qquad \quad- 6\beta (q_{i-1} - q_{i}) \Big[3q_i - 3q_{i-1} - q_{i+1} + q_{i-2} + \beta \big[ (q_{i} - q_{i+1})^3 - 2 (q_{i-1} - q_{i})^3 + (q_{i-2} - q_{i-1})^3 \big] \Big\} \\ \gamma_{i+2}& = - 2 \Big[ 1 + 3\beta (q_{i+2} - q_{i+1})^2 \Big] \Big[ 1 + 3 \beta ( q_{i+1} - q_{i})^2 \Big] \\ \gamma_{i-2}& = - 2 \Big[ 1 + 3\beta(q_{i-1} - q_{i-2})^2 \Big] \Big[ 1 + 3 \beta ( q_{i} - q_{i-1})^2 \Big] \end{split}. \end{equation} | (A.32) |
2. The 1D Klein-Gordon chain model
The operator e^{\tau L_{BZ}} of the Klein-Gordon chain [Eq. (A.27)] is
\begin{equation} \begin{split} e^{ \tau L_{BZ}}: \left\{ \begin{array}{rl} q_i' & = q_i \\ p_i' & = \big\{ - q_i(1 + q_i^2) + k(q_{i+1} + q_{i-1} -2 q_i ) \big\} \tau + p_i\\ \delta q_i' & = \delta q_i \\ \delta p_i'& = \big\{ - [ 1 + 3q_i^2 + 2k ]\delta q_i + k \delta q_{i+1} + k \delta q_{i-1} \big\}\tau + \delta p_i \end{array} \right. \end{split}. \end{equation} | (A.33) |
The corresponding corrector Hamiltonian K [Eq.(A.7)] is written as
\begin{equation} K(\boldsymbol{ q}) = \sum\limits_{i = 1}^{N}\left[ q_i(1 + q_i^2) + k(2 q_i - q_{i+1} - q_{i-1} ) \right]^2\ , \end{equation} | (A.34) |
while e^{\tau L_{KZ}} takes the form
\begin{equation} \begin{split} e^{ \tau L_{KZ}}: \left\{ \begin{array}{rl} q_i' & = q_i \\ p_i' & = 2 \Big\{ \big[ - q_i(1 + q_i^2) + k( q_{i+1} + q_{i-1} - 2 q_i) \big]\big[ 1 + 3 q_i^2 + 2k \big] \\ &\quad + k\big[ q_{i-1}(1 + q_{i-1}^2) - k( q_{i} + q_{i-2} - 2 q_{i-1}) \big]\\ & \quad + k\big[ q_{i+1}(1 + q_{i+1}^2) - k( q_{i+2} + q_{i} - 2 q_{i+1})\big] \Big\}\tau + p_i \\ \delta q_i' & = \delta q_i \\ \delta p_i' & = \big\{ \gamma_i \delta q_{i} + \gamma_{i+1} \delta q_{i+1} + \gamma_{i+2} \delta q_{i+2} + \gamma_{i-1} \delta q_{i-1} + \gamma_{i-2} \delta q_{i-2} \big\}\tau + \delta p_i \end{array} \right. \end{split}, \end{equation} | (A.35) |
with
\begin{equation} \begin{split} \gamma_{i}& = -2 \Big\{ \big[ 1 + 3 q_i^2 + 2k \big]^2 + 6 q_i \big[ q_i(1 + q_i^2) + k(2 q_i - q_{i+1} - q_{i-1} ) \big] +2k^2 \Big\} \\ \gamma_{i+1}& = 2k \big[ 2 + 4k + 3 q_i^2 + 3q_{i+1}^2 \big] \\ \gamma_{i-1}& = 2k \big[ 2 + 4k + 3 q_i^2 + 3q_{i-1}^2 \big] \\ \gamma_{i+2}& = -2 k^2 \\ \gamma_{i-2}& = -2 k^2 \end{split}. \end{equation} | (A.36) |
3. The XY model of a Josephson junctions array
The operator e^{\tau L_{BZ}} for the potential of Eq. (A.28) is
\begin{equation} \begin{split} e^{ \tau L_{BZ}}: \left\{ \begin{array}{rl} q_i' & = q_i \\ p_i' & = E_J \big[ \sin(q_{i+1}-q_i)+\sin(q_{i-1} - q_i) \big] \tau+ p_i\\ \delta q_i' & = \delta q_i \\ \delta p_i' & = \big\{ - E_J \big[\cos(q_{i+1}-q_i)+\cos(q_i-q_{i-1}) \big] \delta q_i \\ &\qquad E_J \big[ \cos(q_{i+1}-q_i) \big]\delta q_{i+1} + E_J \big[ \cos(q_i-q_{i-1}) \big] \delta q_{i-1} \big\}\tau + \delta p_i \end{array} \right. \end{split}. \end{equation} | (A.37) |
In this case the corrector Hamiltonian K of Eq.(A.7) becomes
\begin{equation} K(\boldsymbol{ q}) = \sum\limits_{i = 1}^{N} E_J^2 \left[ \sin(q_{i+1}-q_i)+\sin(q_{i-1} - q_i) \right]^2, \end{equation} | (A.38) |
and the operator e^{\tau L_{KZ}} is given by the following set of equations
\begin{equation} \begin{split} e^{ \tau L_{KZ}}: \left\{ \begin{array}{rl} q_i' & = q_i \\ p_i' & = E_J^2 \Big\{ 2 \big[ \sin(q_{i+1}-q_i)+\sin(q_{i-1} - q_i) \big] \cdot \big[ \cos(q_{i+1}-q_i) + \cos(q_{i-1} - q_i) \big] \\ & -2 \big[ \sin(q_{i+2}-q_{i+1})+\sin(q_{i} - q_{i+1}) \big] \cdot \cos(q_{i} - q_{i+1}) \\ & - 2 \big[ \sin(q_{i}-q_{i-1})+\sin(q_{i-2} - q_{i-1}) \big] \cdot \cos(q_{i} - q_{i-1}) \Big\}\tau + p_i \\ \delta q_i' & = \delta q_i \\ \delta p_i' & = \big\{ \gamma_i \delta q_{i} + \gamma_{i+1} \delta q_{i+1} + \gamma_{i+2} \delta q_{i+2} + \gamma_{i-1} \delta q_{i-1} + \gamma_{i-2} \delta q_{i-2} \big\}\tau + \delta p_i \end{array} \right. \end{split} \end{equation} | (A.39) |
with
\begin{equation} \begin{split} \gamma_{i}& = E_J^2 \Big\{ - 4 \cos(2 (q_{i+1}-q_i)) - 4 \cos(q_{i-1} -2 q_i+ q_{i+1}) - 4 \cos(2 (q_{i-1}-q_i)) \\ & \qquad + 2 \sin(q_{i+2}-q_{i+1}) \sin(q_{i} - q_{i+1}) + 2 \sin(q_{i-2} - q_{i-1}) \sin(q_{i} - q_{i-1}) \Big\} \\ \gamma_{i+1}& = E_J^2 \Big\{ 2 \cos(q_{i-1} -2 q_i+ q_{i+1}) + 4 \cos(2 (q_{i+1}-q_i)) + 2 \cos(q_{i+2}-2 q_{i+1}+q_{i}) \Big\} \\ \gamma_{i-1}& = E_J^2 \Big\{ 2 \cos(q_{i-1} -2 q_i+ q_{i+1}) + 4\cos(2 (q_{i}-q_{i-1})) + 2 \cos(q_{i-2} - 2 q_{i-1}+ q_{i}) \Big\} \\ \gamma_{i+2}& = E_J^2 2 \cos(q_{i+2}-q_{i+1}) \cos(q_{i} - q_{i+1}) \\ \gamma_{i-2}& = E_J^2 2 \cos(q_{i-2} - q_{i-1}) \cos(q_{i} - q_{i-1}) \end{split}. \end{equation} | (A.40) |
We thank S. Flach for useful discussions. C.D. and M.T. acknowledge financial support from the IBS (Project Code No. IBS-R024-D1). Ch.S. and B.M.M. were supported by the National Research Foundation of South Africa (Incentive Funding for Rated Researchers, IFFR and Competitive Programme for Rated Researchers, CPRR). B.M.M. would also like to thank the High Performance Computing facility of the University of Cape Town (http://hpc.uct.ac.za) and the Center for High Performance Computing (https://www.chpc.ac.za) for the provided computational resources needed for the study of the two DDNLS models, as well as their user-support teams for their help on many practical issues.
All authors declare no conflict of interest in this paper.
1. | Celia Martínez-Lázaro, Marco Antonio Taneco-Hernández, Ramón Reyes-Carreto, Cruz Vargas-De-León, Existence of a Conserved Quantity and Stability of In Vitro Virus Infection Dynamics Models with Absorption Effect, 2019, 2019, 1748-670X, 1, 10.1155/2019/2954041 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_C | Integrator | n | Steps | \tau | T_C |
ABA864 | 4 | 15 | 0.6 | 88 | SRKN_{14}^a | 6 | 29 | 0.45 | 160 |
ABAH864 | 4 | 17 | 0.55 | 115 | SRKN_{11}^b | 6 | 23 | 0.35 | 177 |
SABA_2Y6 | 6 | 29 | 0.575 | 167 | s11SABA_26 | 6 | 45 | 0.3 | 536 |
ABA864Y6 | 6 | 43 | 0.625 | 202 | s19SABA_28 | 8 | 77 | 0.45 | 594 |
s9SABA_26 | 6 | 37 | 0.575 | 205 | s9ABA82\_6 | 6 | 73 | 0.35 | 607 |
FR4 | 4 | 7 | 0.14 | 228 | s15SABA_28 | 8 | 61 | 0.35 | 611 |
SBAB_2Y6 | 6 | 29 | 0.5 | 233 | SABA_2Y6 | 6 | 29 | 0.14 | 683 |
SABA_2K | 4 | 9 | 0.3 | 234 | ABA864 | 4 | 15 | 0.08 | 717 |
ABA82Y4 | 4 | 25 | 0.375 | 240 | s19ABA82\_8 | 8 | 153 | 0.65 | 773 |
s11SABA_26 | 6 | 45 | 0.65 | 247 | s9SABA_26 | 6 | 37 | 0.16 | 779 |
SABA_2Y4 | 4 | 13 | 0.18 | 265 | ABA864Y6 | 6 | 43 | 0.16 | 791 |
ABA82 | 2 | 5 | 0.125 | 278 | s15ABA82\_8 | 8 | 121 | 0.475 | 838 |
ABA82Y6 | 6 | 57 | 0.675 | 283 | ABA82Y6 | 6 | 57 | 0.2 | 841 |
s15SABA_28 | 8 | 61 | 0.65 | 339 | s11ABA82\_6 | 6 | 89 | 0.275 | 941 |
SABA_2 | 2 | 5 | 0.07 | 347 | SBAB_2Y6 | 6 | 29 | 0.12 | 965 |
s19SABA_28 | 8 | 77 | 0.775 | 356 | ABAH864 | 4 | 17 | 0.055 | 1013 |
SBAB_2Y4 | 4 | 13 | 0.18 | 358 | SABA_2Y8\_D | 8 | 61 | 0.175 | 1223 |
FR4Y6 | 6 | 19 | 0.21 | 366 | ABA82Y8\_D | 8 | 121 | 0.25 | 1575 |
s9ABA82\_6 | 6 | 73 | 0.575 | 369 | SABA_2Y4Y6 | 6 | 37 | 0.07 | 1701 |
s11ABA82\_6 | 6 | 89 | 0.675 | 382 | FR4Y6 | 6 | 19 | 0.45 | 1787 |
SBAB_2 | 2 | 5 | 0.07 | 387 | ABA82Y4Y6 | 6 | 73 | 0.125 | 1932 |
SABA_2Y4Y6 | 6 | 37 | 0.3 | 394 | ABA82Y4 | 4 | 25 | 0.0375 | 2156 |
ABA82Y4Y6 | 6 | 73 | 0.525 | 405 | SBAB_2Y4Y6 | 6 | 37 | 0.065 | 2239 |
SABA_2Y8\_D | 8 | 61 | 0.525 | 408 | SABA_2K | 4 | 9 | 0.03 | 2344 |
SBAB_2K | 4 | 9 | 0.2 | 416 | SABA_2KY6 | 6 | 19 | 0.09 | 2465 |
s19ABA82\_8 | 8 | 153 | 1.15 | 439 | FR4 | 4 | 7 | 0.01 | 2597 |
SABA_2KY6 | 6 | 19 | 0.4 | 535 | SABA_2Y4 | 4 | 13 | 0.018 | 2654 |
SBAB_2Y4Y6 | 6 | 37 | 0.275 | 553 | SBAB_2Y4 | 4 | 13 | 0.018 | 3156 |
s15ABA82\_8 | 8 | 121 | 0.775 | 618 | SABA_2Y8\_A | 8 | 61 | 0.06 | 3570 |
ABA82Y8\_D | 8 | 121 | 0.6 | 656 | SBAB_2K | 4 | 9 | 0.02 | 4167 |
SABA_2Y8\_A | 8 | 61 | 0.225 | 1090 | ABA82Y8\_A | 8 | 121 | 0.07 | 5624 |
LF | 2 | 3 | 0.018 | 1198 | ABA82 | 2 | 5 | 0.00125 | 27796 |
ABA82Y8\_A | 8 | 121 | 0.225 | 1749 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 31409 |
SABA_2 | 2 | 5 | 0.0007 | 34595 | |||||
SBAB_2 | 2 | 5 | 0.0007 | 39004 | |||||
LF | 2 | 3 | 0.0002 | 95096 | |||||
TIDES | - | \delta=10^{-16} | 0.05 | 232785 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_C | Integrator | n | Steps | \tau | T_C |
ABA864 | 4 | 15 | 0.6 | 77 | SRKN_{14}^a | 6 | 29 | 0.475 | 156 |
ABAH864 | 4 | 17 | 0.55 | 101 | SRKN_{11}^b | 6 | 23 | 0.35 | 179 |
SABA_2Y6 | 6 | 29 | 0.575 | 183 | s11SABA_26 | 6 | 45 | 0.3 | 480 |
ABA864Y6 | 6 | 43 | 0.625 | 195 | s19SABA_28 | 8 | 77 | 0.45 | 596 |
s9SABA_26 | 6 | 37 | 0.575 | 214 | s9ABA82\_6 | 6 | 73 | 0.35 | 611 |
ABA82Y4 | 4 | 25 | 0.375 | 224 | ABA864 | 4 | 15 | 0.08 | 613 |
s11SABA_26 | 6 | 45 | 0.65 | 227 | s15SABA_28 | 8 | 61 | 0.35 | 618 |
FR4 | 4 | 7 | 0.14 | 231 | SABA_2Y6 | 6 | 29 | 0.14 | 676 |
SABA_2K | 4 | 9 | 0.3 | 241 | s19ABA82\_8 | 8 | 153 | 0.65 | 760 |
SBAB_2Y6 | 6 | 29 | 0.5 | 255 | ABA864Y6 | 6 | 43 | 0.16 | 811 |
SABA_2Y4 | 4 | 13 | 0.18 | 270 | s15ABA82\_8 | 8 | 121 | 0.475 | 828 |
ABA82 | 2 | 5 | 0.125 | 280 | s9SABA_26 | 6 | 37 | 0.16 | 838 |
ABA82Y6 | 6 | 57 | 0.675 | 285 | ABA82Y6 | 6 | 57 | 0.2 | 937 |
SBAB_2Y4 | 4 | 13 | 0.18 | 316 | SBAB_2Y6 | 6 | 29 | 0.12 | 964 |
s15SABA_28 | 8 | 61 | 0.65 | 329 | ABAH864 | 4 | 17 | 0.055 | 1023 |
s19SABA_28 | 8 | 77 | 0.775 | 336 | s11ABA82\_6 | 6 | 89 | 0.275 | 1062 |
FR4Y6 | 6 | 19 | 0.21 | 337 | SABA_2Y8\_D | 8 | 61 | 0.175 | 1230 |
SBAB_2K | 4 | 9 | 0.2 | 366 | FR4Y6 | 6 | 19 | 0.45 | 1577 |
s9ABA82\_6 | 6 | 73 | 0.575 | 373 | ABA82Y8\_D | 8 | 121 | 0.25 | 1613 |
SABA_2 | 2 | 5 | 0.07 | 392 | ABA82Y4Y6 | 6 | 73 | 0.125 | 1702 |
SBAB_2 | 2 | 5 | 0.07 | 393 | ABA82Y4 | 4 | 25 | 0.0375 | 2113 |
ABA82Y4Y6 | 6 | 73 | 0.525 | 398 | SABA_2Y4Y6 | 6 | 37 | 0.07 | 2159 |
SABA_2Y8\_D | 8 | 61 | 0.525 | 407 | SBAB_2Y4Y6 | 6 | 37 | 0.065 | 2310 |
s11ABA82\_6 | 6 | 89 | 0.675 | 415 | SABA_2K Y6 | 6 | 19 | 0.09 | 2380 |
s19ABA82\_8 | 8 | 153 | 1.15 | 431 | FR4 | 4 | 7 | 0.01 | 2615 |
SABA_2Y4Y6 | 6 | 37 | 0.3 | 444 | SABA_2K | 4 | 9 | 0.03 | 2651 |
SBAB_2Y4Y6 | 6 | 37 | 0.275 | 533 | SABA_2Y4 | 4 | 13 | 0.018 | 3016 |
SABA_2K Y6 | 6 | 19 | 0.4 | 540 | SBAB_2Y4 | 4 | 13 | 0.018 | 3585 |
s15ABA82\_8 | 8 | 121 | 0.775 | 598 | SABA_2Y8\_A | 8 | 61 | 0.06 | 3647 |
ABA82Y8\_D | 8 | 121 | 0.6 | 629 | SBAB_2K | 4 | 9 | 0.02 | 3663 |
SABA_2Y8\_A | 8 | 61 | 0.225 | 952 | ABA82Y8\_A | 8 | 121 | 0.07 | 5691 |
LF | 2 | 3 | 0.018 | 1059 | ABA82 | 2 | 5 | 0.00125 | 31576 |
ABA82Y8\_A | 8 | 121 | 0.225 | 1787 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 31709 |
SABA_2 | 2 | 5 | 0.0007 | 39333 | |||||
SBAB_2 | 2 | 5 | 0.0007 | 45690 | |||||
LF | 2 | 3 | 0.0002 | 106653 | |||||
TIDES | - | \delta=10^{-16} | 0.05 | 225565 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s11\mathcal{ABC}6 | 6 | 45 | 0.115 | 3395 | s19\mathcal{ABC}8 | 8 | 77 | 0.09 | 7242 |
s9\mathcal{ABC}6 | 6 | 37 | 0.095 | 3425 | s17\mathcal{ABC}8 | 8 | 69 | 0.08 | 7301 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.07 | 3720 | s11\mathcal{ABC}6 | 6 | 45 | 0.025 | 15692 |
SS864S | 4 | 17 | 0.05 | 6432 | s9\mathcal{ABC}6 | 6 | 37 | 0.02 | 16098 |
\mathcal{ABC}Y4 | 4 | 13 | 0.0125 | 10317 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 18408 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.015 | 35417 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.002 | 258891 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.008 | 40109 | TIDES | - | \delta = 10^{-16} | 0.05 | 419958 |
\mathcal{ABC}S4 | 4 | 21 | 0.00085 | 267911 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.0002 | 320581 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s11\mathcal{ABC}6 | 6 | 45 | 0.4 | 1132 | s19\mathcal{ABC}8 | 8 | 77 | 0.3 | 2184 |
s9\mathcal{ABC}6 | 6 | 37 | 0.285 | 1147 | s17\mathcal{ABC}8 | 8 | 69 | 0.225 | 2632 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.2 | 1308 | s11\mathcal{ABC}6 | 6 | 45 | 0.1 | 4137 |
SS864S | 4 | 17 | 0.265 | 1365 | s9\mathcal{ABC}6 | 6 | 37 | 0.075 | 4462 |
\mathcal{ABC}Y4 | 4 | 13 | 0.055 | 2354 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.065 | 8528 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.105 | 4965 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 14998 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.04 | 8091 | TIDES | - | \delta = 10^{-16} | 0.05 | 420050 |
\mathcal{ABC}S4 | 4 | 21 | 0.02 | 9774 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.0055 | 11700 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s9\mathcal{ABC}6 | 6 | 45 | 0.105 | 13914 | s17\mathcal{ABC}8 | 8 | 77 | 0.075 | 36528 |
s11\mathcal{ABC}6 | 6 | 37 | 0.125 | 14000 | s19\mathcal{ABC}8 | 8 | 69 | 0.08 | 38270 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.08 | 15344 | s9\mathcal{ABC}6 | 6 | 45 | 0.0235 | 65287 |
\mathcal{ABC}Y4 | 4 | 13 | 0.025 | 23030 | s11\mathcal{ABC}6 | 6 | 37 | 0.0275 | 67314 |
SS864S | 4 | 17 | 0.085 | 23887 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.008 | 140506 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.03 | 77424 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 218704 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.0165 | 87902 | |||||
\mathcal{ABC}S4 | 4 | 21 | 0.0065 | 132713 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.005 | 157694 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s11\mathcal{ABC}6 | 6 | 45 | 0.1515 | 11443 | s19\mathcal{ABC}8 | 8 | 77 | 0.135 | 20952 |
s9\mathcal{ABC}6 | 6 | 37 | 0.11 | 13408 | s17\mathcal{ABC}8 | 8 | 69 | 0.0875 | 28966 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.0775 | 14607 | s11\mathcal{ABC}6 | 6 | 45 | 0.0335 | 50301 |
SS864S | 4 | 17 | 0.0915 | 15564 | s9\mathcal{ABC}6 | 6 | 37 | 0.024 | 58187 |
\mathcal{ABC}Y4 | 4 | 13 | 0.0215 | 25898 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.009 | 150045 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.035 | 64423 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 166145 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.01375 | 102580 | |||||
\mathcal{ABC}S4 | 4 | 21 | 0.005 | 185615 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.00155 | 198534 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_C | Integrator | n | Steps | \tau | T_C |
ABA864 | 4 | 15 | 0.6 | 88 | SRKN_{14}^a | 6 | 29 | 0.45 | 160 |
ABAH864 | 4 | 17 | 0.55 | 115 | SRKN_{11}^b | 6 | 23 | 0.35 | 177 |
SABA_2Y6 | 6 | 29 | 0.575 | 167 | s11SABA_26 | 6 | 45 | 0.3 | 536 |
ABA864Y6 | 6 | 43 | 0.625 | 202 | s19SABA_28 | 8 | 77 | 0.45 | 594 |
s9SABA_26 | 6 | 37 | 0.575 | 205 | s9ABA82\_6 | 6 | 73 | 0.35 | 607 |
FR4 | 4 | 7 | 0.14 | 228 | s15SABA_28 | 8 | 61 | 0.35 | 611 |
SBAB_2Y6 | 6 | 29 | 0.5 | 233 | SABA_2Y6 | 6 | 29 | 0.14 | 683 |
SABA_2K | 4 | 9 | 0.3 | 234 | ABA864 | 4 | 15 | 0.08 | 717 |
ABA82Y4 | 4 | 25 | 0.375 | 240 | s19ABA82\_8 | 8 | 153 | 0.65 | 773 |
s11SABA_26 | 6 | 45 | 0.65 | 247 | s9SABA_26 | 6 | 37 | 0.16 | 779 |
SABA_2Y4 | 4 | 13 | 0.18 | 265 | ABA864Y6 | 6 | 43 | 0.16 | 791 |
ABA82 | 2 | 5 | 0.125 | 278 | s15ABA82\_8 | 8 | 121 | 0.475 | 838 |
ABA82Y6 | 6 | 57 | 0.675 | 283 | ABA82Y6 | 6 | 57 | 0.2 | 841 |
s15SABA_28 | 8 | 61 | 0.65 | 339 | s11ABA82\_6 | 6 | 89 | 0.275 | 941 |
SABA_2 | 2 | 5 | 0.07 | 347 | SBAB_2Y6 | 6 | 29 | 0.12 | 965 |
s19SABA_28 | 8 | 77 | 0.775 | 356 | ABAH864 | 4 | 17 | 0.055 | 1013 |
SBAB_2Y4 | 4 | 13 | 0.18 | 358 | SABA_2Y8\_D | 8 | 61 | 0.175 | 1223 |
FR4Y6 | 6 | 19 | 0.21 | 366 | ABA82Y8\_D | 8 | 121 | 0.25 | 1575 |
s9ABA82\_6 | 6 | 73 | 0.575 | 369 | SABA_2Y4Y6 | 6 | 37 | 0.07 | 1701 |
s11ABA82\_6 | 6 | 89 | 0.675 | 382 | FR4Y6 | 6 | 19 | 0.45 | 1787 |
SBAB_2 | 2 | 5 | 0.07 | 387 | ABA82Y4Y6 | 6 | 73 | 0.125 | 1932 |
SABA_2Y4Y6 | 6 | 37 | 0.3 | 394 | ABA82Y4 | 4 | 25 | 0.0375 | 2156 |
ABA82Y4Y6 | 6 | 73 | 0.525 | 405 | SBAB_2Y4Y6 | 6 | 37 | 0.065 | 2239 |
SABA_2Y8\_D | 8 | 61 | 0.525 | 408 | SABA_2K | 4 | 9 | 0.03 | 2344 |
SBAB_2K | 4 | 9 | 0.2 | 416 | SABA_2KY6 | 6 | 19 | 0.09 | 2465 |
s19ABA82\_8 | 8 | 153 | 1.15 | 439 | FR4 | 4 | 7 | 0.01 | 2597 |
SABA_2KY6 | 6 | 19 | 0.4 | 535 | SABA_2Y4 | 4 | 13 | 0.018 | 2654 |
SBAB_2Y4Y6 | 6 | 37 | 0.275 | 553 | SBAB_2Y4 | 4 | 13 | 0.018 | 3156 |
s15ABA82\_8 | 8 | 121 | 0.775 | 618 | SABA_2Y8\_A | 8 | 61 | 0.06 | 3570 |
ABA82Y8\_D | 8 | 121 | 0.6 | 656 | SBAB_2K | 4 | 9 | 0.02 | 4167 |
SABA_2Y8\_A | 8 | 61 | 0.225 | 1090 | ABA82Y8\_A | 8 | 121 | 0.07 | 5624 |
LF | 2 | 3 | 0.018 | 1198 | ABA82 | 2 | 5 | 0.00125 | 27796 |
ABA82Y8\_A | 8 | 121 | 0.225 | 1749 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 31409 |
SABA_2 | 2 | 5 | 0.0007 | 34595 | |||||
SBAB_2 | 2 | 5 | 0.0007 | 39004 | |||||
LF | 2 | 3 | 0.0002 | 95096 | |||||
TIDES | - | \delta=10^{-16} | 0.05 | 232785 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_C | Integrator | n | Steps | \tau | T_C |
ABA864 | 4 | 15 | 0.6 | 77 | SRKN_{14}^a | 6 | 29 | 0.475 | 156 |
ABAH864 | 4 | 17 | 0.55 | 101 | SRKN_{11}^b | 6 | 23 | 0.35 | 179 |
SABA_2Y6 | 6 | 29 | 0.575 | 183 | s11SABA_26 | 6 | 45 | 0.3 | 480 |
ABA864Y6 | 6 | 43 | 0.625 | 195 | s19SABA_28 | 8 | 77 | 0.45 | 596 |
s9SABA_26 | 6 | 37 | 0.575 | 214 | s9ABA82\_6 | 6 | 73 | 0.35 | 611 |
ABA82Y4 | 4 | 25 | 0.375 | 224 | ABA864 | 4 | 15 | 0.08 | 613 |
s11SABA_26 | 6 | 45 | 0.65 | 227 | s15SABA_28 | 8 | 61 | 0.35 | 618 |
FR4 | 4 | 7 | 0.14 | 231 | SABA_2Y6 | 6 | 29 | 0.14 | 676 |
SABA_2K | 4 | 9 | 0.3 | 241 | s19ABA82\_8 | 8 | 153 | 0.65 | 760 |
SBAB_2Y6 | 6 | 29 | 0.5 | 255 | ABA864Y6 | 6 | 43 | 0.16 | 811 |
SABA_2Y4 | 4 | 13 | 0.18 | 270 | s15ABA82\_8 | 8 | 121 | 0.475 | 828 |
ABA82 | 2 | 5 | 0.125 | 280 | s9SABA_26 | 6 | 37 | 0.16 | 838 |
ABA82Y6 | 6 | 57 | 0.675 | 285 | ABA82Y6 | 6 | 57 | 0.2 | 937 |
SBAB_2Y4 | 4 | 13 | 0.18 | 316 | SBAB_2Y6 | 6 | 29 | 0.12 | 964 |
s15SABA_28 | 8 | 61 | 0.65 | 329 | ABAH864 | 4 | 17 | 0.055 | 1023 |
s19SABA_28 | 8 | 77 | 0.775 | 336 | s11ABA82\_6 | 6 | 89 | 0.275 | 1062 |
FR4Y6 | 6 | 19 | 0.21 | 337 | SABA_2Y8\_D | 8 | 61 | 0.175 | 1230 |
SBAB_2K | 4 | 9 | 0.2 | 366 | FR4Y6 | 6 | 19 | 0.45 | 1577 |
s9ABA82\_6 | 6 | 73 | 0.575 | 373 | ABA82Y8\_D | 8 | 121 | 0.25 | 1613 |
SABA_2 | 2 | 5 | 0.07 | 392 | ABA82Y4Y6 | 6 | 73 | 0.125 | 1702 |
SBAB_2 | 2 | 5 | 0.07 | 393 | ABA82Y4 | 4 | 25 | 0.0375 | 2113 |
ABA82Y4Y6 | 6 | 73 | 0.525 | 398 | SABA_2Y4Y6 | 6 | 37 | 0.07 | 2159 |
SABA_2Y8\_D | 8 | 61 | 0.525 | 407 | SBAB_2Y4Y6 | 6 | 37 | 0.065 | 2310 |
s11ABA82\_6 | 6 | 89 | 0.675 | 415 | SABA_2K Y6 | 6 | 19 | 0.09 | 2380 |
s19ABA82\_8 | 8 | 153 | 1.15 | 431 | FR4 | 4 | 7 | 0.01 | 2615 |
SABA_2Y4Y6 | 6 | 37 | 0.3 | 444 | SABA_2K | 4 | 9 | 0.03 | 2651 |
SBAB_2Y4Y6 | 6 | 37 | 0.275 | 533 | SABA_2Y4 | 4 | 13 | 0.018 | 3016 |
SABA_2K Y6 | 6 | 19 | 0.4 | 540 | SBAB_2Y4 | 4 | 13 | 0.018 | 3585 |
s15ABA82\_8 | 8 | 121 | 0.775 | 598 | SABA_2Y8\_A | 8 | 61 | 0.06 | 3647 |
ABA82Y8\_D | 8 | 121 | 0.6 | 629 | SBAB_2K | 4 | 9 | 0.02 | 3663 |
SABA_2Y8\_A | 8 | 61 | 0.225 | 952 | ABA82Y8\_A | 8 | 121 | 0.07 | 5691 |
LF | 2 | 3 | 0.018 | 1059 | ABA82 | 2 | 5 | 0.00125 | 31576 |
ABA82Y8\_A | 8 | 121 | 0.225 | 1787 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 31709 |
SABA_2 | 2 | 5 | 0.0007 | 39333 | |||||
SBAB_2 | 2 | 5 | 0.0007 | 45690 | |||||
LF | 2 | 3 | 0.0002 | 106653 | |||||
TIDES | - | \delta=10^{-16} | 0.05 | 225565 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s11\mathcal{ABC}6 | 6 | 45 | 0.115 | 3395 | s19\mathcal{ABC}8 | 8 | 77 | 0.09 | 7242 |
s9\mathcal{ABC}6 | 6 | 37 | 0.095 | 3425 | s17\mathcal{ABC}8 | 8 | 69 | 0.08 | 7301 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.07 | 3720 | s11\mathcal{ABC}6 | 6 | 45 | 0.025 | 15692 |
SS864S | 4 | 17 | 0.05 | 6432 | s9\mathcal{ABC}6 | 6 | 37 | 0.02 | 16098 |
\mathcal{ABC}Y4 | 4 | 13 | 0.0125 | 10317 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 18408 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.015 | 35417 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.002 | 258891 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.008 | 40109 | TIDES | - | \delta = 10^{-16} | 0.05 | 419958 |
\mathcal{ABC}S4 | 4 | 21 | 0.00085 | 267911 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.0002 | 320581 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s11\mathcal{ABC}6 | 6 | 45 | 0.4 | 1132 | s19\mathcal{ABC}8 | 8 | 77 | 0.3 | 2184 |
s9\mathcal{ABC}6 | 6 | 37 | 0.285 | 1147 | s17\mathcal{ABC}8 | 8 | 69 | 0.225 | 2632 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.2 | 1308 | s11\mathcal{ABC}6 | 6 | 45 | 0.1 | 4137 |
SS864S | 4 | 17 | 0.265 | 1365 | s9\mathcal{ABC}6 | 6 | 37 | 0.075 | 4462 |
\mathcal{ABC}Y4 | 4 | 13 | 0.055 | 2354 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.065 | 8528 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.105 | 4965 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 14998 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.04 | 8091 | TIDES | - | \delta = 10^{-16} | 0.05 | 420050 |
\mathcal{ABC}S4 | 4 | 21 | 0.02 | 9774 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.0055 | 11700 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s9\mathcal{ABC}6 | 6 | 45 | 0.105 | 13914 | s17\mathcal{ABC}8 | 8 | 77 | 0.075 | 36528 |
s11\mathcal{ABC}6 | 6 | 37 | 0.125 | 14000 | s19\mathcal{ABC}8 | 8 | 69 | 0.08 | 38270 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.08 | 15344 | s9\mathcal{ABC}6 | 6 | 45 | 0.0235 | 65287 |
\mathcal{ABC}Y4 | 4 | 13 | 0.025 | 23030 | s11\mathcal{ABC}6 | 6 | 37 | 0.0275 | 67314 |
SS864S | 4 | 17 | 0.085 | 23887 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.008 | 140506 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.03 | 77424 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 218704 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.0165 | 87902 | |||||
\mathcal{ABC}S4 | 4 | 21 | 0.0065 | 132713 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.005 | 157694 |
E_r \approx 10^{-5} | E_r \approx 10^{-9} | ||||||||
Integrator | n | Steps | \tau | T_{\text{C}} | Integrator | n | Steps | \tau | T_{\text{C}} |
s11\mathcal{ABC}6 | 6 | 45 | 0.1515 | 11443 | s19\mathcal{ABC}8 | 8 | 77 | 0.135 | 20952 |
s9\mathcal{ABC}6 | 6 | 37 | 0.11 | 13408 | s17\mathcal{ABC}8 | 8 | 69 | 0.0875 | 28966 |
\mathcal{ABC}Y6\_A | 6 | 29 | 0.0775 | 14607 | s11\mathcal{ABC}6 | 6 | 45 | 0.0335 | 50301 |
SS864S | 4 | 17 | 0.0915 | 15564 | s9\mathcal{ABC}6 | 6 | 37 | 0.024 | 58187 |
\mathcal{ABC}Y4 | 4 | 13 | 0.0215 | 25898 | \mathcal{ABC}Y8\_D | 8 | 61 | 0.009 | 150045 |
\mathcal{ABC}S4Y6 | 6 | 49 | 0.035 | 64423 | DOP853 | 8 | \delta=10^{-16} | 0.05 | 166145 |
\mathcal{ABC}Y4Y6 | 6 | 37 | 0.01375 | 102580 | |||||
\mathcal{ABC}S4 | 4 | 21 | 0.005 | 185615 | |||||
\mathcal{ABC}2 | 2 | 5 | 0.00155 | 198534 |