Investigating Chaos by the Generalized Alignment Index (GALI) Method
aa r X i v : . [ n li n . C D ] J u l Investigating Chaos by the Generalized Alignment Index(GALI) Method
Henok Tenaw Moges niversity of Cape Town
Faculty of Science
Department of Mathematics and Applied Mathematics
Investigating Chaos by theGeneralized Alignment Index (GALI)Method
Author:
Henok Tenaw Moges
Supervisor:
A/Prof Haris SkokosMarch, 2020 bstract
One of the fundamental tasks in the study of dynamical systems is the discrimination betweenregular and chaotic behavior. Over the years several methods of chaos detection have beendeveloped. Some of them, such as the construction of the system’s Poincar´e Surface ofSection, are appropriate for low-dimensional systems. However, an enormous number ofreal-world problems are described by high-dimensional systems. Thus, modern numericalmethods like the Smaller (SALI) and the Generalized (GALI) Alignment Index, which canalso be used for lower-dimensional systems, are appropriate for investigating regular andchaotic motion in high-dimensional systems. In this work, we numerically investigate thebehavior of the GALIs in the neighborhood of simple stable periodic orbits of the well-knownFermi-Pasta-Ulam-Tsingou lattice model. In particular, we study how the values of theGALIs depend on the width of the stability island and the system’s energy. We find thatthe asymptotic GALI values increase when the studied regular orbits move closer to the edgeof the stability island for fixed energy, while these indices decrease as the system’s energyincreases. We also investigate the dependence of the GALIs on the initial distribution of thecoordinates of the deviation vectors used for their computation and the corresponding anglesbetween these vectors. In this case, we show that the final constant values of the GALIs areindependent of the choice of the initial deviation vectors needed for their computation.
Plagiarism declaration
I know the meaning of plagiarism and certify that this dissertation is my work, based on mystudy and I have acknowledged all the materials and resources used in the preparation.
Henok Tenaw Moges cknowledgments First and foremost, I would like to express my sincere gratitude to my supervisor, A/ProfHaris Skokos, for his continuous support and guidance. I am so grateful to have a supervisorwho cares about my work. I would also like to thank A/Prof Thanos Manos for severaluseful discussions. I will forever be thankful to my colleagues in the ‘Nonlinear Dynamicsand Chaos’ group; Bob Senyange, Many Manda, Malcolm Hillebrand, and the former groupmember Chinenye Ani. Thank you, guys. I would like to pay my special regards to theDepartment of Mathematics and Applied Mathematics of the University of Cape Town (UCT)as well. This work was fully funded by the Woldia University and the Ministry of Scienceand Higher Education of Ethiopia and partially supported by the UCT International andRefugee Grant 2019. Finally, I would like to thank the High-Performance Computing facilityof UCT and the Center for High-Performance Computing for providing the computationalresources needed for this work. 2 ontents
Abstract 1Acknowledgments 2Abbreviations 4Introduction 51 Hamiltonian mechanics 7 ist of abbreviations • dof : degrees of freedom • FPUT : Fermi-Pasta-Ulam-Tsingou • GALI : Generalized Alignment Index • IC : Initial Condition • LE : Lyapunov Exponent • mLE : maximum Lyapunov Exponent • N D : N -Dimensional • PSS : Poincar´e Surface of Section • SALI : Smaller Alignment Index Method • SI : Symplectic Integrator • SVD : Singular Value Decomposition 4 ntroduction
Dynamical Systems theory attempts to understand, or at least describe, changes over timethat occur in physical and artificial systems. It consists of a set of subaltern theories that canbe applied to solve complex problems in numerous fields, such as weather predictability [1],stochastic processes [2] and evolutions of astronomical systems [3]. Because of the remarkableapplicability of this theory, it has become one of the top topics in modern scientific research.In particular, Hamiltonian systems of many degrees of freedom (dof) have been broadlyused for a general study and a better understanding of energy transport and equipartitionphenomena (see for example [4–8] and references therein). Equipartition phenomena areassociated with the dynamical nature of motion fostered in the respective phase space namely,the regular or chaotic evolution of the system’s orbits [6]. A chaotic dynamical system is onein which nearby initial conditions (ICs) lead to very different states as time evolves. Thismeans that small variations in ICs can lead to large, and frequently unpredictable, variationsin final states.A fundamental aspect in studies of dynamical systems is the identification of chaoticbehavior, both locally, i.e. in the neighborhood of individual orbits and globally, meaningfor large samples of ICs. There are several ways to identify the chaotic behavior of a givennonlinear system. The rapid and efficient detection of the regular or chaotic nature of motionin many dof systems has been a very active research topic over the years. The most commonlyused method to characterize chaos is the computation of the Lyapunov Exponents (LEs)[9–13]. In general, LEs are asymptotic measures characterizing the average rate of growthor shrink of small perturbations to the orbits of dynamical systems. The development ofefficient algorithms for the calculation of LEs [11, 12] (see also [13] for a recent review on thistopic) led to a gallery of chaos detection methods (see e.g. [14–17]). The most commonlyemployed chaos indicators were recently reviewed in [17].This study focuses on the Generalized Alignment Index (GALI) method [18], anefficient and fast chaos detection technique that can be used effectively for the descriptionof the behavior of multidimensional systems. As its name would suggest, the GALI is ageneralization of the Smaller Alignment Index (SALI) [19]. The main advantages of the GALImethod are its ability to distinguish between regular and chaotic motion more quickly thanother techniques and to determine the dimensionality of the torus on which the quasi-periodicmotion occurs. In addition, it can predict slow diffusion in multidimensional Hamiltoniansystems [20]. A concise overview of the theory and numerical applications of the SALI andthe GALI methods can be found in [21].Here we aim to understand the GALI’s prediction efficiency and performance when theIC of an orbit moves gradually from a regular to a chaotic motion region in Hamiltoniansystems, by studying in detail the variation of the asymptotic values of the indices. Thistransition can take place either by changing the total energy of the system or by movingthe regular orbit’s IC towards the edge of a stability island. In addition, we follow the time5volution of the deviation vectors, which are small perturbations of the IC and are needed forthe computation of the GALIs, and examine if different initial distributions of the coordinatesof these vectors affect the final GALI value. Our investigation was performed in phasespace regions around some simple, stable, periodic orbits of the Fermi-Pasta-Ulam-Tsingou(FPUT) model [22, 23], which describes a chain of harmonic oscillators coupled throughnonlinear interactions. The dynamics of this system has been studied extensively in the lastdecades (see e.g. [24, 25] and references therein) and this model is considered nowadays as aprototypical, multidimensional, nonlinear system.Moreover, an investigation of the behavior of the GALIs for regular orbits of coupledarea-preserving standard mappings [26] was carried out. More details on this mapping canbe found in [27]. In some sense, our study completes some previous works on the GALImethod [18, 20] where other aspects of the index have been investigated, for example, itsbehavior for periodic orbits [28].The thesis is organized as follows:1. Chapter 1 provides a general overview of the general theory of Hamiltonian mechanicsand chaotic dynamics. The Hamiltonian formulation of finite-dimensional systems isalso presented.2. In Chapter 2, we present the symplectic integration methods used for our computationsalong with the introduction of several chaos detection techniques, namely the Poincar´eSurface of Section (PSS), the LEs, as well as the SALI and the GALI methods.Particular emphasis is given to the properties and the numerical computation of theGALI method. In order to illustrate the behavior of the LEs, as well as the SALIand the GALI methods for regular and chaotic motion, simple Hamiltonian modelsand symplectic mappings are used: two-dimensional (2D) and four-dimensional (4D)mappings, as well as the 2D H´enon-Heiles system [29]. Moreover, the behavior of theseindicators for dissipative dynamical systems is briefly discussed.3. In Chapter 3, aspects of the behavior of the GALI method in the case ofmultidimensional conservative Hamiltonian systems are investigated in detail. Mostof the presented numerical simulations are performed using one of the most classicalsystems of statistical physics: the FPUT model [22]. The equations of motion andvariational equations of this model are presented. The asymptotic behavior of the GALIfor regular orbits in the neighborhood of two simple periodic orbits of the β − FPUTmodel is analyzed. An investigation of the dynamical changes induced by the increaseof the coupling between the 2D mappings and their influence on the behavior of theGALI using a coupled standard mapping [21, 27] is also undertaken.4. In Chapter 4, we finally summarize and discuss the results and findings of our study.6 hapter 1Hamiltonian mechanics
The phase space of a dynamical system is the collection of all possible states of the system inquestion. Each state represents a complete snapshot of the system at some moment in time.The dynamical evolution of the system is governed by rules which transform one point in thephase space, representing the state of the system “now”, into another point representing thestate of the system one time step “later”.There are two different classes of dynamical systems:1. Discrete dynamical systems. They are described by recurrence relations iteratedmappings or sets of difference equations x n +1 = f ( x n ) , (1.1)where f is a set of n functions and x n is the state vector x at the discrete time t = n , n ∈ Z .2. Continuous dynamical systems. They are described by differential equations˙ x = dxdt = f ( x ( t )) , (1.2)where x ( t ) ∈ R m is a vector of state variables, f : R m −→ R m , m ∈ N is avector field, and ˙ x is the time-derivative, which we can write as dxdt . We mayregard (1.2) as describing the evolution in continuous time t of a dynamical systemwith finite-dimensional state x ( t ) of dimension m . In component form, we write x = ( x , . . . , x m ), f ( x ) = ( f ( x , . . . , x m ) , . . . , f m ( x , . . . , x m )) and the system (1.2) isgiven as ˙ x = f ( x , . . . , x m ) , (1.3)˙ x = f ( x , . . . , x m ) , ...˙ x m = f m ( x , . . . , x m ) . Notice that, here we do not use notations such as ~x to explicitly denote vectors. In addition, we alsowrite vectors as row or column matrices. Lets briefly discuss the definition of chaos.
Chaos theory in mechanics and mathematics,studies the apparently random or unpredictable behavior in systems governed bydeterministic laws. A more accurate term, deterministic chaos , suggests a paradox becauseit connects two notions that are familiar and commonly regarded as incompatible [30]. Thefirst is that of randomness or unpredictability, for example in the trajectory of a molecule ina gas or the voting choice of a particular individual in a population. Usually, randomnessis considered more apparent than real, arising from the ignorance of several causes in anexperiment. In other words, it is commonly believed that the world is unpredictable becauseit is too complicated. The second notion is that of deterministic behaviors/laws, like forexample the ones appearing in the motion of a pendulum or a planet, which have been usedsince the time of Isaac Newton, exemplifying the success of science in predicting the evolutionof cases which are initially considered complex.
Devaney’s definition of chaos : we present here a formal definition of chaos followingthe presentation of [31]. Let ( χ, d ) be a metric space, where χ is a set which could consist ofvectors in R N and d : χ × χ → R is the distance or metric function. A function f : χ → χ iscalled chaotic if and only if it satisfies the following three conditions:1. f has sensitive dependence on ICs. This means that ∃ δ > U and for any point x ∈ U , there exists a point y ∈ U such that d ( f k ( x ) , f k ( y )) > δ forsome k ∈ N , where f k denotes k successive application of f . The positive real number δ is called a sensitivity constant and it only depends on the space χ and function f .2. f is topologically transitive. Thus, for any two open sets U and V there exists k suchthat f k ( U ) ∩ V = ∅ .3. The set of periodic points of f is dense. A point x is called periodic if f k ( x ) = x forsome k ≥ Newton’s second law gives rise to systems of second-order differential equations in R n andso to a system of first-order equations in R n , i.e. in an even-dimensional space. If the forcesare derived from a potential function, the equations of motion of the mechanical system havemany special properties, most of which follow from the fact that the equations of motioncan be derived from a Hamiltonian formulation. The Hamiltonian formalism is the naturalmathematical structure in which to develop the theory of conservative dynamical systems.8HAPTER 1. HAMILTONIAN MECHANICSIn a Hamiltonian system with N dof, the dynamics is derived from a Hamiltonianfunction H ( p, q, t ) , (1.4)where q = ( q , . . . , q N ) and p = ( p , . . . , p N ) are respectively the system’s canonicalcoordinates and momenta. An orbit in the 2 N -dimensional (2 N D) phase space S of thissystem is defined by the state vector x ( t ) = ( q ( t ) , q ( t ) , . . . , q N ( t ) , p ( t ) , p ( t ) , . . . , p N ( t ) , t ) , (1.5)where x i = q i , x i + N = p i , and i = 1 , , . . . , N . In many cases, the Hamiltonian does not explicitly depend on time and the energy E = H ( q, p ) , (1.6)is conserved along trajectories, i.e. it is a constant of motion [32]. Examples of such systemsare the pendulum, the harmonic oscillator, dynamical billiards or the motion of a particlewith mass m in a potential V ( q ) described by the Hamiltonian H ( p, q ) = p m + V ( q ) . One might ask, when is the energy conserved? The answer to this question is given bythe Noether‘s theorem [33]. This theorem applies to any action in classical mechanics weatherit is described in the Lagrangian or Hamiltonian formulation. Focusing on Hamiltoniansystems, let us consider a general action I described by I [ p i , q i ] = Z dt ( p i ˙ q i − H ( p, q )) , (1.7)and its Poisson bracket [ F, G ] = ∂F∂q i ∂G∂p i − ∂F∂p i ∂G∂q i , (1.8)where F and G are general functions of the generalized coordinates q i and p i , for i =1 , , . . . , N . The system’s equations of motion are˙ q i = ∂H∂p i = [ q i , H ] , ˙ p i = − ∂H∂q i = [ p i , H ] , The time derivative of a function G ( p, q, t ) of the canonical variables, which couldexplicitly depend on time t , is expressed in terms of Poisson brackets as dG ( p, q, t ) dt = ∂G∂q i ˙ q i + ∂G∂p i ˙ p i + ∂G∂t , = ∂G∂q i ∂H∂p i − ∂G∂p i ∂H∂q i + ∂G∂t , = [ G, H ] + ∂G∂t . (1.9)9HAPTER 1. HAMILTONIAN MECHANICSLet us assume that we are interested in a quantity Q ( q, p, t ) which is conserved under thedynamics. Then Equation (1.9) indicates that dQ ( p, q, t ) dt = 0 , ⇒ [ Q, H ] + ∂Q∂t = 0 . (1.10)In many cases of practical interest, the quantity Q has no explicit time dependence and thenEq. (1.10) reduces to having the Poisson bracket of Q with the Hamiltonian H equal to zero,i.e. dQdt = 0 ⇒ [ Q, H ] = 0. For example in the case of the free motion of a unit mass particleon the plane, which is described by the Hamiltonian H = 12 ( p x + p y ), the angular momentum Q = xp y − yp x is a constant of motion as [ H, Q ] = − p x p y + p y p x = 0.If the Hamiltonian function (1.4) of a system does not explicitly depend on time thenit is an integral of motion, i.e. its value is conserved H = const , because from Eq. (1.10) we see that dHdt = 0. This behavior can be interpreted as follows:under stable conditions, performing an experiment today or tomorrow one expects to get thesame results. In other words, for autonomous dynamical systems (i.e. systems which do notexplicitly depend on time t ) the system’s evolution is independent of the choice of the initialtime value t . Consider a continuous autonomous Hamiltonian system (1.6) with N dof described by 2 N variables p i , q i , i = 1 , . . . , N . The time evolution of an orbit of this system is governed byHamilton equations of motion˙ x = f ( x ) = (cid:20) ∂H∂p − ∂H∂q (cid:21) T = J N · D H , (1.11)where J N = (cid:20) N I N − I N N (cid:21) , with I N and 0 N being the identity and zero N × N matrices, respectively, and D H = (cid:20) ∂H∂q ∂H∂q . . . ∂H∂q N ∂H∂p ∂H∂p . . . ∂H∂p N (cid:21) . The time evolution of an initial deviation vector at time t , v ( t ) = δx ( t ) =( δq ( t ) , δp ( t )) = ( δq ( t ) , . . . , δq N ( t ) , δp ( t ) , . . . , δp N ( t )), from a given orbit with ICs x ( t ) = ( q ( t ) , p ( t )) is defined by the so-called variational equations [13]:˙ v ( t ) = (cid:2) J N D H ( x ( t )) (cid:3) · v ( t ) , (1.12)where D H ( z ( t )) is the 2 N × N Hessian matrix with elements D H ( x ( t )) i,j = ∂ H∂x i ∂x j (cid:12)(cid:12)(cid:12) x ( t ) , (1.13)10HAPTER 1. HAMILTONIAN MECHANICSwith i, j = 1 , , . . . , N . Eq. (1.12) can be considered as the Hamilton equations of motionof the so-called tangent dynamics Hamiltonian H v of (1.6) H v ( δq, δp ) = 12 N X i =1 δp i + 12 N X i,j =1 D H ( q ( t )) ij δq i δq j . (1.14)Let us now consider the 2 N D area-preserving symplectic mapping (1.1) of discretetime. The evolution of x deviation vector v n at time t = n , n ∈ Z related to a reference orbit x n , is the so-called the tangent map v n +1 = M n · w n , (1.15)where M n = ∂f ( x n +1 ) ∂x n +1 = ∂f ∂x ∂f ∂x . . . ∂f ∂x N ∂f ∂x ∂f ∂x . . . ∂f ∂x N ... ... ... ∂f N ∂x ∂f N ∂x . . . ∂f N ∂x N , (1.16)is the system’s Jacobian matrix.The evolution of the deviation vectors is needed for the computation of the chaosindicators. For a continuous Hamiltonian system, this evolution is done along thesimultaneous integration of the equations of motion. In the case of symplectic mappings,the evolution of the deviation vectors is computed by simultaneously iterating the mapping(1.1) and the tangent map (1.15). Any mapping f ∈ ( R N , R N ) is area-preserving if ∀ M ∈ R N , d ( f ( M )) = d ( M ), where d ( M ) is the N Dmeasure of M . hapter 2Numerical techniques and models Symplectic integrators (SIs) are numerical schemes aiming to determine the solution ofHamilton equations of motion, preserving at the same time the Hamiltonian system’sunderlying symplectic structure [34–38]. An advantage of SIs is that their applicationtransforms the numerical integration of the Hamilton equations of motion into the applicationof a symplectic mapping. The use of SIs is facilitated by the presence of a separableHamiltonian function, i.e. when the whole Hamiltonian system can be written as a sumof Hamiltonian terms, whose solution is known explicitly.Let us present a general way of constructing explicit SIs for the separable Hamiltonian H ( q, p ) = T ( p ) + V ( q ) , (2.1)where T ( p ) is the system’s kinetic energy and V ( q ) is the potential energy, following anapproach based on Lie algebraic notion, Eq. (1.11) can be simply expressed as (see e.g. [35]) dxdt = L H x, (2.2)where L H is a differential operator defined by the Poisson bracket L H f = [ f, H ]. The solutionof this set of equations, for ICs x (0)= x is formally given as x = e tL H x . (2.3)In the common case of Eq. (2.1) the Hamiltonian function can be split into twointegrable parts as H ( q, p ) = A ( p ) + B ( q ), with A ( p ) being the kinetic energy T ( p ), which isa function of only the momenta p i , and B ( q ) being the potential energy V ( q ) depending onlyon the coordinates q i . A symplectic scheme for integrating (2.2) from time t to t + τ , with τ being the integration step, consists of approximating the operator e τL H by an integratorof j steps involving products of operators e c i τL A and e d i τL B , i = 1 , , ..., j , which are exactintegrations over times c i τ and d i τ of the integrable Hamiltonian functions A and B . c i and d i being carefully chosen constants in order to improve the accuracy of the integrationscheme. Thus, a SI approximates the action of the operator of e τL H by a product of the form e τL H = j Y i =1 e ( c i τ ) L A e ( d i τ ) L B + O ( τ j +1 ) , (2.4)12HAPTER 2. NUMERICAL TECHNIQUES AND MODELSwhere c i and d i are constants such that P ji =1 c i = P ji =1 d i = 1 and j ∈ N is the so-called orderof the integrator . Each operator e c i τL A and e d i τL B corresponds to a symplectic mapping, andconsequently the product appearing on the right-hand side of Eq. (2.4) is also a symplecticmapping. Various approaches have been developed over the years in order to determine thevalues of the coefficients c i and d i resulting to schemes of different orders (see e.g. [34, 36–42]and references therein). A basic second order SI can be written in the form S nd ( τ ) = e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A , (2.5)like for example: the so-called leapfrog method having three steps (i.e. number of applicationsof the simple operator e c i τL A and e d i τL A ) with constants c = c = 0 . d = 1. Thus e τL H = e ( τ/ L A e τL B e ( τ/ L A . (2.6)Let us now present some other 2nd order SIs like the 5 step SABA2 SI with compositionconstants c = − √ , c = √ and d = e τL H = e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A , (2.7)and the SBAB2 scheme with constants c = , d = and d = [42] e τL H = e ( d τ ) L A e ( c τ ) L B e ( d τ ) L A e ( c τ ) L B e ( d τ ) L A . (2.8)Another SI of this kind is the 9 step ABA82 scheme [36] e τL H = e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A , (2.9)where c = − √ √ , c = √ √ − √ − √ , c = √ − √ ,d = − √ , d = + √ . A 4th order SI can be obtained by a symmetric repetition (product) of 2nd order SIs (2.5)in the form S th ( τ ) = S nd ( r τ ) S nd ( r τ ) S nd ( r τ ) , (2.10)where r and r are two real adequately determined constants. The construction can lead tothe integrator developed by Forest and Ruth [39] involving 7 step e τL H = e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A e ( d τ ) L B e ( c τ ) L A , (2.11)with composition constants c = c = − / ) , c = c = − / − / ) ,d = d = − / , d = − / − / . c i , d i , i = 1 , ,
3. Higher order coefficients for thisfamily of SIs is found in [42].Quite efficient fourth order schemes named ABA864 and ABAH864 having respectively15 and 17 steps were developed in [36]. The corresponding coefficients c i , d i , i = 1 , , , Once a 4th order SI is found, it is easy to obtain a 6th order scheme using the 4th order andimplementing the same composition process, i.e. S th ( τ ) = S nd ( e τ ) S nd ( e τ ) S nd ( e τ ) . (2.12)An example of this construction is the SI developed in [34] which has 19 steps.More generally, if a SI of order 2 n , S n , is already known, a SI of order (2 n + 2) can beobtained through the composition S n +2 ( τ ) = S n ( z τ ) S n ( z τ ) S n ( z τ ) , (2.13)with z = − / (2 n +1) − / n +1) , and z = 12 − / n +1) . An efficient numerical technique for visualizing the behavior of a dynamical system is theso-called PSS (see for e.g. [6, 29]) which is named after Henri Poincar´e. According to thismethod the dynamics in the phase space of a high-dimensional system is understood byobserving the behavior induced by the flow on a particular section of the phase space. Inparticular, the dynamics is represented by the successive intersections of orbits with the PSSwhen this section is crossed in the same direction. In this way, a mapping is defined (seeFig. 2.1). The created Poincar´e mapping is a discrete dynamical system which represents thecontinuous flow of the original dynamical system.14HAPTER 2. NUMERICAL TECHNIQUES AND MODELSFigure 2.1: The Poincar´e mapping P evolves point x onto point P ( x ) (taken from [43]).Let us discuss this in more detail by considering an N D dynamical system defined bythe set of ordinary differential equations ˙ x = f ( x ) . (2.14)Let S be an ( N − S should cross it and should not evolve parallel to it. The Poincar´emapping is a mapping P : S → S , which is obtained by following trajectories from oneintersection after the other. Let us denote by x k the k th intersection of an orbit with thePSS, k ∈ N , and define the Poincar´e mapping as x k +1 = P ( x k ) . (2.15)If x ∗ is a fixed point of the mapping, the trajectory starting at this point comes backafter some iterations m (equivalently after some time T of the original dynamical system).Indicating that this is a T periodic orbit of the original system (2.14). This fixed pointcorresponds to an m -periodic orbit of the Poincar´e mapping (2.15).An efficient numerical approach to determine the Poincar´e mapping of a dynamicalsystem was proposed by H´enon [44]. Let us discuss this approach in more detail. Given an N D autonomous dynamical system of the form of Eq. (2.14), dx dt = f ( x , . . . , x N ) ,dx dt = f ( x , . . . , x N ) , ... dx N dt = f N ( x , . . . , x N ) , (2.16)the PSS S can be defined by a function of the form S ( x , x , . . . , x N ) = 0 . (2.17)15HAPTER 2. NUMERICAL TECHNIQUES AND MODELSThen (2.17) implicitly defines the Poincar´e mapping of dimension ( N − x N − a = 0 where a is a constant number. In order to findthe crossing points of a trajectory with the PSS defined by x N = a , H´enon introduced in [44]a simple idea using a system of differential equation which is equivalent to the system (2.17)by using x N as the independent integration variable instead of time tdtdx = 1 f ,dx dx = f f , ... dx N dx = f N f . (2.18) The LEs were introduced by Lyapunov [9] and they measure the exponential divergenceof nearby orbits in phase space. LEs are quantities which measure the system’s sensitivedependence on ICs, i.e. the rate at which the information is lost as the system evolves intime. The maximum Lyapunov Exponent (mLE), σ (see [13] and references therein), is theaverage rate of divergence (or convergence) of two neighboring trajectories in the phase spaceof a dynamical system and it is given by σ = lim t →∞ t ln k d x Φ t v k , (2.19)where the operator d x Φ t maps the deviation vector v from the tangent space at one point ofthe trajectory to the tangent space of the next point along the orbit. Practically the mLEcan be computed as the limit for t −→ ∞ of the quantity χ ( t ) = 1 t ln k v ( t ) kk v (0) k , (2.20)which is usually called finite-time mLE [13]. In (2.20) v (0) and v ( t ) are deviation vectorsfrom a given orbit, at times t = 0 and t >
0, respectively, and k . k denotes any norm of avector. Taking the limit as time goes to infinity, we have σ = lim t →∞ χ ( t ) . (2.21)We note that the direct application of Eq. (2.20) faces some computational problems [11]as in general, it results in an exponential increase of k v ( t ) k in the case of chaotic orbit.More specifically, with increasing t this norm rapidly exceeds the possibilities of ordinarynumerical computations. Luckily, this can be avoided by taking into account the linearity ofthe evolution operator d x Φ t and the composition law of these mappings [45]. So the quantity σ can be computed as σ = lim t →∞ kτ k X i =1 ln( α i ) . (2.22)16HAPTER 2. NUMERICAL TECHNIQUES AND MODELSwhere α i = k v ( iτ ) kk v (0) k , τ is a small time interval of the integration where time t = kτ , k =1 , , . . . . Then after time τ the deviation vector is renormalized in order to avoid overflowproblems. The mLE, in the case of autonomous Hamiltonian systems, is positive for chaoticorbits whereas for regular orbits it goes to zero by following a power law σ ∝ t − [11]. The spectrum of LEs is another useful tool for estimating the stability and chaos of 2 N Ddynamical systems. While the information from the mLE σ can be used to determine theregular ( σ = 0) or chaotic ( σ >
0) nature of orbits, the knowledge of part, or of the wholeset of LEs, σ , σ , . . . , σ N , provides additional information on the underlying dynamics andthe statistical properties of the system as they can be used for example to measure the fractaldimension of strange attractors in dissipative systems [13]. Following [13] the ordering of thespectrum of LEs is given by σ ( x ) ≥ σ ( x ) ≥ . . . ≥ σ N ( x ) ≥ − σ N ( x ) ≥ . . . ≥ − σ ( x ) ≥ − σ ( x ) . (2.23)The sum s of all LEs (i.e. N P i =1 σ i ) measures the contraction rate of volumes in the phasespace. In the so-called dissipative systems, s <
0, meaning that volumes visited by generictrajectories shrink exponentially. In autonomous Hamiltonian systems, s = 0, i.e. volumesare preserved. Similarly the sum of LEs for area-preserving symplectic mappings is zero N X i =1 σ i ( x ) = 0 . (2.24)The LEs come in pairs of values having opposite signs σ i ( x ) = − σ N − i +1 ( x ) , i = 1 , . . . , N. (2.25)In addition, at least one pair of LEs is by default zero σ N ( x ) = − σ N +1 ( x ) = 0 . (2.26)Other vanishing exponents may signal the existence of additional constants of motion apartfrom the Hamiltonian itself.For 2 N D Hamiltonian systems, the computation of the first p LEs, with 1 < p ≤ N ofan orbit is performed using the so-called standard method in [46]. This method involves thetime evolution of p initial linearly independent and orthonormal deviation vectors, with anew set of orthonormal vectors, obtained by the Gram-Schmidt orthonormalization process,replacing the evolved deviation vectors, along with the simultaneous integration of theequations of motions. 17HAPTER 2. NUMERICAL TECHNIQUES AND MODELS The idea behind the introduction of a simple, fast and efficient chaos indicator, such as theSALI [19, 21] was the need to overcome the slow limit convergence of σ ( t ) (2.21). Instead ofestimating the average rate of exponential growth, i.e. σ ( t ), SALI uses the possible alignmentof any two normalized deviation vectors to identify the chaotic nature of orbits. In order tocompute the SALI, we follow the evolution of the orbit and two deviation vectors, v ( t ) and v ( t ). The SALI is computed at any unit time by [19]SALI( t ) = min n k ˆ v ( t ) + ˆ v ( t ) k , k ˆ v ( t ) − ˆ v ( t ) k o , (2.27)where t is either continuous or discrete time and ˆ v and ˆ v are unit vectors given by therelation ˆ v i = v i k v i k , i = 1 , . (2.28)Note that 0 ≤ SALI( t ) ≤ √
2. SALI = 0 means the two deviation vectors, v ( t ) and v ( t )align (i.e. both being either parallel or antiparallel) and SALI = √ The asymptotic behavior of the SALI for regular and chaotic motion is given as follows:1. In the case of regular motion SALI does not become zero [47]. The two deviationvectors, v and v , tend to fall on the tangent space of the torus, following a t − timeevolution and having in general two different directions. Thus, in this case, SALI attainsa constant positive value, i.e. SALI ∝ const . (2.29)2. On other hand, for chaotic orbits, the deviation vectors align in the direction definedby the mLE and the SALI tends exponentially fast to zero at a rate which is related tothe difference of the two largest LEs, σ and σ as discussed in [48], i.e.SALI( t ) ∝ e − ( σ − σ ) t . (2.30)Rather than evaluating SALI (2.27), in order to see if the two deviation vectors arealigned or not we can use the wedge product of these vectors [21] k ˆ v ∧ ˆ v k = k ˆ v + ˆ v k · k ˆ v − ˆ v k , (2.31)which represents the area of the parallelogram formed by the two deviation vectors. Then k ˆ v ∧ ˆ v k → k ˆ v ∧ ˆ v k 6 = 0 corresponds to regular motion. If we takemore than two deviation vectors, say v , v , . . . , v k , 2 ≤ k ≤ N , then the wedge product ofthe corresponding unit vectors evaluates the volume of the parallelepiped formed by thesedeviation vectors, which leads to the introduction of the GALI Method.18HAPTER 2. NUMERICAL TECHNIQUES AND MODELS Consider the 2 N D phase space of a conservative dynamical system represented by eithera Hamiltonian flow of N dof or a 2 N D symplectic mapping. In order to study whether anorbit is chaotic or not, we examine the asymptotic behavior of k initially linearly independentdeviations from this orbit, denoted by vectors v , v , . . . , v k , 2 ≤ k ≤ N . Thus, we follow theorbit, using the Hamilton equations of motion or the mapping equations and simultaneouslywe solve the corresponding variational equations or the related tangent map to study thebehavior of deviation vectors from this orbit.The Generalized Alignment Index of order k (GALI k ) represents the volume of thegeneralized parallelogram defined by the evolved k unit deviation vectors at any given timeand it is determined as the norm of the wedge (or exterior) product of these vectors [18]GALI k ( t ) = k ˆ v ∧ ˆ v ∧ ˆ v ∧ · · · ∧ ˆ v k k . (2.32)In the case that the number of deviation vectors k exceeds the dimension 2 N ofthe system’s phase space, by definition the vectors will be linearly dependent and thecorresponding volume (the value of GALI k , for k > N ) will be zero. Let us consider N dof Hamiltonian systems. For regular orbits, if we start with k ≤ N linearly independent initial deviation vectors, then the deviation vectors eventually fall onthe N D tangent space of the torus [18]. In this case, the asymptotic GALI value will bepractically constant. Whereas, if we start with
N < k ≤ N linearly independent initialdeviation vectors, then the asymptotic GALI value will be zero since some deviation vectorswill eventually become linearly dependent as again all of them will fall on the N D tangentspace of the torus. The general behavior of the GALI k for regular orbits lying on an N Dtorus is given by [18] GALI k ( t ) ∝ constant if 2 ≤ k ≤ N,t − k − N ) if N < k ≤ N. (2.33)We note that, the behavior of GALI k for regular orbits on lower-dimensional tori i.e. k Dtorus with k < N is given by [28].GALI k ( t ) ∝ t − ( k − if 2 ≤ k ≤ N − ,t − N if k = 2 N. (2.34)On the other hand, for chaotic and unstable periodic orbits all deviation vectors alignin the direction defined by the mLE and the value of GALI k decays exponentially fast to zerofollowing a rate which depends on the values of several LEs [18, 28]GALI k ( t ) ∝ e − [( σ − σ )+( σ − σ )+ ··· +( σ − σ k )] t , (2.35)where σ ≥ σ ≥ · · · ≥ σ k are approximations of the k largest LEs of the orbit. In particularwhen k = 2, GALI tends to zero exponentially fast as follows:GALI ( t ) ∝ e − ( σ − σ ) t , (2.36)19HAPTER 2. NUMERICAL TECHNIQUES AND MODELSwhich is similar to the behavior of the SALI shown in (2.30). This is due to the fact that theSALI is equivalent to the GALI [18]. The behavior of the GALI for regular orbits in 2Dmappings requires special attention. Since in this case the motion exists on a 1D torus, thetwo deviation vectors tend to fall on the tangent space of this torus, which is again 1D. Thusthe two vectors will eventually become linearly dependent as in the case of chaotic orbits,but this happens with a different time rate given by [21]GALI ( n ) ∝ n . (2.37)where n is the number of iterations.An efficient way of computing the value of GALI k is through the Singular ValueDecomposition (SVD) of the matrix A = (cid:0) ˆ v ˆ v . . . ˆ v k (cid:1) = ˆ v , ˆ v , . . . ˆ v k, ˆ v , ˆ v , . . . ˆ v k, ... ... ...ˆ v , N ˆ v , N . . . ˆ v k, N , (2.38)having as columns the coordinates of the k unitary vectors ˆ v i ( t ) = v i ( t ) k v i ( t ) k = (ˆ v i, , ˆ v i, , . . . , ˆ v i, N )[20]. Following [20] we see that the value of GALI k can be given byGALI k ( t ) = p det [ A ( t ) · A T ( t )] . (2.39)where ‘ det ( A )’ denotes the determinant of a matrix A . The matrix product in (2.39) is a k × k symmetric positive definite matrix AA T = < ˆ v , ˆ v > < ˆ v , ˆ v > < ˆ v , ˆ v > . . . < ˆ v , ˆ v k >< ˆ v , ˆ v > < ˆ v , ˆ v > < ˆ v , ˆ v > . . . < ˆ v , ˆ v k > ... ... ... ... < ˆ v , ˆ v k > < ˆ v , ˆ v k > < ˆ v , ˆ v k > . . . < ˆ v k , ˆ v k > , (2.40)where each element of i th row and j th column is the inner product of the unit deviationvectors ˆ v i and ˆ v j , so that < ˆ v i , ˆ v j > = cos θ ij , i, j = 1 , , . . . , k, (2.41)where θ ij is the angle between vectors ˆ v i and ˆ v j . Thus, the matrix of Eq. (2.40) can bewritten as AA T = θ cos θ . . . cos θ k cos θ θ . . . cos θ k ... ... ... ...cos θ k cos θ k cos θ k . . . . (2.42)The value of the GALI k can be computed through the norm of the wedge product ofthe k deviation vectors defined in (2.32) asGALI k = X ≤ i
Let us start our presentation by considering a simple area-preserving 2D symplectic mappingthe so-called standard mapping [46]: x ′ = x + x ,x ′ = x − K sin( x + x ) , (2.46)where K is a real parameter and x ′ and x ′ denote the values of coordinates x x after oneiteration of the motion. We note that all coordinates have (mod 1), i.e. 0 ≤ x i < i = 1 , K are shown in Figure 2.2. More specifically, Fig. 2.2a shows that the phase space is practicallyfilled with stability islands for a small value of the parameter K = 0 .
5. When we increasethe value K to K = 1 . K = 5, the chaotic sea dominates as we seein Fig. 2.2c. (a) x x (b) x x (c) x x Figure 2.2: Phase space portraits of the standard mapping (2.46) for (a) K = 0 .
5, (b) K = 1 .
2, and (c) K = 5. The portraits are produced by 10 iterations of several ICs.Different colors are used to help distinguishing between regular (blue points) and chaotic(red points) orbits.The tangent map (1.15) of the 2D mapping (2.46) which is needed in order to compute22HAPTER 2. NUMERICAL TECHNIQUES AND MODELSchaos indicators is given by δx ′ = δx + δx ,δx ′ = δx − K ( δx + δx ) cos( x + x ) . (2.47)We consider several orbits in the 2D mapping (2.46) for K = 1 . . The obtainedresults are shown in Fig. 2.3. In order to to illustrate the behavior of the chaos indicators, weconsider three different regular and chaotic orbits from Fig. 2.2b. The phase space portraitsof the standard mapping (2.46) for these orbits is shown in Fig. 2.3a.Fig. 2.3b shows that the mLE goes to zero for the regular orbits following an evolutionwhich is proportional to the power law n − (which is indicated by the dashed straight line)whereas mLE eventually saturates to a constant value for chaotic orbits. Fig. 2.3c illustratesthe time evolution of the SALI for the same regular and chaotic orbits. In the case of regularorbits, the evolution of the SALI is proportional to the theoretical prediction n − (2.37),while for chaotic orbits the SALI goes to zero exponentially fast following the exponentialdecay in (2.30). If we compare the time evolution of the mLE and SALI for chaotic orbits ofFigs. 2.3b and 2.3c we can notice that SALI discriminates the chaotic orbits very fast.23HAPTER 2. NUMERICAL TECHNIQUES AND MODELS [b][b][y] [y] [y][g][g] [p][k] [r] (a) x x -6-5-4-3-2-1 0 0 1 2 3 4 5 6 [b] [y] [g][p] [k] [r] (b) l og χ ( n ) log n Slope = -1 -12-9-6 -3 0
0 1 2 3 4 5 6 [b][y] [g][p][k] [r] (c) l og [ S A L I( n ) / S A L I( )] log n Slope = -2
Figure 2.3: (a) Phase space portraits of the standard mapping (2.46) for three differentregular orbits with ICs (0 . , .
12) (blue ‘[b]’), (0 . , .
22) (green ‘[g]’) and (0 . , . . , .
4) (black ‘[k]’), (0 . , .
6) (red ‘[r]’)and (0 . , .
7) (purple ‘[p]’) for K = 1 .
2. We see the evolution of the mLEs and the SALIsof these orbits, with respect to the number of iteration n , respectively in panels (b) and (c).In (b) σ ( n ) tends to zero following a power law n − for regular orbits and becomes constantfor chaotic ones. In (c) SALI( n ) tends to zero following a power law n − for regular orbitsand it goes to zero exponentially fast for chaotic orbits. The dotted straight lines in (b) and(c) correspond respectively to functions proportional to n − and n − .Similar behaviors are also observed in the case of the 4D mapping [45, 50, 51] x ′ = x + x ,x ′ = x − K sin( x + x ) − µ [1 − cos( x + x + x + x )] ,x ′ = x + x ,x ′ = x − K sin( x + x ) − µ [1 − cos( x + x + x + x )] , (2.48)Consisting of two coupled 2D standard mappings with parameters K and K , which arecoupled through a term whose strength is defined by the parameter µ . In this case allcoordinates x i , i = 1 , . . . , π ), i.e. − π ≤ x i < π .Figures 2.4 to 2.6 show 2D projection of the 4D mapping for different values of µ byfixing the parameters K = 0 . K = 0 .
1. In each case we iterated the mapping (2.48)by using five different ICs of the form ( x , x , x , x ) = ( c i , , c i , i = 1 , . . . ,
5, where c = 0 . c = 2, c = 2, c = 2 . c = 3. From the results of Fig. 2.4 we see that theconsidered orbits are mainly regular for µ = 0 except for the orbit with x = x = 3 which24HAPTER 2. NUMERICAL TECHNIQUES AND MODELSshows a weak chaotic behavior. As the value of the coupling parameter increases to µ = 10 − (Fig. 2.5) and µ = 10 − (Fig. 2.6) more orbits become chaotic. -1 0 1 -3 -2 -1 0 1 2 3 [b] [g][y][k] [r] (a) x x [b] [g][y][k] [r] (b) x x Figure 2.4: Two-dimensional projection of orbits of the 4D mapping (2.48) for µ = 0 withICs (0 . , , . ,
0) (blue ‘[b]’), (1 , , ,
0) (green ‘[g]’), (2 , , ,
0) (yellow ‘[y]’), (2 . , , . , , , , . ) (red ‘[r]’) on (a) the x x plane and (b) the x x plane. -1 0 1 -3 -2 -1 0 1 2 3 [b] [g][y][k] [r] (a) x x [b] [g][y][k] [r] (b) x x Figure 2.5: Similar to Fig. 2.4 but for µ = 10 − .25HAPTER 2. NUMERICAL TECHNIQUES AND MODELS -1 0 1 -3 -2 -1 0 1 2 3 [b] [g][y][k] [r] (a) x x [b] [g][y][k] [r] (b) x x Figure 2.6: Similar to Fig. 2.4 but for µ = 10 − .The tangent map of system (2.48) is given by δx ′ = δx + δx ,δx ′ = δx − K ( δx + δx ) cos( x + x ) − µ ( δx + δx + δx + δx ) sin( x + x + x + x ) ,δx ′ = δx + δx ,δx ′ = δx − K ( δx + δx ) cos( x + x ) − µ ( δx + δx + δx + δx ) sin( x + x + x + x ) , (2.49)Let us consider the orbits of the 4D mapping (2.48) with ICs x = 0 . , x = 0 , x = 0 . x = 0 as a representative of a regular (blue points in Fig. 2.6) and x = 3 , x = 0 , x = 3and x = 0 for a chaotic orbit (red points in Fig. 2.6). In Fig. 2.7a we plot the four LEs forthe regular and the chaotic orbit. For regular orbit, the LEs go to zero following an evolutionwhich is proportional to the power law n − . For the chaotic orbit, the LEs come in pairs ofopposite values, i.e. χ = − χ and χ = − χ . For this reason, we plot χ , χ , | χ | , | χ | . Thefour finite time LEs eventually saturate to a constant value for the chaotic orbit. In Fig. 2.7bwe see that the time evolution of the SALI for the regular orbit leads to a constant valuewhereas the SALI goes to zero exponentially fast for the chaotic orbit.26HAPTER 2. NUMERICAL TECHNIQUES AND MODELS -6-5-4-3-2-1 0 1 0 1 2 3 4 5 6 [b][r] (a) l og χ ( n ) , χ ( n ) , χ ( n ) , χ ( n ) log n Slope = -1 -12-9-6-3 0
0 1 2 3 4 5 6 [b][r] (b) l og [ S A L I( n ) / S A L I( )] log n Figure 2.7: Evolution of (a) the four LEs and (b) the SALIs with respect to the numberof iteration n of the 4D mapping (2.48) with µ = 10 − , K = 0 . K = 0 . x = 0 . , x = 0 , x = 0 . x = 0 (blue ‘[b]’ curves) and thechaotic orbit with IC, x = 3 , x = 0 , x = 3 and x = 0 (red ‘[r]’ curves). In (a) the LEsof regular orbit tend to zero following a power law n − while they attain a constant positivevalue for chaotic orbits. We note that for the chaotic orbit we plot the quantities χ , χ , | χ | , | χ | in order to avoid the logarithm of negative values. In (b) SALI( n ) reach a non-negativeconstant value for the regular orbit and it goes to zero exponentially fast for the chaotic one.The dashed straight line in (a) corresponds to a function proportional to n − . Axes in bothpanels are in logarithmic scale.Furthermore, in Fig. 2.8 we computed the evolution of the GALIs for these two orbitsof Fig. 2.7. In Fig. 2.8a we plot the evolution of GALI , GALI and GALI for the regularorbit. We see that GALI eventually saturates to a non-zero constant value while GALI and GALI go to zero following respectively the asymptotic power law decays t − and t − .In Fig. 2.8b we see the evolution of GALI , GALI and GALI for the chaotic orbit. In thiscase, the GALIs go to zero following an exponential decay which is proportional to exp[ − σ n ],exp[ − σ n ] and exp[ − σ n ] for σ = 0 . -12-8-4 0
1 2 3 4 5 6 (a)
GALI GALI GALI l og [ GA L I k ( n ) / GA L I k ( )] log n Slope = -2 Slope = -4 -12-8-4 0
0 200 400 600 800 (b)
GALI GALI GALI n Slope = - σ /ln(10) Slope = -2 σ /ln(10)Slope = -4 σ /ln(10) Figure 2.8: Evolution of the GALI k ( n ), k = 2 , , n ofthe 4D mapping (2.48) with µ = 10 − , K = 0 . K = 0 . x = 0 . , x = 0 , x = 0 . x = 0 and (b) the chaotic orbit with IC x = 3 , x = 0 , x = 3and x = 0. The dashed straight lines in (a) correspond to functions proportional to n − and n − while in (b) they represented functions proportional to exp[ − σ n ], exp[ − σ n ] andexp[ − σ n ] for σ = 0 . The H´enon-Heiles system is a prototypical 2D Hamiltonian dynamical system initially usedin [29] to investigate the motion of a star in a simplified galactic potential. The star isassumed to move on the galactic plane with coordinates x and y . The Hamiltonian functionof this system is H = 12 ( p x + p y ) + 12 ( x + y ) + x y − y , (2.50)where p x and p y are the conjugate momentum. The system’s equations of motion are˙ x = p x , ˙ y = p y , ˙ p x = − x (1 + 2 y ) , ˙ p y = y − x − y. (2.51)In Fig. 2.9 we present the PSS of the H´enon-Heiles system (2.50) with H = 0 . x = 0 and p x ≥
0. There ordered orbits correspond to closedsmooth curves, while the chaotic orbits are represented by scattered points.28HAPTER 2. NUMERICAL TECHNIQUES AND MODELS -0.6-0.4-0.2 0 0.2 0.4 0.6-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 (a) P y y Figure 2.9: The PSS of the H´enon-Heiles system (2.50) with energy H = 0 .
125 for x = 0, p x ≥
0. Orbits were integrated up to a final time t = 10 ,
000 time units using the SABA2 SIscheme with an integration time step of τ = 0 . δ ˙ x = δp x ,δ ˙ y = δp y ,δ ˙ p x = − (cid:2) (1 + 2 y ) δx + 2 xδy (cid:3) ,δ ˙ p y = − (cid:2) xδx + (1 + 2 y ) δy (cid:3) . (2.52)In our simulations we simultaneously integrated Eqs. (2.51) and (2.52) using the SABA2 SIscheme with an appropriate integration time step which keep the relative energy error below10 − . We compute the mLEs and the SALIs for three different regular orbits which are takenfrom the islands of stability around periodic orbits of different periods (Fig. 2.10a). TheirICs are (0 , − . , . , − . , . , − . , , . , . , − .
25) for the regular orbitsand (0 , . , . , . , − . , . , − . , − . , . , − .
03) for the chaoticones. Fig. 2.10b shows that the mLEs go to zero for regular orbits following an evolutionwhich is proportional to the power law t − , while they eventually attain constant positivevalues for chaotic orbits. In Fig. 2.10c we present the time evolution of the SALIs for the sameregular and chaotic orbits. In the case of regular orbits, the SALIs obtain non-zero constantvalues. On the other hand, the SALIs go to zero exponentially fast for chaotic orbits. Similarto the discrete system of Section 2.3.1.1 (Figs. 2.3 and 2.7). The SALIs discriminate chaoticorbits very fast. An orbit x ( t ), in the system’s phase space, is called periodic if x ( t + T ) = x ( T ), for a non-zero constantnumber T -0.6-0.4-0.2 0 0.2 0.4 0.6-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 [b] [o][g][p] [k] [r] (a) P y y -5-4-3-2-1 0 0 1 2 3 4 5 6 [b][o][g] [p][k] [r] (b) l og χ ( t ) log t Slope = -1 -12-9-6 -3 0
0 1 2 3 4 5 6 [b] [o][g] [p][k] [r] (c) l og [ S A L I( t ) / S A L I( )] log t Figure 2.10: (a) The PSS ( x = 0 , p x ≥
0) of the 2D system (2.50) with energy H = 0 . y, p x , p y ) ( − . , . , − .
41) ‘blue [b]’, (0 . , − . ,
0) ‘green [g]’,and (0 . , . , − .
25) ‘orange [o]’ and chaotic orbits with ICs (0 . , . , .
14) ‘black [k]’,( − . , . , − .
2) ‘red [r]’, and ( − . , . , − .
03) ‘purple [p]’. The time evolution of themLEs and the SALIs are respectively shown in (b) and (c). In (b) the plotted line correspondsto a function proportional to t − .In Fig. 2.11 we also study the behavior of the GALIs for the regular orbit with IC x =0 ., y = 0 . , p x = 0 .
495 and p y = 0, and the chaotic orbit with IC x = 0 , p x = − . , x = 0 . x = 0. The regular orbit lies on a 2D torus and therefore its GALI eventually saturatesto a non-zero constant value, whereas the GALI and the GALI go to zero respectivelyfollowing the asymptotic power law decays t − and t − (Fig. 2.11a). The GALIs for thechaotic orbit are shown in Fig. 2.11b. All of them go to zero exponentially fast following theexponential decays which are described in Eq. (2.35).30HAPTER 2. NUMERICAL TECHNIQUES AND MODELS -10-5 0
1 2 3 4 5 6 (a)
GALI GALI GALI l og [ GA L I k ( t ) / GA L I k ( )] log t Slope = -2 Slope = -4 -10-5 0
150 300 450 (b)
GALI GALI GALI t Slope = - σ /ln(10) Slope = -2 σ /ln(10)Slope = -4 σ /ln(10) Figure 2.11: The time evolution of GALI k ( t ), k = 2 , , x =0 , y = 0 . , p x = 0 .
495 and p y = 0 and (b) the chaotic orbit with IC x = 0 , p x = − . , x =0 .
421 and x = 0 for the system (2.50) with energy H = 0 . t − and t − while in (b) they correspond tofunctions proportional to exp[( − σ ) t ], exp[( − σ ) t ] and exp[( − σ ) t ] for σ = 0 . x -axis is linearwhile the y -axis is in logarithmic scale.Overall all the methods we implemented can correctly capture the dynamical natureof the studied orbits. GALI k , with higher-order k , is the fastest in discriminating betweenregular and chaotic motion whereas the mLE is the least efficient technique. For example,consider the evolution of the GALI and the mLE for the 4D mapping (2.48) and 2DHamiltonian (2.51). The exponential decay of the GALI for the chaotic orbit with IC x = 3 , x = 0 , x = 3 and x = 0 (black curve in 2.8b) of (2.48) reaches the value 10 − around time t = 10 . , while the mLE (red curve in 2.7a) clearly indicates that the orbit ischaotic around time t = 10 . . Also, the exponential decay of the GALI for the chaotic orbitwith IC x = 0 , p x = − . , x = 0 .
421 and x = 0 (black curve in 2.11b) of (2.51) reachesthe value 10 − around time t = 10 . , while the mLE (purple curve in 2.10b) shows that theorbit is chaotic around time 10 . Based on our results and discussion, we can conclude thatthe GALI method is the most effective among the methods we have considered. Even thoughthe convergence of the LEs is theoretically guaranteed, its rate of convergence can becometoo slow, while the PSS is impractical for systems with many dof. A dynamical system is called non-conservative if its phase space volume changes over time.The determinant of the system’s Jacobian matrix is an important tool to explore theproperties of the dynamic, in particular to describe the change in volume. If this determinantis smaller than one, then the volume of phase space decreases over time [52].For dissipative systems, in the limit of t → ∞ the dynamics converges to a set ofvalues for various ICs, the so-called attractor . A dynamical system can have several differentattractors. We should note that the concept of an attractor is only relevant when discussingthe asymptotic behavior of orbits.In this section, we compute various chaos indicator of two dissipative systems. The aimis to briefly investigate the behavior of the GALIs for these systems. Our results suggest31HAPTER 2. NUMERICAL TECHNIQUES AND MODELSthat further work is certainly required. However, this study provides a good starting pointfor further investigations. First we consider a mapping which was introduced by H´enon [53] based on concepts ofstretching and folding of areas in the phase space. The H´enon mapping takes a point ( x, y )in the plane and maps it to a new point ( x ′ , y ′ ) according to the rule x ′ = y + 1 − αx ,y ′ = βx. (2.53)According to [53], the dynamics converges to a strange attractor for parameter values α = 1 . β = 0 .
3. This attractor is seen in Fig. 2.12 -1.2-0.6 0 0.6 1.2-0.4 -0.2 0 0.2 0.4 x y Figure 2.12: Phase plot of the strange attractor of the H´enon mapping (2.53) with α = 1 . β = 0 .
3. An IC (0 . , .
5) is iterated 10 times in order to produce this plot.The rate of expansion of an area in the map’s phase space is given by the determinantof the Jacobian matrix of the mapping (2.53) J = (cid:12)(cid:12)(cid:12)(cid:12) det (cid:20) − αx n β (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) = | β | . (2.54)Thus for β = 1 the mapping is area preserving but when β <
1, as we have done for the case β = 0 .
3, the mapping is dissipative and phase space areas contract.There is an equally simple (as the H´enon mapping (2.53)) 2D quadratic mappingadmitting several attractors [54]. This mapping, has two parameters, α and β , and it isgiven by x ′ = 1 − αy + βx,y ′ = x. (2.55)32HAPTER 2. NUMERICAL TECHNIQUES AND MODELSNote that this mapping is conservative for values of the parameter β ∈ ( − . , . δx ′ = − βδx − αy δy,δy ′ = δx. (2.56)We compute the LEs as well as the SALI of the mapping (2.55) of a representativeregular and a chaotic orbit. For the regular orbit we use the parameters ( α, β ) = (0 . , .
1) andthe IC x = 0 . y = 0 .
5. In Fig. 2.13a we plot the evolution of the two LEs of this orbit. Wesee that the LEs eventually attain negative values for this regular orbit, i.e. χ = − . χ = − . -0.135-0.132-0.129 -0.126
1 2 3 4 5 6 (a) χ ( n ) , χ ( n ) log n 0 0.06 0.12 0.18
1 2 3 4 5 6 (b) l og [ S A L I( n ) / S A L I( )] log n Figure 2.13: The evolution of (a) the two LEs and (b) the SALI with respect to the number ofiteration n of the mapping (2.55) for the regular orbit with IC ( x, y ) = (0 . , . α = 0 . β = 0 . α, β ) = (0 . , .
6) for the mapping (2.55) andcompute the chaos indicators for the chaotic orbit with IC ( x, y ) = (0 . , . χ = 0 . χ = 0 . (a) χ ( n ) , χ ( n ) log n -12-9-6-3 0
0 1 2 3 4 (b) l og [ S A L I( n ) / S A L I( )] log n Figure 2.14: Similar to Fig. 2.13 but for a chaotic orbit with IC ( x, y ) = (0 . , . α = 0 . β = 0 . Similar behaviors are observed for a more complicated mapping model given by [55] x n +1 = − Dx n − + 2 x n " C + BS Q − exp ( − r n ) r n + Eξ i ,y n +1 = − Dy n − + 2 y n " C + BS Q − exp ( − r n ) r n + Eη i , (2.57)where r n = x n + y n , n = 0 , , , . . . and Q , Q , B ∈ R , 0 < D < ξ i , η i are randomly generated noise variables, within the interval ( − . , +0 . c i = cos 2 πQ i , S i = sin 2 πQ i , i = 1 ,
2. The corresponding tangent map is given by δx n +1 = − Dδx n − + 2 C δx n + BS Q " x n + δx n ) 1 − exp ( − r n − p n ) r n + 2 p n − x n − exp ( − r n ) r n ,δy n +1 = − Dδy n − + 2 C δy n + BS Q " y n + δy n ) 1 − exp ( − r n − p n ) r n + 2 p n − y n − exp ( − r n ) r n , (2.58)where p n = x n δx n + y n δy n . Mapping (2.57) is conservative when D = 1.We also studied the time evolution of the mLE and the GALI k , k = 2 , , χ = − . x , y , x , y ) =( √ . , , , √ .
1) when the mapping parameters are set to E = 1 . Q . σ = 0 . D = 0 . B = 0 . x , y , x , y ) = ( √ . , , , √ .
1) when the mapping parameters are set to E = 1 . Q . σ = 0 . D = 1 . B = 0 . χ = 0 . -0.06-0.04-0.02 0 0.02 0.04
1 2 3 4 5 6 χ ( n ) log n -4-2 02 3 4 5 6 (b) l og [ GA L I k ( n ) / GA L I k ( )] log n GALI-2GALI-3GALI-4
Figure 2.15: Evolution of (a) the mLE and (b) the GALI k ( n ), k = 2 , , n of the mapping (2.57) for a regular orbit with IC ( x , y , x , y ) =( √ . , , , √ .
1) when the mapping parameters are set to E = 1 . Q . σ = 0 . D = 0 . B = 0 .
1 2 3 4 5 6 (a) χ ( n ) log n -12-8-4 0 0 1 2 3 4 (b) l og [ GA L I k ( n ) / GA L I k ( )] log n GALI-2GALI-3GALI-4
Figure 2.16: Similar to Fig. 2.15 but for a chaotic orbit with IC ( x , y , x , y ) =( √ . , , , √ .
1) for D = 1 . hapter 3Behavior of the GALI for regularmotion in multidimensionalHamiltonian systems The FPUT model [22, 23] describes a multidimensional lattice system and it is related tothe famous paradox of the appearance of a recurrent behavior in a nonlinear system. It istraditionally called the ( α + β ) − FPUT model and represents a one-dimensional chain of N identical particles with nearest-neighbor interactions (see Fig. 3.1). . . . x j . . . x j + . . . N Figure 3.1: Schematic model of the FPUT model: identical masses that can move only inone spatial dimension are coupled by non-linear springs. x j is the displacement from theequilibrium position of the j th mass. The boundary conditions are fixed.The ( α + β ) − FPUT Hamiltonian is H N = 12 N X j =1 p j + N X j =0 (cid:20)
12 ( x j +1 − x j ) + 13 α ( x j +1 − x j ) + 14 β ( x j +1 − x j ) (cid:21) , (3.1)where x j is the displacement of the j th particle from its equilibrium position and p j is thecorresponding conjugate momentum. In our study we impose fixed boundary condition tosystem (3.1), i.e. x = x N +1 = p = p N +1 = 0. We also chose to work with the β − FPUTsystem (i.e. α = 0 in (3.1)) H N = 12 N X j =1 p j + N X j =0 (cid:20)
12 ( x j +1 − x j ) + 14 β ( x j +1 − x j ) (cid:21) , (3.2)because it does not allow particles to escape.36HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSIn order to use a SI scheme for the integration of the system’s equations of motionand variational equations (using the so-called tangent map method [35, 56, 57]) we split theHamiltonian in two integrable parts H N = A + B in which A is a function of the momenta p j and B is a function of the coordinates q j : A = 12 N X j =1 p j , and B = N X j =0 (cid:20)
12 ( x j +1 − x j ) + 14 β ( x j +1 − x j ) (cid:21) . (3.3)Then the corresponding equations of motion for the Hamiltonian functions A and B ,for 1 ≤ j ≤ N , are ˙ X = L AZ X : ˙ x j = p j ˙ p j = 0˙ δx j = δp j ˙ δp j = 0 , (3.4)and˙ X = L BZ X : ˙ x j = 0˙ p j = x j +1 − x j + x j − − β ( x j +1 − x j ) + β ( x j − x j − ) ˙ δx j = 0˙ δp j = (1 + 3 β ( x i − x i − ) ) δx i − − (2 + 3 β ( x i − x i − ) − β ( x i +1 − x i ) ) δx i +(1 + 3 β ( x i +1 − x i ) ) δx i +1 . (3.5)These two sets of differential equations can be easily integrated. The operators e τL A z and e τL B z which propagate the set of ICs x j , p j , δx j and δp j at time t to x ′ j , y ′ j , δx ′ j and δp ′ j time t + τ ( τ being the integration time step), for 1 ≤ j ≤ N , are e τLAz X = x ′ j = x j + τ p j p ′ j = p i δx ′ j = δx j + τ δp j δp ′ j = δp j , (3.6)and e τLBz X = x ′ j = x j p ′ j = p j + τ (cid:2) x j +1 − x j + x j − − β ( x j +1 − x j ) + β ( x j − x j − ) (cid:3) δx ′ j = δx j δp ′ j = δp j + τ (cid:2) (1 + 3 β ( x i − x i − ) ) δx i − − (2 + 3 β ( x i − x i − ) − β ( x i +1 − x i ) ) δx i +(1 + 3 β ( x i +1 − x i ) ) δx i +1 (cid:3) . . (3.7) Computational considerations
Here we discuss some practical aspects of our numerical simulation. In order to computethe time evolution of GALI k of a 2 N dof Hamiltonian system, we need k initially linearlyindependent deviation vectors. In most considered cases the coordinates of these vectors are37HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSnumbers taken from a uniform random distribution in the interval [ − . , . n v = 10 sets of initial vectors. Therandom choice of the initial vectors leads to different GALI k (0) values, i.e. the value ofthe GALI k at the beginning of the evolution. Thus, to fairly and adequately compare thebehavior of the indices for different initial sets of vectors we normalize the GALIs’ evolutionby registering the ratio GALI k ( t )/GALI k (0), i.e. we measure the change of the volume definedby the k deviation vectors with respect to the initially defined volume. Another option isto start the evolution of the dynamics by considering a set of k orthonormal vectors so thatGALI k (0) = 1. As we will see later on, both approaches lead to similar results, so in ourstudy we will follow the first procedure, unless otherwise specifically stated.In our investigation, we implement an efficient fourth-order symplectic integrationscheme the so-called ABA864 method, to integrate the equations of motion and thevariational equations of the β − FPUT system (3.2). In our numerical simulations, we typicallyintegrate the β − FPUT system (3.2) up to a final time t f of 10 time units. By adequatelyadjusting the used integration time step we always keep the absolute relative energy error below 10 − .The scope of our work is to quantitatively analyze the chaotic and regular behavior oforbits of the multidimensional β − FPUT model (3.2). In Chapter 2, we described in detailthe behavior of the mLE, and the spectrum of LEs, as well as of the SALI and the GALI for2D and 4D systems. These behaviors should hold for multidimensional systems. In order toillustrate the behavior of these indices let us take the β -FPUT model (3.2) with N = 8 dofand consider some typical regular and chaotic orbits. More specifically we consider the regularorbit R with total energy H = 0 . − x = x = − x = x = − x = x = 0 . x = 0, x = 0 . p = p = p = 0, − p = p = − p = p = − p = 0 .
08, and β = 1and the chaotic orbit C with total energy H = 10 .
495 and ICs x = x = 2, x = x = 1, x = x = 0 . x = x = 0 . p i = 0, i = 1 , , . . . ,
8, and β = 1 . t − whereas the σ ( t ) remainspractically constant for the chaotic orbit (red curve), which is exactly what we expect from achaotic trajectory. In Fig. 3.2b we see that SALI goes to zero exponentially fast following thelaw of (2.30), i.e. SALI( t ) ∝ exp [ − ( σ − σ ) t ] for σ = 0 .
170 and σ = 0 . The relative energy error, RE , is defined as RE = (cid:12)(cid:12)(cid:12) H c − H H (cid:12)(cid:12)(cid:12) , where H c is the computed Hamiltonianvalue at the current time unit and H is the initial Hamiltonian value. -6-5-4-3-2-1 0 0 1 2 3 4 5 6 7 [b] [r] (a) l og χ ( t ) log t Slope = -1 -12-9-6-3 0
0 1 2 3 4 5 6 7 [r] [b] (b) l og [ S A L I( t ) / S A L I( )] log t Slope = exp[- (σ −σ ) ] Figure 3.2: Time evolution of (a) the mLE and (b) the SALI of the regular orbit R (bluecurves denoted by [b]) and the chaotic orbit C (red curves denoted by [r]) of the β − FPUTsystem (3.2) for N = 8. The dotted curve indicates a function proportional to (a) t − (and(b) e − ( σ − σ ) t for σ = 0 .
170 and σ = 0 . R and C . In the case of the regular orbit (Fig. 3.3a) all 16 LEs approach zero followingan evolution proportional to t − , while in the case of the chaotic orbit (Fig. 3.3b) the setof LEs consist of pairs of values having opposite signs. Positive LEs are plotted in greenwhile negative values in red. Moreover, two of the exponents, χ and χ (yellow curves inFig. 3.3b), eventually become zero. We will use the values of these 16 exponents later onwhen we will discuss the results of Fig. 3.5. -7-6-5-4-3-2-1 0 0 1 2 3 4 5 6 7 (a) l og χ i ( t ) log t Slope = -1 -0.2 0 0.2 0 1 2 3 4 5 6 [g][y][r] (b) χ i ( t ) log t Figure 3.3: Time evolution of the 16 LEs, χ i , i = 1 , . . . ,
16 of the β − FPUT system (3.2) with N = 8 for (a) the regular orbit R and (b) the chaotic orbit C . The plotted straight linein (a) corresponds to a function proportional to t − . In (b) labels indicate green ‘[g]’, yellow‘[y]’ and red ‘[r]’ curves. Axes in (a) are in logarithmic scale, while in (b) the x - and y -axisis in logarithmic and linear scale, respectively.Let us now study the behavior of the GALIs for these two orbits. The regular orbit R lies on an 8 D torus and therefore GALI up to GALI eventually oscillate around a non-zeroconstant value which decreases with increasing order k (Fig. 3.4a), while the GALIs withorder greater than the dimension of the torus go to zero following the asymptotic power39HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSlaw decays in accordance with Eq. (2.33) (Fig. 3.4b). On the other hand, the GALIs of thechaotic orbit C follow exponential decays in accordance with Eq. (2.35) (Fig. 3.5). -12-9-6-3 0 1 2 3 4 5 6 (a) GALI GALI GALI GALI GALI GALI GALI l og [ GA L I k ( t ) / GA L I k ( )] log t -40-30-20-10 0
3 4 5 6 (b)
GALI GALI GALI GALI GALI log t slope=-2 slope=-4 slope=-8 slope=-12slope=-16 Figure 3.4: Time evolution of the GALI k for (a) k = 2 , . . . , k = 9 , . . . ,
16 of theregular orbit R of the β − FPUT system (3.2). The plotted straight lines in (b) correspondto functions proportional to t − , t − , t − , t − , and t − . -24-16-8 0
0 300 600 900 (a)
GALI GALI GALI GALI GALI l og [ GA L I k ( t ) / GA L I k ( )] log t Slope=- (σ −σ ) /ln(10)Slope=- (2σ −σ −σ ) /ln(10)Slope=- (3σ −σ −σ −σ ) /ln(10)Slope=- (4σ −σ −σ −σ −σ ) /ln(10)Slope=- (5σ −σ −σ −σ −σ −σ ) /ln(10)
0 50 100 150 (b)
GALI GALI GALI GALI GALI t Slope=- (7σ −σ −σ −σ −σ −σ −σ ) /ln(10)Slope=- (9σ −σ −σ −σ −σ −σ ) /ln(10)Slope=- (11σ −σ −σ −σ ) /ln(10)Slope=- (13σ −σ ) /ln(10)Slope=- (16σ ) /ln(10) Figure 3.5: Time evolution of the GALI k for (a) k = 2 , . . . , k = 8 , , , , C of the β − FPUT system (3.2). The plotted straight lines correspondto functions that follow the asymptotic decay laws in Eq. (2.35) for σ = 0 . σ = 0 . σ = 0 . σ = 0 . σ = 0 . σ = 0 .
042 and σ = 0 . For Hamiltonian ows and symplectic maps, a random small perturbation of a stable periodicorbit in general leads to regular motion. Hence, finding a stable periodic orbit of the β − FPUTsystem (3.2) and performing small perturbations to it will allow us to study the behavior ofGALI for several regular orbits.The eigenvalues and eigenvectors of the so-called monodromy matrix M ( T ) dictate thestability of the periodic orbit. M ( T ) is the fundamental solution matrix of the variationalequations of the periodic orbit evaluated at a time equal to one period T of the orbit (seee.g. [58]). This matrix is symplectic, and its columns are linearly independent solutions ofthe equations that govern the evolution of deviation vectors from the periodic orbit, i.e. thevariational equations. Since the Hamiltonian system (3.2) is conservative, the mappingdefined by the monodromy matrix is volume-preserving, i.e. the determinant of M ( T ) isone. In addition, the product of the matrix eigenvalues, which is equal to the determinant,is also one.If all the eigenvalues of M ( T ) are on the unit circle in the complex plane, then thecorresponding periodic orbit is stable otherwise it is unstable. The instability can be ofdifferent types (more discussion about this can be found in [58] and references therein) buthere we are only interested in the stability of periodic orbits. So we do not pay much attentionto the types of instability periodic orbits have when they become unstable.In an N D autonomous Hamiltonian system, it turns out that two eigenvalues are alwaysequal to λ = 1 [58]. In practice, the remaining 2( N −
1) eigenvalues define the stability ofthe periodic orbit. In addition, we expect the eigenvalues of M ( T ) computed at any pointof the stable orbit to remain the same since the stability of the orbit does not change alongthe orbit. Thus, we can reduce our investigation to a 2( N − N −
1) eigenvalues, none of which is by default λ = 1.In our analysis, we investigate the behavior of the GALIs for regular orbits in theneighborhood of two simple periodic orbits (SPOs) of the β − FPUT system (3.2), which werefer to as SPO1 and SPO2. In the remaining part of our study, we set β = 1. The existenceand dynamics of these SPOs for the β − FPUT system were discussed in [59, 60].
The first SPO we study is called SPO1 in [60] and it is obtained by considering the β − FPUTlattice (3.2) with N being an odd integer so that all particles at even-numbered positionsare kept stationary at all times, while the odd-numbered particles are always displacedsymmetrically to each other. For example, when N = 7 the 2nd, 4th and 6th particlesare fixed, i.e. x = x = x = 0 and the odd numbered particles, i.e. 1st, 3rd, 5th and 7thare equally displaced but in opposite directions, i.e. x = − x = x = − x (see Fig. 3.6). Figure 3.6: Example of the SPO1 with N=741HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSIn general, the SPO1 for N dof is given by x j ( t ) = 0 , x j − ( t ) = − x j +1 ( t ) = x ( t ) , (3.8)for j = 1 , , . . . , N − . Using condition (3.8) in the system’s equations of motion we endup with a second order nonlinear differential equation for variable x ( t ) which describes theoscillations of all moving particles of the SPO1. In particular, for j = 1 , , , . . . , N we have¨ x j ( t ) = − x ( t ) − βx ( t ) , (3.9)while x j ( t ) = 0, j = 2 , , , . . . , ( N −
1) for all times.In Fig. 3.7 we see the arrangement of the eigenvalues of the monodromy matrix of theSPO1 orbit of Hamiltonian (3.2) with N = 11 when H N /N = 0 . H N /N = 0 . Im( λ ) - 0.5 - -0.5 (a) Re( λ ) Im( λ ) - 0.5 - -0.5 (b) Re( λ ) Figure 3.7: Representation of the eigenvalues λ i , i = 1 , , . . . ,
20, of the monodromy matrixof the SPO1 (3.8) of the β − FPUT system (3.2) on the complex plane with N = 11 for(a) the stable SPO1 with H N /N = 0 . x (0) = − . H N /N = 0 . x (0) = − . H N (3.2)the orbit is stable, but it becomes unstable when the energy increases beyond a certainthreshold H cN . For instance in Fig. 3.7 the transition from stability to instability occurs at H cN /N ≈ . N = 11 particles. In [59, 60] the authors computed the energy densitythreshold H cN /N for different number of particles N , and they found that H cN /N decreasesas the value of N increases following an asymptotic law H cN /N ∝ N − . Fig. 3.8 illustratesthe destabilization energy threshold for SPO1 with respect to increasing N .42HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS -3 -2 -1 H N / N N Slope = -1
Figure 3.8: The energy per particle threshold of the first destabilization H cN /N of the SPO1 ofthe β − FPUT system (3.2) for β = 1. The straight line corresponds to a function proportionalto N − .The main objective of our analysis is to study the behavior of GALIs for regular orbitsthat are located in the neighborhood of the stable periodic orbits. In order to find theseregular orbits we initially start with the stable SPO1 having ICs x i (0), p i (0), i = 1 , , . . . , N and energy density H N /N . Then we perturb this orbit to obtain a nearby orbit with ICs X i (0), P i (0) while keeping the total energy fixed. Practically, in order to find the perturbedorbit, we add a small random number to the positions of the stationary particles of thestable SPOs, i.e. we have x j (0) + w j , j = 2 , , , . . . , ( N −
1) where w j s are small real randomnumbers. Then we change the position of one of the initially moving particles to keep thesame energy value H N .The phase space distance D between these two orbits, x i (0), p i (0) and X i (0), P i (0) isgiven by D = ( N X i =1 (cid:2) ( x i (0) − X i (0)) + ( p i (0) − P i (0)) (cid:3)) / . (3.10)As we have already shown in Sec. 2.2, in multidimensional Hamiltonian systems, thevalue of the GALI k with 2 ≤ k ≤ N remains practically constant for regular orbits lying onan N D torus and these constant values decrease as the order of the index k grows [20]. Inorder to improve our statistical analysis we evaluate the average value of GALI k ( t ) / GALI k (0)over a set of n v = 10 different random initial deviation vectors and denote this quantity as h GALI k ( t ) / GALI k (0) i . We stop our integration when the values of the GALIs more or lesssaturate showing some small fluctuations. Then we estimate the asymptotic GALI k value,denoted by GALI k , by finding the mean value of h GALI k ( t ) / GALI k (0) i over the last n t = 20recorded values and the related error is estimated through the standard deviation of theaveraging process.Fig. 3.9 shows the time evolution of h GALI k ( t ) / GALI k (0) i for a regular orbit withdistance D = 0 .
12 from the stable SPO1 and H N /N = 0 .
01. The gray area around thesecurves represents one standard deviation. We observe that GALIs of higher orders convergeto lower values and need more time to settle to these values. In Fig. 3.9b we see more clearlyhow these final asymptotic values, GALI k , decrease with increasing k .43HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS -20-15-10-5 0 2 3 4 5 6 7 8 (a) GALI GALI GALI GALI GALI GALI l og < GA L I k ( t ) / GA L I k ( ) > log t
2 4 6 8 10 12 | | | (b) GA L I k k Figure 3.9: (a) The time evolution of the averaged GALI, h GALI k ( t ) / GALI k (0) i for k =2, 3,5, 7, 9 and 11, of a regular orbit close to the stable SPO1 with H N /N = 0 .
01 and D = 0 . n v = 10 sets of random initial deviation vectors. (b) Estimationof the asymptotic GALI k value, GALI k with respect to the order k for the results of (a). Theshaded gray area around each curve in (a) and the error bars in (b) denote one standarddeviation.Similar results are obtained for a regular orbit in the neighborhood of SPO1 for N = 21(Fig. 3.10). Again we observe that the almost constant final values of the GALIs drasticallydecrease as the order of the index k increases from GALI to GALI . Also, GALI startssaturating at around t = 10 but GALI becomes constant around t = 10. -60-45-30-15 0 2 3 4 5 6 7 8 (a) GALI GALI GALI GALI GALI GALI GALI l og < GA L I k ( t ) / GA L I k ( ) > log t
2 4 6 8 10 12 14 16 18 20 22 | | | | | | | | (b) GA L I k k Figure 3.10: Similar to Fig. 3.9 but for a regular orbit in the vicinity of the stable SPO1 with H N /N = 0 .
04 and D = 0 . N . We do that by perturbing the ICs of the stable periodic orbit by a small number,i.e. by increasing D (3.10). Fig. 3.11 shows the outcome of this process. There we see howthe behavior of GALIs change with increasing distance D from the stable SPO1, when weonly consider one set of initial deviation vectors. In particular, we see the evolution of threedifferent GALIs, namely GALI (Fig. 3.11a), GALI (Fig. 3.11b), and GALI (Fig. 3.11c) forseveral orbits in the neighborhood of the stable SPO1 with H N /N = 0 .
01 and D = 0 . D = 0 . D = 0 . D = 0 . D = 0 .
22 and D = 0 .
4. All values of GALIs for the smallerdistance D = D go to zero asymptotically following the power law of Eq. (2.33), i.e. GALI with t − [Fig. 3.11a], GALI with t − [Fig. 3.11b], and GALI with t − [Fig. 3.11c]. Thebehavior of the GALIs for the perturbed orbit with the smallest distance from the SPO1, D = D , is still similar to what we expect for a stable periodic orbit. When we increase D the asymptotic behavior of the GALIs start deviating from the above mentioned powerdecay law, and eventually GALI k becomes constant for D = D i , i = 2 , . . . ,
5. Finally, forvery large D values the perturbed orbits will be outside the stability island. This happensfor D = D which leads to chaotic motions. In this case, GALIs go to zero exponentiallyfast. Note that the values of the GALIs for the regular orbit grow with increasing distance D from the stable SPO1. In Fig. 3.12a we see the values of GALI k with respect to increasing D when E = H N /N is kept constant to H N /N = 0 .
01, for regular orbits around the SPO1orbit. We chose an appropriate value for D in order to avoid a transition to chaotic motion. -10-8-6-4-2 0 1 2 3 4 5 6 (a) D D D D D D l og [ GA L I ( t ) / GA L I ( )] log (t) Slope=-1
1 2 3 4 5 6 (b) D D D D D D l og [ GA L I ( t ) / GA L I ( )] log (t) Slope=-3 -40-30-20-10 0 1 2 3 4 5 6 (c) D D D D D D l og [ GA L I ( t ) / GA L I ( )] log (t) Slope=-10
Figure 3.11: Time evolution of (a) GALI ( t ) / GALI (0), (b) GALI ( t ) / GALI (0), and (c)GALI ( t ) / GALI (0) for orbits with H N /N = 0 .
01 and distances D = 0 . D = 0 . D = 0 . D = 0 . D = 0 .
22 and D = 0 . β − FPUT system(3.2) with N = 11. The straight lines correspond to functions proportional to (a) t − , (b) t − and (c) t − .There is another way to move from stability to instability, that is by increasing theenergy densities E = H N /N . We have seen in Fig. 3.8 that as long as we are below thedestabilization energy of SPO1, the periodic orbit is stable. For instance, in the case of N = 11 the SPO1 is stable below H cN /N ≈ . k as a function of the GALI’s order k for regular orbits in the neighborhood of thestable SPO1 for increasing energy densities E = H N /N . In particular, we consider regular45HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSorbits with D = 0 .
12 for E = 0 . D = 0 . E = 0 .
04 and with D = 0 .
01 for E = 0 . E = 0 .
16. Therefore, the GALI k values decrease as we increase the energy density forregular orbits around the stable SPO. -40-30-20-10 0
2 4 6 8 10 | | | | | | | (a) GA L I k k D D D D
2 4 6 8 10 (b) k E E E E Figure 3.12: Estimation of the asymptotic GALI values, GALI k , as a function of their order k for regular orbits in the vicinity of the stable SPO1 of the β − FPUT system (3.2) with N = 11. (a) Regular orbits for E = H N /N = 0 .
01 and distances D = 0 . D = 0 . D = 0 .
08 and D = 0 .
12. (b) Regular orbits with E = 0 .
01 and D = 0 . E = 0 . D = 0 . E = 0 .
12 and D = 0 .
01, and E = 0 .
16 and D = 0 .
01. The error bars inboth panels denote one standard deviation. Note that in this case the destabilization energydensity of the SPO1 is H cN /N ≈ . k for regular orbits in two ways. Weshowed that the asymptotic GALI k values change when we increase the distance D of thestudied orbit from the stable SPO1 for constant energy H N (Fig. 3.12a). We also saw that,the asymptotic value of GALI k decreases with increasing energy densities inside the stabilityisland [Fig. 3.12b].For example, for the stable SPO1 with N = 11 and E = 0 .
01 regular motion occursup to D = 0 . D value is chaotic.Fig. 3.13 shows an approximation of the size of the stability island around the SPO1 orbitwith different energy densities H N /N by finding the largest distance, denoted by D m , inwhich regular motion is observed for N = 11.46HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS D m H N /N Figure 3.13: Estimation of the extent in phase space of the stability island around SPO1: Themaximum value of D , denoted by D m , for which regular motion occurs, as a function of theenergy density H N /N of the β − FPUT Hamiltonian (3.2) with N = 11. The destabilizationof the SPO1 takes place at H cN /N ≈ . N = 21. In Fig. 3.14a we seethe asymptotic, GALI k values for regular orbits in the vicinity of the stable SPO1 of the β − FPUT system (3.2) for E = H N /N = 0 .
002 with respect to increasing distances D = 0 . D = 0 . D = 0 .
08 and D = 0 .
12, while in Fig. 3.12b we see how GALI k changes for E = 0 .
002 and D = 0 . E = 0 .
006 and D = 0 . E = 0 .
01 and D = 0 .
09, and E = 0 . D = 0 . -60-40-20 0 4 8 12 16 20 | | | | | | | | | (b) GA L I k k D D D D
4 8 12 16 20 (a) k E E E E Figure 3.14: Similar to Fig. 3.12 but for regular orbits in the vicinity of the stable SPO1 ofsystem (3.2) for N = 21 (a) with H N /N = 0 . D = 0 . D = 0 . D = 0 . D = 0 .
12, and (b) E = 0 .
002 and D = 0 . E = 0 .
006 and D = 0 . E = 0 . D = 0 .
09, and E = 0 .
04 and D = 0 .
03. The error bars in both panels denote onestandard deviation. Note that in this case the destabilization energy density of the SPO1 is H cN /N ≈ . So far our analysis was based on results obtained in the neighborhood of one SPO type,i.e. SPO1. We will now show that these findings are quite general as they remain similarfor another SPO. In particular, to illustrate that a similar analysis is performed for regularorbits in the neighborhood of what was called SPO2 in [60].The SPO2 of the β − FPUT system (3.2) with N = 5 + 3 m , m = 0 , , , . . . , particles,corresponds to an arrangement where every third particle remains always stationary and thetwo particles in between move in opposite directions, i.e. x j ( t ) = 0 , j = 1 , , , . . . , N − ,x j ( t ) = − x j +1 ( t ) = x ( t ) , j = 1 , , , . . . , N − . (3.11)For example, for N = 8, particles located in the 3rd and 6th position do not move, x = x = 0, while the remaining consecutive particles are moving symmetrically and in oppositedirections, i.e. x = − x = x = − x = x = − x (Fig. 3.15). Figure 3.15: Example of the SPO2 with N=8As in the case of the SPO1, we can obtain a single differential equation for the timeevolution of the SPO2 when substituting (3.11) in the equations of motion of the β − FPUTsystem (3.2). In particular we get¨ x j ( t ) = − x ( t ) − βx ( t ) , (3.12)for the moving particles, while x j ( t ) = 0 for the stationary particles with j = 3 , , , . . . , N − λ i , i = 1 , , . . . ,
20, ofthe monodromy matrix of the SPO2 of Hamiltonian (3.2) with N = 11. In Fig. 3.16a alleigenvalues are on the unit circle, which means that the SPO2 with H N /N = 0 .
01 and x (0) = − . x = − . p >
0. In Fig. 3.16b, there are four eigenvalues outside the unit circle, thusthe SPO2, with H N /N = 0 .
02 and x (0) = − . x = − . p >
0. The critical energy density for whichthe SPO2 encounters its first transition from stability to instability is much smaller that theone seen for SPO1. In particular, it is H cN /N ≈ . Im( λ ) - 0.5 - -0.5 (a) Re( λ ) Im( λ ) - 0.5 - -0.5 (b) Re( λ ) Figure 3.16: Similar to Fig. 3.7 but for 20 eigenvalues of the monodromy matrix of theSPO2 (3.11) of the β − FPUT Hamiltonian (3.2) with N = 11 for (a) the stable SPO2 with H N /N = 0 .
01 and x (0) = − . H N /N = 0 .
02 and x (0) = − . H cN /N and the dimension N of the system(Fig. 3.17). The first destabilization energy density of SPO2 follows a slower asymptoticdecrease compared to SPO1, i.e. H cN /N ∝ N − [59]. -5 -4 -3 -2 -1 5 50 100 300 H N / N N Slope = -2
Figure 3.17: Similar to Fig. 3.8 but for the SPO2. The straight dashed line corresponds to afunction proportional to N − .The estimation of the size of the stability island around the stable SPO2 with N = 11shows a similar behavior to what we observed for the same N in the case of the SPO1(Fig. 3.13). In Fig. 3.18 we report the maximum value of D (i.e. D m ) for which regularmotion occurs as a function of the energy density for the stable SPO2 of (3.2). Note thatthe energy axis is smaller than the one of Fig. 3.13 since the critical energy density for whichthe SPO2 becomes unstable for the first time is H cN /N ≈ . D m H N /N Figure 3.18: Similar to Fig. 3.13 but for the stable SPO2 of the β − FPUT Hamiltonian (3.2)with N = 11.Similarly to Fig. 3.12 in Fig. 3.19a we see the dependence of GALI k on the order k for increasing D and constant H N /N , whereas in Fig. 3.19b we have similar results but forincreasing H N /N values for regular orbits close to the stable SPO2 with N = 11. -30-20-10 0
2 4 6 8 10 | | | | | | | (a) GA L I k k D D D D
2 4 6 8 10 (b) k E E E E Figure 3.19: Similar to Fig. 3.12 but for regular orbits in the vicinity of the stable SPO2 ofthe β − FPUT system (3.2) with N = 11. Regular orbits with ICs (a) E = H N /N = 0 .
001 and D = 0 . D = 0 . D = 0 .
08 and D = 0 .
12, (b) E = 0 .
001 and D = 0 . E = 0 . D = 0 . E = 0 .
008 and D = 0 .
04, and E = 0 .
012 and D = 0 . k values for regular orbits around the stable SPO2 for N = 20. Note that the destabilization threshold H cN /N ≈ . N = 20is much smaller than in the case with N = 11.50HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS -60-40-20 0 4 8 12 16 20 | | | | | | | | | | (a) GA L I k k D D D D
4 8 12 16 20 (b) k E E E E Figure 3.20: Similar to Fig. 3.19 but for N = 20. Regular orbits for (a) E = H N /N = 0 . D = 0 . D = 0 . D = 0 .
08 and D = 0 .
12, (b) E = 0 . D = 0 . E = 0 . D = 0 . E = 0 .
002 and D = 0 .
04, and E = 0 .
004 and D = 0 . k values increase when the regularorbit approaches the edge of the stability island (Figs. 3.19a and 3.20a) while it decreaseswhen this orbit moves towards the destabilization energy ( Figs. 3.19b and 3.20b). Let us now turn our attention to the properties of the deviation vectors needed for thecomputation of the GALIs. Regular motions take place on a torus and the time evolution ofall the initially linearly independent deviation vectors brings them on the tangent space ofthis torus, having in general different directions. On the other hand, in the case of chaoticorbits, the deviation vectors gradually align to the direction defined by the mLE followingthe exponential law given in Eq. (2.35).We start our study by noting that the values of the GALIs are practically independentof the choice of the initial deviation vectors. To illustrate this property, we considerdifferent sets of initial deviation vectors whose coordinates are drawn from three differenttypes of probability distributions and compute the corresponding GALIs. In particular, weconsider a uniform distribution in the interval [ − . , .
5] (green curves in Fig. 3.21), a normaldistribution (blue curves in Fig. 3.21) with mean 0 and standard deviation 1 in the interval[ − . , . h GALI k ( t ) / GALI k (0) i for a regular orbit in the neighborhoodof the stable SPO1 of Hamiltonian (3.2) for N = 11 and H N /N = 0 .
01. The average is doneover n v = 10 sets of normalized unit (denoted by [u]) initial deviation vectors. Fig. 3.21bshows a similar computation but over n v = 10 sets of orthonormalized (denoted by [o])deviation vectors, i.e. GALI k (0) = 1. In Fig. 3.21c we compare asymptotic values, GALI k ,obtained from the results of Figs. 3.21a and 3.21b. The curves are more or less the same.This indicates that the choice of the initial deviation vectors does not affect the values of the51HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSGALIs. The three curves corresponding to the different initial distribution of the vectors’coordinates practically overlap both in Fig. 3.21a and 3.21b which is a clear indication thatthe time evolution of the GALI k is similar, irrespective of the used distribution. -20-15-10-5 0
2 4 6 8 (a)
GALI GALI GALI GALI GALI GALI l og < GA L I k ( t ) / GA L I k ( ) > log t
2 4 6 8 (b)
GALI GALI GALI GALI GALI GALI l og < GA L I k ( t ) > log t
2 4 6 8 10 | | | | | (c) GA L I k k (a)u(a)n(a)e(b)u(b)n(b)e Figure 3.21: The average over different initial sets of deviation vectors time evolution of (a) h GALI k ( t ) / GALI k (0) i and (b) h GALI k ( t ) i for k = 2, 3, 5, 7, 9, 11, for a regular orbit inthe neighborhood of the stable SPO1 ( D = 0 . H N /N = 0 .
01) of Hamiltonian (3.2) with N = 11. The coordinates of the initially (a) unit [(b) orthonormalized] deviation vectors werechosen from a uniform [u] (green curves), normal [n] (blue curves) and exponential [e] (redcurves) distributions. In the legend of panel (c), ‘(a)’ and ‘(b)’ denotes the first and secondpanel, and ‘u’, ‘n’, and ‘e’ denote uniform, normal and exponential distributions respectively.Furthermore, we analyze the distribution of the coordinates of the vectors ˆ v j and thecorresponding angles between these vectors θ ij (seen in Eq. (2.42)). In particular, we considerthe two cases N = 11 and N = 21. In Fig. 3.22a we show vector coordinates distributionsobtained from the evolution of GALI k for a regular orbit of Hamiltonian (3.2) with N = 11which is close to the stable SPO1, with H N /N = 0 .
01 and D = 0 .
12. When we follow thetime evolution of GALI k , k = 2 , . . . ,
11 using unit or orthonormal initial deviation vectors,the coordinates of the k deviation vectors and the corresponding angles on the torus have thesame distribution. We note that for GALI k with high order k we obtain better statistics as wehave in our disposal a large number of vectors and corresponding angles. Thus, considering,for example, GALI k with the largest possible order, i.e. GALI , we obtain a fine picture of thedistributions. More specifically, in order to follow the evolution of GALI we use 11 initialdeviation vectors with 22 ×
11 = 242 coordinates and 11! / (9!2!) = 55 angles. In Fig. 3.22awe plot the probability density distributions of the coordinates of the unit deviation vectorsneeded for the evaluation of GALI for a regular orbit close to the SPO1 of the Hamiltonian(3.2) with N = 11. The distributions are generated from the coordinates of 11 sets ofdeviation vectors over 10 snapshots when the GALI has reached its asymptotic value. Asimilar behavior is observed when increasing the number of particles. The time evolution ofGALI for regular orbits of (3.2) with N = 21 close to the stable SPO1 H N /N = 0 . D = 0 .
04 using 21 initial unit [u] and orthonormal [o] deviation vectors led to the creation ofFig. 3.22b where we plot the distributions obtained from the coordinates of 10 sets of vectorsas in Fig. 3.22a. The two curves, in both Figs. 3.22a and 3.22b practically overlap.52HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS ^ (a) P r ob a b ilit y D e n s it y D i s t r i bu ti on v i uo ^ (b) v i uo Figure 3.22: Probability density distributions of the coordinates ˆ v i of (a) 11, and (b) 21initially unit (blue curves - ‘u’) and orthonormalized (black curves - ‘o’) deviation vectorsfor a regular orbit of Hamiltonian (3.2) close to the stable SPO1, with (a) H N /N = 0 . D = 0 .
12 for N = 11, and (b) H N /N = 0 . D = 0 .
04 for N = 21.We further investigated how the distribution of the angles between each deviationvectors θ ij behaves and a behavior similar to that seen in Fig. 3.22 is obtained. Fig. 3.23a[Fig. 3.23b] displays the probability density distributions of the angles corresponding to thedeviations in Fig. 3.22a [Fig. 3.22b]. The curves are almost the same, with distributions ofthe angles corresponding to vectors whose coordinates are drawn from the initial unit [u] andorthonormal [o]. Π /2 Π ^ (a) P r ob a b ilit y D e n s it y D i s t r i bu ti on θ ij uo Π /2 Π ^ (b) θ ij uo Figure 3.23: Probability density distributions of the angles θ ij corresponding to the deviationvectors in Fig. 3.22 for (a) N = 11 and (b) N = 21.Based on the results from Fig. 3.21 and Fig. 3.22 we can argue that the distribution ofthe coordinates of the initial deviation vectors does not affect the asymptotic behavior of theGALI. Therefore, we can, for example, use initial unit deviation vectors whose coordinatesare drawn from a uniform distribution for our analysis.53HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSSo far we have investigated the behavior of the distribution of the coordinates ofdeviation vectors when the vectors have fallen on the tangent space of the torus, but it isinteresting to also investigate the dynamics of the deviation vectors for the whole evolutionof the regular orbit. In order to do so, we follow the evolution of GALI (Fig. 3.24a) for theregular orbit close to the stable SPO1 of (3.2) for N = 11, with H N /N = 0 .
004 and D = 0 . : 0 ≤ log t < . : 1 . ≤ log t < : 3 ≤ log t < . :4 . ≤ log t < : 6 ≤ log t ≤ . ) the distribution of the coordinates of the deviation vectors hasa sharply peaked shape with a high concentration in the middle. -20-15-10-5 0
0 2 4 6 8 I I I I I (a) l og [ GA L I k ( t ) / GA L I k ( )] log t ^ (b) P r ob a b ilit y D e n s it y D i s t r i bu ti on v i I I I I I -1 -0.5 0 0.5 1 ^ (c) v i I I I I I Figure 3.24: (a) The time evolution of the GALI for a regular orbit in the neighborhood ofthe stable SPO1 orbit of Hamiltonian (3.2) for N = 11 with H N /N = 0 .
004 and D = 0 . , I , I ,I and I of (a). (c) Similar to (b) but for a different random set of initial deviation vectors.Finally we evaluated the asymptotic coordinate distributions of the deviation vectorsused to compute the GALIs for regular orbits in the neighborhood of the stable SPO1Fig. 3.25a and Fig. 3.26a and SPO2 Fig. 3.25b and Fig. 3.26b orbits. In our computationswe use 10 sets of 11 initially linearly independent unit deviation vector. In Fig. 3.25 wesee the dependence of the final coordinate distribution of these vectors when the distance D from the SPO is increased (for fixed H/N values) for the stable SPO1 (Fig. 3.25a) and SPO254HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS ^ (a) P r ob a b ilit y D e n s it y D i s t r i bu ti on v i D D D D ^ (b) v i D D D D Figure 3.25: The probability density distributions of the coordinates of 11 unit deviationvectors used for the time evaluation of GALI for regular orbits close to the stable (a) SPO1and (b) SPO2 of the β − FPUT Hamiltonian (3.2) with N = 11. Regular orbits have (a) E = H N /N = 0 .
01 and distances D (3.10) from the stable SPO1 D = 0 . D = 0 . D = 0 .
12 and D = 0 .
22 and (b) E = H N /N = 0 .
001 and distances D from the stable SPO2 D = 0 . D = 0 . D = 0 .
08 and D = 0 . D and E , the distributions havea peaked shape with a concentration in their middle. As the regular orbits get closer to theboundary of the stability island (increase D in Fig. 3.25 and increase E in Fig. 3.26) thedistributions become more sharply peaked with a very high concentration in their centers.55HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMS ^ (a) P r ob a b ilit y D e n s it y D i s t r i bu ti on v i E E E E ^ (b) v i E E E E Figure 3.26: Similar to Fig. 3.25 but for regular orbits with (a) E = 0 .
001 and D = 0 . E = 0 .
01 and D = 0 . E = 0 . D = 0 . E = 0 .
175 and D = 0 . E = 0 .
001 and D = 0 . E = 0 .
004 and D = 0 . E = 0 .
008 and D = 0 .
04, and E = 0 .
012 and D = 0 . So far we have discussed the behavior of the GALIs for the β − FPUT Hamiltonian model(3.2) which has a continuous time t . In particular, we have studied the asymptotic behaviorof GALIs for regular orbits in the neighborhood of the stable SPOs when their distance D from them increases. Here, we will implement the same methodology for regular orbitsof a multidimensional area-preserving mapping in order to investigate the generality of ourfindings.As we already discussed in Sec. 2.3, the standard mapping describes a universal, genericof area-preserving mapping with divided phase space in which the integrable islands ofstability are surrounded by chaotic regions. A 2 N D system of coupled mappings is: x ′ j = x j + y ′ j ,y ′ j = y j + K j π sin (cid:0) πx j (cid:1) − µ π ( sin h π (cid:0) x j +1 − x j (cid:1)i + sin h π (cid:0) x j − − x j (cid:1)i) , (3.13)where ( ′ ) indicates the new values of variables after one mapping iteration and K j isa dimensionless parameter that influences the degree of chaos, while µ j is the couplingparameter between neighboring mappings where j = 1 , , . . . , N . All values have (mod 1),i.e. 0 ≤ x j <
1, 0 ≤ y j < x = x N and x N +1 = x .For this coupled standard mapping model, the behavior of the GALIs for a given ICdepends on the parameters K j , as K j strongly influences the regular or chaotic nature of themapping’s orbits. Here, the objective is to show that what we have observed so far for thebehavior of the GALIs in the case of the β − FPUT model (Sec. 3.2) and in particular theresults of Figs. 3.9 and 3.11, hold also for the discrete system (3.13). In other words, weinvestigate the behavior of the GALIs when the transition of periodic orbits from regular to56HAPTER 3. BEHAVIOR OF THE GALI FOR REGULAR MOTION INMULTIDIMENSIONAL HAMILTONIAN SYSTEMSchaotic motion happens for the mapping (3.13) using the same analysis. In order to do this,we start from the center of the big island ( x j , y j ) = (0 . ,
0) [27], j = 1 , , . . . , N for a smallvalue of K j and coupling µ and we locate a stable periodic orbit. The location of this islandis shown in Fig. 3.27 for the simple case of the 2D mapping (2.46) x x Figure 3.27: Phase space portraits of the standard mapping (2.46) with K j = 2 and µ = 10 − .As we have seen in Fig. 3.11 when the stable orbit is very close to the SPO, the GALIstend to zero following the power law (2.34) (curve D in Fig. 3.11). Then as the orbit’sIC is moved further and further away from the SPO the GALIs asymptotically tend to apositive value (curves D to D in Fig. 3.11) while when the orbit’s IC is outside the stabilityisland the motion becomes chaotic and the GALIs will go to zero exponentially fast as ithappens for the curve D in Fig. 3.11. All these behaviors can also be observed for themapping, Fig. 3.28 clearly illustrate these results. There we see the behavior of the GALIsfor increasing distances from the center of the island of stability. In particular, the evolutionof the GALI for different orbits starting in the neighborhood of the stable periodic orbit of(3.13) with parameters K j = 2 and µ = 10 − , and distances D = 0, D = 0 . D = 0 . D = 0 . D = 0 . D = 0 . for the periodicorbit ( D = 0) goes to zero asymptotically following the power law t − (2.33). Then, as thevalue of D is increased the GALI gradually starts to saturate to a constant value ( D to D ).Finally, the perturbed orbit becomes unstable and the GALI goes to zero exponentially fast( D = 0 . -6-4-2 0
0 1 2 3 4 5 6 D D D D D D l og [ GA L I ( t ) / GA L I ( )] log n Slope=-1
Figure 3.28: The time evolution of GALI ( t ) / GALI (0) for orbits of the mapping (3.13) with K j = 2 and µ = 10 − , and distances D = 0, D = 0 . D = 0 . D = 0 . D = 0 . x j , y j ) = (0 . , j = 1 , , . . . ,
20 of the stability island for 1 set of initialdeviation vectors. The straight line correspond to a function proportional to n − .In Fig. 3.29 we follow the evolution of the GALI k in a similar way to Figs. 3.9 and3.10. From the results of Fig. 3.29a we see that the GALI k with higher order k needs moretime to become constant. For example, GALI asymptotically saturates around t = 10 . while GALI saturates faster around t = 10 time units. From Fig. 3.29b we see that the finalasymptotic value GALI k drastically decrease as the order of the index k increases. -30-24-18-12-6 0 1 2 3 4 5 6 (a) GALI GALI GALI GALI GALI GALI GALI l og < GA L I k ( t ) / GA L I k ( ) > log n
4 8 12 16 20 | | | | | (b) GA L I k k Figure 3.29: (a) The time evolution of the average, over n v = 10 sets of random initialdeviation vectors, h GALI k ( t ) / GALI k (0) i for k =2, 6, 9, 12, 15, 18 and 20, for a regular orbitof the mapping (3.13) with N = 20. The orbit’s ICs are given in [61] (b) Estimation of theasymptotic GALI values of panel (a), GALI k , as a function of the order k . The shaded grayarea around each curve in panel (a) and the error bars in panel (b) denote one standarddeviation. 58 hapter 4Summary and discussion In this work, we performed several numerical investigations of multidimensional dynamicalsystems by implementing several chaos detection techniques. In the first chapters of thethesis, we gave a brief overview of the basic notions of Hamiltonian dynamics along with apresentation of some basic numerical techniques for investigating chaos. The creation of thePoincar´e surface of section is an important tool to understand the behavior of dynamicalsystems by representing trajectories of the full 2 N D phase space by an object in a 2( N − k (GALI k ) represents thevolume of a parallelepiped formed by k initially linearly independent deviation vectors ofunit length.For the computation of these chaos indicators, a symplectic integration scheme wasused to follow the time evolution of the equations of motion and variational equations ofthe Hamiltonian. In our presentation, we applied these indicators on the 2D H´enon-HeilesHamiltonian, 2D and 4D symplectic mappings as well as on simple dissipative systems.The main results of our work was presented in Chapter 3, where we thoroughlyinvestigated the behavior of the GALI chaos indicator for regular motions in multidimensionalHamiltonian systems. In the case of stable periodic orbits, GALI k goes to zero following thepower law described in [28], while it tends to zero exponentially fast following a rate thatdepends on the values of different LEs for unstable orbits. In our study, more attentionwas given to regular motions that lie on N D tori, i.e. GALI k , k = 2 , , . . . , N asymptoticallyattain a non zero constant value. To analyze these asymptotic values of GALIs we consideredseveral regular orbits in the vicinity of two basic simple periodic orbits of the β − FPUT model(3.2), namely, SPO1 and SPO2 orbits [60], for various numbers of the system’s degrees offreedom.The main results of our study are the following: the asymptotic GALI values, GALI k ,59HAPTER 4. SUMMARY AND DISCUSSION k = 2 , , . . . , N for regular orbits on N D tori of multidimensional Hamiltonian models dependon: • The order k of the index. The constant values of the GALIs decrease as k increases(Figs. 3.9, 3.10, and 3.29). • The phase space distance D of the regular orbit from the nearby stable periodic orbit.The asymptotic GALI values increase when we increase the distance D for a fixedenergy H N (Figs. 3.12a and 3.20a). • The orbit’s energy H N . The GALI k values decrease when we approach thedestabilization energy (Figs. 3.12b and 3.20b).In addition, we computed the evolution of GALIs using initially linearly independentdeviation vectors generated from different random distributions (uniform, normal andexponential). We choose the deviation vectors in two ways: firstly, by normalizing each vectorand secondly, by setting the angles between vectors to be right angles (i.e. orthonormalvectors). In all cases, we observed that the volume of the parallelepiped formed by thedeviation vectors does not depend on the distribution from which the vectors were createdand the way we choose our vectors (Fig. 3.21). This means that the values of GALIs donot depend on the type of initial deviation vectors but rather on the initial conditions.Furthermore, we showed that the shape of the coordinate distributions of the deviationvectors and the angles between them when the GALIs have reached their asymptotic valuesis also independent of the way the initial deviation vectors were produced (Fig. 3.22 andFig. 3.23).Finally, in order to corroborate the generality of our results, we evaluated the asymptoticconstant values of GALIs for regular orbits of the 2 N D multidimensional area-preservingmapping (3.13) and showed that similar results to the ones seen for the β − FPUT modelwere also obtained. 60HAPTER 4. SUMMARY AND DISCUSSION
Dissemination of the results of this work
Parts of this work were included in:1. Paper in a peer review journal • “On the behavior of the Generalized Alignment Index (GALI) method for regularmotion in multidimensional Hamiltonian systems”, Moges H.T. , Manos T. andSkokos Ch.: 2019, Nonlin. Phenom. Complex Syst. (in press), preprint version: https://arxiv.org/abs/2001.00803
2. Conference presentations • “Investigation of Chaos by the Generalized Alignment Index (GALI) Method”,Poster presentation in the International Conference on Mathematical Modeling ofComplex Systems, Pescara (Italy), July 3-11, 2019. • “Investigation of Chaos by the Generalized Alignment Index (GALI) Method”,Poster presentation in the 12th Annual University of Cape Town/StellenboshFaculty of Science Postgraduate Symposium, Stellenbosch (South Africa),September 11, 2019. • “Investigating the properties of regular motion in multidimensional nonlinearlattices by the Generalized Alignment Index (GALI) method”, Oral presentationin the 62nd Annual Congress of the South African Mathematical Society, CapeTown (South Africa), December 2-4, 2019.61 ibliography [1] E. Lorenz, The essence of chaos (Adam Hilger, Bristol, 1993).[2] V. S. Anishchenko, V. V. Astakhov, and T. E. Vadivasova,
Nonlinear Dynamics ofChaotic and Stochastic Systems. (Springer, s.l., Berlin, Heidelberg, 2007).[3] G. Contopoulos,
The development of nonlinear dynamics in astronomy.
Found. Phys. , 89 (2001).[4] R. S. MacKay, and J. D. Meiss, Hamiltonian Dynamical Systems (University ofWashington Press, Seattle, 1987).[5] S. Wiggins,
Global Bifurcations and Chaos: Analytical Methods (Springer-Verlag, NewYork, 1988).[6] A. J. Lichtenberg, and M. A. Lieberman,
Regular and Chaotic Dynamics (2nd edition)Springer-Verlag, Berlin (1992).[7] C. Sim´o, C. ed,
Hamiltonian Systems with Three or More Degrees of Freedom , (NatoScience Series C, Springer Netherlands, 1999).[8] T. Bountis, and Ch. Skokos,
Complex Hamiltonian Dynamics (Springer-Verlag, Berlin2012).[9] A. M. Lyapunov,
The general problem of the stability of motion (Taylor and Francis,London 1992). English translation from the French: Liapounoff, A.: Probl`emeg´en´eralde la stabilit´e du mouvement. Annal. Fac. Sci. Toulouse , 203–474(1907). The Frenchtext was reprinted in Annals Math. Studies Vol. 17 Princeton Univ. Press (1947). Theoriginal was published in Russian by the Mathematical Society of Kharkov in 1892.[10] V. I. Oseledec, A multiplicative ergodic theorem. Characteristic Ljapunov, exponents ofdynamical systems.
Trans. Moscow Math. Soc. , 197 (1968).[11] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Lyapunov CharacteristicExponents for smooth dynamical systems and for Hamiltonian systems; a method forcomputing all of them. Part 1: Theory.
Meccanica , 9 (1980).[12] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Lyapunov CharacteristicExponents for smooth dynamical systems and for Hamiltonian systems; A method forcomputing all of them. Part 2: Numerical application.
Meccanica , 21 (1980).[13] Ch. Skokos, The Lyapunov characteristic exponents and their computation. In Dynamicsof Small Solar System Bodies and Exoplanets, Lect. Notes Phys. , 63 Springer,(2010). 62IBLIOGRAPHY[14] N. P. Maffione, L. A. Darriba, P. M. Cincotta, and C. M. Giordano, A comparison ofdifferent indicators of chaos based on the deviation vectors: application to symplecticmappings.
Cel. Mech. Dyn. Astron. , 285 (2011).[15] L. A. Darriba, N. P. Maffione, P. M. Cincotta, C. M. Giordano,
Comparative study ofvariational chaos indicators and ODE’numerical integrators.
Int. J. Bifurcation Chaos , 1230033 (2012).[16] D. D. Carpintero, N. Maffione, and L. Darriba, LP-VIcode: A program to compute asuite of variational chaos indicators.
Astronomy and Computing , 19 (2014).[17] Ch. Skokos, G. Gottwald and J. Laskar (eds.), Chaos detection and predictability (Lecture Notes in Physics, Vol. 915, Springer Berlin Heidelberg, 2016).[18] Ch. Skokos, T. C. Bountis, and Ch. Antonopoulos, Geometrical properties of localdynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method.Physical. D , 30 (2007).[19] Ch. Skokos, Alignment Indices: A new simple method for determining the ordered orchaotic nature of orbits. J. Phys. A: Mathematical and General, , 10029 (2001).[20] Ch. Skokos, T. Bountis, and Ch. Antonopoulos, Detecting chaos, determining thedimensions of tori and predicting slow diffusion in fermipastaulam lattices by thegeneralized alignment index method. Eur. Phys. J. Spec. Top. , 5 (2008).[21] Ch. Skokos, and T. Manos, The Smaller (SALI) and the Generalized (GALI) alignmentindices: Efficient methods of chaos detection.
Lect. Notes Phys. , 129 (2016).[22] E. Fermi, P. Pasta, S. Ulam, and M. Tsingou,
Studies of the nonlinear problems. I , LosAlamos Report
LA-1940 (1955).[23] J. Ford,
The Fermi-Pasta-Ulam problem: Paradox turns discovery.
Phys. Rep. , 271(1992).[24] D. K. Campbell, P. Rosenau, and G. M. Zaslavsky (eds.),
Focus Issue: The’Fermi-Pasta-Ulam’ Problem-The First Fifty Years , Chaos The Fermi-Pasta-Ulam Problem. A Status Report (Lecture Notes inPhysics, Vol. 728, Springer Berlin Heidelberg, 2008).[26] H. Kantz, and P. Grassberger,
Internal Arnold diffusion and chaos thresholds in coupledsymplectic maps. , J. Phys. A , L127 (1988).[27] T. Bountis, T. Manos, and H. Christodoulidi, Application of the GALI method tolocalization dynamics in nonlinear systems.
J. Comp. Appl. Math. , 17 (2009).[28] T. Manos, Ch. Skokos, and Ch. Antonopoulos,
Probing the local dynamics of periodicorbits by the generalized alignment index (GALI) method.
Int. J. Bifurcation Chaos ,1250218 (2012).[29] M. H´enon and C. Heiles, The applicability of the third integral of motion: some numericalexperiments. . The Astronomical Journal , 73 (1964).63IBLIOGRAPHY[30] E. Gregersen, Math Explained: The Britannica Guide to Analysis and Calculus (Britannica Digital Learning, New York 2011).[31] R. Devaney,
Chaotic Dynamical System (Addison-Wesley, New York 1989).[32] H. Korsch, and H. Jodl,
Chaos: A Program Collection for the PC. (Springer-Verlag,New York 1989).[33] M. Banados, I. Reyes, A short review on Noether’s theorems, gauge symmetries andboundary terms. International Journal of Modern Physics D, , 1630021 (2016).[34] H. Yoshida, Construction of higher order symplectic integrators.
Phys. Lett. A , 262(1990).[35] Ch. Skokos, and E. Gerlach,
Numerical integration of variational equations.
Phys. Rev. E , 036704 (2010).[36] S. Blanes, F. Casas, A. Farres, J. Laskar, J. Makazaga, and A. Murua, New familiesof symplectic splitting methods for numerical integration in dynamical astronomy.
App. Num. Math. , 58 (2013).[37] B. Senyange, and Ch. Skokos, Computational efficiency of symplectic integrationschemes: application to multidimensional disordered KleinGordon lattices.
Eur. Phys. J. Spec. Top. , 625 (2018).[38] C. Danieli, B. Many Manda, T. Mithun, and Ch. Skokos,
Computational efficiencyof numerical integration methods for the tangent dynamics of many-body Hamiltoniansystems in one and two spatial dimensions.
Math. Engin. , 447 (2019).[39] E. Forest and R.D. Ruth , Fourth-Order Symplectic Integration.
Physica D , 105(1989).[40] M. Suzuki , Fractal decomposition of exponential operators with applications tomany-body theories and Monte Carlo simulations.
Phys. Lett. A , 319 (1990).[41] E. Forest ,
Geometric integration for particle accelerators.
J. Phys. A , 5351 (2006).[42] J. Laskar and P. Robutel , High order symplectic integrators for perturbed Hamiltoniansystems.
Cel. Mech. Dyn. Astr. , 39 (2000).[43] Wikipedia, Poincar´e map, Wikipedia, the Free Encyclopedia, https://en.wikipedia.org/wiki/Poincar%C3%A9_map .Accessed Sep 2019.[44] M. H´enon, On the numerical computation of Poincar´e maps . Physica D: NonlinearPhenomena, : 412 (1982).[45] C. Froeschl´e,
Numerical study of a four-dimensional mapping.
Astronomy andAstrophysics, , 172 (1972).[46] B. Chirikov, A universal instability of many-dimensional oscillator systems.
Physicsreports, 52(5): 263 (1979). 64IBLIOGRAPHY[47] Ch. Skokos, Ch. Antonopoulos, T. C. Bountis, and M. N. Vrahatis,
How does the SmallerAlignment Index (SALI) distinguish order from chaos?
Prog. Theor. Phys. Supp. ,439 (2003).[48] Ch. Skokos, Ch. Antonopoulos, T. C. Bountis, and M. N. Vrahatis, Detecting orderand chaos in Hamiltonian systems by the SALI method, J. Phys. A: Mathematical andGeneral, , 6269 (2004).[49] Wikipedia, Singular value decomposition, Wikipedia, the Free Encyclopedia, https://en.wikipedia.org/wiki/Singular_value_decomposition . Accessed Feb2019.[50] G. Contopoulos, and A. Giorgilli Bifurcations and complex instability in a 4-dimensionalsymplectic mapping.
Meccanica, : 19 (1988).[51] Ch. Skokos, G. Contopoulos, and C. Polymilis Structures in the phase space of fourdimensional symplectic map.
Celest. Mech. Dyn. Astron, : 223 (1997).[52] H. Wen, A review of the H´enon map and its physical interpretations.
School of PhysicsGeorgia Institute of Technology, Atlanta, GA, U.S.A., : 30332 (2014).[53] M. H´enon,
A two-dimensional mapping with a strange attractor.
Gordon and Breach,Science Publishers, S.A., : 181 (1986).[54] Z. Elhadj, and J.C.Sprott, A minimal 2-D quadratic map with quasi-periodic route tochaos.
Inter. J. Bifur and Chaos, : 1023 (2009).[55] T. Bountis, Analytical and numerical studies of 4-D mapping models of colliding electronbeams.
Communications in Mathematical Physics, : 67 (1976).[56] E. Gerlach, and Ch. Skokos, Comparing the efficiency of numerical techniques for theintegration of variational equations .
Discr. Cont. Dyn. Sys.-Supp. 2011, 475 (2011).[57] E. Gerlach, S. Eggl, and Ch. Skokos,
Efficient integration of the variational equations ofmulti-dimensional Hamiltonian systems: Application to the Fermi-Pasta-Ulam lattice.
Int. J. Bifurcation Chaos , 1250216 (2012).[58] Ch. Skokos, On the stability of periodic orbits of high dimensional autonomousHamiltonian systems.
Physica D , 155 (2001).[59] Ch. Antonopoulos, T. Bountis, and Ch. Skokos,
Chaotic Dynamics of Ndegree of FreedomHamiltonian Systems.
Int. J. Bifurcation Chaos , 1777 (2006).[60] Ch. Antonopoulos, and T. Bountis, Stability of simple periodic orbits and chaos in aFermi-Pasta-Ulam lattice.
Phys. Rev. E, , 056206 (2006).[61] IC is chosen using the procedure in [27] with µ = 10 − and K j ’s in triplets of − . , − . , − .
55, for j = 1 , , . . . ,
20 and theresults are x = ( − . , . , . , − . , − . , . , . , . , . , − . , . , − . , . , − . , − . , − . , − . , − . , . , . y = (0 . , − . , − . , . , . , − . , . , − . , − . , . , . , . , − . − . , . , . , . , . , − . , − ..