The Lyapunov Characteristic Exponents and their computation
We present a survey of the theory of the Lyapunov Characteristic Exponents (LCEs) for dynamical systems, as well as of the numerical techniques developed for the computation of the maximal, of few and of all of them. After some historical notes on the first attempts for the numerical evaluation of LCEs, we discuss in detail the multiplicative ergodic theorem of Oseledec O_68 (), which provides the theoretical basis for the computation of the LCEs. Then, we analyze the algorithm for the computation of the maximal LCE, whose value has been extensively used as an indicator of chaos, and the algorithm of the so–called ‘standard method’, developed by Benettin et al. BGGS_80b (), for the computation of many LCEs. We also consider different discrete and continuous methods for computing the LCEs based on the QR or the singular value decomposition techniques. Although, we are mainly interested in finite–dimensional conservative systems, i. e. autonomous Hamiltonian systems and symplectic maps, we also briefly refer to the evaluation of LCEs of dissipative systems and time series. The relation of two chaos detection techniques, namely the fast Lyapunov indicator (FLI) and the generalized alignment index (GALI), to the computation of the LCEs is also discussed.
Keywords:Lyapunov exponents; Multiplicative ergodic theorem; Numerical techniques; Dynamical systems; Chaos; Variational equations; Tangent map; Chaos detection methods
For want of a nail the shoe was lost.
For want of a shoe the horse was lost.
For want of a horse the rider was lost.
For want of a rider the battle was lost.
For want of a battle the kingdom was lost.
And all for the want of a horseshoe nail.
For Want of a Nail (proverbial rhyme)
One of the basic information in understanding the behavior of a dynamical system is the knowledge of the spectrum of its Lyapunov Characteristic Exponents (LCEs). 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 nonstationary solutions of ordinary differential equations Lyapunov_1892 (), and has been widely employed in studying dynamical systems since then. The value of the maximal LCE (mLCE) is an indicator of the chaotic or regular nature of orbits, while the whole spectrum of LCEs is related to entropy (Kolmogorov–Sinai entropy) and dimension–like (Lyapunov dimension) quantities that characterize the underlying dynamics.
By dynamical system we refer to a physical and/or mathematical system consisting of a) a set of real state variables , whose current values define the state of the system, and b) a well–defined rule from which the evolution of the state with respect to an independent real variable (which is usually referred as the time ) can be derived. We refer to the number of state variables as the dimension of the system, and denote a state using the vector , or the matrix notation, where denotes the transpose matrix. A particular state corresponds to a point in an –dimensional space , the so–called phase space of the system, while a set of states parameterized by is referred as an orbit of the dynamical system.
Dynamical systems come in essentially two types:
Continuous dynamical systems described by differential equations of the form
with dot denoting derivative with respect to a continuous time and being a set of functions known as the vector field.
Discrete dynamical systems or maps, described by difference equations of the form
with being a set of functions and denoting the vector at a discrete time (integer).
Let us now define the term chaos. In the literature there are many definitions. A brief and concise presentation of them can be found for example in LC_04 (). We adopt here one of the most famous definitions of chaos due to Devaney (Devaney_1989, , p. 50), which is based on the topological approach of the problem.
Let be a set and a map on this set. We say that is chaotic on if
has sensitive dependence on initial conditions.
is topologically transitive.
periodic points are dense in .
Let us explain in more detail the hypothesis of this definition.
has sensitive dependence on initial conditions if there exists such that, for any and any neighborhood of , there exist and , such that , where denotes successive applications of .
Practically this definition implies that there exist points arbitrarily close to which eventually separate from by at least under iterations of . We point out that not all points near need eventually move away from under iteration, but there must be at least one such point in every neighborhood of .
is said to be topologically transitive if for any pair of open sets there exists such that .
This definitions implies the existence of points which eventually move under iteration from one arbitrarily small neighborhood to any other. Consequently, the dynamical system cannot be decomposed into two disjoint invariant open sets.
From Definition 1 we see that a chaotic system possesses three ingredients: a) unpredictability because of the sensitive dependence on initial conditions, b) indecomposability, because it cannot be decomposed into noninteracting subsystems due to topological transitivity, and c) an element of regularity because it has periodic points which are dense.
Usually, in physics and applied sciences, people focus on the first hypothesis of Definition 1 and use the notion of chaos in relation to the sensitive dependence on initial conditions. The most commonly employed method for distinguishing between regular and chaotic motion, which quantifies the sensitive dependence on initial conditions, is the evaluation of the mLCE . If the orbit is chaotic. This method was initially developed at the late 70’s based on theoretical results obtained at the end of the 60’s.
The concept of the LCEs has been widely presented in the literature from a practical point of view, i. e. the description of particular numerical algorithms for their computation F_84b (); ER_85 (); GPL_90 (); LichtenbergL_1992 (); DV_95 (). Of course, there also exist theoretical studies on the LCEs, which are mainly focused on the problem of their existence, starting with the pioneer work of Oseledec O_68 (). In that paper the Multiplicative Ergodic Theorem (MET), which provided the theoretical basis for the numerical computation of the LCEs, was stated and proved. The MET was the subject of several theoretical studies afterwards Ra_79 (); R_79b (); JPS_87 (); W_93 (). A combination of important theoretical and numerical results on LCEs can be found in the seminal papers of Benettin et al. BGGS_80a (); BGGS_80b (), written almost 30 years ago, where an explicit method for the computation of all LCEs was developed.
In the present report we focus our attention both on the theoretical framework of the LCEs, as well as on the numerical techniques developed for their computation. Our goal is to provide a survey of the basic results on these issues obtained over the last 40 years, after the work of Oseledec O_68 (). To this end, we present in detail the mathematical theory of the LCEs and discuss its significance without going through tedious mathematical proofs. In our approach, we prefer to present the definitions of various quantities and to state the basic theorems that guarantee the existence of the LCE, citing at the same time the papers where all the related mathematical proofs can be found. We also describe in detail the various numerical techniques developed for the evaluation of the maximal, of few or even of all LCEs, and explain their practical implementation. We do not restrict our presentation to the so–called standard method developed by Benettin et al. BGGS_80b (), as it is usually done in the literature (see e. g. F_84b (); ER_85 (); LichtenbergL_1992 ()), but we include in our study modern techniques for the computation of the LCEs like the discrete and continuous methods based on the singular value decomposition (SVD) and the QR decomposition procedures.
In our analysis we deal with finite–dimensional dynamical systems and in particular with autonomous Hamiltonian systems and symplectic maps defined on a compact manifold, meaning that we exclude cases with escapes in which the motion can go to infinity. We do not consider the rather exceptional cases of completely chaotic systems and of integrable ones, i. e. systems that can be solved explicitly to give their variables as single–valued functions of time, but we consider the most general case of ‘systems with divided phase space’ (Contopoulos_2002, , p. 19) for which regular111Regular orbits are often called ordered orbits (see e. g. (Contopoulos_2002, , p. 18)). (quasiperiodic) and chaotic orbits co–exist. In such systems one sees both regular and chaotic domains. But the regular domains contain a dense set of unstable periodic orbits, which are followed by small chaotic regions. On the other hand, the chaotic domains contain stable periodic orbits that are followed by small islands of stability. Thus, the regular and chaotic domains are intricately mixed. However, there are regions where order is predominant, and other regions where chaos is predominant.
Although in our report the theory of LCEs and the numerical techniques for their evaluation are presented mainly for conservative systems, i.e. system that preserve the phase space volume, these techniques are not valid only for such models. For completeness sake, we also briefly discuss at the end of the report the computation of LCEs for dissipative systems, for which the phase space volume decreases on average, and for time series.
We tried to make the paper self–consistent by including definitions of the used terminology and brief overviews of all the necessary mathematical notions. In addition, whenever it was considered necessary, some illustrative examples have been added to the text in order to clarify the practical implementation of the presented material. Our aim has been to make this review of use for both the novice and the more experienced practitioner interested in LCEs. To this end, the reader who is interested in reading up on detailed technicalities is provided with numerous signposts to the relevant literature.
Throughout the text bold lowercase letters denote vectors, while matrices are represented, in general, by capital bold letters. We also note that the most frequently used abbreviations in the text are: LCE(s), Lyapunov Characteristic Exponent(s); –LCE, Lyapunov Characteristic Exponent of order ; mLCE, maximal Lyapunov Characteristic Exponent; –mLCE, maximal Lyapunov Characteristic Exponent of order ; MET, multiplicative ergodic theorem; SVD, Singular Value Decomposition; PSS, Poincaré surface of section; FLI, fast Lyapunov indicator; GALI, generalized alignment index.
The paper is organized as follows:
In Section 2 we present the basic concepts of Hamiltonian systems and symplectic maps, emphasizing on the evolution of orbits, as well as of deviation vectors about them. In particular, we define the so–called variational equations for Hamiltonian systems and the tangent map for symplectic maps, which govern the time evolution of deviation vectors. We also provide some simple examples of dynamical systems and derive the corresponding set of variational equations and the corresponding tangent map.
Section 3 contains some historical notes on the first attempts for the application of the theoretical results of Oseledec O_68 () for the actual computation of the LCEs. We recall how the notion of exponential divergence of nearby orbits was eventually quantified by the computation of the mLCE, and we refer to the papers where the mLCE or the spectrum of LCEs were computed for the first time.
The basic theoretical results on the LCEs are presented in Section 4 following mainly the milestone papers of Oseledec O_68 () and Benettin et al. BGGS_80a (); BGGS_80b (). In Section 4.1 the basic definitions and theoretical results of LCEs of various orders are presented. The practical consequences of these results on the computation of the LCEs of order 1 and of order are discussed in Sections 4.2 and 4.3 respectively. Then, in Section 4.4 the MET of Oseledec O_68 () is stated in its various forms, while its consequences on the spectrum of LCEs for conservative dynamical systems are discussed in Section 4.5.
Section 5 is devoted to the computation of the mLCE , which is the oldest chaos indicator used in the literature. In Section 5.1 the method for the computation of the mLCE is discussed in great detail and the theoretical basis of its evaluation is explained. The corresponding algorithm is presented in Section 5.2, while the behavior of for regular and chaotic orbits is analyzed in Section 5.3.
In Section 6 the various methods for the computation of part or of the whole spectrum of LCEs are presented. In particular, in Section 6.1 the standard method developed in SN_79 (); BGGS_80b (), is presented in great detail, while the corresponding algorithm is given in Section 6.2. In Section 6.3 the connection of the standard method with the discrete QR decomposition technique is discussed and the corresponding QR algorithm is given, while Section 6.4 is devoted to the presentation of other techniques for computing few or all LCEs, which are based on the SVD and QR decomposition algorithms.
In Section 7 we briefly refer to various chaos detection techniques based on the analysis of deviation vectors, as well as to a second category of chaos indicators based on the analysis of the time series constructed by the coordinates of the orbit under consideration. The relation of two chaos indicators, namely the fast Lyapunov indicator (FLI) and the generalized alignment index (GALI), to the computation of the LCEs is also discussed.
Although the main topic of our presentation is the theory and the computation of the LCEs for conservative dynamical systems, in the last section of our report some complementary issues related to other types of dynamical systems are concisely presented. In particular, Section 8.1 is devoted to the computation of the LCEs for dissipative systems, while in Section 8.2 some basic features on the numerical computation of the LCEs from a time series are presented.
Finally, in the appendix 9 we present some basic elements of the exterior algebra theory in connection to the evaluation of wedge products, which are needed for the computation of the volume elements appearing in the definitions of the various LCEs.
2 Autonomous Hamiltonian systems and symplectic maps
In our study we consider two main types of conservative dynamical systems:
Continuous systems corresponding to an autonomous Hamiltonian system of degrees (D) of freedom 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 a vector
with , , . The time evolution of this orbit is governed by the Hamilton equations of motion, which in matrix form are given by
with , , and
Matrix has the following block form
with being the identity matrix and being the matrix with all its elements equal to zero. The solution of (2) is formally written with respect to the induced flow as
Symplectic maps of dimensions having the form
A symplectic map is an area–preserving map whose Jacobian matrix
The state of the system at the discrete time is given by
where , times.
2.1 Variational equations and tangent map
Let us now turn our attention to the (continuous or discrete) time evolution of deviation vectors from a given reference orbit of a dynamical system. These vectors evolve on the tangent space of . We denote by the linear mapping which maps the tangent space of at point onto the tangent space at point , and so we have with
where , are deviation vectors with respect to the given orbit at times and respectively.
We underline that equations (8) represent 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 . The solution of (8) can be written as
where is the so–called fundamental matrix of solutions of (8), satisfying the equation
In the case of the symplectic map (4) the evolution of a deviation vector , with respect to a reference orbit , is given by the corresponding tangent map
Thus, the evolution of the initial deviation vector is given by
with satisfying the relation
2.2 Simple examples of dynamical systems
As representative examples of dynamical systems we consider a) the well–known 2D Hénon–Heiles system HH_64 (), having the Hamiltonian function
with equations of motion
and b) the 4–dimensional (4d) symplectic map
with parameters , and . All variables are given (mod ), so , for . This map is a variant of Froeschlé’s 4d symplectic map F_72 () and its behavior has been studied in CG_88 (); SCP_97 (). It is easily seen that its Jacobian matrix satisfies equation (5).
2.3 Numerical integration of variational equations
When dealing with Hamiltonian systems the variational equations (8) have to be integrated simultaneously with the Hamilton equations of motion (2). Let us clarify the issue by looking to a specific example. The variational equations of the 2D Hamiltonian (14) are
This system of differential equations is linear with respect to , , , , but it cannot be integrated independently of system (15) since the and variables appear explicitly in it. Thus, if we want to follow the time evolution of an initial deviation vector with respect to a reference orbit with initial condition , we are obliged to integrate simultaneously the whole set of differential equations (15) and (17).
A numerical scheme for integrating the variational equations (8), which exploits their linearity and is particularly useful when we need to evolve more than one deviation vectors is the following. Solving the Hamilton equations of motion (2) by any numerical integration scheme we obtain the time evolution of the reference orbit (3). In practice this means that we know the values for , , where is the integration time step. Inserting this numerically known solution to the variational equations (8) we end up with a linear system of differential equations with constant coefficients for every time interval , which can be solved explicitly.
which is a linear system of differential equations with constant coefficients and thus, easily solved. In particular, equations (18) can by considered as the Hamilton equations of motion corresponding to the Hamiltonian function
with , , and
2.4 Tangent dynamics of symplectic maps
In the case of symplectic maps, the dynamics on the tangent space, which is described by the tangent map (11), cannot be considered separately from the phase space dynamics determined by the map (4) itself. This is because the tangent map depends explicitly on the reference orbit .
3 Historical introduction: The early days of LCEs
Prior to the discussion of the theory of the LCEs and the presentation of the various algorithms for their computation, it would be interesting to go back in time and see how the notion of LCEs, as well as the nowadays taken for granted techniques for evaluating them, were formed.
The LCEs are asymptotic measures characterizing the average rate of growth (or shrinking) of small perturbations to the orbits of a dynamical system, and their concept was introduced by Lyapunov Lyapunov_1892 (). Since then they have been extensively used for studying dynamical systems. As it has already been mentioned, one of the basic features of chaos is the sensitive dependence on initial conditions and the LCEs provide quantitative measures of response sensitivity of a dynamical system to small changes in initial conditions. 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. Therefore, the presence of positive LCEs is a signature of chaotic behavior. Usually the computation of only the mLCE is sufficient for determining the nature of an orbit, because guarantees that the orbit is chaotic.
Characterization of the chaoticity of an orbit in terms of the divergence of nearby orbits was introduced by Hénon and Heiles HH_64 () and further used by several authors (e. g. FL_70 (); F_70 (); F_72 (); SF_73 (); CF_75 (); CDGS_76 ()). In these studies two initial points were chosen very close to each other, having phase space distance of about , and were evolved in time. If the two initial points were located in a region of regular motion their distance increased approximately linearly with time, while if they were belonging to a chaotic region the distance exhibited an exponential increase in time (Figure 1).
Although the theory of LCEs was applied to characterize chaotic motion by Oseledec O_68 (), quite some time passed until the connection between LCEs and exponential divergence was made clear BGS_76 (); P_77 (). It is worth mentioning that Casartelli et al. CDGS_76 () defined a quantity, which they called ‘stochastic parameter’, in order to quantify the exponential divergence of nearby orbits, which was realized afterwards in BGS_76 () to be an estimator of the mLCE for .
So, the mLCE was estimated for the first time in BGS_76 (), as the limit for of an appropriate quantity , which was obtained from the evolution of the phase space distance of two initially close orbits. In this paper some nowadays well–established properties of were discussed, like for example the fact that tends to zero in the case of regular orbits following a power law , while it tends to nonzero values in the case of chaotic orbits (Figure 2).
The same algorithm was immediately applied for the computation of the mLCE of a dissipative system, namely the Lorenz system NS_77 ().
The next improvement of the computational algorithm for the evaluation of the mLCE was introduced in CGG_78 (), where the variational equations were used for the time evolution of deviation vectors instead of the previous approach of the simultaneous integration of two initially close orbits. This more direct approach constituted a significant improvement for the computation of the mLCE since it allowed the use of larger integration steps, diminishing the real computational time and also eliminated the problem of choosing a suitable initial distance between the nearby orbits.
In BGGS_78 () a theorem was formulated, which led directly to the development of a numerical technique for the computation of some or even of all LCEs, based on the time evolution of more than one deviation vectors, which are kept linearly independent through a Gram–Schmidt orthonormalization procedure (see also BG_79 ()). This method was explained in more detail in SN_79 (), where it was applied to the study of the Lorenz system and was also presented in BFS_79 (), where it was applied to the study of an D Hamiltonian system with varying from 2 to 10.
The theoretical framework, as well as the numerical method for the computation of the maximal, some or even all LCEs were given in the seminal papers of Benettin et al. BGGS_80a (); BGGS_80b (). In BGGS_80b () the complete set of LCEs was calculated for several different Hamiltonian systems, including four and six dimensional maps. In Figure 3
we show the results of BGGS_80b () concerning the 3D Hamiltonian system of CGG_78 (). The importance of the papers of Benettin et al. BGGS_80a (); BGGS_80b () is reflected by the fact that almost all methods for the computation of the LCEs are more or less based on them. Immediately the ideas presented in BGGS_80a (); BGGS_80b () were used for the computation of the LCEs for a variety of dynamical systems like infinite–dimensional systems described by delay differential equations F_82 (), dissipative systems ER_85 (), conservative systems related to Celestial Mechanics problems F_84 (); F_85 (), as well as for the determination of the LCEs from a time series WSSV_85 (); SS_85 ().
4 Lyapunov Characteristic Exponents: Theoretical treatment
In this section we define the LCEs of various orders presenting also the basic theorems which guarantee their existence and provide the theoretical background for their numerical evaluation. In our presentation we basically follow the fundamental papers of Oseledec O_68 () and of Benettin et al. BGGS_80a () where all the theoretical results of the current section are explicitly proved.
We consider a continuous or discrete dynamical system defined on a differentiable manifold . Let denote the state at time of the system which at time was at (see equations (3) and (6) for the continuous and discrete case respectively). For the action of over two successive time intervals and we have the following composition law
The tangent space at is mapped onto the tangent space at by the differential according to equation (7). The action of is given by equation (9) for continuous systems and by equation (12) for discrete ones. Thus, the action of on a particular initial deviation vector of the tangent space is given by the multiplication of matrix for continuous systems or for discrete systems with vector . From equations (9) and (12) we see that the action of over two successive time intervals and satisfies the composition law
This equation can be written in the form
where is the matrix corresponding to . We note that since we get and . A function satisfying relation (23) is called a multiplicative cocycle with respect to the dynamical system .
Let be a measure space with a normalized measure such that
for . Suppose also that a smooth Riemannian metric is defined on . We consider the multiplicative cocycle corresponding to and we are interested in its asymptotic behavior for . Since, as mentioned by Oseledec O_68 (), the case is analogous to the case , we restrict our treatment to the case , where time is increasing. In order to clarify what we are practically interested in let us consider a nonzero vector of the tangent space at . Then the quantity
is called the coefficient of expansion in the direction of . If
we say that exponential diverge occurs in the direction of . Of course the basic question we have to answer is whether the characteristic exponent (also called characteristic exponent of order 1)
We will answer this question in a more general framework without restricting ourselves to multiplicative cocycles. So, the results presented in the following Section 4.1, are valid for a general class of matrix functions, a subclass of which contains the multiplicative cocycles which are of more practical interest to us, since they describe the time evolution of deviation vectors for the dynamical systems we study.
4.1 Definitions and basic theorems
Let be an matrix function defined either on the whole real axis or on the set of integers, such that , for each time the value of function is a nonsingular matrix and the usual 2–norm of 222The 2–norm of an matrix A is induced by the 2–norm of vectors, i. e. the usual Euclidean norm , by and is equal to the largest eigenvalue of matrix .. In particular, we consider only matrices satisfying
with a suitable constant.
Considering a matrix function as above and a nonzero vector of the Euclidian space the quantity
is called the 1-dimensional Lyapunov Characteristic Exponent or the Lyapunov Characteristic Exponent of order 1 (1-LCE) of with respect to vector .
For simplicity we will usually refer to 1–LCEs as LCEs.
We note that the value of the norm does not influence the value of . For example, considering a vector , with a nonzero constant, instead of in Definition 4, we get the extra term (with denoting the absolute value) in equation (26) whose limiting value for is zero and thus does not change the value of . More importantly, the value of the LCE is independent of the norm appearing in equation (26). This can be easily seen as follows: Let us consider a second norm satisfying the inequality
for some positive real numbers , , and for all vectors . Such norms are called equivalent (see e.g. (HornJ_1985, , §5.4.7)). Then, by the above–mentioned argument it is easily seen that the use of norm in (26) leaves unchanged the value of . Since all norms of finite dimensional vector spaces are equivalent, we conclude that the LCEs do not depend on the chosen norm.
Let , be a set of linearly independent vectors in , be the subspace generated by all and be the volume of the –parallelogram having as edges the vectors . This volume is computed as the norm of the wedge product of these vectors (see Appendix 9 for the definition of the wedge product and the actual evaluation of the volume)
Let also be the volume of the initial –parallelogram defined by all , since is the identity matrix. Then the quantity
is called the coefficient of expansion in the direction of and it depends only on and not on the choice of the linearly independent set of vectors. Obviously for an 1–dimensional subspace the coefficient of expansion is . If the limit
exits it is called the characteristic exponent of order p in the direction of .
Considering the linearly independent set , and the corresponding subspace of as above, the p-dimensional Lyapunov Characteristic Exponent or the Lyapunov Characteristic Exponent of order p (p-LCE) of with respect to subspace is defined as
Similarly to the case of the 1–LCE, the value of the initial volume , as well as the used norm, do not influence the value of .
From (25) and the Hadamard inequality (see e. g. O_68 ()), according to which the Euclidean volume of a –parallelogram does not exceed the product of the lengths of its sides, we conclude that the LCEs of equations (26) and (27) are finite.
From the definition of the LCE it follows that
for any two vectors and with , while the Hadamard inequality implies that if , is a basis of then
where is the determinant of matrix .
It can be shown that for any the set of vectors is a vector subspace of and that the function with , takes at most different values, say
For the subspaces
with and if and only if for . So in descending order each LCE ‘lives’ in a space of dimensionality less than that of the preceding exponent. Such a structure of linear spaces with decreasing dimension, each containing the following one, is called a filtration.
A basis , of is called normal if attains a minimum at this basis. In other words, the basis , is a normal basis if
where , is any other basis of .
A normal basis , is not unique but the numbers depend only on and not on the particular normal basis and are called the LCEs of function . By a possible permutation of the vectors of a given normal basis we can always assume that .
Let , be a normal basis of and , with , , the LCEs of these vectors. Assume that value , appears exactly times among these numbers. Then is called the multiplicity of value and the collection is called the spectrum of LCEs.
In order to clarify the used notation we stress that , are the (possibly nondistinct) LCEs, satisfying , while , represent the (), different values the LCEs have, with .
The matrix function is called regular as if for each normal basis