Strange attractors in coupled reaction-diffusion cells

Physica 5D (IQS2) 258-272 North-Holland Publishing Company

STRANGE ATTRACTORS IN COUPLED REACTION-DIFFUSION CELLS Igor S C H R E I B E R and Milog MAREK 13W,parlmenr of Chemical En,~ineering. Prague Institute oF Chemical Technology. 166 26 Prague 6. Czechoslovakia Received 15 July 1981 Revised 29 January 1982

A numerical study of two identical reaction cells with diffusion coupling has shown that the structure of motion in the system in pr/nciple agrees with results reported for variety of dynamic syxtems. When the characteristic parameter is varied, alternating sub-intervals of stable periodic (P~) and stable aperiodic (Aj) solutions appear. The sub-intervals are ,:onnected by inlervals, where tangent bifurcations and infinite sequences of subharmonic bifurcations occur. Feigenbavm relation holds for the studied sequence of subharmonic bifurcations. Aperiodic (chaotic) states are characterized :',y a complete set of one-dimensional Lyapunov exponents, by power spectra, and by corresponding Poincar6 maps. The ,,pectra of Lyapunov exponents are of the type ( + , 0 , - , - ) , and show that the topological dimension of the chaotic attractor is lwo. The power spectra are of two different types (a) the spectra with sharp peaks ~ocated above the broad-band noise, showing statistical phase coherence of the attractor, (b) the flat spectra showing only broad-band noise, corresponding to phase incoherent attractor. The phase coherence is present close after every point of accumulation. Phase incoherence arises when the strange attractors contain unstable periodic orbits with different topology. The relations of bifurcated stable and unstable periodic solutions (computed by means of continuation techniques) to the structure of the strange attractor is discussed and Poincar~ maps are used to illustrate tie dependence of the structure of the attractors on the value of the characteristic parameter.

1. I[ntroduction

~ oupled cell,: with reaction and mutual mass exchange are often used as model systems for description of processes in living cells, tissues, in compartmental models and various other dis'.ributed systems with chemical reactions and transport (including various forms of chemical reactors) [I]. A system of coupled well-stirred reaction cells with mutual mass exchange through common walls is used for experimentation in our research group for a number of years. The time course of concentrations of reaction components is monitored and various steady states and time dependent regimes are observed [2]. When Zhabotinskii reaction (oxidation ~.~f malonic acid by acidic bromate with ceric/cerous ions catalyst) is used in the coupled cells, coexisting two steady states in two cells [3a] and various combinations of steady states (Turing stnJ,ctures) in the system of up to seven coupled ceils can be observed [3b]. The obser-

0167-27~9/82/000~0000/$02.75 © 1982 North-Holland

ved time-dependent regimes are either periodic (various synchronized regimes) or aperiodic [4]. In these experiments the course of oscillations in the cells have been observed in dependeace on the varying intensity of mutual mass exchange between the cells. The mass exchange can be described as a linear diffusion coupling. However, the models of detailed reaction mechanism of the experimentally studied Zhabotinskii reaction system are complex and resulting systems of differential equations are of a high order and stiff. Hence simple models are studied first. We discuss the results of analysis of dynamic behaviour of two identical coupled reaction cells with the Brusselator model. Recently there appeared a large number of papers on the subject of strange attractors. In the case of coupled oscillators turbulence can be due to coupling, i.e. a set of decoupled oscillators does not reveal chaotic behaviour. The most simple coupled system exhibiting chaotic behaviour is nonlinear oscillator subjected to an

Y. 8chreibcr and M. MarekI Strang¢ attractors in coupled cells

external periodic force [5,6, 7]. Chaotic states may also occur, w h e n two or more oscillators are coupled together [8,9, 10]. In m ~st of the above systems, t infmite ~sequence cations, However, turbulence are still p0ssible. F o r example metastable chaos [11] and intermRtent chaos [12] are known from thp ~Lorenz mc~el [13], stochastization of quasiperiodic motion was described in the studies of conservative systems etc. 2. Strange attractors in coupled Brusselatorsnumerical results In this numerical study we demonstrate that bifurcation phenomena leading to a strange attractor and chaotic states in two identical coupled nonlinear oscillators with linear diffusion coupling are very similar to the bifurcation phenomena which are well known from the g~ublished studies of discrete and continuous dynamic systems e,g, refs. 14-1"]. Coexistence of attractors for given value of the parameter is al~:o illustrated. Three methods are used for characterization of periodic and chaotic states. Power spectra a-e computed (fast Fourier transform is used [18]), all one-dimensional Lyapunov characteristic exponents are determined [19, 20] and Poincar6 maps are studied. A strange attractor is characterized by power spectra with high noise level, po,*itive value of at least one of the characteristic exponents or by a Poincar6 map with a Cantor set-like structure. 2.1. Model

Two Brusselators* with linear diffusion coupling are described by the following system of ordinary differential equations: Jtt = A - ( B + l)xt + xl2yt+ DI(x2- x0, * Tyson 110] was the first to observe complicated motions in a modified Brusselator model.

259

~I ffi B x l - x~Yl

+ D:(Y2- Yl), ~2 ffi A - (8 + l)x2 + x~y: + DI(x~- xg, ~ 2 - 8x~-x~y:

+ D:(y,- y:).

(1)

It is well known [21] that one Brusselator can have one stationary state (x = A, y = B / A ) , which is stable for B < 1 + A 2, and one periodic solution, bifurcating from the stationary state at B = 1 + A 2. The periodic solution is stable for B > 1 + A 2. Both solutions correspond to homogeneous solutions of the system (1). The homogeneous solutions satisfy the condition xt -- x2, y t - Y2. It means that all homogeneous solutions of the system (1) are located in the (two-dimensionaD diagonal plane xt = x2, yt ---- Y2 of the R' phase space (xt, Yt, X2, Yz)- However, the stability o f homogeneous stationary H $ and periodic'. P , solutions of (1) can be different than is the stability of comparable solutions of the single Brusselator. Nonhomogeneous solutions (xt ~ x2, y~ ~ Y0 of the system (1) can also exist. The values of parameters in our numerical study are fixed at A = 2, B - 5.9, Dt/D: -- q - O. !. We follow behaviour of the system in dependence on the value of D~ or equivalently D2 (bifurcation parameters). Structure and stability of stationary solutions, existence of subharmonic bifurcations and existence of a strange attractor have been discussed in [4b]. Here we discuss in detail the structure of solutions of (1) for the values of the bifurcation parameter DtE(0.9, 1.5). Global situation is sketched in fig. 1. Unstable homogeneous stationary solutions HS (xtfx2fA, Yt=Y2=BIA) exist for all values of DI. Other two existing stationary solutions NS, N S are nonhomogeneous. 'Fhey ...,.,.,,, mlrror images with respect t(, the s vs diagonal plane, i.e. N S =(xl, y~, x2, . ~) and N"S = (x], y~, x~, yl). Their stability is the same and both are stable for DI <0.a542. At D! = 0.9542 two symmetric pe;~odic solutions P, P bifurcate; the solutions are stablte on the interval DtE(0.9542, 1.1720). These period~,c solutions are also symmetric with respect to the diagonal

260

I. Schreiber and M. MarektStrange atlractors in coupled cells

/

~

/,,J

$~,e ~ o e o e O ° v ~ ° ~ ° ° ° ° ° ° ° ° ° O O o o ~ _

~"

,.,

oo-' . . . . . . . . .

' i,"

,~x-eooec'oo.~c o o o o o ~ e e ~ o ~ o o e O0 0

0

Oo

.-'/'

!

:;~ 0

o0~oo~O0@@'OO0

o00% O0000

C~0 0 0 0 0

O0'z

"'/~

!

"

. . . . . . . . . .

e o o~ e eee ~eee •

~i

,

~

___z-"_

oo~oo--__

_,, _ , , ' : ~ . , ~ , ~ ~ o o'

"" 'i- A

~ , ,,_,..~,a_................. -r ........................

Fig. I. S c h e m a t i c picture of s t e a d y s t a t e s a n d basic p e r i o d i c ~olufion,, of eqs. (i) for D ~ , ( 0 . 9 , 1.5). A = 2, B = 5.9. D,/I): : ,! :- 0.1: stable s t e a d y state; unstable s t e a d y state: O ~ O s t a b l e periodic orbit; ~? ." ~ ,? u n s t a b l e periodic orbit.

plane: the phase shift between both oscillators is negligible. We call such solutions "in phase" solutions. We study just o . e of them but in the following discussion we have to keep in mind that nonhomogeneous periodic (and any other solutions of the system (I)) occur in pairs symmetric with respect to the diagonal plane ~v~ith the exception of anti-phase periodic solutions, which are self-symmetric). After the loss of stability both P and /5 pass through the limit points and finally both disappear at the point (D~ = 1.17.19) where the homogeneous periodic solution P . changes its stabil-

ity (see section 2.4 for more details). The homogeneous periodic solution PH is unstable (with respect to, inhomogeneous perturbations) for D~ < 1.1219 and stable for D, > 1.1219. In cb ~.,tic region, Dt E (1.1933, 1.2o35), other perio~tic and aperiodic solutic~ exist. These are shown in more detail in fig. 2. Because nonhomogeneous attractors (periodic or strange) exist for D~ > 1.1219, multiplicity of attractors and hysteresis in the dependence of solution on the parameter Dt occurs. Hence to approach various attractors different initial conditions have to be used. For example: (a) Initial condition (2.1, 2.9, 0, 0) leads to PH. (b) Initial condition (2.1, 2.9, 1.9, 3.0) leads to a nonhomogeneous attractor NA. (c) Initial condition (I.9, 3.0, 2.1, 2.9) leads to the symmetric NA. The above initial conditions are used for rill values ,,f D,. Usual numerical integration procedure (Runge-Kutta-Merson routine), speciall:¢ constructed numerical algorithm for continuation of periodic solutions in dependence on the parameter D~, D: and construction of Poincar,~ map are used to cha, acterize existing attractors. In most cases !nitial conditions (b) were used for integration: after transients have died away (i.e. usu'.ily after hundred or more oscillations), either Poincar6 maps o: chosen projections of

A3

i ~.l,q

_..~___.

1.18

3 b_ b !

W t

!.22

P2

.

i

CO~0

.

.

1.20

.

II

:

'

p.

C¢,~ 3

''

1.22

3 !

1.24

C;

1.26

D1

~-'i~, 2. Detailed piclure of the region /)= E~ ( I. 16. 1.27) t'~ - region o f a sta! ',c p e r i o d i c .~olmion w i l h k fixed points in the Poincar6 m;,r, A ~. A : . . . . . . 4, - region~ of a p e r i o d i c behaviour

L Schreiber and M. MareklStrange atlractors in coup ed cells

trajectories were plotted. The results are schematically presented in fig. 2. In the interesting range of the parameter D~ E (!. 16, 1.27), periodic Pk and strange attractors SA alternate. Periodic omits Pk can be differentiated by a number of fixed points in a properly chosen Poincar6 map; Pk will denote periodic orbit with k fixed points (for PH is k = 1). The region A where the chaotic motion can occur (D! E(1.1933, 1.2635)) is divided by the windows of periodic attractors. The number of periodic windows is high, but only some of them are relatively larger. Considering the larger periodic windows, we have divided the region A into several parts A~ (j - I, 2, 3, 4, 5). Every interval Aj still contains a large number of periodic windows, however, tiley are very narrow and thus hardly detectable in observations of real physical systems. Sequences of bifurcations of new periodic solutions with double period appea." very often in the windows of periodic solutions. The rate of convergence of these sequences to the limit values D , - C~ can be cotapared with the universal value suggested b,., Feigenbaum for large class of maps on a unit interval [22]. Such a comparison for the first sequence from the left in fig. 2 is discussed later. The basic informations on the attractors shown in fig. 2 are listed in table I. 2.2. Power spectra Fast Fourier transfcrmation [18] of the discrete time series obtained by numerical integration of the system (I) was used to obtain power spectra. Every time series consisted of 214 samples of the time dependence of one of the coordinates xd, y, (it approximately corresponds to 400 oscillations): the samples were taken after transients have died away. Let us discuss the properties of t0ower spectra in the region of DIE(I.16, 1.27). Every subharmonic bifurcation P ~ : , ~ P k , , . , corresponds to the appearence of sharp peak at the half of t h original frequency. The difference between the mean level of the noise and the highest

261

peaks of the power spectra is 9-10 orders of magnitude. When the first accumulation point is reached, the situation changes. The sharp peaks are still present, but the p_oise level increases by 4 - 6 orders. It indicates the presence of a strange attractor: If we move away from the accumulation point, a reverse b~furcation sequence occurs [23]. Periodic components with lowest frequency subsequently disappear from the power spectrum. The final state corresponds to initial state of the forward sequence of bifurcations of periodic orbits, i.e. one periodic component (together with its harmonics) is present in the spectrum above the broad-band noise. The frequency of the periodic component corresponds to the frequency of the basic periodic solution Pk but it is a little shifted. Examples of computed power spectra are given in figs. 3a-g. The spectra 3a. b correspond to the reverse bifurcation sequence appearing after the first point of accumulation C~ (i.e. it is located in the first chaotic region A I in fig. 2). The periodic component with the lowest frequency appearing in the spectrum 3a corresponds to the second subharmonic (i.e. to the P4 orbit in the forward bifurcation sequence). The spectrum 3b arises after two reverse bifurcations from 3a. The spectra 3d and 3f contain basic periodic component corresponding to a P~ orbit or to a P.~ orbit, respectively. The spectrum 3d (DI = 1.2250) corresponds to a final state reached after the reverse bifurcation sequence to the right of the periodic window PW~t. has taken place. This spectrum corresponds to the trajectory starting from the initial conditions (b) or (c). However, if we choose other initial conditions (for example an initial condition located on the periodic orbit P, in the window PW~) we obtain power spectrum corresponding t¢, the periodic solution PI2. This solution belongs to R the sequence of periodic orbits in PW~. Five different attractors thus exist for D! = !.2250. One symmetric pair of strange attractors SA, S-"A, one pair of periodic orbits P~2, P~2 and one homogeneous periodic solution Pa. In most cases the observed reverse bifur-

_¢,_

I. Schreiber and M. Marekl $1range attractors in coupled ceils TABLE I List of attractors in fig. 2 Range of Dt

Attractor. remark~

i.172-1.189 1,190-LI92 1.1925- I. 1930 1.19311.1933- 1.2050

P: attractor P+ attractor Ps attractor P,~ attractor chaotic region A,: contains many narrow periodic window containing, for example, attractors PI:, Pro. P.,, Ps first larger window of periodic solutions PW~, contains progression P;, Ps, Pt:,

1.205 !- !.2068

P24.

I ...'~127- i .2157

!~2~8 !.22M

i .~~ h8

- -

1.2240

!....50-1 n~

1.2250

1.2420-1.2440 1.2635

• •

second larger window of periodic solutions PW, contains only Ps (r:o subharmonic solutions) two windows of periodic solutions very close together a~e ~:ontained in this region left window of periodic solutions PW L, contains sequence (from the left to the right) P~, P~, P t : . . . right window of periodic solutions PW~: the sequence P.~. P~, Pu is now from the right to the left: attractors in the space between Pg,'[ and PW~ are mostly chaotic, however, for Dt = 1.2245 attractor Pts exists for this value of Dt exist simultaneously P,, (from PW$) and a st, ange attractor SA, which probably corresponds to a final state of reverse bifurcation sequence arising to the right from PW~+; several narrow windows with P~ orbits exists to the left from PW~ (i.e. in A0 and to the right from PW$ (i.e. in A+) last larger window PW: of periodic solutions with attractors P:, P,, P, . . . . end poin~t F of the region of existence of SA; windows of periodic solutions containing e.g. basic attractors P~, P~, P:, P4, P~ still exist above this value of D~. Each of these windows contains a sequence of bifurcated periodic, solutions with double period

Attracting periodic homogeneous solution Pu exists for all values of O1> 1 . 1 2 1 9 . . . C,tti~,tl " ' : ' - 5~qu~iqces

+_t._ tar, c place on a far smaller

int,+rval of Dt than the corresponding forward bifu:cation sequences. The ,~orward and reverse bifurcation sequences can mutually superimpose; this probably occurs in the reighbourhood of the "'third" window (i.e. the windows PW3L and PW~) because we can observe power spectra

with periodic components also outside of the above "window". The presence of a periodic component in the spectrum means that the attractor is statistically phase coherent [23, 24], i.e. a :+et of originally close points on the attractor will expand mainly in the transversal direction to ,'he flow (phase

L Schreiber and M. Marek/Stran~eattractors in coupled cells

I°g PI

I

o

10g F' l

I

o

4 ........................................

P

log

b

l

.....

l ........

l

1"-" . . . . . . .

o.,, b4

.......

log P

d

-4

'

!

4-! .................................................................. tog P i

e

/

- 4 4

I

.

.

.

.

.

-i

.

.

.

.

.

.

....

.

.

I .................. ~

.

.

.

.

.

.

t

f i

I

0~

-4 ................... l ...................

0

r ............

o.2

Fig. 3. Power spectra of various strange attr actor s - log P -f [Hz]. A = 2 , B - S . 9 , DI/D_~ -- 0.1. (a) D l = 1.1940- phase

4 7 ................................................................... t o g P~ c

t°g Pt

~

coherent SA; (b) D~= 1.20-phase coherent SA: (c) D,-1.21 - phase incoherent SA; (d) D~= 1.225- phase coherent SA; (e) D~ = 1.235- phase incoherent SA; (f) Dt--!.245pha,,e coherent 5A; (g) Dt = 1.263-phase incoherent $A.

-4.1

4 7.

263

0.2

1 .................

0.4

1 .......................

0.6

f

,-~.~ 0.8

L~.I

shift between individual points is changed only little by the flow). Hence every strange attractor arising in the reverse bifurcation sequence is phase coherent. The power spectrum 3f illustrates the fact that phase coherent attractors are not limited to the neighbourhood of larger windows of per,iodic solutions. This power spectrum corresponds to ~he phase coherent strange attractor which is related to a p e : : . l i c orbit P~. A number of power spectra simf'~r to the spectrum shown in fig. 3f can be foullY. Their existence indicates the presence of a window of periodic solutions in a close neighbourhood. There is a possibility to connect the phase coherence with the existence of the Net of unstable period;~ solutions (and their unstable invariant manifolds) which are the continuations of solutions arising in forward seq~lence of subharmonic bifurcations. This was suggested in [15] for a diffeomorphism and may be relevant also in the system (I). We now follow the changes in the power spectra after the reverse bifurcation sequence behind the points C~, C~, C~ has occurred. When DI is increased a little, the sharp peaks in the spectrum either completely disappear or become far broader and lower. The power spectra are then similar to those shown in figs. 3c, e, g. These P~t spectra cover most of the chaotic regions A2, A~, A4 and A~. In the chaotic region A~ such a flat I:ower spectra occur for

2(',-1

!. Schrt,',:,er and M'. Marekt S t n m g e attractors in cottpled cells

D, "-1.201. The stran~,e attractor corresponding to ,he flat power spe,:trum is phase incoherent [2-11. The flow on the attractor is mixing both in the transversal and longitudinal directions and thus the phase shift between the laoints which were originally cl<~se together varies considerabl~~. The tin'e course ef oscillations ef .x,, y~, i = 1,2 for [), whi,'h c~rresponds to a phase inct;heren~ attractor can exhibit irregular alterna',ion of oscillations similar to various periodic orbits. For example, in fig. 4 we can clearly differentmte two-pe~ak and three-peak oscillations. Thi~ supports the idea that the pha~.e incoherence depends on the existence of unstable periodic trajectories with qualitatively different periods Fr (i.e. periods which are not approximz~,te powers of 2. for example T~,, and T,,), o~ generally on the differences in the topology of unstable periodic orbits. This is discuqsed in detail later.

O:

t

1

,

................. ~...............

X~

oi o

"W. . . . . . . .

60

t-'-

80

..............

--t

100

Fig. 4. T i m e c o u r s e o f o s c i l l a t i o n s of t h e p h a s e v a r i a b l e s in t h e first r e a c t i o n cell. x~(t), y~(t). I r r e g u l a r a l t e r n a t i o n o f t w o and three peak cycles. A=2, B = 5 . 9 , DdDt =0.1. DI = 1.21. Initial c o n d i t i o n s {xt. Yr, x,., y.,) = (2.1, 2.9, 1.9, 3.0).

T A B I . E It ( ) n e - d i m e n s i o n a l L y a p u n o v e x p o n e n t ~ Xk a n d t h e H a u s d o r f f d i m e n s i o n D of s o m e s t r a n g e attrac':ors o f t h e s y s t e m (I)

2.3. Lyapunc,r c h a r a c t e r i s t i c e x p o n e n t s ",.Ue have computed 1 2 5 ] k-di:mensional l_y:~tpunov e:,laonents k (k = 1, 2, 3, 4 in our ,::tse} [26] b~, tbe method similar to tha! proposed in ll9! or in [210]. Each number k k exp~esses a,,,~rage exponential rate of increase of :he volume ef k-dimensional parailelepiped by the action of the flow. The numerical method used ,,~ives m:~.ximal values k,,~,k with probability one ( k can attain at most (~) different values). All one-dimensional exponents kk (k = I, 2, 3, 4) can be computed as k:=k I

fl~, , t $. '

k~=k t

I'll ; I ~

-k ~ I

Illll ~ •

k = ~~

'l

3,4.

()no-dimensional Lyapunov exponents kk describe the rate of divergence of the neighbt, aring Tra.iectories in four independent directions. We have chosen for illustration several examples from chaotic regions A~, A , and ~,~. The results are given in table 11. The values of kk are accurate., within 5'% rel. A better convergem-e of

!.20 1.21 I "•_ 6

I).(178 0 107 0.200

0.IX) 0.04) 0.00

- 2.68 - 2.6'0 - ~. . 8 .~

- 30.87 - 30.67 -32. !7

2.005 2.(X}fi .'~.011

characteristic e x p o n e n t s requires far higher integration times. The spectra of the Lyapunov exponents are of the type ( + , 0, , ), i.e. they correspond to a strange attractor. As only one kk is positive, the attractor can be approximaled by a sheet and the topological dimension of the attractor is two [27]. In reality the strange attractor has a fractal dimension (higher than two), whic.~. can be evaluated directly fron't one-dimensional Lyapunov exponents [28]. This quantity is also given in table II. The fractal dimension increases with increasing Dt which means that the thickness of the sheet-like surface of the strange attractor also increases.

L $chreiber and M. Marek/Strange attractors in coupled cells

2.4. Birth and development of periodic orbits in

dependence on Dt The presence of a number of periodic "windows" in chaotic region A (but also outside of it) is a very interesting phenomenon, particularly with respect to behaviour of experimental systems. In diffeomorphisms thL:e windows arise by the tangent bifurcation [24]; periodic solutions bifurcate in pairs, one solution is stable and the other one unstable. Numerical examples of existence of similar behaviour in dynamic system described by a set of ODE's have also appeared [29,30]. AnG,ther type of bifurcation corresponds to the btfurcafion of stable periodic orbit with double period. Examples of such a behaviour are also well known both for discrete and continuous dynamic systems. This type of the bifurcation can repeat in an infinite sequence with universal metric properties [22]. A conjecture was made for a forced Van der Pol oscillator that an infinite number of different one-parametric families of periodic orbits exist (each of the periodic orbits being stable in a certain interval of the parameter). This implies an infinite number of tangent bifurcations. An algorithm for continuation of peJiodic solutions in dependence on the parameter [31, 25] wa~ used to deterr~ine chosen families of periodic orbits. This algorithm is based ,m the shooting method combined with a modified routine for continuation of the dependence of the stationary state on a parameter [32], and can be used for the construction of the dependence of both stable and unstable periodic orbits on DI. The method requires an integration of variational equations of the system (I). We thus have an information on the behaviour of linearized traject,~ries in the neighbourhood of the given periodic orbit in the form of FIoquet multipliers and eigenvec~ors of the monodromy matrix. One eigenwlue of the monodromy matrix is always equal to one, the remaining eigenvalues goveP n the stability of the orbit with respect to

265

small perturbations. If at least one eigenvalue lies outside the unit circle in the complex plane, the solution is unstable. The computed dependences cf periodic orbits o n the value ,)f the parameter D, are shown in fig. 5. An "average period" ~p~ = TpJk has been chosen as a r~;r~resentative characteristic of the orbit P~ (f~g given Dr). The time ~p, can be taken as an average time necessary to make one loop in the phase space or as an average time required to reach the next fixed point in the Poincar6 map. We have followed the families PI, P2 and P4 orbits (i.e. first three terms of the infinite bifurcation sequence ending at C'~, P3, Ps and/)6 orbits and, finally, P , orbits. Qualitative behaviour of all Pk families is in principle similar, and therefore we describe only parametric dependencies of the orbits P,, P~, P~ and P~. The branch Pm starts in the point of Hopf bifurcation (DI = 0.9542...) and is stable up to the point D! = 1.1720 . . . . Here one eigenvalue of the monodromy matrix, let us denote it p, leaves the unit circle at -1. It indicates, tha~, periodic orbit with double period P2 bifurcates at this point [33]. The branch PI becomes unstable with an eigenvalue p_ ~ - 1 (this branch is denoted PI_). The periodic orbit P1- becomes first more unstable with the increasing value of the parameter Dt (i.e. p_ increases in its absolute value) and then its instability diminish~ (the absolute value of P- starts to decrease) and fi~.ally O- enters the unit circle again at p- = -1. This happens for D~ = 1.4721 . . . . In the interval Dt (1.4721..., 1.4724...) is the periodic orbit (Pts) stable. The branch Pts, joins at the limit point Dt = LPI = 1.4724 with another branch of the periodic solution P,. This whol~ branch, let us denote it PI. is unstable with one eigenvalue of the monodromy matr~ (~,÷) greater than 1. The branch PI+ exists in the interval D I E (1.1219, 1.4724); on the left ooundary of its existence the solution Pt÷ coalesce with its mirror image 131÷in another limit point PB. At the point of coalescence it holds xl(t) = x2(t), yt(t) =

I. Schreiber and M. Marek/ Strange atlractors in coupletl cells

266

5.0 ~&_

.... _

_

P8

Tp i"

/

I

/

~,>.

. t'"

#,,,.

,' ~,.5-

/

/ / / 2<,I / ;/ :'/

I

t

i.-

/,

/f

/

p

/

. . . . . -f . . . . . "-t. .+'~"~ .-~ ¢, ..

j'p

~,,~-~p

,, •..~. ..... s . . . . " I .... 2.;,,.-~..:.y

"" ~ ~.'.~,- :..,.:.....'-"z: .. >-,." ..... k- --"~

t.

.'iL

."

. •" ..." p

.- %,:,

/

z,.O.

.,

/"

/

I~'

..' .,'

.. ."F2

6-.,. ...p

•

.

.. .'"

,.- , / p 2,

!.........

- "...~..::.."

T"

!

1,1

1.2

r

...................

1 ..........................................

1.3

"r- ....................................

1.t.

D1

Fig 5 The dependence of the "average period" ~,, = Ti,Jk of some families of periodic orbits of eqs. (!) on the parameter D,, A-: 2. B = 5.9. D i D : = q = 0.1. Full l i n e s - stable periodic orbits; dotted or dashed l i n e s - u n s t a b l e periodic orbits. At every point denoted by o a new periodic solution with double period bifurcates (with the exception of the "pitchfork bifurcation", the poin~ t Bk No bifurcations occur at other intersections.

y:(t). This solution is thus at the same time the homogeneous solution PH and hence a branch of PH must pass through this limit point. The computations show that to the left from the bifurcation point PB is the homogeneous solution unstable ( P ~ ) a n d to the right it is stable (PH+). Such bifurcations occur in symmetric systems ("pitchfork" bifurcation). Qualitative changes on the P~ branch starting from the Hopf bifurcation can be scheqaatically written as:

Hopf

'P,s

' P!

"

'

PIs LP~

PB ~

the iefl one (Dr= 1.1720...) we can follow the family P~ to go through iwo limit points and reach the right bifurcation point ( D t 1.4721...). Both bifurcation points on the P= branch are thus joined by one family of P2 orbits. The behaviour in the neighbourhood of the limit points on the P2 family of solutions is similar to the behaviour in th: neighbourhood of the limit point on the Pt family of solutions. We can again describe the sequen,-e of qualitative changes on the P: branch as P:s

~ P2---'-~ P2s

PI,

On the branch of PL orbits exist two points,

where the P: orbit bifilrcates. If we start at

~LP: LPz ~

P2+ P : s - - - + P 2 - - - ' * P2s

L~5ch~iber and M. MareklStrangeanractors in couplcdceils

26 ~

birth tions. Two pairs o f theqe ~ " ~ .... .... point., a r e joined by two continuous families of. P4 solutions~ (each

Similar situation can b e l o b s e r v e d for the branches of P3 orbits, The only difference is that P~ orbits form a closed loop. Schematically we can describe the changes of stability by the diagram

LP,

P3+

P.~s

P.~-

) P~s...~ L p ~ ~ L p ~

P.I÷ LPs ~

" P.~s

' ' P.~-

I ; P.~s

Similarly as it is in the case of P2 orbits, two pairs of the bifurcation points are joi.3ed by two families of P6 orbits, each of the P6 families have two limit points. It seems probable that four families of PI2 orbits again repeat the behaviour etc. The larger from the two closed loops of the P5 orbits behaves similarly as P~ loop, but the behaviour of the smaller loop is far simpler, It can be schematically described as

t the

'~

i

. . . .

TABLE III The rate of convergence 8. of the first subharmonic sequence. 8. =

(D~'- D~'-')/(D~'÷'- D~). D~ is Ps+

Pss ~"

No subharmonic bifurcations occur on this branch. It follows fr,~m 4he results of computations that in the range of values of DI, Dt E (1.0, 1.5) two types of bifurcations exist:

by

the first subh~monic sequence starting the orbit PI. The estimates of values of the p a r a m e t e r DI corresponding ....to bifurcation points have been obtained from two orbits Ph located very close t o the left and t o the right from t h e bifurcation point. A linear interpolation with respect to the eigenvalue p which passes through the point - 1 on the unit circle was used. It appears that the results shown in table III confirm the Feigenbaum's prediction of the asymptotic value 8 =~.669 . . . . The first accumulation point C~ will then correspond to the value of DI = 1.1932616 . . . . We have not verified the rates of convergence for other subharmonic sequences, but it seems reasonable to assume that accumulation points of other subharmonic bifurcation sequences will also exist. Then every periodic window will be

the bifurcation point at which the bifurcation from P2. to P2.', occui's

n

D~'

0

1.1720000

t

1.1893.093

5.s22

2 3 4

1.1923881 1.1930756 1.1932219

4.347 4.699

i

8.

i

-

•

26S

I. Schreiber and M. MarekI Strange attractors in coupled cells

limited by a tangent bifurcation on erie side and by an accumulation point on the other side. In the chaotic region A, i.e. in the range of the values of D, from the interval C~ ~ 1.1932616< D~ < F ~ 1.2635, the strange attractors appear on both sides of periodic windows. The chaotic t~,ehaviour close to the point of the tangent bifurcation cab be of the intermittent type [12]; the strange attractors close to the accumulation point are phase coherent with the structure of the power spectrum determined by the reverse bifurcation sequence. When sup~.rposition of several windows of periodic solutions occurs (or the windows are very close together as it is in the case of PW~ and PW~) we can expect an accumulation of taagent bifurcations (cf. limit points on the families of P, orbits in fig. 5) and hence also an accumulation of other periodic windows. The intervals between these narrow windows will be titled by phase coherent strange attractors. This explains the presence ol phase coherent attractors on the left from PW~ and on the right from

PW . A large number of periodic solutions still exist above the point F. where the strange attractor disappears (see table I). In all tl:ese windows is the direction of the subsequent subharmonic bifurcations from the right to the left. The last possible window is boul,ded by the limit point LP,: for the values of DI > LP, only unique periodic solution P , exists. Finally a remark on phase incoherent attractors. First window of periodic solutions Pk... with an odd k (k = 5) which has been found begins at Dt = i.201 . . . . The P~ orbits have both different "average pe,iod" from the "average periods" of the orbits P~ .o. k = I. n = 0, i, 2 . . . . . see fig. 5, and are arranged in a differellt way iv~ the phase space. For Dt~ 1.201 ... the power spectra become verv fiat. It indicates, that the P~ orbits (a,ld for higher values of D, also P~ orbits) ~ffect the flow on strange attractor and its topolt,gy in such a way that longitudinal mixing and phase incoherence arises.

2.5. Poincar~ maps

Let us consider a map Z:R3--,R 3, which relates two subsequent intersections of a solution trajectory of the system (1) with codimension one hyperplane 2 transversal to the flow (intersections in one direction are considered). This map is often called the Poincar6 map. If Z has a k-point periodic orbit, the system (1) has a periodic solution P~. Each of the k points of the periodic orbit of Z remains fixed under Z ~ and the stability of Pk is determined by the eigenvalues of DZ~ k. These eigenvalues are equal to the eigenvalues of the monodromy matrix, cf. section 2.4 (the eigenvalue equal to one is missing). We set 2: y : = 2.5, the R ~ space has coordinates (x,, Yl, x2) and we consider a projection into the plane (xl, x2) for simplicity. Subsequent iterations of Z produce in the case of a strange attractor a map which is a representation of a strange attractor. Several examples of such maps corresponding to strange attractors are shown in fig. 6a-e. We can discern between forward and reverse bifurcation sequences in Poincar6 maps similarly as in the case of power spectra. A forward bifurcation sequence P~.:~--, Pk:~.,, n = 0 , I, 2, 3 . . . . . is reflected in the Poincar6 map by subsequent occurrence of k • 2" fixed points (n = 0, I, 2 , . . . ) . Corresponding reverse sequence can be observed as gradual fusion of k . 2" islands representing, a strange attractor. The fusion occurs in a pairwise manner and k . 2 " fixed points of now unstable peric,Jic orbit Pk :~ are located approximately in the middle of the.,;e islands. The final state has k islands. When reverse sequences overlap, the number of islands can change. This is probably the case of the Poincar6 map in fig. 6a which has two islands instead of four (as could be expected from corresponding spectrum 3a). The Poincarg map in fig. 6b (corresponding to spectrum 3b) has one (more or less) continuous island (,k = 1). The Poincar6 maps 6c, d, e belong to phase

!. Schreiber and M. Marekl Strange attractors in coupled cells

269

ties o f two-dimensional H6non map [ 15, 34]

,if

°

T(x, y) = (l + y - ax ~, bx).

(2)

\

incoherent attractors. Here only one island appears in contrast to the Poincar~ ,naps of the pha,e coherent strange attractors (with an exception of the Poincar6 map shown in fig. 6b all ~hase coherent strange attractor Poincar~

Our map Z is three-dimensional and is directly connected with the s y s t e m ( l ) which makes the analysis more complicated, However, we c a n make use of ideas developed in the course of analysis of the H6non map to explain disappearance of the strange attractor at D, = 1.2635. In the case of the H6non's map it is known [35, 15] that the basic role in the process of formation of attractors and their domains of attractions play two hyperbolic fixed points P., P_ of the map Z and their stable and unstable invariant manifolds. Because both fixed points have a character of saddles, they have one stable W s and one unstable W U manifold (both are one-dimensional). The points in W~+ N wS+ of the manifolds of the point P+ form a set of horn,clinic points which have the structure of a Cantor set. This set is then responsible for the appearance of the strange attractor. The stable manifol(~ W,sp_ of the second fixed point P_ forms a boundary of the domain of attraction of the strange attractor. The structure of a strange attractor and various bifurcation phenomena (fusion of islands, appearance of attracting n-periodic orbits in chaotic region) can be explained by means of stable and unstable manifolds of fixed points of the map T". A strange attractor formed by iterations of ,' disappears when outer heteroclinic points appear, i.e. when W s_ N W~. ~ ~. Appro.dmately, the strange attractor loses its attractive character at the moment when it touches the boundary of its own ,~omain of

m a r s harP- m a r ~

tattti

'I I

....

3

T

4

I"-- . . . .

I". . . . . - " T - - - ' - - I

5

6

~

I [

.,e.l

em~"

,."

I"- ..... I-- .... -'t-- .... ~.....

7x~8

4

S

7x~ 8

G

\

'!

1-

\

]

41

° •

"

I

3

4

5

G

3

7x~ 8

t,

5

fi

7x~ 8

......................;, l X

%,

ii. . . . . . . . . . . . . P "~,.

'~.

•

3

4

S

6

! , i

7x 1 8

Fig. 6. Projection of the Poincar6 map into the plane (x~. x,). Figs. 6a. b. c, e correspond to figs. 3a. b. c, g. respectively. Increasing thickness of the sheet can be seen in figs. 6d, e. A = 2, B -- 5.9, D J D , = 0 . 1 . (a) O~ = 1.194; (b) D~ = 1.20; (c) DI = 1.21; (d) Dt ffi 1.24; (e) DI = 1.263.

kl~nd~).

From figs. 6c, d, e it can be seen that the strange attractor is stretched on its upper e~sd when DI is increased. This phenomenon leads to appearance of heteroclinic points and the strange attractor is destroyed. Poincar~ maps arising in iterations of Z have mP,~y properties qualitatively similar to proper-

~ttrehot~nn

t4,1w t i,t~ i h

Our map Z has in the i~lterval D, E (1.1219, 1.4724) five fixed points corresponding to periodic sonutions P,, PI.,/'1-, Pm. 1~-. Because of the symmetry in t.he system (!) we need not follow fixed points Pt., /31- (here v,e keep the same notation both ior closed orbits of (1) and for correspo~ din8 fixed points of Z, even though

. ,0

I. Schreiber and M. MareklStrange attractors in coupled cells

the symmetry of (1) does no~: hold for Z owing to our specific choice of 2;). F:om the previous section we know that whenever a closed orbit P~ is unstable, lhe corresponding monodromy matrix has exactly one e~genvalue located 9:ttside and two eigenvalues located inside the ~mit circle. The same holds for the matrix DZ~,~. This implies that every unstable fixed point P~_. of Z ~ is hyperbolic and has one-dimensional unstable manifoid W~ and two-dimensional stable manifold W~,. For example, all computed fixed points P~. for D~ = 1.21 (with an exception of P~_), i.e. P~_., P : (P:. does not exist for this valoe of Dr), P~.,, P~- and Pn are shown in fig. 7. If this figure is compared with fig. 6c, which represents the strange attractor, we can follow evident connections between the set P~_. (with an exception of P~.~ and the strange attractor. Hence we can expect that there exists a set of homoclinic points w S, n W ~ and that the strange attractor is a closure of unstable manifolds of hyperbolic reriodic orbits [36. 6, t5]. The point P~,, which is located between the

6

point PH and SA is the key to explanation of the disappearance of the s~range attractor. If Z behaves in a similar rr~anner as T, then the stable manifold wS~. fo~ms a boundary of the region of attraction of SA. The solutions outside of the region of attraction then would ~ attracted to Pn, eventually to a symmetric SA. We have verified the idea about WS~÷ locally in the neighbourhood of Pt+. Iterations of Z with initial conditions located on the edges, corners and walls of a small cube centered at Pt÷ have been followed, cf. fig. 8. Successive iterations of Z from one part of the cube approached SA, the iterations from the second one approached PH, see fig. 8. The boundary between both parts of the cube is just wS~÷. We can observe in fig. 7 that for D - 1.21 is the point Pt÷ still far away from SA. When the value of D~ is increased, Ps+ and SA approach each other until SA touches vW sp,, i.e. heteroclinic points wS,~ n W~_ appear. In this way a "gap" arises in the boundary of attraction of SA. Any "Jrbit will eventually escape through the gap and be attracted to P,. The time interval spent in the chaotic motion can be, however, very long. This happens just at the point F = Di I.,.63..

• P

V~ u

o

o

÷ r

o~AO

°

•

\ sA \

r

2

~,

~ . . . . .

6

!

........

X~

Fig. 7. Fixed points of Z k (k = 1, 2, 3, 4, 5) for Dt = |.21. PLaint PH is stable, t h e o t h e r points are u n s t a b l e . Points P t , P , P~ , P , . , P~ , P~, P~. are parts o f t h e SA, see fig. 6c, The ~oint P~. I;e~ in the b o u n d a r y of a t t r a c t i o n o f the SA.

Fig. 8. Stable manifold W~,I, of the point Pt, for D i - 1.21. Points in the vicinity of PI, tend under iterations of Z to PH from one part of the cube and to the SA from the other one.

L $chreiber and M. MareklStrange attraclors in coupled cells

271

Because heteroclinic ,points a p p e a r in, the vicinity Of P t + i n our case, we ~can v e ~ the above mechanism in the following way, We choose some points on SA for the critical value

3. Discussion

these points p a s s : see fig, 9a, b)i.e, heteri Dynamic behaviour of the system(i)slightly above F is similar as in the Lorenz model [37]. Starting from initial conditions (b) the system oscillates chaotically for some time until it reaches the "gap" in wSl+ and is attracted to homogeneous periodic solution Pn. The closer we are to F the longer is the "chaotic phase". This behaviour is known as metastable chaos [;!]

sequence~ of subsequent: period doubling bifurcations occurs and chaos a r i s e s a f t e r t h e accumulation point (C~ = 1'1932616). When D! is decreased, metastable chaotic oscillations occur. The metastable chaos changes to stable chaotic motion when outer heteroclinic points disappear (F = D! = 1.2635). In the chaotic region D~ ~ (C~, F) periodic and strange attractors alternate. The transition between them can be via a tangent bifurcation or via a cascade of subharmonic bih~rcations. The flow on strange attractors is sensitively dependent on initial conditions and can be either statistically phase coheren! or incoherent. Above the point F there still exist a number of periodic windows with cascades of subharmonic bifurcations. The type of chaotic and periodic states studied here can be called "in phase" type because the phase shift between both oscillators is very small. This seems to happen for a relatively stronger coupling be'~ween cells while other types of chaotic and periodic states are possible when weak coupling is considered. The study of these states will be a subject of another communication.

a)

O1 < F

/ " , . '•QI? /

b)

f

O~>

/wl.

t"

F

/

-,

\~,S

-

f

/~

We have studied in more detail a system which is of interest in many fields- a system of

"q

Acknowledgements Fig. 9. Mutual position of the stable manifold wS,, ,tad the point S on the upper edge of the SA at the critical value Dt = F = 1.2635. (a) Dt < F, the point $ forms a part of the SA; (b) DI > F, the point $ have crossed through the eb,.,~plane W~,, and the SA vanishes via a metastable .naos.

Discussions with A. Kli~, M. Kubl~ek, M. Holodniok and the referee's comments are gratefully acknowledged.

272

I. Schreiber and M MureklStrange attraclors in coupled cells

References

[111 J. Yorke and E. Yorke, J. Stat. Phys. 21 (1979) 263.

[1] !. Prigogine and R. Lefever, J. Chem. Phys 48 (1968) 2¢,3. T, Pav)ldis, Biological Oscillators: Their Mathematical Analy,~is (Academic Press, New York, 1973), R, Aris, The Mathematical Theory of Diffusion and Reaction ;n Permeable Calalysts (Clarendon Press, Oxford, ! )75). P, C. Fit,-, Mathematical Aspects of Reacting and Diffusing Systems (Springer, New York, 1979). 121 M. Marek, in Synergetics-Far from Equilibrium, A. Pacault and C. Vidai, eds. (Springer, New York, 1979). M. Marek and E. Svobodov:i, Sci. Papers of the Prague Institute of Chemical Technology K 10 (1976) 199. M. Filipovh, Cascade of Reactors with Mutual Mass Exchange, Thesis, 1978, Prague Institute of Chemical Technology. M. Marek and I. Stuch|, Biophys. Chem. 3 (1975) 24. M. Marek and E. Svobodovfx, Biophys. Chem. 3 (1975) 263. 13a] 1. Stuchi and M. Marek, J. Chem. Phys., to appear August !%2. [3hi i. Stuchl and M. Marek, J. Chem. Phys., to apgear September 1982. [4al i. Stuchl and M. Marek, preprint, CHISA Congress, Prague 1978. [4b] I. Schreiber, M. KubR:ek and M. Marek, ia New Approaches to Nonlinear Problems in Dynamics, P. J. Hoi,~es ed. (SIAM, Philadelphia, !~J80) p. 496. [5} K. Tomita and T. Kai, Prog. Theor Phys. 61 (1979) 54; J. Stat. Phys. 21 (1979) 65; M. L. t/artwright and J. E. l.ittlewood, J. l.ond, Math. Soc. 20 (1945) 180. [6,] Y. Ueda. Ch. Hayashi, N. Akamatsu, Trans. IECE .l~pan. 56-A (1973) 218. Y Ueda. J. Stat. Phys. 20 (1979) 181. [7!. ~. Guckenheimer, Physica ID (1980) 227. [8] (). E. Rossler. Z. Naturforsch. 31a (1976) 1168. [9] C. R. Kennedy and R. Aris, in New Approaches to Nonlinear Problems in Dynamics, P. J. Holmes, ed. (SIAM, Philadelphia, 1980) ~. 21 I. [10] H Fujisaka and T. Yamada, Phys. Lett. 66A (1978) 450: Prog. Theor. Phys. Supk~l. 64 (1978) 269; J. J. Tyson, J. Chem. Phys. 58 (1973) 3919.

L. Kaplan and J. Yorke, Commun. Math. Phys. 67 (1979) 93, [12] P. Manneville and Y. Pomeau, Physica ID (1980) 219. [131 E. N. Lorenz, J. Arm. Sci. 20 (1963) 130, [,~41 R. M. May, Nature 261 (1976) 459. [;5] C. Sire6, J. Star. Phys. 21 (1979)465. [161 I. Shimada and T. Nagashima, Prog. Theor. Phys. 59 (1978) 1033. T. Nagashima, Prog. Theor. Phys. Suppl. 64 (1978) 368. [lS] R. J. HiBgins, Am. J. Phys. 44 (1976) 766. W. Cooley and J. W. Tukey, Math. Comput. 19 (1965) 297. [19] I. Shimada and T. Nagashima, ProM. Theor. Phys. 61 (1979) 1605. [2o1 G. Benettin, C. Froeschle and J. P. Scheideeker, Phys. Rev. AI9 (1979) 2454. [21] G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems (Wiley, New York, 1977). [221 M. J. Feigenbaum, J. Star. Phys. 19 (1978) 25; J. Star. Phys. 21 (1979) 669. [23] E. N. Lorenz, Ann. N.Y. Acad. Sci. 357 (1980) 282. [241 D. Farmer et al,, Ann. N.Y. Acad. Sci. 357 (1980) 453. [251 1. Schreiber, A. KIf& M. Kubf6ek and M. Marek, preprint CHISA Congress, Prague 1981, 1261 V. I. Oseledec. Trans. Moscow. Math. Soc. 19 (1968) 497. [27] N. H. Packard et al., Phys. Rev. Lett. 45 (1980) 712 I281 H, Mori, ProM. Theor. Phys. 63 (1980) 1044. [29] V. Franceschini and C, Tebaldi, J. Stat. Phys. 21 (1979) 707. I30l V. Franceschini, J. Stat. Phys. 22 (1980~ 397. !311 M. Holodniok and M. Kubi~,ek, to be published; M. Kubf~,ek and M. Marek, Computational Methods in Bifurcation Theory and Dissipative Structures, l,ect. Notes in Comp. Physics, ¢Springer, Berlin, New York 1982). [321 M. Kubf~ek, ACM TOMS 2 (1976) 98. [331 G. Iooss, Bifurcation of Maps and Applications (NorthHolland, Amsterdam, 1979). [3q M. Henon, Commun. Math. Phys. 50 (1976) 69. [351 H. Daido, Prog. Theor. Phys. 63 (1980) 1190. [36] S. E. Ne::'~:use, Ann. N.Y. Acad. Sci. 357 (1980) 292. [371 V. S. Aframovich, V. V. Bykov and L. P. Shilnikov, Dokl. Akad. Nau*: SSSR 234 (1977) 366.