# Lyapunov Exponents and Strange Attractors in Discrete and Continuous Dynamical Systems

Lyapunov Exponents and Strange Attractors in Discrete and Continuous Dynamical Systems

Jo Bovy Jo.Bovy@student.kuleuven.ac.be Theoretical Physics Project September 11, 2004

Contents

1 Introduction 2 Overview 3 Discrete and continuous dynamical tems 3.1 The Logistic Map . . . . . . . . . . 3.2 The H? enon Map . . . . . . . . . . 3.3 The Lozi Map . . . . . . . . . . . . 3.4 The Zaslavskii Map . . . . . . . . . 3.5 The Lorenz System . . . . . . . . . 3.6 The R¨ ossler System . . . . . . . . sys. . . . . . . . . . . . 2 2 5.2

5.1.2 5.1.3 5.1.4 5.1.5 5.1.6

4 Lyapunov Exponents 4.1 De?nition and basic properties . . . . 4.2 Constraints on the Lyapunov exponents 4.3 Calculating the largest Lyapunov exponent - method 1 . . . . . . . . . . . 4.4 Calculating the largest Lyapunov exponent - method 2 . . . . . . . . . . . 4.4.1 Maps . . . . . . . . . . . . . . 4.4.2 Continuous systems . . . . . . 4.5 Calculating the other Lyapunov exponents . . . . . . . . . . . . . . . . . . . 4.6 Numerical Results . . . . . . . . . . . 5 Strange Attractors 5.1 Dimensions and de?nition of a strange attractor . . . . . . . . . . . . . . . . . 5.1.1 Topological dimension . . . . .

2 2 3 5.3 4 4 6 Conclusion 17 4 5 A Oseledec’s multiplicative ergodic theorem 18 A.1 The theorem . . . . . . . . . . . . . . 18 5 A.2 A measure on the attractors . . . . . . 19 6 7 B Computer Programs 19 7 8 8 8 9 9 12 12 12 1

Box-counting dimension . . . . Correlation dimension . . . . . Kaplan-Yorke dimension . . . . De?nition of a strange attractor Concerning dimensions from experimental data . . . . . . . Algorithms for calculating dimensions 5.2.1 Box-counting dimension . . . . 5.2.2 Correlation dimension . . . . . Results . . . . . . . . . . . . . . . . . .

12 13 13 14 14 14 14 15 15

3

DISCRETE AND CONTINUOUS DYNAMICAL SYSTEMS

2

1

same algoritmes can be used on a large variety of dynamical systems. Therefore I shall use the techHaving already introduced a chaotic system (the niques that I develop later in this text not only on Lorenz system) in a previous paper [23] and hav- the Lorenz system, but on a selection of the most iming already studied chaos qualitatively in that paper, portant and famous non-linear dynamical systems. In I will, in this paper, try to obtain some quantita- this section I will introduce these dynamical systems tive results concerning chaos. Therefor I will intro- and list there basic properties. I begin by describing duce some more dynamical systems, discrete systems some discrete systems, i.e. two-dimensional maps. I (meaning they have discrete time steps) as well as will also introduce a second three-dimensional concontinuous systems. Discrete systems are a lot eas- tinuous dynamical system in addition to the Lorenz ier to handle than continuous systems. They can be system: the R¨ ossler system. chaotic even in less than two dimensions (which is precluded for continuous systems by the Poincar? eBendixon theorem). And their solution doesn’t in- 3.1 The Logistic Map volve solving di?erential equations. The Logistic map is obtained by replacing the LogisI will then study two aspects of these chaotic sys- tic equation for population growth with a quadratic tems. Firstly I will concentrate on measuring the recurrence relation. This model was ?rst published chaos in the system. This is done by introducing by Pierre Verhulst [27, 28]. Lyapunov exponents, which measure the exponential The Logistic equation is divergence of nearby trajectories. When a Lyapunov dx exponents is positive, we will say that the system is = rx(1 ? x) (1) chaotic. dt All these systems also show a strange attractor for The Logistic map is given by certain parameter values. We will calculate the dimensions of these attractors and see that the dimenxn+1 = rxn (1 ? xn ) (2) sions don’t have to be an integer. This fact is the reason why we call them strange. with r a positive constant sometimes known as the “biotic potential”. We can understand this equation as follows [5]: as2 Overview suming constant conditions every year, the population (of insects for example) at year n uniquely deterFirstly I will describe some dynamical systems, both mines the population at year n + 1. Therefor we have discrete as continuous, and give their basic proper- a one-dimensional map. Say that there are zn insects ties. Then I will introduce the Lyapunov exponents at year n and that every insect lays, on the average, r and give methods for their calculation. Next I shall eggs, each of which hatches at year n + 1. This would move on to strange attractors, give their de?nition yield a population at year n + 1 of zn+1 = rzn . When and explain how I used these de?nitions to write pro- r > 1, this yields an exponentially increasing populagrams that calculate their dimension. tion. When the population is too large however, the insects will exhaust there food supply as they eat and grow, and not all insects will reach maturity. Hence 3 Discrete and continuous dy- the average number of eggs laid will become less than r as zn is increased. A simple assumption is then that namical systems the number of eggs laid per insect decreases linearly The techniques for calculating Lyapunov exponents with the insect population, r[1 ? (zn /z )], where z and dimensions of strange attractors are not speci?c is the population at which the insects exhaust their to the dynamical system under investigation. The food supply. We then have the one-dimensional map

Introduction

3

DISCRETE AND CONTINUOUS DYNAMICAL SYSTEMS

3

zn+1 = rzn [1 ? (zn /z )]. Dividing both sides by z and The H? enon equations are letting x = z/z , we obtain the Logistic map (2). In general, this equation cannot be solved in closed xn+1 = 1 ? ax2 (7) n + yn form and it is capable of very complicated behavior. yn+1 = bxn (8) Many aspects of this equation and it’s chaotic behavior can however be studied exactly. We’ll take This map is invertible, with the inverted map given r ∈ [0, 4] to keep all xn in the interval [0, 1] [10]. The by Jacobian is yn+1 dxn+1 xn = (9) =| r(1 ? 2xn ) | (3) J= b dxn 2 y (10) yn = xn+1 ? 1 + a n+1 b2 The map is stable at a point x0 if J (x0 ) < 1. To ?nd the ?xed points of the map we have to This makes the calculation of basins of attraction solve the equation (dropping the subscripts on x for much simpler. We shall however not look at these convenience) basins of attraction (see for instance [8, 14]). x = rx(1 ? x) (4) The Jacobian matrix J of the H? enon map is so the ?xed points are x1 = 0 and x2 = 1 ? 1 r . The Jacobian evaluated at these points gives is J (x1 ) =| r |

(1) (1) (1)

?2ax 1 b 0

The determinant of this Jacobian is b. By a basic (6) result of linear algebra, the factor by which an area grows under a linear transformation is given by the so the ?rst ?xed point is stable until r = 1, where absolute value of the determinant of the matrix repreit becomes unstable. The second ?xed point is un- senting the transformation. Locally we can linearize enon transformation, so a small area near a stable for r < 1 and stable for 1 < r < 3. At the H? r = 3 both ?xed points are unstable. At this value point P = (x, y ) is reduced by the factor given by for r period doubling occurs. The system begins a the absolute value of the determinant of the derivaperiod-doubling route to chaos that is completed at tive (the Jacobian matrix) of the transformation at r = 3.5699456. At this value the system becomes that point. This gives |det(J )| = |b|, a constant value, chaotic. This is the value we will use to calculate the not dependent on the location of P . Lyapunov exponent and dimension of the strange atThe ?xed points of the H? enon map are tractor that occurs. √ ?1 + b ± 1 ? 2b + b2 + 4a (x, y ) = (1, b) (11) 2a 3.2 The H? enon Map

(1) J (x2 )

(5)

=| 2 ? r |

H? enon introduced this map as a simpli?ed version of the Poincar? e map of the Lorenz system [25]. The Poincar? e map of a system is the map which relates the coordinates of one point at which the trajectory crosses a Poincar? e plane to the coordinate of the next (in time) crossing point. The existence of this map is assumed to be the consequence of the uniqueness of the solution of the equation of the dynamical system.

When we calculate the eigenvalues of the Jacobian matrix J , we can see when these ?xed points are stable. The eigenvalues are λ(1,2) = ?ax ± a2 x2 + b (12)

For a = 1.4 and b = 0.3 (the values for which the H? enon map has a strange attractor), these eigenval-

3

DISCRETE AND CONTINUOUS DYNAMICAL SYSTEMS

4

ues for the ?xed points are λ1 = ?0.092029562 λ2 = 0.6915136742 λ1 = ?0.8079567198 λ2 = 0.1559463222 so both ?xed points are not-stable.

(2) (2) (1) (1)

(again a = 1.7, b = 0.5) are

(1) λ1 (1)

=

λ2 = λ1 = λ2 =

(2) (2)

√ 489 ?17 ? 20 20 √ 489 ?17 + 20 √ 20 489 17 ? 20 √20 489 17 ? 20 20

= ?1.95566722 = 0.25566722

= ?0.25566722 = 1.95566722

3.3

The Lozi Map

so both points are not-stable.

The Lozi map is a simpli?cation of the H? enon map. The quadratic term in xn (?ax2 The Zaslavskii Map n ) is replaced by 3.4 ?a|xn |. The Lozi map is then given by The Zaslavksii map [13] is given by the following xn+1 = 1 ? a|xn | + yn (13) equations yn+1 = bxn (14) xn+1 = xn + ν + ayn+1 (mod1) yn+1 = cos(2πxn ) + e

?r

(19) (20)

We can write the Jacobian matrix using the Heaviside step-function Θ(x) = 0 x<0 1 x≥0 (15)

yn

To compute the Jacobian, we write the equations as follows xn+1 = xn + ν + a cos(2πxn ) + e?r yn (mod1) (21) yn+1 = cos(2πxn ) + e?r yn The Jacobian is then (22)

The Jacobian J then becomes ?a 2Θ(x) ? 1 b 1 0

The determinant is still b. 1 ? 2πa sin(2πxn ) ae?r This map has a strange attractor we shall study ?2π sin(2πxn ) e?r for the parameter values a = 1.7, b = 0.5. The ?xed points for these parameter values are The determinant is equal to e?r , so areas shrink by a factor of e?r every iteration (r > 0). 5 5 The Zaslavskii map shows a strange attractor for (x, y )1 = ( , ) (16) 11 22 ν = 400, r = 3, a = 12.6695. ?5 ?5 (x, y )2 = ( , ) (17) 6 12

3.5

The Lorenz System

The eigenvalues of the Jacobian matrix are I have already studied the Lorenz system [7] extensively in a previous project [23]. I studied the basic λ(1,2) (18) properties of the Lorenz system (?xed points, stability and basins of attraction). So the reader should and the eigenvalues evaluated at the ?xed points consult this for an overview of the basic properties. ?a 2Θ(x) ? 1 ± = 2 √ a2 + 4b

4

LYAPUNOV EXPONENTS

5

Other good references are [15, 31]. I will just state the Lorenz equations ˙ = p(Y ? X ) X ˙ = rX ? XZ ? Y Y ˙ = XY ? bZ Z

Setting the parameter c = 5.7 (we shall study the strange attractor that occurs at this parameter value) the eigenvalues are λ1 = 0.0970008560175134871 + i0.995193491034748634 λ2 = 0.0970008560175134871 ? i0.995193491034748634 λ3 = ?5.68697550703502762 and λ1 = ?0.459615167119897806e ? 5 + i5.42802593149083901 λ2 = ?0.459615167119897806e ? 5 ? i5.42802593149083901 λ3 = 0.192982987303353531

(23)

I shall study the Lorenz strange attractor for the parameter values p = 10, b = 8/3, r = 28.

3.6

The R¨ ossler System

The R¨ ossler system [26] is one of the most simple so both ?xed points are not-stable. chaotic continuous systems. It is arti?cially designed solely with the purpose of creating a simple model for a chaotic strange attractor. The system of di?erential 4 Lyapunov Exponents equations is As with everything else in Physics, and especially in x ˙ = ?(y + z ) (24) Theoretical Physics, we don’t content ourselves with y ˙ = x + ay (25) just a qualitative picture of chaos. In this section I’ll introduce a quantitative measure of chaos, the z ˙ = b + xz ? cz (26) celebrated Lyapunov Exponents. Introducing a quantitative measure of chaos is imwith a,b and c adjustable paramters of which a and portant for several reasons. Most importantly it alb are usely ?xed at a = 0.2, b = 0.2. This system is lows us to de?ne exactly what we mean by chaos. obviously simpler than the Lorenz system, because it When we only have a qualitative picture of chaos, contains only one non-linear term. It is however not everybody has a di?erent opinion of what chaos is. the simplest chaotic ?ow, since the simplest dissipa... One just looks for instance at the picture in phase 2 tive ?ow has been argued to be x +Ax ¨?x ˙ +x = 0 [17] space and decides wether or not he ?nds it chaotic. (This equation can easily be converted to a system of Science would not have come as far as it is now if ?rst order di?erential equations using the following we had contented ourselves with such qualitative picsubstitutions: x ˙ = y, x ¨ = z ). tures. So introducing a measure of chaos allows us to The divergence of the ?ow, ?· f (writing the system rigourously de?ne what we mean by chaos. ? y ˙ x ˙ ?z ˙ ˙ = f) = ? as y ?x + ?y + ?z = a + x ? c. Note that this Having this measure of chaos allows us to go furis not a constant, hence the shrinking/expansion of ther and compare di?erent systems. We can de?ne volumes is not uniform over phase space. what we mean by saying that one system is more The ?xed points are chaotic than another system. Thus we can compare the chaoticness of a system with di?erent parameter √ c ± c2 ? 4ab (x, y, z ) = (a, ?1, 1) (27) values, or the chaoticness of two completely di?erent 2a systems. In the ?rst subsection of this section I will de?ne The Jacobian matrix is and discuss the Lyapunov exponents. In the next ? ? subsections I will discuss how I calculated the Lya0 ?1 ?1 ? 1 a punov exponents and state the results I got for the 0 ? maps described in the previous section. z 0 x?c

4

LYAPUNOV EXPONENTS

6 decomposed in the orthonormal eigenvectors ej of Hn

4.1

De?nition and basic properties

n Since we want to measure chaotic behavior, which we u0 = aj ej (31) intuitively de?ned as sensitivity on initial conditions, j =1 we shall look at the evolution of a small displacement of a initial condition x0 (with corresponding orbit xn and the Lyapunov exponent of such an initial disn = 0, 1, 2, . . .). First take a map M. If we consider placement can be written in function of the eigenvalan in?nitesimal displacement from x0 in the direction ues of Hn (x0 ), since of a tangent vector y0 , the evolution of the tangent n vector, given by n u? · H ( x ) · u = a2 (32) 0 0 j exp[2nλjn (x0 )] 0 yn+1 = DM(xn ) · yn (28) j =1

determines the evolution of the in?nitesimal displacement of the orbit from the unperturbed orbit xn (DM is just the Jacobian matrix). So |yn |/|y0 | is the factor by which the in?nitesimal displacement grows or shrinks. From (28) we have yn = DMn (x0 ) · y0 , where DMn (x0 ) = DM(xn?1 ) · DM(xn?2 ) · . . . · DM(x0 ). We then de?ne the Lyapunov exponent1 for initial condition x0 and initial orientation of the in?nitesimal displacement given by u0 = y0 /|y0 | as λ(x0 , u0 ) = lim 1 ln(|yn |/|y0 |) n 1 = lim ln|DMn (x0 ) · u0 | n→∞ n

n→∞

(29)

An important question now is wether the limits in (29) exist. This existence is guaranteed by Oseledec’s multiplicative ergodic theorem under very general circumstances. This theorem is stated in appendix A with some of its consequences. One of its consequences is that the Lyapunov exponents are the same for almost every x0 with respect to the natural measure on the attractor (this natural measure is also described in appendix A). Now I will move on to de?ne Lyapunov exponents in continuous time systems. All the above considerations remain valid when we replace Eq. (29) by λ(x0 , u0 ) = lim 1 ln(|yn |/|y0 |) n 1 = lim ln|O(x0 , t) · u0 | n→∞ n

n→∞

(33) Note that these exponents are (for now) dependent on the initial condition. If the dimension of the map is N , then there will dy(t) be N Lyapunov exponents, since there are N inde- where dt = F(x(t)) · y(t), x0 = x(0), y0 = y(0), u0 = y(0)/|y(0)|, and O(x0 , t) is the matrix solution pendent initial displacement directions. Since of the equation 1 λ(x0 , u0 ) λn (x0 , u0 ) ≡ ln|DMn (x0 ) · u0 | dO/dt = DF(x(t)) · O (34) n 1 n = ln|u? with initial condition 0 · H (x0 ) · u0 | 2n (30) O(x0 , 0) = I where Hn = [DMn ]? DMn , and ? denotes the This equation is called the variational equation. To transpose, we get N eigenvalues of Hn (Hjn ) each calculate the Lyapunov exponents we’ll have to solve corresponding to one Lyapunov exponent (λjn = an additional system of di?erential equations. (2n)?1 ln Hjn ). Every initial displacement u0 can be Now that we have de?ned Lyapunov exponents, we can de?ne what we mean by saying that a system is 1 In the literature one encounters often reference to Lyapunov numbers ?. These are given in terms of the Lyapunov chaotic. We say that a system is chaotic when it exponents λ by ? = exp(λ). has at least one strictly positive Lyapunov exponent.

4

LYAPUNOV EXPONENTS

7

When a system has a positive Lyapunov exponent, small disturbances will give rise to exponential divergence. So this de?nition is in accordance with the qualitative picture we had. Now how can we represent what these di?erent Lyapunov exponents mean. When we consider a small ball (say in?nitesimal, radius dr) of initial conditions around an initial condition x0 , and then evolve every initial condition inside this ball for n iterates. The ball will have evolved into an ellipsoid . In the limit of large time the Lyapunov exponents give the time rate of exponential growth or shrinking of the principal axes of the evolving ellipsoid. The axes will be (approximately) given by the expressions a1 = exp(nλ1 )dr and a2 = exp(nλ2 )dr.

tonomous, it will follow the same trajectory as the not-disturbed trajectory. These trajectories clearly do not diverge from each other, on the average they will keep at the same distance of each other. So in this direction the Lyapunov exponent is zero. So in every continuous dynamical system one of the Lyapunov exponents will be zero (except for systems with a ?xed point). Applying this to the Lorenz system, which has a constant divergence equal to ?(1 + b + p), we have the following two relations λ1 + λ2 + λ3 = ?(1 + b + p) λ2 = 0 (35) (36)

4.2

Constraints on the Lyapunov exponents

so again we could only calculate the largest Lyapunov exponent (λ1 ), and deduce the other two from these two relations. We can’t use this on the R¨ ossler system however because its divergence isn’t a constant.

We can deduce some constraints on the Lyapunov exponents of a dynamical system. Using these constraints we should only calculate the largest Lyapunov exponent for all but one of the dynamical systems I described in section 3. I shall however in the following sections calculate all Lyapunov exponents and use them to check the validity of these relations. In systems where the area-reducing factor (or volume reducing factor in three dimensional systems) is constant, we can derive a relation between the Lyapunov exponents. Since exp( i λi ) is equal to the area-reduction, this must be equal to the determinant of the Jacobian (for maps) or to exp(? · f ) (? · f is the divergence of the ?ow for continuous dynamical systems). Applying this to the H? enon system, we see that λ1 + λ2 = ln(|b|). We could only calculate the largest exponent and derive the second from this relation. For a continuous dynamical system we can look at disturbances in the direction of the velocity vector of the trajectory (see page 718 in [14]). Consider a trajectory (x(t), y (t), z (t)) and a perturbed trajectory (? x(t), y ?(t), z ?(t)) with initial condition x ?(0) = x(δ ), y ?(0) = y (δ ), z ?(0) = z (δ ) with δ > 0. So the trajectory starts on the reference trajectory, and, if we assume the system to be au-

4.3

Calculating the largest Lyapunov exponent - method 1

As pointed out in appendix A, the largest Lyapunov exponent is the easiest one to calculate. If we start by choosing a direction vector randomly, we can decompose it as in (31). It will most probably have a non-zero component in the direction of the eigenvector of the largest exponent. Since in the cases we consider only one Lyapunov exponent is positive, evolution will be dominated by the largest exponent, as seen from (32). We can then calculate this exponent by just straightforwardly using our qualitative understanding (in the next subsection I will give a better method based on the di?erential map). We have to calculate λ1 = lim lim

n→∞ E0 →0

1 ln(|En |/|E0 |) n

(37)

with E0 an initial disturbance. We can write |En | |En | |En?1 | |E1 | = ··· |E0 | |En?1 | |En?2 | |E0 | substituting (38) in (37), we get λ1 = lim lim 1 n→∞ E0 →0 n

n

(38)

ln(

k=1

|Ek | ) |Ek?1 |

(39)

4

LYAPUNOV EXPONENTS

8

We will however in our algorithm renormalize our error to some chosen value . For one-dimensional maps such as the Logistic map, we essentially haven’t got any choice of direction of the error. For twoor three-dimensional systems we have to choose the direction of the error. For two-dimensional maps we’ll de?ne an error of size on (x, y ) (? x, y ?) = (x + cos φ, y + sin φ) (40)

4.4.1

Maps

For maps this is not such a hard task. We just calculate the Jacobian matrix, and this gives us the differential map. This map tells us how an in?nitesimal small error gets ampli?ed. The algorithm is then as follows 1. We choose an initial point an let the map iterate this for say 100 times. We do this to let transients die out. 2. We consider for an arbitrary angle φ the direction (cos φ, sin φ) (generalizing this to three dimensions is obvious). 3. We compute the next point of the map H (x, y ) and the transformed error E (x, y ) using the differential map. So we get the point E (x, y ) = DM(x, y ) · (cos φ, sin φ)? . 4. The error has increased by d = E (x, y ) . We add log d to an accumulator. 5. We renormalize the direction vector to a unit vector (E (x, y ) = E (x, y )/d). Then we go back to step 3 using the new point and the new error. 6. The result is once again the accumulator divided by the number of iterations. 4.4.2 Continuous systems

with φ an arbitrary angle. For three-dimensional systems this becomes on (x, y, x) (? x, y ?, z ?) = (x + sin φ cos σ, y + sin φ sin σ, z + cos φ) (41) The algorithm2 now works as follows (it is described for maps, but it is easily adjusted to continuous systems): 1. We choose an initial point an let the map iterate this for say 100 times. We do this to let transients die out. 2. We compute the perturbed point according to (40) or (41). 3. We iterate both points and compute the distance d between them 4. We add log d/ to an accumulator 5. We renormalize the error 6. We iterate steps 3-5

The algorithm is essentially the same as the one described in the previous section. We only have to ad7. The result is the average of the log di / , thus the just the manner of computing the transformed direcaccumulator divide by the number of iterations tion vector. To ?nd this transformed direction vector, we have to solve the variational equation. Also in continuous systems, we have to choose the time 4.4 Calculating the largest Lyapunov interval over which each iteration reaches. When we consider discrete maps, we have always iterated the exponent - method 2 map one discrete time step. When considering conThe program outlined in the previous subsection isn’t tinuous maps we choose to let the system evolve for quite accurate, because we should actually take the one timestep (such that t n+1 = tn + 1). So when we limit → 0. We can take this limit by considering choose an integration step of 1/N in our numerical the di?erential map. solver, we have iterate this for N times. After these 2 this algorithm and the ones in the following section are N times we look at our new direction vector, and written along the lines given in [14] pg.710 measure its ampli?cation.

4

LYAPUNOV EXPONENTS

9

4.5

Calculating the other Lyapunov have to be careful about the renormalization however. Because the ?rst Lyapunov is the largest, the exponents

area will shrink very quickly, and the parallellogram spanned will become more and more a line in the direction of the ?rst eigenvector. Thus we shall renormalize the two transformed vectors by Gram-Schmidt orthogonalization, so that they once again form a orthonormal pair. I shall give an overview of the algorithm I used for a two-dimensional map 1. We choose an initial point an let the map iterate this for say 100 times. We do this to let transients die out. 2. We consider for an arbitrary angle φ the direction (cos φ, sin φ) and the direction (? sin φ, cos φ) (generalizing this to three dimensions is obvious). 3. We compute the next point of the map H (x, y ) and the transformed errors E1 (x, y ), E2 (x, y ) using the di?erential map. So we get the points E1 (x, y ) = DM(x, y ) · (cos φ, sin φ)? and E2 (x, y ) = DM(x, y ) · (? sin φ, cos φ)? . 4. The area has increased/shrunk by d = det |E1 (x, y ) E2 (x, y )|. We add log d to an accumulator. 5. We renormalize the directions using GramSchmidt orthogonalization. Then we go back to step 3 using the new point and the new errors. 6. The result is once again the accumulator divided by the number of iterations. The adjustments necessary for continuous systems should be clear.

Because of theorem 2 in appendix A.1, when calculating the second (or third or higher) Lyapunov exponent, we should start with a vector that is orthogonal to the eigenvector of the ?rst Lyapunov exponent (or in the case of higher Lyapunov exponents, we should start with a vector orthogonal to all eigenvectors, belonging to the Lyapunov exponents bigger than the one we want to calculate). This is however not an easy task. For continuous systems, which we have to solve numerically, the situation is even worse because due to numerical integration errors, components from the eigenvector of the largest Lyapunov exponent will become not-zero, and dominate the further evolution. So we have to think of another way of calculating the remaining Lyapunov exponents. I will give the description of this method for the second Lyapunov exponent, but one can easily see from this exposition how the algorithm works for arbitrary Lyapunov exponents. The idea is to compute not the second Lyapunov exponent, but the sum of the ?rst and second Lyapunov exponent, and then compute the second out of our knowledge of the ?rst (for the third for instance we shall have to compute the sum of the ?rst three exponents and then calculate the third from our knowledge of the ?rst and the second). How are we going to calculate the sum of the ?rst and second exponent? By looking at the evolution of two orthonormal error directions, and look at the ampli?cation of the area they span (we shall immediately use the di?erential map to calculate the transformed errors). When assuming that the original vectors have components in the directions of the eigenvectors of the ?rst and the second eigenvalue, the area will behave as An ≈ exp[n(λ1 + λ2 )]A0 (when the vectors also have components in the direction of the other eigenvectors, the largest two will dominate the growth). Thus we have λ1 + λ2 = lim 1 ln(An /A0 ) n→∞ n

4.6

Numerical Results

I calculated all Lyapunov exponents of the systems (42) from section 3 using 10000 iterations and using the di?erential map method. The results are stated in So the algorithm doesn’t change much. We have to table 1. I have also included graphs showing the calculate two transformed directions, and the area convergence towards the Lyapunov exponent for the that they span will be the ampli?cation factor. We largest Lyapunov exponent of the H? enon system and

4

LYAPUNOV EXPONENTS

10

the Lorenz system. In both cases we see that the convergence is very quick. We immediately see that all results are in agreement with the relations we discussed in section 4.2. The small deviations on the exponents that should be zero (for instance for the Lorenz system) can actually be used to estimate the error we make in calculating the Lyapunov exponents. So using this and taking in consideration the ?uctuations we see when we zoom in on the convergence-?gure 1, I would estimate the error for the Lorenz system to be ±0.0005. For the R¨ ossler system the error will be of the same magnitude. And since we don’t make numerical integration errors in the two-dimensional maps, we expect the value for the Lorenz system to be an upper bound on the error of the two-dimensional maps. So I shall use this value for all systems. I will now shortly comment on the individual results. Logistic map The Logistic map is not chaotic for the parameter value r = 3.5699456. We see that it is on the verge of getting chaotic and this is in agreement with the value one ?nds at which the perioddoubling ends. Setting r = 3.9, we see that for this value the Logistic map is chaotic. H? enon map We see that the H? enon map is chaotic. The sum of its two Lyapunov exponents (?1.204) is in very good agreement with its area reducing factor (0.3, we should take the logarithm of this number in order to get the sum of the exponents, ln 0.3 = ?1.20397). Lozi map The Lozi map is also chaotic, and it is even more chaotic than the H? enon map. The sum of its Lyapunov exponents (-0.6931) is in very good agreement with its theoretical value (ln 0.5 = ?0.693147). Zaslavskii map The Zaslavskii map is very chaotic (largest Lyapunov exponent 3.6865). The sum of the Lyapunov exponents should equal the parameter r, and it does so very nicely.

Lorenz system The Lorenz system is chaotic for two of the parameter values I studied. For r = 148 it is not chaotic. In my previous project [23] we saw that the Lorenz system has a periodical attractor at this parameter value. The calculated Lyapunov exponents con?rm this : the largest is zero for the same reason as stated in section 4.2, along the periodical orbit there is on the average no divergence or convergence of nearby trajectories, so one of the exponents has to be zero. Only when we have a ?xed point, none of the exponents will be zero. So this orbit must be a periodical orbit (since we have only three possibilities: ?xed point, periodical orbit or chaotic attractor). In all three cases under study the three Lyapunov exponents add up nicely to ?(1 + p + b) = ?41/3 = ?13.6666 . . .. R¨ ossler system We see that the R¨ ossler system is only slightly chaotic. Adding the Lyapunov exponents has no use, since the divergence of the ?ow isn’t a constant. The sum however does agree nicely with the average of the divergence over the attractor.

4

LYAPUNOV EXPONENTS

11

Logistic Logistic (r = 3.9) H? enon Lozi Zaslavskii Lorenz (r=28) Lorenz (r=148) Lorenz (r=142) R¨ ossler

λ1 ± 0.0005 -0.001 0.4945 0.4189 0.4721 3.6865 0.9051 2.39E-04 1.2533 0.0696

λ2 ± 0.0005

λ3 ± 0.0005

i

λi

-1.6229 -1.1652 -6.6865 8.12E-05 -0.4271 9.20E-05 -2.11E-04

-14.5718 -13.2398 -14.9201 -5.3928

-1.204 -0.6931 -3 -13.667 -13.667 -13.667 -5.3234

Table 1: Lyapunov Exponents. Parameters are as in section 3, unless otherwise indicated.

(a) H? enon

(b) Lorenz (r = 28)

Figure 1: Convergence of ?rst Lyapunov exponent.

5

STRANGE ATTRACTORS

12

5

Strange Attractors

5.1

5.1.1

Dimensions and de?nition of a strange attractor

Topological dimension

The systems described in section 3 all show a strange attractor for certain parameter values. I will de?ne what a strange attractor is in the next subsection. In investigating the properties of these strange attractors, I will focus on their dimension. We shall shortly see that dimension is a broader notion then one might think. The dimension doesn’t have to be an integer for example. Now what is the meaning of the dimension of a strange attractor? In a dissipative system, almost all initial values will eventually settle on an attractor. When we know the dimension of this attractor, we could say that the degrees of freedom the dynamical system has, is essentially this dimension, and this could be signi?cantly less than the dimension of the underlying phase space. In the systems under study here this isn’t very spectacular (twodimensional phase space with 1.5-dimensional attractor for example), but we must not forget that these systems are idealizations which don’t occur in real life. Processes in real life can for example be described by a system of partial di?erential equations, which typically has a in?nite number of degrees of freedom. When we see that these systems exhibit a ?nite-dimensional strange attractor, we can look for a small set of variables which describe the system in its attractor state. Also in experimental situations, one often doesn’t know the dimension of the phase space he’s working in, since one doesn’t know exactly all the variables contributing to a phenomenon. So the dimension of the strange attractor a experimental system may show, is the only thing one knows about the degrees of freedom of the system.

The dimension we are all familiar with is called the topological dimension. A set of disconnected points (by this I mean that they are not in?nitely close together as on a line, which is also just a set of points) has dimension zero. Curves have dimension one, surfaces have dimension two. The topological dimension of an object is always an integer. As a non-trivial example I mention the Cantor set (which one gets by dividing the unit interval in three equal pieces, throwing away the middle piece and iterate this procedure on the two remaining pieces). This set is just a set of points (however uncountable) and thus has a topological dimension zero. 5.1.2 Box-counting dimension

One can de?ne for every object something called its Hausdor?-Besicovitch dimension. The de?nition of this dimension is however rather intricate. It can be found in the appendix of [1]. For most systems it coincides with the Box-counting dimension. This dimension is de?ned as follows: Consider ’boxes’ (in one dimension this would be intervals, in two dimensions squares, in three dimensions cubes, and so on) of side R, then cover the object with these boxes, and then we count the number N (R) of boxes necessary to contain all points of the object. As we let the size of these boxes get smaller, we expect the N (R) to increase. The box-counting is then de?ned as the number Df that satis?es N (R) = lim kR?Df

R→0

(43)

where k is a proportionality constant. We ?nd Db For all these reasons, strange attractors have made thus by taking the logarithm of both sides, before their way into the whole o? the scienti?c world. They taking the limit are used to describe a large variety of systems, for log N (R) log k example in immunology [32] or biology [30]. Df = lim ? + (44) R→0 log R log R Firstly I shall give some explanation on the dimensions involved and then use these dimensions to de?ne Since k is a constant the last term will go to zero, a strange attractor. Then I will describe how I used when R goes to zero, so actually we have these de?nitions to calculate the dimensions of the log N (R) Df = ? lim (45) systems under study and state my results. R→0 log R

5

STRANGE ATTRACTORS

13

In practice we shall verify this law (43) by computing N (R) for an appropriate scaling region, and then plotting ? log N (R) against log R. By using linear regression we can see how well this law is satis?ed and the Df is then given by the slope of the ?tted line. Is this a good de?nition for a dimension? Well it is easy to see that for all Euclidean objects (a point, a line, a plane, . . .) it just gives the topological dimension. I.e. we always need one box to cover a point, so the limit (45) will be zero. In case of a line segment of length L we need N (R) = L/R boxes to cover it. So the limit will be one. 1 ∞ }n=2 . This is just a Consider however the set { n countable set of disconnected points, so one would expect its dimension to be zero. However, I calculated the box-counting dimension in the following 1 way. Consider intervals of length R(R ?1) . This rep1 1 resents the distance between the points R and R? 1. 1 All points i for which i < R will be separated from 1 each other by an amount greater than R(R ?1) . So we need for each of these points an individual interval. This gives us already R ? 2 intervals. To cover the 1/R remaining points we need 1/R (R?1) = R ? 1 intervals. Thus N (R) = 2R ? 3 and we get

1982 Grassberger and Procaccia suggested a new way to de?ne a dimension [29]. Start with a long-time series on the attractor N N {Xi }i=1 ≡ {X(t + iτ )}i=1 where τ is an arbitrary but ?xed time increment. We then de?ne the correlation to be C (r) = lim 1 N →∞ N (N ? 1)

N

Θ(r ? |Xi ? Xj |)

i,j =1 i=j

(47) where Θ is the Heaviside function as de?ned in (15) and we use for example the Euclidean norm (one could also take the maximum norm to speed up the calculations). We then assume that for small r C (r) behaves as follows C (r) ∝ rDc (48)

Dc is then called the correlation dimension. We calculate it using the same procedure as in the previous subsection, choosing an appropriate scaling interval and ?tting a line on the plot of log C (r) against log r. This dimension is sensitive to the distribution of the points on the attractor, since crowded regions will yield an higher correlation. When the distribution of points on the attractor is uniform, the correlation log 2R ? 3 1 log N (R) dimension equals the box-counting dimension. Oth= lim = Df = ? lim 1 R→0 log R(R ? 1) R→0 log erwise it is smaller. We could compare box-counting 2 R(R?1) (46) and correlation dimension with average and variance which obviously isn’t zero. So this can be re- in statistics, the average also doesn’t care much about 3 garded as a failure of the de?nition of the box- the distribution of the values . counting dimension. It can be shown however that the Hausdor?-Besicovitch dimension of this set is 5.1.4 Kaplan-Yorke dimension zero. So the Hausdor?-Besicovitch dimension hasn’t got this problem. Kaplan and Yorke proposed a dimension based on the Lyapunov exponents of the system. Let us rank the Lyapunov exponents from the largest λ1 to the 5.1.3 Correlation dimension smallest λd . Let j be the largest integer such that The box-counting dimension, however relatively easy λ1 + λ2 + . . . + λj > 0, then the Kaplan-Yorke dimento understand, isn’t easy to calculate. In higher di3 actually we should compare the correlation dimension with mensional spaces there is so much more space, so we the second moment, since the ?rst moment, the variance in need much more boxes, and thus much more points statistics, can be compared with yet another dimension, the on our attractor, which are mostly calculated by a information dimension. This dimension is part of a bunch computer or come from experimental data. And cal- of generalized dimensions which can be compared with the culating these points can take a lot of time. So in arbitrary moments of a distribution in statistics.

5

STRANGE ATTRACTORS

14

Now there are theorems (Takens, 1981; Ma? ne, 1981 which state that we can embed see for instance [21]) j i=1 λi a one-dimensional time series in a high dimensional DKY = j + (49) ?λj +1 space (typically twice the Hausdor? dimension), thus obtaining a good projection of the attractor, which Kaplan and Yorke conjectured that, for a two- has the same properties as the one which is under dimensional mapping, the box-counting dimension investigation. Df equals the Kaplan-Yorke dimension DKY [24]. We construct this from our one-dimensional time This was subsequently proven to be true in 1982. A series {x } N as follows: i i=1 later conjecture held that the Kaplan-Yorke dimension is generically equal to another dimension called xi = (xi , xi+tL , xi+2tL , . . . , xi+(d?1)tL ) the information dimension, which is also closely related to the box-counting dimension and the correla- where tL is some appropriate time lag (which is not tion dimension. This conjecture is partially veri?ed so easy to choose [16]) and d is the embedding dimenby Ledrappier (1981). (This paragraph partially from sion. Calculating dimensions and also Lyapunov expo[9]). nents from experimental data is however much more di?cult than from the systems we study. One should 5.1.5 De?nition of a strange attractor for example use enough points in order to get a corWe can now de?ne what we mean by a strange attrac- rect result [22, 6]. tor. We shall say that an attractor is strange when its Hausdor? dimension strictly exceeds the topolog5.2 Algorithms for calculating dimenical dimension [1] (in dynamical systems we speak of sions strange attractors, otherwise these objects are called fractals ). So a strange attractor is an attractor with 5.2.1 Box-counting dimension a non-integer Hausdor? dimension. For two-dimensional attractors the box-counting dimension is relatively easy to calculate. For three5.1.6 Concerning dimensions from experidimensional systems such as the Lorenz system, this mental data calculation already becomes much more complicated Before turning to the description of the algorithms (for reasons stated in section 5.1.3). used to calculate dimensions, I would like to say As a ?rst attempt I just made an algorithm purely something about experimental situations. As I said in from the de?nition. I divided the region of phase the introduction to strange attractors, one wants to space in squares (in the two-dimensional case) and ?nd the dimension of strange attractors observed in looked for every square wether a point of the attracan experimental situation. Now when we look at the tor was in it or not. Repeating this for smaller and de?nitions of the dimensions, we see that we have to smaller squares, I could then ?nd the box-counting know in what phase space we’re working to be able dimension from linear regression described in section to calculate the dimensions. When calculating the 5.1.2. Realizing that I was doing a lot of work that box-counting dimension we have to use boxes. Now wasn’t necessary (in every step I checked for every these boxes have the dimension of the phase space square wether a point was in it), I programmed a we’re working in (this is also partially solved in the recursive algorithm. This recursive algorithm begins Hausdor?-Besicovitch dimension). In calculating the with the region the attractor is contained within. If correlation dimension we also use a time series made a point is in it (which obviously is true) it gets split up of points in our phase space. In experimental sit- up in for smaller squares. Then the same algorithm uations however we can only measure one, or a few does the same for these squares. So it checks for the variables. four squares wether a point is in it and if so it splits

sion is de?ned to be

5

STRANGE ATTRACTORS

15

it further up, until some bottom level, which can be speci?ed, is reached. These algorithms are very slow. They also use a enormous amount of memory: applying them on the Lorenz system didn’t work because my computer didn’t have enough working memory to hold the information of all the boxes. I then came up with a better idea. Instead of looking at all the squares, I just looked at all the calculated points of the attractor. So I ‘hypothetically’ divided the region of phase space where the attractor lies (by ‘hypothetically’ I mean that the computer didn’t have to do this explicitly). Then I looked for every point in which square it lies and then assigned to that point a number, uniquely determined by the square. Doing this for all the calculated points, I then had an array of numbers. Of course when two points are within the same square, they will get the same number assigned to them. So this array contains doubles. Getting rid of these doubles leaves us with an array which length is the number of squares needed to cover the whole attractor. This program works a thousand times as fast as the other (this is just an estimation, not a rigourously checked statement, it isn’t however an exaggerated guess).

already calculated the Lyapunov exponents of all systems). For the two-dimensional maps I calculated the box-counting dimension. For the three-dimensional continuous systems I calculated the correlation dimension. I can immediately note that all results are independent of the initial value. For the calculation of the box-counting dimension I checked this for 160 or 200 di?erent initial values. For the correlation dimension I only checked ?ve initial values, because of the computing time required. All the initial values gave the same result, with, especially for the maps, very small variance. For the box-counting dimension of the twodimensional maps I would estimate the error to be ±0.01 since this is how much the dimension changes when I change the scaling interval without getting a bad ?t. For the Lorenz system and the R¨ ossler system I would estimate this error to be larger. The scaling interval I used here is very small (for the Lorenz system 0.1:0.025:0.2). This gives the best ?t, but changing the interval doesn’t make the ?t much worse, the slope however changes signi?cantly. So the error should be at least ±0.05 and maybe even ±0.1. All results (including the Lyapunov exponents) are given in table 2. In each case the Kaplan-Yorke di5.2.2 Correlation dimension mension is in good agreement with our results. When Calculating the correlation dimension was pretty consulting the literature [20, 29, 2, 3] we see that all straightforward. I just implemented the de?nition dimensions are in rather good agreement with the (47). The only thing to note here is that I ad- literature results. justed the de?nition a little, to prevent autocorrelation. Calculating the points of a continuous system yields that for small stepsize, successive points will automatically be close together. So they will surely be correlated, exaggerating the correlation. So instead of just i = j , I demanded |j ? i| > k for some suitable k (I used 50).

5.3

Results

Since the amount of computer time necessary to get good dimension estimates was enormous, I calculated for every strange attractor only one dimension (for every system I also calculated the Kaplan-Yorke dimension, which obviously wasn’t much work since I

5

STRANGE ATTRACTORS

16

H? enon (a=1.4, b=0.3) Lozi (a=1.7, b=0.5) Logistic (r = 3.5699456) Logistic (r = 3.9) Zaslavskii (ν = 400, r =3, a = 12.6695) Lorenz (p=10, b=8/3, r=28) Lorenz (p=10, b=8/3, r=148) Lorenz (p=10, b=8/3, r=142) R¨ ossler (a=0.2, b=0.2, c=5.7)

Df ± 0.01 1.251 1.365 0.5462 0.9118 1.479

Dc ± 0.05

λ1 0.4189 0.4721 -0.001 0.4945 3.6865

λ2 -1.6229 -1.1652

λ3

DKY 1.258118 1.405166

-6.6865

1.551335

2.0519

0.9051 2.39E-04 1.2533

8.12E-05 -0.4271 9.20E-05 -2.11E-04

-14.5718 -13.2398 -14.9201 -5.3928

2.062

1.9729

0.0696

2.013

Table 2: Results

REFERENCES

17

6

Conclusion

[5] Edward Ott. Chaos in dynamical systems. Cambridge University Press, 1993.

When calculating the Lyapunov exponents and di[6] Edward Ott, Tim Sauer, and James A. Yorke, mensions of the strange attractors of the systems I editors. Coping With Chaos. Wiley Series in introduced in section 3, we see that all methods work Nonlinear Science. John Wiley & Sons, Inc., fairly good for all systems. Only calculating the box1994. counting dimension of the three-dimensional systems didn’t work properly, but this has more to do with [7] E.N. Lorenz. Deterministic non-periodic ?ow. J. the limited computer time than with a failure of the Atmos. Sci., 20:130–141, 1963. algorithm. I also noticed clearly that doing calculations with two-dimensional discrete maps was a lot [8] Eric W. Weisstein. H? enon map. From easier and quicker than doing the same calculations Mathworld–A Wolfram Web Resource. http: for three-dimensional continuous systems, which was //mathworld.wolfram.com/HenonMap.html. to be expected. One can ?nd my ?nal results in table 2. One task [9] Eric W. Weisstein. Kaplan-yorke conjecfor the future could be ?lling in the blanks of this ture. From MathWorld–A Wolfram Web table (note that some are blank because they have Resource. http://mathworld.wolfram.com/ no meaning for a speci?c system, like the second and Kaplan-YorkeConjecture.html. third Lyapunov exponent of the Logistic map). More Logistic map. From recently attempts have already been made to general- [10] Eric W. Weisstein. 7 Mathworld–A Wolfram Web Resource. http:// ize some of these results. Investigating up to 4 ×10 mathworld.wolfram.com/LogisticMap.html. low-dimensional, low order polynomial maps and ordinary di?erential equations, J.C. Sprott has found that typically a few percent of them is chaotic [18]. [11] Eric W. Weisstein. Lozi map. From Mathworld– A Wolfram Web Resource. http://mathworld. By investigating the strange attractors which these wolfram.com/LoziMap.html. chaotic maps show, he found a relation between the dimension of the phase space and the dimension of [12] Eric W. Weisstein. Zaslavskii map. the strange attractor [19]. From Mathworld–A Wolfram Web ReOne could go further on this road and thus obtain source. http://mathworld.wolfram.com/ a good estimator of the dimension of a strange atZaslavskiiMap.html. tractor. [13] G.M. Zaslavskii. The simplest case of a strange attractor. Physics Letters A, 69:145–147, 1978.

References

[1] Benoit B. Mandelbrot. Fractals: Form, Chance and Dimension. W.H. Freeman and Company.

[14] Heinz-Otto Peitgen, Hartmut J¨ urgens, and Dietmar Saupe. Chaos and Fractals: New Frontiers of Science. Springer-Verlag, 1992.

[2] David A. Russell, James D Hanson, and Edward [15] Robert C. Hilborn. Chaos and Nonlinear DyOtt. Dimension of strange attractors. Physical namics. Oxford University Press, 1994. Review Letters, 45(14), 1980. [16] H.S. Kim, R. Eykholt, and J.D. Salas. Nonlinear [3] Divakar Viswanath. The fractal property of the dynamics, delay times and embedding windows. lorenz attractor. Physica D, 190:115–128, 2004. Physica D, 127:48–60, 1999. [4] Edward N. Lorenz. The essence of chaos. UCL Press, 1993. [17] J. C. Sprott. Simplest dissipative chaotic ?ow. Physics Letters A, 228:271–274.

A

OSELEDEC’S MULTIPLICATIVE ERGODIC THEOREM

18

[18] J. C. Sprott. How common is chaos. Physics Letters A, 173:21–24, 1993.

[19] J. C. Sprott. Predicting the dimension of strange attractors. Physics Letters A, 192:355–360, [31] Colin Sparrow. The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors. Springer1994. Verlag, 1982. [20] J. C. Sprott. Improved correlation dimension calculation. International Journal of Bifurcation [32] Thomas Schall. Fractalkine - a strange attractor in the chemokine landscape. Immunology Today, and Chaos, 11(7):1865–1880, 2001. 18(4):147, 1997. [21] J.-P. Eckmann and D. Ruelle. Ergodic theory of chaos and strange attractors. Reviews of Modern [33] V. I. Oseledec. A multiplicative ergodic theorem, ljapunov characteristic numbers for dynamical Physics, 57(3):617–656, 1985. systems. Trudy Mosk. Obsch., 19, 1969. [22] J.-P. Eckmann and D. Ruelle. Fundamental limitations for estimating dimensions and lyapunov exponents in dynamical systems. Physica D, 56, 1992. [23] Jo Bovy and Mark Cox. Chaos en het lorenz systeem. Theoretisch Projectwerk, Eerste licentie Natuurkunde, Katholieke Universiteit Leuven, 2004.

[30] Sakire Pogun. Are attractors ’strange’, or is life more complicated than the simple laws of physics. BioSystems, 63:101–114, 2001.

A

Oseledec’s multiplicative ergodic theorem

In this appendix I will state Oseledec’s multiplicative [24] Kaplan, J.L. and Yorke, J. A. Functional Difergodic theorem [33] and it’s consequences concerning ferential Equations and Approximations of Fixed Lyapunov exponents. I have taken the statement of Points, page 204. Springer-Verlag, 1979. this theorem from [21]. [25] M. H? enon. A two-dimensional mapping with a strange attractor. Commun. Math. Phys., 50:69– A.1 The theorem 77, 1976. Theorem 1 (Continuous-time Multiplicative Er[26] O.E. R¨ ossler. An equation for continuous chaos. godic Theorem). Let ρ be a probability measure on a space M , and f : M → M a measure preserving Physics Letters A, 57:397–398, 1976. map such that ρ is ergodic. Let also T : M → the [27] P.-F. Verhulst. Recherches math? ematiques sur la m × m matrices be a measurable map such that loi d’accroissement de la population. Nouv. mm. de l’Academie Royale des Sci. et Belles-Lettres ρ(dx) log+ T (x) < ∞, de Bruxelles, 18:1–41, 1845.

+ [28] P.-F. Verhulst. Deuxi` eme m? emoire sur la loi where log (u)=max(0,log(u)). De?ne the matrix n n?1 x) · · · T (f x)T (x). Then, for ρ-almost d’accroissement de la population. Mm. de Tx = T (f all x, the following limit exists: l’Academie Royale des Sci., des Lettres et des Beaux-Arts de Belgique, 20:1–32, 1847. n? n 1 / 2 n lim (Tx Tx ) = Λx (50) n→∞ [29] Peter Grassberger and Itamar Procaccia. Charn? n acterisation of strange attractors. Physical Re- (We have denoted by Tx the adjoint of Tx , and n? n taken the 2nth root of the positive matrix Tx Tx ) view Letters, 50:346–349, 1983.

B

COMPUTER PROGRAMS

19

The logarithms of the eigenvalues of Λx are called characteristic exponents. These are just the Lyapunov exponents as de?ned in (29). They are ρalmost everywhere constant. So the Lyapunov exponents are not dependent on the initial value, only in a subset of ρ-measure zero can they be di?erent. Let λ(1) > λ(2) > · · · be the characteristic exponents, no longer repeated by multiplicity; we call m(i) ( i) the multiplicity of λ(i) . Let Ex be the subspace of Rm corresponding to the eigenvalues ≤ exp λ(i) of (1) (2) Λx . Then Rm = Ex ? Ex ? · · · and the following holds Theorem 2 For ρ-almost all x,

frequencies the natural measures of the boxes. For a typical x0 in the basin of attraction, the natural measure of a typical box Ci is ?i = lim η (Ci , x0 , T ) T →∞ T (52)

where η (Ci , x0 , T ) is the amount of time the orbit originating from x0 spends in Ci in the time interval 0 ≤ t ≤ T.

B

Computer Programs

The implementation of the algorithms given in the text were all written as m-?les in Matlab. The source 1 code for all these programs can be found on my webn lim log Tx u = λ ( i) (51) site http://m0219684.kuleuven.be (click on Matlab n→∞ n Programs). Here you can also ?nd some pictures of (i+1) ( i) . In particular, for all vectors u the strange attractors. if u ∈ Ex \Ex (2) that are not in the subspace Ex , the limit is the (1) largest characteristic exponent λ .

When we randomly choose a vector u, it will most probably not lie in one of the special subspaces ( i) Ex i > 1. So this theorem then says that we will get the largest Lyapunov exponent when calculating the limit. To get the ith exponent, we have to carefully (i+1) ( i) . This select a vector u such that u ∈ Ex \Ex makes the computation of the other Lyapunov exponents very di?cult.

A.2

A measure on the attractors

In order to apply these theorems to the strange attractors we encounter in the dynamical systems under study in this project, we need to make sure that there is a measure on these attractors. This measure can be de?ned as follows. We can cover the attractor with a grid of boxes (as in computing the box-counting dimension) and then look at the frequency with which typical orbits visit the various boxes covering the attractor in the limit that the orbit length goes to in?nity. When these frequencies are the same for all initial conditions in the basin of attraction of the attractor except for a set of Lebesgue measure zero, then we call these