- Email: [email protected]
Predicting pattern formation in coupled reaction—diffusion systems
Chemical Engheering Science. Vol. 44. No. 9, pp. 1893-1901. Printed m Great Britain.
1989.
am-2soV/nv s3.00+0.0 0 IV89 Pergamon Press plc
PREDICTING PATTERN FORMATION IN COUPLED REACTION-DIFFUSION SYSTEMS IOANNIS Department of Chemical (Received
G. KEVREKIDIS and HARRY S. BROWN Engineering, Princeton University, Princeton, NJ 08544, U.S.A. 10 November
1988; accepted 2 February
1989)
Abstract-Coupled reaction4iffusion equations are known to exhibit a wealth of multiple coexisting stationary solution patterns as the characteristic length of the system grows. We describe and implement a technique which allows us. by studymg only stationary solution branches at small system lengths, to predict the complex structure of steady-state bifurcations occurring at large system lengths without actually computing this structure. This technique is applicable to arbitrary isothermal or nonisothermal systems of coupled reaction-StTusion equations with periodic or no-flux boundary conditions. To illustrate it we use a standard literature example (the Brussellator), where. by computing only the first two solution branches, we accurately predict the steady-state bifurcation structure reported up to much larger system lengths. This technique also provides a compact way ofdescribing and comparing stationary pattern formation in a large class of systems, extending beyond coupled reaction-diffusion equations. We demonstrate this by comparing stationary pattern formation for our test problem (the Brussellator), with the formation of complex surface wave patterns in thin liquid film flow.
INTRODUCTION A number of nonlinear reaction-diffusion systems can have more than one steady-state solution; some of these solutions may be spatially nonuniform. This property can be used in modeling a variety of pattern fcirmation phenomena in biochemical and chemical contexts, from morphogenesis and cell differentiation in embryonic development to surface concentration patterns in catalytic kinetics. This property is also closely related to the emergence of spatial nonuniformities (asymmetric steady states) in systems of interacting catalyst particles. Many possible models of interest exist, giving rise to different reaction rate functions and different bifurcation diagrams as the important parameters (the kinetic rate constants, the diffusion coefficients and the characteristic system size) vary. In particular, it is well known by numerical studies on a number of specific models [from, for example, Babloyantz Hiernaux (1975), and Herschkowitz-Kaufman (1975) and Kubicek et aZ. (1978) to, more recently, Duncan and Eilbeck (1987)] that the steady-state bifurcation diagrams for isothermal, coupled reaction-diffusion systems with noflux boundary conditions often possess a bewildering complexity, which increases as the characteristic system size increases. This complexity results from the interplay between an intrinsic pattern size (influenced by the reaction rate form, the values of the kinetic constants and the diffusion coefficients) and the characteristic system size, imposed by the boundary conditions. While the results of this interplay are indeed model-dependent, one would expect the general mechanisms of pattern formation to be reasonably robust to changes in the models. In this paper we describe and implement a technique which allows us, by studying stationary solution 1893
branches at only small often complex structure much larger system sizes this structure. This can
tational behavior.
system sizes, to predict the of bifurcations occurring at without actually computing lead to significant compu-
savings in accurately describing the system The technique is applicable to arbitrary
nonisothermal systems of coupled equations with no-flux or periodic Its general applicability, boundary conditions. coupled with the ease of implementation and the resulting computational savings, make it a useful tool in predicting and understanding the detailed steadystate bifurcation behavior of any specific model. We believe, however, that its usefulness goes beyond mere computational savings in a particular model study. Since an entire complex pattern-forming bifurcation process can be reduced to the study of a few branches, isothermal or reaction4ffusion
this technique provides a compact way of describing and comparing pattern formation across several models. This is an efficient approach to quickly and accurately assessing the effect of changes in the model to the predicted steady-state behavior over large intervals of parameter values without actually performing detailed computations. It may also prove helpful in understanding the degree to which this apparently complicated process of pattern formation is robust and model-independent. This paper is structured as follows. In Section 2 we describe the conceptual underpinnings of the method, and present the computational steps involved in its implementation. In Section 3 we ilhistrate this method by applying it to a popular literature example: the “Brussellator”. The main reason for choosing this as a “typical” example to illustrate the applicability of the method is that many chemical engineers are probably familiar with its behavior, since it has been the object
IOANNIS
1894
G.
and
KEVREKIDIS
of intense numerical calculations in several scientific papers, and it has also appeared as an illustrative example in monographs [e.g. Kubicek and Marek (1983)]. A different example (i.e. one with a rational polynomial rate expression, like the ones resulting from enzyme or Langmuir-Hinshelwood kinetics) could serve equally well to illustrate the method. in this example, we are abte to accurately predict the secondary bifurcation structure observed and reported in the literature for large characteristic system sizes, by studying only the first two steady-state branches at small system sizes. Finally, in Section 4, we use the method to discuss the robustness of the pattern formation process in different models. This is accomplished by comparing the common features of the bifurcation diagram for our test example-the Brussellator-with the formation of complex surface wave patterns in thin liquid film flow. 2. THE METHOD Mnf ivation
To motivate our computational procedure we consider for simplicity a set of two coupled reaction-diffusion equations, occurring in a one-dimensional medium. The steady-state equations are
D”D,, +$I(% 0) = 0
(1)
where u and v represent concentrations of the two species, D, and D, are diffusion coefficients, and f (u, u) and g(u, o) are the reaction rate forms. Let us for the moment consider periodic boundary conditions [u(O) =u(L), u(O)=v(L)] and suppose that system (1) has, for some set of parameter values and some system size L, a certain “primary” solution u(x). v(x) (0 G x < L). It is then obvious that if we double the system size and
31
I 02
S. BROWN
extend the steady-state pattern two-fold [u(x + I_) =u(x), V(X+ L)= v(x)] we obtain a soIution of the same equations with the same boundary conditions for the new system size 2L. This holds in general for the n-fold extension of this first solution, which will be a solution of the same equations with the same type of boundary conditions when the system size in nL. This by no means implies that the nth replication of the “primary” solution in the only solution of eqs (1) when the system size in nL. Neither does it imply that the nth replication of the “primary” solution has the same stability characteristics for system size nL that the “primary” solution has at system size L. This is easily seen upon inspection of the bifurcation diagram reported in the literature for the “Brussellator’‘-see Fig. 1. At L= o! the “primary” branch A, of spatially nonuniform stationary solutions bifurcates from a spatially flat solution. This implies then that at L = 20~ the two-fold replica A, of A, bifurcates from the spatially flat solution, that at L=3a the three-fold replica A, of A, bifurcates from the spatially flat solution etc. We observe, however, that, while soon after its birth, the “primary” branch A, is the only spatially nonuniform stationary solution, its two-fold replica A, over most ofits interval of existence coexists with several other spatially nonuniform stationary solutions. Another important observation is the following: the only two steady-state bifurcations observed on the diagram for the “primary” branch AI are its birth (at L=a) and its “death” when it collides with the A, branch at L=c. By the simple “replication” arguments, this means that the two-fold replica of A, will be “born” at L = 21x,and it will “die” colliding with the four-fold replica of A, at L = 2c, which is indeed o,bserved. This was observed, for example, by Kubicek et al. (1978), and is related to what they termed “composition of solutions”. We notice, however, that
I 04
I 06 Length
Fig. 1. Steady-state periodic boundary
HARRY
I 08
I I
L
bifurcation diagram with respect to the system size L for the “Brussellator” with conditions up to L= 1 for A=2, B=4.6, D,= 1.6 x 10-j, D,=8 x 10m3.Only solutions symmetric about x = L/2 are shown.
Predicting pattern formation in coupled reaction-diffusion the two-fold
replica
also connects
with several
other
branches at several (four, to be exact) “secondary” bifurcation points, for which there is no counterpart on the “primary” branch A,. Furthermore, it is obvious even from the part of the diagram shown, that the four-fold replica A, of the A, branch will have even more secondary bifurcations, giving birth to a large number of so-called “mixed-mode” branches. Our goal in this paper is to present a simple technique for accurately predicting, as the system size L grows, all these secondary symmetry-breaking bifurcations off of all the n-fold replications of a given “primary” branch by only studying the “primary” solution branch, which in itself does not exhibit them. This technique is generally applicable in systems modeling pattern formation with periodic boundary conditions [see, for example, Scovel et al. (1988)]. What makes it particularly interesting in the case of reaction-diffusion systems, where no-flux boundary conditions are often the reasonable choice, is the following simple observation. Consider a solution u(x). u(x) of eqs (1) with no-flux (Neumann) boundary conditions at x = 0 and x = L/2: u,=u,=O If we extend
at x=0
this solution
and x=-.
between
L 2
L/2 and L so that
we obtain a solution of system (1) with periodic boundary conditions at x=0 and x= L. This means that u[Z stationary solutions of system (1) with Neumann boundary conditions for system size L/2 are also stationary solutions of eqs (1) with periodic boundary conditions for system size L. The converse is, of course, not necessarily true, i.e. there may exist solutions of eqs (1) with periodic boundary conditions that are not necessarily also solutions of the Neumann problem at half the system size. The method we describe can, of course, be used to study and predict the bifurcations of these latter “periodic” but LLnonNeumann” solutions when they exist. This simple property of the coupled reaction-diffusion equations, however, allows us to use the method, which hinges on periodic boundary conditions, also in the case of Neumann boundary conditions for reaction+liffusion systems. Indeed, all the stationary solution branches shown in Fig. 1, which have been computed with periodic boundary conditions for various system sizes, have Neumann counterparts at half the system size shown.
L = 2a, and collides
systems
1895
with A, at L = 2~). It is straightfor-
ward to predict that these bifurcations will occur, and the parameter values where they will occur, by simply knowing their “counterparts” for the primary branch A,. As a matter of fact, it is obvious that one can predict the “replications” of these bifurcations for the three-fold replication of A,, the branch A,. This branch will be born at L= 3a, and will collide with A, at L=3c, beyond the values of L shown. Similar results hold for all other n-fold replications of A,. There exist, however, four. parameter values at which A, exhibits steady-state bifurcations for which A, has no counterpart. We are interested in the possibility of the predicting these “secondary” bifurcations by only studying the primary branch. When, For practical computational purposes, we discretize a problem like eqs (I), say through finitedifference or pseudo-spectral techniques, the approximation of its steady-state solutions reduces to the solution of severai (say N) nonlinear simultaneous algebraic equations. These solutions can be followed using standard continuation techniques [e.g. AUTO (Doedel 1981, 1986)] as the system parameters (here the system size L) vary. To accurately calculate steadystate bifurcation points one searches for parameter values for which an eigenvalue of an N x N square matrix (the linearization of the N nonlinear algebraic equations around their solution) crosses zero. Such points are then further examined to determine whether a steady-state bifurcation does indeed occur there. “Secondary” bifurcations, occurring off of replications of a primary branch, but having no counterpart on the primary branch itself, become easier to understand and predict in the following simple fashion. Consider a steady-state solution on the primary branch for some parameter value (say L ~0.2) [Fig. 2(a)) (remember that “half” of this solution satisfies the Neumann boundary conditions for half the system size). Then at L=O.4 the two-fold replication of this solution [Fig. 2(c)] will be again a solution, lying on the branch A,. A sready state of problem (1) with periodic boundary conditions and system size L can be considered as a limit cycle (with period L) of the following set of ordinary differential equations (ODES) (in the spatial variable x):
du, dx
*
4 dx
4
The technique and its implementation In Fig. 1 the two-fold replication A, of the primary branch AL does possess the replications of the steadystate bifurcations that A L exhibits (i.e. A, is “born” at
where (u,, ul, uj, u,)~=(u, u,, u, u,)~. A projection of this limit cycle on the (u,, ul) plane (i.e. u and u,) is
IOANNIS G. KEVREKIDIS and HARRY S. BROWN
1896
2fb) ;
El 0
-0
4
”
branch is now apparent. Suppose that, for some value of L, the steady solution of eqs (1) lying on A, haswhen viewed as a limit cycle of eqs (2 ja Floquet multiplier at the nth root of unity (say for example at the square root of unity, - 1, or at a fourth root of unity, &-i). Then, at n times that system size, n_L, the nth replication of the primary branch A, will have-when viewed as a limit cycle of (2) -a Floquet multiplier exactly at 1, and that is a potential bifurcation point of limit cycles for (2), or steady states of (Itpotential because one still has to verify that a bifurcation will indeed occur there. Consider for example all j&r ‘*secondary” bifurcation points observed on branch A,. They occur at precisely twice the system sizes for which the corresponding solution on A, had-when viewed as a limit cycle of (2kFloquet multipliers at
-1. The implementation
Fig. 2. Spatial profile of a solu&ion at L = 0.2 [Fig. 2(a)] and its representation as a limit cycle in x with period L=O.Z projected on the u-u, plane [Fig. 2(b)]. Figure 2(c) shows a two-fold replica of the same solution at LzO.4. Species u: solid lines; species ~1:dotted lines. plotted in Fig. 2(b). Moving along the solution curve as x varies between 0 and L=O.2 in Fig. 2(a) is equivalent to moving the phase point [u(x), u,(x)] once around the limit cycle in Fig. 2(b). It is now obvious that the two-fold replication of this solution when L=O.4 is exnctly the same limit cycle, only the phase point goes around it twice. This is the crucial observation on which the technique is based: steady states of system (1) with periodic boundary conditions at system size L are Ziinit cycles of the. set of ODES (2) steady-state bifurcations of with period I.. Therefore, problem (I) to other steady-states can be viewed as bifurcations of limit cycles of problem (2) to other limit cycles. A necessary condition for bifurcations of limit cycle branches to other limit cycles (with nearby periods) to occur is that an eigenvalue of the state transition matrix (the so-called Floquet or characteristic multipliers) must cross unity [see, for example, Iooss and Joseph (1980)]. This means that, at parameter values L for which the branch A, does exhibit a steady-state bifurcation (while branch A, dots not), rotating around a limit cycle twice (solution on branch AZ) yields a Floquet multipher at unity, while rotating around the limit cycle only once (solution on branch A, at half the system size) does not. There exists, however, a very simple relation between the Floquet multipliers obtained when we go around a limit cycle once, and those we obtain when we go around the same limit cycle twice. Namely, going around twice will syuarr the Floquet multipliers (similarly, going around 3 times will raise them to the third power, and going around n times will raise them to the nth power). The technique for predicting the secondary bifurcations by only using information on the primary
of this technique is therefore very simple, and can be appended with minimal effort to any continuation code currently used: after a solution on the primary branch (with periodic boundary conditions) has been computed-through any discretization method-we simply integrate the variational equations of system (2) around this solution once. We will obtain a (2n) x (2n) state transition matrix (where n is the number of coupled reaction-diffusion equations, in our test case n = 2). Its 2n eigenvalues (in our test case 4) are the Floquet multipliers; one of them has to always be at unity [see, for example, Iooss and Joseph (1980)]. When one of the remaining multipliers goes through an nth root of unity in the complex plane, then-at n times that system size the nth replication of the primary branch will potentially have a “secondary” bifurcation. These bifurcations give rise to the so-called “mixed-mode” branches.
3. THE EXAMPLE The equations We will illustrate the method by applying it to a system of two coupled reaction-diffusion equations, the “Brussellator”: U,,+;CA+U%-_(B+ U
l)u]=O
I u,+~(Bu-u%)=O.
In particular, we will examine the steady-state bifurcation diagram obtained by numerical computations for A=2, B=4.6, D,=1.6x lo-’ and D,=8x 10e3 with periodic boundary conditions. This diagram (the one presented in Fig. I) may be more familiar to the
reader in the form it (1978) or Kubicek Fig. 3(c)]. The reader after our discussion
was presented by Kubicek et al. and Marek (1983) [see also should
however keep in mind, 2, that the branches reported here at a given system size have no-flux counterparts, reported by Kubicek et al., at half that system size. To further facilitate the comparison, we
in Section
1897
Predicting pattern formation in coupted reaction-diffusion systems
a
2a
I 02
Length
0
I 01
L
I
I
I
0.2
0.4
0.6
Length
13
I 0.4
03
I a*
1 I
L
Fig. 3. (a) Two “primary” branches. At specially labeled points the principal Floquet multipliers are located at nth roots of unity on the unit circle (n = 1 is symbolized by a, n = 2 by b, n = 3 by t, n = 4 by q etc). The same symbols preceded by a B denote such occurrences on the second branch B,.(b) Secondary bifurcations on replicated branches that can be predicted based on Floquet multiplier information on the primary branches. (c) Same bifurcation diagram in the format reported by Kubicek et al. (1978): instead of the L* norm of the solution,
u (0) is plotted
have only included in our diagram solutions of the problem with periodic boundary conditions that do possess no-flux counterparts. As we will discuss, there exist bifurcations to periodic solutions that do not possess no-flux counterparts. We will perform Floquet “primary” calculations on the two multiplier
vs t.
branches, A, and E,, and from that information we will predict all other secondary bifurcations in the diagram and several others that will be observed at values of L far beyond the ones reported here. There exists only one more technicality we must discuss before examining the bifurcation diagram in
IOANNIS G. KEVREKIDIS and HARRY S. BROWN
1898
Fig. 3. In our test problem, and every problem involving two coupled reaction diffusion equations, we will calculate four Floquct multipliers [since we have four ODES in system (2)]. We observe computationallyand it is possible to show theoretically-that, for the stationary solutions we discuss in Fig. 3, two of these four multipliers (let us call them E., and ia) will always be fixed at unity, and the remaining two, & and E.,, which we will call the principal Floquet multiphers (PFMs) will vary as the system size (and the entire solution) varies, while their product always remains at 1. This means that they may either be both positive or both negative such that A,= I/>.,, or that they are complex conjugate with modulus 1 (i.e. they will lie on the unit circle in the complex plane). Since this fact is not germane to the technique, the reader may simply consider it as a computational observation. Finally, in what follows, we have tried IO use a convenient mnemonic notation for the parameter values at which the PFMs lie at roots of unity, and for the predicted bifurcations. Points on a primary branch, for which the PFMs lie at an nth root of unity are marked b, (for square roots), ti (for cube roots), qi (for fourth rootstwhere i is an index to discriminate between them. The corresponding predicted bifurcation points on replications of the primary branch are marked with rhe same symbol preceded by an R, (R stands for replication and k denotes a k-fold replication of the primary branch). For example, b, denotes the first time (L) on the primary branch A 1 for which we observe a PFM at the square root (b) of unity. This predicts a bifurcation off of the two-fold replication (R,) of A, at twice (I) the system size (RA). The primary
branch
A, (Spatially At L= LX = 0.159689 the first nontrivial varying) branch bifurcates off of the flat solution U(Y) = A, u(x) = B/A. Figure 3(a) shows an enlargcmcnt of the bifurcation diagram of Fig. 1 between L-=0.1 and L=O.5, obtained through a spectral discretixation and a modification of the bifurcation code AUTO written by Doedel (1981). Figure 3(c) is the same bifurcation diagram reported [as Kubicek et al. (1978) did] by plotting u(O) vs the bifurcation parameter L. At L = o! =0.159689 the two principal Floquet multipliers lie at i 3,4 = - 0.6428 f 0.76601, or at an angle of approximately 130” (see Fig. 4, and continue referring to it). As L grows, these two multipliers move along the unit circle in the complex plane, and they go through the fifth root of unity at L = 0.16527 (point f,), which will then mean that the fifth replica of this branch, marked A, in the diagram, will have a symmetry-breaking bifurcation at L = 0.826 % 5 x 0.16527 (R,f,). Indeed, that is a bifurcation point of the fifth replica A, of A,. From now on, we will only consider bifurcations of replicated solutions for k,<4 (up to fourth roots of unity), since higher replication branches lie mostly outside the diagram of Fig. 1. Subsequently the PFMs come to a double - 1 at L =0.17404 (point b,)-see again Fig. 4-which pre-
Fig. 4. Structure and direction of motion of the principal Floquet multipliers as the parameter L changes at the points labelled
on the first primary
branch
A, of Fig. 3 (a).
diets a symmetry-breaking bifurcation for branch A, at L=O.3481 (point R,b, where branches A, and A, collide-and further replica at R,b, where branch A, runs into branch A, in the same fashion that branch A, runs into branch A, at R,b,). After that, the PFMs become real and split on the real axis, one inside and one outside the unit circle. They return to a double - 1 at L=O.18363 (point b2), which predicts a symmetrybreaking bifurcation for branch A, at L = 2 x 0.18363 = 0.36726 (point R,b2, bifurcation of branch A, into a “new” secondary branch which we will call I?-with a They then go further replica at R,b,, L=O.7344). through the third root of unity [&, = - 0.5 f 0.5iJ3, or exp ( + i2z/3), 0 = 120” (see Fig. 4)] at L = 0.19423 (point t2), with a predicted symmetry breaking bifurcation on branch A,
Predicting pattern formation in coupled reaction+Iiffusion which again there is no bifurcation in the subspace of functions with the property that u(L/2 - x) = u(L/2 +x). At that point we do, however, have a bifurcation to an asymmetric (around L/2) periodic branch, which has no Neumann counterpart at half the system size. They now go again around the unit circle. They pass through the fourth root of unity at point q3 (L = 0.3409, symmetry-breaking bifurcation on the A, branch at L = 4 x 0.3409 = 1.3634, point R,q,, off scale in the diagram). They pass through a cube root of unity (point t,, L = 0.34 1146, symmetry-breaking bifurcation at L= 3 x 0.341146= 1.023438 off the A, branch, point R,t,, again off scale). They hit the double-negative one at L=O.341382 (point b3), with a symmetry-breaking bifurcation for the A, branch at L =0.68276 (point R,b,). They split on the real axis (both negative, one inside and one outside the unit circle) and they return to a bifurdouble ~ 1 at L = 0.35430 (point b,, predicted cation at L= 0.7086 off of branch A,, point R,b,). They go through the cube root of unity at L=O.35359 off of the A, branch at (point t,, predicted bifurcation L = 1.0608, off scale in the Fig. 4) and through the fourth root of unity at L~O.35259 (point q4, predicted bifurcation off of A, at L=1.4104, off scale). They finally come to a quadruple 1 at point u,-R,b,, where the A, branch merges with A,.
The primary
branch 9, This branch also bifurcates from the flat solution, at L=O.442243 (point p) and, like the A, branch, it eventually connects with the A, branch. The principal Floquet multipliers on this branch start at point fi at as is predicted 0~83.0”) E.,,,=0.1216f0.9926i(angle from linear analysis around the flat steady state. They then start moving around the unit circle (as L decreases) and they reach 90” (a fourth root of unity) at L bifurcation =0.441443 (point Bq,, with a replicated of the replicated branch B, at L= 1.7658, off scale). They reach 120” at L=O.43778 (point Bt,, with a symmetry-breaking bifurcation of the replicated branch B, at L= 1.31334, again off scale). They come to a double ~ 1 at L=O.431876 (point Bb,), which predicts a symmetry-breaking bifurcation off of the replicated branch B, at L = 0.863753 (point R,Bb,, birth of branch Dl). They split on the negative real axis (one inside and one outside the unit circle) and return to - 1 at L=O.42645 (point Bb,, with replicated symmetry-breaking bifurcation at L=O.85290, point R,Bb,, death of branch DI). They traverse the unit circle, going through the third and fourth roots of symmetryunity (points Bt2 and Bq2, with replicated breaking bifurcations off of branches B, and B,, respectively, off scale in the graph). They come to a quadruple multiplier at 1 at point Bu, (L=O.4089) split on the real axis and come back again to a quadruple multiplier at 1 at point Bu, (L=O.37624). As was the case with the A, branch, no bifurcation is observed for one of these quadruple 1 points, while an asymmerric branch, not depicted in the diagram, is
systems
1899
expected at the other quadruple 1 point. The PFMs subsequently start moving on the unit circle, and now a ‘&new” phenomenon is observed: instead of moving monotonically, they reach a “maximum” angle at approximately 24” (LzO.37157, point S) and then return towards a quadruple 1, where they arrive at L =0.367108, point Bu,, which has been predicted as a symmetry-breaking bifurcation off of branch A, (point L=2xO.l8355) (remember, at L=O.18355, R,b,, point b, branch A, had two multipliers at - 1). There is of course a replicated (not symmetry-breaking) bifurcation at R,Bu,= R,b,, when branch B, colludes with branch A, in precisely the same fashion that B, collides with A,.
4. DISCUSSlON
We presented a computational approach that allows the prediction of secondary bifurcations (without actually calculating these bifurcations) off of replications of a “primary” branch of stationary solutions for a set of coupled reaction diffusion equations. This is accomplished through the computation of the Floquet multipliers of the limit cycle of a system of equations like eqs (2), corresponding to stationary solutions of the reaction-diffusion system with periodic boundary conditions. This computation is performed only around the primary branch, but it allows directly the prediction of the bifurcations off of its replications over a much wider interval of system sizes. As we discussed in Section 2, our calculations can also be used to predict the corresponding structure of the bifurcation diagram for no-flux boundary conditions, and, more generally, for boundary conditions that can be extended so as to be periodic. Obviously, the method can be used for larger sets of coupled reaction-diffusion equations. For example, for three coupled reactiondiffusion equations we would have six ODES in set (2), and therefore six Floquet multipliers to monitor on every primary branch of stationary solutions. In nonisothermal cases, since the equation for temperature is formally similar to a reaction-diffusion equation, the same technique can again be applied. The method is currently applicable in one spatial dimension. It can obviously be extended to systems in two dimensions with periodic boundary conditions in one of these two dimensions. While this approach will predict the various bifurcation points on replicated branches, like the birth of the branch C, A, and A, at R2b3 ofTof A,, as well as its connecting collision with A, at R,t, and its end on A, at R,b, (all by just looking at the A, branch), it will not tell us that there will indeed exist one single branch connecting all these three bifurcation points. We are currently working, in collaboration with Professor J. C. Eilbeck, on a two-parameter study of a reaction-diffusion system. The purpose is to examine codimension-two steady-state bifurcations in this setting. Another fruitful direction for research would be to combine this approach with results from studies of bifurcations with symmetry-like the local struc-
1900
IOANNIS
G.
KEVREKIDIS
and HARRY S. BROWN
Bifurcoton
parameter
0
Fig. 5. Bifurcation diagram of the Kuramotc+Sivashinsky PDE with periodic boundary conditions and its replicated bifurcations (the parameter dl is proportional to the square of the system length). Notice the striking similarity of the secondary bifurcation patterns (e.g. the same number of symmetry-breaking bifurcations on the second replica of the “primary” branch etc.).
ture in the neighborhood of a point like RJtZ, analyzed in the work of Golubitsky and Krupa (Krupa, 1988). It is obvious that, through this method, a large amount of information on the location and nature of many steady-state bifurcations over a large interval of parameter values, is “compressed” in the study of the Floquet multipliers on a couple of “primary” branches. It is therefore straightforward to compare the steady-state bifurcation behavior of different models-as the system size varies-by simply comparing the scenario of movement of the Floquet multipliers on these primary branches. One of the most startling observations of this work came from comparing Fig. 3 of the Brussellator with Fig. 5, from an analysis-in one spatial dimension-of the steady solutions of the Kuramoto-Sivashinsky (K-S) equation as the system size varies [see Scovel et al. (1988) and Kevrekidis et al. (1988)]: u, + 4%X,
+ G(u,, + ruu, = 0.
This is an amplitude equation describing thin film and interfacial instabilities as well as wavy @ame fronts. The structure of the two diagrams is almost identical: there is no branch comparable to B, here, but we again have the equivalents of A ,, A2 etc. “A,” again go again 4 times hits “A ?“. Along “A I” the multipliers around the unit circle (predicting exactly four extra bifurcations off of branch “A,“!). The branch “C,“, connecting “A,” and “A, ” is going through strikingly similar bifurcations interacting with almost the same branches as in the reaction diffusion system. The two systems (the Brussellator and the K-S) are not related in any apparent fashion: they model different physical systems, and their nonlinearities are completely different. The ‘most significant common feature is the periodic boundary conditions--lt appears that many
features of pattern formation and the interplay between intrinsic pattern size and characteristic system size is indeed quite robust to system changes. We believe that this method will prove a valuable computational tool in predicting the detailed structure of the steady-state bifurcation diagrams of any coupled reaction diffusion system .as the system size varies. It is also worth mentioning that, if we nondimensionaiize equations like eqs (1) with the system size, this size appears in the same dimensionless group as the diffusion coefficients. It is therefore possible that, under some conditions, while this method is principally applicable in predicting bifurcation structures with respect to characteristic system size, it may also be used to predict bifurcations with respect to other parameters, entering the same dimensionless groups the system size affects. In our opinion, however, this approach should be particularly useful in comparing the complex bifurcation structures obtained by different models and the pattern formation mechanisms across different physical systems [e.g. Sirovich and Newton (1986) and Zufiria (198711. Acknowledgements-This work was_ partially supported through NSF grants Nos CBT-8707090 and EET-8717787. The inspiration and assistance of Dr James C. Scovel of Los Alamos National Laboratory, as well as discussions with Professor J. C. Eilbeck and his group are gratefully acknowledged. NOTATION A,
bj
B
concentration of species A and B (constants) a point on the first branch where the principal Floquet multipliers are at the square root of unity (= - 1) (the subscript indicates that this occurred during the jth time that the multipliers traverse the unit circle)
1901
Predicting pattern formation in coupled reactionditLus.ion systems
Dw Du ‘1
i
L qj
&I tj
u
” X
Greek u
d 0
diffusion coefficients for species u and D, respectively a point on the first branch where the principal Floquet multipliers are at the fifth root of unity (subscript is the same as for bj) square root of - I length of reactor, typical pattern size a point on the first branch where the principal Floquet multipliers are at the fourth root of unity (= f i) (subscript is the same as for bj) kth replica of a primary solution [ R,u(x) = u(kx)] a point on the first branch where the principal Floquet multipliers are at the cube root of unity [ =exp (+2ni/3)] (subscript is the same as for bj) concentration of a reacting chemical species model (also interfacial in the Brussellator film height in the Kuramoto-Sivashinsky equation) concentration of a reacting chemical species spatial variable
letters
bifurcation parameter for the Kuramoto-Sivashinsky equation Floquet multiplier argument of a principal Floquet muttiplier REFERENCES
Babloyantz, A. and Hiemaux differentiation and generation
J.. 1975, Models for cell of polarity in diffusion-
governed 637-657.
morphogenetic
fields. Bull. math.
Biol.
37,
Doedel, E. J., 198t, AUTO: a program for the bifurcation analysis of autonomous systems. Gong. Num. 30,265-284; also, 1986, AUTO 86 User Manual. Caltech, Pasadena, CA. Duncan, K. and Eilbeck, J. C., 1987, Numerical studies of symmetry-breaking bifurcations in reactiondiffusion systems, in Prbceedings of the international Workshop “Biomathematics and Related Computational Problems”(Edited by L. M. Ricciardi). Reidel, Dordrecht. Herschkowitz-Kaufman, M., 1975, Bifurcation analysis of
‘nonlinear reaction-diffusion solutions and comparison
equations--II. Steady state with numerical simulations.
Bull. math. Bid. 37, 589-636. Iooss, G. and Joseph, D. D., 1980, Elementary Stability and Bifurcation Theory. Springer, New York. Kevrekidis, 1. G., Nicolaenko, B. and Scovel, J. C., 1988, Back in the saddle again: a computer-assisted study of the Kuramoto-Sivashinsky equation. SIAM J. appl. Math. (in press). Krupa, M., 1988, Bifurcations of critical group orbits. Ph.D.
dissertation, University of Houston. Kubicek, M. and Marek, M., 1983, Computational Methods in Bifurcation Theory and Dissipative Structures. Springer
Series in Computational Physics, Springer, New York. Kubicek, M., Ryzler, V. and Marek, M., 1978, Spatial structures in a reactiondiffusion systemdetaited analysis of the “Brussellator”. Eiophys. Chem. 8, 253-246. Scovel, C.. Kevrekidis, I. G. and Nicotaenko, B., 1988, Scaling laws and the prediction of bifurcations in systems modeling pattern formation. Phys. Letr. A 130, 73-80. Sirovich, L. and Newton, P. K.. 1986, Periodic solutions of the Ginzburg-Landau equation. Physica 21D, 115-125. Zufiria, J. A., 1987,’ Symmetry breaking in periodic and solitary gravity-capillary waves on water of finite depth.
JFM 184, 183-206.
1989.
am-2soV/nv s3.00+0.0 0 IV89 Pergamon Press plc
PREDICTING PATTERN FORMATION IN COUPLED REACTION-DIFFUSION SYSTEMS IOANNIS Department of Chemical (Received
G. KEVREKIDIS and HARRY S. BROWN Engineering, Princeton University, Princeton, NJ 08544, U.S.A. 10 November
1988; accepted 2 February
1989)
Abstract-Coupled reaction4iffusion equations are known to exhibit a wealth of multiple coexisting stationary solution patterns as the characteristic length of the system grows. We describe and implement a technique which allows us. by studymg only stationary solution branches at small system lengths, to predict the complex structure of steady-state bifurcations occurring at large system lengths without actually computing this structure. This technique is applicable to arbitrary isothermal or nonisothermal systems of coupled reaction-StTusion equations with periodic or no-flux boundary conditions. To illustrate it we use a standard literature example (the Brussellator), where. by computing only the first two solution branches, we accurately predict the steady-state bifurcation structure reported up to much larger system lengths. This technique also provides a compact way ofdescribing and comparing stationary pattern formation in a large class of systems, extending beyond coupled reaction-diffusion equations. We demonstrate this by comparing stationary pattern formation for our test problem (the Brussellator), with the formation of complex surface wave patterns in thin liquid film flow.
INTRODUCTION A number of nonlinear reaction-diffusion systems can have more than one steady-state solution; some of these solutions may be spatially nonuniform. This property can be used in modeling a variety of pattern fcirmation phenomena in biochemical and chemical contexts, from morphogenesis and cell differentiation in embryonic development to surface concentration patterns in catalytic kinetics. This property is also closely related to the emergence of spatial nonuniformities (asymmetric steady states) in systems of interacting catalyst particles. Many possible models of interest exist, giving rise to different reaction rate functions and different bifurcation diagrams as the important parameters (the kinetic rate constants, the diffusion coefficients and the characteristic system size) vary. In particular, it is well known by numerical studies on a number of specific models [from, for example, Babloyantz Hiernaux (1975), and Herschkowitz-Kaufman (1975) and Kubicek et aZ. (1978) to, more recently, Duncan and Eilbeck (1987)] that the steady-state bifurcation diagrams for isothermal, coupled reaction-diffusion systems with noflux boundary conditions often possess a bewildering complexity, which increases as the characteristic system size increases. This complexity results from the interplay between an intrinsic pattern size (influenced by the reaction rate form, the values of the kinetic constants and the diffusion coefficients) and the characteristic system size, imposed by the boundary conditions. While the results of this interplay are indeed model-dependent, one would expect the general mechanisms of pattern formation to be reasonably robust to changes in the models. In this paper we describe and implement a technique which allows us, by studying stationary solution 1893
branches at only small often complex structure much larger system sizes this structure. This can
tational behavior.
system sizes, to predict the of bifurcations occurring at without actually computing lead to significant compu-
savings in accurately describing the system The technique is applicable to arbitrary
nonisothermal systems of coupled equations with no-flux or periodic Its general applicability, boundary conditions. coupled with the ease of implementation and the resulting computational savings, make it a useful tool in predicting and understanding the detailed steadystate bifurcation behavior of any specific model. We believe, however, that its usefulness goes beyond mere computational savings in a particular model study. Since an entire complex pattern-forming bifurcation process can be reduced to the study of a few branches, isothermal or reaction4ffusion
this technique provides a compact way of describing and comparing pattern formation across several models. This is an efficient approach to quickly and accurately assessing the effect of changes in the model to the predicted steady-state behavior over large intervals of parameter values without actually performing detailed computations. It may also prove helpful in understanding the degree to which this apparently complicated process of pattern formation is robust and model-independent. This paper is structured as follows. In Section 2 we describe the conceptual underpinnings of the method, and present the computational steps involved in its implementation. In Section 3 we ilhistrate this method by applying it to a popular literature example: the “Brussellator”. The main reason for choosing this as a “typical” example to illustrate the applicability of the method is that many chemical engineers are probably familiar with its behavior, since it has been the object
IOANNIS
1894
G.
and
KEVREKIDIS
of intense numerical calculations in several scientific papers, and it has also appeared as an illustrative example in monographs [e.g. Kubicek and Marek (1983)]. A different example (i.e. one with a rational polynomial rate expression, like the ones resulting from enzyme or Langmuir-Hinshelwood kinetics) could serve equally well to illustrate the method. in this example, we are abte to accurately predict the secondary bifurcation structure observed and reported in the literature for large characteristic system sizes, by studying only the first two steady-state branches at small system sizes. Finally, in Section 4, we use the method to discuss the robustness of the pattern formation process in different models. This is accomplished by comparing the common features of the bifurcation diagram for our test example-the Brussellator-with the formation of complex surface wave patterns in thin liquid film flow. 2. THE METHOD Mnf ivation
To motivate our computational procedure we consider for simplicity a set of two coupled reaction-diffusion equations, occurring in a one-dimensional medium. The steady-state equations are
D”D,, +$I(% 0) = 0
(1)
where u and v represent concentrations of the two species, D, and D, are diffusion coefficients, and f (u, u) and g(u, o) are the reaction rate forms. Let us for the moment consider periodic boundary conditions [u(O) =u(L), u(O)=v(L)] and suppose that system (1) has, for some set of parameter values and some system size L, a certain “primary” solution u(x). v(x) (0 G x < L). It is then obvious that if we double the system size and
31
I 02
S. BROWN
extend the steady-state pattern two-fold [u(x + I_) =u(x), V(X+ L)= v(x)] we obtain a soIution of the same equations with the same boundary conditions for the new system size 2L. This holds in general for the n-fold extension of this first solution, which will be a solution of the same equations with the same type of boundary conditions when the system size in nL. This by no means implies that the nth replication of the “primary” solution in the only solution of eqs (1) when the system size in nL. Neither does it imply that the nth replication of the “primary” solution has the same stability characteristics for system size nL that the “primary” solution has at system size L. This is easily seen upon inspection of the bifurcation diagram reported in the literature for the “Brussellator’‘-see Fig. 1. At L= o! the “primary” branch A, of spatially nonuniform stationary solutions bifurcates from a spatially flat solution. This implies then that at L = 20~ the two-fold replica A, of A, bifurcates from the spatially flat solution, that at L=3a the three-fold replica A, of A, bifurcates from the spatially flat solution etc. We observe, however, that, while soon after its birth, the “primary” branch A, is the only spatially nonuniform stationary solution, its two-fold replica A, over most ofits interval of existence coexists with several other spatially nonuniform stationary solutions. Another important observation is the following: the only two steady-state bifurcations observed on the diagram for the “primary” branch AI are its birth (at L=a) and its “death” when it collides with the A, branch at L=c. By the simple “replication” arguments, this means that the two-fold replica of A, will be “born” at L = 21x,and it will “die” colliding with the four-fold replica of A, at L = 2c, which is indeed o,bserved. This was observed, for example, by Kubicek et al. (1978), and is related to what they termed “composition of solutions”. We notice, however, that
I 04
I 06 Length
Fig. 1. Steady-state periodic boundary
HARRY
I 08
I I
L
bifurcation diagram with respect to the system size L for the “Brussellator” with conditions up to L= 1 for A=2, B=4.6, D,= 1.6 x 10-j, D,=8 x 10m3.Only solutions symmetric about x = L/2 are shown.
Predicting pattern formation in coupled reaction-diffusion the two-fold
replica
also connects
with several
other
branches at several (four, to be exact) “secondary” bifurcation points, for which there is no counterpart on the “primary” branch A,. Furthermore, it is obvious even from the part of the diagram shown, that the four-fold replica A, of the A, branch will have even more secondary bifurcations, giving birth to a large number of so-called “mixed-mode” branches. Our goal in this paper is to present a simple technique for accurately predicting, as the system size L grows, all these secondary symmetry-breaking bifurcations off of all the n-fold replications of a given “primary” branch by only studying the “primary” solution branch, which in itself does not exhibit them. This technique is generally applicable in systems modeling pattern formation with periodic boundary conditions [see, for example, Scovel et al. (1988)]. What makes it particularly interesting in the case of reaction-diffusion systems, where no-flux boundary conditions are often the reasonable choice, is the following simple observation. Consider a solution u(x). u(x) of eqs (1) with no-flux (Neumann) boundary conditions at x = 0 and x = L/2: u,=u,=O If we extend
at x=0
this solution
and x=-.
between
L 2
L/2 and L so that
we obtain a solution of system (1) with periodic boundary conditions at x=0 and x= L. This means that u[Z stationary solutions of system (1) with Neumann boundary conditions for system size L/2 are also stationary solutions of eqs (1) with periodic boundary conditions for system size L. The converse is, of course, not necessarily true, i.e. there may exist solutions of eqs (1) with periodic boundary conditions that are not necessarily also solutions of the Neumann problem at half the system size. The method we describe can, of course, be used to study and predict the bifurcations of these latter “periodic” but LLnonNeumann” solutions when they exist. This simple property of the coupled reaction-diffusion equations, however, allows us to use the method, which hinges on periodic boundary conditions, also in the case of Neumann boundary conditions for reaction+liffusion systems. Indeed, all the stationary solution branches shown in Fig. 1, which have been computed with periodic boundary conditions for various system sizes, have Neumann counterparts at half the system size shown.
L = 2a, and collides
systems
1895
with A, at L = 2~). It is straightfor-
ward to predict that these bifurcations will occur, and the parameter values where they will occur, by simply knowing their “counterparts” for the primary branch A,. As a matter of fact, it is obvious that one can predict the “replications” of these bifurcations for the three-fold replication of A,, the branch A,. This branch will be born at L= 3a, and will collide with A, at L=3c, beyond the values of L shown. Similar results hold for all other n-fold replications of A,. There exist, however, four. parameter values at which A, exhibits steady-state bifurcations for which A, has no counterpart. We are interested in the possibility of the predicting these “secondary” bifurcations by only studying the primary branch. When, For practical computational purposes, we discretize a problem like eqs (I), say through finitedifference or pseudo-spectral techniques, the approximation of its steady-state solutions reduces to the solution of severai (say N) nonlinear simultaneous algebraic equations. These solutions can be followed using standard continuation techniques [e.g. AUTO (Doedel 1981, 1986)] as the system parameters (here the system size L) vary. To accurately calculate steadystate bifurcation points one searches for parameter values for which an eigenvalue of an N x N square matrix (the linearization of the N nonlinear algebraic equations around their solution) crosses zero. Such points are then further examined to determine whether a steady-state bifurcation does indeed occur there. “Secondary” bifurcations, occurring off of replications of a primary branch, but having no counterpart on the primary branch itself, become easier to understand and predict in the following simple fashion. Consider a steady-state solution on the primary branch for some parameter value (say L ~0.2) [Fig. 2(a)) (remember that “half” of this solution satisfies the Neumann boundary conditions for half the system size). Then at L=O.4 the two-fold replication of this solution [Fig. 2(c)] will be again a solution, lying on the branch A,. A sready state of problem (1) with periodic boundary conditions and system size L can be considered as a limit cycle (with period L) of the following set of ordinary differential equations (ODES) (in the spatial variable x):
du, dx
*
4 dx
4
The technique and its implementation In Fig. 1 the two-fold replication A, of the primary branch AL does possess the replications of the steadystate bifurcations that A L exhibits (i.e. A, is “born” at
where (u,, ul, uj, u,)~=(u, u,, u, u,)~. A projection of this limit cycle on the (u,, ul) plane (i.e. u and u,) is
IOANNIS G. KEVREKIDIS and HARRY S. BROWN
1896
2fb) ;
El 0
-0
4
”
branch is now apparent. Suppose that, for some value of L, the steady solution of eqs (1) lying on A, haswhen viewed as a limit cycle of eqs (2 ja Floquet multiplier at the nth root of unity (say for example at the square root of unity, - 1, or at a fourth root of unity, &-i). Then, at n times that system size, n_L, the nth replication of the primary branch A, will have-when viewed as a limit cycle of (2) -a Floquet multiplier exactly at 1, and that is a potential bifurcation point of limit cycles for (2), or steady states of (Itpotential because one still has to verify that a bifurcation will indeed occur there. Consider for example all j&r ‘*secondary” bifurcation points observed on branch A,. They occur at precisely twice the system sizes for which the corresponding solution on A, had-when viewed as a limit cycle of (2kFloquet multipliers at
-1. The implementation
Fig. 2. Spatial profile of a solu&ion at L = 0.2 [Fig. 2(a)] and its representation as a limit cycle in x with period L=O.Z projected on the u-u, plane [Fig. 2(b)]. Figure 2(c) shows a two-fold replica of the same solution at LzO.4. Species u: solid lines; species ~1:dotted lines. plotted in Fig. 2(b). Moving along the solution curve as x varies between 0 and L=O.2 in Fig. 2(a) is equivalent to moving the phase point [u(x), u,(x)] once around the limit cycle in Fig. 2(b). It is now obvious that the two-fold replication of this solution when L=O.4 is exnctly the same limit cycle, only the phase point goes around it twice. This is the crucial observation on which the technique is based: steady states of system (1) with periodic boundary conditions at system size L are Ziinit cycles of the. set of ODES (2) steady-state bifurcations of with period I.. Therefore, problem (I) to other steady-states can be viewed as bifurcations of limit cycles of problem (2) to other limit cycles. A necessary condition for bifurcations of limit cycle branches to other limit cycles (with nearby periods) to occur is that an eigenvalue of the state transition matrix (the so-called Floquet or characteristic multipliers) must cross unity [see, for example, Iooss and Joseph (1980)]. This means that, at parameter values L for which the branch A, does exhibit a steady-state bifurcation (while branch A, dots not), rotating around a limit cycle twice (solution on branch AZ) yields a Floquet multipher at unity, while rotating around the limit cycle only once (solution on branch A, at half the system size) does not. There exists, however, a very simple relation between the Floquet multipliers obtained when we go around a limit cycle once, and those we obtain when we go around the same limit cycle twice. Namely, going around twice will syuarr the Floquet multipliers (similarly, going around 3 times will raise them to the third power, and going around n times will raise them to the nth power). The technique for predicting the secondary bifurcations by only using information on the primary
of this technique is therefore very simple, and can be appended with minimal effort to any continuation code currently used: after a solution on the primary branch (with periodic boundary conditions) has been computed-through any discretization method-we simply integrate the variational equations of system (2) around this solution once. We will obtain a (2n) x (2n) state transition matrix (where n is the number of coupled reaction-diffusion equations, in our test case n = 2). Its 2n eigenvalues (in our test case 4) are the Floquet multipliers; one of them has to always be at unity [see, for example, Iooss and Joseph (1980)]. When one of the remaining multipliers goes through an nth root of unity in the complex plane, then-at n times that system size the nth replication of the primary branch will potentially have a “secondary” bifurcation. These bifurcations give rise to the so-called “mixed-mode” branches.
3. THE EXAMPLE The equations We will illustrate the method by applying it to a system of two coupled reaction-diffusion equations, the “Brussellator”: U,,+;CA+U%-_(B+ U
l)u]=O
I u,+~(Bu-u%)=O.
In particular, we will examine the steady-state bifurcation diagram obtained by numerical computations for A=2, B=4.6, D,=1.6x lo-’ and D,=8x 10e3 with periodic boundary conditions. This diagram (the one presented in Fig. I) may be more familiar to the
reader in the form it (1978) or Kubicek Fig. 3(c)]. The reader after our discussion
was presented by Kubicek et al. and Marek (1983) [see also should
however keep in mind, 2, that the branches reported here at a given system size have no-flux counterparts, reported by Kubicek et al., at half that system size. To further facilitate the comparison, we
in Section
1897
Predicting pattern formation in coupted reaction-diffusion systems
a
2a
I 02
Length
0
I 01
L
I
I
I
0.2
0.4
0.6
Length
13
I 0.4
03
I a*
1 I
L
Fig. 3. (a) Two “primary” branches. At specially labeled points the principal Floquet multipliers are located at nth roots of unity on the unit circle (n = 1 is symbolized by a, n = 2 by b, n = 3 by t, n = 4 by q etc). The same symbols preceded by a B denote such occurrences on the second branch B,.(b) Secondary bifurcations on replicated branches that can be predicted based on Floquet multiplier information on the primary branches. (c) Same bifurcation diagram in the format reported by Kubicek et al. (1978): instead of the L* norm of the solution,
u (0) is plotted
have only included in our diagram solutions of the problem with periodic boundary conditions that do possess no-flux counterparts. As we will discuss, there exist bifurcations to periodic solutions that do not possess no-flux counterparts. We will perform Floquet “primary” calculations on the two multiplier
vs t.
branches, A, and E,, and from that information we will predict all other secondary bifurcations in the diagram and several others that will be observed at values of L far beyond the ones reported here. There exists only one more technicality we must discuss before examining the bifurcation diagram in
IOANNIS G. KEVREKIDIS and HARRY S. BROWN
1898
Fig. 3. In our test problem, and every problem involving two coupled reaction diffusion equations, we will calculate four Floquct multipliers [since we have four ODES in system (2)]. We observe computationallyand it is possible to show theoretically-that, for the stationary solutions we discuss in Fig. 3, two of these four multipliers (let us call them E., and ia) will always be fixed at unity, and the remaining two, & and E.,, which we will call the principal Floquet multiphers (PFMs) will vary as the system size (and the entire solution) varies, while their product always remains at 1. This means that they may either be both positive or both negative such that A,= I/>.,, or that they are complex conjugate with modulus 1 (i.e. they will lie on the unit circle in the complex plane). Since this fact is not germane to the technique, the reader may simply consider it as a computational observation. Finally, in what follows, we have tried IO use a convenient mnemonic notation for the parameter values at which the PFMs lie at roots of unity, and for the predicted bifurcations. Points on a primary branch, for which the PFMs lie at an nth root of unity are marked b, (for square roots), ti (for cube roots), qi (for fourth rootstwhere i is an index to discriminate between them. The corresponding predicted bifurcation points on replications of the primary branch are marked with rhe same symbol preceded by an R, (R stands for replication and k denotes a k-fold replication of the primary branch). For example, b, denotes the first time (L) on the primary branch A 1 for which we observe a PFM at the square root (b) of unity. This predicts a bifurcation off of the two-fold replication (R,) of A, at twice (I) the system size (RA). The primary
branch
A, (Spatially At L= LX = 0.159689 the first nontrivial varying) branch bifurcates off of the flat solution U(Y) = A, u(x) = B/A. Figure 3(a) shows an enlargcmcnt of the bifurcation diagram of Fig. 1 between L-=0.1 and L=O.5, obtained through a spectral discretixation and a modification of the bifurcation code AUTO written by Doedel (1981). Figure 3(c) is the same bifurcation diagram reported [as Kubicek et al. (1978) did] by plotting u(O) vs the bifurcation parameter L. At L = o! =0.159689 the two principal Floquet multipliers lie at i 3,4 = - 0.6428 f 0.76601, or at an angle of approximately 130” (see Fig. 4, and continue referring to it). As L grows, these two multipliers move along the unit circle in the complex plane, and they go through the fifth root of unity at L = 0.16527 (point f,), which will then mean that the fifth replica of this branch, marked A, in the diagram, will have a symmetry-breaking bifurcation at L = 0.826 % 5 x 0.16527 (R,f,). Indeed, that is a bifurcation point of the fifth replica A, of A,. From now on, we will only consider bifurcations of replicated solutions for k,<4 (up to fourth roots of unity), since higher replication branches lie mostly outside the diagram of Fig. 1. Subsequently the PFMs come to a double - 1 at L =0.17404 (point b,)-see again Fig. 4-which pre-
Fig. 4. Structure and direction of motion of the principal Floquet multipliers as the parameter L changes at the points labelled
on the first primary
branch
A, of Fig. 3 (a).
diets a symmetry-breaking bifurcation for branch A, at L=O.3481 (point R,b, where branches A, and A, collide-and further replica at R,b, where branch A, runs into branch A, in the same fashion that branch A, runs into branch A, at R,b,). After that, the PFMs become real and split on the real axis, one inside and one outside the unit circle. They return to a double - 1 at L=O.18363 (point b2), which predicts a symmetrybreaking bifurcation for branch A, at L = 2 x 0.18363 = 0.36726 (point R,b2, bifurcation of branch A, into a “new” secondary branch which we will call I?-with a They then go further replica at R,b,, L=O.7344). through the third root of unity [&, = - 0.5 f 0.5iJ3, or exp ( + i2z/3), 0 = 120” (see Fig. 4)] at L = 0.19423 (point t2), with a predicted symmetry breaking bifurcation on branch A,
Predicting pattern formation in coupled reaction+Iiffusion which again there is no bifurcation in the subspace of functions with the property that u(L/2 - x) = u(L/2 +x). At that point we do, however, have a bifurcation to an asymmetric (around L/2) periodic branch, which has no Neumann counterpart at half the system size. They now go again around the unit circle. They pass through the fourth root of unity at point q3 (L = 0.3409, symmetry-breaking bifurcation on the A, branch at L = 4 x 0.3409 = 1.3634, point R,q,, off scale in the diagram). They pass through a cube root of unity (point t,, L = 0.34 1146, symmetry-breaking bifurcation at L= 3 x 0.341146= 1.023438 off the A, branch, point R,t,, again off scale). They hit the double-negative one at L=O.341382 (point b3), with a symmetry-breaking bifurcation for the A, branch at L =0.68276 (point R,b,). They split on the real axis (both negative, one inside and one outside the unit circle) and they return to a bifurdouble ~ 1 at L = 0.35430 (point b,, predicted cation at L= 0.7086 off of branch A,, point R,b,). They go through the cube root of unity at L=O.35359 off of the A, branch at (point t,, predicted bifurcation L = 1.0608, off scale in the Fig. 4) and through the fourth root of unity at L~O.35259 (point q4, predicted bifurcation off of A, at L=1.4104, off scale). They finally come to a quadruple 1 at point u,-R,b,, where the A, branch merges with A,.
The primary
branch 9, This branch also bifurcates from the flat solution, at L=O.442243 (point p) and, like the A, branch, it eventually connects with the A, branch. The principal Floquet multipliers on this branch start at point fi at as is predicted 0~83.0”) E.,,,=0.1216f0.9926i(angle from linear analysis around the flat steady state. They then start moving around the unit circle (as L decreases) and they reach 90” (a fourth root of unity) at L bifurcation =0.441443 (point Bq,, with a replicated of the replicated branch B, at L= 1.7658, off scale). They reach 120” at L=O.43778 (point Bt,, with a symmetry-breaking bifurcation of the replicated branch B, at L= 1.31334, again off scale). They come to a double ~ 1 at L=O.431876 (point Bb,), which predicts a symmetry-breaking bifurcation off of the replicated branch B, at L = 0.863753 (point R,Bb,, birth of branch Dl). They split on the negative real axis (one inside and one outside the unit circle) and return to - 1 at L=O.42645 (point Bb,, with replicated symmetry-breaking bifurcation at L=O.85290, point R,Bb,, death of branch DI). They traverse the unit circle, going through the third and fourth roots of symmetryunity (points Bt2 and Bq2, with replicated breaking bifurcations off of branches B, and B,, respectively, off scale in the graph). They come to a quadruple multiplier at 1 at point Bu, (L=O.4089) split on the real axis and come back again to a quadruple multiplier at 1 at point Bu, (L=O.37624). As was the case with the A, branch, no bifurcation is observed for one of these quadruple 1 points, while an asymmerric branch, not depicted in the diagram, is
systems
1899
expected at the other quadruple 1 point. The PFMs subsequently start moving on the unit circle, and now a ‘&new” phenomenon is observed: instead of moving monotonically, they reach a “maximum” angle at approximately 24” (LzO.37157, point S) and then return towards a quadruple 1, where they arrive at L =0.367108, point Bu,, which has been predicted as a symmetry-breaking bifurcation off of branch A, (point L=2xO.l8355) (remember, at L=O.18355, R,b,, point b, branch A, had two multipliers at - 1). There is of course a replicated (not symmetry-breaking) bifurcation at R,Bu,= R,b,, when branch B, colludes with branch A, in precisely the same fashion that B, collides with A,.
4. DISCUSSlON
We presented a computational approach that allows the prediction of secondary bifurcations (without actually calculating these bifurcations) off of replications of a “primary” branch of stationary solutions for a set of coupled reaction diffusion equations. This is accomplished through the computation of the Floquet multipliers of the limit cycle of a system of equations like eqs (2), corresponding to stationary solutions of the reaction-diffusion system with periodic boundary conditions. This computation is performed only around the primary branch, but it allows directly the prediction of the bifurcations off of its replications over a much wider interval of system sizes. As we discussed in Section 2, our calculations can also be used to predict the corresponding structure of the bifurcation diagram for no-flux boundary conditions, and, more generally, for boundary conditions that can be extended so as to be periodic. Obviously, the method can be used for larger sets of coupled reaction-diffusion equations. For example, for three coupled reactiondiffusion equations we would have six ODES in set (2), and therefore six Floquet multipliers to monitor on every primary branch of stationary solutions. In nonisothermal cases, since the equation for temperature is formally similar to a reaction-diffusion equation, the same technique can again be applied. The method is currently applicable in one spatial dimension. It can obviously be extended to systems in two dimensions with periodic boundary conditions in one of these two dimensions. While this approach will predict the various bifurcation points on replicated branches, like the birth of the branch C, A, and A, at R2b3 ofTof A,, as well as its connecting collision with A, at R,t, and its end on A, at R,b, (all by just looking at the A, branch), it will not tell us that there will indeed exist one single branch connecting all these three bifurcation points. We are currently working, in collaboration with Professor J. C. Eilbeck, on a two-parameter study of a reaction-diffusion system. The purpose is to examine codimension-two steady-state bifurcations in this setting. Another fruitful direction for research would be to combine this approach with results from studies of bifurcations with symmetry-like the local struc-
1900
IOANNIS
G.
KEVREKIDIS
and HARRY S. BROWN
Bifurcoton
parameter
0
Fig. 5. Bifurcation diagram of the Kuramotc+Sivashinsky PDE with periodic boundary conditions and its replicated bifurcations (the parameter dl is proportional to the square of the system length). Notice the striking similarity of the secondary bifurcation patterns (e.g. the same number of symmetry-breaking bifurcations on the second replica of the “primary” branch etc.).
ture in the neighborhood of a point like RJtZ, analyzed in the work of Golubitsky and Krupa (Krupa, 1988). It is obvious that, through this method, a large amount of information on the location and nature of many steady-state bifurcations over a large interval of parameter values, is “compressed” in the study of the Floquet multipliers on a couple of “primary” branches. It is therefore straightforward to compare the steady-state bifurcation behavior of different models-as the system size varies-by simply comparing the scenario of movement of the Floquet multipliers on these primary branches. One of the most startling observations of this work came from comparing Fig. 3 of the Brussellator with Fig. 5, from an analysis-in one spatial dimension-of the steady solutions of the Kuramoto-Sivashinsky (K-S) equation as the system size varies [see Scovel et al. (1988) and Kevrekidis et al. (1988)]: u, + 4%X,
+ G(u,, + ruu, = 0.
This is an amplitude equation describing thin film and interfacial instabilities as well as wavy @ame fronts. The structure of the two diagrams is almost identical: there is no branch comparable to B, here, but we again have the equivalents of A ,, A2 etc. “A,” again go again 4 times hits “A ?“. Along “A I” the multipliers around the unit circle (predicting exactly four extra bifurcations off of branch “A,“!). The branch “C,“, connecting “A,” and “A, ” is going through strikingly similar bifurcations interacting with almost the same branches as in the reaction diffusion system. The two systems (the Brussellator and the K-S) are not related in any apparent fashion: they model different physical systems, and their nonlinearities are completely different. The ‘most significant common feature is the periodic boundary conditions--lt appears that many
features of pattern formation and the interplay between intrinsic pattern size and characteristic system size is indeed quite robust to system changes. We believe that this method will prove a valuable computational tool in predicting the detailed structure of the steady-state bifurcation diagrams of any coupled reaction diffusion system .as the system size varies. It is also worth mentioning that, if we nondimensionaiize equations like eqs (1) with the system size, this size appears in the same dimensionless group as the diffusion coefficients. It is therefore possible that, under some conditions, while this method is principally applicable in predicting bifurcation structures with respect to characteristic system size, it may also be used to predict bifurcations with respect to other parameters, entering the same dimensionless groups the system size affects. In our opinion, however, this approach should be particularly useful in comparing the complex bifurcation structures obtained by different models and the pattern formation mechanisms across different physical systems [e.g. Sirovich and Newton (1986) and Zufiria (198711. Acknowledgements-This work was_ partially supported through NSF grants Nos CBT-8707090 and EET-8717787. The inspiration and assistance of Dr James C. Scovel of Los Alamos National Laboratory, as well as discussions with Professor J. C. Eilbeck and his group are gratefully acknowledged. NOTATION A,
bj
B
concentration of species A and B (constants) a point on the first branch where the principal Floquet multipliers are at the square root of unity (= - 1) (the subscript indicates that this occurred during the jth time that the multipliers traverse the unit circle)
1901
Predicting pattern formation in coupled reactionditLus.ion systems
Dw Du ‘1
i
L qj
&I tj
u
” X
Greek u
d 0
diffusion coefficients for species u and D, respectively a point on the first branch where the principal Floquet multipliers are at the fifth root of unity (subscript is the same as for bj) square root of - I length of reactor, typical pattern size a point on the first branch where the principal Floquet multipliers are at the fourth root of unity (= f i) (subscript is the same as for bj) kth replica of a primary solution [ R,u(x) = u(kx)] a point on the first branch where the principal Floquet multipliers are at the cube root of unity [ =exp (+2ni/3)] (subscript is the same as for bj) concentration of a reacting chemical species model (also interfacial in the Brussellator film height in the Kuramoto-Sivashinsky equation) concentration of a reacting chemical species spatial variable
letters
bifurcation parameter for the Kuramoto-Sivashinsky equation Floquet multiplier argument of a principal Floquet muttiplier REFERENCES
Babloyantz, A. and Hiemaux differentiation and generation
J.. 1975, Models for cell of polarity in diffusion-
governed 637-657.
morphogenetic
fields. Bull. math.
Biol.
37,
Doedel, E. J., 198t, AUTO: a program for the bifurcation analysis of autonomous systems. Gong. Num. 30,265-284; also, 1986, AUTO 86 User Manual. Caltech, Pasadena, CA. Duncan, K. and Eilbeck, J. C., 1987, Numerical studies of symmetry-breaking bifurcations in reactiondiffusion systems, in Prbceedings of the international Workshop “Biomathematics and Related Computational Problems”(Edited by L. M. Ricciardi). Reidel, Dordrecht. Herschkowitz-Kaufman, M., 1975, Bifurcation analysis of
‘nonlinear reaction-diffusion solutions and comparison
equations--II. Steady state with numerical simulations.
Bull. math. Bid. 37, 589-636. Iooss, G. and Joseph, D. D., 1980, Elementary Stability and Bifurcation Theory. Springer, New York. Kevrekidis, 1. G., Nicolaenko, B. and Scovel, J. C., 1988, Back in the saddle again: a computer-assisted study of the Kuramoto-Sivashinsky equation. SIAM J. appl. Math. (in press). Krupa, M., 1988, Bifurcations of critical group orbits. Ph.D.
dissertation, University of Houston. Kubicek, M. and Marek, M., 1983, Computational Methods in Bifurcation Theory and Dissipative Structures. Springer
Series in Computational Physics, Springer, New York. Kubicek, M., Ryzler, V. and Marek, M., 1978, Spatial structures in a reactiondiffusion systemdetaited analysis of the “Brussellator”. Eiophys. Chem. 8, 253-246. Scovel, C.. Kevrekidis, I. G. and Nicotaenko, B., 1988, Scaling laws and the prediction of bifurcations in systems modeling pattern formation. Phys. Letr. A 130, 73-80. Sirovich, L. and Newton, P. K.. 1986, Periodic solutions of the Ginzburg-Landau equation. Physica 21D, 115-125. Zufiria, J. A., 1987,’ Symmetry breaking in periodic and solitary gravity-capillary waves on water of finite depth.
JFM 184, 183-206.