- Email: [email protected]
Stability and pattern formation in a coupled arbitrary order of autocatalysis system
Available online at www.sciencedirect.com
Applied Mathematical Modelling 33 (2009) 884–896 www.elsevier.com/locate/apm
Stability and pattern formation in a coupled arbitrary order of autocatalysis system q Li Zhang *,1, San Yang Liu Department of Applied Mathematics XiDian University, Xi’an, ShannXi 710071, PR China Received 19 June 2007; received in revised form 9 December 2007; accepted 11 December 2007 Available online 23 December 2007
Abstract Spatiotemporal structures arising in two identical cells, each governed by arbitrary order autocatalator kinetics and coupled via the diffusive interchange of a reactant, are discussed. The stability of two homogeneous steady states is obtained with the use of linear stability analysis. By studying the linearized equations, it is found that two steady states, in the uncoupled and coupled system respectively, may give rise to the possibility of bifurcations to spatially nonuniform pattern forms. Further information about Turing bifurcation solutions close to these bifurcation points are obtained by weakly nonlinear theory. It is seen that the coupling leads to bifurcations not present in the uncoupled system which give rise to locally stable nonuniform pattern forms. Finally the stability of the equilibrium points of the amplitude equation is discussed by weakly nonlinear theory, with the bifurcation branches about small coupled system with 0 < a 1 and large coupling for a 1. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Reaction diffusion system; Stability; Bifurcation; Coupling; Pattern
1. Introduction The feature common to nearly all isothermal oscillatory reactions is autocatalysis. Autocatalysis is shared not only by the Belousov–Zhabotinskii reaction, a wide range of halide-based oscillations, and the arsenite plus iodate reaction but by numerous enzyme systems and by chain-branching reactions in the gas phase [1]. All mechanistic schemes include autocatalytic steps (or autocatalytic combinations of elementary steps). However, even the simplified forms of these schemes are still complicated, and there is still need to look at even simpler prototypes in a logically ordered way. Accordingly, we consider throughout steps with the overall stoichiometry A ! B. We may also take this as written to represent a reaction satisfying the rate law = ka. When the same change is catalyzed by a species Y and has a rate satisfying the expression kay, we may q *
1
Foundation items: National Natural Science Foundation of China (60574075). Corresponding author. Tel./fax: +86 02988202861. E-mail address: [email protected] (L. Zhang). Currently working towards the Ph.D. degree in applied mathematics at Xidian University.
0307-904X/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.12.013
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
885
represent this as A þ Y ! B þ Y . Autocatalysis corresponds to catalysis by the product B itself. Two exemplary cases [2,3] span the whole range of behavior likely to be encountered: quadratic autocatalysis rate ¼ k q ab cubic autocatalysis rate ¼ k c ab2 These may be represented by the ‘chemical’ equations A þ B ! 2B;
ðaÞ
A þ 2B ! 3B:
ðbÞ
and Normally these will arise from a subscheme of elementary steps. For instance, the arsenite- iodate reaction produces iodate ions with an experimentally determined rate law of the form 2
þ d½I =dt ¼ k 1 þ k 2 ½I ½I ½IO 3 ½H :
This does not, however, require us to believe that it requires a triple encounter as an elementary step. It is important to note (a) that reversibility makes little difference to behavior until the equilibrium constant ð½B=½AÞeq becomes small and (b) that less simple stoichiometry A þ nB ! ðn þ mÞB can readily be incorporated if required [4,5]. In a series of recent papers [6–8], we have looked at the various phenomena which can arise in a simple prototype chemical reaction due to chemical coupling. Cubic autocatalysis has been studied in many papers, for example, the review lecture of Gary [9]. In particular, Hill et al. [10] has shown the formation of stable patterns in a reaction diffusion system based on the cubic autocatalator A þ 2B ! 3B; B ! C, with the reaction taking place within a closed region, the reactant A being replenished by the slow decay of precursor P via the reaction P ! A. In [11], the spatiotemporal structures that can arise in two identical cells, each governed by cubic autocatalator kinetics and coupled via the diffusive interchange of the autacatalyst B, are discussed. In [12], the initiation and propagation of reaction diffusion travelling waves in two regions, coupled together by a linear diffusive interchange across a semipermeable membrane is considered.In [13], a three variable Turing system with two kinds of cells and a diffusive chemical is studied, and the formation and maintenance of fish skin patterns with multiple pigment cells are considered. In [14], pattern formation in biological systems for which the tissue on which the spatial pattern resides is growing at a rate which is itself regulated by the diffusible chemicals that establish the spatial pattern, is investigated. In the present context, we wish to allow a more general approach in which arbitrary order of autocatalysis can be encompassed, namely A þ mB ! ðm þ 1ÞB;
rate ¼ kabm :
This simple mth order power dependence on the autocatalyst concentration plays the role of the usual exponential dependence of the reaction rate on the temperature rise DT in the systems, which can be represented somewhat schematically as A þ T ! T þ DT ;
rate ¼ kae/DT ;
where / ¼ E=RT 2ref and T ref is some suitable reference temperature at which k is evaluated. In this paper, we consider the prototype chemical reaction scheme based on arbitrary order autocatalator P ! Aðrate k 0 pÞ;
ðdÞ
A þ mB ! ðm þ 1ÞBðrate k 1 abm Þ;
ðeÞ
B ! Cðrate k 2 bÞ;
ðfÞ
where a and b are the concentrations of the reactant A and the autocatalyst B and the k i are rate constants (both reactions are assumed to be isothermal). Here we suppose that the reaction is taking place inside a closed system with reactant A now being replenished by the slow decay of a precursor P, via the simple isothermal step (d), where p is the concentration of reactant P.
886
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
In [15], a model based on an autocatalytic, two-step reaction mechanism including two ionic components (of the same charge) and two non-ionic component, where both reactions are of second order overall, is considered when an electric field is applied to the system. In [16] the effects of electric fields on the reaction fronts that arise in a system governed by an autocatalytic reaction and a complex reaction between the autocatalyst and a complex agent are considered. In [17], Hubbard et al. investigated pattern formation in a coupled system of reaction diffusion equations in two spatial dimensions. These equations arise as a model of isothermal chemical autocatalysis with termination in which the orders of autocatalysis with termination, m and n, respectively, are such that 1 < n < m. In [18], the propagation of planar reaction diffusion waves in the isothermal autocatalytic system A þ mB ! ðm þ 1ÞB; rate kabm , is considered. Attention is paid to the case when the diffusion coefficient of reactant A is much less than that of the autocatalyst; the case when A is completely immobilized is discussed in detail. In the present paper, we consider the reaction scheme (d),(e),(f) taking place in a closed vessel without stirring. This allows for spatial variations in the reactant concentrations throughout the vessel arising via the combined mechanisms of reaction and molecular diffusion. The reaction scheme take place in two identical adjacent cells, coupled by the diffusive interchange of the reactant A. The equations governing the evolution in concentration of the reactants in the vessel are 8 2 oa1 > > ¼ k ooxa21 þ l a1 bm1 þ aða2 a1 Þ; > > ot > > > > 2 ob1 > > < ¼ kd ooxb21 þ a1 bm1 b1 ; ot ð1Þ 2 > oa2 m o a > 2 > ¼ k ox2 þ l a2 b2 þ aða1 a2 Þ; > > ot > > > > > : ob2 ¼ kd o2 b22 þ a2 bm b2 ; 2 ox ot where m > 1; m 2 N , ða1 ; b1 Þ and ða2 ; b2 Þ are the dimensionless concentrations in cells 1 and 2, respectively, a is the dimensionless coupling parameter, and d is the diffusion rate of reactant A and B. We impose zero-flux boundary condition namely oai ¼ 0; ox
obi ¼ 0; ox
on x ¼ 0; 1 ði ¼ 1; 2Þ:
ð2Þ
It is assumed that the cells are sufficiently thin so that transverse diffusion across them can be considered to be instantaneous. Thus we have two identical regions, divided by a semipermeable membrane which allows the passage of reactant A only, with the same reaction taking place in each region. We find solution is more complicated than that described in [11]. We know that there are some spatially uniform solution branches, one of which is the symmetric solution and the other spatially uniform asymmetric solutions.These solution will primarily be discussed. Here we find that there are bifurcations to pattern forms (diffusion-driven instabilities) on both the symmetric and asymmetric solution branches. These patterns emerge as pitchfork bifurcations from the appropriate spatially uniform solution. We find that the system (1) has the three spatially uniform steady states, a1s ¼ a2s ¼ l1m ; b1s ¼ b2s ¼ l
ð3Þ
and l þ ð2lÞ1m ; a l 1m a2s ¼ þ ð2lÞ ; a a1s ¼
a2s ¼ ð2lÞ1m ; a1s ¼ ð2lÞ
1m
;
b1s ¼ 0;
b2s ¼ 2l;
ð4aÞ
b2s ¼ 0;
b1s ¼ 2l:
ð4bÞ
The steady state (4a) and (4b)are identical. In the Section 2, we discuss primarily the stability of the spatially ~>l ~M l ¼ lm is bifurcation parameuniform steady state (3) and (4). The steady state (3) is stable if l c ðaÞ (~ m b ~>l ~c (~ ter)and the steady state (4) is stable for l l ¼ ð2lÞ is bifurcation parameter). These steady states which,
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
887
in the without stirred system, are stable for some range of parameter values, may in the presence of diffusion, exhibit diffusion-driven (Turing) instability. In Sections 3–5, the local behavior of Turing instability is obtained via weakly nonlinear theory. 2. Linear stability analysis Here we consider the temporal stability of the uniform state (3) and (4) when subject to a disturbance of small amplitude. We suppose that the disturbance is imposed at t ¼ 0 and can be represented by the initial conditions ai0 ðxÞ; ai ðx; 0Þ ¼ ais þ e 0 < x < 1; ði ¼ 1; 2Þ: ð5Þ bi ðx; 0Þ ¼ bis þ e bi0 ðxÞ; where ðais ; bis Þ are steady state and ð ai0 ; bi0 Þ are bounded functions, 0 < e 1 is a measure of the amplitude of the initial disturbance. With e 1 we look for a solution of Eq. (1) together with conditions (2) and (5) in the form ai ðx; tÞ ¼ ais þ eai0 ðx; tÞ þ ; bi ðx; tÞ ¼ bis þ ebi0 ðx; tÞ þ ;
ði ¼ 1; 2Þ:
ð6Þ
where ai0 , bi0 are Oð1Þ as e ! 0. After substitution from (6) into (1), expanding and neglecting terms Oðe2 Þ, we obtain the following linear equations governing the perturbed quantities ai0 ; bi0 , 8 2 oa10 > > ¼ k ooxa210 ðbm1s þ aÞa10 ma1s bm1s b10 þ aa20 ; > > ot > > > > 2 ob > 10 > < ¼ kd ooxb210 þ bm1s a10 þ ðma1s bm1 1Þb10 ; 1s ot 0 < x < 1; t > 0 ð7Þ > oa20 o2 a20 m m > > ¼ k ox2 ðb2s þ aÞa20 ma2s b2s b20 þ aa10 ; > > ot > > > > ob20 2 > : ¼ kd ooxb220 þ bm2s a20 þ ðma2s bm1 1Þb20 : 2s ot Further substitution from (6)into (2) and (5) leads to the initial and boundary conditions in terms of ai0 and bi0 as ( ai0 ðxÞ; bi0 ðx; 0Þ ¼ bi0 ðxÞ; 0 < x < 1; ði ¼ 1; 2Þ ai0 ðx; 0Þ ¼ ð8Þ oai0 obi0 ¼ ox ¼ 0; on x ¼ 0; 1; ði ¼ 1; 2Þ: ox The form of the boundary condition suggests looking for a solution of Eq. (7)in the form 8 1 P > > Ain ext cosðnpxÞ; < ai0 ðx; tÞ ¼ n¼0 ði ¼ 1; 2Þ: 1 P > > : bi0 ðx; tÞ ¼ Bin ext cosðnpxÞ;
ð9Þ
n¼0
After substituting expansions (9) into (7), we see that there are non-trivial solutions provided that x satisfies the equation detðMÞ ¼ 0; where 2
x þ k 2n þ bm1s þ a 6 bm1s 6 M ¼6 4 a 0
ma1s bm1 1s 2 x þ dk n þ 1 ma1s bm1 1s
a 0
0 0
0 0
x þ k 2n þ bm2s þ a bm2s
ma2s bm1 2s x þ dk 2n þ 1 ma2s bm1 2s
3 7 7 7; 5
888
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
with k 2n ¼ kn2 p2 ðn ¼ 0; 1; 2; ; 1Þ. On expanding this necessary condition, we find that x satisfies the characteristic equation, Þðn2 g2 þ ma2s b2m1 Þ a2 g1 g2 ¼ 0; ðn1 g1 þ ma1s b2m1 1s 2s k 2n
dk 2n
bmis
ð10Þ
mais bm1 is ;
þ1 ði ¼ 1; 2Þ. There are essentially two cases to conwhere ni ¼ x þ þ þ a; gi ¼ x þ sider, in which ais ; bis ði ¼ 1; 2Þ correspond to the steady state (3) or either of the steady states (4). 2.1. The neutral curves corresponding to the symmetric steady state In the case of the steady state (3), the characteristic Eq. (10) factorizes to give f ðxÞgðxÞ ¼ 0; where f ðxÞ ¼ x2 þ p1 ð~ l; k n Þx þ q1 ð~ l; k n Þ; gðxÞ ¼ x2 þ p2 ð~ l; k n Þx þ q2 ð~ l; k n Þ and ~ þ 1 m; l; k n Þ ¼ ðd þ 1Þk 2n þ l p1 ð~
q1 ð~ l; k n Þ ¼ ð1 þ dk 2n Þ~ l þ k 2n ðdk 2n þ 1 mÞ;
~ þ 1 m þ 2a; l; k n Þ ¼ ðd þ 1Þk 2n þ l p2 ð~
q2 ð~ l; k n Þ ¼ ð1 þ dk 2n Þ~ l þ ðk 2n þ 2aÞðdk 2n þ 1 mÞ;
l ~ ¼ l1m with l ¼ lm (the scale of initial concentration of reaction A and autocatalyst B)is bifurcation parameter. We define thatf ðxÞ ¼ 0 corresponds to the equation for the linear stability of the uncoupled system, gðxÞ ¼ 0 is that of the coupled system. We label the roots of f ðxÞ ¼ 0 as x 1 , and those of gðxÞ ¼ 0 as x2 . These roots are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð11Þ xi ¼ pi p2i 4qi ði ¼ 1; 2Þ: 2 The local stability of the steady state is then determined by examining the behavior of x 1 and x2 in the first ~Þ plane. In [11] it was seen that the determination of the neutral curves was a necessary quadrant of the ðk n ; l precursor to determining where the spatially uniform steady state bifurcated to stable spatially nonuniform ~Þ steady states. The neutral curve is subjected to x i ¼ 0 and Reðxi Þ ¼ 0 in the first quadrant of the ðk n ; l plane. We obtain the function which the neutral curve we label L1 , corresponding to f ðxÞ ¼ 0, is 8 < m 1 ðd þ 1Þk 2n ; 0 6 k n 6 k H n; qffiffiffiffiffiffiffi ~c1 ðk n Þ ¼ k2 ðm1dk2 Þ ð12Þ l H n : n ; ; k n 6 k n 6 m1 d 1þdk 2 n
where by
~H Þ ðk H n ;l
þ ~Þ plane where the curves Reðx is the point on the ðk n ; l 1 Þ ¼ 0 and x1 ¼ 0 intersect and is given
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 ðd 1Þ m þ 4dm þ ðm 2Þd m ; ¼ 2 2d ~H ¼ m 1 ðd þ 1Þk H2 l n :
k H2 n
The first of these branches 0 6 k n 6 k H n corresponds to Reðx1 Þ ¼ 0, while the second corresponds to x1 ¼ 0. qffiffiffiffiffiffiffi . The neutral curves L1 attains its maximum value when Note that L1 intersects the k n axis at k n ¼ m1 d pffiffiffiffi m1 k 2max ¼ ; d pffiffiffi m1Þ2 ~ is consider 0 < d < ð m1 , the corresponding value of l pffiffiffiffi ð m 1Þ2 ~ð1Þ : l ¼ max d
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
889
We now consider the neutral curves which we label L2 corresponding to gðxÞ ¼ 0 is given by ~ c2 ¼ m 1 2a ð1 þ dÞk 2n l
ð13Þ
and ~c2 ¼ l
ðk 2n þ 2aÞðm 1 dk 2n Þ ; 1 þ dk 2n
ð14Þ
respectively. The latter curve attains its maximum value when pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mð1 2adÞ 1 2 k max ¼ d m1 ~ is provided that ad < 2m , the corresponding value of l pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð m 1 2adÞ ~ð2Þ : l max ¼ d
ð15Þ
If ad > 12, then (14) is strictly decreasing. , In order to discuss the neutral curve L2 , we consider three ranges of a separately, when a > m1 2 ~ þ 1 m þ 2a is strictly positive for all (positive) parameter values. Therefore, the p2 ð~ l; k n Þ ¼ ðd þ 1Þk 2n þ l l; k n Þ < 0Þ. neutral curve L2 is made up of the only branch (14). On L2 , q2 ¼ 0 ðxþ 2 ¼ 0 and x2 ¼ p 2 ð~ þ Any point above the neutral curve L2 have both x2 and x2 are real and negative, while just below L2 , one m1 ~>l ~ð2Þ point has at least xþ 2 is real and positive and x2 is real and negative. Therefore, when a > 2 , if l max , any point in the regions are Reðx2 Þ < 0, the leading order terms ðai0 ; bi0 Þ ði ¼ 1; 2Þ in (6) decay exponentially ~ 0, the terms in (6) are grow exponentially as t ! 1, and the system will diverge from , the k 0 bifurcation first occurs at the uniform state (3) when perturbed. We also note that, when a > m1 2 ~ ¼ 2aðm 1Þ, corresponding to the point at which the unstable uniform state (3) emerge. l ~Þ plane where < a < m1 , the branch (13) exists but lies below (14) and there no point on the ðk n ; l When m1 2m 2 the branch (13) and (14) intersect. It dose not form part of the neutral curve because one point has at least xþ 2 m1 is positive for parameter values lying on it. Therefore, when m1 < a < ,the neutral curve consist of only 2m 2 , we find that the neutral curve is made up two parts (14). Finally, if a takes values in the range 0 < a < m1 2m is a section of the curve (13) and the second a section of (14). The curves (13) and (14) intersect in qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k 2 ¼ 1 ðm dm þ 2dÞ m2 ð1 dÞ2 þ 4mdð1 2adÞ ; ð16Þ n 2d2 with ad < 12 and ~ ¼ m 1 2a ð1 þ dÞk 2 : l n The neutral curve L2 is therefore given by 8 < m 1 2a ðd þ 1Þk 2n ; ~c2 ðk n Þ ¼ m1dk2 2 l n : ðk n þ 2aÞ; 1þdk 2 n
0 6 k n 6 k n ; qffiffiffiffiffiffiffi k n 6 k n 6 m1: d
ð17Þ
~>l ~mc ¼ maxfm 1 2a; l ~ð2Þ Therefore, when l max g, any point in this plane have Reðx2 Þ < 0 or x2 < 0, so any ~ 0 or x2 > 0, the uniform steady state (3) is unstable. On the other hand, if a > 2 , the curve L2 lies . For values of a satisfying above L1 and the stability of the steady state (3) is controlled by the signs of x 2 ^ 0 < a < m1 , we find that the first part of L lies above L , while 0 < k < k , where 1 2 n n 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ^k n ¼ 1 4d ðm 1Þð2a 1Þ ðm dðm 2 þ 2aÞÞ ; ðm dðm 2 þ 2aÞÞ 2d2 the curves L1 and L2 intersect in ^k n . The neutral curves are given in Fig. 1.
890
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
Fig. 1. Sketches of the neutral curves L1 and L2 .
In order to see whether pattern formation can occur we must examine those points at which the curves k n ¼ constant ðn ¼ 0; 1; 2 Þ intersect the neutral curves L1 and L2 . At each neutral mode one the possibility of the emergence of a spatially nonuniform solution (providing n–0). As we have already seen, pattern formation may arise through a Turing instability for certain ranges of parameter values. Our primary concern is to establish whether diffusion driven instability can occur: a necessary condition is that the value of bifur~ satisfies l ~>l ~M cation parameter l c ðaÞ where ( ~ð2Þ a > m1 l max ; 2 ~M ð18Þ l c ðaÞ ¼ ~ð2Þ maxfm 1; l 0 < a < m1 : max g; 2 ~ð2Þ Therefore, if a < m1 a necessary condition for Turing instability is give by l max > ðm 1Þ together with 2 ad < 12 ð1 m1 Þ, i.e. pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 1 2ad > dðm 1Þ , Turing instability is possible provided that 0 < d < dmin , where Therefore, when a < m1 2 pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi m 1ð m 1 2aÞ dmin ¼ m 1 þ 2a pffiffiffi m and dmin ð12Þ ¼ m1 . Thus In 0 < a < 12, dmin is a monotonically decreasing function of a with dmin ð0Þ ¼ mþ12 m1 m the effect of coupling is to increase the range of values of the diffusion coefficient ratio d available for possible pattern formation. Finally we note that the curve L2 lies above the curve L1 (which is the neutral curve for the uncoupled sys~ as the bifurcation parameter, the primary bifurcations from the spatially symmetric tem). Thus regarding l steady state (3) to stable pattern forms will arise as a result of the coupling and consequently not be seen on the uncoupled system. 2.2. The neutral curves corresponding to the asymmetric steady state In this section, we consider the stability of the asymmetric steady state (4) have the characteristic Eq. (10) given by ðx þ dk 2n þ 1Þðx3 þ a1 x2 þ a2 x þ a3 Þ ¼ 0; where ~ þ 2a þ 1 m þ ð2 þ dÞk 2n ; a1 ¼ l ~ða þ 1 þ ðd þ 1Þk 2n Þ þ ð1 þ 2dÞk 4n þ 2ððd þ 1Þa þ 1 mÞk 2n þ 2að1 mÞ; a2 ¼ l ~ðdk 4n þ ð1 þ adÞk 2n þ aÞ þ k 2n ðdk 4n þ ð1 þ 2ad mÞk 2n þ 2að1 mÞÞ; a3 ¼ l ~ ¼ ba2s2s ¼ ð2lÞm is bifurcation parameter. with l
ð19Þ
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
891
We label the solutions of the cubic equation in (19) as x1 ; x2 and x3 , and employ the Routh–Hurwitz conditions to find regions of parameter space where Reðx1 ; x2 ; x3 Þ < 0. The Routh–Hurwitz conditions for the real parts of the roots to be negative, and hence for the steady state to be stable, are a1 > 0; a3 > 0; a1 a2 a3 > 0:
ð20Þ
The possibility of primary bifurcations to steady spatially nonuniform solutions is given by the neutral curves which is obtained by letting a1 a2 a3 ¼ 0 and a3 ¼ 0. For given m, we may plot the curves a1 a2 a3 ¼ 0 and ~Þ plane. a3 ¼ 0; a1 ¼ 0 in the ðk n ; l On the boundary curves we can evaluate the eigenvalues and we find that if a3 ¼ 0 when a1 > 0 and a1 a2 a3 > 0 the one root of the cubic is zero and the other two have negative real parts. Crossing this curve corresponds to the bifurcation to spatially nonuniform steady state solutions. If a1 a2 a3 ¼ 0 when a1 > 0 and a3 > 0 then one root is negative and the other two form an imaginary pair, and crossing this curve may lead to the bifurcation to unsteady spatially nonuniform solutions. Diffusion-driven instability to pattern formation is possible if the first of the above conditions is realized, which is possible only if the maximum of ~ given by (21) the curve a3 ¼ 0 lie above the value of l ~bc l
1 ½mð1 þ 3aÞ 1 4a 2a2 þ ¼ 2ð1 þ aÞ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðma þ 1 mÞ þ 4a2 ða2 þ ma þ m 1Þ:
ð21Þ
~ is monotone decreasing in the ðk n ; l ~Þ plane and thus any neutral model with n P 1 and lying on The curve l ~l ~ bc giving the possibility of primary bifurcations ranges of values of k n over which the a3 ¼ 0 curve has l to steady spatially nonuniform solutions. Finally, we note that these bifurcations occur at a value of ~¼l ~c ðk n Þ, where l ~c ¼ l
k 2n ð2aðm 1Þ þ ðm 1 2adÞk 2n dk 4n Þ dk 4n þ ð1 þ adÞk 2n þ a
~bc . ~c > l provided that l 3. Perturbation analysis for 0 < a 1 on the symmetric branch ~¼l ~b corresponds to the uncoupled bifurcation We suppose that the mode k b is unstable and assume that l point, that is ~b ¼ l
k 2b ðm 1 dk 2b Þ ; for some positive interge b 1 þ dk 2b
ð22Þ
~ given by, when 0 < a 1, we introduce perturbation of bifurcation parameter l a; ~¼l ~b þ l l
ð23Þ 1
is Oð1Þ, we expand the solution in powers of a2 , and introduce a long timescale s ¼ at, writing where l 8 1 3 a1 ðx; t; sÞ ¼ l1m þ p1 ðx; t; sÞa2 þ p2 ðx; t; sÞa þ p3 ðx; t; sÞa2 þ Oða2 Þ; > > > > > < b1 ðx; t; sÞ ¼ l þ q ðx; t; sÞa12 þ q ðx; t; sÞa þ q ðx; t; sÞa32 þ Oða2 Þ; 1 2 3 ð24Þ > a ðx; t; sÞ ¼ l1m þ r ðx; t; sÞa12 þ r ðx; t; sÞa þ r ðx; t; sÞa32 þ Oða2 Þ; > 2 1 2 3 > > > : 1 3 b2 ðx; t; sÞ ¼ l þ s1 ðx; t; sÞa2 þ s2 ðx; t; sÞa þ s3 ðx; t; sÞa2 þ Oða2 Þ: 1
On substituting the above expressions (22)–(24) into the governing Eq. (1), we find that the solution at Oða2 Þ is
892
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
p1
¼ A1 ðsÞ
q1 where
d1 d2
d1
cosðbpxÞ;
d2
r1
s1
d1 ¼ A2 ðsÞ cosðbpxÞ; d2
m k2 þ~ , b lb
¼ p2
and amplitudes A1 and A2 are functions of the slow variable s. At OðaÞ, we obtain ! d1 A21 D 1 4dk 2b A21 D 1 ¼ cosðbpxÞ; þ B1 ðsÞ cosð2bpxÞ þ 2DðbÞ 4k 2b 2~ lb 0 q2 d2 ! r2 d1 A22 D 1 4dk 2b A22 D 1 ¼ cosðbpxÞ; þ B2 ðsÞ cosð2bpxÞ þ 2 2DðbÞ 4k b 2~ lb 0 s2 d2
1
~bm d 22 ½m 1þ where B1 and B2 are the function of s which will be determined at higher order, D ¼ m2 l 3 d1 2 2 2 ~b ; DðbÞ ¼ m~ 2 d2 l lb ðm 1 4dk b Þð~ lb þ 4k b Þ. At Oða Þ, we obtain the equation of ðp3 ; q3 Þ and ðr3 ; s3 Þ,
where
" Mb ¼
!
p1 r1 ¼ G1 þ Mb þ 1 q3 0 1 r3 r 1 p1 þ ¼ G2 þ Mb s3 0 1 p3
op1 os ; oq1 os ! or1 os ; os1 os
1
2
d ~b k dx 2 l
m
~b l
d kd dx 2 þ ðm 1Þ
2
ð25Þ ð26Þ
# :
with G1 ¼ G1 ðp1 ; q1 ; p2 ; q2 Þ; G2 ¼ G2 ðr1 ; s1 ; r2 ; s2 Þ. To obtain equation for the amplitudes A1 and A2 of the leading order solution, we invoke the Fredholm alternative for the solution of Eqs. (25) and (26). After some algebra, we obtain, ( dA 1 ¼ E½A1 ð lb þ GA21 Þ RðA1 A2 Þ; ds ð27Þ dA2 ¼ E½A2 ð lb þ GA22 Þ RðA2 A1 Þ; ds pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi ð1þdk 2b Þ2 2 2 l >0 where E > 0 for m2þ2 m 4 < dk 2b < m 1 and E < 0 for 0 < dk 2b < m2þ2 m 4, R ¼ ð~l m~ 2 2 > 0, b ¼ mk 2b b þk b Þ and G is the Landau constant is given by G¼
2
~bm d 22 l
m d1 1 þ 4dk 2b DðbÞ 2 2 ~b þ mðm 1Þk b m~ þ G1 ; m1þ2 l lb ð1 þ dk b Þ m~ lb d2 4 2DðbÞ 2~ lb
~b dd 12 þ mðm1Þðm2Þ with G1 ¼ 3mðm1Þ : For given m, we may plot the curves G ¼ GðhÞðh ¼ dkb2 p2 Þ. The steady l 8 8 state of amplitude Eq. (27) is given by ðaÞ A1s ¼ A2s ¼ 0; lb ; ðbÞ A21s ¼ A22s ¼ G b 2R l ; A2s ¼ A1s ; ðcÞ A21s ¼ G qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðdÞ A21s ¼
b Rl
ð lb RÞ2 4R2 2G
;
A22s ¼
b Rl
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð lb RÞ2 4R2 2G
:
We give the stability conditions for these steady states for E > 0 and similarities discuss them for E < 0. The b > 2R. For given parameter values such that G > 0, the second steady state (b) is steady state (a) is stable if l b < R, the third (c) is stable if l b < 2R, (d) is unstable for all values of l . On the other hand, for stable if l G < 0, all the non-trivial steady states are unstable. Therefore if G > 0, the stable pattern formation bifurcation solutions are given by the steady state (b) and (c).
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
893
4. The analysis of the strong coupling In this section, we discuss the behavior of the nonhomogeneous solution branches away from their bifurcation points, for a 1,as [19]. We note that the solution develops on two timescale: a fast timescale given by T ¼ at and the longer timescale t. The equations are rewritten in terms of the independent variables T and t and the concentrations are expanded in powers of a1 as follows ai ðx; t; T Þ ¼ ai0 ðx; t; T Þ þ a1 ai1 ðx; t; T Þ þ Oða2 Þ; ði ¼ 1; 2Þ ð28Þ bi ðx; t; T Þ ¼ bi0 ðx; t; T Þ þ a1 bi1 ðx; t; T Þ þ Oða2 Þ and the time derivative are accordingly replaced using o o o ¼ þa : ot ot oT
ð29Þ
Substitution from (28) and (29) into (1), we obtain the equation of leading order ða10 ; a20 Þ at OðaÞ ( oa 10 ¼ a20 a10 ; oT oa20 oT
¼ a10 a20 ;
which has the solution a10 ¼ A1 ðx; tÞ þ A2 ðx; tÞe2T ; a20 ¼ A1 ðx; tÞ A2 ðx; tÞe2T ; where Ai ðx; tÞði ¼ 1; 2Þ are undetermined functions which satisfy the initial and boundary conditions. At OðaÞ, the leading order for ðb10 ; b20 Þ satisfied ob10 ob20 ¼ ¼ 0; oT oT which has the solution b10 ¼ B1 ðx; tÞ; b20 ¼ B2 ðx; tÞ; where B1 ðx; tÞ and B2 ðx; tÞ are again functions which satisfy the initial and boundary conditions. We find that, after eliminating the secular terms that arise at Oð1Þ, the equations reduce to 8 > oA1 o 2 A1 > > k ¼ l 12 A1 ðBm1 þ Bm2 Þ; > 2 > ot ox > > < oB1 o 2 B1 ð30Þ kd ¼ A1 Bm1 B1 ; 2 > ot ox > > > 2 > > > : oB2 kd o B2 ¼ A1 Bm2 B2 : ot ox2 The concentrations of reactant A in the two cells equilibrate to A1 ðx; tÞ on the fast timescale, whereas the concentrations of autocatalyst B evolve on the slower timescale and are governed by the Eq. (30) above. Eq. (30) have three spatially uniform steady state: ðA1s ; B1s ; B2s Þ ¼ ðð2lÞ
1m
; 0; 2lÞ;
ðasÞ
ðA1s ; B1s ; B2s Þ ¼ ðð2lÞ
1m
; 2l; 0Þ;
ðbsÞ
1m
ðA1s ; B1s ; B2s Þ ¼ ððlÞ
; l; lÞ:
ðcsÞ
In the absence of diffusion, pffiffiffiwe 2 find that the steady state (cs) is unstable and that the steady states (as,bs) are m ~ > maxfm 1; ð md1Þ gð~ l ¼ ð2lÞ Þ.This is as we would expect from letting a ! 1 in the results of stable if l the linear stability analysis. Therefore we conclude that, in the limit as a ! 1, patterns can bifurcate only from the asymmetric steady states (as,bs), where
894
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
~¼ l
2k 2n ðm 1 dk 2n Þ ; 1 þ dk 2n
which is precisely the limit we obtain by allowing a ! 1 in a3 ¼ 0. 5. Bifurcation analysis on the asymmetric branches In this section, we analyse the forms of the spatially nonuniform steady solutions close to their bifurcation points on the asymmetric solution branches (4). We suppose that the steady state (4) loses stability such that one real (simple) eigenvalue passes through zero. This corresponds to a point on the a3 ¼ 0 part of the neutral ~ given by curve discussed earlier. In this case the steady state loses stability at the value of l ~b ¼ l
k 2b ½2aðm 1Þ ð1 m þ 2adÞk 2b dk 4b ; dk 4b þ ð1 þ adÞk 2b þ a
for some positive interge b:
ð31Þ
We take 0 < e 1 as measure of the perturbation and introduce the small perturbation of bifurcation parameter via e2 ; ~¼l ~b þ l l
ð32Þ
¼ 1. We look for a solution of (1) in the form where l 8 > a1 ðx; t; sÞ ¼ la þ ð2lÞ1m þ ep1 þ e2 p2 þ e3 p3 þ Oðe4 Þ; > > > < b ðx; t; sÞ ¼ 0 þ eq þ e2 q þ e3 q þ Oðe4 Þ; 1 1 2 3 1m 2 3 4 > a ðx; t; sÞ ¼ ð2lÞ þ er þ e r > 2 1 2 þ e r3 þ Oðe Þ; > > : 2 3 4 b2 ðx; t; sÞ ¼ 2l þ es1 þ e s2 þ e s3 þ Oðe Þ;
ð33Þ
with boundary conditions opi oqi ori osi ¼ ¼ ¼ ¼ 0ði ¼ 1; 2; 3; 4Þ on x ¼ 0; 1 ox ox ox ox
ð34Þ
Substitution these Eqs. (31)–(34) into (1), we find that qi ¼ 0; ði ¼ 1; 2; 3Þ, and that pi ; ri ; si satisfy the following equations: at OðeÞ 0 1 p1 B C Lb @ r1 A ¼ 0; ð35Þ s1 where 0
2
d k dx 2 a
B Lb ¼ B @a 0 At Oðe2 Þ: 0
p2
1
2
0
0
1
0
a d ~b k dx 2 a l
m
~b l
d kd dx 2 þ ðm 1Þ
2
C C: A
1
h i 1 m1 C B C B lbm s21 þ m~ Lb @ r2 A ¼ @ 1 A mðm 1Þ~ lbm r1 s1 : 1 s2
At Oðe3 Þ 0
p3
1
0
0
ð36Þ
1
h i m1 1 m2 B C B C r1 þ m~ Lb @ r3 A ¼ @ 1 A l lbm ðr1 s2 þ r2 s1 Þ þ 2mðm 1Þ~ lbm s1 s2 þ mðm 1Þ~ lbm r1 s21 : 1 s3
ð37Þ
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
895
Without lose of generally, the solution of the first Eq. (35) is 0 1 0 1 d1 p1 B C B C @ r1 A ¼ A@ d 2 A cosðbpxÞ; s1
d3
where A is a measure of the amplitude of the solution which we shall determine. The solution of (36) is found to be 3 2 0 1 0 1 0 1 b1 p2 1 7 6 B C B C B C @ r2 A ¼ A2 4b0 @ 1 A þ @ b2 A cosð2bpxÞ5; s2 b3 0 where b0 ¼
D Da D 4Dk 2b ð2k 2b þ aÞ; ; b1 ¼ ð1 þ 4dk 2b Þ; b2 ¼ ð4k 2b þ aÞð1 þ 4dk 2b Þ; b3 ¼ ~ ~ ~ 2~ lb 2DðbÞ 2DðbÞ 2DðbÞ 1
~ ~b dd 23 ; DðbÞ ~b ðk 2b þ aÞ m~ with D ¼ m~ lbm d 23 ½m 1 þ l ¼ ðm 1 dk 2b Þ½k 2b ðk 2b þ 2aÞ þ l lb ðk 2b þ aÞ. We are now in a position to apply a solvability condition to (37) and hence determine A. 2
~mb ll ; A ¼ md 23 H ðk; a; dÞ 2
where
! 3ðm 1Þd 2 d 3 ðm 1Þ ð4k 2b þ aÞð1 þ 4dk 2b Þ 2k 2b d 2 ð2k 2b þ aÞ d 3 ~b ~b þ H¼ l þ : l ~ ~ 4 2~ lb 4 Dð2bÞ Dð2bÞ
For given m, we may plot the curves H ¼ H ðhÞðh ¼ dkb2 p2 Þ. The sign of H determines whether the bifurcation ¼ 1 and the bifurcation occurs in l ~ 0 then l at l ¼ 1 and the bifurcation is subcritical. ical, whereas if H < 0 then l 6. Conclusion We have shown that the coupling together of two identical cells, each governed by cubic autocatalytic kinetics, via the diffusive interchange of the reactant A, can give rise to stable pattern forms which primarily bifurcate from both nonhomogeneous spatial structures. The case when the coupling is through A is significantly different from that when the coupling is through the autocatalyst B. In the former case there are new spatially uniform steady states in addition to those which occur in the uncoupled system. We give a necessary condition that there are bifurcations from the symmetric and asymmetric branches. When the bifurcation is from the ~ (~ symmetric steady state, it always occurs at larger values of l l ¼ lm is bifurcation parameter) than for the corresponding uncoupled system; the primary bifurcation is due to the coupling. The stable pattern forms will arise as a result of the coupling and consequently not be seen on the spatially uniform steady states which occur in the uncoupled system are possible. The primary bifurcation corresponds to that in the uncoupled system, and any new stable spatially nonuniform steady states arise via secondary bifurcation from a lower solution branch. In this paper, we found that spatially nonuniform solutions lost stability via Hopf bifurcations leading to spatially nonuniform oscillatory behaviour and the conversion step became dominant and that any such oscil~ and had small basins of attraction. There were ranges of l ~ for lations existed only over very small ranges of l which b1 –b2 and the oscillations in the two cells and different amplitudes and were out of phase. In the case only period one oscillations were found. we did not find any of the secondary bifurcation sequences leading to the complex time dependent behavior that was reported for the spatially uniform system.
896
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
References [1] P. Gray, S.K. Scott, Sustained oscillations and other exotic patterns of behavior in isothermal reactions, J. Phys. Chem. 89 (1985) 22– 32. [2] P. Gray, S.K. Scott, Autocatalytic reaction in the isothermal continuous stirred tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38 (1983) 29–35. [3] P. Gray, S.K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor: qscillations and instabilities in the system A þ 2B ! 3B; B ! C, Chem. Eng. Sci. 39 (1984) 1087–1098. [4] K.F. Lin, Multiplicity, stability and dynamics for isothermal autocatalytic reaction in CSTR, Chem. Eng. Sci. 36 (1981) 1447–1452. [5] B.F. Gray, S.K. Scott, Multistability and sustained oscillations in isothermal, open system, cubic autocatalysis and the influence of competitive reactions, J. Chem. Soc., Faraday Trans. 81 (1) (1985) 1563–1567. [6] J. Billingham, D.J. Needham, The development of travelling waves in quadratic and cubic autocatalysis with unequal diffusion rates III. Large time development in quadratic autocatalysis, Quart. Appl. Math. I.2 (1992) 343–372. [7] J.H. Merkin, D.J. Needham, S.K. Scott, Coupled reaction-diffusion waves in an isothermal autocatalytic chemical system, IMA J. Appl. Math. 50 (1993) 43–76. [8] S.M. Collier, J.H. Merkin, S.K. Scott, Multistability, oscillations and travelling waves in a product-feedback autocatalator model II. The initiation and propagation of travelling waves, Phil. Trans. R. Soc. Lond. A 349 (1994) 389–415. [9] P. Gray, Review lecture: instabilities and oscillations in chemical reactions in closed and open system, in: Proceedings of the Royal Society, A415 (1988) pp. 1–29. [10] R. Hill, J.H. Merkin, D.J. Needham, Stable pattern and standing wave formation in a simple isothermal cubic autocatalytic reaction scheme, J. Eng. Math. 29 (1995) 413–436. [11] R. Hill, J.H. Merkin, Patten formation in a coupled cubic autocatalator system, IMA. J. Appl. Math. 53 (1994) 265–322. [12] M.J. Metcalf, J.H. Merkin, Reaction diffusion waves in coupled isothermal autocatalytic chemical systems, IMA J. Appl. Math. 51 (1993) 269–298. [13] Koichiro Vriu, Yoh Iwasa, Turing pattern formation with kinds of cells and a diffusive chemical, Bull. Math. Biol. 69 (2007) 2515– 2536. [14] A.A. Neville, P.C. Matthews, H.W. Byrne, Interactions between pattern formation and domain growth, Bull. Math. Biol. 68 (2006) 1975–2003. [15] J.H. Merkin, H. Sevcikova, D. Snita, The effect of an electric field on the local stoichiometry of front waves in an ionic chemical system, IMA. J. Appl. Math. 64 (2000) 157–188. [16] J.H. Merkin, H. Sevcikova, The effects of a complexing agent on travelling waves in autocatalytic systems with applied electric fields, IMA. J. Appl. Math. 70 (2005) 527–549. [17] M.E. Hubbard, J.A. Leach, J.C. Wei, Pattern formation in a 2D simple chemical system with general orders of autocatalysis and decay, IMA. J. Appl. Math. 70 (2005) 723–747. [18] M.J. Metcalf, J.H. Merkin, S.K. Scott, Oscillating wave fronts in isothermal chemical systems with arbitrary powers of autocatalysis, Proc. Royal Soc. Lond. A 447 (1994) 155–174. [19] R. Hill, J.H. Merkin, The effects of coupling on pattern formation in a simple autocatalytic system, IMA. J. Appl. Math. 54 (1995) 257–281.
Applied Mathematical Modelling 33 (2009) 884–896 www.elsevier.com/locate/apm
Stability and pattern formation in a coupled arbitrary order of autocatalysis system q Li Zhang *,1, San Yang Liu Department of Applied Mathematics XiDian University, Xi’an, ShannXi 710071, PR China Received 19 June 2007; received in revised form 9 December 2007; accepted 11 December 2007 Available online 23 December 2007
Abstract Spatiotemporal structures arising in two identical cells, each governed by arbitrary order autocatalator kinetics and coupled via the diffusive interchange of a reactant, are discussed. The stability of two homogeneous steady states is obtained with the use of linear stability analysis. By studying the linearized equations, it is found that two steady states, in the uncoupled and coupled system respectively, may give rise to the possibility of bifurcations to spatially nonuniform pattern forms. Further information about Turing bifurcation solutions close to these bifurcation points are obtained by weakly nonlinear theory. It is seen that the coupling leads to bifurcations not present in the uncoupled system which give rise to locally stable nonuniform pattern forms. Finally the stability of the equilibrium points of the amplitude equation is discussed by weakly nonlinear theory, with the bifurcation branches about small coupled system with 0 < a 1 and large coupling for a 1. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Reaction diffusion system; Stability; Bifurcation; Coupling; Pattern
1. Introduction The feature common to nearly all isothermal oscillatory reactions is autocatalysis. Autocatalysis is shared not only by the Belousov–Zhabotinskii reaction, a wide range of halide-based oscillations, and the arsenite plus iodate reaction but by numerous enzyme systems and by chain-branching reactions in the gas phase [1]. All mechanistic schemes include autocatalytic steps (or autocatalytic combinations of elementary steps). However, even the simplified forms of these schemes are still complicated, and there is still need to look at even simpler prototypes in a logically ordered way. Accordingly, we consider throughout steps with the overall stoichiometry A ! B. We may also take this as written to represent a reaction satisfying the rate law = ka. When the same change is catalyzed by a species Y and has a rate satisfying the expression kay, we may q *
1
Foundation items: National Natural Science Foundation of China (60574075). Corresponding author. Tel./fax: +86 02988202861. E-mail address: [email protected] (L. Zhang). Currently working towards the Ph.D. degree in applied mathematics at Xidian University.
0307-904X/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.12.013
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
885
represent this as A þ Y ! B þ Y . Autocatalysis corresponds to catalysis by the product B itself. Two exemplary cases [2,3] span the whole range of behavior likely to be encountered: quadratic autocatalysis rate ¼ k q ab cubic autocatalysis rate ¼ k c ab2 These may be represented by the ‘chemical’ equations A þ B ! 2B;
ðaÞ
A þ 2B ! 3B:
ðbÞ
and Normally these will arise from a subscheme of elementary steps. For instance, the arsenite- iodate reaction produces iodate ions with an experimentally determined rate law of the form 2
þ d½I =dt ¼ k 1 þ k 2 ½I ½I ½IO 3 ½H :
This does not, however, require us to believe that it requires a triple encounter as an elementary step. It is important to note (a) that reversibility makes little difference to behavior until the equilibrium constant ð½B=½AÞeq becomes small and (b) that less simple stoichiometry A þ nB ! ðn þ mÞB can readily be incorporated if required [4,5]. In a series of recent papers [6–8], we have looked at the various phenomena which can arise in a simple prototype chemical reaction due to chemical coupling. Cubic autocatalysis has been studied in many papers, for example, the review lecture of Gary [9]. In particular, Hill et al. [10] has shown the formation of stable patterns in a reaction diffusion system based on the cubic autocatalator A þ 2B ! 3B; B ! C, with the reaction taking place within a closed region, the reactant A being replenished by the slow decay of precursor P via the reaction P ! A. In [11], the spatiotemporal structures that can arise in two identical cells, each governed by cubic autocatalator kinetics and coupled via the diffusive interchange of the autacatalyst B, are discussed. In [12], the initiation and propagation of reaction diffusion travelling waves in two regions, coupled together by a linear diffusive interchange across a semipermeable membrane is considered.In [13], a three variable Turing system with two kinds of cells and a diffusive chemical is studied, and the formation and maintenance of fish skin patterns with multiple pigment cells are considered. In [14], pattern formation in biological systems for which the tissue on which the spatial pattern resides is growing at a rate which is itself regulated by the diffusible chemicals that establish the spatial pattern, is investigated. In the present context, we wish to allow a more general approach in which arbitrary order of autocatalysis can be encompassed, namely A þ mB ! ðm þ 1ÞB;
rate ¼ kabm :
This simple mth order power dependence on the autocatalyst concentration plays the role of the usual exponential dependence of the reaction rate on the temperature rise DT in the systems, which can be represented somewhat schematically as A þ T ! T þ DT ;
rate ¼ kae/DT ;
where / ¼ E=RT 2ref and T ref is some suitable reference temperature at which k is evaluated. In this paper, we consider the prototype chemical reaction scheme based on arbitrary order autocatalator P ! Aðrate k 0 pÞ;
ðdÞ
A þ mB ! ðm þ 1ÞBðrate k 1 abm Þ;
ðeÞ
B ! Cðrate k 2 bÞ;
ðfÞ
where a and b are the concentrations of the reactant A and the autocatalyst B and the k i are rate constants (both reactions are assumed to be isothermal). Here we suppose that the reaction is taking place inside a closed system with reactant A now being replenished by the slow decay of a precursor P, via the simple isothermal step (d), where p is the concentration of reactant P.
886
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
In [15], a model based on an autocatalytic, two-step reaction mechanism including two ionic components (of the same charge) and two non-ionic component, where both reactions are of second order overall, is considered when an electric field is applied to the system. In [16] the effects of electric fields on the reaction fronts that arise in a system governed by an autocatalytic reaction and a complex reaction between the autocatalyst and a complex agent are considered. In [17], Hubbard et al. investigated pattern formation in a coupled system of reaction diffusion equations in two spatial dimensions. These equations arise as a model of isothermal chemical autocatalysis with termination in which the orders of autocatalysis with termination, m and n, respectively, are such that 1 < n < m. In [18], the propagation of planar reaction diffusion waves in the isothermal autocatalytic system A þ mB ! ðm þ 1ÞB; rate kabm , is considered. Attention is paid to the case when the diffusion coefficient of reactant A is much less than that of the autocatalyst; the case when A is completely immobilized is discussed in detail. In the present paper, we consider the reaction scheme (d),(e),(f) taking place in a closed vessel without stirring. This allows for spatial variations in the reactant concentrations throughout the vessel arising via the combined mechanisms of reaction and molecular diffusion. The reaction scheme take place in two identical adjacent cells, coupled by the diffusive interchange of the reactant A. The equations governing the evolution in concentration of the reactants in the vessel are 8 2 oa1 > > ¼ k ooxa21 þ l a1 bm1 þ aða2 a1 Þ; > > ot > > > > 2 ob1 > > < ¼ kd ooxb21 þ a1 bm1 b1 ; ot ð1Þ 2 > oa2 m o a > 2 > ¼ k ox2 þ l a2 b2 þ aða1 a2 Þ; > > ot > > > > > : ob2 ¼ kd o2 b22 þ a2 bm b2 ; 2 ox ot where m > 1; m 2 N , ða1 ; b1 Þ and ða2 ; b2 Þ are the dimensionless concentrations in cells 1 and 2, respectively, a is the dimensionless coupling parameter, and d is the diffusion rate of reactant A and B. We impose zero-flux boundary condition namely oai ¼ 0; ox
obi ¼ 0; ox
on x ¼ 0; 1 ði ¼ 1; 2Þ:
ð2Þ
It is assumed that the cells are sufficiently thin so that transverse diffusion across them can be considered to be instantaneous. Thus we have two identical regions, divided by a semipermeable membrane which allows the passage of reactant A only, with the same reaction taking place in each region. We find solution is more complicated than that described in [11]. We know that there are some spatially uniform solution branches, one of which is the symmetric solution and the other spatially uniform asymmetric solutions.These solution will primarily be discussed. Here we find that there are bifurcations to pattern forms (diffusion-driven instabilities) on both the symmetric and asymmetric solution branches. These patterns emerge as pitchfork bifurcations from the appropriate spatially uniform solution. We find that the system (1) has the three spatially uniform steady states, a1s ¼ a2s ¼ l1m ; b1s ¼ b2s ¼ l
ð3Þ
and l þ ð2lÞ1m ; a l 1m a2s ¼ þ ð2lÞ ; a a1s ¼
a2s ¼ ð2lÞ1m ; a1s ¼ ð2lÞ
1m
;
b1s ¼ 0;
b2s ¼ 2l;
ð4aÞ
b2s ¼ 0;
b1s ¼ 2l:
ð4bÞ
The steady state (4a) and (4b)are identical. In the Section 2, we discuss primarily the stability of the spatially ~>l ~M l ¼ lm is bifurcation parameuniform steady state (3) and (4). The steady state (3) is stable if l c ðaÞ (~ m b ~>l ~c (~ ter)and the steady state (4) is stable for l l ¼ ð2lÞ is bifurcation parameter). These steady states which,
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
887
in the without stirred system, are stable for some range of parameter values, may in the presence of diffusion, exhibit diffusion-driven (Turing) instability. In Sections 3–5, the local behavior of Turing instability is obtained via weakly nonlinear theory. 2. Linear stability analysis Here we consider the temporal stability of the uniform state (3) and (4) when subject to a disturbance of small amplitude. We suppose that the disturbance is imposed at t ¼ 0 and can be represented by the initial conditions ai0 ðxÞ; ai ðx; 0Þ ¼ ais þ e 0 < x < 1; ði ¼ 1; 2Þ: ð5Þ bi ðx; 0Þ ¼ bis þ e bi0 ðxÞ; where ðais ; bis Þ are steady state and ð ai0 ; bi0 Þ are bounded functions, 0 < e 1 is a measure of the amplitude of the initial disturbance. With e 1 we look for a solution of Eq. (1) together with conditions (2) and (5) in the form ai ðx; tÞ ¼ ais þ eai0 ðx; tÞ þ ; bi ðx; tÞ ¼ bis þ ebi0 ðx; tÞ þ ;
ði ¼ 1; 2Þ:
ð6Þ
where ai0 , bi0 are Oð1Þ as e ! 0. After substitution from (6) into (1), expanding and neglecting terms Oðe2 Þ, we obtain the following linear equations governing the perturbed quantities ai0 ; bi0 , 8 2 oa10 > > ¼ k ooxa210 ðbm1s þ aÞa10 ma1s bm1s b10 þ aa20 ; > > ot > > > > 2 ob > 10 > < ¼ kd ooxb210 þ bm1s a10 þ ðma1s bm1 1Þb10 ; 1s ot 0 < x < 1; t > 0 ð7Þ > oa20 o2 a20 m m > > ¼ k ox2 ðb2s þ aÞa20 ma2s b2s b20 þ aa10 ; > > ot > > > > ob20 2 > : ¼ kd ooxb220 þ bm2s a20 þ ðma2s bm1 1Þb20 : 2s ot Further substitution from (6)into (2) and (5) leads to the initial and boundary conditions in terms of ai0 and bi0 as ( ai0 ðxÞ; bi0 ðx; 0Þ ¼ bi0 ðxÞ; 0 < x < 1; ði ¼ 1; 2Þ ai0 ðx; 0Þ ¼ ð8Þ oai0 obi0 ¼ ox ¼ 0; on x ¼ 0; 1; ði ¼ 1; 2Þ: ox The form of the boundary condition suggests looking for a solution of Eq. (7)in the form 8 1 P > > Ain ext cosðnpxÞ; < ai0 ðx; tÞ ¼ n¼0 ði ¼ 1; 2Þ: 1 P > > : bi0 ðx; tÞ ¼ Bin ext cosðnpxÞ;
ð9Þ
n¼0
After substituting expansions (9) into (7), we see that there are non-trivial solutions provided that x satisfies the equation detðMÞ ¼ 0; where 2
x þ k 2n þ bm1s þ a 6 bm1s 6 M ¼6 4 a 0
ma1s bm1 1s 2 x þ dk n þ 1 ma1s bm1 1s
a 0
0 0
0 0
x þ k 2n þ bm2s þ a bm2s
ma2s bm1 2s x þ dk 2n þ 1 ma2s bm1 2s
3 7 7 7; 5
888
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
with k 2n ¼ kn2 p2 ðn ¼ 0; 1; 2; ; 1Þ. On expanding this necessary condition, we find that x satisfies the characteristic equation, Þðn2 g2 þ ma2s b2m1 Þ a2 g1 g2 ¼ 0; ðn1 g1 þ ma1s b2m1 1s 2s k 2n
dk 2n
bmis
ð10Þ
mais bm1 is ;
þ1 ði ¼ 1; 2Þ. There are essentially two cases to conwhere ni ¼ x þ þ þ a; gi ¼ x þ sider, in which ais ; bis ði ¼ 1; 2Þ correspond to the steady state (3) or either of the steady states (4). 2.1. The neutral curves corresponding to the symmetric steady state In the case of the steady state (3), the characteristic Eq. (10) factorizes to give f ðxÞgðxÞ ¼ 0; where f ðxÞ ¼ x2 þ p1 ð~ l; k n Þx þ q1 ð~ l; k n Þ; gðxÞ ¼ x2 þ p2 ð~ l; k n Þx þ q2 ð~ l; k n Þ and ~ þ 1 m; l; k n Þ ¼ ðd þ 1Þk 2n þ l p1 ð~
q1 ð~ l; k n Þ ¼ ð1 þ dk 2n Þ~ l þ k 2n ðdk 2n þ 1 mÞ;
~ þ 1 m þ 2a; l; k n Þ ¼ ðd þ 1Þk 2n þ l p2 ð~
q2 ð~ l; k n Þ ¼ ð1 þ dk 2n Þ~ l þ ðk 2n þ 2aÞðdk 2n þ 1 mÞ;
l ~ ¼ l1m with l ¼ lm (the scale of initial concentration of reaction A and autocatalyst B)is bifurcation parameter. We define thatf ðxÞ ¼ 0 corresponds to the equation for the linear stability of the uncoupled system, gðxÞ ¼ 0 is that of the coupled system. We label the roots of f ðxÞ ¼ 0 as x 1 , and those of gðxÞ ¼ 0 as x2 . These roots are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð11Þ xi ¼ pi p2i 4qi ði ¼ 1; 2Þ: 2 The local stability of the steady state is then determined by examining the behavior of x 1 and x2 in the first ~Þ plane. In [11] it was seen that the determination of the neutral curves was a necessary quadrant of the ðk n ; l precursor to determining where the spatially uniform steady state bifurcated to stable spatially nonuniform ~Þ steady states. The neutral curve is subjected to x i ¼ 0 and Reðxi Þ ¼ 0 in the first quadrant of the ðk n ; l plane. We obtain the function which the neutral curve we label L1 , corresponding to f ðxÞ ¼ 0, is 8 < m 1 ðd þ 1Þk 2n ; 0 6 k n 6 k H n; qffiffiffiffiffiffiffi ~c1 ðk n Þ ¼ k2 ðm1dk2 Þ ð12Þ l H n : n ; ; k n 6 k n 6 m1 d 1þdk 2 n
where by
~H Þ ðk H n ;l
þ ~Þ plane where the curves Reðx is the point on the ðk n ; l 1 Þ ¼ 0 and x1 ¼ 0 intersect and is given
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 ðd 1Þ m þ 4dm þ ðm 2Þd m ; ¼ 2 2d ~H ¼ m 1 ðd þ 1Þk H2 l n :
k H2 n
The first of these branches 0 6 k n 6 k H n corresponds to Reðx1 Þ ¼ 0, while the second corresponds to x1 ¼ 0. qffiffiffiffiffiffiffi . The neutral curves L1 attains its maximum value when Note that L1 intersects the k n axis at k n ¼ m1 d pffiffiffiffi m1 k 2max ¼ ; d pffiffiffi m1Þ2 ~ is consider 0 < d < ð m1 , the corresponding value of l pffiffiffiffi ð m 1Þ2 ~ð1Þ : l ¼ max d
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
889
We now consider the neutral curves which we label L2 corresponding to gðxÞ ¼ 0 is given by ~ c2 ¼ m 1 2a ð1 þ dÞk 2n l
ð13Þ
and ~c2 ¼ l
ðk 2n þ 2aÞðm 1 dk 2n Þ ; 1 þ dk 2n
ð14Þ
respectively. The latter curve attains its maximum value when pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mð1 2adÞ 1 2 k max ¼ d m1 ~ is provided that ad < 2m , the corresponding value of l pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð m 1 2adÞ ~ð2Þ : l max ¼ d
ð15Þ
If ad > 12, then (14) is strictly decreasing. , In order to discuss the neutral curve L2 , we consider three ranges of a separately, when a > m1 2 ~ þ 1 m þ 2a is strictly positive for all (positive) parameter values. Therefore, the p2 ð~ l; k n Þ ¼ ðd þ 1Þk 2n þ l l; k n Þ < 0Þ. neutral curve L2 is made up of the only branch (14). On L2 , q2 ¼ 0 ðxþ 2 ¼ 0 and x2 ¼ p 2 ð~ þ Any point above the neutral curve L2 have both x2 and x2 are real and negative, while just below L2 , one m1 ~>l ~ð2Þ point has at least xþ 2 is real and positive and x2 is real and negative. Therefore, when a > 2 , if l max , any point in the regions are Reðx2 Þ < 0, the leading order terms ðai0 ; bi0 Þ ði ¼ 1; 2Þ in (6) decay exponentially ~
0 6 k n 6 k n ; qffiffiffiffiffiffiffi k n 6 k n 6 m1: d
ð17Þ
~>l ~mc ¼ maxfm 1 2a; l ~ð2Þ Therefore, when l max g, any point in this plane have Reðx2 Þ < 0 or x2 < 0, so any ~
890
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
Fig. 1. Sketches of the neutral curves L1 and L2 .
In order to see whether pattern formation can occur we must examine those points at which the curves k n ¼ constant ðn ¼ 0; 1; 2 Þ intersect the neutral curves L1 and L2 . At each neutral mode one the possibility of the emergence of a spatially nonuniform solution (providing n–0). As we have already seen, pattern formation may arise through a Turing instability for certain ranges of parameter values. Our primary concern is to establish whether diffusion driven instability can occur: a necessary condition is that the value of bifur~ satisfies l ~>l ~M cation parameter l c ðaÞ where ( ~ð2Þ a > m1 l max ; 2 ~M ð18Þ l c ðaÞ ¼ ~ð2Þ maxfm 1; l 0 < a < m1 : max g; 2 ~ð2Þ Therefore, if a < m1 a necessary condition for Turing instability is give by l max > ðm 1Þ together with 2 ad < 12 ð1 m1 Þ, i.e. pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 1 2ad > dðm 1Þ , Turing instability is possible provided that 0 < d < dmin , where Therefore, when a < m1 2 pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi m 1ð m 1 2aÞ dmin ¼ m 1 þ 2a pffiffiffi m and dmin ð12Þ ¼ m1 . Thus In 0 < a < 12, dmin is a monotonically decreasing function of a with dmin ð0Þ ¼ mþ12 m1 m the effect of coupling is to increase the range of values of the diffusion coefficient ratio d available for possible pattern formation. Finally we note that the curve L2 lies above the curve L1 (which is the neutral curve for the uncoupled sys~ as the bifurcation parameter, the primary bifurcations from the spatially symmetric tem). Thus regarding l steady state (3) to stable pattern forms will arise as a result of the coupling and consequently not be seen on the uncoupled system. 2.2. The neutral curves corresponding to the asymmetric steady state In this section, we consider the stability of the asymmetric steady state (4) have the characteristic Eq. (10) given by ðx þ dk 2n þ 1Þðx3 þ a1 x2 þ a2 x þ a3 Þ ¼ 0; where ~ þ 2a þ 1 m þ ð2 þ dÞk 2n ; a1 ¼ l ~ða þ 1 þ ðd þ 1Þk 2n Þ þ ð1 þ 2dÞk 4n þ 2ððd þ 1Þa þ 1 mÞk 2n þ 2að1 mÞ; a2 ¼ l ~ðdk 4n þ ð1 þ adÞk 2n þ aÞ þ k 2n ðdk 4n þ ð1 þ 2ad mÞk 2n þ 2að1 mÞÞ; a3 ¼ l ~ ¼ ba2s2s ¼ ð2lÞm is bifurcation parameter. with l
ð19Þ
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
891
We label the solutions of the cubic equation in (19) as x1 ; x2 and x3 , and employ the Routh–Hurwitz conditions to find regions of parameter space where Reðx1 ; x2 ; x3 Þ < 0. The Routh–Hurwitz conditions for the real parts of the roots to be negative, and hence for the steady state to be stable, are a1 > 0; a3 > 0; a1 a2 a3 > 0:
ð20Þ
The possibility of primary bifurcations to steady spatially nonuniform solutions is given by the neutral curves which is obtained by letting a1 a2 a3 ¼ 0 and a3 ¼ 0. For given m, we may plot the curves a1 a2 a3 ¼ 0 and ~Þ plane. a3 ¼ 0; a1 ¼ 0 in the ðk n ; l On the boundary curves we can evaluate the eigenvalues and we find that if a3 ¼ 0 when a1 > 0 and a1 a2 a3 > 0 the one root of the cubic is zero and the other two have negative real parts. Crossing this curve corresponds to the bifurcation to spatially nonuniform steady state solutions. If a1 a2 a3 ¼ 0 when a1 > 0 and a3 > 0 then one root is negative and the other two form an imaginary pair, and crossing this curve may lead to the bifurcation to unsteady spatially nonuniform solutions. Diffusion-driven instability to pattern formation is possible if the first of the above conditions is realized, which is possible only if the maximum of ~ given by (21) the curve a3 ¼ 0 lie above the value of l ~bc l
1 ½mð1 þ 3aÞ 1 4a 2a2 þ ¼ 2ð1 þ aÞ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðma þ 1 mÞ þ 4a2 ða2 þ ma þ m 1Þ:
ð21Þ
~ is monotone decreasing in the ðk n ; l ~Þ plane and thus any neutral model with n P 1 and lying on The curve l ~
k 2n ð2aðm 1Þ þ ðm 1 2adÞk 2n dk 4n Þ dk 4n þ ð1 þ adÞk 2n þ a
~bc . ~c > l provided that l 3. Perturbation analysis for 0 < a 1 on the symmetric branch ~¼l ~b corresponds to the uncoupled bifurcation We suppose that the mode k b is unstable and assume that l point, that is ~b ¼ l
k 2b ðm 1 dk 2b Þ ; for some positive interge b 1 þ dk 2b
ð22Þ
~ given by, when 0 < a 1, we introduce perturbation of bifurcation parameter l a; ~¼l ~b þ l l
ð23Þ 1
is Oð1Þ, we expand the solution in powers of a2 , and introduce a long timescale s ¼ at, writing where l 8 1 3 a1 ðx; t; sÞ ¼ l1m þ p1 ðx; t; sÞa2 þ p2 ðx; t; sÞa þ p3 ðx; t; sÞa2 þ Oða2 Þ; > > > > > < b1 ðx; t; sÞ ¼ l þ q ðx; t; sÞa12 þ q ðx; t; sÞa þ q ðx; t; sÞa32 þ Oða2 Þ; 1 2 3 ð24Þ > a ðx; t; sÞ ¼ l1m þ r ðx; t; sÞa12 þ r ðx; t; sÞa þ r ðx; t; sÞa32 þ Oða2 Þ; > 2 1 2 3 > > > : 1 3 b2 ðx; t; sÞ ¼ l þ s1 ðx; t; sÞa2 þ s2 ðx; t; sÞa þ s3 ðx; t; sÞa2 þ Oða2 Þ: 1
On substituting the above expressions (22)–(24) into the governing Eq. (1), we find that the solution at Oða2 Þ is
892
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
p1
¼ A1 ðsÞ
q1 where
d1 d2
d1
cosðbpxÞ;
d2
r1
s1
d1 ¼ A2 ðsÞ cosðbpxÞ; d2
m k2 þ~ , b lb
¼ p2
and amplitudes A1 and A2 are functions of the slow variable s. At OðaÞ, we obtain ! d1 A21 D 1 4dk 2b A21 D 1 ¼ cosðbpxÞ; þ B1 ðsÞ cosð2bpxÞ þ 2DðbÞ 4k 2b 2~ lb 0 q2 d2 ! r2 d1 A22 D 1 4dk 2b A22 D 1 ¼ cosðbpxÞ; þ B2 ðsÞ cosð2bpxÞ þ 2 2DðbÞ 4k b 2~ lb 0 s2 d2
1
~bm d 22 ½m 1þ where B1 and B2 are the function of s which will be determined at higher order, D ¼ m2 l 3 d1 2 2 2 ~b ; DðbÞ ¼ m~ 2 d2 l lb ðm 1 4dk b Þð~ lb þ 4k b Þ. At Oða Þ, we obtain the equation of ðp3 ; q3 Þ and ðr3 ; s3 Þ,
where
" Mb ¼
!
p1 r1 ¼ G1 þ Mb þ 1 q3 0 1 r3 r 1 p1 þ ¼ G2 þ Mb s3 0 1 p3
op1 os ; oq1 os ! or1 os ; os1 os
1
2
d ~b k dx 2 l
m
~b l
d kd dx 2 þ ðm 1Þ
2
ð25Þ ð26Þ
# :
with G1 ¼ G1 ðp1 ; q1 ; p2 ; q2 Þ; G2 ¼ G2 ðr1 ; s1 ; r2 ; s2 Þ. To obtain equation for the amplitudes A1 and A2 of the leading order solution, we invoke the Fredholm alternative for the solution of Eqs. (25) and (26). After some algebra, we obtain, ( dA 1 ¼ E½A1 ð lb þ GA21 Þ RðA1 A2 Þ; ds ð27Þ dA2 ¼ E½A2 ð lb þ GA22 Þ RðA2 A1 Þ; ds pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi ð1þdk 2b Þ2 2 2 l >0 where E > 0 for m2þ2 m 4 < dk 2b < m 1 and E < 0 for 0 < dk 2b < m2þ2 m 4, R ¼ ð~l m~ 2 2 > 0, b ¼ mk 2b b þk b Þ and G is the Landau constant is given by G¼
2
~bm d 22 l
m d1 1 þ 4dk 2b DðbÞ 2 2 ~b þ mðm 1Þk b m~ þ G1 ; m1þ2 l lb ð1 þ dk b Þ m~ lb d2 4 2DðbÞ 2~ lb
~b dd 12 þ mðm1Þðm2Þ with G1 ¼ 3mðm1Þ : For given m, we may plot the curves G ¼ GðhÞðh ¼ dkb2 p2 Þ. The steady l 8 8 state of amplitude Eq. (27) is given by ðaÞ A1s ¼ A2s ¼ 0; lb ; ðbÞ A21s ¼ A22s ¼ G b 2R l ; A2s ¼ A1s ; ðcÞ A21s ¼ G qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðdÞ A21s ¼
b Rl
ð lb RÞ2 4R2 2G
;
A22s ¼
b Rl
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð lb RÞ2 4R2 2G
:
We give the stability conditions for these steady states for E > 0 and similarities discuss them for E < 0. The b > 2R. For given parameter values such that G > 0, the second steady state (b) is steady state (a) is stable if l b < R, the third (c) is stable if l b < 2R, (d) is unstable for all values of l . On the other hand, for stable if l G < 0, all the non-trivial steady states are unstable. Therefore if G > 0, the stable pattern formation bifurcation solutions are given by the steady state (b) and (c).
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
893
4. The analysis of the strong coupling In this section, we discuss the behavior of the nonhomogeneous solution branches away from their bifurcation points, for a 1,as [19]. We note that the solution develops on two timescale: a fast timescale given by T ¼ at and the longer timescale t. The equations are rewritten in terms of the independent variables T and t and the concentrations are expanded in powers of a1 as follows ai ðx; t; T Þ ¼ ai0 ðx; t; T Þ þ a1 ai1 ðx; t; T Þ þ Oða2 Þ; ði ¼ 1; 2Þ ð28Þ bi ðx; t; T Þ ¼ bi0 ðx; t; T Þ þ a1 bi1 ðx; t; T Þ þ Oða2 Þ and the time derivative are accordingly replaced using o o o ¼ þa : ot ot oT
ð29Þ
Substitution from (28) and (29) into (1), we obtain the equation of leading order ða10 ; a20 Þ at OðaÞ ( oa 10 ¼ a20 a10 ; oT oa20 oT
¼ a10 a20 ;
which has the solution a10 ¼ A1 ðx; tÞ þ A2 ðx; tÞe2T ; a20 ¼ A1 ðx; tÞ A2 ðx; tÞe2T ; where Ai ðx; tÞði ¼ 1; 2Þ are undetermined functions which satisfy the initial and boundary conditions. At OðaÞ, the leading order for ðb10 ; b20 Þ satisfied ob10 ob20 ¼ ¼ 0; oT oT which has the solution b10 ¼ B1 ðx; tÞ; b20 ¼ B2 ðx; tÞ; where B1 ðx; tÞ and B2 ðx; tÞ are again functions which satisfy the initial and boundary conditions. We find that, after eliminating the secular terms that arise at Oð1Þ, the equations reduce to 8 > oA1 o 2 A1 > > k ¼ l 12 A1 ðBm1 þ Bm2 Þ; > 2 > ot ox > > < oB1 o 2 B1 ð30Þ kd ¼ A1 Bm1 B1 ; 2 > ot ox > > > 2 > > > : oB2 kd o B2 ¼ A1 Bm2 B2 : ot ox2 The concentrations of reactant A in the two cells equilibrate to A1 ðx; tÞ on the fast timescale, whereas the concentrations of autocatalyst B evolve on the slower timescale and are governed by the Eq. (30) above. Eq. (30) have three spatially uniform steady state: ðA1s ; B1s ; B2s Þ ¼ ðð2lÞ
1m
; 0; 2lÞ;
ðasÞ
ðA1s ; B1s ; B2s Þ ¼ ðð2lÞ
1m
; 2l; 0Þ;
ðbsÞ
1m
ðA1s ; B1s ; B2s Þ ¼ ððlÞ
; l; lÞ:
ðcsÞ
In the absence of diffusion, pffiffiffiwe 2 find that the steady state (cs) is unstable and that the steady states (as,bs) are m ~ > maxfm 1; ð md1Þ gð~ l ¼ ð2lÞ Þ.This is as we would expect from letting a ! 1 in the results of stable if l the linear stability analysis. Therefore we conclude that, in the limit as a ! 1, patterns can bifurcate only from the asymmetric steady states (as,bs), where
894
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
~¼ l
2k 2n ðm 1 dk 2n Þ ; 1 þ dk 2n
which is precisely the limit we obtain by allowing a ! 1 in a3 ¼ 0. 5. Bifurcation analysis on the asymmetric branches In this section, we analyse the forms of the spatially nonuniform steady solutions close to their bifurcation points on the asymmetric solution branches (4). We suppose that the steady state (4) loses stability such that one real (simple) eigenvalue passes through zero. This corresponds to a point on the a3 ¼ 0 part of the neutral ~ given by curve discussed earlier. In this case the steady state loses stability at the value of l ~b ¼ l
k 2b ½2aðm 1Þ ð1 m þ 2adÞk 2b dk 4b ; dk 4b þ ð1 þ adÞk 2b þ a
for some positive interge b:
ð31Þ
We take 0 < e 1 as measure of the perturbation and introduce the small perturbation of bifurcation parameter via e2 ; ~¼l ~b þ l l
ð32Þ
¼ 1. We look for a solution of (1) in the form where l 8 > a1 ðx; t; sÞ ¼ la þ ð2lÞ1m þ ep1 þ e2 p2 þ e3 p3 þ Oðe4 Þ; > > > < b ðx; t; sÞ ¼ 0 þ eq þ e2 q þ e3 q þ Oðe4 Þ; 1 1 2 3 1m 2 3 4 > a ðx; t; sÞ ¼ ð2lÞ þ er þ e r > 2 1 2 þ e r3 þ Oðe Þ; > > : 2 3 4 b2 ðx; t; sÞ ¼ 2l þ es1 þ e s2 þ e s3 þ Oðe Þ;
ð33Þ
with boundary conditions opi oqi ori osi ¼ ¼ ¼ ¼ 0ði ¼ 1; 2; 3; 4Þ on x ¼ 0; 1 ox ox ox ox
ð34Þ
Substitution these Eqs. (31)–(34) into (1), we find that qi ¼ 0; ði ¼ 1; 2; 3Þ, and that pi ; ri ; si satisfy the following equations: at OðeÞ 0 1 p1 B C Lb @ r1 A ¼ 0; ð35Þ s1 where 0
2
d k dx 2 a
B Lb ¼ B @a 0 At Oðe2 Þ: 0
p2
1
2
0
0
1
0
a d ~b k dx 2 a l
m
~b l
d kd dx 2 þ ðm 1Þ
2
C C: A
1
h i 1 m1 C B C B lbm s21 þ m~ Lb @ r2 A ¼ @ 1 A mðm 1Þ~ lbm r1 s1 : 1 s2
At Oðe3 Þ 0
p3
1
0
0
ð36Þ
1
h i m1 1 m2 B C B C r1 þ m~ Lb @ r3 A ¼ @ 1 A l lbm ðr1 s2 þ r2 s1 Þ þ 2mðm 1Þ~ lbm s1 s2 þ mðm 1Þ~ lbm r1 s21 : 1 s3
ð37Þ
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
895
Without lose of generally, the solution of the first Eq. (35) is 0 1 0 1 d1 p1 B C B C @ r1 A ¼ A@ d 2 A cosðbpxÞ; s1
d3
where A is a measure of the amplitude of the solution which we shall determine. The solution of (36) is found to be 3 2 0 1 0 1 0 1 b1 p2 1 7 6 B C B C B C @ r2 A ¼ A2 4b0 @ 1 A þ @ b2 A cosð2bpxÞ5; s2 b3 0 where b0 ¼
D Da D 4Dk 2b ð2k 2b þ aÞ; ; b1 ¼ ð1 þ 4dk 2b Þ; b2 ¼ ð4k 2b þ aÞð1 þ 4dk 2b Þ; b3 ¼ ~ ~ ~ 2~ lb 2DðbÞ 2DðbÞ 2DðbÞ 1
~ ~b dd 23 ; DðbÞ ~b ðk 2b þ aÞ m~ with D ¼ m~ lbm d 23 ½m 1 þ l ¼ ðm 1 dk 2b Þ½k 2b ðk 2b þ 2aÞ þ l lb ðk 2b þ aÞ. We are now in a position to apply a solvability condition to (37) and hence determine A. 2
~mb ll ; A ¼ md 23 H ðk; a; dÞ 2
where
! 3ðm 1Þd 2 d 3 ðm 1Þ ð4k 2b þ aÞð1 þ 4dk 2b Þ 2k 2b d 2 ð2k 2b þ aÞ d 3 ~b ~b þ H¼ l þ : l ~ ~ 4 2~ lb 4 Dð2bÞ Dð2bÞ
For given m, we may plot the curves H ¼ H ðhÞðh ¼ dkb2 p2 Þ. The sign of H determines whether the bifurcation ¼ 1 and the bifurcation occurs in l ~
896
L. Zhang, S.Y. Liu / Applied Mathematical Modelling 33 (2009) 884–896
References [1] P. Gray, S.K. Scott, Sustained oscillations and other exotic patterns of behavior in isothermal reactions, J. Phys. Chem. 89 (1985) 22– 32. [2] P. Gray, S.K. Scott, Autocatalytic reaction in the isothermal continuous stirred tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38 (1983) 29–35. [3] P. Gray, S.K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor: qscillations and instabilities in the system A þ 2B ! 3B; B ! C, Chem. Eng. Sci. 39 (1984) 1087–1098. [4] K.F. Lin, Multiplicity, stability and dynamics for isothermal autocatalytic reaction in CSTR, Chem. Eng. Sci. 36 (1981) 1447–1452. [5] B.F. Gray, S.K. Scott, Multistability and sustained oscillations in isothermal, open system, cubic autocatalysis and the influence of competitive reactions, J. Chem. Soc., Faraday Trans. 81 (1) (1985) 1563–1567. [6] J. Billingham, D.J. Needham, The development of travelling waves in quadratic and cubic autocatalysis with unequal diffusion rates III. Large time development in quadratic autocatalysis, Quart. Appl. Math. I.2 (1992) 343–372. [7] J.H. Merkin, D.J. Needham, S.K. Scott, Coupled reaction-diffusion waves in an isothermal autocatalytic chemical system, IMA J. Appl. Math. 50 (1993) 43–76. [8] S.M. Collier, J.H. Merkin, S.K. Scott, Multistability, oscillations and travelling waves in a product-feedback autocatalator model II. The initiation and propagation of travelling waves, Phil. Trans. R. Soc. Lond. A 349 (1994) 389–415. [9] P. Gray, Review lecture: instabilities and oscillations in chemical reactions in closed and open system, in: Proceedings of the Royal Society, A415 (1988) pp. 1–29. [10] R. Hill, J.H. Merkin, D.J. Needham, Stable pattern and standing wave formation in a simple isothermal cubic autocatalytic reaction scheme, J. Eng. Math. 29 (1995) 413–436. [11] R. Hill, J.H. Merkin, Patten formation in a coupled cubic autocatalator system, IMA. J. Appl. Math. 53 (1994) 265–322. [12] M.J. Metcalf, J.H. Merkin, Reaction diffusion waves in coupled isothermal autocatalytic chemical systems, IMA J. Appl. Math. 51 (1993) 269–298. [13] Koichiro Vriu, Yoh Iwasa, Turing pattern formation with kinds of cells and a diffusive chemical, Bull. Math. Biol. 69 (2007) 2515– 2536. [14] A.A. Neville, P.C. Matthews, H.W. Byrne, Interactions between pattern formation and domain growth, Bull. Math. Biol. 68 (2006) 1975–2003. [15] J.H. Merkin, H. Sevcikova, D. Snita, The effect of an electric field on the local stoichiometry of front waves in an ionic chemical system, IMA. J. Appl. Math. 64 (2000) 157–188. [16] J.H. Merkin, H. Sevcikova, The effects of a complexing agent on travelling waves in autocatalytic systems with applied electric fields, IMA. J. Appl. Math. 70 (2005) 527–549. [17] M.E. Hubbard, J.A. Leach, J.C. Wei, Pattern formation in a 2D simple chemical system with general orders of autocatalysis and decay, IMA. J. Appl. Math. 70 (2005) 723–747. [18] M.J. Metcalf, J.H. Merkin, S.K. Scott, Oscillating wave fronts in isothermal chemical systems with arbitrary powers of autocatalysis, Proc. Royal Soc. Lond. A 447 (1994) 155–174. [19] R. Hill, J.H. Merkin, The effects of coupling on pattern formation in a simple autocatalytic system, IMA. J. Appl. Math. 54 (1995) 257–281.