Numerical integration of variational equations

Numerical integration of variational equations

Ch. Skokos Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, D-01187 Dresden, Germany    E. Gerlach Lohrmann Observatory, Technical University Dresden, D-01062 Dresden, Germany
August 17, 2019

We present and compare different numerical schemes for the integration of the variational equations of autonomous Hamiltonian systems whose kinetic energy is quadratic in the generalized momenta and whose potential is a function of the generalized positions. We apply these techniques to Hamiltonian systems of various degrees of freedom, and investigate their efficiency in accurately reproducing well-known properties of chaos indicators like the Lyapunov Characteristic Exponents (LCEs) and the Generalized Alignment Indices (GALIs). We find that the best numerical performance is exhibited by the ‘tangent map (TM) method’, a scheme based on symplectic integration techniques which proves to be optimal in speed and accuracy. According to this method, a symplectic integrator is used to approximate the solution of the Hamilton’s equations of motion by the repeated action of a symplectic map , while the corresponding tangent map , is used for the integration of the variational equations. A simple and systematic technique to construct is also presented.

45.10.-b, 05.45.-a, 02.60.Cb

I Introduction

Numerical integration is very often the only available tool for investigating the properties of nonlinear dynamical systems. Different numerical techniques Hairer_etal_93 ; Hairer_etal_02 have been developed over the years which permit the fast and accurate time evolution of orbits in such systems.

Of particular interest are the so-called ‘symplectic integrators’ which are numerical methods specifically aimed at advancing in time the solution of Hamiltonian systems with the aid of symplectic maps (see for example (Hairer_etal_02, , Chapt. VI), SI_Ham and references therein). Another challenging numerical task in conservative Hamiltonian systems is to discriminate between order and chaos. This distinction is a delicate issue because regular and chaotic orbits are distributed throughout phase space in very complicated ways. In order to address the problem several methods have been developed, which can be divided into two major categories: the ones based on the study of the evolution of deviation vectors from a given orbit, like the computation of the maximal Lyapunov Characteristic Exponent (mLCE) S10 , and those relying on the analysis of the particular orbit itself, like the frequency map analysis of Laskar FMA .

Other chaos detection methods, belonging to the same category with the evaluation of the mLCE, are the fast Lyapunov indicator (FLI) FLI and its variants Barrio , the smaller alignment index (SALI) S01b_SABV03_SABV04 and its generalization, the so-called generalized alignment index (GALI) SBA07 ; SBA08 , and the mean exponential growth of nearby orbits (MEGNO) MEGNO . The computation of these indicators require the numerical integration of the so-called variational equations, which govern the time evolution of deviation vectors.

The scope of this paper is to present, analyze and compare different numerical methods for the integration of the variational equations. In our study we consider methods based on symplectic and non-symplectic integration techniques. The integration of the variational equations by non-symplectic methods is straightforward since one simply has to integrate these equations simultaneously with the equations of motion. This approach requires in general, more CPU time than schemes based on symplectic integration techniques for the same order of accuracy and integration time step. For this reason we focus our attention on methods based on symplectic schemes, explaining in detail their theoretical foundation and applying them to Hamiltonian systems of different numbers of degrees of freedom.

The numerical solution of the variational equations obtained by the various integration schemes studied are used for the computation of the spectrum of the Lyapunov Characteristic Exponent (LCEs) and the GALIs. We chose to compute these two chaos indicators among the indices based on the evolution of deviation vectors, because the computation of the mLCE is the elder and most commonly employed chaos detection technique, while the computation of the whole spectrum of LCEs and GALIs requires the evolution of more than one deviation vector and thus is strongly influenced by inaccuracies of the integration procedure. We investigate the numerical efficiency of the different integration methods by comparing the CPU times they require for the computation of the LCEs and the GALIs, as well as their accuracy in reproducing well-known properties of these chaos indicators. In particular, we check whether the set of computed LCEs consists of pairs of values having opposite signs, and if the time evolution of GALIs follows specific theoretically predicted laws.

The paper is organized as follows: after introducing the concept of variational equations in the next section, we describe in sections III and IV the LCEs and the GALIs respectively, which are the two chaos indicators we use in our study. Then, in section V we give the basic properties of symplectic integrators. Section VI is devoted to the detailed description of several numerical schemes for the integration of the variational equations of Hamiltonian systems. Applications of these schemes to regular and chaotic orbits of systems with two or more degrees of freedom are presented in section VII, where also the efficiency of each technique is discussed. Finally, in section VIII, we summarize the results and present our conclusions, while in the appendix the explicit expressions of the various integration methods for the Hénon-Heiles system are given.

Ii The variational equations

Let us consider an autonomous Hamiltonian system of degrees of freedom (D) having a Hamiltonian function


where and , are the generalized coordinates and conjugate momenta respectively. An orbit in the -dimensional phase space of this system is defined by the vector


with , , . The time evolution of this orbit is governed by the Hamilton’s equations of motion, which in matrix form are given by


with , , and

with denoting the transpose matrix. Matrix has the following block form

with being the identity matrix and being the matrix with all its elements equal to zero.

An initial deviation vector from an orbit evolves in the tangent space of according to the so-called variational equations


with being the Hessian matrix of Hamiltonian (1) calculated on the reference orbit , i. e.

Equations (4) are a set of linear differential equations with respect to , having time dependent coefficients since matrix depends on the particular reference orbit, which is a function of time .

In the present paper we consider autonomous Hamiltonians of the form


with being the potential function. The Hamilton’s equations of motion (3) become


while the variational equations (4) of this system take the form


with , , and


Thus, the tangent dynamics of Hamiltonian (5) is represented by the time dependent Hamiltonian function


which we call the ‘tangent dynamics Hamiltonian’ (TDH), and whose equations of motion are exactly the variational equations (7).

Iii The Lyapunov Characteristic Exponents

The LCEs are asymptotic measures characterizing the average rate of growth (or shrinking) of small perturbations to the solutions of a dynamical system. Their concept was introduced by Lyapunov when studying the stability of non-stationary solutions of ordinary differential equations L1892 , and has been widely employed in studying dynamical systems since then. A detailed review of the theory of the LCEs, as well as of the numerical techniques developed for their computation can be found in S10 .

The theory of LCEs was applied to characterize chaotic orbits by Oseledec O68 , while the connection between LCEs and exponential divergence of nearby orbits was given in BGS76 ; P77 . For a chaotic orbit at least one LCE is positive, implying exponential divergence of nearby orbits, while in the case of regular orbits all LCEs are zero or negative. Therefore, the computation of the mLCE is sufficient for determining the nature of an orbit, because guarantees that the orbit is chaotic.

The mLCE is computed as the limit for of the quantity


often called finite time mLCE, where , are deviation vectors from a given orbit, at times and respectively, and denotes the norm of a vector. So, we have


If the energy surface defined by (1) is compact, it has been shown that this limit is finite, independent of the choice of the metric for the phase space and converges to for almost all initial vectors O68 ; BGGS80a ; BGGS80b . tends to zero in the case of regular orbits following a power law BGS76


while it tends to nonzero values in the case of chaotic orbits.

An D Hamiltonian system has (possibly non-distinct) LCEs, which are ordered as . In BGGS78 a theorem was formulated, which led directly to the development of a numerical technique for the computation of all LCEs, based on the time evolution of many deviation vectors, kept linearly independent through a Gram-Schmidt orthonormalization procedure. The theoretical framework, as well as the corresponding numerical method for the computation of all LCEs (usually called the ‘standard method’), were given in BGGS80a ; BGGS80b . According to this method all other LCEs , etc., apart from the mLCE obtained from (11), are computed as the limits for of some appropriate quantities , etc., which are called the finite time LCEs (see BGGS80b ; S10 for more details). We note that throughout the present paper, whenever we need to compute the values of the LCEs, we apply the discrete QR-decomposition technique (NumRec, , Sect. 2.10), which is a variation of the standard method (see Sect. 6.3 of S10 for more details).

It has been shown in BGGS80a that in the case of an autonomous Hamiltonian flow, the set of LCEs consists of pairs of values having opposite signs


In addition, since the Hamiltonian function is an integral of motion, at least two LCEs vanish, i. e.


while the presence of any additional independent integral of motion leads to the vanishing of another pair of LCEs.

Iv The Generalized Alignment Index

The GALI is an efficient chaos detection technique introduced in SBA07 as a generalization of a similar indicator called the smaller alignment index (SALI) S01b_SABV03_SABV04 . The method has been applied successfully for the discrimination between regular and chaotic motion, as well as for the detection of regular motion on low dimensional tori to different dynamical systems SBA08 ; GALI_appl .

The GALI of order () is determined through the evolution of initially linearly independent deviation vectors , . The time evolution of each deviation vector is governed by the variational equations (7). Each evolved deviation vector is normalized from time to time, having its norm equal to 1, in order to avoid overflow problems, but its direction is left intact. Then, according to SBA07 , is defined to be the volume of the -parallelogram having as edges the unitary deviation vectors , . This volume is equal to the norm of the wedge product of these vectors, and is given by


From this definition it is evident that if at least two of the deviation vectors become linearly dependent, the wedge product in (15) becomes zero and the vanishes.

Expanding the wedge product (15) into a sum of determinants and studying the asymptotic behavior of those who vary the slowest in time, it is possible to show analytically the following SBA07 : in the case of a chaotic orbit all deviation vectors tend to become linearly dependent, aligning in the direction defined by the mLCE and tends to zero exponentially following the law


where are approximations of the first largest Lyapunov exponents. On the other hand, in the case of regular motion on an -dimensional torus, all deviation vectors tend to fall on the -dimensional tangent space of this torus. Thus, if we start with general deviation vectors they will remain linearly independent on the -dimensional tangent space of the torus, since there is no particular reason for them to become aligned. As a consequence is different from zero and remains practically constant for . On the other hand, tends to zero for , since some deviation vectors will eventually become linearly dependent, following a particular power law which depends on the dimensionality of the torus and the number of deviation vectors. The behavior of for regular orbits lying on -dimensional tori is given by


If the regular orbit lies on a low dimensional torus, i. e. an -dimensional torus with then remains practically constant and different from zero for and tends to zero for following particular power laws (see SBA08 for more details).

In order to compute the value of we consider the matrix having as columns the coordinates of the unitary deviation vectors , , , with respect to the usual orthonormal basis , , …, of the -dimensional tangent space and perform the Singular Value Decomposition (SVD) of this matrix. Then, as it was shown in SBA08 , is equal to the product of the singular values , of matrix , i. e.


V Symplectic integrators

Let us discuss in some detail how we can integrate the equations of motion (3) of a general Hamiltonian (1) by a symplectic integration scheme, focusing our attention on a particular family of integrators presented in LR01 . Defining the Poisson bracket of functions , by comment1 :


the Hamilton’s equations of motion (3) take the form


where is the differential operator defined by . The solution of Eq. (20), for initial conditions , is formally written as


Let us assume that the Hamiltonian function can be split into two integrable parts as . A symplectic scheme for integrating equations (20) from time to time consists of approximating, in a symplectic way, the operator by an integrator of steps involving products of operators and , , which are exact integrations over times and of the integrable Hamiltonians and . The constants , , which in general can be positive or negative, are chosen to increase the order of the remainder of this approximation. So , are actually symplectic maps acting on the coordinate vector . Therefore the integration of equations (20) over one time step , which evolves the initial coordinate vector to its final state , is represented by the action on of a symplectic map produced by the composition of products of and . In this context several symplectic integrators of different orders have been developed by various researchers SI_methods ; ML95 .

In LR01 the families of SBAB (and SABA) symplectic integrators, which involve only forward (positive) integration steps were introduced. These integrators were adapted for the integration of perturbed Hamiltonians of the form , where both and are integrable and is a small parameter. A particular integrator SBAB (), or SABA (), involves steps, i. e. applications of products of and , and is of order with respect to the integration step . This means that by using these integrators, we are actually approximating the dynamical behavior of the real Hamiltonian by a Hamiltonian , i. e. we introduce an error term of the order .

The accuracy of the () integrator can be improved when the commutator term comment2 leads to an integrable system, as in the common situation of being quadratic in momenta and depending only on positions . In this case, two corrector terms of small backward (negative) steps can be added to the integrator


A similar expression is valid also for . The value of constant is chosen in order to eliminate the dependence of the remainder which becomes of order . The SBAB (SABA) integrators have already proved to be very efficient for the numerical study of different dynamical systems LR01 ; si_appl1 ; si_appl2 . We note that several authors have used commutators for improving the efficiency of symplectic integrators (e. g. Chin97 ; OMF_23 ).

Setting we can apply the SBAB (SABA) integration schemes for the integration of Hamiltonian (5), since this Hamiltonian can be written as , with


being both integrable. The maps , , which propagate the set of initial conditions at time , to their final values at time , for the Hamiltonian functions and (23) are




respectively. For Hamiltonian (5) the corrector term is given by


which is a function of only the coordinates and thus easily integrated as


In Appendix A.1 we give the explicit formulas of equations (24), (25) and (27) for the Hénon-Heiles system (54).

Vi Numerical integration of variational equations

In this section we present several numerical schemes for the integration of the variational equations, considering both non-symplectic techniques and methods based on symplectic integrators. The latter schemes are quite general and any symplectic integrator can be used for their implementation. In our study we consider an efficient fourth order symplectic integrator, the Chin97 ; LR01 , which has an extra degree of complexity with respect to integrators composed of products of maps , and , since it requires the application of the corrector term (26).

vi.1 Non-symplectic schemes

In order to follow the evolution of a deviation vector, the variational equations (7) have to be integrated simultaneously with the Hamilton’s equations of motion (6), since matrix depends on the particular reference orbit , which is a solution of equations (6). Any non-symplectic numerical integration algorithm can be used for the integration of the whole set of equations (6) and (7).

In our study we use the DOP853 integration method which has been proven to be very efficient. The DOP853 integrator foot1 is an explicit non-symplectic Runge-Kutta integration scheme of order 8, based on the method of Dormand and Price (see (Hairer_etal_93, , Sect. II.5)). Two free parameters, and , are used to control the numerical performance of the method. The first one defines the time span between two successive outputs of the computed solution. After each step of length the values of LCEs (GALIs) are computed and the deviation vectors are orthonormalized (normalized). For the duration of each step , the integrator adjusts its own internal time step, so that the local one-step error is kept smaller than the user-defined threshold value . For DOP853 the estimation of this local error and the step size control is based on embedded formulas of orders 5 and 3.

vi.2 Integration of the tangent dynamics Hamiltonian

Another approach to compute the evolution of deviation vectors is to initially integrate the Hamilton’s equations of motion (6), in order to obtain the time evolution of the reference orbit , and then to use this numerically known solution for solving the equations of motion of the TDH (9), which are actually the variational equations (7).

In practice one numerically solves the Hamilton’s equations of motion (6) by any (symplectic or non-symplectic) integration scheme to obtain the values at , , where is the integration time step of these orbits. Of course, the accuracy of the particular numerical scheme used for the construction of the time series will affect the quality of the numerical solution of the variational equations, regardless of the numerical scheme used for solving them. Having computed the values different methods can be applied for approximating the solution of the variational equations, which will be discussed in the following sections.

vi.2.1 TDH with piecewise constant coefficients

One method is to approximate the actual time dependent TDH (9) by a Hamiltonian with piecewise constant coefficients. This means to assume that the coefficients of (9) are constants equal to for the time interval . These constants are determined by the values of the orbit’s coordinates and are known, since we know the time series . Thus, for each time interval we end up with a quadratic form Hamiltonian function , whose equations of motion form a linear system of differential equations with constant coefficients.

The Hamiltonian can be integrated by any symplectic or non-symplectic integration scheme, or can be explicitly solved by performing a canonical transformation to new variables , , so that the transformed Hamiltonian becomes a sum of uncoupled 1D Hamiltonians, whose equations of motion can be integrated immediately. To this end, let be the eigenvalues and , the unitary eigenvectors of the constant matrix . Then matrix , having as columns the eigenvectors , defines a canonical change of variables , , which gives the diagonal form


The equations of motion of are then easily solved.

In our study we use the same symplectic integrator () both for obtaining the time series and for integrating the quadratic form Hamiltonian in the time interval . We name this approach the TDHcc method (cc: constant coefficients). An alternative approach is to compute the exact solution of the equations of motion of (whose piecewise constant coefficients are obtained by the symplectic integration of the orbit using the scheme) by transforming it to a system of uncoupled harmonic oscillators through the canonical transformation induced by matrix . This approach is called the TDHes method (es: exact solution).

In general, the transformation matrix is determined for each time interval by solving numerically the eigenvalue problem


a procedure which could become computationally very time consuming, especially for systems with many degrees of freedom. On the other hand, in some simple low dimensional cases, like for example the Hénon-Heiles system (54), the transformation matrix can be determined analytically (see Appendix A.2.1).

vi.2.2 Integration of the TDH in an extended phase space

Instead of approximating (9) by a quadratic form having constant coefficients for each time interval , we can explicitly treat as a time dependent Hamiltonian. This time dependency is due to the fact that the coefficients of are functions of the orbit’s coordinates . Like in the previous approach, we consider the time series to be known from the numerical integration of the Hamilton’s equations (6).

The D time dependent Hamiltonian can be transformed to a time independent Hamiltonian with an extra degree of freedom by considering the time as an additional coordinate (see for example (LichtenbergL_1992, , Sect. 1.2b)). For this purpose, we add to the Hamilton’s equations of motion of the equations


Then we set and as an additional coordinate and momentum respectively, i. e. , , and define the new Hamiltonian


where and are respectively the new coordinates and momenta. The flow in the -dimensional extended phase space of the D Hamiltonian is parameterized by a ‘new’ time such that , which does not appear explicitly in the functional form of (31). The set of equations (7) and (30) are the Hamilton’s equations of motion of .

The dynamics of the D TDH (9) is equivalent to that of the D Hamiltonian


This Hamiltonian can be easily integrated by any symplectic integration scheme, since it can be split into two integrable parts


The maps , , which propagate the set of initial conditions at time , to their final values at time are


The corrector term of the SBAB and SABA integration schemes


is a function of only the coordinates and thus easily integrated


The explicit expressions of these maps for the Hénon-Heiles system (54) are given in Appendix A.2.2.

From equations (34), (35) and (37) we see that time is changed only by the act of operator . On the other hand, operators and require the knowledge of positions at specific times for the evaluation of the partial derivatives of and . We also note that for all these operators the last equation for can be neglected, since the knowledge of its value does not influence the evolution of the other quantities, and consequently the solution of the variational equations (7).

Since the coordinates of the orbit are known only at specific times , , one is restricted to use integration schemes that require the knowledge of at exactly these times. Such a scheme is, for example, the integrator


(which is practically the well-known Störmer/Verlet or leap-frog method) with . The right operator which acts first, requires the knowledge of , while the left operator needs the values of , because the time value has changed from to by . Note that the integrator


requires the knowledge of for the application of . This second order integration scheme could be used with , leading in general to a less accurate algorithm compared to (38), which is also a second order integrator but uses a smaller time step . For is in general more efficient to apply the integration scheme


which was initially derived in ML95 . This integrator needs the known values , and .

The above integration schemes can also be combined with a corrector step, since (37) does not change the time values, and acts before and after the main body of the integrator (see equation (22)), when has values for which we know the coordinates . We refer to this technique as the TDHeps method (eps: extended phase space). For the numerical applications of the TDHeps method (presented in Sect. VII) we use the fourth order integrator both for the integration of the variational equations and for the computation of the orbit.

Higher order SBAB or SABA integrators cannot be used in this framework, because they require the knowledge of at non equidistant time values, different from . In order to apply such schemes one could initially compute the solution of equations (6) also at these specific times (e. g. by interpolation), but this would lead to a cumbersome, complex, time consuming, and consequently inefficient scheme.

vi.3 The tangent map (TM) method

The set of equations (6) and (7) can be considered as a unified set of differential equations


where is a vector formed by the phase space vector and the deviation vector , and is the differential operator of the whole system. In analogy to equation (21), the solution of system (41) for an initial condition can be formally written as . We describe now how symplectic integrators can be used to obtain this solution.

First of all, let us note that equations (41) cannot be considered as the Hamilton’s equations of motion of some generalized Hamiltonian function. If such a Hamiltonian existed, and could be split into two integrable parts, any symplectic integrator could be used for finding the solution of system (41). Since this is not the case, we follow a different approach to achieve this goal. In section V the integration of the equations of motion of Hamiltonian (5) over one integration time step was split into steps over appropriate time intervals , , where the dynamics was determined either by Hamiltonian or (23). During these intermediate steps the tangent dynamics of the system is governed by the variational equations


for , and by


for . Therefore, for each intermediate step of the symplectic integration scheme the dynamics of the phase and the tangent space is governed by the set of equations




for Hamiltonians and (23) respectively, with and being the corresponding differential operators.

These sets of equations are immediately solved, leading to maps


Obviously the first two equations of maps , are exactly maps (24) and (25), respectively.

Thus, any symplectic integration scheme used to solve the Hamilton’s equations of motion (6), which involves the successive application of maps (24), (25), can also be used for the simultaneous integration of the variational equations (7), i. e. for solving the set of equations (41), by replacing maps , with maps (46) and (47) respectively. This statement is a specific application of a more general result which is stated for example in LR01 : Symplectic integration schemes can be applied to first order differential systems that can be written in the form , where , and are differential operators for which the two systems and are integrable. The system of differential equations (41) belongs to this category since it can be split into the integrable systems (44) and (45).

Let us discuss this splitting in more detail. The system (41) can be written as


with , , and being a vector with coordinates


Then the dynamics of any general variable is given by