Bifurcations and strange attractors in a climate related system Текст научной статьи по специальности «Математика»

r

dx

DIFFERENTIAL EQUATIONS AND

CONTROL PROCESSES N 1, 2005 Electronic Journal, reg. N P23275 at 07.03.97

dt

http://www. neva. ru/journal e-mail: diff@osipenko. stu. neva. ru

Ordinary differential equations

BIFURCATIONS AND STRANGE ATTRACTORS IN A CLIMATE RELATED SYSTEM

Niklas Lundstrom

Department of Mathematics, Lulea Univeristy of Technology,

Lulea, SWEDEN, e-mail: niklas.lundstrom@ltu.se

Abstract.

R. V. Bekryaev derived a system for a horizontally baroclinic atmosphere consisting of six ordinary differential equations. We prove dissipativity and find estimates for the location of the global attractor. The evolution of a complicated attractor is analysed with a Poincare map showing difficult bifurcation behaviour. Investigations in bifurcation diagrams show a rich dynamical behaviour including a lot of known complicated bifurcations, of which a fold-Hopf bifurcation is examined in detail. Finally, we give some theory about the Lya-punov spectra and present a method for determining the exponents.

1 Introduction

The equations for atmospheric flow are one of the most studied dynamical systems of the last century. Both their practical importance and their mathematical richness have attracted much attention. R. V. Bekryaev at the Arctic and Antarctic Research Institute of Saint Petersburg derived a model for a horizontally baroclinic atmosphere. In the derivation of the system a mean approximation of the atmosphere was used, this means that the atmosphere was treated as a thin layer at some height above the surface of the Earth. The Navier-Stokes equation, the equation for heat conduction, the continuity equation and a Galerkin procedure were then used to obtain the continuous system consisting of six ordinary quadratic differential equations given by

d X Y f G\

= -BUT - PrX - BeG + 2.2GH + (Uv - H)— + [A + B + — ]Z

or 50 \ 50/

d G F ( X\

— = Bex - PrG - 2.2XH + (H - U„)- + [C - -J Z

^ = _ p + (2L_£W (1)

dr \ P 160 7 V160 PJ

d F

= Qf ~ UTG — F + (Uu — H)Y + GZ = Qy + UTX + (H - UU)F - Y - XZ ^-GF+XY-PZ

dr

dY

dr

dZ_ dr

Here t is a dimensionless time, A, B, Be, C, P, Pr, QF, QY, UT and Uu are parameters of which P > 0 and Pr > 0. We will refer to Bekryaev [1] for a derivation of the system and explanations to physical relations to the parameters and variables. More theory about dynamic meteorology can be found in [2]. Note that, in the three-dimensional modification of the Bekryaev model, G = F = H = 0, a linear change of variables, X = X, Y = and Z = allows us to reduce the system to the famous Lorenz equations

dX

— = PrY - PrX or

d YP Pe2 -

= -XZ + X-y

d t 2000Pr

or 5

where Pe is the Pecle number. This system was found by Edward N. Lorenz in 1963, see [3]. In 1984 Lorenz published a modification of his system, referred to as the Lorenz-84 model, see [4]. This system has been of great interest for mathematicians and yields

^ = -aX -Y2 - Z2 + aF

or

d Y

— = -Y + XY - bXZ + G

or

d Z

— = -Z + bXY + XZ

or

where a and b are parameters. In a later paper [5], Lorenz pointed out that F and G should be allowed to vary periodically, and in 2002, H. Broer, C. Simo and R. Vitolo published an article [6] where they examined this system for

F F(1 + e cos(^t))

G ^ G(1 + e cos(^t))

Except for the Bekryaev system, there is another six-dimensional system examined by L. Van Veen [7] in 2003, with quite analogous behaviour as the Bekryaev system. For example, both systems give rise to fold-Hopf bifurcations and the dynamics seems to take place on three-dimensional invariant manifolds in both systems. Also, L. Van Veen found a reduction of his six-dimensional system to the Lorenz-84 model. It is interesting that all of these systems have so much in common, even if they represent different models.

In Section 2 we will prove the existence of a global attractor by showing that the Bekryaev system (1) is dissipative by finding a suitable Lyapunov function.

After this proof, seven of the parameters are given values found in [1], and are from now on fixed when the system is investigated by changing the three parameters A, B and QF. The values are

Be = -8.5242, C = 0, P =16/5, Pr = 1 (2)

Qy = 0, Ut = 620.15, Uu = 42.467

Estimates for the bounds of the global attractor are found for various values of the parameter QF for the case when A = 2 and B = 0.

In Section 3 one attains knowledge about some non-wandering sets that occur when changing QF for A = 2 and B = 0. The system is simulated using a Runge-Kutta method programmed in C-code. A complicated attractor that arises close to a fold bifurcation is discussed and analysed using a Poincare map, showing difficult bifurcation behaviour. The Lyapunov spectrum is also determined for this attractor at different values of QF and compared to the results from the Poincare map.

Bifurcations in the plane of the parameters A and B will be analysed in Section 4 using the software package AUTO 97, which can be found on the internet at [8]. The fold and Hopf bifurcation curves found in the AB-plane and different regions and their corresponding set of equilibria are discussed. Two bifurcations of codimension two are discovered, the Bogdanov-Takens and the fold-Hopf bifurcation. The fold-Hopf bifurcation will be studied in more detail. An approximation of the locally defined three-dimensional center manifold will be found and the system is transformed into a normal form for fold-Hopf bifurcations derived in [10], from which the behaviour near the bifurcation can be analysed. Some of the results are then verified by further simulations with AUTO. In the end of the section some attractors are shown, one of which shows an invariant torus.

The definition of the Lyapunov spectrum and some facts about the exponents can be found in Section 5, and a method for computing the spectrum is presented in which one augment the dynamical system with an orthonormal frame and a Lyapunov vector. The implementation is also described and some results are shown.

2 Dissipativity

In this section dissipativity will be proved for the Bekryaev system (1). Before giving the definition of dissipativity, we introduce the evolution operator which determines the state xt of a system at time t, provided the initial state x0 is known. Assume that for each t > 0 a single-valued map ^ is defined in the phase space X,

^ : X ^ X

which transforms an initial state x0 £ X into some state xt £ X at time t > 0 given by

xt = ^xo

We are now ready to give the following definition

Definition 2.1 An n-dimensional dynamical system is said to be dissipative if there exists a bounded region S £ Rn such that V x0 £ Rn 3 to > 0 such that Vt > t0, (£>tx0 £ S.

In words, the above definition means that, in a dissipative dynamical system there exists a region such that all trajectories enter this region after sufficiently long time, and nothing can escape from this region. This means that if we can prove dissipativity, we have proved the existence of a global attractor. It is therefore of great interest knowing if a dynamical system is dissipative or not.

2.1 Proof of dissipativity

The dissipativity of the Bekryaev system is given by the following.

Theorem 2.1 The Bekryaev system (1) is dissipative.

Proof

The idea of the proof is to find a proper Lyapunov function V, and then prove that it is possible to take a constant r such that all regions where V is increasing in time are totally contained in the ellipsoid formed by V = r. This means that outside the ellipsoid V = r the Lyapunov function V is decreasing, and therefore, trajectories outside V = r tend inward this ellipsoid. Due to symmetries in the Bekryaev system (1), we get rid of some improper terms in the time derivative by considering as Lyapunov function

V = ^X2 + G2 + y (H - Uu)2 + F2 + Y2 + {Z - UT)2^J (3)

The level surface to (3) given by V = r forms an ellipsoid centered at the point X = G = F = Y = 0, H = Uu, Z = UT. By taking the derivative with respect to time t of the Lyapunov function (3), and substituting the derivatives from (1), we obtain

V Vlinear + Vquad + Vcross (4)

where Vlinear, Kuadl and Vcross are given by

VLear = -BUTX + fpPrUuH + (V + B)^ + Qf^J F

flQCUu ~hp

+ + Qy ] Y + PUTZ

Vquad = -PrX2 - PrG2 - fpPrH2 -F2-Y2- PZ2

16 1

VLss = -(A + B)—HF - —HY + CGZ + (A + B)XZ (5)

Since P > 0 and Pr > 0, the part Vqiuad is always negative and dominates Viinear far from the origin. The problem is the cross terms in V'cross, which can be positive even far away from the origin. To get rid of these terms, we consider another function Vfyz in the FYZ-subspace given by

Vfyz = \{F2 + Y2 + {Z- UT)2) (6)

Taking the derivative of (6) with respect to t gives

v^yz = fqf - f2 + yqy - y2 - pz2 + utpz (7)

Clearly, expression (7) is a quadratic form only in the variables F, Y and Z. Further, outside of the sphere V'YZ = 0, we have V'YZ < 0 due to the negative signs of the second order terms. This means that we can take a constant r such that outside the sphere VFYZ = r the function VFYZ is decreasing in time. This means that we have shown dissipativity in the FYZ-subspace.

Observe that all terms in (5) involve one of the variables F, Y or Z. Therefore, due to the dissipativity in the FYZ-subspace, bounds for the cross terms can be found by replacing F, Y and Z by sufficiently large constants, say

Pmax, Ymax and Zmax such that

1 F \ — Fmax 1 Y 1 — Ymax 1 Z \ — Zmax

From this and the triangle inequality we can find a bound for expression (5) by writing

16 16 VLss < 5P I (A + B)HF I +— I CHY I + I CGZ | + | (A + B)XZ | <

16 16

— | (A + B)H | Fmax + — | CH | ymax+ | CG | Zmax+ | (A + B)X | Zmax 5P 5P

Thus, VCross is less than a function including only linear terms, and it follows that the derivative V' given by (4) becomes negative far away from the origin. Therefore, it is possible to take a constant r such that all regions where the Lyapunov function given by (3) is increasing in time are totally contained in the ellipsoid given by the level surface V = r, and we have proved the theorem.

2.2 Locating the attractor for a special case

In this section, estimates for the location of the global attractor will be found at different QF values for the case when A = 2, B = 0. The rest of the parameters

are fixed to the values given in (2). This case was studied by Bekryaev in [1] and will be analysed further in Section 3. By giving QF a value, we will be able to determine the smallest possible value of the constant r defined in Section 2.1, say rs, giving the smallest boundary of the global attractor. To do so, consider the derivative with respect to time t of the Lyapunov function given by (3). In this case it takes the simpler form

V' = -X2 - G2 - F2 - Y2 + 2(ZX - HF + UuF)

256 16

■(H2 + HUu) - — (Z2 + ZUt) + FQf (8)

25 5

Expression (8) is negative outside the ellipsoid V' = 0. In order to find rs, we will determine a point on both surfaces, that is, the ellipsoids V' = 0 and V = rs, of which the corresponding normals nV and nV' are parallel. From this we have the equations

n V = a n V' a > 0

V = r (9)

V' = 0

where the normals to the level surfaces can be calculated from

nV = ±VV and nV' = ±VV'

The first equation in (9) is six-dimensional, and since we have eight unknown, this is enough for finding solutions with the software package MAPLE. In this case however, we obtained three real solutions, which is natural because the ellipsoids may touch several times with parallel normals when rs is increased from zero. But, choosing the solution that gives the largest rs ensures that the ellipsoid V' = 0 is totally contained in the ellipsoid V = rs. To illustrate the result, we look for the maximum value each of the variables X, G, H, F, Y, Z can attain on the global attractor. From the Lyapunov function (3) we find by setting all except one variable equal to zero the following bounds for the variables

■V^fl < X,G,F,Y < V^s

4 4

UT ~ VZTs < Z <UT+ v/2r\

The value of rs is determined for fixed parameter values. Thus, to illustrate how the above intervalls change varying QF, we repeat the calculations for the 401 different parameter values, Qf = —3000, —2990,..., 1000. The result is shown in figure 1, the curves plotted are the functions /i(Qf) = Ut + \/^s{Qf) and /2{Qf) = — \/2rs{Qf)- Thus, on the global attractor, all variables are bounded between these curves.

Figure 1: The functions fi(QF) and f2(Qf) which constitute the bounds for all variables on the global attractor for -3000 < QF < 1000.

3 Analysis in the parameter QF

In this section, different non-wandering sets that occur in the Bekryaev system (1) when changing the parameter QF will be studied. The parameter values are A = 2 and B = 0, the rest of the parameter values are given in (2). In Section 2 it was proved that the Bekryaev system has a global attractor, and estimates for the bounds of this attractor were found especially for this case, see Section 2.2. With these parameter values, the Bekryaev system yields

Y f C

X' = -X- BeG + 2, 2GH + {Uu - H)— + 2 + —

5U \ 50

F 1

G = BeX -G- 2, 2 XH + (H — Uu)---XZ

50 50

50

Z

(10)

F' = qf - utC - F + (uu - H) Y + CZ

Y' = UTX + (H - Uu)F - Y - X

Z' = -GF + XY--Z

5

. d -X

Here we have used the simpler notation X = ^ and so on. 3.1 Finding equilibria

In this section we will find equilibria to system (10) by solving the equations obtained by requiring the derivatives to be zero. The software package MAPLE was used in order to accomplish the following calculations. Unfortunately, the package was not able to solve the system of equations immediately, therefore symmetries in the system were used to make simplifications. By taking the following linear combinations

of the equations in (10) we obtain three new equations which contain less quadratic terms than the previous ones. These expressions yield, after substituting the derivatives from (10),

V' = 160H' - Z'

W' = 50X' - F'

S' = 50C' — Y'

V' = -512H - 100F + —Z

5

W' = -50(X + BeC) + 110CH + 100Z - Qf + UtG + F S' = 50(BeX - C) - 110XH - UTX + Y

Solving V' = 0 for H and substituting H into the above expressions W', S', and into F', Y' and Z' from (10) gives the five differential equations

(95 F Z \

-—--— + 100Z - Qf + UTG + F

190 160/

(25F Z \

128 "TeoJ + y

(25F Z \

T28 " 160 ) +GZ

( 25F Z \

y' = + F---- FUu - Y - XZ

V 128 160 J 16

Z' = -GF + XY--Z

5

After solving Y' = 0 for Y, and S' = 0 for G, we still have three expressions, W', F', Z', each involving the variables X,Z and F. In W' = 0 and F' = 0, X is found from first order equations, therefore one easily obtain two expressions for X as functions of F and Z, say Xi and X2. In the following step, substitute one of these expressions into Z'. By plotting the curves X2 = X1 and Z' = 0 for Qf = 0 in the same plot, one can graphically determine intervals small enough for MAPLE to find the intersections of the curves.

For Qf = 0, we have found apart from the unstable equilibrium at the origin (Fo), two equilibria, one of which is unstable (Fi), and one stable (F2). Unfortunately, we cannot with this method show that these three equilibria are the only existing ones.

3.2 Analysing equilibria varying QF

The software package AUTO was used in order to examine the behaviour of the three equilibria F0, Fi and F2 when changing QF. In figure 2, the distance from the origin to all three equilibria and the maximum distance from the origin to a limit cycle marked P are plotted versus QF. F0 moves slowly when Qf changes, this unstable equilibrium has four contracting directions and no bifurcation occurs. The two other equilibria, Fi and F2 are created in a fold bifurcation at QF = 200. For large QF values we have the stable limit cycle P as an attracting set, the cycle is shown for QF = 300 in figure 5; Section 3.3. This cycle loses stability in a Neimark-Sacker bifurcation at QF = 185, and disappears in a subcritical Hopf bifurcation for QF = 172 at equilibrium Fi.

Figure 3 illustrates the behaviour of Fi in each coordinate. The fold bifurcations are marked by dotted lines and the Hopf bifurcations by broken lines. The first branch of equilibrium Fi (Fia) exists in the interval -61.0 < QF < 200, the number of contracting directions changes from five to three in the sub-critical Hopf bifurcation at QF = 172. Moreover, at QF = -10.3 we have a fold bifurcation where two more equilibria, one stable (Fib) and one unstable (Fic) are created. Thus, from Fi, there are three equilibria in the interval -61.0 < Qf < -10.3. In this interval the stable equilibrium Fib loses two contracting directions in a subcritical Hopf bifurcation at QF = -28.6, where also an unstable cycle is created. Further, the branches Fib and Fia meet and disappear in another fold bifurcation at QF = -61.0. Thus, for QF < -61.0 only Fic remains and changes slowly.

F2 consists of one stable equilibrium (F2a) when -14.4 < QF < 200. Similar to Fi , two more equilibria, one unstable branch (F2b) and one stable branch (F2c) arise in a fold bifurcation at QF = -14, 4. The stable branch F2c loses two contracting directions in a supercritical Hopf bifurcation at QF = -14, 7, where a stable limit cycle occurs. Further, F2a and F2b are destroyed in a fold bifurcation at Qf = -50.7, and for QF < -50.7 only F2c remains and changes slowly. Figure 4 illustrates the existence of equilibria in different intervals from F2. Table 3.1 summarizes the number of equilibria existing in different intervals of QF.

1400

1200

1000

400

200

0

-1000

-500

500

0

Q

F

Figure 2: The distance from the origin to the three equilibria F0, Fi, F2 and the maximum distance from the origin to a limit cycle marked P for -1000 < QF < 500.

14

12

10 F1a

* 8

6 Fb

Fc

4

-200 -100 0 100 200

qf

qf

Figure 3: The variation of Fi in all coordinates for -200 < QF < 200, dotted lines show the location of fold bifurcations and broken lines indicate the location of Hopf bifurcations.

Figure 4: The variation of F2 in all coordinates for -200 < QF < 200, dotted lines show the location of fold bifurcations and the broken line indicates the location of a Hopf bifurcation.

Table 3.1: Equilibria in different intervals of QF. The number of contracting directions of each equilibrium are given by the number in each square.

Qf No F0 Fla Fib Flc F2a F2b F2c

(1000,200) 1 4 - - - - - -

(200,172) 3 4 5 - - 6 - -

(172,-10.3) 3 4 3 - - 6 - -

(-10.3,-14.4) 5 4 3 6 5 6 - -

(-14.4,-14.7) 7 4 3 6 5 6 5 6

(-14.7,-28.6) 7 4 3 6 5 6 5 4

(-28.6,-50.7) 7 4 3 4 5 6 5 4

(-50.7,-61.0) 5 4 3 4 5 - - 4

(-61.0,-3000) 3 4 - - 5 - - 4

3.3 Attractors at different QF values

In this section we will look for different attractors that occur in system (10) when changing QF and taking a starting point close to the origin. For large QF values, the stable cycle P discussed in Section 3.2 is the attracting set, see figure

5. By decreasing the parameter, cycle P disappears, and the stable equilibrium F2a takes over.

H Y

Figure 5: The limit cycle P of system (10) for QF = 300.

A more complicated attractor (CA) occurs when F2a destroys in the fold bifurcation at Qf = -50, 7. Figure 6 shows this attractor for QF = -300. The Z-coordinate for the five different parameter values QF = -300, -450, -650, -1000 and QF = -2500 is shown in figure 7 to give an overview of how the attractor changes with the parameter. For QF = -300, the attractor CA clearly shows one repeating laminar and one repeating more complicated part. In fact, the laminar part becomes longer by increasing the parameter until we have the stable equilibrium. In order to explain this behaviour, consider the following. Figure 8 illustrates the location of the attracting set for QF = -50.7 and Qf = -51.0. In the first case, trajectories are attracted by the stable equilibrium F2a, and in the latter, we have the attractor CA. In figure 8 one can see that the laminar part lies close to the position of the stable equilibrium, and we may guess, due to this and figure 7, that the attracting behaviour of F2a remains even for small QF values, where the equilibrium does not exist. To see that this is the case, the local linearization for the equilibrium F2a are evaluated at Qf = -50.7, that is, we calculate the eigenvalues and eigenvectors to the Jacobian matrix of system (10) at that point. Now, consider the eigenvector corresponding to the eigenvalue which is zero in the fold bifurcation where F2a disappears. By comparing this eigenvector to the direction of trajectories in the laminar part of CA for QF = -51.0, we conclude that trajectories leave the

laminar part in a direction close to this eigenvector. Therefore, it is clear that the laminar part occurs due to the earlier existing stable equilibrium. Moreover, the behaviour of the more complicated parts of the attractor CA remains for some parameter values larger than QF = -50.7. This can be seen by choosing a starting point at each side of equilibrium F2b, (which has five contracting directions, cf table 3.1) in the repelling direction. From both starting points, the trajectory arrives F2a, and one of these orbits behaves similar to the more complicated part of CA for nearby QF values. From this we conclude that the attractor CA creates when the branches F2a and F2b meet and disappear in the fold bifurcation.

Figure 6: The attractor CA of system (10) for QF = -300. In this case the attractor shows a limit cycle.

3.4 Analysing CA with a Poincare map

In this section we will examine the behaviour of the attractor CA when changing Qf, that is, we will investigate whether CA shows some periodic cycles or chaotic behaviour. For doing this, a Poincare map1 Pqf , depending on QF is defined.

PQf : R5 ^ R5

1Theory about Poincare maps can be found in [9], [10] and [11].

1000

500 -

(a)

50 55 60 65 70 75 80 85 90 95 100

10001-1-1-1-1-1-1-1-1-1-

AiL .»ill ___.AMLL .AALL .diiL ___M

500 - -

(b)

50 55 60 65 70 75 80 85 90 95 100

1000

500 -

0

- -(c)

55 60 65 70 75 80 85 90 95 100

1000 500 0 1000 500 0

(d)

55 60 65 70 75 80 85 90 95 100

(e)

100

Figure 7: The Z-coordinate of attractor CA when (a) QF = -300, (b) QF = -450, (c) QF -650, (d) Qf = -1000, (e) QF = -2500.

200 400 600 800

1000 1200 t

1600 1800

Figure 8: The distance from the origin to the stable equilibria F2a at QF = -50.7 as a broken line, and to the attractor CA at QF = -51.0 as a solid curve.

3.4.1 Construction of the Poincare map

Consider a five-dimensional plane nearly orthogonal to the direction of the trajectory when it passes the middle of the laminar part of the attractor CA at Qf = -51.0. The Poincare map VqF will be defined by this plane, and when investigating periodicity in the laminar part, we also need another condition

Z

1000

900

800

700

600

500

1400

for Pqf , namely, that trajectories cross the plane almost perpendicular. This is to avoid all intersections from other parts of CA. With this definition of the Poincare map, it is clearly not well defined for all QF values, because the behaviour of the laminar part changes with QF, but it is enough for investigating the attractor to approximately QF = -400.

To find the plane, new coordinates (x, g, h, f, y, z) are introduced such that system (10) has, in these coordinates, the stable equilibrium F2a at the origin for Qf = -50.7. We now have system (10) written in the form

x = K(x), x e R6 (11)

where k forms the right hand side of the equations and k(0) = 0. In the next step, system (11) is transformed into the eigenbasis to the Jacobian matrix for F2a at Qf = -50.7. This is done using diagonalization in the following way. Let e1,... ,e6 denote the set of eigenvectors to the Jacobian matrix, where e1 is the eigenvector corresponding to the eigenvalue which is zero at the fold bifurcation where F2a disappears. Then we put

P =[e1,...,ee] and x = [x, g,h, f,y, z]

T

Now, new coordinates u = [a, b, c, d, k, l] are defined by

x = Pu

System (11) can now be transformed into u coordinates by

PU = K(PU)

giving

u = P-1 K(PU) (12)

In system (12), variable a corresponds to the direction of the eigenvector e1, therefore, a will change quite fast during the laminar part of CA, where the rest of the variables, b, c, d, k, l will stay almost constant. The five-dimensional plane seeked for can now be found as a = 0.

3.4.2 Implementation and results

To see when a trajectory crosses the plane a = 0, we check if

alal-l < 0, i = 1, 2,3,... (13)

where i symbolizes the number of steps taken by the C-programme. In words, it is checked if the variable a has changed sign after every iteration. To constrain Pqf to intersections only from the laminar part, we define the norm N in the variables b,c,d,k,l by

N = b? + c? + d? + k? +i? Further, let the difference AN be given by

AN = N - Ni-1

We demand that ANi is bounded by a sufficiently small number e, that is

AN < e (14)

By testing, we found a proper value of e, giving only the intersections from the laminar part of CA. If (13) and (14) simultaneously hold after iteration i, the program checks if

ai <5 or ai-1 < 5 (15)

where 5 is a sufficiently small number determining the accuracy. If (15) does not hold, the programme reverses one iteration and takes half the step size and so on, until (15) is satisfied.

We show the result of 350 implementations for the parameter values QF = -51, -52,..., -400. By comparing the time t at different intersections with the plane, we found the proper value e = 0.001. Further, the accuracy parameter was 5 = 10-7. In each simulation, 2 • 107 loops were carried out, (corresponds to integrating to about t = 3 • 104) before starting to analyse the Poincare map to avoid transients. Then we let the simulations go on until 100 points were found

in the map. With this method it took about 17 minutes to analyse one parameter value, which means that the whole simulation needed about 100 hours. For QF = -51, -52,..., -296, we have a fixed point in the Poincare map, which means that the laminar part of CA occurs periodically with periodicity one. When the parameter decreases further, the attractor jumps between different periodicities and more complicated behaviour due to unknown bifurcations. An overview of this behaviour is illustrated in figure 9, which shows the results for Qf = -280, -281,..., -400. In figure 9(a), black lines are plotted for QF values corresponding to detected periodic solutions. In figure 9(b), solid lines indicate the location of periodic cycles with periodicity two, and the broken line indicates a periodic cycle having periodicity four.

«I

3

2

(a)

-400 -380 -360 -340 -320

Q_

3

(b)

2

-380 -360 -340 -320

Q

Figure 9: Strange bifurcation behaviour for QF = -280, -281,..., -400. (a) A black line is plotted for detected periodic attractors. (b) Solid lines indicate periodicity two, and the broken line indicates periodicity four.

When decreasing the parameter from QF = -51, the first more complicated result is obtained for QF = -306. The Poincare map for this case together with the result for QF = -305, a fixed point, are shown in figure 10. Here, solid squares are plotted for the case QF = -305, and dots for QF = -306. Figure 11 shows the Poincare map for QF = -324, a cycle with periodicity two, and for Qf = -359, a cycle with periodicity four. The Z-coordinate of the attractor CA for the above cases QF = -305, -306, -324 and QF = -359 is shown in figure 12.

To verify whether there are chaos in the more complicated attractors detected above or not, one can, in each case, determine the corresponding Lya-

punov spectrum. In Section 5, we present some theory about Lyapunov spectra and a method for determining the spectrum. Calculations have been carried out for QF = -300, -301,..., -400. The result indicates that there are strange attractors in the cases where no period was detected in the Poincare map. In table 5.1; Section 5.4, we show the resulting Lyapunov spectrum for the attractor CA at some parameter values between QF = -3000 and QF = -300, from which we conclude that there are chaos for QF values below approximately Qf = -600.

-0.1

-0.05 0

b

0.05

0.65 0.7 c

-1.56 -1.54 -1.52 -1.5 -1.48 -1.46 d

-0.45 -0.4 -0.35 -0.3 -0.25 e

Figure 10: Poincare map for QF = -305 as solid squares, and QF = -306 as dots. In each case 100 intersections are plotted. The first case shows a fixed point while the second case shows a more complicated behaviour.

0 0.05

b

0.75 0.8 0.85 0.9 0.95 c

-1.9 -1.8 -1.7 d

-0.4 -0.35

Figure 11: Poincare map for QF = -324 as dots, and QF = -359 as solid squares. In each case 100 intersections are plotted. The first case shows a two periodic cycle while the second case shows a four periodic one.

4 44 44 44 4 <

j[-irr IT W ^ IT tTT W W ^ w ^

(a)

(b)

(c)

(d)

Figure 12: Z-coordinate for the cases (a) QF = -305, showing a periodic cycle with periodicity one. (b) Qf = -306, more complicated attractor. (c) QF = -324, a two periodic cycle, (d) Qf = -359, a four periodic cycle.

e

Z

0

20

40

60

80

800

600

400

200

0

20

40

60

80

t

4 Bifurcations in the AB-plane

In this section we will investigate the fold and Hopf bifurcations2 that occur when varying the parameters A and B for QF = 0 and the remaining parameters to the values given in (2). In this case, the Bekryaev system (1) takes the form

d X Y f G\

— = -BUT X- BeG + 2.2GH + (Uv - H)— + [A + B + — ]Z ot 50 \ 50/

d G F 1

— = BeX -G- 2.2 XH + (H - Uu) — + —XZ ot 50 50

d F

= -UTG — F + (Uu — H)Y + GZ = UTX + (H - Uu)F - Y - XZ

dr

dY

dr

M^-GF + XY-^Z

Ot 5

One can see that system (16) is symmetric in such a way that the substitution

A = -A, B = -B, X = -X, G = -G, F = -F, Y = -Y

leaves it invariant. This means that we only need to consider variations in some values of the parameters A and B to understand the parameter plane.

4.1 The 20 regions

From Section 3 we know that system (16) has three equilibria for A = 2 and B = 0. The software package AUTO was used in order to follow these equilibria under variations in the free parameters A and B. By doing this, several bifurcations occurred that could be followed by the program in order to find the corresponding bifurcation curves in the plane of the parameters A and B. This analysis contains the fold and Hopf bifurcations, limit cycles are not analysed with exception of the fold-Hopf case studied in Section 4.2, however, some

2Theory about fold and Hopf bifurcations can be found in [9], [10] and [11].

stable limit cycles that occur will be illustrated in Section 4.3. Figure 14 and 15 show that the interesting part of the bifurcation diagram are fairly close to the origin in the AB-plane, therefore, an overview of all bifurcation curves and regions are shown in figure 13, which has an "unreal" scale to enable us to show all 20 regions bounded by the bifurcation curves in the same figure.

Figure 13: Overview in an unreal scale of the bifurcation diagram in the AB-plane. The Roman numerals symbolize the number of equilibria existing in each region. Solid curves indicate fold bifurcations and broken curves indicate Hopf bifurcations.

The bifurcation diagram will be described from figure 13, but we will also refer to the real bifurcation diagram illustrated in figure 14 - 17. One need to mention that nothing ensures that all existing equilibria and bifurcations are found in this analysis.

Far from the origin in the AB-plane the four regions Ia, Ib, Ic and 11 Ia dominate, and in region Ib, we have only one unstable equilibrium. At the Hopf bifurcation curves h2 and h3, this equilibrium becomes stable, therefore, we

have one stable equilibrium in the regions Ia and Ic. The curves marked /1, meeting at the cusp bifurcation point c1, denote a fold bifurcation. By crossing /1 from region Ia, two equilibria are created, one stable and one unstable. Thus, in region III'a and IIIa, we have three equilibria, two stable and one unstable. Moreover, close to the fold bifurcation curve /1 we have a Hopf bifurcation curve marked h1. These curves bound the narrow regions III'a, IIIh, V', IIIC, V and V'. The equilibria in these regions differ from their corresponding regions to the right with respect to the number of contracting directions on the unstable equilibrium created at the fold bifurcation curve /1. When crossing the Hopf bifurcation curve h1 to the left, the number of contracting directions increase from three to five, and simultaneously an unstable limit cycle is created.

From the fact that the stable equilibrium in region Ia loses stability when entering region Ih due to the Hopf bifurcation at the curve h2, we find that in region IIIh, we have one stable and two unstable equilibria. In region III2, bounded by the fold bifurcation curves marked /2, meeting at the cusp bifurcation c2, there are three unstable equilibria. By crossing the fold bifurcation curve /2 from IIIh one enters region V2. Here we have five equilibria, the unstable one from Ih together with the equilibria created at /1 and /2 respectively. Thus, there are one stable and four unstable equilibria in this region. In region IIIC, we have the same set of equilibria as in region IIIh, that is, one stable and two unstable ones.

For the parameter values A = 0.23114, B = -0.023035, the Hopf bifurcation curve h3 collides with the fold bifurcation curve /3 and disappears. At this point, we have two zero eigenvalues which means a Bogdanov-Takens (fold-fold) bifurcation (marked B-T in the figures). Figure 16 shows the location of this bifurcation. Moreover, the fold curves /1, joint at the cusp bifurcation c1 and the Hopf curve h1 intersect at A = 1.2426, B = -0.016664. Here we have one zero and a pair of pure imaginary eigenvalues which means a fold-Hopf bifurcation (marked F-H in the figures). Figure 17 shows the location of this bifurcation, which will be investigated further in Section 4.2. On the curve segments /3, meeting at cusp point c3, one stable and one unstable equilibrium arise if the curve is crossed above the Bogdanov-Takens bifurcation, otherwise, two unstable equilibria are created. Therefore, in region IIId, we have one stable and two unstable equilibria. Since the stable equilibria in Ic loses two contracting directions at the supercritical Hopf bifurcation represented by the

curve h3, there are three unstable equilibria in region III'd, coexisting with the stable limit cycle created at the Hopf bifurcation. This cycle is shown for some parameter values in region III'd in figure 24; Section 4.3.

In the regions Vd and Ve, we have the equilibria from region IIId together with the two equilibria created at the fold bifurcation marked /1. This means that there are two stable and three unstable equilibria in region Vd, and due to the Hopf bifurcation at the curve h1, there are one stable and four unstable equilibria in region Ve. Table 4.1 summarizes the number of stable equilibria in each region, also, it tells about the number of contracting directions of the equilibria. To understand the set of equilibria in region Vf, see Section 4.2.5.

From table 4.1 it is clear that all equilibria found in this analysis have one feature in common, they all have at least three contracting directions. This is also true for the equilibria discovered in Section 3, cf table 3.1; Section 3.2. Moreover, all limit cycles found also have this property, and from Section 5.4, it follows that all Lyapunov spectra determined have at least three strongly negative exponents. From these observations one may suspect that the dynamics in the Bekryaev system (1) take place on a three-dimensional invariant manifold.

A

Figure 14: Bifurcation diagram in the AB-plane.

A

Figure 15: Magnification of part of figure 14.

0.6 \ _ f2 1 1 1 1 1 1 1 / / f -'l >h1

0.4 - m2 /' _ / ./

0.2 V2 _

0 -0.2 f3

B-T ^ IIId f3

-0.4 I " c

-0.6

-0.8 1 \ 1 1 c3 1 1 1 1 1 1

0 0.5 1 1.5 2 2.5 3 3.5 4

A

Figure 16: Magnification showing the location of the Bogdanov-Takens bifurcation.

A

Figure 17: Magnification showing the location of the fold-Hopf bifurcation.

Table 4.1: The number of equilibria and their type in each of the 20 regions bounded by the fold and Hopf bifurcation curves.

Region No Stable Unstable Contracting directions

Ia 1 1 0 6

h 1 0 1 4

Ic 1 1 0 6

Uh 3 0 3 4, 4, 3

Ilia 3 2 1 6, 6, 3

Hl'a 3 2 1 6, 6, 5

Uh 3 1 2 6, 4, 3

mi 3 1 2 6, 5, 4

Ulc 3 1 2 6, 4, 3

Hl'c 3 1 2 6, 5, 4

Uh 3 1 2 6, 5, 4

Hl'd 3 0 3 5, 4, 4

Hie 3 0 3 5, 4, 4

V2 5 1 4 6, 4, 4, 3, 3

Vi 5 1 4 6, 5, 4, 4, 3

Vd 5 2 3 6, 6, 5, 4, 3

Vi 5 2 3 6, 6, 5, 5, 4

ve 5 1 4 6, 5, 4, 4, 3

VI 5 1 4 6, 5, 5, 4, 4

Vf 5 2 3 6, 6, 5, 5, 4

4.2 The fold-Hopf bifurcation

In this section we will start by presenting some theory from [10] about center manifolds that allow one to reduce the dimension of a given dynamical system near a local bifurcation. The three-dimensional dynamical system restricted to the center manifold for the fold-Hopf bifurcation at A = 1.2426, B = -0.016664 will be approximated with up to second order terms. Then this system is transformed into a normal form for fold-Hopf bifurcations, and this normal form is then analysed giving the local bifurcation diagram.

4.2.1 Center manifold theorems

Consider a continuous dynamical system defined by

x = f (x), x e Rn (17)

where f is sufficiently smooth and f (0) = 0. Let the eigenvalues of the Jacobian matrix A evaluated at the equilibrium point x0 = 0 be Ai, A2,..., An. Suppose that there are eigenvalues with zero real part, and assume that there are n+ eigenvalues (counting multiplicities) with Re A > 0 , n0 eigenvalues with Re A = 0 and n_ eigenvalues with Re A < 0. Let Tc denote the linear eigenspace corresponding to the union of the n0 eigenvalues on the imaginary axis and let ^ denote the flow3 associated with (17). Under the assumptions stated above, the following theorem holds.

Theorem 4.1 (Center Manifold Theorem) There is a locally defined smooth n0-dimensional invariant manifold Wc of (17) that is tangent to Tc at x = 0. Moreover, there is a neighborhood U of x0 = 0, such that if ^x e U for all t > 0 (t < 0), then ^x ^ Wc for t ^ to (t ^ -to).

Definition 4.1 The manifold Wc is called the center manifold.

3A definition of ¡p1 can be found in the beginning of Section 3.

In its eigenbasis4, system (17) can be written as

ü = Bu + g(u, v) (18)

v = Cv + h(u, v)

where u G Rn0, v G Rn++n-, B is an n0 x no matrix with all its eigenvalues on the imaginary axis, while C is an (n+ + n-) x (n+ + n-) matrix with no eigenvalue on the imaginary axis. Functions g and h have Taylor expansions starting with at least quadratic terms. Further, the center manifold Wc of system (18) can be locally represented as a graph of a smooth function

Wc = {(u, v) : v = V(u)}

where V : Rn0 ^ Rn++n-, and due to the tangent property of Wc, V(u) = O(yuy2). Figure 18 illustrates a two-dimensional center manifold in a three-dimensional dynamical system as the graph of a function v = V(u).

Figure 18: Center manifold as the graph of a function v = V(u).

4Recall that the eigenbasis is a basis formed by all (generalized) eigenvectors of A (or their linear combinations if the corresponding eigenvalues are complex). Actually, the basis used in the following may not be the true eigenbasis. Any basis in the noncritical eigenspace is allowed. In other words, the matrix C may not be in real canonical (Jordan) form.

The following theorem ends this section.

Theorem 4.2 (Reduction Principle) System (18) is locally topologically equivalent near the origin to the system

Notice that the equations for u and v are uncoupled in (19). This means that only the first equation, which is the restriction of (18) to its center manifold, has to be analysed to understand the dynamics locally near the bifurcation. In Section 4.2.2, we will show how to determine this equation for the fold-Hopf bifurcation.

4.2.2 Approximating the system restricted to Wc

The following calculations have been carried out using the software package MAPLE. As a first step, it is verified by the eigenvalues to the Jacobian matrix that n+ = 0, n0 = 3 and n_ = 3, which means that the center manifold to the fold-Hopf bifurcation is three-dimensional and stable.

We now define new parameters a = (a1 ,a2) and new coordinates (x,g, h,f, y, z) in system (16) such that the fold-Hopf bifurcation occurs at the parameter values a = 0 and at the origin in the new coordinates. These parameter values are called the critical parameter values and the point is referred to as the critical point. With these coordinates and parameters, system (16) takes the form

where k forms the right hand side of the equations and k(0) = 0. In the next step, diagonalization is used to transform system (20) into the eigenbasis of the Jacobian matrix at the critical point. This is done in the following way, let e1,... ,e6 denote the set of eigenvectors to the Jacobian matrix, where ei,e2 and e3 are the eigenvectors corresponding to the three eigenvalues on the imaginary axis. Then we put

U = Bu + g(u, V(u)) v = Cv

(19)

x = k(x, a), x e R6, a e R

2

(20)

P =[ei,...,ee], x = [x, g, h, f, y, z]

T

and define new coordinates w = [u\, u2, u3, v\, v2, v3]T by

x = Pw

The system X = k(x, a) can now be transformed into w coordinates by the transformation

PW = k(Pw, a)

which gives

w = P-1k(Pw, a) (21)

By writing equation (21) as two vector equations in u and v respectively, we find

u = a(a) + B(a)u + g(u, v, a), u e R3 v = b(a) + C(a)v + h(u, v, a), v e R

3 (22)

where a(0) = b(0) = 0. By setting a = 0 in (22) we obtain the 3 x 3 matrices B, C and the two functions g and h in the representation (18) of our system in the eigenspace to the Jacobian matrix at the critical point.

Our next task is to find the function v = V(u), where V : R3 ^ R3, locally representing the center manifold Wc. As mentioned before, Wc is tangent to the critical eigenspace Tc at the origin, and therefore, to obtain a second order approximation of the function V(u), we only need to include second order terms. We start by making the following approach, including all possible second order terms terms

vt = aU + bui + CiU2 + diUiU2 + etuiu3 + fiu2u3 i = 1, 2, 3 (23)

To determine the 18 coefficients abi,..., f\, i = 1, 2,3, we substituted (23) and the derivative of (23) into equation (18). Then equations for the coefficients could be verified and solved.

An approximation of the restriction of (18) to its center manifold given by the first equation in (19) is now determined, and we will follow the theory from Chapter 8 in [10] to determine the normal form in Section 4.2.3.

4.2.3 The normal form

In Chapter 8 in [10], the following lemma is derived for fold-Hopf bifurcations.

Lemma 4.1 Suppose that a three-dimensional system

x = f (x, a), x e R3, a e R2 (24)

with smooth f, has at a = 0 the equilibrium x = 0 with eigenvalues

Ai(0) = 0, Л2,з(0) = ± iwo, w > 0

Let

(ZP.1) g200 = 0, where g200 is defined by (29).

(ZP.2) g011 = 0, where g011 is defined by (29).

(ZP.3) E(0) = 0, where E(0) can be calculated using (31).

(ZP4) the map a ^ (7(a),^(a))T is regular at a = 0.

Then, by introducing a complex variable, making smooth and smoothly parameter dependent transformations, reparametrizing time (reversing it if E(0) < 0), and introducing new parameters, one can bring system (24) into the following form

С = A + С2 + s I С I2 +O(|| (C,C,C) II4) (25)

c = (ft + iw1)z + (0 + W)ic + С2Z + O(|| (i,c,c) II4)

where С e R, Z e C are new variables; в1, в2 are new parameters; 0 = 0(в), $ = $(в),w1 = w1(e) are smooth real-valued functions; w1(0) = 0; and

s = sign (^200^011) = ±1 (26)

m = Rehio (27)

^200

Only s and 0(0) are important in what follows. Assume that (ZP.5) 0(0) = 0

Four different cases of the fold-Hopf bifurcation can occur, established by the signs of the parameters s and 0(0). The two cases corresponding to s 0(0) < 0 are much more complex. To determine s and 0(0), we follow parts of the derivation of the normal form (25). During these calculations we also verify the conditions (ZP.1), (ZP.2), (ZP.3) and (ZP.5). Condition (ZP.4) is not verified here, this because there are some printing mistakes in [10], making it difficult to understand how y(a) is defined. Therefore, we suppose that (ZP.4) is satisfied, and apply Lemma 4.1 on the three-dimensional system restricted to the center manifold. The system is given by the first equation in (19), that is

u = Bu + g(u,V(u)), u e R3 (28)

where g(u,V(u)) = O(|| u ||2). Matrix B yields

B

0 0 0 0 0 -^0 0 ^0 0

where u0 = 46.529

and has the simple eigenvalues

Ax = 0, À2,3 = ± i^o

System (28) will now be transformed into a complex form. Let q0 e R3 and qi e C3 be the eigenvectors to B corresponding to the eigenvalues Ai = 0 and A = respectively, that is

Bqo = Aiqo, Bqi = Aqi

giving

qo =

1 0 0

and q\ =

0

i 1

Moreover, the adjoint eigenvectors p0 e R3 and p1 e C3 can be defined by

B Tpo = Aipo, BT pi = Api

which gives

Po =

and pi = —¿= v2

0

-1

These normalized eigenvectors has the property

(qo, po) = (qi, pi

The following orthogonality properties simultaneously hold due to the Fredholm Alternative Theorem

(qo, pi) = (qi, po) = 0 Now any real vector u can be represented as

u = ygo + Z qi + z qi

with

y = (u, Po > z = (u, Pi

In the coordinates y G R and z G C system (28) reads

y = g(y,z,z ) z = w0z + h(y, z, z )

where

g(y, z, z ) = (F(yqo + zqi + zqi), po> h(y, z, Z) = (F(yqo + zqi + Zqi),pi

are smooth functions of y, z, z whose Taylor expansions start with quadratic terms and are given by

5 As usual, (v, w) = v\w\ + v2w2 + v3w3 for two complex vectors v, w G C3.

g{y,z,z) = ~àû\gikiyJzk~zl ^

j+k+l>2 j ' ' '

and

1

jlklll

h(y,z,z) = Y, ^W**^ (3°)

j+k+l>2'

From expansions (29) and (30), the 14 coefficients needed to calculate s, 0(0) and E(0) are verified. Table 4.2 shows these coefficients from which we directly get the conditions (ZP.1) and (ZP.2).

Table 4.2: The 14 complex constants needed to calculate s, 0(0) and E(0).

Poll -0.00468 ^101 —0.0564 — 0.0192«

9no -0.0260 - 0.0207« ^110 0.0131 - 0.0972«

9in -0.0000461 ^002 0.000404 + 0.00253«

9020 0.0242 - 0.00494« ^020 -0.0000280 + 0.0209«

#200 -0.000108 ^200 0.197-0.131«

#300 -0.0000920 ^021 -0.00000679 - 0.0000876«

^011 -0.00221 + 0.00890« ^210 0.000220 - 0.000400«

E(0) is given by the expression

T?in\ Iß f TT I i f ReH02i g300 gm\ H021g20ö\ /Q1x

E{ 0) = -Re #210 + hno----1-------(31)

2 V V g0ii g200 g0ii j 2g0ii j

where G300, G m, H210 and H021 are calculated from

6

C300 = #300 --Ini(giio/i2oo)

^0

Gm = gm - — (2Im(#no/ioii) + Im(#020^i0i)) ^0

i 2 _

^210 = ^210 H--(^200(^020 — 2^110)— I hioi I —/1011^200)

^0

if i i ^

Hq2\ = ^021 H--^011^020 — 77^020^101 ~ 2 | /¿oil | — 77 | ^002 |

V 2 3

This gives E(0) = -0.0054650. Thus, condition (ZP.3) is satisfied. Finally, from (26) and (27) we establish that in the normal form given by (25) the parameters are s = 1 and 0(0) = -121.860, from which we also have the condition (ZP.5) satisfied. This means that we have one of the more complicated cases of fold-Hopf bifurcation.

4.2.4 Bifurcation diagram of the normal form

In coordinates (£,p, with Z = pe^, system (25) without O(|| • ||4) terms can be written as

£ = A + £2 + sp2

p = p(A + 0£ + £2) (32)

= +

the first two equations of which are independent of the third one. The equation for ^ describes a rotation around the £-axis with almost constant angular velocity ^ « for | £ | small. Thus, to understand the bifurcations in (32), one needs to study only the planar system for (£, p) with p > 0 given by

£ = ßi + £2 + sp2 p = p(ß2 + 0£ + £2)

2N (33)

The bifurcation diagram of (33) corresponding to the case s = 1, 0 < 0 is shown in figure 19. The system can have between zero and three equilibria in a small neighborhood of the origin for || A || small. Two equilibria with p = 0 exists for A < 0 and are given by

El = and E2 = (v^A, o)

These equilibria appear at the fold bifurcation on the line

S = {(A,&): A = 0}

The line S has two branches, S+ and S-, corresponding to A > 0 and A < 0 respectively. Crossing the branch S+ gives rise to an unstable node and a saddle, while passing through S- implies a stable node and a saddle. The node can bifurcate further, namely, a nontrivial equilibrium with p > 0,

ES = + o(/32), ]j--s{p1 + § + <H/%)))

that appears at the bifurcation curve

The nontrivial equilibrium E3 is a stable focus if A < 0, and an unstable focus if A > 0. For parameter values belonging to the line

T = {(A, A) : A = 0, A < 0}

the nontrivial equilibrium has a nondegenerate Hopf bifurcation and a unique unstable limit cycle exists for nearby parameter values. The cycle coexists with the two trivial equilibria Ei;2 which are saddles. Under parameter variation, the cycle can approach a heteroclinic cycle formed by the separatrices of the saddles, its period tends to infinity and the cycle disappears. The heteroclinic cycle appears along a curve orginating at ft = 0 and having the representation

p = {(A, A) : A = + A < o j

Now we can use the obtained bifurcation diagram for (33) to reconstruct bifurcations in the three-dimensional truncated normal form (32) by "suspension" of the rotation in ^ around the £-axis. The equilibria Ei;2 correspond to equilibrium points of (32). Therefore, the curve S is a fold bifurcation curve where two equilibria appear. Equilibrium E3 corresponds to a limit cycle in (32) of the same stability as E3. The curve H, at which a small cycle bifurcates from an equilibrium, corresponds to a Hopf bifurcation in (32). Moreover, the limit cycle corresponds to an invariant torus. Therefore, the Hopf bifurcation curve T describes a Neimark-Sacker bifurcation of the cycle, at which it gains stability and an unstable torus appears "around" it.

Figure 19: Bifurcation diagram for fold-Hopf bifurcation for the case s = 1, 0(0) < 0 in the ^-plane. In the phase portraits £ is on the horizontal axis and p on the vertical axis.

Finally, the curve P corresponds to a "sphere" in (32).

Adding higher order terms to the truncated system (32) will result in a nonequivalent bifurcation diagram. This because the spherelike surface that appears for parameter values on the curve P is an extremely degenerate structure, which disappears when adding higher order terms. Therefore, the torus cannot approach the "sphere", since it simply does not exist, and must therefore disappear before. Instead, system (32) may have near the curve P, in addition to local bifurcation curves, a bifurcation set corresponding to global bifurcations (heteroclinic tangencies, homoclinic orbits) and bifurcations of long-periodic limit cycles (folds and period-doubling cascades).

4.2.5 Comparing to simulated bifurcation diagram

By checking the stability of the two equilibria created when crossing the fold bifurcation curve / at each side of the fold-Hopf bifurcation we verified the branches S+ and S_. See figure 20, which shows the bifurcation diagram in the AB-plane close to the critical parameter values A = 1.2426, B = -0.016664.

Figure 20: Bifurcation diagram in the AB-plane close to the fold-Hopf bifurcation.

When moving from region IIId counter clockwise around the fold-Hopf bifurcation, we have the following. By crossing the curve /i entering region Ve, two unstable equilibria are created, which means that we have crossed branch S+. One of the equilibria has three (E2) and one has four (Ei) contracting directions. E2 gains two contracting directions when crossing the curve hi, which corresponds to branch H+, and an unstable limit cycle occurs which corresponds to the nontrivial equilibrium E3. By following this limit cycle in the negative A-direction, we found the Neimark-Sacker bifurcation. The corresponding bifurcation curve marked ti splits region Ve' into two regions, Ve'i and Ve'2. At the Neimark-Sacker bifurcation the cycle becomes stable and an unstable torus is created. Thus, the curve ti corresponds to the line T defined in Section 4.2.4. The stable limit cycle coexisting with the unstable torus in region Ve'2 is illustrated in figure 21. This cycle disappears by crossing the curve

c.

-0.018

1.2 1.21 1.22 1.23 1.24 1.25 1.26 1.27 1.28 1.29 1.3

A

hi into region Vf and simultaneously, equilibrium E becomes stable. From this we also conclude that in region Vf, we have two stable and three unstable equilibria. Finally, by crossing the fold bifurcation curve fi back into region IIId, the equilibria Ei and E2 destroys.

Figure 21: Limit cycle existing in region V' for A = 1.23, B = -0.0133.

4.3 Some attractors in the AB-plane

In this section, we show some attractors to system (16) that occur for different values of the parameters A and B. Figure 22 and 24 illustrate two limit cycles, one in region Ib and one in III'd, created at the supercritical Hopf bifurcations represented by the bifurcation curves h2 and h3 respectively. The following attractors are produced by choosing a point close to the origin as an initial point. Figure 23 shows a limit cycle in region Ib for A = B = 0, in which the variables H and Z are constant at the values H = 3.6397 and Z = 582.3547. Figure 25 and 26 illustrates limit cycles in region III2 for A = 2, B = 1, and in Ib for A = 0, B = 15 respectively. Recall that due to the symmetry in system (16), we obtain similar attractors by changing signs of A and B.

Figure 27 shows a more complicated periodic attractor that clearly reminds about the behaviour of the attractor CA, discussed in Section 3. In this case it is verified by calculating Lyapunov spectra that the attractor becomes chaotic by decreasing the parameter B, see Section 5 for how the spectra were determined.

Figure 28 illustrates a limit cycle for A = 15, B = 18. By increasing both parameters to A = 20, B = 25, the attractor shows a quasi-periodic behaviour, and the attracting set is an invariant torus, see figure 29. This was verified by calculating the corresponding Lyapunov spectrum, having two zero exponents.

-150 -100 -50 0 50 100

600 550 M 500 450 400

420 Y

440

>- 420

400

X

F

600

550

N 500

450

400

400

440

H

Figure 23: Limit cycle in region Ib for A = B = 0, the variables H and Z are fixed to the values

H = 3.6397 and Z = 582.3547.

Figure 24: Limit cycle in region III'd close to the Hopf bifurcation curve h3 for A = 2, B = -0.3536.

Figure 25: Limit cycle in region III2 for A = 2, B = 1.

Figure 26: Limit cycle in region Ib for A = 0, B = 15.

Figure 27: More complicated attractor for A =15, B = -5.

5 The Lyapunov spectrum

In this section we begin by presenting a definition and some properties of the Lyapunov spectrum. Then a method for computing the exponents is shown, and finally, we describe how to implement this method on a computer and show some results. The theory presented here has been found in [12] and [13], and we refer to these papers for more details.

5.1 The Lyapunov spectrum defined

The Lyapunov spectrum is a striking characterization of an n-dimensional dynamical system. It associates a set of n real values to each orbit of the system which describes exponential instabilities of infinitesimal deviations from the orbit. Moreover, for an ergodic6 dynamical system, the spectrum is independent of which orbit you choose. In more detail, consider a continuous n-dimensional dynamical system given by

x = v(x), v,x e Rn

For an initial condition x(0) = x0, we integrate the system to obtain a corresponding orbit

x(t) = ^txo

To examine the stability of this orbit, we look at the evolution of a nearby orbit x(t) + u(t). In the first step, we linearize the equations of motion in u to obtain

u = J (x(t)) u(t) (34)

where J(x) is the Jacobian matrix at the point x. By integrating (34) along the orbit we obtain the tangent map

u(t) = MX0 (t)uo where the time dependent n x n matrix Mx0 (t) is given by

6Theory about ergodic dynamical systems can be found in [14].

The exponential instabilities of a trajectory are now reflected in the eigenvalue spectrum to the matrix Mx0 (t), or rather, since the Lyapunov exponents are related to the modulus of the eigenvalues, the spectrum of the symmetric product

M X0 (t)Mxo (t)

The eigenvalues of this matrix are real and positive and we order them in the following way

M?(t) >M2(t) > ... >A(t) >0

From the analysis above one realizes that these eigenvalues are dependent of the initial point x0 chosen. To overcome this we have a theorem by Oseledec giving that the limit

Xk = lim \ log fik(t), k = 1,...,n (35)

t

is independent of, and exists for almost every initial point x0. Thus, by taking an arbitrary initial point and calculating the above limits, with probability 1 you will get its unique Lyapunov spectrum, {Ai > A2 > ... > An}.

5.2 Properties of the Lyapunov spectrum

From the definition of the spectrum it can be proved that it is independent of the choice of coordinate system. Moreover, any continuous time dependent dynamical system without a fixed point will have at least one zero exponent, corresponding to the linear changing magnitude of a principal axis tangent to the flow. Axes that are on the average expanding correspond to positive exponents, and contracting axes correspond to negative exponents. Also, for a dissipative dynamical system there will be at least one negative exponent. Further, an attractor for a dissipative system with one or more positive exponents is said to be a strange, or a chaotic attractor. From this discussion, we conclude that, for example in a three-dimensional continuous dissipative dynamical system, we can have the following spectra; (+, 0, —), a strange attractor; (0,0, —),

a torus (quasi-periodic behaviour); (0, —, —), a limit cycle; and (—, —, —), a fixed point.

5.3 Determining the Lyapunov spectrum

From the numerical point of view, the above description is insufficient because

rjl

the matrix Mxo Mx0 becomes singular rather fast since its eigenvalues separate exponentially in time, if not all exponents are equal. This makes it difficult to measure the spectrum, and in this section we shall present a method for doing this in which we augment the given dynamical system with an orthonormal frame and a Lyapunov vector. The method applies to any finite-dimensional dynamical system.

We define a time-dependent orthonormal k-frame to be the set of k(k < n) orthonormal n-dimensional vectors

s(t) = ei(t),..., efe(t) (e*, ej) = bl3 i < i,j < k (36)

where (.,.) is the usual Euclidean product in Rn. Using this frame, we let the matrix elements of the Jacobian matrix J be given by

Jim = (ei, Jem), 1,m < k

which depend on time both through the Jacobian and the frame. Further, we introduce a stability parameter ^ > 0, and define the symmetric stabilized matrix elements Lmm and Lim as

Lmm — Jmm + @ ((em, em) 1) , m < k

and

Lim = Jim + Jmi + 2^(ef, em), I = m, 1,m < k

Finally, let A = (A1(t),..., Ak(t)} be a k-dimensional real vector. The augmented dynamical system is now given by the following set of differential equations, of which the first two are vector equations

x = v(x), v,x e Rn

em = Jem — ^ ei Lim m = 1,...,k (37)

l<m

Am Jmm m' 1, . . . , k

For system (37) we have the following theorem.

Theorem 1 Let x0 be an initial point for which the associated Lyapunov spectrum (cf equation 35) Ai > A2 > ... > An exists. Set A(t = 0) = 0. Choosing the stability parameter ^ > — Ak, then for almost any (i.e. with probability 1 when choosing randomly) initial frame s(t = 0) the time evolution of the dynamical system (37) yields

lim -Am{t) = Am m = 1,..., k

t—t

Thus, by following a trajectory of the augmented system we obtain almost surely the k first exponents in the Lyapunov spectrum for the given orbit. The condition on the stability parameter is satisfied, for example by setting

maxi|e||=i(—(e, Je)) (38)

where the maximum is over all unit-length vectors e and over the relevant region of phase space. Dynamically, this corresponds to finding the strongest local contraction. For details and a proof of the above theorem, see [12].

5.4 Implementations and results

In order to implement the augmented dynamical system (37), we use the same algorithm as in the rest of the simulations on the Bekryaev system, that is, a Runge-Kutta method in C-code. When calculating the complete spectrum, we get k = n = 6, which means that the augmented system to implement gets a dimension of 48. As an initial frame, we chose the orthonormal set of vectors

ei (j 1 j 1 < < 6

^ 0 i = j " "

Also, we start at a point on the attractor of which we want to calculate the spectrum. Further, in addition to the main system to simulate, the program needs the Jacobian matrix explicitly to be able to set up the equations for x, em and Am in system (37). To choose the stability parameter 3 such that condition (38) is satisfied, we started with 3 = 0, and at each step in the simulation we check the condition

3 < max(-Jmm) m = l,...,k (39)

If (39) holds we let 3 = max(-Jmm) + l. After completed one simulation one may fix 3 as some value greater than the value obtained above to speed up future simulations.

Figure 30 shows the result of a simulation of the Bekryaev system (l) with the parameters A = 2, B = 0, QF = -2500 and the remaining parameters as in (2). We have the attractor CA discussed in Section 3, the spectrum obtained yields (+, 0, —, —, —, —) which indicates a strange attractor. In this case we had the stability parameter 3 = 416. The elapsed time for this simulation was approximately 20 hours. In table 5.1, we present the result of some more calculations of the Lyapunov spectrum for the attractor CA. All spectra obtained for the Bekryaev system l include at least three strongly negative exponents.

Table 5.1: Lyapunov spectra for the attractor CA at different values of QF.

Qf Lyapunov spectra Qf Lyapunov spectra

-300 (o,-,-,-,-,-) -600

-325 -650

-350 -700

-375 -800

-400 -900

-450 -1000

-500 -2500

-550 -3000

3.5

S 3

0

500

-0.5

CO

* -1 -1.5

^v——

-1.5

^ -2

0

500

-2.5

-3

<cz

500

l5 -3.5

500

500

1000

1000

1000

1000

1000

1500

2000

1500

2000

1500

2000

1500

2000

1500

2000

2500

2500

2500

2500

2500

3000

3000

3000

3000

3000

3500

3500

3500

3500

3500

3500

Figure 30: The six Lyapunov exponents from a single run of the Bekryaev system (1) for the parameters A = 2, B = 0 and Qf = —2500.

References

[1] R.V. Bekryaev, Violation of exponential divergence of trajectories in a system of hydrodynamic type with a chaos, 1994, St. Petersburg: Allerton Press, Inc. BRAS Physics / Supplement, Physics of Vibrations: Vol. 58, No.1: 49-56.

[2] James R. Holton, An Introduction to Dynamic Meteorology, 1992, San Diego, Academic Press, ISBN 0-12-354355-X.

[3] Edward N. Lorenz, Deterministic Nonperiodic Flow, 1963, Journal of the Atmospheric sciences, Vol 20: 130-141.

[4] Edward N. Lorenz, Irregularity: A fundamental property of the atmosphere, 1984, Tellus 36A:98-110.

[5] Edward N. Lorenz, Can chaos and intransitivity lead to interannual variability?, 1990, Tellus 42A: 378-389.

0

0

0

t

[6] H. Broer, C. Simo and R. Vitolo, Bifurcations and strange attractors in the Lorenz-84 climate model with seasonal forcing, 2002, Nonlinearity, Vol 15: 1205-1267.

[7] L. Van Veen, Baroclinic flow and the Lorenz-84 model, 2003, International Journal of Bifurcation and Chaos, Vol 13, No. 8: 2117-2139.

[8] http://indy.cs.concordia.ca/auto/main.html

[9] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical systems, and Bifurcations of Vector Fields, 1983, New York: Springer-Verlag. ISBN 0-387-90819-6.

[10] Yuri A. Kuznetsov, Elements of Applied Bifurcation Theory, 1995, New York: Springer-Verlag. ISBN 0-387-94418-4.

[11] Stephen Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, 1990, New York: Springer-Verlag. ISBN 0-387-97003-7.

[12] F. Christiansen and H. H. Rugh, Computing Lyapunov spectra with continuous Gram-Schmidt orthonormalization, 1997, Germany: Department of Mathematics, University of Warwick, Coverty CV4 7AL, UK.

[13] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Determining lyapunov exponents from a time series, 1984, USA: Department of physics, University of Texas, Austin, Texas 78712.

[14] A. Katok, B. Hasselblatt, Introduction to Modern Theory of Dynamical Systems, 1995, Cambridge: Cambridge University Press. ISBN 0-52134187-6.