Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
WWaldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
Franz Waldner a* , Rainer Klages ba Physics Institute, University of Zurich ,Winterthurerstr. 190, CH-8057 Switzerland b Queen Mary University of London, School of Mathematical Sciences, Mile End Road, London E1 4NS, UK _____________________________________________________________________________________A B S T R A C T______________________________________________________________________________________
The stability analysis introduced by Lyapunov and extended by Oseledec is an excellent tool to describe the character of nonlinear n- dimensional flows by n global exponents if these flows are stable in time. However, there are two main shortcomings: (a) The local exponents fail to indicate the origin of instability where trajectories start to diverge. Instead, their time evolution contains a much stronger chaos than the trajectories, which is only eliminated by integrating over a long time. Therefore, shorter time intervals cannot be characterized correctly, which would be essential to analyse changes of chaotic character as in transients. (b) Moreover, although Oseledec uses an n dimensional sphere around a point x to be transformed into an n dimensional ellipse in first order, this local ellipse has yet not been evaluated. The aim of this contribution is to eliminate these two shortcomings. Problem (a) disappears if the Oseledec method is replaced by a frame with a ‘constraint’ as performed by Rateitschak and Klages (RK) [Phys. Rev. E 65 036209 (2002)]. The reasons why this method is better will be illustrated by comparing different systems. In order to analyze shorter time intervals, integrals between consecutive Poincaré points will be evaluated. The local problems (b) will be solved analytically by introducing the symmetric ‘Jacobian deformation ellipsoid’ and its orthogonal submatrix, which enable to search in the full phase space for extreme local separation exponents. These are close to the RK exponents but need no time integration of the RK frame. Finally, four sets of local exponents are compared: Oseledec frame, RK frame, Jacobian deformation ellipsoid and its orthogonal submatrix. PACS:
Keywords:
Lyapunov exponents, stability analysis,, dynamical instability, Jacobian deformation ellipsoid.________________________________________________________________________________ * Corresponding author: Fax +41 44 635 57 04, Tel. +41 44 923 46 27
E-mail address : [email protected]; [email protected]; [email protected]
1. Introduction
Lyapunov exponents [1] are a well-established tool [2-9] to analyse the type of chaos displayed by a trajectory x ( t ) as a solution of a nonlinear n- dimensional equation in the phase space R { x },d x / d t = X ( x ) (1) . Their evaluation is based on proofs by Oseledec [2] that a vector u to any arbitrarily chosen neighbouring point of x (except the instantaneous direction along d x as discussed later) will rotate to the direction of integrated extreme expansion. After an initial period, this local direction obtained through integrated rotation is an inherent feature of each point on the trajectory. Moreover, this local uniqueness applies to the set of all n directions of a frame orthogonalised subsequently after each time step. The integration for infinite time t → ∞ of the corresponding local exponents results in a set of n Lyapunov exponents assessing the character of chaos in terms of dynamical instability of the given dynamical system. The fact that an initial frame rotates after a short time to a locally unique orientation has fascinated many authors. This fascination includes also the corresponding local exponents which has, probably, prevented one to ask: Is the local exponent leading to the largest Lyapunov exponent already locally assessing the strongest expansion, i.e. the strongest divergence of neighbouring trajectories, thus identifying the places where trajectories ‘break away’, the origin of instability and chaos? Even a quick look at this aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited exponent for the Lorenz model reveals that this is not the case resulting in misleading results for shorter time intervals. The origin of this discrepancy is inherent to the idea of Lyapunov and Oseledec: They considered a sphere of neighbouring points. Each vector u not orthogonal to the direction d x of the trajectory contains a component parallel to d x . This component measures only the expansion along the trajectory, hence the acceleration, which is not connected to instability. Only the component orthogonal to d x analyzes the change of the distance to neighbouring trajectories and can indicate where trajectories start to diverge. This feature is obscured by the admixing of the acceleration. Since in many cases the average acceleration is zero, the integration for long times cancels this contribution, thus only the value for divergence remains.These shortcomings can easily be avoided by using a method introduced by Rateitschak and Klages (RK) [10]: Select as an initial direction u the direction of d x , thus u is parallel to the trajectory. Certainly, this direction will remain parallel to the trajectory for all times. The integral of the corresponding exponent will only describe the average acceleration which for many cases is zero, i.e., if the flow remains bounded and does neither shrink nor blow up. This exponent is not related to the problem of instability. However, all remaining directions of the initial frame, subsequently orthogonalised and their rotation integrated, are describing the divergence of neighbouring trajectories. Their corresponding exponents are the important numbers for analysing the nature of chaos. After an initial period much shorter than for the Oseledec method these directions are unique to the point x ( t ). The first in the orthogonalizing procedure of the corresponding local exponents is the largest one quantifying where trajectories diverge. Moreover, integrals for shorter time intervals analyse the chaotic character of the trajectories thus discriminating between intervals of strong and weak or no chaos. This problem of finding coordinate systems in which trivial eigendirections are eliminated has already been addressed by Eckhardt and Wintgen in 1991 [11]. However, they focussed on periodic orbits in conservative two degree of freedom Hamiltonian systems for which they eliminated the two trivial neutral directions along the orbit and perpendicular to it on the energy shell by also requiring certain smoothness properties. Their Hamiltonian method was reviewed and further amended by Gaspard for calculating local Lyapunov exponents and local stretching rates [12]. Other simple numerical methods for computing local Lyapunov exponents and stretching rates have been proposed and tested by Dellago and Hoover [13] and by Rateitschak and Klages (RK) [10]. Both Oseledec and RK methods will be illustrated for the Lorenz [14] and the Rössler [15] model by comparing them with each other and by also applying them to shorter time intervals for transient chaos. It is one of the main points of this contribution to show that the method using a set orthogonal to the flow is more adequate to describe the chaotic behaviour than the Oseledec method using general directions, thus explaining why Rateitschak and Klages [10] found much improved results introducing this method for complex chaotic behaviour. Jacobian deformation ellipsoid
In order to understand the definition and use of what later on we call ‘Jacobian deformation ellipsoid’ we first briefly review the origins of Lyapunov instability analysis:(i) The basic idea of Lyapunov exponents is to follow the evolution of points close to the points x ( t ) on a trajectory. First, these neighbouring points are chosen on an n- dimensional sphere around the starting point x ( t o ). Then this sphere is continuously deformed during the time evolution. After a sufficiently long integration time τ this deformed object is analysed. From the largest expansion direction the largest Lyapunov exponent is evaluated. Starting with the direction of this largest expansion, further orthogonal directions are used to measure their expansions, which complete the set of n Lyapunov exponents. (ii) In principle, the evolution of the deformation should be numerically computed for an infinite number of points on the initial sphere, a tremendous task even for large computers. (iii) However, the works of Lyapunov [1] and Oseledec [2] propose a well-established, much less elaborate method: Arbitrarily defined at the start, choose an orthogonal frame of only n directions, follow their evolution - orthogonalised always in the same order by the Gram-Schmidt orthonormalisation method before each integration step - and measure their expansion, thus disregarding their rotation. The Lyapunov exponents are then obtained as the averages of the logarithms of these expansions.It would be interesting to test to which extent (iii) corresponds to (i) in case (ii) could be treated with less numerical effort. It is the aim of this section to propose a simpler method. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
Before using the computer for doing (ii), let’s go back to the 19 th century of Jacobi. The main idea was then: Think of a reduction of the problem and retain only the essentials. Jacobi [16] used only the first derivative to define his matrix J at a point x, J ( x ) = d X / d x (2) . He realised that the sum of the diagonal elements of J measures the local rate of change of the volume of a sphere around x in the limit of infinitesimally small radius. However, nothing has been said then about the deformation of this sphere. Although finally global quantities should be evaluated, we will first consider ‘instantaneous local’ quantities defined by the deformation between a fixed time t and t+ d t, starting as a sphere at t .Applying the 19 th century method to the deformed sphere, what is the first approximation? This has been shown already by Farmer et al. [4] with a picture of an ellipse in two dimensions. Green and Kim [7] describe a general ellipsoid for n dimensions, which is continuously deformed in time but remains always an ellipsoid. The problem (i) is then solved by using the n principal axes and their corresponding expansions for the final ellipsoid. Obviously, these principal axes are orthogonal in accord with (iii), since the ellipsoid has inversion symmetry. An n- dimensional ellipsoid is described by a symmetric n x n matrix E , with orthogonal eigenvectors as principal axes and real eigenvalues, and is determined by n ( n +1)/2 parameters. At this point it seems worthwhile to note that E yields the radius for the infinite number of possible directions as diagonal elements ( E U ) ii of E U , after transforming to a new coordinate system with the unitary matrix U and its transpose U T as E U = U T E U . The matrix E has the same size as J . However, J is in principle not symmetric causing complex eigenvalues and non-orthogonal complex eigenvectors. Furthermore, its diagonal elements are connected with the rate of change of the radius of a sphere, not with the radius of the ellipsoid described by E . Therefore, why not try to symmetrise J while keeping the diagonal elements by defining a symmetric matrix S = ½ ( J + J T ) with orthogonal eigenvectors and real eigenvalues in order to describe the rates of change of the radii when the sphere is transformed in first order into an ellipsoid? This paper aims to convince the reader that S is exactly describing in first order the ‘instantaneous local rates of exponential stretching ratio’ – in short ‘local exponents’ - for all possible directions, although only n numbers are involved in J . This symmetric result has almost been found by Greene and Kim [7] in their eqs. (26)-(28). They showed that their instantaneous local expansion exponents λ U along any orthogonal set of directions U are given by the diagonal elements K ii of K = U T J U (note that in general J and K are not symmetric). Simply reduce this equation for only one direction u to the scalar product λ u =( u , J u ). This form can be derived in a direct geometrical way providing probably one of the simplest approaches to define local Lyapunov exponents. The derivation is so short that it will be sketched here in the introduction, as follows: The local exponent λ for a dynamical system given by eq.(1) can be defined as λ = (1/Δ t ) ln ( r ), with r = ( d +Δ d )/ d = 1 + Δ d / d being the stretching ratio for the length d of a vector u pointing from a point x ( t ) towards a neighbouring point. Using the Jacobian matrix J , the new vector after time Δ t is found in first order to be u (t+Δ t ) = u ( t )+ J u ( t ) Δ t. This new vector has rotated and changed its length to d +Δ d. It is now important to eliminate the rotation by projecting the change Δ u = J u Δ t onto u by using the scalar product Δ d = ( u , J u ) Δ t / │u │ . valid for Δ t →0. Putting this into r = 1 + Δ d / d and expanding ln(1+ α ) ≈ α for│α│<<1 results in λ = (1/Δ t ) ln ( r ) ≈ (1/Δ t ) ( u , J u ) Δ t / │u │ . The simple result λ =( u , J u ) is then found if the vector u has been normalised to unity. Note that the `logarithmisation’ is no longer visible in the form for λ after the above approximation has been applied. Furthermore, note that this result has the correct dimension of reciprocal time as needed for a Lyapunov exponent defined above as ‘the rate of exponential stretching ratio’. Interestingly, the form λ u =( u , J u ) shows that u and – u give the same value implying inversion symmetry. In more detail, the elements J ik of the matrix J appear in pairs ( J ik + J ki ) in the scalar product ( u , J u ). Hence J can be replaced by the symmetric matrix S = ½ ( J + J T ), with J T denoting the transposed matrix. This gives the same value for λ u = ( u , S u ) as λ u = ( u , J u )=( u , J T u ). Going now back to the above form K = U T J U of Green and Kim [7] and replacing there J by S , the diagonal elements T ii of the new matrix T = U T S U are equal to K ii . In contrast, this new T is a symmetric matrix corresponding to the symmetrised K with T = ½ ( K + K T ). Furthermore, the transformation of J with K = U T J U could be interpreted as transforming J into the new reference system U and the result might be denoted by J U . Similarly, T = U T S U transforms S into S U = U T S U with its diagonal elements ( S U ) ii as instantaneous local exponents in the U directions. Thus S is the generating form containing all deformations of the ellipsoid and will be called ‘Jacobian deformation ellipsoid’ aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited (not to be confused with ‘deformed ellipsoid E ' which describes directly the radii, not the rates of their changes). Note that using the symmetry argument is a consequence of the reduction of the deformation to the lowest order. In reality, the deformations are more complicated than being captured in lowest order. As a by-product, the largest eigenvalue of the main principal axis of S gives the extreme rate of expansion of the sphere, its eigenvectors rarely being parallel to any of the frames in Section 1.1 The above first step defines the instantaneous local deformation. Before doing the second step, let us recall here three relevant features of the ‘Jacobian deformation ellipsoid’ S : First, the matrix S = ½ ( J + J T ) is obviously symmetric by definition. Second, the `logarithmisation’ is not visible explicitly since the approximation ln(1+ α ) ≈ α for small α has been used; nevertheless, S describes the exponential rates of instantaneous local stretching. Third, S yields this rate for all directions around the point x , although only n numbers of the elements of the local Jacobian matrix J ( x ) are needed to define the matrix S .The second step now consists of evaluating global quantities both for schemes (i) and (iii) above in order to test them in the long-time limit. In both cases, global exponents are defined as final deformations after integrating all deformations during a time interval τ=t final – t initial for τ→ ∞.Before going into further detail, two important questions have to be answered: Have `deformations of previous deformations' to be evaluated? And have rotations of the ellipsoids to be incorporated? The answers to both questions are illustrated in figure 1 of Benettin et al. [3]: After each time step the new deformation does not account for the previous deformation, that is, the deformation of the previous deformation is of higher order and can be neglected. Thus instantaneous local quantities can be averaged if integration is replaced by small but finite time steps as in numerical work. As far as rotation is concerned, each new start in Benettin's figure 1 is from a rotated direction. Since S does not contain any information about the rotation, it has to be incorporated separately, thus at each time step the transformed S U ( t ) = U T ( t ) S ( t ) U ( t ) has to be evaluated. Note that for U ( t ) any of the two orthonormalised frames of section 1.1 could be used. Furthermore, note that S U ( t ) has usually non-zero off-diagonal elements implying that the diagonal elements are no eigenvalues. Now, on the basis of the new Jacobian deformation ellipsoid S , differences between (i) and (iii) will be described. The evaluation according to (i) has to be done in two steps: First, average all local matrices S U ( t ) along a typical trajectory. Obviously, each of the n elements of S U ( t ) has to be averaged separately. Hence this averaging results in a final symmetric matrix S final . In a second step, only for this final matrix S final the eigenvalues and the eigenvectors have to be evaluated corresponding to global Lyapunov exponents and Lyapunov directions, respectively. The recipe for performing (iii) is much simpler: The averages of the n diagonal elements ( S U ) ii ( t ) of the local matrix S U ( t ) along a typical trajectory are the global Lyapunov exponents. Again obviously, each of the n diagonal elements of ( S U ) ii ( t ) has to be averaged separately. Comparing both methods, the results of (iii) are already incorporated in the final matrix S final as diagonal elements. Hence, to test the equivalence between (i) and (iii) it is sufficient to compare the values of the diagonal elements of S final with its eigenvalues. This test has been performed numerically for the ‘stable’ Lorenz chaos. There are small but distinct deviations for the free running rotation of Oseledec. The test is successful for the constraint rotation of RK.So far the literature has mainly focused on local and global Lyapunov exponents. Both quantities have been defined above within a new approach, and tests of this concept will be described later on in this paper. So what else? Previous concepts of Lyapunov instability are only appropriate for ‘stable’ chaos, where a trajectory does not change its character in time. They are not suitable to analyse transients, crises, or continuous changes of parameters in time in equations of motion. However, an adequate method is easy to find: Instead of only considering ‘global quantities’ defined for an infinite interval of time τ , try a series of successive finite time intervals τ n , each starting at t n and ending at t n+1 . It is essential to define the successive times t n such that the resulting data correspond to the character of the trajectory. To give an example, for a ‘stable’ Lorenz chaos time intervals τ n between successive Poincaré points are convenient, which have been chosen such that each interval describes a loop on the left or on the right hand side of the strange attractor. Again both methods (i) and (iii) are compared for different frames of rotation. Only the RK frame produces rather small but distinct deviations of the values for (i) and (iii). These discrepancies are caused by residual acceleration parts parallel to the flow dx . It is easy to eliminate all these directions by forming the n -1 dimensional subspace S ┴(2… n ) orthogonal to the flow of each instantaneous local S U ( t ) with U as the RK frame, a straightforward evaluation. For (i) the eigenvalues of the averaged subspace are then compared with (iii) as the averages of the exponents orthogonal to d x of RK aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited in the same time interval. Indeed, both values are equal within numerical accuracy for all time intervals.The method (iii) based on a series of finite intervals is then successfully applied to a transient Lorenz chaos. It distinguishes well between an initial period of weak chaos, an intermediate period of strong chaos, and a final spiralling onto a stable fixed point. Also here the RK frame is superior to the Oseledec frame. As a by-product, each instantaneous local subspace S ┴(2… n ) - to be calculated at any point x of the phase space R { x } – has n -1 eigenvalues. The largest eigenvalue within this subspace accounts for the extreme rate of divergence between neighbouring trajectories, not obscured by partial acceleration. This yields a novel indicator for extreme local divergence, which can be used to find ‘hot spots’ of maximum local dynamical instability in the whole phase space.These applications demonstrate that the new ‘Jacobian deformation ellipsoid’ is a powerful tool, which furthermore opens a pedestrian approach to defining both local and global Lyapunov exponents. The novelty of this article is that it compares four types of local directions and exponents with each other. The first two are strictly local to x . They use only the knowledge of the local Jacobian J and the value of d x , there is no need to evaluate any trajectory by integration: 1. The instantaneous local extreme expansion exponent of a local sphere and its corresponding direction are found as the largest eigenvalue and eigenvector of S . Note that this yields the maximal possible value of the exponent for all four different methods, and the direction of this maximal deformation is rarely parallel to the direction obtained by the other three methods. 2. The extreme exponent for the instantaneous local divergence of neighbouring trajectories measured orthogonal to d x and its direction are found as the largest eigenvalue and eigenvector of the subset S ┴ .Note that this exponent is the maximal possible value of the largest exponent for divergence in what we will call the W frame of RK as described in method 4 below.The other two types are evaluated by method (iii). and need integrations both of a trajectory and the directions of the frames.3. The standard method of Oseledec (O) with an integrated free running frame called V frame produces n instantaneous local exponents and directions.4. The method of RK using an integrated W frame, where the first direction is always constrained to be parallel to d x thus assessing the instantaneous local acceleration parallel to d x (needs no integration). The remaining instantaneous local exponents are orthogonal to d x describing the divergence of neighbouring trajectories, with their corresponding directions found by integration of the W frame. A test of how well the exponents of methods 3. and 4. are correlated with the exponents of methods 1. and 2.
A heuristic test uses Poincaré points and the corresponding local exponents. The distribution of these local exponents is tested as a function of the distance between the corresponding Poincaré points. The Lorenz model reveals that the O method has a much larger chaotic character than the RK method. This confirms again that the RK method is more efficient than the O method. Only the RK method is adequate to evaluate meaningful Lyapunov exponents, which furthermore indicate where the instabilities start locally.
2. Treating a local point x of the phase space Local expansion of a vector u to a neighbouring point of x ( t ) during a time interval Δ t will be described geometrically. In first order, using the Jacobian matrix J u ( t +Δ t ) = u ( t ) + J ( t ) u ( t ) Δ t = u ( t ) + Δ u (3) the new vector u ( t +Δ t ) will have rotated and changed its length from d to d +Δ d . If only the change Δ d of the length is considered, the rotation can be eliminated by projecting Δ u onto u by using the scalar product Δ d ≈ ( u , J u ) Δ t / │ u │ valid for Δ t→
0. With the ratio r of stretching r = ( d +Δ d ) / d = 1 + Δ d / d the local exponent λ u for the direction of u is found according to the form, see Greene and Kim [7], λ u = (1/Δ t ) ln ( r ) (4) giving λ u = (1/Δ t ) ln ( 1 + Δ d / d ) = (1/Δ t ) ln [ 1 + ( u , J u ) Δ t / │ u │ ]. Using ln(1+α) ≈ α for │α│ << 1, the result is the scalar product λ u = ( u , J u ) , (5) if the vector u has been normalized to unity. Greene and Kim [7] published in their eqs. (26-28) the same result in form of the diagonal elements K ii aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited λ i = K ii K = V T J V (6) for an orthonormalised reference frame V and the Jacobian matrix J . Note that in the scalar product ( u, J u ) the values J ik and J ki appear in pair sums ( J ik + J ki ). This symmetry is enforced by the special structure of the scalar product ( u , J u ): on both sides there occurs the same vector u . Therefore, the same result is obtained for a symmetrised Jacobian matrix S constructed by adding the transpose J T . S = ½ ( J + J T ) (7) λ u = ( u , S u ) (8) In more detail, the Jacobi matrix J can be decomposed into J = S + D with the symmetric matrix S defined by eq. (7) and the anti-symmetric matrix D with D ik = - D ki defined by the difference D = ½ ( J – J T ). The scalar product of eq. (5) thus results in ( u , J u )=( u , S u ) + ( u , D u ). The second term ( u , D u ) is zero, because the same vector on both sides of the scalar product enforces terms which add up to zero due to the anti-symmetry of D . Therefore, the reduction of ( u , J u ) to ( u , S u ) as performed in eq. (8) is justified. Furthermore, note that ( u , D u )=0 implies that the vector D u is orthogonal to u describing only a rotation of u . At first glance, for capturing the rotation of u in eq. (3) it is thus tempting to reduce the term J u therein to D u . However, the vector S u can have any direction thus incorporating both elongation and rotation before it is projected onto u by the scalar product ( u , S u ) to find only elongation. Hence, for describing the full rotation of u it is essential to keep both terms ( S+D ) u in eq. (3).The symmetrised Jacobian matrix S will pave the way to identify the Jacobian deformation ellipsoid as demonstrated in the next section. A local n dimensional sphere around a point x of eq. (1) will be rotated and deformed after a time interval Δ t into a complicated geometrical object. Only in linear approximation its form can be described by a symmetric n dimensional ellipsoid with orthogonal principal axes of the principal deformation exponents. Here, a straightforward simple way to find this ellipsoid will be described: The expansion exponent λ u in any direction u is found by the scalar product ( u , S u ), see eq. (8). Defining an orthonormal set U and its transpose U T , the transformed matrix S U = U T SU can be evaluated. It is now important to note that the scalar product of eq. (8) implies that only the diagonal elements ( S U ) ii are equal to the local expansion exponents λ i corresponding to the n orthogonal normalised directions u i in U . Note further that the off-diagonal elements ( S U ) ik with i≠k can have any non-zero value. Let us first consider the special case that all off-diagonal elements ( S U ) ik with i≠k are zero. In this case, the diagonal elements ( S U ) ii = λ i could be eigenvalues of S . Since S is by definition symmetric, its eigenvalues α i are real and the corresponding eigenvectors a i are both real and orthogonal. It is then obvious from eq. (8) that the eigenvalues α i can indeed be local expansion exponents. By expanding the arbitrary vector u in eq.(8) into the basis of eigenvectors of S , in complete analogy to the expectation value problem of a quantum mechanical operator [17], one concludes that the eigenvalues α i are the extreme local expansion exponents, denoted by λ i(e) . This implies that the eigenvectors a i correspond to the extreme vectors u i(e) of U (e) , which are orthogonal by definition. They transform the matrix S into its diagonal eigenvalue form S diag =( U (e) ) T S U (e) . The vectors u i (e) thus define the principal axes of the ellipsoid into which the local sphere is deformed.Note that in the general case of arbitrary directions of U the corresponding local expansions exponents λ i =( S U ) ii are not eigenvalues of S . In summary, according to eq. (8) the matrix S can be considered as the generator of deformations in all possible local directions including the principal axes which yield the extreme values [18]. It is therefore justified to call S the n -dimensional ‘Jacobian deformation ellipsoid’. Note that S does not describe the deformed sphere. It accounts for the rate of expansion (positive values) or contraction (negative values) in any direction; hence, this ‘deformation ellipsoid’ can have values of both signs. It will be convenient to order the exponents α according to their values α k > α k+ , thus the largest first, and the eigenvectors a k , written as columns in the matrix A x accordingly, its components expressed in the frame of R { x }.The Jacobian deformation ellipsoid can be written as a symmetric n dimensional tensor T , its explicit form depending on the frame of reference. The simplest way is in the local frame of principal axes, the matrix T diag with only the exponents α k in the diagonal. Its trace as the sum of the exponents is clearly the rate of change of the volume of the sphere around x . The more convenient form T x would be expressed in the reference frame of R {x}, to be found by back transformation, resulting, obviously, in the form of S expressed usually in the reference frame of R { x }, with the trace unchanged as the rate of change of the volume, aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited T x = A x T diag A xT = S x . (9) An alternative frame of reference will be described in the next section. ┴(2,…n) for divergence Although the Jacobian deformation ellipsoid contains all local information about neighbouring points, it is essential to find the principal exponents and their directions orthogonal to the flow. Only these directions indicate true divergence, because they are not disturbed by partial acceleration.First the reason will be explained why these directions are important. Then an arbitrary local frame of reference is introduced serving to find the local frame of extreme divergence. The ellipsoid is transformed into this frame. After its truncation to the n -1 dimensional subspace orthogonal to the flow d x , the new extreme exponents and principal directions are found by solving the n -1 dimensional eigenvalue problem. Although these procedures are straightforward, the matrix operations will be described in more detail in order to facilitate their implementation in programs for Lyapunov exponents. The Jacobian deformation tensor will be described now by S x ( x ) at the point x , the subscript denotes that the matrix is written in the coordinate system R { x }. S x ( x ) describes the rate of stretching along all possible directions u to a neighbouring point of any of the points x of the phase space without the need to execute any integration of eq. (1). Whereas most of these points belong to a neighbouring trajectory, there is one point exactly on the trajectory through x . This point is found for Δ t→ X of eq. (1) which describes the flow through x . Defining a flow unit vector f ║ = X /│ X │, the corresponding expanding exponent φ ║ φ ║ = ( f ║ , S x f ║ ) (10) does not describe any divergence of a neighbouring trajectory. Instead, it is the relative acceleration (d v /d t )/ v related only to the change of velocity v along the trajectory through x [7]. Hence, any vector u can be decomposed into a vector sum u = c ║ f ║ + c ┴ g ┴ , where g ┴ is the appropriate unit vector in the space orthogonal to f ║ . Therefore, u has its expansion exponent λ u composed of the exponent φ ║ related only to acceleration and λ ┴g describing only the divergence of the neighbouring trajectory line in the g ┴ direction, λ u = c ║ ( f ║ , S f ║ ) + c ┴ ( g ┴ , S g ┴ ) = c ║ φ ║ + c ┴ λ ┴g (11) Hence, if the sole interest is to find locally the largest divergence of neighbouring trajectory lines, and not to be disturbed by interference with acceleration, the n -1 dimensional subspace orthogonal to the flow direction f ║ has to be constructed. This will be performed using an arbitrary local orthonormal set of reference F with the first column as the unit vector f ║ . Then, construct the remaining n -1 vectors f k by permutation of the components of f ║ . Finally, use a Gram-Schmidt procedure to make the f k orthogonal as columns of the matrix F x ( x ) expressed in the frame of reference R { x } (denoted by subscript x ) at the local point x .The Jacobian deformation ellipsoid S x is transformed into the F frame by S F = F xT S x F x (12) The first row and column of S F describe the expansion along the flow. The remaining n -1 dimensional reduced square submatix S ┴ (2,…, n ) contains all expansions in the orthogonal subspace. The reduced square submatix S ┴ (2,… n ) is symmetric and has n -1 eigenvalues β ┴k and eigenvectors b ┴k as the principal perpendicular local exponents and directions, respectively, again both to be ordered according to their values β ┴k > β ┴k+ , the largest first, and B is the n -1 dimensional matrix containing the ordered vectors b ┴ as columns.It is worthwhile to construct a local reference frame { h } as the matrix H with the first vector as flow direction f ║ and the remaining directions as local perpendicular extreme expansion directions b ┴ k . In the F frame, the first row and column are zero except H F =1. The remaining submatrix is filled with the matrix B . This frame H F can be transformed into the reference frame of R { x } by the arbitrary matrix F x H x = F x H F (13) aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited Although the exponents β ┴k were evaluated already by solving the eigenvalue problem of the reduced submatrix S ┴ (2,… n ) , a general relation as a numerical control could be written using the frame matrix H x in a matrix product with the diagonal elements (…) ii giving φ ║ as the first, and β ┴k=1,…,n-1 as the remaining numbers φ ║ = ( H xT S x H x ) β ┴k =(1,…,n-1) = ( H xT S x H x ) (k+1),(k+1) (14) Before a specific trajectory is evaluated, it seems worthwhile to explore the phase space by producing an n dimensional map of the principal exponent of local divergence β ┴ k=1 in order to find 'hot' regions of large divergence or 'cool' regions of missing divergence. However, applying this procedure to the Lorenz attractor the interpretation is not trivial; there are strong ‘hot’ and very ‘cool’ regions well outside the strange attractor. Examples will be given later. Moreover, changing the parameters in eq. (1) could result in a very different behaviour, more, stronger, less or no 'hot spots' in the whole phase space. To test this would be very elaborate by evaluating various trajectories. The procedure described here is much faster.
3. Exponents following a trajectory
A specific trajectory x ( t ) is found after choosing a starting point x = x ( t =0) by integrating eq. (1) starting at x . In the spirit of Lyapunov, n local exponents could be evaluated if a specific orthogonal set of directions is defined for each x ( t ). The average of these local exponents will for t →∞ lead to the Lyapunov exponents Λ k characterizing the type of trajectory.Therefore, the problem of finding these local directions is essential. First, the well established method of Oseledec (O) without a constraint will be shortly described. Then the new method of Rateitschak and Klages (RK) with their constraint will be introduced. A comparison and a test will how later on that only the second method should be used, implying a fundamental change for the description of instability. According to Oseledec [2], any arbitrarily chosen orthogonal frame V at the origin x rotated according to eq. (1), then orthogonalised and normalised after each time step Δ t , will become a unique local frame V [ x ( t )] for each point x ( t ) of the trajectory after an initial transient time τ transV . The rotation can be evaluated by applying eq. (1) to neighbouring points or by using eq. (3) as proposed by Greene and Kim [7]. The idea is that the first direction, which is never adjusted by the orthogonalising process, will turn to the direction of strongest divergence from the trajectory, independently from its starting direction in V . The instantaneous local expansion exponents λ Vk can be found by eq. (5) or (8), and their averages will be the global Lyapunov exponents Λ Vk for t→∞ . The first exponent Λ Vk= will be the largest.At this point it is interesting to note that and Meier [19] found in their numerical analysis of the angles of the V frame that the direction corresponding to the smallest (‘most negative’) local exponent is always nearly orthogonal to the flow. Small deviations are probably due to finite step integration and numerical limitations. Rateitschak and Klages (RK) [10] introduced a new concept for a local frame W . The rotation can be made by eq. (3), and the first vector w k =1 = f ║ is always set parallel to X /│ X │. The remaining directions are then orthogonalised always in the same order. Also here, after a transient time τ transW a unique frame W [ x ( t )] for each point of the trajectory x ( t ) will be defined. The local expansion exponents λ Wk can be found by eq. (5) or (8) by using W instead of V , and their time averages will be the global exponents Λ Wk for t→∞ . Now, the first Λ Wk= = Λ ║ will not be the largest. It has a different function: it measures the mean acceleration of the flow, which is zero for many chaotic models. The remaining Λ Wk> = Λ ┴ Wk all describe only divergence of neighbouring trajectory lines, hence the second Λ Wk= will be the largest. A numerical test with the Lorenz model [14] confirmed that the third direction of the V frame is not only always nearly orthogonal to the flow [20], but also always nearly parallel to the third direction of the W frame, with deviations of the order of the deviations within the V frame. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited Fig. 1.
Lorenz chaos: Difference of the angles between the maximal expansion direction in a local frame flipped at time t =6.8 and the respective direction in the unflipped frame vs. time t . Above: free frame V of O [3]; Below: constraint frame W of RK [10] exhibiting a much shorter recovery time. How fast does an arbitrarily chosen frame approach the locally unique frame directions, described by a transient time τ transient ? This question is analysed in the Lorenz system by measuring the recovery time τ recover of a suitably rotated frame. The tests starts with a frame V at time t . At a time t >> τ transient a new frame V rot ( t ) is constructed, which is rotated by 90 degrees with respect to V ( t ). Both frames are integrated in time, and the difference of the angles between the maximal expansion direction of the flipped V rot ( t ) and the unflipped V ( t ) are plotted in Fig. 1 as a function of time t . Clearly, the recovery time τ recover is much longer for O (top) than for RK (bottom). For RK, the direction of the flow is obviously not rotated at t Therefore, only n -1 directions have to readjust. These recovery times τ recover are an indication for the transient times τ trans after t .In addition, it seems worthwhile to note that the exponent of the acceleration of the flow has zero recovery time, since its direction X /│ X │ is the local value of eq. (1) for each point x of the phase space. The corresponding values of the Lyapunov exponents Λ for t→∞ are numerically tested for three and four dimensions using the Lorenz and the Rössler model, respectively. Fig. 2 (Lorenz), and Fig. 3 (Rössler) display the results as functions of integration time t . At the bottom, the Lorenz x component or the Rössler x component is shown. On top, the free V frame of O is used, below the constraint W frame of RK. Fig. 4 displays in more detail the values between t = 700 and 800 of the first three integrated exponents. Fig. 2.
Lorenz chaos: Bottom: x -component vs. time. Top: The first two exponents integrated vs. time following the free V frame according to O. [3]. Centre: the same integrated exponents vs. time following the constraint W frame according to RK [10]. Fig. 3.
Rössler hyper-chaos: Bottom: x -component vs. time. Top: The first three exponents integrated vs. time following the free V frame according to O. [3]. Centre: the same integrated exponents vs. time following the constraint W frame according to RK [10]. The first two exponents are positive indicating hyper-chaotic behaviour; the third is zero for zero mean acceleration. The final magnitudes are nearly equal and within the fluctuations of both frames, see Fig.4. The fluctuations are more pronounced for the W frame (below), since its recovery time is much shorter, therefore the strong peaks of x (see Fig. 3 bottom) are not so well integrated out as for the V frame (above). Note that the third exponent of the free V frame (above) is slightly below zero although it should be zero according to the theory. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
10 / 24
Fig. 4.
Rössler hyper-chaos: Blow-up of the final time sequence between t =700 and 800 of Fig. 3 . These interesting results confirm that the largest Oseledec exponent Λ Vk= describes the divergence, although its local values are mixed with acceleration. For the Lorenz model the second Λ Vk= and for the Rössler model the third exponent Λ Vk= tend to a zero value, thus describing the mean acceleration Λ Wk= of the RK method by averaging out local acceleration and divergence.Hence, for infinite time and ‘stable’ chaos both methods are equivalent since the contribution of acceleration is cancelled by the integration. However, as shown later, for transients and chaotic regimes changing their character only the RK method is adequate.
4. Comparing local directions and exponents as functions of time t The following local sets of directions have been studied, each with n corresponding local exponents:(1.1.) (A) The principal axes a in A of the ‘full’ Jacobian deformation ellipsoid S , with local exponents α. (1.2.) (B) The ‘constrained’ flow direction f ║ and the extreme divergence directions b ┴ of the ‘reduced’ Jacobian submatrix S ┴ (2,… n ) (in contrast to the ‘full’ S ), both described by H , with local exponents φ ║ and β ┴ . (2.1.) (V). The ‘free’ Oseledec (O) unique local frame V approximated after a transient time τ transV , with local exponents λ V .(2.2.) (W) The ‘constrained’ Rateitschak and Klages (RK) unique local frame W approximated after a transient time τ transW , with local exponents λ ║ W = φ ║ and λ ┴W . The connections between these four local schemes are illustrated in Table 1. free ‘full’ constrained ║and ‘reduced’┴ to d x ║ f ║ strictly local for allpoints x in phase spaceno integration 1.1. A a i α i B b ║ = f ║ , b ┴ k β ║ = φ ║ and β ┴ k local extreme matrix, unit vectorslocal exponentslocal at x ( t ), but after integration alongtrajectory { x ( t k )} k= n V v i λ Vi W w ║ = f ║ , w ┴ k λ ║ W = φ ║ and λ ┴Wk local ‘unique’ matrix, unit vectorslocal exponents Table 1.
Illustration of the relations between the elements in the four local frames. Note that the parallel exponents of 1.2. and 2.2. are equal and both need no integration.
Fig. 5.
Lorenz chaos; Bottom,
Left : y component vs. time t. Right : z component vs. y component showing the two loops. Top: The angle between the principal axis a of the local Jacobian deformation ellipsoid S with the largest exponent and the first axis v of the V frame. Left : vs. time t . Right : vs. y component. Centre: The angle between the axis b ┴ of the local extreme exponent orthogonal to the flow (main axis of the reduced ellipsoid S ┴ (2,…, n ) and the first axis w ┴1 orthogonal to the flow of the W frame. Left : vs. time t . Right : vs. y component.Note that the free V frame (top) has a much more complex relation to the trajectory than the constraint W frame (centre). aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
11 / 24
Obviously, the extreme local exponents of ‘full’ (α) and ‘reduced’ (β) frames are the limiting values of the corresponding exponents of ‘free full’ ( λ V ) and ‘constrained’ ( λ W ) frames, respectively. Deviations between these exponents will be large when the corresponding directions of the different frames are very different. Therefore, angles between the main directions in the different frames for the Lorenz model with { x }={ x,y,z } are shown in Fig. 5 as functions of time t (left) and of the variable y (right), with the variable y ( t ) and z ( y ) of the trajectory at the bottom. On top, the ‘full’ case indicates that the angle between the principal axes ( a ) of the ‘full’ frame A and the main axis ( v ) of the ‘free’ frame V is a complicated function of the trajectory. At the centre of Fig. 5, however, the ‘constrained’ case is shown with the angle between the extreme orthogonal divergence direction and the main orthogonal axis of the frame W . This angle has only a small variation (centre left). It is also somewhat closer (centre right) to the shape of the attractor (bottom right). In both the O and RK method the starting frames V and W are ambiguous. This ambiguity can be removed if the corresponding local frame defined by the Jacobian deformation ellipsoid is used. For the Oseledec frame the set A x refers to the principal deformation directions. For RK the frame H x has as a first vector the flow direction, the remaining ones are the local extreme directions of orthogonal divergence. For RK this setting shortens the transient time until the W x frame is close to the unique local frame, see Fig. 5 centre, whereas A x and V x might be very different, see Fig. 5 top, resulting in a long transient. The Lorenz model will be used to study local exponents as functions of time t . Fig. 6 (bottom) shows the x component as a function of time for a short time interval. Above, local exponents are displayed for the same time interval, the left side for the ‘free’ ‘full’ case, the right side for the ‘constrained’ ‘reduced’ case (see Table 1). For the ‘full’ case, Fig. 6 (left, top) plots the first ‘free’ exponent λ V1 together with the main principal exponent α (fine line) as the local extreme value. The first exponent λ V1 is not at all times the largest and rarely has the extreme possible value of α . The difference ( λ V1 - α ), (left side, second row) has no relation to the x component (bottom). Fig. 6.
Lorenz chaos: Local exponents vs. time t . Bottom: x -component. Top left : first exponent λ V1 (thick line); main extreme exponent α (thin line). Right . first exponent λ W1 and extreme local separation β ┴ orthogonal to flow (nearly same values). Second row: Left : Difference ( λ V1 - α ); right: difference ( λ ┴W1 - β ┴ ) at the same scale as left side. Third row: thin line: main extreme exponent α Left : second exponent λ V2 ; right : second exponent λ W2 = λ ║ W = φ ║ acceleration along the flow. Note: Both local exponents λ V1 and λ V2 (left) consist of a mixing of separation λ ┴W1 and acceleration λ ║ W (right) . The third row, left, displays the second ‘free’ exponent λ V2 , again with α (fine line). Remember that for t →∞ the average of the first exponent is positive and the second average approaches zero. Both local exponents λ V1 and λ V2 are a mixing of divergence and acceleration.For the ‘constrained’ case, Fig. 6 (right), which discriminates between divergence and acceleration, Fig 6 (right, top) shows the main principal exponent orthogonal to the flow β ┴ thus indicating the possible maximum of divergence. The first exponent for divergence λ ┴W1 is very close to that maximum β ┴ k , with the small difference ( λ ┴W1 - β ┴ ) displayed below, related to the small deviation of the relative directions, as shown in Fig. 5 centre. However, the exponent for the local acceleration λ ║ W = φ ║ of Fig. 6 right, third row, can have values nearly twice as large as the exponent λ ┴W1 for divergence, but clearly never exceeds the main principal ‘full’ exponent α (fine line). At this point it is interesting to ask why the global exponents are equal for both frames as aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
12 / 24 demonstrated in Fig. 2, since their local behaviour is quite different. Assuming validity of the conjecture of Moser and Meier [20] that the direction for the smallest exponent is always orthogonal to d x for the V frame and the numerical result that it is parallel to the corresponding direction of the W frame, both smallest exponents are equal locally and globally. Knowing further that one global exponent must be zero for the Lorenz model and that the sum of all three exponents is equal, also the largest global exponents should be equal for the V and W frames.
5. Local exponents as functions of phase space R { x } vs. x component The temporal behaviour of the four local exponents has been extensively discussed in the previous sections. For a better understanding it seems worthwhile to explore their behaviour in the phase space R { x }. The directions of the Oseledec V -frame were already shown in 3-D plots and discussed by Wolf et al. [5] and Green and Kim [7] for a short part of the Lorenz trajectory. Here, longer parts of trajectories changing the loop several times will be displayed.. First the magnitude of the exponents will be shown as a function of the x component. Then the trajectories will be shown projected onto the yz plane if the values of the associated different exponents exceed a certain limit c . Fig. 7 shows on top a time series of the x component of ‘stable’ chaos. Below, left side, the exponents of the ‘free’ V frame are displayed to be compared to the exponents of the ‘constrained’ W frame. In order to show also the variation of the speed, the figures plot points at equal time intervals. First, the largest exponent λ Vi (left) has a rather complex behaviour as compared to the true divergence exponent λ ┴Wk (right). Moreover, close to the x value zero all values of the divergence exponent λ ┴Wk are positive, whereas the ‘free’ exponent λ Vi has a wide spread of positive and negative values in this x range around zero. This range is important, since here there are only trajectories which diverge from one loop to the other loop, as is easily seen on the display on top. The reason for the spread of values of λ V1 in this diverging range is its mixing of divergence with acceleration, which is strongly negative in this range. The local exponent λ ║ W = φ ║ is shown on the second row, right, together with the second exponent λ V2 left, again exhibiting a complex structure with respect to the x component. Fig. 7.
Lorenz chaos: Top: x component vs. time t for a longer time interval. Below: Local exponents vs. x component for the same time interval. Left : V frame. Right : W frame. First row λ V1 and λ W1 , below λ V2 and λ W2 = λ ║ W = φ ║ . Bottom row: λ V1 as a function of λ V2 (left), and λ ┴Wk as a function of λ ║ W = φ ║ (right). How is the relation between the first and the second exponent? In order to indicate this relation, the row at the bottom displays λ V1 as a function of λ V2 (left), and λ ┴Wk as a function of λ ║ W = φ ║ (right). In summary, in Fig. 7 all plots of the ‘free’ V frame (left) show a rather complicated structure compared with the plots of the ‘constrained’ W frame (right). The reason for this will be illustrated in the next figure. Figure 8 (top left) shows the angle of the first vector v of the ‘free’ V frame relative to the flow direction f ║ as a function of the x component. Clearly, this angle has a wide spread between being nearly parallel and nearly antiparallel to the flow direction f ║ . In both these extreme cases the exponent has a strong admixing of acceleration. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
13 / 24
In contrast, the first vector w ┴ of the ‘constrained’ W frame for divergence is by definition always orthogonal to the flow direction. Fig. 8.
Lorenz chaos: Top left: angle between first vector v of the ‘free’ V frame and flow direction f ║ vs. x component. Top right: angle between direction w ┴ of the ‘constrained’ W frame ( w ┴ always orthogonal to flow direction f ║ ) and direction b ┴ of the extreme separation exponent β ┴ vs. x component. Bottom right: exponent λ ┴W1 vs. extreme separation exponent β ┴ Straight diagonal line: limit λ ┴W1 = β ┴ . Bottom left: first ‘free’ exponent λ V1 vs. β ┴ , similar straight line λ V1 = β ┴ .In Fig 8 (top right) the angle of this direction w ┴ within the orthogonal plane relative to the direction b ┴ of the extreme divergence exponent β ┴ is shown to be always small as a function of the x component with a range between -16 to 8 degrees. Therefore, the values of the divergence exponent λ ┴W1 are never far from the value of β ┴ . Fig 8 (bottom right) shows the values of λ ┴W1 as a function of β ┴ with the straight diagonal indicating the limit λ ┴W1 = β ┴ . A similar diagonal λ ┴V1 = β ┴ is plotted in Fig. 8 (bottom left), showing the first ‘free’ exponent λ V1 also as a function of β ┴ . All the excess on the left side as compared to the right side has to be averaged out during integration in time t to finally describe only divergence. It seems worthwhile to compare in Fig. 9 the local exponents of the integrated frames V (left) and W (right) as functions of the strictly local maximum exponent α . Obviously, α is the upper limit, which is attained at specific points of the trajectory for the first (top) and the second (bottom) exponents, indicating that the local axes of the frames coincide with the main principal axis of the Jacobi deformation ellipsoid S at these specific points of the trajectory. Again the left side ( V ) has a more complex structure than the right side ( W ). An interesting question is: Where in the phase space occur the points where the exponents are close to the limiting value? A yz plot reveals that about the same number of such points is widely distributed for the V frames but is concentrated in a few small regions for the W frame. Fig. 9.
Lorenz chaos: Local exponents of the integrated frames V (left) and W (right) vs. strictly local maximum exponent α . Top: first exponents, bottom second exponents. Diagonal lines: exponent = α as a limit.
6. Local exponents and angles as functions of phase space R { x } in the yz plane Fig. 10 and Fig. 11 plot the projection of the trajectory onto the yz -plane when the local exponents exceed c with c = 0 in Fig. 10 and c = 4 in Fig. 11 Again the trajectory is plotted at equal time intervals in order to show the change of the velocity. The small circle denotes the maximal value. On both figures, the left side displays the ‘free’ ‘full’ case, the right side the ‘constrained’ ‘reduced’ case (see table 1). The top rows show the strictly local exponents α (general maximum) and β ┴ (maximum of divergence only). The centre display the largest exponents of the integrated frames V and W as λ V1 (mixing divergence and acceleration) and λ ┴W1 (divergence only), respectively. The bottom rows show the second exponents λ V2 (mixing acceleration and divergence) and λ ║ W = φ ║ (acceleration only). Comparing the ‘free’ exponents of the left side with each other, it is obvious that only α has clear cut borders of the limit c , whereas centre aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
14 / 24 and bottom exhibit a peculiar pattern of the limit. Furthermore, centre and bottom have regions where both exponents are above the limit.
Fig. 10 . Lorenz chaos: Projection of the trajectory at equal time intervals onto the yz -plane when the local exponents exceed c = 0. Left: ‘free’ case, right: ‘constrained case. Top: maximum exponent α and maximum separation exponent β ┴ . Centre: first exponent λ ┴V1 and λ ┴W1 ; bottom: second exponent λ V2 and λ ║ W = φ ║ , respectively. Fig. 11.
Lorenz chaos: Similar to Fig. 10, but exceeding the limit c = 4. The right side of ‘constrained’ exponents shows a very different behaviour. The limits are clear functions of the phase space. Moreover, the top and centre are almost similar, the centre pattern only slightly smaller. The centre and bottom patterns do not coincide. Finally, and most importantly, only the ‘constrained’ case right, top and centre, has large values where the real divergence of the two loops occurs. The discrepancy to the left side is easy to understand: Since acceleration is larger than divergence, and the first exponent λ V1 is a strong mixing of divergence and acceleration, the pattern bottom right λ ║ W = φ ║ is easy to see in the pattern centre left λ V1 , and the maximum (small circle) is nearly at the same place. Hence, only the ‘constrained’ case seems to analyse directly the local chaotic behaviour. The directions of both frames V and W are dependent on integration and not only on the location x . This brings about a dependence on the former sections of the trajectory. Although the exponents are unique for each location x , they are not simple functions of the phase space, since the former section of each location is different. The integration is a nonlinear procedure and, therefore, might have in itself a chaotic behaviour sensitive to small changes in the previous conditions. This implies a chaotic behaviour in addition to the chaotic behaviour of the analysed chaotic trajectory, resulting in the fact that another trajectory nearby might have a very different local exponent caused by a tiny difference at earlier points of that trajectory. However, this additional complexity is very different for the two frames, dependent on their ‘transient times’ as a measure to ‘forget’ earlier sections of the trajectory and on the range of changing angles of the frames. This complex behaviour is also present when angles of selected directions of these frames are plotted, as is performed in the next section. Figure 12 shows angles within certain bounds in the yz plane: left the angle of the first vector of the ‘free’ V frame relative to the flow direction; right the angle between the first orthogonal vector of the ‘constrained’ W frame relative to the direction of the extreme divergence direction orthogonal to the flow. The same angles were shown in Fig.8, top, but only as functions of the x component.The bounds imposed on the left are smaller than 10, between 60 and 120, and larger than 120 degrees, from top to bottom, respectively. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
15 / 24
The bounds on the right are much narrower, since this angle is between -16 and -4, between -2 and 2, and between 4 and 8 degrees, from top to bottom, respectively, with the total range being from -16 to 8 degrees, see Fig. 8, top right. Since
Fig. 12.
Lorenz chaos: Projection of the trajectory at equal time intervals onto the yz -plane when angles are within limits. Left: angle between first vector v of ‘free’ V frame and flow direction f ║ . The limits are smaller than 10, between 60 and 120, and larger than 120 degrees, from top to bottom, respectively. Right: angle between first orthogonal vector w ┴ of ‘constrained’ W frame and direction b ┴ of the extreme separation orthogonal to the flow. The limits are much narrower, since this angle is within the range from -16 to 8 degrees in the Fig. 8 top right: here between -16 and -4, between -2 and 2, and between 4 and 8 degrees from top to bottom, respectively. this range is small and the angle fluctuates on both sides (defined via the zero degree) of the extreme direction, the resulting exponents fluctuate much less than the exponents of the V frame. Furthermore, although the V frame changes sometimes through 90 degrees, see Fig. 12 left centre and Fig. 8 left top, and thus through the plane orthogonal to the flow, its direction within this plane has not to coincide with the orthogonal direction of the W frame at this point. However, since the third directions of both frames have been found to be nearly parallel, this plane is not far from the second direction of the W frame.
7. Local divergence in the phase space R { x } outside trajectories Since the ‘constrained’ exponent β ┴ of true divergence is only a function of the position x in the phase space and does, therefore, not depend on a trajectory, its value can be evaluated directly in the full phase space R { x }. These values are shown in Fig. 13 as a mesh plot for the xy plane at z values 60 (top left) and 20 (bottom left). The strange Lorenz attractor is about along the diagonal x ≈ y where the values β ┴ are low, but with a higher pass between the sections of the two loops. A simple local indicator of the curvature of trajectories in the phase space can be found easily: For the evaluation of the local acceleration, the new vector J f ║ after the time interval d t is projected by the scalar product ( f ║ , J f ║ ) onto the flow direction f ║ . The absolute magnitude of the new vector minus the absolute magnitude of the projection is a measure for the deviation at t+d t from the direction f ║ at t of the local trajectory and hence a measure of the curvature of the trajectory, d = │ J f ║ │- │ ( f ║, J f ║ ) │ (15) Fig. 13.
Lorenz chaos: The ‘constrained’ exponent β ┴ of true separation is only a function of the position x (independent of any trajectory) in the phase space R { x }. Its values are plotted as a mesh on an xy plane at z values 60 (top) and 20 (bottom) on the left side. The local measure d, eq. (15), of the curvature along the trajectory through x is plotted on the right side for the same z values. This ‘deviation’ d is shown in Fig. 13 on the right side. In more detail, the Jacobian matrix J can be written as sum of the symmetric matrix S and the anti-symmetric D : J = S + D , see sect. . The evaluated separate action of these matrices onto d reveals a very different behaviour: The symmetric S causes a strongly varying peculiar pattern in the region of the strange attractor, but its action is nearly negligible outside the attractor. In contrast, the anti-symmetric D creates a rather smooth structure, small at the attractor, but increasingly large with larger distance from the aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
16 / 24 attractor. The combined action of these different sources are visible in Fig. 13 on the right side as varying peculiar pattern along the diagonal ( S ), and smooth increase outside ( D ).
8. Why integrating exponents for short intervals between Poincaré points ?
At the beginning of analysing chaos, the main interest was to find general numbers defined in the limit of t→∞ , which presupposes a ‘stable’ chaos. Later on, also chaotic motions were analysed where the strength of the chaotic behaviour was changing. Therefore, shorter time intervals were used. The large variation of the local exponents in time made the resulting numbers strongly dependent on the limits where the time intervals start and end. In order to test different methods, a chaotic behaviour with strongly changing character would be welcome. Such a ‘transient’ chaos will be described in the next section.
Until now, the Lorenz chaos with parameters ( σ,ρ,β ) = (16,40,4) was ‘stable’ with repelling unstable fixed points at X =0, where trajectories starting nearby would spiral out, join the well-known double loop strange attractor and stay there in theory to infinity, in practice until the build-up of computing errors will be too high. The top of Figure 14 displays a section of the ‘stable’ chaos where the trajectories are plotted on the yz plane at equal time intervals to show the variation of the speed as in earlier plots. A change of the parameter ρ in the Lorenz equation (d y /d t = - x z + ρ x – y ) from 40 to 28.165 replaces the unstable fixed points by attracting stable fixed points. Figure 14 bottom displays a peculiar ‘transient’ Lorenz trajectory starting with several loops on the right side with increasing radius, followed by a chaotic interval with loops on both sides, and a final decay spiralling on the left side to a stable fixed point. p between Poincaré points for Lorenz trajectories For the Lorenz chaos, integrating over one loop would be reasonable. Since there are no closed loops, well chosen Poincaré points will be used to define the start and end of a single loop. Poincaré points will be defined here when a trajectory is crossing a plane at y = c p ‘from above’, i.e. for decreasing values of y . The time interval τ p is then defined between two consecutive Poincaré points. The plane parameter c p is chosen at a level where the trajectories are about normal to the Poincaré plane. Figure 14 top shows this plane ( c p =40) as a Fig. 14.
Lorenz model: ‘stable chaos’ (above) and ‘chaotic transient to a fixed point’ (bottom). Trajectories are projected onto the yz plane at equal time intervals to indicate the changes in velocity. The lines indicate the position of the Poincaré planes z = c p =const. with c p =40 (above) and c p =27.165 (below). Over each time interval between consecutive Poincaré points local exponents are integrated. line in the yz plot for a ‘stable chaos’ with instable fixed points, whereas Fig. 14 bottom displays a similar line ( c p =27.165) for stable attracting fixed points with the ‘transient’ trajectory described above. The integration between Poincaré points is therefore performed over one loop. The results are shown in the next sections for the ‘free’ and ‘constraint’ case, first for the ‘stable’ Lorenz chaos. These results will be compared with the integration of all possible local directions, which are evaluated with the local Jacobian deformation ellipsoid. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
17 / 24
9. The idea of integrating all directions of the local sphere
It seems worthwhile to recall here the basic idea behind the definition of Lyapunov exponents, which consists of integrating along a trajectory the subsequent deformations of the surface of a sphere, surrounding a given initial condition, along all direction in phase space. At the final time, the principal axes and eigenvalues of this deformed ellipsoid are considered to be the relevant characteristics of the problem, thus only n exponents are sufficient to describe the type of dynamical instability. It was proved by Oseledec that only the knowledge of the deformations of a frame of n orthogonal directions is necessary to find these final axes. The evaluation of the deformation with the Jacobian ellipsoid S enables to test this established fact numerically. It is straightforward to evaluate S [ x ( t )] along a trajectory x ( t ). However, for calculating Lyapunov exponents from it the rotations of the directions after each time step d t still have to be eliminated. This important fact is not too obvious, since the expression J u of eq. (3) implies elongation and rotation. However, in order to find instantaneous local exponents describing only elongation, the rotation must be projected out by using the scalar product of eq. (5). Moreover, because this scalar product has the same vector on both sides, it enforces the symmetrisation of the Jacobian matrix J in form of the symmetric Jacobian deformation ellipsoid S . Therefore, a rotation of S has to be performed after each time step d t along the trajectory. In principle, each direction should be rotated according to eq. (1) separately. However, as a crude oversimplification all orientations might be rotated equally according to the rotation of an n dimensional orthogonal frame evaluated with the equation of motion eq. (1). Along these lines, transformations S into the frames V or W will be used. Again, the average of values computed at discrete time steps will approximate the continuous time integration. Let the interval τ be divided into m parts Δ t. With x k = x ( t + k Δ t ), and first for the frame V, the resulting arithmetic average S ( τ ) V can be written as S (τ) V = (1/ τ )∑ mk= V T ( x k ) S x ( x k ) V ( x k ) (16) Secondly, in complete analogy the average S (τ) W is defined for the RK frame W with the first vector constrained along the flow direction, S (τ) W = (1/ τ )∑ mk= W T ( x k ) S x ( x k ) W ( x k ) (17) Thirdly, another approach consists of considering only all the directions orthogonal to the flow as described by the reduced Jacobian ellipsoid S ┴ (2,… n ) , obviously to be transformed by a reduced matrix W ┴ (2,… n ) , resulting in the average S ┴ ( τ ) , S ┴W (τ) = (1/ τ )∑ mk= W ┴T ( x k ) S x ( x k ) W ┴ ( x k ) (18) At the final time of the integration, for each of the above averages the eigenvalues are evaluated and compared with the integrals of the exponents of the corresponding frames. For convenience, the eigenvalues will be ordered with the largest first for the principal axes of the averaged sums of the ellipsoids. In contrast, the averaged sums of the local exponents of the frames are ordered according to the order during the orthogonal normalization process, which might result in a reversed order, such that the value for the first exponent is smaller than for the second exponent.
10. Integration of all directions and integration of n local exponents for ‘stable’ chaos Figure 15 shows for ‘stable’ Lorenz chaos, see the variable x ( t ) at the bottom, the resulting first and second exponents, evaluated as eigenvalues of the final Poincaré integration matrix S (τ) V (top) of eq. (16) , S (τ) W (second row) of eq. (17), and S ┴ W ( τ ) (third row) of eq. (18), see symbols o and *, respectively. They are compared to the corresponding integrals of the first two exponents, evaluated as averages of the instantaneous local exponents in the directions of the frames V and W , see symbols x and + in the figure , respectively . For the first case S (τ) V (top), large discrepancies are visible. For the second case S (τ) W (second row), there exist rather small deviations. Only for the third case S ┴ W ( τ ) (third row), the eigenvalues of the integral of the reduced deformation ellipsoid give the same values as the integral of the instantaneous local exponents of the W ┴ frame. The discrepancies might be caused by the oversimplification of rotating all orientations equally. However, the most important result is that the RK method yields values that are equivalent to the method based on computing aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
18 / 24 eigenvalues by integrating all directions in the n -1 space orthogonal to the flow, even when these orientations are rotated equally but within the orthogonal subspace. The second exponent is by definition equal for RK and S ║ W ( τ ) . The small but distinct discrepancies between the full S (τ) W (second row) and S ┴ W ( τ ) (third row) are easy to understand: Already at each point x the eigenvectors of the reduced matrix S ┴ and the direction of the flow d x only accidentally coincide with the principal axes of the full local S Fig. 15. ‘Stable’ Lorenz chaos; Bottom: x ( t ). Above: Integrals of deformations between consecutive Poincaré points and subsequent eigenvalues for S (τ) V (top) , S (τ) W (second row) and S ┴W ( τ ) (third row) as o and *, respectively, together with the integrals of the first two exponents of the frames V and W as x and + , respectively. Lines are only guides to the eye. . The same comparison is made after a long integration time τ when the sums are close to the asymptotic values. First, using eq. (16), the ellipsoid S (τ) V of the Oseledec type V frame of free rotations has the diagonal elements (1.59, 0.002, -22.59). According to the conjecture for Lyapunov exponents, the global ellipsoid should have principal axes parallel to the directions of the V frame. Equally important, the eigenvalues of S (τ) V should be equal to the global exponents evaluated by the integration of the instantaneous local exponents in the local directions of the frame V, resulting in the values (1.38, 0.002, -22.38). Indeed, there is a small but distinct difference for two values. Moreover, the directions of the principal axes of S (τ) V deviate by about 5 degrees from the directions of the V frame. Secondly, using eq. (17) the ellipsoid S (τ) W has the eigenvalues (1.37, 0.0001, -22.37), exactly as the corresponding global exponents of the W frame (1.37, 0.0001, -22.37). Moreover, the corresponding directions differ by less than 1 degree. Note that for short intervals as described in the previous section, for S (τ) W the corresponding values are slightly different. Hence, for the global values these small discrepancies are levelled out for long integration time τ . At this point it seems interesting to discuss the question: Are the global ‘Lyapunov eigenvectors’ fixed in time? The answer is yes and no , depending on the observer. For a static observer at the origin of the coordinate system of R { x }, the eigenvectors are fixed only if no rotating system U ( t ) is used. In reality, these eigenvectors are rotating as fast as U ( t ) itself. Only a rotating observer at the origin of U ( t ) sees the Lyapunov eigenvectors fixed in time. Moreover, although the final ellipsoid is fixed accordingly for this rotating observer, he would recognise that even after a very large interval τ the instantaneous local ellipsoid is still deforming as rapidly as after a short time.In conclusion, the approach by Lyapunov as already been described in sect. 1.2 (iii) proposes that the global exponents evaluated only for n local orthogonal directions are the eigenvalues of the global deformation ellipsoid, which contains all directions. This idea has been tested numerically. However, it is only confirmed if the constrained frame W of RK is used.
11. Poincaré integrals for ‘transient’ chaos with a sudden flip of the frames V frame Figure 16 bottom shows the x component vs. time t for the ‘transient’ chaos described above. On top, the Poincaré integrals for the first two exponents of the V frame are plotted as o and +, respectively. The centre shows the angle of the first axis of the V frame relative to the flow direction. The frame starts with an arbitrary orientation. Then, at t = 7 the frame is suddenly aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
19 / 24 rotated by 90 degrees. After this flip the reorientation to being nearly parallel to the flow is rather slow. During the short strongly chaotic regime the angle flips up and down. Long after the beginning of the spiralling down to the fixed point the angle reorients again between t ≈ 22 and 27. This delayed motion causes a virtually chaotic pattern of the exponents, see top panel within approximately the same time interval, similar to Fig. 16.
Lorenz ‘chaotic transient to a fixed point’; Bottom: x ( t ). Centre: Angle of the first axis of the V frame relative to the flow direction. This angle is artificially rotated by 90 degrees at t = 7. Top: Poincaré integrals of the first two exponents of the V frame denoted by the symbols o and +, respectively. Lines are only guides to the eye. Note the ‘artificial bump’ between t ≈ 22 and 27 as explained in the text .the reaction at the beginning and after the flip. Furthermore, note that during the strong chaotic behaviour the integral of the second exponent is sometimes larger than the first exponent. Clearly, the Oseledec V frame has problems with delayed reorientation, which could give wrong exponents for finite time intervals. RK W frame Figure 17 bottom displays the x component of the same ‘transient’ chaos as in Fig. 16. Here the frame is also suddenly rotated by 90 degrees at t =7, as shown in the centre for the angle between the first orthogonal direction of the W frame and the direction of the local orthogonal direction of extreme divergence. This flip changes the exponents only for the integral over one loop, as is shown on top for the second exponent, since the recovery time of the angle is small, and one direction along the flow remains fixed according to the definition of the ‘constrained’ W frame. This implies that the direction for the first exponent is independent of the orientation of the other directions of the W frame, so the flip at t =7 does not affect the first exponent, see + in Fig. 17 top panel, Fig. 17 . Lorenz ‘chaotic transient to a fixed point’; Bottom: x ( t ). Centre: Angle between the first orthogonal direction of the W frame and the direction of the local orthogonal direction of extreme separation. The W frame is artificially rotated by 90 degrees at t = 7. Top: Poincaré integrals of the first two exponents of the W frame denoted by the symbols o and +, respectively. The dots plotted above the o symbols are the integrals of the local extreme orthogonal exponents, which are independent of orientation and flipping of the W frame. It is also interesting to compare the RK Poincaré integrals with the integrals of the local extreme orthogonal exponents evaluated as largest eigenvalues of the local reduced Jacobian ellipsoid, see the dots above the symbols for the second exponent in Fig. 17. Since these integrals are for local exponents, they are independent of the orientation of any frame. Their values are always larger than the RK values, but they provide also an excellent description of the trajectory. They are small during the last part of the spiralling dynamics to a fixed point but still slightly positive there, which indicates chaotic behaviour. However, the negative integral of the exponent along the acceleration simultaneously indicates non-chaotic behaviour. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited
20 / 24
Hence the combination of both the local extreme exponent with the exponent along the acceleration leads already to a good estimate of the type of regime, chaotic or non-chaotic. Furthermore, both exponents have zero transient time, they are correct already at the staring time t , Since no frame V or W is needed, the ambiguity of a starting frame, their integration and their transient times are avoided.
12. Quantitative analysis using Poincaré point correlations
Qualitatively, the previous figures seem to favour the RK method: The structure of the time-dependent dynamics of the exponents is simpler than the one of the established Oseledec method. Hence, it seems worthwhile to quantitatively compare both methods. There are ways to test the complexity of their local exponents, see Chlouverakis et al. [20]. However, the local exponents themselves are merely tools to analyse the chaos of a given nonlinear dynamics. The question is therefore rather: How well are the given trajectories assessed by the local exponents? How well are the exponents correlated to the trajectories? How can the degree of the associated complexity be evaluated? A method similar to an established method to compare the volatility of financial values will be constructed, namely the evaluation of the standard deviation of the distribution of the values.
It seems sufficient to analyse only a typical set of points in a reduced n-1 dimensional subspace. For Lorenz chaos, an xy -plane will be defined by z =constant. Pairs of neighbours of the points in this plane are selected and their distance r is evaluated. Then, the absolute difference Δ λ of the corresponding local exponents is evaluated as a function of the distance r of these neighbouring points i,k . Δ λ ik ( r ) = │ λ i – λ k │ with r = ( x i - x k ) + ( y i - y k ) (19) Finally, plots of the strictly local exponents are compared to the exponents of the V and W frames which incorporate integrals over previous times. This comparison can be performed in a quantitative way and results in a measure of how well the V and W frames correspond to the trajectories to be analysed. The strength of chaotic behaviour is manifest in the ‘pseudo’-irregularity of the differences Δ λ ik ( r ), which is measurable by its standard deviation, after a possible underlying regular behaviour has been eliminated (we call it ‘pseudo’-irregular, since all data points are regularly determined by eq. (1), including the rotations of the frames).The following quantitative procedure will be used to evaluate the distribution of the data depending on the function Δ λ ik ( r
21 / 24 considered. A plot of these Poincaré points in the Poincaré xy- plane reveals a particular structure, see Fig. 18 (top): There are two branches of p (-) points around the centre, one with y>x, denoted by c1, and the other with y Lorenz chaos: xy plot of Poincaré points at z =40. Top: all points. Bottom: blowup of the central sections for decreasing z values. Upper left branch c1: symbol o denotes the sequence ℓ→ℓ(left torus), symbol + denotes the sequence ℓ→r(right torus). Lower right branch c2: * for r→r and x for r→ℓ. goes to ℓ or r by analysing the time series p (+) ( t i-1 ), p (-) ( t i ), p (+) ( t i+1 ). The result is striking: the points c1 display only the sequences ℓ→ℓ or ℓ→r plotted in Fig. 18 (bottom) as circles ○ and x signs, respectively. Moreover, these two possibilities are separated in the Poincaré plane. Similarly, the c2 points consist only of an r→r or r→ℓ sequences, plotted as stars * and + signs, respectively. This symmetry allows restricting the analysis to the c1 branch. At this point an interesting question can be asked: Where on the trajectory could additional noise most easily change the dynamical behaviour? As an example, perturb the trajectory at z =40 where the ℓ→ℓ sequence is close to the ℓ→r sequence, thus the circles ○ and x are very close, such that the trajectory jumps from one to the other sequence. At this position a noise-induced transition as described by Gassmann [21] would be most easily possibleTo provide an overview, the local exponents of the c1 points are plotted as a function of the x component in Fig. 19. The separation between the ℓ→ℓ and ℓ→r sequences is marked by a vertical line. The strictly local exponents are: extreme expansion α (top left), extreme orthogonal divergence β (top right). The second exponent λ ║ W of the W frame of local acceleration is also a strictly local exponent (bottom right). The central row displays the first exponents λ V and λ ┴ W , the bottom row the second exponents λ V and λ ║ W of the V frame (left) and of the W frame (right), respectively. It is interesting that all three plots on the right side (extreme orthogonal divergence β and W frame) share the same feature: The values of the exponents all differ for the two sequences, the barrier is marked by a horizontal line indicating a complete correlation to the behaviour of remaining on the same left loop, or changing the loop, which certainly implies a larger divergence. This feature is missing for both exponents of the V frame (left centre and left bottom) showing a large spreading. Fig. 19. Lorenz chaos: Local exponents of c1 points ℓ→ℓ (o) and ℓ→r (+) of Fig. 18 vs. the x component. Vertical lines separate o from +. Left: top extreme expansion α i , centre λ V1 , bottom λ V2 . Right: top extreme orthogonal expansion β ┴ , centre λ ┴W1 , bottom λ ║ W = φ ║ . Horizontal lines separate values of o from +. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited 22 / 24 Fig. 20. Lorenz chaos, PPCD analysis: Absolute difference │Δ λ ik ( r )│of the corresponding local exponents λ vs. radius r Lorenz chaos: The same PPCD data points as shown in Fig. 20, but all plotted on the same scale to indicate the differences of the magnitude of the V frame (left centre and bottom) to the W frame (right centre) and to the strictly local exponents. Therefore, the quantitative factors f defined by eqs. (21) and (22) are f V ≈ 100 (23) f W ≈ 1 (24) Although the heuristic PPCD analysis provides a strong oversimplification by using only a small selected portion of the data, its outcome is surprisingly clear. In order to visualize this final result, Fig. 20 is potted again as Fig. 21, but now with the same scale for all plots. Obviously, only the Oseledec V frame (left centre and bottom) is adding a strong chaotic complexity to the complexity of the trajectories to be analyzed. 13. Conclusions This work suggests that the evaluation of Lyapunov exponents should be revisited. The introduction of the symmetric deformation Jacobian ellipsoid and its submatrix orthogonal to the flow allows the direct determination of the principal exponents and of the extreme local exponent for diverging trajectories for every point in the phase space without the need to integrate along a specific trajectory to find the local aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited 23 / 24 ‘Lyapunov directions’ according to the procedure of Oseledec. Moreover, a fundamentally different approach evaluates the divergence of nearby trajectories by avoiding a mixing with local acceleration. This is performed by adding a simple constraint to the Oseledec frame as already done in Ref. [10]: The first vector remains fixed along the flow direction d x , and any further vectors are subsequently orthogonalised after each rotation, which implies that only the divergence and not any partial acceleration are assessed at every point on the trajectory. This avoids the additional complexity introduced by the Oseledec method, and the local exponents are in accord with the local divergence of the trajectories. The main problem of dynamical instability: Remain neighbouring trajectories in a tube? is thus solved. The largest exponent indicates locally already where trajectories diverge and where the origins of instability occur, Only the evaluation of the Jacobian deformation ellipsoid is treated analytically, the remaining conclusions are based on studying numerical examples. Hence, a formal revision of the original idea of Lyapunov to evaluate the divergence of nearby trajectories instead of considering neighbouring points remains to be carried out.We finally remark that there might exist interesting crosslinks between our work and the very recent active field of calculating and understanding so-called Lyapunov modes in interacting many-particle systems, see Ref. [22] for a short review and further references therein. Lyapunov modes refer to the eigenmodes associated with the spectrum of Lyapunov exponents which are closest to zero, projected onto the single particles from which they originate. They were found to form interesting spatio-temporal periodic patterns. The methods developed in our paper, based on solving the eigenvalue problem for the local Jacobian deformation ellipsoid in a suitable local coordinate system, could possibly serve for developing alternative techniques of computing such Lyapunov modes. Using our methods, it might also be interesting to check for spatio-temporal structures in the corresponding distribution of local Lyapunov exponents in such systems. Acknowledgements F.W. would like to thank P.F. Meier, H.R. Moser, and E.P. Stoll for valuable help. References [1] A.M. Lyapunov, doctoral dissertation on stability of motion, St. Petersburg (Russian) (1892), Ann. Faculté des sciences, Université de Toulouse, 2 ser. (1907), Ann. Math. Study 17 (1977) (Princeton). [2] V.I. Oseledec, A multiplicative ergodic theorem, Lyapunov characteristic numbers for dynamical systems, Trudy Moscow Mat. Obsc. 19 (1968) 179-210. English translation Trans. Moscow Math. Soc. 19 (1968) 197-221. [3] G. Benettin, L. Galgani and J-M. Strelcyn, Kolmogorov entropy and numerical experiments, Phys. Rev. A 14 (1976) 2338-2345.[4] J.D. Farmer, E. Ott, and J.A Yorke, The dimension of chaotic attractors, Physica 7D (1983) 153-180. [5] A. Wolf, J.B. Swift, H.L. Swinney and J.A. Vastano, Determining Lyapunov exponents from a time series, Physica 16D (1985) 285-317, and references therein. [6] I. Goldhirsdch, P.-L. Sulem, and S.A. Orszag, Stability and Lyapunov stability of dynamical systems: a differential approach and a numerical method, Physica 27D (1987) 311-337.[7] J. M. Greene and J.-S. Kim, The calculation of Lyapunov spectra, Physica 24D (1987) 213-225. [8] J.P. Eckmann and D. Ruelle, Ergodic theory of chaos and strange attractors, Rev. Mod. Phys. 57 (1985) 617-656, and references therein. [9] Ch. Skokos, The Lyapunov characteristic exponents and their computation, Lect. Notes Phys. 790 (2010) 63-135. [10] K. Rateitschak and R. Klages, Lyapunov instability for a periodic Lorentz gas thermostated by deterministic scattering, Phys. Rev. E (2002) 036209/1-11. [11] B. Eckhardt and D. Wintgen, Indices in classical mechanics. J. Phys. A: Math. Gen. 24 (1991) 4335-4348. [12] P. Gaspard, Chaos, Scattering and Statistical Mechanics, Cambridge Nonlinear Science Series (No.9) Cambridge 1998, 1-495. [13] Ch. Dellago and Wm. G. Hoover, Are local Lyapunov exponents continuous in phase space? Phys. Lett. A 268 (2000) 330-334. aldner, Klages Jacobian deformation ellipsoid and Lyapunov stability analysis revisited 24 / 24 [14] E.N. Lorenz, Deterministic nonperiodic flow, J. Atmos. Sci. 230 (1963) 130-141.[15] O.E. Rössler, An equation for hyperchaos, Phys. Lett. 71A (1979) 155.[16] C.G.J. Jacobi (1830), Gesammelte Werke, 7 vol., G. Reimer, Berlin (1881-1889). [17] E. Fick, Einführung in die Grundlagen der Quantentheorie. Aula-Verlag, Wiesbaden, 6 thth