- Email: [email protected]
Pattern formation in the Gray–Scott model
Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121 www.elsevier.com/locate/na
Pattern formation in the Gray–Scott model Je& S. McGough∗ , Kyle Riley Department of Mathematics and Computer Science, South Dakota School of Mines and Technology, 501 E. St. Joseph St., Rapid City, SD 57701, USA Received 23 August 2002; accepted 20 December 2002
Abstract We investigate an elliptic system that arises in cubic autocatalytic reactions known as the Gray–Scott model. Complicated patterns were reported by Pearson in a numerical study of this system. We produce the bifurcation analysis to support the existing numerical evidence for patterns. Speci3cally bifurcation results and C 2 bounds for nonuniform steady states are derived. ? 2003 Elsevier Ltd. All rights reserved. MSC: 35K55; 35K57; 35J60 Keywords: Partial di&erential equations; Reaction di&usion equations; Gradient bounds; Gray–Scott model; Chemical kinetics; Autocatalytic reactions
1. Introduction In the last several decades, chemical kinetics has produced a variety of phenomenon that have translated into challenging mathematical problems. A classical example is seen in the waves of the Belousov–Zhabotinskii reaction. Other examples have been produced that are not quite as complicated and require fewer species interactions, but still yield very interesting behavior. One aspect of these systems is that they do not involve thermal transfer as an essential part of the interaction. Since the 1950s, fundamental studies in reaction kinetics have focused on nonisothermal systems, i.e. where thermal feedback is a critical element. In 1968, Selkov described a particular autocatalytic model of glycolysis. The version of this model due to Gray and Scott [7] is investigated below. Gray–Scott wanted to provide the same foundation for isothermal autocatalytic systems, i.e. chemical feedback. This model becomes the basis of our paper. ∗
Corresponding author. Tel.: +1-605-394-2471; fax: +1-605-394-6078. E-mail address: [email protected] (J.S. McGough).
1468-1218/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S1468-1218(03)00020-8
106
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
In the paper by Pearson [13], very complicated patterns were described for a parabolic system known as the Gray–Scott system. Pearson did a thorough numerical study of the system and found a complex structure in the solutions. However, Pearson used a simple integration scheme and left open the question regarding numerical artifacts. We are able to con3rm the results from Pearson and show that more robust numerical schemes produce the same results. We also provide a bifurcation analysis giving the existence of nonuniform solutions for the steady-state problem, i.e. the elliptic system of Gray–Scott. The Gray–Scott model is ut = d1 Iu + F(1 − u) − uv2 ;
x ∈ ;
vt = d2 Iv − (F + k)v + uv2 ;
x ∈ ;
(1) 2
with boundary data: @u=@ = @v=@ = 0 where = [0; 2:5] , and F; k; d1 ; d2 are constants. Similar systems to the Gray–Scott system have been described in the literature, such as, the Glycolysis model [1,17], or the Brusselator model [2–5,15]. These systems can be collected into a general form ut = d1 Iu + a1 u + b1 v − f(u; v) + g1 (x; u; v);
x ∈ ;
vt = d2 Iv + a2 u + b2 v + f(u; v) + g2 (x; u; v);
x∈
(2)
with zero Neumann boundary data: @u=@ = @v=@ = 0. The parameters for the models are chosen as • Gray–Scott: a1 = −F ¡ 0, b1 = 0, a2 = 0, b2 = −(F + k) ¡ 0, f = uv2 , g1 = F, g2 = 0. • Brusselator: a1 = 0, b1 = B, a2 = 0, b2 = −(B + 1) ¡ 0, f = uv2 , g1 = 0, g2 = A. • Glycolysis: a1 = −k ¡ 0, b1 = 0, a2 = k, b2 = −1, f = uv2 , g1 = , g2 = 0. In a forthcoming paper we develop a general class of problems which cover the three systems above. This general class of problems will share some common properties in their steady-state solutions. The analysis of the Gray–Scott model will lend some insight into the behavior of the general system (2). The next section presents the analysis of the Gray–Scott model. We 3nish with some discussion of numerical results and compare our numerical results with the results given in the Pearson paper. Our results are summarized in the following: Theorem. There exists nonuniform smooth steady solutions to the Gray–Scott system (1) with d1 ; d2 ¿ d ¿ 0, F; k ¿ 0. The Gray–Scott system (1) has smooth bounded solutions and has a rich structure of steady-state solutions. There exists both uniform and nonuniform steady states. 2. Autocatalytic models Gray–Scott began with isothermal autocatalytic systems in the continuously Kowing, well-stirred, tank reactor (CSTR). In this model, isothermal implies the reaction takes
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
107
place under constant temperature, autocatalytic means the catalyst is also the product, and continuously Kowing corresponds to an open system. The well-stirred assumption involves systems that have uniform transport of the reactants. The latter assumption will be relaxed to obtain the general model. There are several well-known examples of autocatalysis as isothermal feedback: Belousov–Zhabotinskii reaction, chorite–iodide– malonic acid reaction, arsenite–iodate reaction, and enzyme systems [4,5,11,12,15]. The reactions listed above are still rather complicated. To reduce the problem for a detailed study, Gray and Scott focused on a model with the overall stoichiometry: A → B, using the reaction rate law: ka (a = [A]). When the reaction is catalyzed by some species Y (the catalyst): A + Y → B + Y, the reaction rate is kay. Our focus is on cubic autocatalysis: A + 2B → 3B, with rate=kab2 . The mass balance for this reaction is da = −k1 ab2 + kf (a0 − a); dt db = k1 ab2 − k2 b + kf (b0 − b): dt
(3)
The system in (3) can be recast into dimensionless quantities. De3ne u as the dimensionless reactant concentration, v the dimensionless catalyst concentration, t the dimensionless time, F a dimensionless “Kow” parameter (1/(residence-time)), k a dimensionless catalyst lifetime and b0 is set to zero. This yields the CSTR model: du = −uv2 + F(1 − u); dt dv (4) = uv2 − (F + k)v: dt One goal of this paper is to analyze the unstirred, or nonuniform, variant of the Gray–Scott model. This nonuniform version of the model will be designated as CFUR (continuously feed, unstirred, reactor). A functional reaction vessel has been built by Harry Swinney et al. [10] and provides some of the 3rst observational evidence of the CFUR in the lab. The speci3c Gray–Scott CFUR Model is the system (1) given in the introduction. In a paper by Pearson [13], numerical studies of the Gray–Scott model produced some interesting patterns which were proposed as similar to those found by Swinney [10]. Pearson used a simple Euler scheme to integrate the parabolic equations. We will later determine some basic bounds on the solutions which generate the complicated patterns. Additional background on the Gray–Scott equations may be found in [4,7– 9,14,16,18]. The patterns observed in these models are thought to be examples of the Turing process. 3. Gray–Scott model dynamics In this section we present an overview of the dynamics found in the Gray–Scott model. Our approach is to examine the local problem (setting the di&usion constants
108
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Fig. 1. Nullclines for the Gray–Scott model.
to zero), determine the steady states, and examine the spatially independent dynamics. We use this to then bootstrap up to study the nonuniform steady states which will be the solution to the related elliptic system. The steady states of the local dynamics (setting d1 and d2 to zero) give us the uniform steady-state solutions. The point (1; 0) is a uniform steady state and will be labelled p1 . If F ¿4(F + k)2 , or the interior of the parabolic region is de3ned by F = 12 [(1=4 − 2k) ± 1=16 − k], then there are two further restpoints which are on the intersection of the nullclines (F + k)v + Fu = F and uv = F + k (see Fig. 1). Solving these two equations gives 1 (F + k)2 1± 1−4 and u1; 2 = 2 F F (F + k)2 v1; 2 = 1∓ 1−4 : 2(F + k) F We will call p2 = (u1 ; v1 ) = (u+ ; v− ) (lower right) and p3 = (u2 ; v2 ) (upper left). This line intersects (1; 0) and (0; F=(F + k)). We will later make use of the fact that v1;2 2 = F is equivalent to F = 4(F + k)2 . This is the set of points for which there are exactly two restpoints for the vector 3eld. One last computation gives that v2 6 F(1 ∓ ) where 0 ¡ ¡ 1 and this bounds v in restpoint p2 : v2 6 F. The intersection of the four half spaces, G1 (u; v) = u ¿ 0; G2 (u; v) = v ¿ 0;
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
109
G3 (u; v) = −(u − c1 ) ¿ 0; G4 (u; v) = −((a1 + a2 )u + (b1 + b2 )v − c2 ) ¿ 0; de3nes a bounded invariant region for c1 , c2 suQciently large. Thus, the local dynamics is bounded in the invariant region de3ned by the intersection of the half planes Gi , i = 1; 2; 3; and 4. Setting the di&usion terms in (1) to zero results in ut = −uv2 + F(1 − u); vt = uv2 − (F + k)v: The intersection of the four half spaces, G1 (u; v) = u ¿ 0; G2 (u; v) = v ¿ 0; G3 (u; v) = −(u − 1) ¿ 0; G4 (u; v) = −(u + v − ) ¿ 0; de3nes a bounded invariant region for ¿ 1 + F=k. The point (1; 0) is always stable. The other two restpoint’s stability are given by the eigenvalues of the following matrix: −(F + v2 ) −2uv −(F + v2 ) −2(F + k) A= = : v2 v2 −(F + k) + 2uv (F + k) 2 The determinant is given by det(A)=(v −F)(F +k) and the trace Tr(A)=−(v2 −k). The eigenvalues are = Tr(A)=2 ± Tr(A)2 − 4 det(A)=2. The point where there is a unique restpoint also yields det(A) = 0. Thus, (F; k) = (1=16; 1=16) yields (u∗ ; v∗ ) = (1=2; 1=4). The restpoints split and travel up/down the line (F + k)v + Fu = F as one enters the region. The restpoint p2 is a saddle point. This follows from recalling that v2 6 F which forces the determinant to be positive and expression under the square root to be larger than Tr(A)2 . Restpoint p3 cannot be a saddle point, but it may have stable, unstable, or possibly spiral, node structure. A Hopf bifurcation occurs if one traverses the correct line in √ the F; k plane. For this, the trace of the linearization must be zero, v = k, and the √ determinant must be nonnegative, v ¿ F. Using the 3rst steady-state equation, one obtains u = F=(k + F):√This value may be plugged into the equation (F + k)v + Fu = F giving (F + k)2 = F k. We may compare this line (the Hopf line) to the lower branch of the steady-state bifurcation line. The Hopf line lies above and intersects the steady-state bifurcation line at (0; 0) and the turning point (1=16; 1=16) (Fig. 2).
110
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121 0.25 Static Bifurcation Stable Nodes Stable Spirals Hopf Bifurcation
0.2
F
0.15
0.1
0.05
0 0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
k
Fig. 2. The stability of the upper-left restpoint, p3 .
4. Stability of uniform solutions To determine if uniform solutions are stable, we study the linearized problem for (2). In other words, we want to see if any modes have exponentially growing terms associated with them. To do this we take the FrechSet derivative of the operator de3ning the Kow, and look for exponentially growing eigenfunctions of the space operator. With reuse of the variables u and v, and taking g1 and g2 to be constant, the linearization of (2) is ut = d1 Iu + (a1 − fu )u + (b1 − fv )v;
x ∈ ;
vt = d2 Iv + (a2 + fu )u + (b2 + fv )v;
x ∈ :
(5)
Separation of the time variable (u = e t , v = e t ) gives d1 I + (a1 − fu − ) + (b1 − fv ) ;
x ∈ ;
d2 I + (a2 + fu ) + (b2 + fv − ) ;
x ∈ :
(6)
To simplify the analysis, let c1 = a1 − fu , c2 = b1 − fv , c3 = a2 + fu , c4 = b2 + fv . De3ne !n to be the eigenfunctions on , and #n ¿ 0 the eigenvalues for I!n = −#n !n :
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
111
We may then split the eigenvector from the eigenfunction: ( ; ) = ($1 ; $2 )!n and thus we reduce to $1 c1 − d1 #n $1 c2 $1 = = ; M $2 c3 c4 − d 2 # n $2 $2 Det(M ) = (d1 + d2 )#n2 − (d2 c1 + d1 c4 )#n + c1 c4 − c2 c3 : Since d1 and d2 are positive, this quadratic is positive for large values of #n . If c1 c4 − c2 c3 ¡ 0 then we get the existence of a positive root. Hence, we get the appearance of nonuniform solutions. This is discussed in greater detail in Section 5. For the Gray–Scott system (1), the associated eigenvalue problem is d1 I − (vs2 + F + ) − 2us vs = 0; d2 I + vs2 + (2us vs − F − k − ) = 0: For the uniform solution (1; 0), this reduces to d1 I − (F + ) = 0; d2 I − (F + k + ) = 0: Following the separation of variables and Fourier analysis outlined above, we now work through the stability analysis for the Gray–Scott system. For (1; 0), we note that the eigenvalues are 1 = −F − d1 #n ¡ 0
and
2 = −F − k − d2 #n ¡ 0
for F ¿ 0. It follows that (1,0) is a stable uniform solution. For the remaining uniform solutions (if they exist), let −2(F + k) −(d1 #n + vs2 + F) M= −d2 #n + F + k vs2 and then the eigenvalues, , may be determined by det(M − I ) = 0: If Tr(M ) = −#n (d1 + d2 ) − vs2 + k and det(M ) = (d1 #n + vs2 + F)(d2 #n − F − k) + 2(F + k)vs2 = d1 d2 #n2 + (d2 (vs2 + F) − d1 (F + k))#n + (F + k)(vs2 − F) then the solutions are Tr(M ) ± Tr(M )2 − 4 det(M ) : = 2 As #n gets large, the structure of M becomes diagonally dominant, which forces the corresponding eigenvalues to be real valued and negative. The term
112
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121 0.25 Static Bifurcation Turing Bifurcation Hopf Bifurcation
0.2
F
0.15
0.1
0.05
0 0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
k
Fig. 3. Location of the unstable modes.
D = Tr(M )2 − 4 det(M ) is strictly increasing in #n (since dD=d#n = 2(d1 − d2 )2 # + constant). Thus, the magnitude of the imaginary component is strictly decreasing. One may conclude that if at # = 0 the eigenvalues are real then they are always real. Dynamical instabilities occur at the lowest modes when the eigenvalues are complex valued at # = 0. A static bifurcation can occur at higher modes, which will occur in our subsequent analysis. The boundary of the region in parameter space where (F; k) : sup R( (F; k; #)) = 0 ; #¿0
de3nes a region of instability that extends the region de3ned by the Hopf bifurcation. This curve is not the same as the curve for Hopf bifurcations. It does intersect the origin, but not the turning point of the static bifurcation curve. One may verify this by plugging into the matrix M and computing the eigenvalues. The 3rst case is = 0 and for the second case ≈ 0:021 for # ≈ 1550. The maximum of may be computed directly. Solving the equations det(M (#)) = 0 and [det(M (#))]# = 0, we obtain from the second equation # = [d1 (F + k) − d2 (v2 + F)]=2d1 d2 . Using the 3rst equation one arrives at [d1 (F + k) − d2 (v2 + F)]2 − 4d1 d2 (F + k)(v2 − F) = 0; which clearly has no solution when v2 ¡ F. First, it is informative to determine where this curve intersects the static bifurcation curve, F − 4(F + k)2 = 0. Since v2 = F, the line k = 0. One would (d1 − 2d2 )F + k = 0, which for the Gray–Scott model is √ also like to determine intersections with the Hopf curve, F k − (F + k)2 = 0. This is
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
113
0.005 Dispersion curves 0
-0.005
R (λ)
-0.01
-0.015
-0.02
-0.025
-0.03 0
1000
2000
3000 µ
4000
5000
6000
0.02 Dispersion curves 0.015 0.01 0.005
R (λ)
0 -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 0
1000
2000
3000
4000
5000
6000
µ
Fig. 4. Top: Family of instability curves for F = 0:038 where k = 0:057 is the lowest curve with increments of Ik = 0:0003. Bottom: Family of instability curves for F = 0:063 and k = 0:06149 (increments of 0:00025).
done numerically for d1 = 2(10−5 ) and d2 = 10−5 giving f = 0:0471, and k = 0:06056. Fig. 3 provides the stability information for d1 = 2(10−5 ) and d2 = 10−5 . Fig. 4 graphs a family of curves for 3xed F and increasing k. The curves are uniformly increasing in k. This gives one typical bifurcation scenario; one that the
114
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
uniform mode becomes unstable 3rst. The x-axis is really a plot of the eigenvalues for the continuous operator; which on the present scale appear dense in the interval. For the second graph in Fig. 4, the uniform mode is not the 3rst mode to become unstable. For the following, we take to be a rectangle with l1 ; l2 to be the length scales in the x; y directions. The eigenvalues for the Laplacian on the square (with zero Neumann data) are k2 m2 #km = + + ; k = 0; 1; 2; : : : ; and m = 0; 1; 2; : : : : l21 l22 For 3xed #, the collection of k and m lie on an ellipse with major and minor axes de3ned by (l1 #=+)2 and (l2 #=+)2 . Pearson studies the case where l1 = l2 = l, and so the curve corresponds to a circle. Because it is symmetric about the m = k line, any eigenvalue #km for which k = m has #km = #mk , in other words, it occurs in pairs. Thus if m = k does not produce the #km , the multiplicity of #km is even. To get eigenvalues of odd multiplicity, we need then that k = m produces an eigenvalue. Restricting, and √ relabelling, to this set in the Pearson example, yields #m = m+ 2=l. 5. A priori bounds for Gray–Scott elliptic problem To prove the existence of nonuniform solutions, a bifurcation approach is used, which is similar to the approach in [3]. There are at most three steady states. The solution (1; 0) is always stable and if it exists, the steady state p2 is unstable for the local dynamics. We focus our attention on the third uniform solution, p3 , which is denoted (us ; vs ) in the following analysis. For the reader’s convenience, system (1) is reproduced here. The functions u and v satisfy d1 Iu + F(1 − u) − uv2 = 0;
x ∈ ;
d2 Iv − (F + k)v + uv2 = 0;
x ∈ ;
(7)
with boundary data: @u=@ = @v=@ = 0 on @ . This system is shifted about the steady state i.e., let u˜ = u − us and v˜ = v − vs and so we obtain d1 Iu˜ − u˜v˜2 − 2u˜vv ˜ s − uv ˜ s2 − us v˜2 − 2us vv ˜ s − F u˜ = 0; 2
d2 Iv˜ + u˜v˜ + 2u˜vv ˜ s+
uv ˜ s2
2
x ∈ ;
+ us v˜ + 2us vv ˜ s − (F + k)v˜ = 0;
x ∈ ;
(8)
with @u=@ = @v=@ = 0 on @ . We drop the tildes for now and recast the problem as an integral equation. Choose Banach spaces X and Y where
@f @g 2; Y = (f; g) | f; g ∈ C ( ); = = 0 on @ and @n @n X = {(f; g) | f; g ∈ C ( )}
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
for some 0 ¡ ¡ 1. Next de3ne operators, L, M , and N , such that 2 −vs − F + d1 L(u; v) = [ − d1 (Iu − u); −d2 (Iv − v)]; M = vs2 and
N (u; v) =
−uv2 − 2uvvs − us v2 uv2 + 2uvvs + us v2
115
−2(F + k)
F + k + d2
:
Thus the elliptic problem in (1) may be re-expressed as L(u; v) = M (u; v) + N (u; v): To apply the standard bifurcation theorems, we must be able to invert the operator L. To do this, one must examine the invertibility of the operator (0 − I ) with zero Neumann boundary data: Iu − u = f;
· ∇u = 0;
x ∈ ; x ∈ @ :
(9)
This is a strictly elliptic operator with c 6 0, we can invoke Theorem 6.31 from [6], U for each f ∈ C ( ). U This proves and 3nd that Eq. (9) above has a unique u ∈ C 2; ( ) that the operator L above is invertible and we obtain (u; v) = L−1 [M (u; v) + N (u; v)] ≡ T (u; v): Hence, T : X → Y is a bounded operator. Y is compactly contained in X so T : X → X is a compact operator. Solving (u; v) = T (u; v) is then equivalent to solving (1). Let K(u; v) = L−1 M (u; v) and G(u; v) = L−1 N (u; v). The associated linear eigenvalue problem is K( ; u; v) − (u; v) = 0 for nontrivial u, v. In this analysis, the parameter represents one of the constants: F, k, d1 , or d2 . Let eigenvalues for 0 − I on with Neumann data be represented by #n and let !m be the associated eigenfunctions. We assume that is a domain for which #n are isolated and that there is an in3nite subset of #n , #nm , for which the algebraic multiplicity of #n is 3nite and odd and that K is Fredholm Index zero. This assumption is valid for the square domain [19]. We will drop the double subscript and refer to this subsequence as #m . This gives rise to the associated algebraic eigenvalue problem, to determine the value of that solves 2 −vs − F + d1 − d1 #m −2(F + k) = 0: det vs2 F + k + d 2 − d2 # m The solutions to the above will be denoted m± . Normally we will choose = F and so m± = Fm± . The nontrivial solution is then K( m± ; c1± !m ; c2± !m ) − (c1± !m ; c2± !m ) = 0: De3ne w = (u; v) and the operator equation is written as w − T (w) = 0:
116
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
Theorem 1. The point m± is a bifurcation point for w − T (w) = 0. Standard bifurcation results assure bifurcation of nonuniform solutions from the uniform steady state in a neighborhood of m+ or m− . For some interval | − m± | ¡ 4, we have w( ) which solves w − T (w) = 0. This gives existence of nonuniform solutions. To determine global behavior we need estimates on solution bounds. Lemma 2. Let u and v satisfy (7) and assume that d1 ; d2 ; k; F ¿ 0, then u and ∇u are L2 ( ). Add the two equations in (7) and integrate, we get F u + (F + k) v = F| |:
Assuming F is positive, we have that |u| 6 | | and |v| 6 | |:
Integrating the second equation in (7) generates (F + k) v= uv2 ;
so then
uv2 6 (F + k)| |:
Multiply the 3rst equation in (7) by u and integrate d1 |∇u|2 + u2 v 2 + F u2 = F u 6 F| |:
We then obtain u2 6 | |;
and
|∇u|2 6
uv 6 F| |;
F| | d1
u2 v2 6 F| |:
Lemma 3. Let u and v satisfy (7) and assume that d1 ; d2 ; k; F ¿ 0, then v and ∇v are L2 ( ). Multiply the second equation in (7) by u and integrate −d2 (∇u · ∇v) + u2 v2 − (F + k) uv = 0;
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
117
which gives F| | |∇u · ∇v| 6 : d2 Multiply the 3rst equation in (7) by v and integrate (∇u · ∇v) + uv3 + F uv = F v d1
and we obtain that F| | uv3 6 (d1 + d2 ): d2 Multiply the second equation in (7) by v and integrate F| | (d1 + d2 ): v2 = uv3 6 |∇v|2 + (F + k) d2 d2 Finally we obtain | | v2 6 (d1 + d2 ); d2
|∇v|2 6
F| | (d1 + d2 ): d22
Lemma 4. Let u and v satisfy (7) and assume that d1 ; d2 ; k; F ¿ 0, then Iu is L2 ( ). Again, summing the equations, multiplying by Iu and integrating yields d1 (Iu)2 + d2 (Iu)(Iv) + F |∇u|2 + (F + k) (∇u · ∇v) = 0;
which gives
F| | (Iu)2 6 d2
(Iu)(Iv)
+ : d1 d2 We then estimate
1
2 2
[uv − F(1 − u)][(F + k)v − uv ] : d2 (Iu)(Iv) =
d 1 Expanding this out
1
d2 (Iu)(Iv) = 2 [(F + k)uv3 − u2 v4 − F(F + k)v d1
+ F(F + k)uv + Fuv2 − Fu2 v2 ]
: We may drop the second, third and sixth terms and estimate the remainder by
F(F + k)(1 + F)| | F(F + k)| | + (d1 + d2 ): d2
(Iu)(Iv)
6 d2 d2 d
1
1 2
118
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
Hence the bound is F(F + k)(1 + F)| | F(F + k)| | F| | (Iu)2 6 + (d1 + d2 ) + : 2 2d d1 d2 d d 1 1 2 Lemma 5. Let u and v satisfy (7) and assume that d1 ; d2 ; F ¿ 0, then Iv is L2 ( ). Manipulating the equations in (7) results in [(F + k)v − F(1 − u)]2 : (d1 Iu + d2 Iv)2 =
Notice we have F(F + k)(1 + F)| | F(F + k)| | (Iv)2 6 + (d1 + d2 ) d1 d22 d1 d32 1 + 2 [(F + k)v − F(1 − u)]2 ; d2 which can be used to obtain F(F + k)(1 + F)| | F(F + k)| | + (d1 + d2 ) (Iv)2 6 d1 d22 d1 d32 +
(F + k)2 | | 2F 2 (F + k)| | F 2 | | (d1 + d2 ) + + : 3 d22 d22 d2
What remains is the case where F =0. This can be performed using the same methods as illustrated above to produce 2 d1 |∇u| + u2 v 2 = 0
and
d2
|∇v|2 + k
v2 =
uv3 :
The 3rst equation implies that u must be constant with u = 0 or v = 0. Taking v = 0 satis3es the latter equation as well. Taking u = 0 in the second equation forces v = 0. For F = 0, there is only the uniform solution (u; v) = (c; 0). Lemma 6. Let u and v satisfy (7) with n = 2. Assume that d1 ; d2 ¿ 0, F; k ¿ 0 and that Iu, Iv are L2 ( ), then I2 u, I2 v is L2 ( ). Following the threads developed above, we focus on the 3rst equation in (7) 2 2 (I u) = [(F + v2 )Iu + 4v(∇u∇v) + 2u|∇v|2 + 2uvIv]2 d1
=
[ − (F + v2 )(F(1 − u) − uv2 ) + 4v(∇u∇v)
+ 2u|∇v|2 + 2uv((F + k)v − uv2 )]2 :
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
119
The right-hand side may be expanded and one obtains powers of u, v, ∇u and ∇v. Mixed terms like uv2 or v(∇u∇v) may be separated via HVolder’s inequality. For n = 2, Sobolev embeddings (7.30 in [6]) provides the required bounds on the terms and then the L2 bound on I2 u. This same analysis may be applied to the second equation or we can continue the thread of summing the two equations and manipulating the expression to gain 2 2 2 (d1 I u + d2 I v) = (FIu + (F + k)Iv)2 :
Expanding out the right side we see that we have integrals over (Iu)2 (Iv)2 (Iu)(Iv) which the 3rst two are bounded by assumption and the later is bounded by applying the HVolder’s inequality. Thus we have 2 2 2 2 2 2 d1 (I u) + 2d1 d2 (I u)(I v) + d2 (I2 v)2 ¡ M:
2
2
Since I u is L ( ), the expression above becomes (I2 v)2 ¡ M − 2d1 d2 (I2 u)(I2 v) d22
¡ M + 2d1 d2
(I2 u)2
¡ M + 2d1 d2 M1
1=2
(I2 v)2
(I2 v)2
1=2
1=2 :
This may be rewritten as 1=2 (I2 v)2 − 2d1 d2 M1 (I2 v)2 ¡ M; d22
or
2
(I v)
2
1=2
d22
2
(I v)
2
1=2
− 2d1 d2 M1 ¡ M
which gives the 3nal estimate. Theorem 7. Let u and v satisfy (7) and assume that d1 ; d2 ¿ d ¿ 0, then, for n = 2, u, v are a priori bounded in C 2 ( ) for F; k ¿ 0. Using Lemma 6, I2 u and I2 v are a priori bounded in L2 . For n = 2, the Sobolev embeddings (7.30 in [6]) give that u and v are bounded in C 2 ( ) in terms F; k ¿ 0. Solution continua remain bounded and must branch from and reconnect to the uniform solutions or persist over a 3xed range of the parameters. This provides the analytical basis for the bifurcation of static nonuniform solutions and Hopf bifurcations of nonstatic solutions. Due to the high mode numbers, direct analysis is not tractable and a numerical approach is necessary.
120
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
Fig. 5. Solution u at t = 20; 000 for F = 0:04 and k = 0:06.
6. Numerical results and conclusion Pearson found patterns when the values for the parameters F and k were chosen in a certain range. This range corresponded to the area between the Hopf curve and the static bifurcation curve, which is the crescent region in Fig. 3. For comparison, we selected the same parameter values from Pearson, but used an implicit algorithm to generate numerical results. The patterns we observed compared reasonably well with Pearson’s results. This body of numerical evidence further supports the conjecture that these types of patterns are true nonuniform solutions. The steady-state solution computed by the ADI routine is then passed to an elliptic solver as an initial guess. The output can then be trusted to solve the discrete problem to a residual less than 10−5 . The worm-like patterns persist agreeing with the analytic bifurcation results. The structure of the patterns is very complex. For the case F = 0:04 and k = 0:06 (Fig. 5), integrating the standard initial condition to obtain a steady-state solution and passing this to the elliptic solver gives a good test case. The question of how the symmetry breaking occurs and the resultant anisotropic phenomenon is still open. The Gray–Scott example that Pearson investigated was shown to generate elaborate patterns using implicit numerical methods and standard bifurcation analysis. Future work on these equations is directed at understanding the exact onset of the nonuniform solutions as a function of some bifurcation parameter. It would also be interesting to understand exactly how certain modes establish overall patterns and know how stable these are with respect to perturbations.
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
121
Acknowledgements The authors would like to thank Prof. Hans Othmer for suggesting the problem and many hours of helpful discussions. References [1] M. Ashkenazi, H.G. Othmer, Spatial patterns in coupled biochemical oscillators, J. Math. Biol. 5 (1978) 305–350. [2] J.F.G. Auchmuty, G. Nicolis, Bifurcation analysis of nonlinear reaction–di&usion equations—I. Evolution equations and the steady state solutions, Bull. Math. Biol. 37 (1975) 323–365. [3] K.J. Brown, F.A. Davidson, Global bifurcation in the brusselator system, Nonlinear Anal. 24 (12) (1995) 1713–1725. [4] I.R. Epstein, Complex dynamical behavior in “simple” chemical systems, J. Phys. Chem. 88 (1984) 187–198. [5] R.J. Field, R.M. Noyes, Oscillations in chemical systems, iv, limit cycle behavior in a model of a real chemical reaction, J. Chem. Phys. 60 (5) (1974). [6] D. Gilbarg, N. Trudinger, Elliptic Partial Di&erential Equations of Second Order, Springer, Berlin, 1983. [7] P. Gray, S.K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38 (1) (1983) 29–43. [8] P. Gray, S. K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor: oscillations and instabilities in the system a + 2b → 3b; b → c, Chem. Eng. Sci. 39 (6) (1984) 1087–1097. [9] P. Gray, S.K. Scott, Sustained oscillations and other exotic patterns of behavior in isothermal reactions, J. Phys. Chem. 89 (1985) 22–32. [10] K.J. Lee, W.D. McCormick, Q. Ouyang, H. Swinney, Pattern formation by interacting chemical fronts, Science 261 (1993) 192–194. [11] J.D. Murray, Mathematical Biology, in: Biomathematics, Vol. 19, Springer, New York, 1993. [12] Q. Ouyang, H.L. Swinney, Transition to chemical turbulence, Chaos 1 (4) (1991) 411–420. [13] J.E. Pearson, Complex patterns in a simple system, Science 261 (1993) 189–192. [14] J.E. Pearson, W. Horsthemke, Turing instabilities with nearly equal di&usion coeQcients, J. Chem. Phys. 90 (3) (1989) 1588–1599. [15] I. Prigogene, R. Lefever, Symmetry breaking instabilities in dissipative systems, ii, J. Chem. Phys. 48 (1968) 1665–1700. [16] S.K. Scott, Isolas, mushrooms and oscillations in isothermal, autocatalytic reaction–di&usion equations, Chem. Eng. Sci. 42 (2) (1987) 307–315. [17] L.A. Segel (Ed.), Mathematical Models in Molecular and Cellular Biology, Cambridge University Press, Cambridge, 1980. [18] J.A. Vastano, J.E. Pearson, W. Horsthemke, H. Swinney, Turing patterns in an open reactor, J. Chem. Phys. 88 (10) (1988) 6175–6181. [19] E. Zeidler, Nonlinear Functional Analysis, Vol. IIb, Springer, New York, 1988.
Pattern formation in the Gray–Scott model Je& S. McGough∗ , Kyle Riley Department of Mathematics and Computer Science, South Dakota School of Mines and Technology, 501 E. St. Joseph St., Rapid City, SD 57701, USA Received 23 August 2002; accepted 20 December 2002
Abstract We investigate an elliptic system that arises in cubic autocatalytic reactions known as the Gray–Scott model. Complicated patterns were reported by Pearson in a numerical study of this system. We produce the bifurcation analysis to support the existing numerical evidence for patterns. Speci3cally bifurcation results and C 2 bounds for nonuniform steady states are derived. ? 2003 Elsevier Ltd. All rights reserved. MSC: 35K55; 35K57; 35J60 Keywords: Partial di&erential equations; Reaction di&usion equations; Gradient bounds; Gray–Scott model; Chemical kinetics; Autocatalytic reactions
1. Introduction In the last several decades, chemical kinetics has produced a variety of phenomenon that have translated into challenging mathematical problems. A classical example is seen in the waves of the Belousov–Zhabotinskii reaction. Other examples have been produced that are not quite as complicated and require fewer species interactions, but still yield very interesting behavior. One aspect of these systems is that they do not involve thermal transfer as an essential part of the interaction. Since the 1950s, fundamental studies in reaction kinetics have focused on nonisothermal systems, i.e. where thermal feedback is a critical element. In 1968, Selkov described a particular autocatalytic model of glycolysis. The version of this model due to Gray and Scott [7] is investigated below. Gray–Scott wanted to provide the same foundation for isothermal autocatalytic systems, i.e. chemical feedback. This model becomes the basis of our paper. ∗
Corresponding author. Tel.: +1-605-394-2471; fax: +1-605-394-6078. E-mail address: [email protected] (J.S. McGough).
1468-1218/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S1468-1218(03)00020-8
106
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
In the paper by Pearson [13], very complicated patterns were described for a parabolic system known as the Gray–Scott system. Pearson did a thorough numerical study of the system and found a complex structure in the solutions. However, Pearson used a simple integration scheme and left open the question regarding numerical artifacts. We are able to con3rm the results from Pearson and show that more robust numerical schemes produce the same results. We also provide a bifurcation analysis giving the existence of nonuniform solutions for the steady-state problem, i.e. the elliptic system of Gray–Scott. The Gray–Scott model is ut = d1 Iu + F(1 − u) − uv2 ;
x ∈ ;
vt = d2 Iv − (F + k)v + uv2 ;
x ∈ ;
(1) 2
with boundary data: @u=@ = @v=@ = 0 where = [0; 2:5] , and F; k; d1 ; d2 are constants. Similar systems to the Gray–Scott system have been described in the literature, such as, the Glycolysis model [1,17], or the Brusselator model [2–5,15]. These systems can be collected into a general form ut = d1 Iu + a1 u + b1 v − f(u; v) + g1 (x; u; v);
x ∈ ;
vt = d2 Iv + a2 u + b2 v + f(u; v) + g2 (x; u; v);
x∈
(2)
with zero Neumann boundary data: @u=@ = @v=@ = 0. The parameters for the models are chosen as • Gray–Scott: a1 = −F ¡ 0, b1 = 0, a2 = 0, b2 = −(F + k) ¡ 0, f = uv2 , g1 = F, g2 = 0. • Brusselator: a1 = 0, b1 = B, a2 = 0, b2 = −(B + 1) ¡ 0, f = uv2 , g1 = 0, g2 = A. • Glycolysis: a1 = −k ¡ 0, b1 = 0, a2 = k, b2 = −1, f = uv2 , g1 = , g2 = 0. In a forthcoming paper we develop a general class of problems which cover the three systems above. This general class of problems will share some common properties in their steady-state solutions. The analysis of the Gray–Scott model will lend some insight into the behavior of the general system (2). The next section presents the analysis of the Gray–Scott model. We 3nish with some discussion of numerical results and compare our numerical results with the results given in the Pearson paper. Our results are summarized in the following: Theorem. There exists nonuniform smooth steady solutions to the Gray–Scott system (1) with d1 ; d2 ¿ d ¿ 0, F; k ¿ 0. The Gray–Scott system (1) has smooth bounded solutions and has a rich structure of steady-state solutions. There exists both uniform and nonuniform steady states. 2. Autocatalytic models Gray–Scott began with isothermal autocatalytic systems in the continuously Kowing, well-stirred, tank reactor (CSTR). In this model, isothermal implies the reaction takes
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
107
place under constant temperature, autocatalytic means the catalyst is also the product, and continuously Kowing corresponds to an open system. The well-stirred assumption involves systems that have uniform transport of the reactants. The latter assumption will be relaxed to obtain the general model. There are several well-known examples of autocatalysis as isothermal feedback: Belousov–Zhabotinskii reaction, chorite–iodide– malonic acid reaction, arsenite–iodate reaction, and enzyme systems [4,5,11,12,15]. The reactions listed above are still rather complicated. To reduce the problem for a detailed study, Gray and Scott focused on a model with the overall stoichiometry: A → B, using the reaction rate law: ka (a = [A]). When the reaction is catalyzed by some species Y (the catalyst): A + Y → B + Y, the reaction rate is kay. Our focus is on cubic autocatalysis: A + 2B → 3B, with rate=kab2 . The mass balance for this reaction is da = −k1 ab2 + kf (a0 − a); dt db = k1 ab2 − k2 b + kf (b0 − b): dt
(3)
The system in (3) can be recast into dimensionless quantities. De3ne u as the dimensionless reactant concentration, v the dimensionless catalyst concentration, t the dimensionless time, F a dimensionless “Kow” parameter (1/(residence-time)), k a dimensionless catalyst lifetime and b0 is set to zero. This yields the CSTR model: du = −uv2 + F(1 − u); dt dv (4) = uv2 − (F + k)v: dt One goal of this paper is to analyze the unstirred, or nonuniform, variant of the Gray–Scott model. This nonuniform version of the model will be designated as CFUR (continuously feed, unstirred, reactor). A functional reaction vessel has been built by Harry Swinney et al. [10] and provides some of the 3rst observational evidence of the CFUR in the lab. The speci3c Gray–Scott CFUR Model is the system (1) given in the introduction. In a paper by Pearson [13], numerical studies of the Gray–Scott model produced some interesting patterns which were proposed as similar to those found by Swinney [10]. Pearson used a simple Euler scheme to integrate the parabolic equations. We will later determine some basic bounds on the solutions which generate the complicated patterns. Additional background on the Gray–Scott equations may be found in [4,7– 9,14,16,18]. The patterns observed in these models are thought to be examples of the Turing process. 3. Gray–Scott model dynamics In this section we present an overview of the dynamics found in the Gray–Scott model. Our approach is to examine the local problem (setting the di&usion constants
108
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Fig. 1. Nullclines for the Gray–Scott model.
to zero), determine the steady states, and examine the spatially independent dynamics. We use this to then bootstrap up to study the nonuniform steady states which will be the solution to the related elliptic system. The steady states of the local dynamics (setting d1 and d2 to zero) give us the uniform steady-state solutions. The point (1; 0) is a uniform steady state and will be labelled p1 . If F ¿4(F + k)2 , or the interior of the parabolic region is de3ned by F = 12 [(1=4 − 2k) ± 1=16 − k], then there are two further restpoints which are on the intersection of the nullclines (F + k)v + Fu = F and uv = F + k (see Fig. 1). Solving these two equations gives 1 (F + k)2 1± 1−4 and u1; 2 = 2 F F (F + k)2 v1; 2 = 1∓ 1−4 : 2(F + k) F We will call p2 = (u1 ; v1 ) = (u+ ; v− ) (lower right) and p3 = (u2 ; v2 ) (upper left). This line intersects (1; 0) and (0; F=(F + k)). We will later make use of the fact that v1;2 2 = F is equivalent to F = 4(F + k)2 . This is the set of points for which there are exactly two restpoints for the vector 3eld. One last computation gives that v2 6 F(1 ∓ ) where 0 ¡ ¡ 1 and this bounds v in restpoint p2 : v2 6 F. The intersection of the four half spaces, G1 (u; v) = u ¿ 0; G2 (u; v) = v ¿ 0;
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
109
G3 (u; v) = −(u − c1 ) ¿ 0; G4 (u; v) = −((a1 + a2 )u + (b1 + b2 )v − c2 ) ¿ 0; de3nes a bounded invariant region for c1 , c2 suQciently large. Thus, the local dynamics is bounded in the invariant region de3ned by the intersection of the half planes Gi , i = 1; 2; 3; and 4. Setting the di&usion terms in (1) to zero results in ut = −uv2 + F(1 − u); vt = uv2 − (F + k)v: The intersection of the four half spaces, G1 (u; v) = u ¿ 0; G2 (u; v) = v ¿ 0; G3 (u; v) = −(u − 1) ¿ 0; G4 (u; v) = −(u + v − ) ¿ 0; de3nes a bounded invariant region for ¿ 1 + F=k. The point (1; 0) is always stable. The other two restpoint’s stability are given by the eigenvalues of the following matrix: −(F + v2 ) −2uv −(F + v2 ) −2(F + k) A= = : v2 v2 −(F + k) + 2uv (F + k) 2 The determinant is given by det(A)=(v −F)(F +k) and the trace Tr(A)=−(v2 −k). The eigenvalues are = Tr(A)=2 ± Tr(A)2 − 4 det(A)=2. The point where there is a unique restpoint also yields det(A) = 0. Thus, (F; k) = (1=16; 1=16) yields (u∗ ; v∗ ) = (1=2; 1=4). The restpoints split and travel up/down the line (F + k)v + Fu = F as one enters the region. The restpoint p2 is a saddle point. This follows from recalling that v2 6 F which forces the determinant to be positive and expression under the square root to be larger than Tr(A)2 . Restpoint p3 cannot be a saddle point, but it may have stable, unstable, or possibly spiral, node structure. A Hopf bifurcation occurs if one traverses the correct line in √ the F; k plane. For this, the trace of the linearization must be zero, v = k, and the √ determinant must be nonnegative, v ¿ F. Using the 3rst steady-state equation, one obtains u = F=(k + F):√This value may be plugged into the equation (F + k)v + Fu = F giving (F + k)2 = F k. We may compare this line (the Hopf line) to the lower branch of the steady-state bifurcation line. The Hopf line lies above and intersects the steady-state bifurcation line at (0; 0) and the turning point (1=16; 1=16) (Fig. 2).
110
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121 0.25 Static Bifurcation Stable Nodes Stable Spirals Hopf Bifurcation
0.2
F
0.15
0.1
0.05
0 0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
k
Fig. 2. The stability of the upper-left restpoint, p3 .
4. Stability of uniform solutions To determine if uniform solutions are stable, we study the linearized problem for (2). In other words, we want to see if any modes have exponentially growing terms associated with them. To do this we take the FrechSet derivative of the operator de3ning the Kow, and look for exponentially growing eigenfunctions of the space operator. With reuse of the variables u and v, and taking g1 and g2 to be constant, the linearization of (2) is ut = d1 Iu + (a1 − fu )u + (b1 − fv )v;
x ∈ ;
vt = d2 Iv + (a2 + fu )u + (b2 + fv )v;
x ∈ :
(5)
Separation of the time variable (u = e t , v = e t ) gives d1 I + (a1 − fu − ) + (b1 − fv ) ;
x ∈ ;
d2 I + (a2 + fu ) + (b2 + fv − ) ;
x ∈ :
(6)
To simplify the analysis, let c1 = a1 − fu , c2 = b1 − fv , c3 = a2 + fu , c4 = b2 + fv . De3ne !n to be the eigenfunctions on , and #n ¿ 0 the eigenvalues for I!n = −#n !n :
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
111
We may then split the eigenvector from the eigenfunction: ( ; ) = ($1 ; $2 )!n and thus we reduce to $1 c1 − d1 #n $1 c2 $1 = = ; M $2 c3 c4 − d 2 # n $2 $2 Det(M ) = (d1 + d2 )#n2 − (d2 c1 + d1 c4 )#n + c1 c4 − c2 c3 : Since d1 and d2 are positive, this quadratic is positive for large values of #n . If c1 c4 − c2 c3 ¡ 0 then we get the existence of a positive root. Hence, we get the appearance of nonuniform solutions. This is discussed in greater detail in Section 5. For the Gray–Scott system (1), the associated eigenvalue problem is d1 I − (vs2 + F + ) − 2us vs = 0; d2 I + vs2 + (2us vs − F − k − ) = 0: For the uniform solution (1; 0), this reduces to d1 I − (F + ) = 0; d2 I − (F + k + ) = 0: Following the separation of variables and Fourier analysis outlined above, we now work through the stability analysis for the Gray–Scott system. For (1; 0), we note that the eigenvalues are 1 = −F − d1 #n ¡ 0
and
2 = −F − k − d2 #n ¡ 0
for F ¿ 0. It follows that (1,0) is a stable uniform solution. For the remaining uniform solutions (if they exist), let −2(F + k) −(d1 #n + vs2 + F) M= −d2 #n + F + k vs2 and then the eigenvalues, , may be determined by det(M − I ) = 0: If Tr(M ) = −#n (d1 + d2 ) − vs2 + k and det(M ) = (d1 #n + vs2 + F)(d2 #n − F − k) + 2(F + k)vs2 = d1 d2 #n2 + (d2 (vs2 + F) − d1 (F + k))#n + (F + k)(vs2 − F) then the solutions are Tr(M ) ± Tr(M )2 − 4 det(M ) : = 2 As #n gets large, the structure of M becomes diagonally dominant, which forces the corresponding eigenvalues to be real valued and negative. The term
112
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121 0.25 Static Bifurcation Turing Bifurcation Hopf Bifurcation
0.2
F
0.15
0.1
0.05
0 0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
k
Fig. 3. Location of the unstable modes.
D = Tr(M )2 − 4 det(M ) is strictly increasing in #n (since dD=d#n = 2(d1 − d2 )2 # + constant). Thus, the magnitude of the imaginary component is strictly decreasing. One may conclude that if at # = 0 the eigenvalues are real then they are always real. Dynamical instabilities occur at the lowest modes when the eigenvalues are complex valued at # = 0. A static bifurcation can occur at higher modes, which will occur in our subsequent analysis. The boundary of the region in parameter space where (F; k) : sup R( (F; k; #)) = 0 ; #¿0
de3nes a region of instability that extends the region de3ned by the Hopf bifurcation. This curve is not the same as the curve for Hopf bifurcations. It does intersect the origin, but not the turning point of the static bifurcation curve. One may verify this by plugging into the matrix M and computing the eigenvalues. The 3rst case is = 0 and for the second case ≈ 0:021 for # ≈ 1550. The maximum of may be computed directly. Solving the equations det(M (#)) = 0 and [det(M (#))]# = 0, we obtain from the second equation # = [d1 (F + k) − d2 (v2 + F)]=2d1 d2 . Using the 3rst equation one arrives at [d1 (F + k) − d2 (v2 + F)]2 − 4d1 d2 (F + k)(v2 − F) = 0; which clearly has no solution when v2 ¡ F. First, it is informative to determine where this curve intersects the static bifurcation curve, F − 4(F + k)2 = 0. Since v2 = F, the line k = 0. One would (d1 − 2d2 )F + k = 0, which for the Gray–Scott model is √ also like to determine intersections with the Hopf curve, F k − (F + k)2 = 0. This is
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
113
0.005 Dispersion curves 0
-0.005
R (λ)
-0.01
-0.015
-0.02
-0.025
-0.03 0
1000
2000
3000 µ
4000
5000
6000
0.02 Dispersion curves 0.015 0.01 0.005
R (λ)
0 -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 0
1000
2000
3000
4000
5000
6000
µ
Fig. 4. Top: Family of instability curves for F = 0:038 where k = 0:057 is the lowest curve with increments of Ik = 0:0003. Bottom: Family of instability curves for F = 0:063 and k = 0:06149 (increments of 0:00025).
done numerically for d1 = 2(10−5 ) and d2 = 10−5 giving f = 0:0471, and k = 0:06056. Fig. 3 provides the stability information for d1 = 2(10−5 ) and d2 = 10−5 . Fig. 4 graphs a family of curves for 3xed F and increasing k. The curves are uniformly increasing in k. This gives one typical bifurcation scenario; one that the
114
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
uniform mode becomes unstable 3rst. The x-axis is really a plot of the eigenvalues for the continuous operator; which on the present scale appear dense in the interval. For the second graph in Fig. 4, the uniform mode is not the 3rst mode to become unstable. For the following, we take to be a rectangle with l1 ; l2 to be the length scales in the x; y directions. The eigenvalues for the Laplacian on the square (with zero Neumann data) are k2 m2 #km = + + ; k = 0; 1; 2; : : : ; and m = 0; 1; 2; : : : : l21 l22 For 3xed #, the collection of k and m lie on an ellipse with major and minor axes de3ned by (l1 #=+)2 and (l2 #=+)2 . Pearson studies the case where l1 = l2 = l, and so the curve corresponds to a circle. Because it is symmetric about the m = k line, any eigenvalue #km for which k = m has #km = #mk , in other words, it occurs in pairs. Thus if m = k does not produce the #km , the multiplicity of #km is even. To get eigenvalues of odd multiplicity, we need then that k = m produces an eigenvalue. Restricting, and √ relabelling, to this set in the Pearson example, yields #m = m+ 2=l. 5. A priori bounds for Gray–Scott elliptic problem To prove the existence of nonuniform solutions, a bifurcation approach is used, which is similar to the approach in [3]. There are at most three steady states. The solution (1; 0) is always stable and if it exists, the steady state p2 is unstable for the local dynamics. We focus our attention on the third uniform solution, p3 , which is denoted (us ; vs ) in the following analysis. For the reader’s convenience, system (1) is reproduced here. The functions u and v satisfy d1 Iu + F(1 − u) − uv2 = 0;
x ∈ ;
d2 Iv − (F + k)v + uv2 = 0;
x ∈ ;
(7)
with boundary data: @u=@ = @v=@ = 0 on @ . This system is shifted about the steady state i.e., let u˜ = u − us and v˜ = v − vs and so we obtain d1 Iu˜ − u˜v˜2 − 2u˜vv ˜ s − uv ˜ s2 − us v˜2 − 2us vv ˜ s − F u˜ = 0; 2
d2 Iv˜ + u˜v˜ + 2u˜vv ˜ s+
uv ˜ s2
2
x ∈ ;
+ us v˜ + 2us vv ˜ s − (F + k)v˜ = 0;
x ∈ ;
(8)
with @u=@ = @v=@ = 0 on @ . We drop the tildes for now and recast the problem as an integral equation. Choose Banach spaces X and Y where
@f @g 2; Y = (f; g) | f; g ∈ C ( ); = = 0 on @ and @n @n X = {(f; g) | f; g ∈ C ( )}
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
for some 0 ¡ ¡ 1. Next de3ne operators, L, M , and N , such that 2 −vs − F + d1 L(u; v) = [ − d1 (Iu − u); −d2 (Iv − v)]; M = vs2 and
N (u; v) =
−uv2 − 2uvvs − us v2 uv2 + 2uvvs + us v2
115
−2(F + k)
F + k + d2
:
Thus the elliptic problem in (1) may be re-expressed as L(u; v) = M (u; v) + N (u; v): To apply the standard bifurcation theorems, we must be able to invert the operator L. To do this, one must examine the invertibility of the operator (0 − I ) with zero Neumann boundary data: Iu − u = f;
· ∇u = 0;
x ∈ ; x ∈ @ :
(9)
This is a strictly elliptic operator with c 6 0, we can invoke Theorem 6.31 from [6], U for each f ∈ C ( ). U This proves and 3nd that Eq. (9) above has a unique u ∈ C 2; ( ) that the operator L above is invertible and we obtain (u; v) = L−1 [M (u; v) + N (u; v)] ≡ T (u; v): Hence, T : X → Y is a bounded operator. Y is compactly contained in X so T : X → X is a compact operator. Solving (u; v) = T (u; v) is then equivalent to solving (1). Let K(u; v) = L−1 M (u; v) and G(u; v) = L−1 N (u; v). The associated linear eigenvalue problem is K( ; u; v) − (u; v) = 0 for nontrivial u, v. In this analysis, the parameter represents one of the constants: F, k, d1 , or d2 . Let eigenvalues for 0 − I on with Neumann data be represented by #n and let !m be the associated eigenfunctions. We assume that is a domain for which #n are isolated and that there is an in3nite subset of #n , #nm , for which the algebraic multiplicity of #n is 3nite and odd and that K is Fredholm Index zero. This assumption is valid for the square domain [19]. We will drop the double subscript and refer to this subsequence as #m . This gives rise to the associated algebraic eigenvalue problem, to determine the value of that solves 2 −vs − F + d1 − d1 #m −2(F + k) = 0: det vs2 F + k + d 2 − d2 # m The solutions to the above will be denoted m± . Normally we will choose = F and so m± = Fm± . The nontrivial solution is then K( m± ; c1± !m ; c2± !m ) − (c1± !m ; c2± !m ) = 0: De3ne w = (u; v) and the operator equation is written as w − T (w) = 0:
116
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
Theorem 1. The point m± is a bifurcation point for w − T (w) = 0. Standard bifurcation results assure bifurcation of nonuniform solutions from the uniform steady state in a neighborhood of m+ or m− . For some interval | − m± | ¡ 4, we have w( ) which solves w − T (w) = 0. This gives existence of nonuniform solutions. To determine global behavior we need estimates on solution bounds. Lemma 2. Let u and v satisfy (7) and assume that d1 ; d2 ; k; F ¿ 0, then u and ∇u are L2 ( ). Add the two equations in (7) and integrate, we get F u + (F + k) v = F| |:
Assuming F is positive, we have that |u| 6 | | and |v| 6 | |:
Integrating the second equation in (7) generates (F + k) v= uv2 ;
so then
uv2 6 (F + k)| |:
Multiply the 3rst equation in (7) by u and integrate d1 |∇u|2 + u2 v 2 + F u2 = F u 6 F| |:
We then obtain u2 6 | |;
and
|∇u|2 6
uv 6 F| |;
F| | d1
u2 v2 6 F| |:
Lemma 3. Let u and v satisfy (7) and assume that d1 ; d2 ; k; F ¿ 0, then v and ∇v are L2 ( ). Multiply the second equation in (7) by u and integrate −d2 (∇u · ∇v) + u2 v2 − (F + k) uv = 0;
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
117
which gives F| | |∇u · ∇v| 6 : d2 Multiply the 3rst equation in (7) by v and integrate (∇u · ∇v) + uv3 + F uv = F v d1
and we obtain that F| | uv3 6 (d1 + d2 ): d2 Multiply the second equation in (7) by v and integrate F| | (d1 + d2 ): v2 = uv3 6 |∇v|2 + (F + k) d2 d2 Finally we obtain | | v2 6 (d1 + d2 ); d2
|∇v|2 6
F| | (d1 + d2 ): d22
Lemma 4. Let u and v satisfy (7) and assume that d1 ; d2 ; k; F ¿ 0, then Iu is L2 ( ). Again, summing the equations, multiplying by Iu and integrating yields d1 (Iu)2 + d2 (Iu)(Iv) + F |∇u|2 + (F + k) (∇u · ∇v) = 0;
which gives
F| | (Iu)2 6 d2
(Iu)(Iv)
+ : d1 d2 We then estimate
1
2 2
[uv − F(1 − u)][(F + k)v − uv ] : d2 (Iu)(Iv) =
d 1 Expanding this out
1
d2 (Iu)(Iv) = 2 [(F + k)uv3 − u2 v4 − F(F + k)v d1
+ F(F + k)uv + Fuv2 − Fu2 v2 ]
: We may drop the second, third and sixth terms and estimate the remainder by
F(F + k)(1 + F)| | F(F + k)| | + (d1 + d2 ): d2
(Iu)(Iv)
6 d2 d2 d
1
1 2
118
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
Hence the bound is F(F + k)(1 + F)| | F(F + k)| | F| | (Iu)2 6 + (d1 + d2 ) + : 2 2d d1 d2 d d 1 1 2 Lemma 5. Let u and v satisfy (7) and assume that d1 ; d2 ; F ¿ 0, then Iv is L2 ( ). Manipulating the equations in (7) results in [(F + k)v − F(1 − u)]2 : (d1 Iu + d2 Iv)2 =
Notice we have F(F + k)(1 + F)| | F(F + k)| | (Iv)2 6 + (d1 + d2 ) d1 d22 d1 d32 1 + 2 [(F + k)v − F(1 − u)]2 ; d2 which can be used to obtain F(F + k)(1 + F)| | F(F + k)| | + (d1 + d2 ) (Iv)2 6 d1 d22 d1 d32 +
(F + k)2 | | 2F 2 (F + k)| | F 2 | | (d1 + d2 ) + + : 3 d22 d22 d2
What remains is the case where F =0. This can be performed using the same methods as illustrated above to produce 2 d1 |∇u| + u2 v 2 = 0
and
d2
|∇v|2 + k
v2 =
uv3 :
The 3rst equation implies that u must be constant with u = 0 or v = 0. Taking v = 0 satis3es the latter equation as well. Taking u = 0 in the second equation forces v = 0. For F = 0, there is only the uniform solution (u; v) = (c; 0). Lemma 6. Let u and v satisfy (7) with n = 2. Assume that d1 ; d2 ¿ 0, F; k ¿ 0 and that Iu, Iv are L2 ( ), then I2 u, I2 v is L2 ( ). Following the threads developed above, we focus on the 3rst equation in (7) 2 2 (I u) = [(F + v2 )Iu + 4v(∇u∇v) + 2u|∇v|2 + 2uvIv]2 d1
=
[ − (F + v2 )(F(1 − u) − uv2 ) + 4v(∇u∇v)
+ 2u|∇v|2 + 2uv((F + k)v − uv2 )]2 :
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
119
The right-hand side may be expanded and one obtains powers of u, v, ∇u and ∇v. Mixed terms like uv2 or v(∇u∇v) may be separated via HVolder’s inequality. For n = 2, Sobolev embeddings (7.30 in [6]) provides the required bounds on the terms and then the L2 bound on I2 u. This same analysis may be applied to the second equation or we can continue the thread of summing the two equations and manipulating the expression to gain 2 2 2 (d1 I u + d2 I v) = (FIu + (F + k)Iv)2 :
Expanding out the right side we see that we have integrals over (Iu)2 (Iv)2 (Iu)(Iv) which the 3rst two are bounded by assumption and the later is bounded by applying the HVolder’s inequality. Thus we have 2 2 2 2 2 2 d1 (I u) + 2d1 d2 (I u)(I v) + d2 (I2 v)2 ¡ M:
2
2
Since I u is L ( ), the expression above becomes (I2 v)2 ¡ M − 2d1 d2 (I2 u)(I2 v) d22
¡ M + 2d1 d2
(I2 u)2
¡ M + 2d1 d2 M1
1=2
(I2 v)2
(I2 v)2
1=2
1=2 :
This may be rewritten as 1=2 (I2 v)2 − 2d1 d2 M1 (I2 v)2 ¡ M; d22
or
2
(I v)
2
1=2
d22
2
(I v)
2
1=2
− 2d1 d2 M1 ¡ M
which gives the 3nal estimate. Theorem 7. Let u and v satisfy (7) and assume that d1 ; d2 ¿ d ¿ 0, then, for n = 2, u, v are a priori bounded in C 2 ( ) for F; k ¿ 0. Using Lemma 6, I2 u and I2 v are a priori bounded in L2 . For n = 2, the Sobolev embeddings (7.30 in [6]) give that u and v are bounded in C 2 ( ) in terms F; k ¿ 0. Solution continua remain bounded and must branch from and reconnect to the uniform solutions or persist over a 3xed range of the parameters. This provides the analytical basis for the bifurcation of static nonuniform solutions and Hopf bifurcations of nonstatic solutions. Due to the high mode numbers, direct analysis is not tractable and a numerical approach is necessary.
120
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
Fig. 5. Solution u at t = 20; 000 for F = 0:04 and k = 0:06.
6. Numerical results and conclusion Pearson found patterns when the values for the parameters F and k were chosen in a certain range. This range corresponded to the area between the Hopf curve and the static bifurcation curve, which is the crescent region in Fig. 3. For comparison, we selected the same parameter values from Pearson, but used an implicit algorithm to generate numerical results. The patterns we observed compared reasonably well with Pearson’s results. This body of numerical evidence further supports the conjecture that these types of patterns are true nonuniform solutions. The steady-state solution computed by the ADI routine is then passed to an elliptic solver as an initial guess. The output can then be trusted to solve the discrete problem to a residual less than 10−5 . The worm-like patterns persist agreeing with the analytic bifurcation results. The structure of the patterns is very complex. For the case F = 0:04 and k = 0:06 (Fig. 5), integrating the standard initial condition to obtain a steady-state solution and passing this to the elliptic solver gives a good test case. The question of how the symmetry breaking occurs and the resultant anisotropic phenomenon is still open. The Gray–Scott example that Pearson investigated was shown to generate elaborate patterns using implicit numerical methods and standard bifurcation analysis. Future work on these equations is directed at understanding the exact onset of the nonuniform solutions as a function of some bifurcation parameter. It would also be interesting to understand exactly how certain modes establish overall patterns and know how stable these are with respect to perturbations.
J.S. McGough, K. Riley / Nonlinear Analysis: Real World Applications 5 (2004) 105 – 121
121
Acknowledgements The authors would like to thank Prof. Hans Othmer for suggesting the problem and many hours of helpful discussions. References [1] M. Ashkenazi, H.G. Othmer, Spatial patterns in coupled biochemical oscillators, J. Math. Biol. 5 (1978) 305–350. [2] J.F.G. Auchmuty, G. Nicolis, Bifurcation analysis of nonlinear reaction–di&usion equations—I. Evolution equations and the steady state solutions, Bull. Math. Biol. 37 (1975) 323–365. [3] K.J. Brown, F.A. Davidson, Global bifurcation in the brusselator system, Nonlinear Anal. 24 (12) (1995) 1713–1725. [4] I.R. Epstein, Complex dynamical behavior in “simple” chemical systems, J. Phys. Chem. 88 (1984) 187–198. [5] R.J. Field, R.M. Noyes, Oscillations in chemical systems, iv, limit cycle behavior in a model of a real chemical reaction, J. Chem. Phys. 60 (5) (1974). [6] D. Gilbarg, N. Trudinger, Elliptic Partial Di&erential Equations of Second Order, Springer, Berlin, 1983. [7] P. Gray, S.K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38 (1) (1983) 29–43. [8] P. Gray, S. K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor: oscillations and instabilities in the system a + 2b → 3b; b → c, Chem. Eng. Sci. 39 (6) (1984) 1087–1097. [9] P. Gray, S.K. Scott, Sustained oscillations and other exotic patterns of behavior in isothermal reactions, J. Phys. Chem. 89 (1985) 22–32. [10] K.J. Lee, W.D. McCormick, Q. Ouyang, H. Swinney, Pattern formation by interacting chemical fronts, Science 261 (1993) 192–194. [11] J.D. Murray, Mathematical Biology, in: Biomathematics, Vol. 19, Springer, New York, 1993. [12] Q. Ouyang, H.L. Swinney, Transition to chemical turbulence, Chaos 1 (4) (1991) 411–420. [13] J.E. Pearson, Complex patterns in a simple system, Science 261 (1993) 189–192. [14] J.E. Pearson, W. Horsthemke, Turing instabilities with nearly equal di&usion coeQcients, J. Chem. Phys. 90 (3) (1989) 1588–1599. [15] I. Prigogene, R. Lefever, Symmetry breaking instabilities in dissipative systems, ii, J. Chem. Phys. 48 (1968) 1665–1700. [16] S.K. Scott, Isolas, mushrooms and oscillations in isothermal, autocatalytic reaction–di&usion equations, Chem. Eng. Sci. 42 (2) (1987) 307–315. [17] L.A. Segel (Ed.), Mathematical Models in Molecular and Cellular Biology, Cambridge University Press, Cambridge, 1980. [18] J.A. Vastano, J.E. Pearson, W. Horsthemke, H. Swinney, Turing patterns in an open reactor, J. Chem. Phys. 88 (10) (1988) 6175–6181. [19] E. Zeidler, Nonlinear Functional Analysis, Vol. IIb, Springer, New York, 1988.