1. No category

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 - jeff.mcgough@sdsmt.edu
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