Efficient integration of the variational equations of multi-dimensional Hamiltonian systems: Application to the Fermi-Pasta-Ulam lattice
We study the problem of efficient integration of variational equations in multi-dimensional Hamiltonian systems. For this purpose, we consider a Runge-Kutta-type integrator, a Taylor series expansion method and the so-called ‘Tangent Map’ (TM) technique based on symplectic integration schemes, and apply them to the Fermi-Pasta-Ulam (FPU-) lattice of nonlinearly coupled oscillators, with ranging from 4 to 20. The fast and accurate reproduction of well-known behaviors of the Generalized Alignment Index (GALI) chaos detection technique is used as an indicator for the efficiency of the tested integration schemes. Implementing the TM technique–which shows the best performance among the tested algorithms–and exploiting the advantages of the GALI method, we successfully trace the location of low-dimensional tori.
From interactions of stars in galaxies to particle beams in high energy accelerators, Hamiltonian mechanics is found at the very heart of modeling and understanding dynamical processes. The necessity to classify evermore complex systems in terms of stability and predictability has lead to a wealth of methods discriminating chaotic from regular behavior (see for example [Skokos, 2010, Sect. 7]). Most of these techniques rely on the study of the time evolution of deviation vectors of a given orbit to discriminate between the two regimes. The time evolution of these vectors is governed by the so-called variational equations.
In Skokos & Gerlach  the ‘Tangent Map’ (TM) technique, an efficient and easy to implement method for the integration of the variational equations of Hamiltonian systems based on the use of symplectic integrators was introduced. In Skokos & Gerlach ; Gerlach & Skokos  the TM method was applied mainly to low-dimensional Hamiltonian systems of 2 and 3 degrees of freedom, and proved to be very efficient and superior to other commonly used numerical schemes, both with respect to its accuracy and its speed.
The scope of the present work is to extend these results by investigating whether the efficiency of the TM method persists also when multi-dimensional Hamiltonian systems are considered. The study of such systems presents a challenging numerical task, which makes the use of fast and accurate numerical tools imperative. In the present paper we use as a toy model the famous Fermi-Pasta-Ulam (FPU-) lattice, which is presented in Sect. 2. In Sect. 3 the different numerical methods we use, i.e. the Generalized Alignment Index (GALI) chaos indicator, the TM method and the SABA family of symplectic integrators, the Taylor series integrator TIDES, and the general-purpose high-accuracy Runge-Kutta integrator DOP853, are presented and their properties are briefly discussed. Then, in Sect. 4.1 the numerical results of the application of these numerical schemes for the integration of variational equations of the FPU system are presented, while in Sect. 4.2 the GALI method is used for locating motion on low-dimensional tori. Finally, in Sect. 5 we summarize our results.
2 The FPU lattice
As a model of a multi-dimensional Hamiltonian system we consider the FPU- lattice Fermi et al. ; Ford ; Chaos Focus Issue , which describes a chain of particles with nearest neighbor interaction. Regarding numerical integration algorithms, the FPU lattice is a very challenging problem, since it exhibits oscillations on largely different time scales. The Hamiltonian of this degrees of freedom (D) system as a function of the momenta and the coordinates is given by
In our study we impose fixed boundary conditions, i. e. , set , and consider models whose number of particles vary from up to . We note that in Skokos et al. ; Skokos & Gerlach  the particular case was studied in detail.
The Hamiltonian (1) can be split into two parts and , which respectively depend only on the momenta and the coordinates, i. e. . The Hamilton’s equations of motion are
with , and denoting the Kronecker delta, which is equal to 1 if and to 0 otherwise. The variational equations governing the time evolution of a deviation vector , that evolves in the tangent space of the Hamiltonian’s phase space are given by
3 Numerical methods
In order to compute the stability of a particular solution of Hamiltonian (1), or in other words of an orbit in the 2-dimensional phase space of the system, the equations of motion (2) have to be integrated together with the variational equations (3). The time evolution of the latter contains information on the stability of the orbit, which can be quantified by using some chaos indicator, for example the GALIs. Different numerical approaches can be used to solve the system of ordinary differential equations given by Eqs. (2) and (3). In this section, after briefly recalling the definition of the GALI and its behavior for regular and chaotic motion, we present an overview of the methods used in the current study. A more detailed description of further possibilities based on symplectic methods can be found in Skokos & Gerlach .
3.1 The Generalized Alignment Index (GALI)
The GALI was originally introduced in Skokos et al.  as an efficient chaos detection method, generalizing a similar indicator called the Smaller Alignment Index (SALI) Skokos ; Skokos et al. [2003, 2004]. The method has been applied successfully to different dynamical systems for the discrimination between regular and chaotic motion, as well as for the detection of regular motion on low dimensional tori Christodoulidi & Bountis ; Skokos et al. ; Bountis et al. ; Manos & Ruffo ; Manos & Athanassoula ; Manos et al. .
For D Hamiltonians the Generalized Alignment Index of order (GALI), , is determined through the evolution of initially linearly independent deviation vectors , which are continually normalized, keeping their directions intact. According to Skokos et al.  GALI is defined as the volume of the -parallelepiped having the unitary deviation vectors , , as edges. GALI is therefore determined through the wedge product of these vectors
while for chaotic orbits GALI tends to zero exponentially following the law Skokos et al. 
where are the first largest Lyapunov characteristic exponents of the orbit.
We chose to apply the GALI method in order to check the accuracy of the numerical integration of variational equations because this index depends on the evolution of an ensemble of deviation vectors. Thus, even small errors in the integration of these vectors are expected to significantly influence the evolution of GALIs. In particular, we focus our attention on regular orbits lying on tori of various dimensions , for which the GALIs exhibit different evolutions depending on and the order of the index. The possible deviations of numerically evaluated GALIs from their expected behaviors (5) will identify the inaccurate integration of the deviation vectors far better than it would be the case using chaotic orbits, for which all GALIs tend to zero exponentially (6). Throughout the paper we denote regular orbits on an -dimensional torus of the D system (1) as , where can range from 2 to .
3.2 The Tangent Map method using symplectic algorithms
Symplectic methods are often the preferred choice when integrating dynamical problems, which can be described by Hamiltonian functions. A thorough discussion of such methods can be found in Hairer et al. . Let us just mention some properties of symplectic integrators which are of interest for our study. Symplectic methods cannot be used with a trivial automated step size control. As a consequence, they are usually implemented with a fixed integration step . Due to their special structure they preserve the symplectic nature of Hamilton’s equations intrinsically, which in turn leads to results that are more robust for long integration times. A side-effect of structure preservation is that the error in energy remains bounded irrespective of the total integration time.
In Skokos & Gerlach  it was shown that it is possible to integrate the Hamilton’s equations of motion and the corresponding variational equations using the TM technique based on symplectic splitting methods. Let us outline the basic idea behind the TM method, which is founded on a general result stated for example in Laskar & Robutel : Symplectic integrators can be applied to systems of first order differential equations , that can be written in the form , where are differential operators defined as and for which the two systems and are integrable. Here are Poisson brackets of functions , defined as:
The set of Eqs. (2) is one example of such a system, since the Hamiltonian (1) can be divided into two integrable parts and with as already noted. A symplectic splitting method separates the equations of motion (2) into several parts, applying either the operator or . These are the equations of motion of the Hamiltonians and , which can be solved analytically, giving explicit mappings over the time step , where the constants are chosen to optimize the accuracy of the integrator. These mappings can then be combined to approximate the solution for a time step . In Skokos & Gerlach  it was shown that the derivatives of these mappings with respect to the coordinates and momenta of the system (the so-called tangent maps) can be used to calculate the time evolution of deviation vectors or, in other words, solve the variational equations (3).
In Laskar & Robutel  a family of symplectic splitting methods called SABA and SBAB was introduced, with being the number of applications of operators and . These integrators were designed to have only positive intermediate steps. Since it is not possible to construct symplectic integrators of order111In this work we call a symplectic integrator to be of order , if it introduces an error of in the approximation of the real solution, with being the integration time step. with this property Suzuki , small negative corrector steps can be added before and after the main body of the integrator to further increase the accuracy. In Sect. 4.1 we test 3 different integrators, namely SABA, SABAC and SBABC, present the results of this comparison and discuss the characteristics of these algorithms when applied to the FPU system.
The fact, that standard adaptive stepping policy is not possible with symplectic integration schemes necessitates an initial assessment of stability for the algorithms used, in order to derive an upper bound on the choice of step-size. Following Hairer et al.  we note that SBAB is equivalent to the well known Størmer-Verlet or leap-frog method, and SABA to its adjoint. In order to perform a linear stability analysis of the FPU- lattice problem for SABA, we introduce normal mode momenta and coordinates as it was done for example in Skokos et al. . The unperturbed () Hamiltonian (1) can then be written as a sum of the so-called harmonic energies , i.e.
where are the corresponding harmonic frequencies. With denoting the Jacobian of the numerical mapping from initial coordinates and momenta to those at the next step for SABA we have
The characteristic polynomial of matrix amounts to . For the eigenvalues to be of modulus one, given the harmonic frequencies (8) satisfy , the maximum admissible step-size for SABA will be . Thus, in our study we always use .
3.3 Taylor series methods-TIDES
The basic idea of the so-called Taylor series methods is to approximate the solution at time of a given -dimensional initial value problem
from the degree Taylor series of at :
For details see e.g. [Hairer et al., 1993, Sect. I.8] and references therein. In this work we call an integrator being of order , when the first neglected term in this Taylor series expansion is of . The computation of the derivatives can be very cumbersome, depending on the structure of , and is done efficiently using automatic differentiation techniques (see for example Barrio ).
In Gerlach & Skokos  two different publicly available implementations of the Taylor method were used and compared regarding their reliability and efficiency in the case of a D Hamiltonian system. One of these integrators, called TIDES222Freely available at http://gme.unizar.es/software/tides. Barrio ; Abad et al. , which showed better performance, will also be used in this work as a representative of methods based on Taylor series expansions. TIDES comes as a Mathematica notebook. After inserting the differential equations one desires to integrate, the notebook generates automatically all the necessary subroutines to compute the given problem. The FORTRAN code produced by TIDES was included into the existing testbed of our work without further modification.
In order to obtain optimal results, the TIDES algorithm is free to choose its order and step size during the whole integration interval. We used one parameter only, the so-called one step accuracy , to control the numerical performance of the algorithm. To be more precise, an integration step from time to is accepted, if the internal accuracy checks estimate that the local truncation error of the solution is less than . If this error is too large, the integrator automatically tries to increase the internal order and/or adjust the step size . After each successful step the deviation vectors are re-normalized.
Let us remark that another elegant way to express Eq. (10) can be achieved using the so called Lie series formalism. Lie series have been rediscovered by Gröbner  and used extensively in the field of dynamical astronomy (e.g. Hanslmeier & Dvorak ; Delva ) to numerically solve ordinary differential equations. Sharing the same ansatz with symplectic maps (Sect. 3.2), Lie series can be used to iterate first order ordinary differential equations of the form as follows:
Note that contrary to Sect. 3.2 no assumptions on the properties of the differential operator are made. Thus, by evaluating the consecutive derivatives and building the corresponding exponential series up to order one is able to follow the trajectory of through phase space. The truncated Lie series’ numerical map is not area-preserving, in general, since the method is non-symplectic. Therefore, the implementation of an adaptive step-size into Lie series algorithms is possible without reservations. This can be achieved via estimates on the size of the remainder of the exponential series (see for example Eggl & Dvorak ).
Algorithms using series expansions become most efficient when combined with recursion relations, where higher order derivatives can be calculated from previous ones. In the course of this work, also extensive tests have been undertaken to adapt the Lie series method to the FPU- lattice. Since for this problem such recursions are not available, we used algebraic manipulation software to compute the successive applications of the operator to , and implemented the generated code into a FORTRAN program. Since this approach proved to be reliable but computationally very expensive, we do not include it in the discussion of our results in Sect. 4.
3.4 General purpose integrators-DOP853
In general, the computation of higher derivatives of functions soon becomes very complicated. Therefore, the methods described in Sect. 3.3 became popular only recently, after automatic differentiation could be performed efficiently by computers. Before that, other methods were developed to approximate Eq. (10). One of these are the so called Runge-Kutta methods (see for example [Hairer et al., 1993, Sect. II.1.1]). An -stage Runge-Kutta method is given as
The real numbers with are chosen to approximate Eq. (10) to the desired order. If one requires further for the integration method will be explicit. Runge-Kutta integrators exist as symplectic as well as non-symplectic variants.
In this work a 12-stage explicit Runge-Kutta integration method called DOP853 is used333DOP853 is freely available from http://www.unige.ch/~hairer/software.html.. This non-symplectic scheme is based on the method of Dormand and Price (see [Hairer et al., 1993, Sect. II.5] for further details). With this integrator we solve the set of differential equations composed of Eqs. (2) and (3) simultaneously. Here we use again the parameter to control the integrator’s overall behavior. For the DOP853 integrator the estimation of the local error and the step size control is based on embedded formulas of orders 5 and 3.
Before investigating how variational equations can be integrated efficiently, one should first clarify the meaning of the term ‘efficiency’. On the one hand, an efficient integration is one that is performed as fast as possible, and on the other hand, this computation should be also done as accurately as possible. Since accuracy always comes at the cost of intense computational efforts –meaning large CPU times– it tends to contradict the first mentioned aspect of efficiency. In addition, since we are especially interested in the computation of chaos indicators, an accurate computation should imply the correct distinction between regular and chaotic orbits of the studied dynamical system. Thus, in this work we consider the integration of variational equations to be efficient, when the computed GALIs are obtained with the least-possible CPU time requirements, achieving at the same time the correct characterization of orbits as regular or chaotic.
In the next section we present a thorough discussion of how variational equations can be integrated numerically using the methods described in Sect. 3. We explain possibilities of estimating the accuracy of such computations and discuss the obtained results with respect to our definition of efficiency. In Sect. 4.2 we apply this knowledge for a more global investigation of the properties of the FPU- lattice. As a final remark we note that all presented computations were performed on an Intel Xeon X3470 with 2.93 GHz, using extended precision (80 bit).
4.1 Efficient integration of variational equations
In most applications it is common to integrate the variational equations together with the equations of motion they are based on. Thus, an important prerequisite to obtain correct stability results is the correct computation of the dynamics of the system itself. For this reason we first discuss the properties of the integrators described in Sect. 3 when applied to some specific orbits of the FPU- lattice. In general, an accuracy estimate of a numerically obtained orbit of a dynamical system can be given via monitoring how well the system’s first integrals (e.g. the total energy, total linear momentum etc.) are conserved. For conservative Hamiltonian systems, like (1), this can be done easily by checking the conservation of the energy itself. The absolute value of the relative errors of the total energy , for different integrators and step sizes , for a orbit over time units are given in Table 4.1. For non-symplectic methods also the one-step accuracy is mentioned.
Firstly, we focus on the application of the TM method with the SABA and SBAB integrators. The corresponding results are given in the first 5 lines in Table 4.1. Comparing the energy conservation between SABA and SABAC one finds that the use of the corrector steps significantly improves for the same step size of . This result can of course be explained by the fact that SABA is an integration scheme of order 2, while SABAC is of order 4. Reducing the step size for SABAC to further improves the energy conservation as expected. We note that this reduction leads to a linear growth by a factor 5 of the required CPU time, as expected. Since SABAC shows the best performance, we use only this integrator for further investigations.
Comparing the results of SABAC with the ones obtained by the non-symplectic methods TIDES and DOP853 one finds that both non-symplectic methods need more CPU time in order to reach the same final value of . If one computes a mean step size for these algorithms, defined as the total integration time divided by the number of accepted steps , one finds that both integrators achieve this final energy error by a larger mean step size, compared to SABAC, which is due to the higher orders of these integrators. While the highest order for DOP853 is fixed to 8, TIDES uses automatic order selection, which explains the larger mean time step for the same value of one-step accuracy .
While monitoring energy conservation serves as a control parameter over the state vector of the system itself, it lacks information on how accurately the corresponding variational equations are solved. If the stability of certain initial conditions is known, one can use the theoretically predicted behaviors (5) of the GALI chaos indicator to estimate the reliability of the numerical computation. Therefore, in the last column of Table 4.1, we provide information on whether the integration was able to identify the orbit correctly as being a regular orbit lying on a 2-dimensional torus.
The time evolution of GALI, , for some of the runs of Table 4.1 is shown in Fig. 1. From Eq. (5) it is known that GALI should be constant for a orbit, while GALI and GALI should decrease proportionally to and respectively. A correct characterization is possible by using the TM method with SABA and a step-size as large as , although the corresponding energy error , is rather high (see first line in Table 4.1). For the non-symplectic methods it is found that leads to a similar error in energy conservation, but is not sufficient for the correct computation of the GALIs (see the first and third columns of Fig. 1). Decreasing the one-step accuracy improves the accuracy, but the integrations become less and less efficient compared to the TM integration as the required CPU time grows. We note that a correct dynamical characterization of the orbit is possible both for the DOP853 and TIDES for .
Using energy conservation as an indicator for the quality of an integration one could argue that the difference in CPU time between SABAC with and DOP853 with is not very significant. While this is true for the orbit, the difference becomes more pronounced, when the number of particles of the FPU- lattice is increased. This is evident in Fig. 2 where we plot, as function of , the ratio of the required CPU times between the TM method with SABAC and the TIDES (blue curves) and DOP853 (red curves) integrators for a regular orbit with initial condition , , for , 8, 12 and 20. If after time units an error in energy conservation of is sufficient, one could use either SABAC with , TIDES with (although for the corresponding GALIs do not exhibit the theoretically expected behavior for the whole time interval), or DOP853 with (see Table 4.1). The ratio of the CPU time of these runs is given in Fig. 2(a). In Fig. 2(b) we show a similar comparison when SABAC with , TIDES with , and DOP853 with are used, which yields at the end of the integration. In this case, the GALIs computed by these methods show the time evolution predicted in Eq. (5).
Figure 2 shows that the efficiency of the TM method improves with increasing when compared to the non-symplectic methods used in this study. While the ratio is 2 between the CPU times of SABAC and TIDES for in Fig. 2(b) it increases up to 7 when using particles in the FPU- lattice. The ratios are generally smaller when comparing the TM method with DOP853, but also here a growing trend can be observed. The results become even more pronounced when a larger error in energy conservation is acceptable, as it is shown in Fig. 2(a).
When analyzing dynamical systems one very often has to follow a trajectory for long time intervals to determine the behavior of a specific initial condition correctly. Especially for weakly chaotic orbits differences to regular orbits are shown only in the late evolution of chaos indicators. In such investigations the TM method also proves to be superior compared to other techniques. As an example we show in Fig. 3 the time evolution of GALIs for a regular orbit with initial condition , , , for the system (1) with , over time units. Further information on these runs are given in Table 4.1.
From Table 4.1 and Fig. 3 it is seen that for the non-symplectic integrators a one-step accuracy of is not sufficient for a correct identification of the regular orbit. Results obtained both by the TIDES and the DOP853 integrators show a deviation from the theoretically predicted behaviors (5) after (see the first and third columns of Fig. 3). To obtain the correct behavior of GALIs over the whole integration interval it is necessary to decrease the one-step accuracy to (see the second and fourth columns of Fig. 3). We note that this decrease in naturally results in a significant increase in CPU time. In contrast, the GALIs obtained via the TM method with require significantly less CPU time, and show the theoretically expected behaviors up to , although the relative energy error is .
While the error in energy conservation grows with time for non-symplectic integrators, it remains bounded for symplectic integrators, as can exemplarily be seen in Fig. 4. Thus, in case it is necessary to integrate beyond , one has to further decrease for the non-symplectic methods, in order to achieve the same final . This, of course, would result in a further increase in the ratio values of the CPU time required by the non-symplectic methods as compared to the TM technique.
4.2 Searching for motion on low-dimensional tori
One of the advantages of the GALI method is its capability to identify motion on low-dimensional tori. For a regular orbit the largest order of its GALIs that eventually remains constant determines the dimension of the torus on which the motion occurs (see Eq. (5)). This ability was verified in Skokos et al. ; Bountis et al. , where some particular orbits on low-dimensional tori were considered.
Since the TM method provides reliable evaluations of the GALIs even for relatively large integration steps, and therefore requires little CPU time, we applied this technique to perform a more global investigation of the FPU- lattice, aiming to trace the location of low-dimensional tori. In particular, we consider the Hamiltonian system (1) with , for which regular motion can occur on an -dimensional torus with , 3, 4. According to Eq. (5) the corresponding GALIs of order will be constant, while the remaining ones will tend to zero following particular power laws. Thus, in order to locate low-dimensional tori we compute the GALI, in the subspace of the system’s phase space, considering orbits with initial conditions , , while is computed to keep the total energy constant at .
Since the constant final values of GALI, , decrease with increasing order (see for example the GALIs with in Fig. 3), we chose to ‘normalize’ the values of GALI, of each individual initial condition, by dividing them by the largest GALI value, , obtained from all studied orbits. Thus, in Fig. 5 we colored each initial condition according to its ‘normalized GALI’ value
In each panel of Fig. 5, large values (colored in yellow or in light red) correspond to initial conditions whose GALI eventually stabilizes to constant, non-zero values. On the other hand, darker regions correspond to small values, which result from power law decays of GALIs.
Consequently, motion on 2-dimensional tori, which corresponds to large final GALI values and small final GALI and GALI values, should be located in areas of the phase space colored in yellow or light red in Fig. 5(a), and in black in Figs. 5(b) and (c). A region of the phase space with these characteristics is for example located in the upper border of the colored areas of Fig. 5. A particular initial condition with , in this region is denoted by a triangle in Fig. 5. This orbit is indeed a regular orbit as we see from the evolution of its GALIs shown in Fig. 6(a). In a similar way, a orbit should be located in regions colored in black only in the plot of Fig. 5(c). An orbit of this type is the one with , denoted by a square in Fig. 5, which actually evolves on a 3-dimensional torus, as only its GALI and GALI remain constant (Fig. 6(b)). Finally, the orbit with , (denoted by a circle in Fig. 5) inside a region of the phase space colored in yellow or light red in all panels of Fig. 5, is a regular orbit on a 4-dimensional torus and its GALI, GALI and GALI remain constant (Fig. 6(c)). It is worth noting, that chaotic motion would lead to very small GALI and values, since all GALIs would tend to zero exponentially, and consequently would correspond to regions colored in black in all panels of Fig. 5. Thus, chaotic motion can be easily distinguished from regular motion on low-dimensional tori.
From the results of Figs. 5 and 6 it becomes evident that the comparison of color plots of ‘normalized GALI’ values can facilitate the tracing of low-dimensional tori. The construction of such plots becomes a very demanding computational task, especially for high-dimensional systems. Thus, the application of the TM method for obtaining such results becomes imperative, since the required computations can be performed very efficiently by this method.
5 Summary and Conclusions
We compared different numerical techniques for the integration of variational equations of multi-dimensional Hamiltonian systems. In particular, we considered the TM method, which uses symplectic integrators for the realization of this task, as well as non-symplectic algorithms, like the general-purpose Runge-Kutta integrator DOP853, and the TIDES algorithm, which relies on Taylor series expansion techniques. These methods were applied to the D Hamiltonian system (1), the FPU- lattice, with varying from to .
We used the numerically obtained solutions of the variational equations for the computation of the GALI chaos indicators. The accurate reproduction of theoretically known behaviors of GALIs was used as a measure of reliability of the numerical techniques tested. In addition, the CPU time required by each method in order to achieve accurate results, was taken into account for the characterization of the efficiency of these algorithms.
The TM method exhibited the best numerical performance in all our simulations, both in accuracy and speed. More specifically, we found that the ratio of the CPU time required by the TIDES and DOP853 algorithms, with respect to the TM method, for correctly characterizing the nature of orbits, increased with increasing (Fig. 2). Thus, the TM method should be preferred over the other two techniques, especially for studies of multi-dimensional systems.
A feature of the TM method, which is of significant practical importance, is that it succeeds in finding the correct GALI behavior, and consequently determines the nature of orbits correctly, even when large integration steps are used, despite the fact that in these cases the energy accuracy is rather low. Therefore, the application of the TM method allows the efficient investigation of the dynamical properties of a large number of initial conditions in feasible CPU times. As an example, we showed in Sect. 4.2 how the TM method can exploit the properties of GALI to efficiently find the location of low dimensional tori in the phase space. Possible applications of this approach could be the tracing of quasiperiodic motion in multi-dimensional systems, when only a few degrees of freedom are excited.
The work of E. G. was financially supported by the DFG research unit FOR584. S. E. would like to acknowledge the support of the Austrian FWF project P20216. Ch. S. was partly supported by the European research project “Complex Matter”, funded by the GSRT of the Ministry Education of Greece under the ERA-Network Complexity Program.
- Abad et al.  Abad, A., Barrio, R., Blesa, F. & Rodriguez, M.  “TIDES: a Taylor series Integrator for Differential EquationS”, preprint 2011.
- Barrio  Barrio, R.  “Performance of the Taylor series method for ODEs/DAEs”, Appl. Math. Comput. 163, 525–545.
- Bountis et al.  Bountis, T., Manos, T. & Christodoulidi, H.  “Application of the GALI method to localization dynamics in nonlinear systems”, J. of Comp. Appl. Math. 227, 17–26.
- Chaos Focus Issue  Campbell D.K., Rosenau P. & Zaslavsky G.M. (eds.) “The Fermi-Pasta-Ulam Problem-The first 50 years”, Chaos, Focus Issue 15.
- Christodoulidi & Bountis  Christodoulidi, H. & Bountis, T.  “Low-dimensional quasiperiodic motion in Hamiltonian systems”, ROMAI Journal 2, 37–44.
- Delva  Delva, M.  “Integration of the elliptic restricted three-body problem with Lie series”, Celestial Mechanics, 34, 145–154.
- Eggl & Dvorak  Eggl, S. & Dvorak, R.  “An Introduction to Common Numerical Integration Codes Used in Dynamical Astronomy”, Lect. Notes in Physics, 790 431–480.
- Fermi et al.  Fermi, E., Pasta, J. & Ulam, S.  “Studies of Nonlinear Problems. I”, Los Alamos Report LA–1940.
- Ford  Ford, J.P.  “The Fermi-Pasta-Ulam Problem: Paradox turned Discovery”, Phys. Rep. 213, 271–310.
- Gröbner  Gröbner, W.  “Die Lie-Reihen und ihre Anwendungen”, Deutscher Verlag der Wissenschaften, 1967 - VI.
- Gerlach & Skokos  Gerlach, E. & Skokos, Ch.  “Comparing the efficiency of numerical techniques for the integration of variational equations” arXiv, 1008.1890.
- Hairer et al.  Hairer, E., Nørsett, S. P. & Wanner, G.  Solving Ordinary Differential Equations. Nonstiff Problems, 2nd Ed. (Springer Series in Comput. Math.).
- Hairer et al.  Hairer, E., Lubich, C. & Wanner, G.  Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, vol. 31, (Springer Series in Comput. Math.).
- Hanslmeier & Dvorak  Hanslmeier, A. & Dvorak, R.  “Numerical integration with Lie-series”, A&A 132, 203–207.
- Laskar & Robutel  Laskar, J. & Robutel, P.  “High order symplectic integrators for perturbed Hamiltonian systems”, Cel. Mech. Dyn. Astr. 80, 39–62.
- Manos & Ruffo  Manos, T., & Ruffo, S.  “Scaling with system size of the Lyapunov exponents for the Hamiltonian Mean Field model”, Transport Theory and Statistical Physics (in press, eprint arXiv:1006.5341).
- Manos & Athanassoula  Manos, T. & Athanassoula, E.  “Regular and chaotic orbits in barred galaxies - I. Applying the SALI/GALI method to explore their distribution in several models”, (submitted, eprint arXiv:1102.1157).
- Pál & Süli  Pál, A. & Süli, Á.  Solving linearized equations of the N-body problem using the Lie-integration method”, MNRAS 381, 1515–1526
- Skokos  Skokos, Ch.  “Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits”, J. Phys. A 34, 10029–10043.
- Skokos  Skokos, Ch.  “The Lyapunov Characteristic Exponents and their computation”, Lect. Notes in Physics, 790 63–135.
- Skokos et al.  Skokos, Ch., Antonopoulos, C., Bountis, T. & Vrahatis, M. N.  “How does the Smaller Alignment Index (SALI) distinguish order from chaos?”, Prog. Theor. Phys. Supp., 150, 439–443.
- Skokos et al.  Skokos, Ch., Antonopoulos, Ch., Bountis, T. & Vrahatis, M. N.  “Detecting order and chaos in Hamiltonian systems by the SALI method”, J. Phys. A 37, 6269–6284.
- Skokos et al.  Skokos, Ch., Bountis, T. & Antonopoulos, Ch.  “Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method”, Physica D 231, 30–54.
- Skokos et al.  Skokos, Ch., Bountis, T. C. & Antonopoulos, Ch.  “Detecting chaos, determining the dimensions of tori and predicting slow diffusion in Fermi–Pasta–-Ulam lattices by the Generalized Alignment Index method” Eur. Phys. J. Sp. T. 165, 5–14.
- Skokos & Gerlach  Skokos, Ch. & Gerlach, E.  “Numerical integration of variational equations”, Phys. Rev. E 82, 036704.
- Manos et al.  Manos, T., Skokos, Ch. & Antonopoulos, Ch.  “Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method”, (submitted, eprint arXiv:1103.0700).
- Suzuki  Suzuki, M.  “General theory of fractal path integrals with applications to many-body theories and statistical physics” J. Math. Phys. 32, 400–407.