Pattern formation in the Gray-Scott model Jeff S. McGough and 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 , TEL: (605) 394-2471, FAX: (605) 394-6078 Contact: J. McGough - [email protected] 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. Specifically bifurcation results and C 2 bounds for non-uniform steady states are derived. Keywords. partial differential equations, reaction diffusion equations, gradient bounds, Gray-Scott model, chemical kinetics, autocatalytic reactions AMS Subject classifications 35K55, 35K57, 35J60 1 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 1950’s, 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. In the paper by Pearson [13], very complicated patterns were described for a parabolic system known as the GrayScott 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 confirm 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 non-uniform solutions for the steady state problem, i.e. the elliptic system of Gray-Scott. The Gray-Scott model is ut = d1 ∆u + F (1 − u) − uv 2 , x ∈ Ω, vt = d2 ∆v − (F + k)v + uv 2 , x ∈ Ω, (1) with boundary data: ∂u/∂ν = ∂v/∂ν = 0 where Ω = [0, 2.5]2 , 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 [15, 5, 2, 4, 3]. These systems can be collected into a general form ut = d1 ∆u + a1 u + b1 v − f (u, v) + g1 (x, u, v), x ∈ Ω, vt = d2 ∆v + 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 = uv 2 , g1 = F , g2 = 0. • Brusselator: a1 = 0, b1 = B, a2 = 0, b2 = −(B + 1) < 0, f = uv 2 , g1 = 0, g2 = A. • Glycolysis: a1 = −k < 0, b1 = 0, a2 = k, b2 = −1, f = uv 2 , 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 finish 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 d 1 , 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 non-uniform steady states. 2 Autocatalytic Models Gray-Scott began with isothermal autocatalytic systems in the continuously flowing, well-stirred, tank reactor (CSTR). In this model, isothermal implies the reaction takes place under constant temperature, autocatalytic means the catalyst is also the product, and continuously flowing 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 1 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 [15, 5, 4, 12, 11]. 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), dt0 db = k1 ab2 − k2 b + kf (b0 − b). dt0 (3) The system in (3) can be recast into dimensionless quantities. Define u as the dimensionless reactant concentration, v the dimensionless catalyst concentration, t the dimensionless time, F a dimensionless “flow” parameter (1/(residencetime)), k a dimensionless catalyst lifetime and b0 is set to zero. This yields the CSTR Model: du = −uv 2 + F (1 − u), dt dv = uv 2 − (F + k)v. dt (4) 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 first observational evidence of the CFUR in the lab. The specific 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, 8, 9, 16, 18, 14]. 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 diffusion constants to zero), determine the steady states, and examine the spatially independent dynamics. We use this to then bootstrap up to study the non-uniform 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 2 point (1, 0) is a uniform and will be i labelled p1 . If F ≥ 4(F + k) , or the interior of the parabolic region h steady state p 1 is defined by F = 2 (1/4 − 2k) ± 1/16 − k , then there are two further restpoints which are on the intersection of the nullclines (F + k)v + F u = F and uv = F + k (see Figure 1). Solving these two equations gives " # " # r r F 1 (F + k)2 (F + k)2 u1,2 = 1± 1−4 , and v1,2 = 1∓ 1−4 . 2 F 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 2 (0, F/(F + k)). We will later make use of the fact that v1,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 field. One last computation givs that v 2 ≤ F (1 ∓ δ) where 0 < δ < 1 and this bounds v in restpoint p2 : v 2 ≤ F . 2 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 Figure 1: Nullclines for the Gray-Scott model. 3 0.8 1 The intersection of the four half spaces, G1 (u, v) = u > 0, G2 (u, v) = v > 0, G3 (u, v) = −(u − c1 ) > 0, G4 (u, v) = −((a1 + a2 )u + (b1 + b2 )v − c2 ) > 0, defines a bounded invariant region for c1 , c2 sufficiently large. Thus, the local dynamics is bounded in the invariant region defined by the intersection of the half planes Gi , i = 1, 2, 3, and 4. Setting the diffusion terms in (1) to zero results in ut = −uv 2 + F (1 − u), vt = uv 2 − (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, defines an 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 + v 2 ) −2(F + k) −(F + v 2 ) −2uv . = A= v2 (F + k) v2 −(F + k) + 2uv The determinant is given by det(A) = (v 2 − F )(F + k) and the trace Tr(A) = −(v 2 − k). The eigenvalues are λ = p 2 Tr(A)/2± Tr(A) − 4 det(A)/2. The point where there is a unique interior restpoint, which gives 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 +F u = F as one enters the region. The restpoint p2 is a saddle point. This follows from recalling that v 2 ≤ 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 non-negative, v ≥ F . Using the first steady state equation, √ one obtains u = F/(k + F ). This value may be plugged into the equation (F + k)v + F u = 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). 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 Frechét derivative of the operator defining the flow, 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 ∆u + (a1 − fu )u + (b1 − fv )v, x ∈ Ω, vt = d2 ∆v + (a2 + fu )u + (b2 + fv )v, x ∈ Ω. (5) Separation of the time variable (u = eλ tα, v = eλ tβ) gives d1 ∆α + (a1 − fu − λ)α + (b1 − fv )β, d2 ∆β + (a2 + fu )α + (b2 + fv − λ)β, 4 x ∈ Ω, x ∈ Ω. (6) To simplify the analysis, let c1 = a1 − fu , c2 = b1 − fv , c3 = a2 + fu , c4 = b2 + fv . Define φn to be the eigenfunctions on Ω, and µn ≥ 0 the eigenvalues for ∆φn = −µn φn . We may then split the eigenvector from the eigenfunction: (α, β) = (Θ1 , Θ2 )φn and thus we reduce to Θ1 c 1 − d 1 µn c2 Θ1 Θ1 M = =λ Θ2 c3 c 4 − d 2 µn Θ2 Θ2 Det(M ) = (d1 + d2 )µ2n − (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 ∆α − (vs2 + F + λ)α − 2us vs β = 0 d2 ∆β + vs2 α + (2us vs − F − k − λ)β = 0. For the uniform solution (1, 0), this reduces to d1 ∆α − (F + λ)α = 0 d2 ∆β − (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 −(d1 µn + vs2 + F ) −2(F + k) M= , vs2 −d2 µn + F + k 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 µ2n + (d2 (vs2 + F ) − d1 (F + k))µn + (F + k)(vs2 − F ) then the solutions are p Tr(M )2 − 4det(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 D = Tr(M )2 − 4det(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. λ= Tr(M ) ± 5 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 <(λ(F, k, µ)) = 0}, µ≥0 defines a region of instability that extends the region defined 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 first case is λ = 0 and for the second case λ ≈ .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 (v 2 + F )]/2d1 d2 . Using the first equation one arrives at [d1 (F + k) − d2 (v 2 + F )]2 − 4d1 d2 (F + k)(v 2 − F ) = 0, which clearly has no solution when v 2 < F . First, it is informative to determine where this curve intersects the static bifurcation curve, F − 4(F + k)2 = 0. Since v 2 = F , (d1 − 2d2 )F + k = 0, which√for the Gray-Scott model is the line k = 0. One would also like to determine intersections with the Hopf curve, F k − (F + k)2 = 0. This is done numerically for d1 = 2(10−5 ) and d2 = 10−5 giving f = 0.0471, and k = 0.06056. Figure 3 provides the stability information for d1 = 2(10−5 ) and d2 = 10−5 . Figure 4 graphs a family of curves for fixed F and increasing k. The curves are uniformly increasing in k. This gives one typical bifurcation scenario; one that the uniform mode becomes unstable first. The x axis is really a plot of the eigenvalues for the continuous operator; which on the present scale, they appear dense in the interval. For the second graph in Figure 4, the uniform mode is not the first 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 s k2 m2 µkm = π + 2 , k = 0, 1, 2, . . . , m = 0, 1, 2, . . . 2 l1 l2 For fixed µ, the collection of k and m lie on an ellipse with major and minor axes defined 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 6= 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. 6 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 k 0.04 0.05 Figure 2: The stability of the upper-left restpoint, p3 . 7 0.06 0.07 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 k 0.04 Figure 3: Location of the unstable modes. 8 0.05 0.06 0.07 0.005 Dispersion curves 0 -0.005 ℜ(λ) -0.01 -0.015 -0.02 -0.025 -0.03 0 1000 2000 3000 µ 0.02 4000 5000 6000 Dispersion curves 0.015 0.01 0.005 ℜ(λ) 0 -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 0 1000 2000 9 3000 µ 4000 5000 6000 Figure 4: Top: Family of instability curves for F = .038 where k = .057 is the lowest curve with increments of ∆k = .0003. Bottom: Family of instability curves for F = .063 and k = .06149 (increments of .00025). 5 Apriori 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 ∆u + F (1 − u) − uv 2 = 0, x ∈ Ω, d2 ∆v − (F + k)v + uv 2 = 0, x ∈ Ω, (7) with boundary data: ∂u/∂ν = ∂v/∂ν = 0 on ∂Ω. This system is shifted about the steady state i.e., let ũ = u − us and ṽ = v − vs and so we obtain d1 ∆ũ − ũṽ 2 − 2 ũṽvs − ũvs 2 − us ṽ 2 − 2 us ṽvs − F ũ = 0, x ∈ Ω, d2 ∆ṽ + ũṽ 2 + 2 ũṽvs + ũvs 2 + us ṽ 2 + 2 us ṽvs − (F + k)ṽ = 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 X = {(f, g) | f, g ∈ C α (Ω)} , ∂n ∂n for some 0 < α < 1. Next define operators, L, M , and N , such that −vs2 − F + d1 L(u, v) = [−d1 (∆u − u), −d2 (∆v − v)] , M = vs2 and N (u, v) = −uv 2 − 2 uvvs − us v 2 uv 2 + 2 uvvs + us v 2 −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 (∆ − I) with zero Neumann boundary data: ∆u − u = f, x ∈ Ω ν · ∇u = 0, x ∈ ∂Ω. (9) This is a strictly elliptic operator with c ≤ 0, we can invoke Theorem 6.31 from [6], and find that equation (9) above has a unique u ∈ C 2,α (Ω) for each f ∈ C α (Ω). This proves 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, d 1 , or d2 . Let eigenvalues for ∆ − 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 infinite subset of µn , µnm , for which the algebraic 10 multiplicity of µn is finite 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 −vs2 − F + d1 − d1 µm −2(F + k) = 0. det vs2 F + k + d 2 − d 2 µm ± ± λ 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. Define w = (u, v) and the operator equation is written as w − T (w) = 0. 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 neighbor− ± hood of λ+ m or λm . For some interval |λ − λm | < , 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 Z Z v = F |Ω|. u + (F + k) F Ω Ω Assuming F is positive, we have that Z Ω |u| ≤ |Ω|, Z and Ω |v| ≤ |Ω|. Integrating the second equation in (7) generates (F + k) so then Z Ω Z v= Ω Z uv 2 , Ω uv 2 ≤ (F + k)|Ω|. Multiply the first equation in (7) by u and integrate Z Z Z Z |∇u|2 + u2 v 2 + F u2 = F u ≤ F |Ω|. d1 Ω From this we obtain and Ω Z Z Ω Ω u2 ≤ |Ω|, uv ≤ F |Ω|, Ω Z Ω Z |∇u|2 ≤ Ω Ω F |Ω| , d1 u2 v 2 ≤ F |Ω|. Lemma 3 Let u and v satisfy (7) and assume that d1 , d2 , k, F > 0, then v and ∇v are L2 (Ω). 11 Multiply the second equation in (7) by u and integrate Z Z Z uv = 0, u2 v 2 − (F + k) −d2 (∇u · ∇v) + Ω Ω Ω which gives Z Ω |∇u · ∇v| ≤ F |Ω| . d2 Multiply the first equation in (7) by v and integrate Z Z Z Z 3 v, uv = F uv + F d1 (∇u · ∇v) + and we obtain that Z Ω uv 3 ≤ Ω Ω Ω Ω F |Ω| (d1 + d2 ) . d2 Multiply the second equation in (7) by v and integrate Z Z Z F |Ω| (d1 + d2 ) . d2 |∇v|2 + (F + k) v2 = uv 3 ≤ d2 Ω Ω Ω Finally we obtain Z Ω v2 ≤ |Ω| (d1 + d2 ) , d2 Z Ω |∇v|2 ≤ F |Ω| (d1 + d2 ) . d22 Lemma 4 Let u and v satisfy (7) and assume that d1 , d2 , k, F > 0, then ∆u is L2 (Ω). Again, summing the equations, multiplying by ∆u and integrating yields Z Z Z Z |∇u|2 + (F + k) (∇u · ∇v) = 0, d1 (∆u)2 + d2 (∆u)(∆v) + F which gives d1 We then estimate Ω Ω Ω Z Ω Z F |Ω| (∆u) ≤ d2 (∆u)(∆v) + . d2 Ω Ω 2 Z Z 1 2 2 uv − F (1 − u) (F + k)v − uv d2 (∆u)(∆v) = . d1 Ω Ω Expanding this out Z Z 1 3 2 4 2 2 2 d2 (∆u)(∆v) = 2 (F + k)uv − u v − F (F + k)v + F (F + k)uv + F uv − F u v . d1 Ω Ω We may drop the second, third and sixth terms and estimate the remainder by Z F (F + k)(1 + F )|Ω| F (F + k)|Ω| d2 (∆u)(∆v) ≤ + (d1 + d2 ). d21 d21 d2 Ω Hence the bound is Z Ω (∆u)2 ≤ F |Ω| F (F + k)(1 + F )|Ω| F (F + k)|Ω| + (d1 + d2 ) + . d21 d21 d2 d1 d2 Lemma 5 Let u and v satisfy (7) and assume that d1 , d2 , F > 0, then ∆v is L2 (Ω). 12 Manipulating the equations in (7) results in Z Z 2 2 [(F + k)v − F (1 − u)] . (d1 ∆u + d2 ∆v) = Ω Ω Notice we have Z F (F + k)(1 + F )|Ω| F (F + k)|Ω| 1 (∆v) ≤ + (d1 + d2 ) + 2 2 3 d d d d d 1 2 1 2 Ω 2 2 Z Ω [(F + k)v − F (1 − u)]2 , which can be used to obtain Z Ω F (F + k)(1 + F )|Ω| F (F + k)|Ω| + (d1 + d2 ) d1 d22 d1 d32 2 2 2F (F + k)|Ω| F 2 |Ω| (F + k) |Ω| (d1 + d2 ) + + . + 3 d2 d22 d22 (∆v)2 ≤ What remains is the case where F = 0. This can be performed using the same methods as illustrated above to produce Z Z d1 u2 v 2 = 0, |∇u|2 + Ω Ω and d2 Z 2 Ω |∇v| + k Z 2 v = Ω Z uv 3 . Ω The first equation implies that u must be constant with u = 0 or v = 0. Taking v = 0 satisfies 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 ∆u, ∆v are L2 (Ω), then ∆2 u, ∆2 v is L2 (Ω). Following the threads developed above, we focus on the first equation in (7) Z Z 2 2 d1 (∆ u) = [(F + v 2 )∆u + 4v(∇u∇v) + 2u|∇v|2 + 2uv∆v]2 Ω Ω = Z Ω [−(F + v 2 )(F (1 − u) − uv 2 ) + 4v(∇u∇v) + 2u|∇v|2 + 2uv((F + k)v − uv 2 )]2 . The right hand side may be expanded and one obtains powers of u, v, ∇u and ∇v. Mixed terms like uv 2 or v(∇u∇v) may be separated via Hölder’s inequality. For n = 2, Sobolev embeddings (7.30 in [6]) provides the required bounds on the terms and then the L2 bound on ∆2 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 Z Z (d1 ∆2 u + d2 ∆2 v)2 = (F ∆u + (F + k)∆v)2 . Ω Ω Expanding out the right side we see that we have integrals over (∆u)2 (∆v)2 (∆u)(∆v) which the first two are bounded by assumption and the later is bounded by applying the Hölder’s inequality. Thus we have Z Z Z d21 (∆2 u)2 + 2d1 d2 (∆2 u)(∆2 v) + d22 (∆2 v)2 < M. Ω Ω 2 Ω 2 Since ∆ u is L (Ω), the expression above becomes Z Z d22 (∆2 v)2 < M − 2d1 d2 (∆2 u)(∆2 v) Ω Ω 13 < M + 2d1 d2 Z (∆2 u)2 Ω 1/2 Z (∆2 v)2 Ω 1/2 < M + 2d1 d2 M1 Z (∆2 v)2 Ω 1/2 . This may be rewritten as d22 or Z Z 2 Ω (∆2 v)2 − 2d1 d2 M1 (∆ v) Ω 2 1/2 " d22 Z Z 2 (∆ v) Ω 2 (∆2 v)2 Ω 1/2 1/2 < M, # − 2d1 d2 M1 < M which gives the final estimate. Theorem 7 Let u and v satisfy (7) and assume that d1 , d2 > d > 0, then, for n = 2, u, v are apriori bounded in C 2 (Ω) for F, k ≥ 0. Using Lemma 6, ∆2 u and ∆2 v are apriori 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 fixed range of the parameters. This provides the analytical basis for the bifurcation of static non-uniform solutions and Hopf bifurcations of non-static solutions. Due to the high mode numbers, direct analysis is not tractable and a numerical approach is necessary. 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 Figure 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 non-uniform 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 = .04 and k = .06 (Figure 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. Acknowledgments. The authors would like to thank Prof. Hans Othmer for suggesting the problem and many hours of helpful discussions. References [1] M. Ashkenazi and H.G. Othmer. Spatial patterns in coupled biochemical oscillators. Jour. Math. Biol., 5:305– 350, 1978. [2] J. F. G. Auchmuty and G. Nicolis. Bifurcation analysis of nonlinear reaction-diffusion equations. - I. Evolution equations and the steady state solutions. Bull. Math. Biol., 37:323–365, 1975. 14 [3] K. J. Brown and F. A. Davidson. Global bifurcation in the brusselator system. Nonlinear Analysis, 24(12):1713– 1725, 1995. [4] Irving R. Epstein. Complex dynamical behavior in “simple” chemical systems. J. Phys. Chem., 88:187–198, 1984. [5] R.J. Field and R.M. Noyes. Oscillations in chemical systems. iv. limit cycle behavior in a model of a real chemical reaction. Jour. Chem. Phys., 60(5), 1974. [6] D. Gilbarg and N. Trudinger. Elliptic Partial Differential Equations of Second Order. Springer Verlag, Berlin, 1983. [7] P. Gray and S. K. Scott. Autocatalytic reactions in the isothermal continuous stirred tank reactor: isolas and other forms of multistability. Chem. Eng. Sci., 38(1):29–43, 1983. [8] P. Gray and 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):1087–1097, 1984. [9] P. Gray and S. K. Scott. Sustained oscillations and other exotic patterns of behavior in isothermal reactions. J. Phys. Chem., 89:22–32, 1985. [10] K. J. Lee, W. D. McCormick, Q. Ouyang, and H. Swinney. Pattern formation by interacting chemical fronts. Science, 261:192–194, 1993. [11] J. D. Murray. Mathematical Biology, volume 19 of Biomathematics. Springer-Verlag, New York, 1993. [12] Q. Ouyang and H.L. Swinney. Transition to chemical turbulence. Chaos, 1(4):411–420, 1991. [13] J. E. Pearson. Complex patterns in a simple system. Science, 261:189–192, 1993. [14] J. E. Pearson and W. Horsthemke. Turing instabilities with nearly equal diffusion coefficients. J. Chem. Phys., 90(3):1588–1599, 1989. [15] I. Prigogene and R. Lefever. Symmetry breaking instabilities in dissipative systems. ii. J. Chem. Phys, 48:1665– 1700, 1968. [16] S. K. Scott. Isolas, mushrooms and oscillations in isothermal, autocatalytic reaction-diffusion equations. Chem. Eng. Sci., 42(2):307–315, 1987. [17] L.A. Segel, editor. Mathematical Models in Molecular and Cellular Biology. Cambridge University Press, 1980. [18] J. A. Vastano, J. E. Pearson, W. Horsthemke, and H. Swinney. Turing patterns in an open reactor. J. Chem. Phys., 88(10):6175–6181, 1988. [19] E. Zeidler. Nonlinear Functional Analysis, volume IIb. Springer Verlag, New York, 1988. 15 Figure 5: Solution u at t = 20000 for F = .04 and k = .06. 16 Figure captions: 1. Nullclines for the Gray-Scott model. 2. The stability of the upper-left restpoint, p3 . 3. Location of the unstable modes. 4. Top: Family of instability curves for F = .038 where k = .057 is the lowest curve with increments of ∆k = .0003. Bottom: Family of instability curves for F = .063 and k = .06149 (increments of .00025). 5. Solution u at t = 20000 for F = .04 and k = .06. 17

© Copyright 2018