Stability and bifurcation of nonconstant solutions to a reaction-diffusion system with conservation of a mass Yoshihisa Morita∗ and Toshiyuki Ogawa†
Abstract We are dealing with a two-component system of reaction-diffusion equations with conservation of a mass. This system is proposed as a conceptual model for cell polarity. Since the system has the conservation of a mass, the steady state problem is reduced to that of a scalar reactiondiffusion equation with a nonlocal term. In particular, we consider the case when the reaction term is linear in one variable. Then the equations are transformed into a similar model of the phase field proposed by Fix and Caginalp, thus the equations allows a Lyapunov function. By investigating the linearized stability of a nonconstant equilibrium solution, we prove that the stability of the equilibrium solution is closely related to that of the nonlocal scalar reaction-diffusion equation. We also exhibit global bifurcation diagrams for equilibrium solutions to specific model equations by numerics together with a normal form near a bifurcation point.
∗
Department of Applied Mathematics and Informatics, Ryukoku University, Seta Otsu 520-2194, Japan, E-mail:
[email protected] † Division of Mathematical Science, Graduate School of Engineering Science, Osaka University, Machikaneyama, Toyonaka, Osaka 560-8531, Japan,
[email protected] AMS 2000 Subject Classification: 35B32, 35B35, 35B40, 35K57. Key words: reaction-diffusion system, conservation of a mass, equilibrium solution, stability, bifurcation.
1
1
Introduction
There are many mathematical models describing pattern formations in nonlinear science. Among other things models of reaction-diffusion equations are very useful for describing the emergence of spatial structures in not only physical and chemical systems but also biological systems. Those have been extensively developed since the discovery of the Turing instability [18] (see also [6], [4], [12], [15], [14], [1] and references therein). In cellular process we can observe molecular accumulations such as chemotaxis and cell polarity. Recently, Otsuji et al [11] proposed a mass-conserved reaction-diffusion system which exhibits the properties of cell polarity. IshiharaOtsuji-Mochizuki [10] also show interesting dynamics of solutions to a similar mass-conserved system together with some formal mathematical analysis. Motivated by their works [10] and [11], we consider the following 2-component reaction-diffusion system: { ut = d1 ∆u + k1 f (u, v), (1.1) vt = d2 ∆v − k2 f (u, v), in a bounded domain, where dj , kj (j = 1, 2) are positive parameters and ∆ stands for the Laplacian. After rescaling the variables as √ k1 t → t, k2 /d2 x → x and putting D := d1 k2 /d2 k1 , we obtain
{ ut = D∆u + f (u, v), τ vt = ∆v − f (u, v)
τ := k1 /k2 (x ∈ Ω),
(1.2)
(1.3)
where Ω is a bounded domain in Rn with smooth boundaries ∂Ω. Henceforth we consider (1.3) with the Neumann boundary conditions ∂u = 0, ∂ν
∂v =0 ∂ν
(x ∈ ∂Ω),
or in Ω = [0, L] with the periodic boundary conditions { u(0, t) = u(L, t), ux (0, t) = ux (L, t), v(0, t) = v(L, t), vx (0, t) = vx (L, t).
(1.4)
(1.5)
The aim of the present article is providing a stability theorem for equilibrium solutions to (1.3) together with global bifurcation diagrams for specific functions of f by numerics. 2
Throughout the paper, we assume f is sufficiently smooth and the zero set {(u, v) : f (u, v) = 0} is nonempty. As for specific cases, we are interested in f (u, v) = − [ f (u, v) = −a1
au + v, u2 + b
(1.6)
] u+v −v , {a2 (u + v) + 1}2
(1.7)
and f (u, v) = a1 (u + v)[(αu + v)(u + v) − a2 ],
(1.8)
where all the coefficients are positive. The case (1.6) is treated in [10] while (1.7) and (1.8) are in [11]. We observe some basic properties for solutions to (1.3). Firstly, by applying the theorem in Henry[9], the equations (1.3) with (1.4) or (1.5) allows a unique time local solution for the initial conditions (u(x, 0), v(x, 0)) = (u0 (x), v0 (x)) ∈ X := H 1 (Ω) × H 1 (Ω),
(1.9)
provided that ∥f (u(·), v(·))∥ < ∞ for (u(·), v(·)) ∈ X, where ∥ · ∥ stands for the L2 -norm. We denote the solution with (1.9) by (u(·, t; u0 , v0 ), v(·, t; u0 , v0 )). In addition, if |fu |+|fv | is uniformly bounded, then the solution map (u0 (·), v0 (·)) 7→ (u(·, t; u0 , v0 ), v(·, t; u0 , v0 )) for t ≥ 0 generates a C 1 semiflow in X. Namely, {S(t)}t≥0 defined by S(t)u0 := (u(·, t; u0 ), v(·, t; u0 )),
u0 := (u0 , v0 )
(1.10)
is a C 1 semiflow in X (see Appendix in this paper). Secondary, we notice conservation of a mass such that the solution (u, v) = (u(x, t), v(x, t)) to (1.3) with the boundary conditions satisfies ∫ ∫ d (u + τ v)dx = (D∆u + ∆v)dx = 0, dt Ω Ω that is, the quantity 1 s := |Ω|
∫ (u0 + τ v0 )dx
(1.11)
Ω
is conserved for any t ≥ 0 as long as the solution is defined. Thirdly, we introduce a new variable w := Du + v. Then the mass conserved system (1.3) is transformed into { ut = D∆u + f (u, w − Du), (x ∈ Ω), τ wt + (1 − τ D)ut = ∆w 3
(1.12)
(1.13)
with or
∂u = 0, ∂ν
∂u =0 ∂ν
(x ∈ ∂Ω),
{ u(0, t) = u(L, t), ux (0, t) = ux (L, t), w(0, t) = w(L, t), wx (0, t) = wx (L, t).
We note that for the case τ D > 1 we can transform (1.3) as { τ vt = ∆v − f ((w − v)/D, v), (x ∈ Ω), wt + (τ D − 1)vt = D∆w
(1.14)
(1.15)
(1.16)
which is the same form as (1.13) if interchanging the role of u and v. Remark 1.1 For a specific f such as f (u, v) = u − u3 + kv,
(1.17)
we notice that (1.13) is a similar to the model equations by Fix [5] and Caginalp [2] if τ1 D < 1. In fact, the phase field model found in [5] and [2] is given as { τ1 φt = ϵ2 ∆φ + W (φ) + cT, (1.18) Tt + 2ℓ φt = κ∆T, where W (φ) = φ − φ3 . Conversely, we can transform (1.18) into (1.1) in what follows. For the case of τ1 κ > ϵ2 , we use )( ) ) ( ( 1 0 φ u = , ϵ2 ℓ 1 T v 2(τ1 κ−ϵ2 ) to obtain
{ 2 ut = τϵ1 ∆u + τ11 f (u, v), f (u, v), vt = κ∆v − 2(τ1ℓκ κ−ϵ2 )
where f (u, v) := u(1 − u2 ) +
cϵ2 ℓ u + cv. 2(τ1 κ − ϵ2 )
On the other hand, if τ1 κ < ϵ2 , we can use )( ) ( ) ( 1 0 φ u = . ϵ2 ℓ 1 T v 2(ϵ2 −τ1 κ) to obtain the similar form.
4
(1.19)
Fourthly, consider the stationary problem of (1.3) { D∆u + f (u, v) = 0, (x ∈ Ω), ∆v − f (u, v) = 0
(1.20)
with the Neumann or the periodic boundary conditions. Let (u∗ (x; s), v ∗ (x; s)) be a solution to (1.20) satisfying ∫ 1 s= (u∗ + τ v ∗ )dx. (1.21) |Ω| Ω Then u∗ is a solution to ( ) ∫ −1 D∆u + f u, −Du + s/τ − (1/τ − D)|Ω| udx = 0,
(1.22)
Ω
with the boundary conditions. In fact, from (1.20) it follows that Du∗ + v ∗ = C
(1.23)
holds for a constant C. Integrating (1.23) and using (1.21), we easily obtain ∫ ∗ ∗ −1 v = −Du + s/τ − (1/τ − D)|Ω| u∗ dx. (1.24) Ω
Thus, a solution of (1.22) together with (1.24) gives a solution (u∗ , v ∗ ) of (1.20), and vice versa. Now we state mathematical results on the solutions to the mass conserved reaction-diffusion system (1.3). By virtue of the studies by Miyasita [13] and Suzuki-Tasaki [17] for the phase field model (see also [16]), if f (u, v) = f0 (u) + kv (k > 0), we notice } ∫ { D kD 2 τk 2 2 |∇u| − F0 (u) + u + (Du + v) dx (1.25) E(u, v) := 2 2 2(1 − τ D) Ω ∫ gives a Lyapunov function, where F0 (u) = f0 (u)du. More precisely, we have the following lemma: Lemma 1.1 Let ∫
f (u, v) = f0 (u) + kv,
(1.26)
and F0 (u) = f0 (u)du, where k is positive constant. Assume −F0 (u) is bounded below and the simiflow S(t) defined by (1.10) is C 1 in X. If 0 < τ D < 1, then E(S(t)u0 ) is non-increasing in t and E(S(t)u0 ) = constant (∀t ∈ R) implies that u0 is an equilibrium. Moreover, the omega limit set of any bounded orbit consists of equilibria, where the bounded orbit implies that S(t)u0 is uniformly bounded above for t ≥ 0. 5
Next we consider the stability of constant equilibrium solutions. Let (u, v) be a constant equilibrium solution to (1.3) and consider the linearized eigenvalue problem { D∆ϕ + fu (u, v)ϕ + fv (u, v)ψ = −λϕ, (1.27) ∆ψ − fu (u, v)ϕ − fv (u, v)ψ = −τ λψ, with the Neumann boundary conditions or the periodic boundary conditions. We let {σj }j≥1 be the eigenvalues of −∆ with (1.4) or (1.5). Those are labeled in increasing order and counting multiplicity. Let {φj }j≥1 be the corresponding normalized eigenfunctions, that is, ∫ ⟨φj , φm ⟩L2 := φj (x)φm (x)dx = δjm . Ω
We note that σ1 = 0. The next lemma is due to [10] and [11]. Lemma 1.2 Assume 0 < fv (u, v) < fu (u, v)/D < fv (u, v)/τ D.
(1.28)
Then for σj (j ≥ 2) satisfying σj + fv (u, v) − fu (u, v)/D < 0, the condition (1.28) allows an eigenvalue λ = λj < 0 of (1.27), where λj is obtained by a negative eigenvalue of −A0 ) ( fv (u, v) −Dσj + fu (u, v) . (1.29) A0 := −fu (u, v)/τ −(σj + fv (u, v))/τ A corresponding eigenfunction to λj is given by ( ) ( ) ϕ ξ = φj (x) , ψ η
(1.30)
where (ξ, η)T is an eigenvector of −A0 to λj . We notice that the constant solution is stable under the homogeneous perturbation provided that (1.28) holds. Indeed, for σ1 = 0, the matrix −A0 has eigenvalues λ = 0, λ = fv (u, v)/τ − fu (u, v) > 0. This zero eigenvalue is ruled out if the mass s is fixed. Therefore, the constant solution is asymptotically stable in the diffusion-free equations under the constraint of the mass. The above lemma implies that the uniform state looses the stability as the diffusion constant D decreases. Hence, we may say that the Turing instability takes places in this system. We also note that this Turing instability is realized for the specific f of (1.6), (1.7) and (1.8) if we choose the parameters appropriately. 6
We consider the next linearized eigenvalue problem for a nonconstant solution (u∗ (·; s), v ∗ (·; s)) to (1.20) ( ) ( ) ( ) ϕ D∆ϕ + fu (u∗ , v ∗ )ϕ + fv (u∗ , v ∗ )ψ ϕ L := − =λ (1.31) ψ {∆ψ − fu (u∗ , v ∗ )ϕ − fv (u∗ , v ∗ )ψ}/τ ψ with the domain DN (L) := DN × DN , DN := {φ ∈ H 2 (Ω) : ∂φ/∂ν = 0 (x ∈ ∂Ω)},
(1.32) (1.33)
or Dp (L) := Dp × Dp , (1.34) 2 Dp := {φ ∈ H (0, L) : φ(0) = φ(L), φx (0) = φx (L) = 0}. (1.35) We note that for the periodic case, ( ∗ ) ( ) ux 0 L = ∗ vx 0 holds. If a nonconstant solution (u∗ (·; s), v ∗ (·; s)) is continuously differentiable with respect to s, then (u∗s , vs∗ ) = (∂u∗ /∂s, ∂v ∗ /∂s) satisfies ∫ (u∗s + τ vs∗ )dx = 1. Ω
Corresponding to (1.22), we also introduce the linearized eigenvalue problem L0 [φ] := −{D∆φ + (fu (u∗ , v ∗ ) − Dfv (u∗ , v ∗ ))φ ∫ ∗ ∗ −1 −(1/τ − D)fv (u , v )|Ω| φdx} = µφ,
(1.36)
Ω
with the domain D(L0 ) = DN or Dp . We see L0 [u∗x ] = 0 in the periodic case. Let w(x) be a solution of k (1.37) L0 [w] = τ with the boundary conditions. The solvability condition for this equation is 1 ∈ R(L0 ) := {g = L0 [ϕ] : ϕ ∈ DK (K = N or p)} If L0 is invertible, then w is uniquely determined. The next theorem is the main result of this article. Theorem 1.3 Let f0 (u) be a C 1 function. Assume (1.26) with k > 0 and 0 < τ D < 1. Let (u∗ (x; s0 ), v ∗ (x; s0 )) be a non-constant solution to (1.20) with the Neumann boundary conditions or the periodic conditions, where s = s0 in (1.21). Then all the eigenvalues of L of (1.31) in the domain (1.32) or (1.34) are real. Moreover, the following holds: 7
(i) If every eigenvalue of L0 of (1.36) in DN is positive, then all the eigenvalues except for zero of L is positive. In addition zero eigenvalue is simple with a corresponding eigenfunction (ϕ, ψ) = (w(x), cD − Dw(x)), where ∫ 1 1 − τD cD := − wdx (1.38) τ τ |Ω| Ω (ii) If zero eigenvalue of L0 of (1.36) in Dp is simple and the remaining eigenvalues are positive, then all the eigenvalues except for zero of L are positive. In addition the eigenspace of zero eigenvalue is the linear hull of (u∗x (·; s0 ), vx∗ (·; s0 )) and (w(·), cD − Dw(·)), where cD is given by (1.38). For the result of Theorem 1.3 (i), we can assert that u∗ (·; s) is continuously differentiable with s in a neighborhood of s0 . As a consequence w(·) = u∗s (·; s0 ) holds. By virtue of Theorem 1.3 we obtain the following dynamical stability of an equilibrium solution: Corollary 1.4 In addition to the assumptions in Theorem 1.3, suppose that the semiflow S(t) defined by (1.10) is C 1 in X = H 1 (Ω) × H 1 (Ω). Let u∗ (·; s0 ) := (u∗ (·; s0 ), v ∗ (·; s0 )). Then the following holds: (i) The case of the Neumann boundary conditions: There exist δ > 0, β > 0 ˜ 0 = (˜ and M > 0 such that if u u0 , v˜0 ) ∈ X satisfies ∥˜ u0 − u∗ (·; s0 )∥X < δ, then the solution (u, v) = S(t)˜ u0 to (1.3)-(1.4) satisfies ∥S(t)˜ u0 − u∗ (·; s˜)∥X ≤ M e−βt ∥˜ u0 − u∗ (·; s˜)∥X , where s˜ is given by 1 s˜ = |Ω|
∫ (˜ u0 (x) + τ v˜0 (x))dx,
(1.39)
Ω
and |˜ s − s0 | → 0 as δ → 0. (ii) The case of the periodic boundary conditions: Let r0 be a positive number. Suppose that for each s ∈ Is := (s0 − r0 , s0 + r0 ), the same condition for L0 as in Theorem 1.3 (ii) holds. Then there exist δ1 > 0 and β1 > 0 such that if ∥˜ u0 − u∗ (·; s0 )∥X < δ1 , then the solution (u, v) = S(t)˜ u0 to (1.3)-(1.5) satisfies ∥S(t)˜ u0 − u∗ (· + c1 ; s˜)∥X = O(e−β1 t ) for a positive number c1 ∈ [0, L), where s˜ ∈ Is is given as (1.39) and (˜ s, c1 ) → (s0 , 0) as δ1 → 0. 8
Remark 1.2 A similar result for the dynamical stability to Corollary 1.4 is given in Suzuki-Tasaki [17] by using the Lyapunov function. They prove the Lyapunov stability of the solution. On the other hand by virtue of Theorem 1.3 we can apply the result of [9] to proving the convergence to a solution with a small shift from s = s0 . Remark 1.3 The assertion of Theorem 1.3 implies that if the solution to (1.22) is stable with some nondegeneracy, then the corresponding solution to (1.20) is also stable. However, it does not assure that an unstable solution to (1.22) provides an unstable one to (1.20). However, if L0 in DN is invertible, then the zero eigenvalue of L is simple (see Lemma 3.1 in §3). Remark 1.4 The linearized stability (by (1.36)) of a nonconstant solution to (1.22) is studied in [17] for the Neumann boundary conditions. Due to their result any stable solution in a cylinder domain must be monotone in the axial direction (in the one-dimensional case, the solution is simply monotone). Moreover, if (1.3)-(1.4) with (1.26) allows a unique constant solution, and if it is unstable, we see that the system of equations has a stable invariant set containing a nonconstant solution with the monotonicity. Indeed, the omega limit set of any orbit consists of equilibrium solutions. (See also [17] for nonconstant solutions bifurcating from a constant one of the one-dimension.) With the aid of some numerics we will give bifurcation diagrams of equilibrium solutions for specific functions f (u, v) in (1.6), (1.7) and (1.8) as the parameter L or s varies. We can observe both supercritical and subcritical bifurcations for the 1-mode solution depending on the parameter values. In addition to the bifurcation of the equilibrium, in the case of (1.8), which is not for (1.26), the numerics tells that a Hopf bifurcation takes place from the 1-mode solution. We also give a mathematical result for a local bifurcation near the bifurcation point at which the constant solution looses the stability (see Theorem 4.1). Finally we give a remark on a characteristic dynamics of solutions to the model equations reported by [11] and [10]. For convenience of explanation, as in [10], we use the function (1.6) in (1.3) with small D and large L under the periodic boundary conditions. By their numerical simulations, adding a small perturbation to a constant steady state leads to the emergence of a spatial wave pattern (by the Turing instability). It grows until several peaks are created. After the growth of the peaks, some peaks go down. Then the distance of neighboring peaks becomes larger because of the middle ones disappear. If the distance is larger enough, it looks that the pattern settles down to a steady state. However, one peak (or more) collapses abruptly. Simultaneously, adjacent peaks grows more. This process continues until a single peak remains. This is quite an interesting long excursion from the Turing instability to a stable single peak. 9
We, however, discuss no mathematical treatment for this transient behavior in the present article. That would be a future study. We organize this paper as follows: In the next section we prove Lemmas 1.1 and 1.2. In §3 we prove Theorem 1.3 together with the Corollary 1.4. In §4 we show numerics for global bifurcation diagrams of equilibrium solutions for specific functions of (1.6), (1.7) and (1.8). Those numerical results are obtained by AUTO. In the same section we also give Theorem 4.1 for a normal form in the center manifold around the bifurcation point where the 1-mode solution emerges. As an Appendix we give a condition on f that allows the C k smoothness of the semiflow S(t) in X.
2
Proof of Lemmas
As mentioned in the introduction the Lyapunov function (1.25) is found in [13], [16] and [17] where the phase field model by Fix and Caginalp is treated. Combining this and a result in Hale[7] we can prove Lemma 1.1. Proof of Lemma 1.1: For a solution (u(x, t), v(x, t)) to (1.3) with (1.4) or (1.5) we compute d E(u(·, t), v(·, t)) dt ∫
τk [D∇u · ∇ut − (f0 (u) − kDu)ut + (Du + v)(Dut + vt )]dx 1 − τD Ω ∫ τk = − [(D∆u + f0 (u) − kDu)ut − (Du + v)(Dut + vt )]dx 1 − τD Ω ∫ τk (Du + v)(Dut + vt )]dx = − [|ut |2 − k(Du + v)ut − 1 − τD Ω ∫ ∫ k 2 = − |ut | dx + (Du + v)(ut + τ vt )dx Ω Ω 1 − τD ∫ ∫ k 2 = − |ut | dx − |∇(Du + v)|2 dx. Ω Ω 1 − τD =
Hence, for τ D < 1 we have d E(u(·, t), v(·, t)) ≤ 0, dt which implies the non-increasing of S(t)u0 in t. By the above computation, we also see d E(u(·, t), v(·, t)) = 0 (∀t ∈ R) dt yields ut = 0, ∇(Du + v) = 0 (∀t ∈ R). 10
Thus we see that E(u(·, t), v(·, t)) = constant(∀t ∈ R) implies (u(·, t), v(·, t)) is an equilibrium solution. Considering the assumption on F0 , we assert that E(u, v) is bounded below. We see that a bounded orbit S(t)u0 (t ≥ 0) is precompact by the smoothing effect of the reaction-diffusion equations. In the sequel we can apply Lemma 3.8.2 of [7] to obtain that the omega limit set for the solution is included in the set of all the equilibrium solutions. This concludes the proof. Lemma 1.2 is proved in [11] and [10]. We give a proof for the reader’s convenience. Proof of Lemma 1.2: Substituting ( ) ( ) ξ ϕ(x) = ϕj (x) ψ(x) η into (1.27) yields
( A0
ξ η
)
( = −λ
ξ η
) .
We obtain the equation for λ as det(A0 + λI) = λ2 − {(D + 1/τ )σj + fv /τ − fu }λ + {Dσj2 + (Dfv − fu )σj )}/τ = 0, where we simply write fu = fu (u, v),
fv = fv (u, v).
Under the condition (1.28), we have (D + 1/τ )σj + fv /τ − fu > 0,
Dfv − fu < 0.
This fact leads to the desired assertion.
3
Proof of Theorem 1.3
We first prove the linearized operator L has no imaginary eigenvalues. When f (u, v) = f0 (u) + kv, (1.31) turns to be ) ( ) ( ) ( ϕ ϕ D∆ϕ + f0′ (u∗ (· : s0 ))ϕ + kψ =λ . (3.1) L =− ψ {∆ψ − f0′ (u∗ (·; s0 ))ϕ − kψ}/τ ψ Put
a(x) := f0′ (u∗ (x; s0 )) − Dk, 11
and χ := Dϕ + ψ. Then (3.1) is transformed as D∆ϕ + a(x)ϕ + kχ + λϕ = 0, ∆χ + τ λχ + λ(1 − τ D)ϕ = 0.
(3.2) (3.3)
We let λ = λ1 + iλ2 ,
ϕ = ϕ1 + iϕ2 ,
χ = χ1 + iχ2 .
We will prove that λ2 must be zero. Multiply ϕ and χ with (3.2) and (3.3) respectively. Then integrating them, we have ∫ (−D|∇ϕ|2 + a(x)|ϕ|2 + kχϕ + λ|ϕ|2 )dx = 0, (3.4) ∫Ω (−|∇χ|2 + τ λ|χ|2 + λ(1 − τ D)ϕχ)dx = 0. (3.5) Ω
∫
Put α :=
∫ (ϕ1 χ1 + ϕ2 χ2 )dx,
(ϕ1 χ2 − ϕ2 χ1 )dx.
β :=
Ω
Ω
Take the imaginary part of (3.4) and write (3.5) the real part and imaginary part separately. Then kβ + λ2 ∥ϕ∥2 = 0, (3.6) { τ λ1 ∥χ∥2 + (1 − τ D)(λ1 α + λ2 β) = ∥∇χ∥2 , (3.7) τ λ2 ∥χ∥2 + (1 − τ D)(λ2 α − λ1 β) = 0, where ∥ · ∥ stands for the usual L2 norm. If β = 0, then λ2 = 0. Hence we may assume β ̸= 0. From (3.7) we see λ2 =
(1 − τ D)β∥∇χ∥2 . ((1 − τ D)α + τ ∥χ∥2 )2 + (1 − τ D)2 β 2
(3.8)
Eliminating β from (3.6) and (3.8) yields ( ) (1 − τ D)∥ϕ∥2 ∥∇χ∥2 1+ λ2 = 0. k{(1 − τ D)α + τ ∥χ∥2 )2 + (1 − τ D)2 β 2 } This implies λ2 = 0, namely, every eigenvalue must be real. Next, letting λ1 be the least eigenvalue of L in (1.31), we show that the assumption of λ1 < 0 leads to a contradiction. Define B := −∆ − τ λ1 From (3.3)
(D(B) = DN or Dp )
χ = λ1 (1 − τ D)B−1 ϕ 12
follows. Hence (3.2)-(3.3) are reduced to D∆ϕ + a(x)ϕ + λ1 (1 − τ D)kB[ϕ] = −λ1 ϕ Multiplying ϕ with this equation and integrate by part yields ∫ ∫ 2 2 {D|∇ϕ| − a(x)|ϕ| − λ1 (1 − τ D)kB[ϕ]ϕ}dx = λ1 |ϕ|2 dx. Ω
(3.9)
(3.10)
Ω
Corresponding to the linearized operator L0 defined by (1.36), consider (∫ )2 ∫ k(1 − τ D) 2 2 Q0 [φ] := {D|∇φ| − a(x)|φ| }dx + φdx . (3.11) τ |Ω| Ω Ω We note that the least eigenvalue µ1 of L0 is obtain by µ1 = inf{Q0 [φ] : φ ∈ D(L0 ), ∥φ∥ = 1}. By the Fourier expansion ϕ(x) =
∑
αj ϕj (x),
αj = ⟨ϕ, ϕj ⟩L2 ,
j≥1
we obtain B[ϕ] =
∑ j≥1
∑ αj αj α1 ϕj = − ϕ1 + ϕj , σj − τ λ1 τ λ1 σ − τ λ1 j≥2 j
thus ⟨B[ϕ], ϕ⟩L2 = −
∑ αj2 α12 + . τ λ1 j≥2 σj − τ λ1
(3.12)
Using (3.10), (3.11) and (3.12) together with ∫ 1 α1 = ϕdx, |Ω|1/2 Ω we obtain
∑
Q0 [ϕ] − k(1 − τ D)λ1
j≥2
αj2 = λ1 ∥ϕ∥2 σj − τ λ1
The assumption of Theorem 1.3 implies that µ1 ≥ 0. We thereby obtain Q0 [ϕ] ≥ 0. Consequently, by the assumption λ1 < 0, we get to −k(1 − τ D)
∑ j≥2
αj2 ≥ ∥ϕ∥2 σj − τ λ1
13
(3.13)
which yields ∥ϕ∥ = 0. This is a contradiction. Hence, λ must be nonnegative. Now we consider the case λ1 = 0. Then (3.2) and (3.3) imply D∆ϕ + a(x)ϕ + kχ = 0, Let χ(x) = c. Then D c= |Ω| Put
∆χ = 0.
∫
1 ϕdx + |Ω| Ω
1 ρ= |Ω|
(3.14)
∫ ψdx. Ω
∫ (ϕ + τ ψ)dx. Ω
From these two equations it follows that ρ 1 − τD c= − τ τ |Ω|
∫ ϕdx. Ω
Thus (3.14) are reduced to (1 − τ D)k L0 [ϕ] = D∆ϕ + a(x)ϕ − τ |Ω|
∫
ρ ϕdx = − . τ Ω
(3.15)
By the assumption of the theorem we can assert that only ϕ = ρw(·) is a solution of (3.15) for the Neumann case while all the solutions for the periodic case are given by ϕ = cu∗x (·; s0 ) + ρw(·). In order to complete the proof we need to show that there is no general eigenfunction of zero eigenvalue of L. We have the following lemma: Lemma 3.1 Under the same assumptions of Theorem 1.3 the generalized eigenspace of zero eigenvalue of L consists of c(w(x), cD −Dw(x))T (c is real number) for the Neumann case, while it consists of the linear combination of (w(x), cD −Dw(x))T and (u∗x (·; s0 ), vx∗ (·; s0 ))T for the periodic cases, where cD is given by (1.38). Proof.
We consider the periodic case. Let ) ( ) ( ) ( ∗ ϕ w(·) ux (·; s0 ) , L = c0 + c1 vx∗ (·; s0 ) ψ cD − Dw(·)
(3.16)
where c0 and c1 are real numbers. We note that by the assumption the null space of L in Dp × Dp consists of the linear hull of (w(·), cD − Dw(·))T and (u∗x (·; s0 ), vx∗ (·; s0 ))T . We prove that there is no solution for any (c0 , c1 ) ̸= (0, 0). This implies the assertion of the lemma. Define the adjoint operator ) ( ) ( p Dpxx + f0′ (u∗ (·; s0 ))p − (1/τ )f0′ (u∗ (·; s0 ))q ∗ . (3.17) L := − (1/τ )qxx + kp − (k/τ )q q 14
(
We solve L
∗
p q
)
( =
0 0
) ,
that is, Dpxx + f0′ (u∗ )(p − q/τ ) = 0,
(1/τ )qxx + k(p − q/τ ) = 0.
(3.18)
These equations are reduced to D(p − q/τ )xx + (f0′ (u∗ ) − Dk)(p − q/τ ) = 0. By the assumption of the simplicity of zero eigenvalue of L0 , we solve this equation as p − q/τ = αu∗x (α ∈ R). Applying this to the second equation of (3.18) yields qxx + τ kαu∗x = 0. After a simple integration, we thereby obtain ( ∫ ) ∫ p = αu∗x + αk − x u∗ dx + x L u∗ dx + β/τ, L 0 0 ) ( ∫ ∫ L x ∗ q = αkτ − u∗ dx + x u dx + β L 0 0 (α, β ∈ R). On the other hand the equation ) ) ( ( ϕ w L0 = c0 , ψ cD − Dw is solvable for nonzero c0 if and only if ∫ L {wp + (cD − Dw)q}dx = 0
(∀(p, q) ∈ N (L∗ )).
(3.19)
(3.20)
(3.21)
0
However, we see that for (p, q) = (1/τ, 1) ∫ 1 L (w/τ + cD − Dw)dx = 1/τ ̸= 0. L 0 This implies that there is no solution to (3.20) for c0 ̸= 0. Next consider ( ) ( ∗ ) ϕ ux , L = c1 vx∗ ψ which is solvable for c1 ̸= 0 if and only if ∫ L (u∗x p + vx∗ q)dx = 0 0
15
(∀(p, q) ∈ N (L∗ )).
(3.22)
(3.23)
Put
∫
x
Y (x) := − 0
x u dx + L ∗
∫
L
u∗ dx.
0
For (p, q) of (3.19) we easily verify ∫ L (u∗x p + vx∗ q)dx 0 ∫ L =α {(u∗x )2 + k(u∗x + τ vx∗ )Y }dx {0∫ L } ∫ L ∗ 2 ∗ ∗ =α (ux ) dx − k(u + τ v )Yx dx . 0
(3.24)
0
We compute as ∫
L
(u∗ + τ v ∗ )Yx dx 0 { } ∫ L ∫ L ∗ ∗ ∗ ∗ = (u + τ v ) −u + (1/L) u dx dx 0 0 ∫ L ∫ L ∗ 2 ∗ ∗ =− ((u ) + τ u v )dx + s u∗ dx, 0
and
∫
(3.25)
0
L
((u∗ )2 + τ u∗ v ∗ )dx 0 ( ))} ( ∫ ∫ L{ 1 − τD L ∗ 1 ∗ 2 ∗ ∗ s− u dx dx = (u ) + τ u −Du + τ L 0 0 )} ( ∫ ∫ L{ 1 − τD L ∗ ∗ 2 ∗ u dx dx = (1 − τ D)(u ) + u s − L 0 0 (∫ L )2 ∫ L ∫ L 1 − τD ∗ 2 ∗ ∗ = (1 − τ D) (u ) dx + s u dx − u dx . (3.26) L 0 0 0 Applying (3.26) to (3.25) yields ∫
L
0
{ ∫ (∫ L )2 } L 1 (u∗ + τ v ∗ )Yx dx = (1 − τ D) − (u∗ )2 dx + u∗ dx . L 0 0
Using
(∫
)2
L
∗
u dx
∫ ≤L
0
we obtain
L
(u∗ )2 dx,
0
∫
L
(u∗ + τ v ∗ )Yx dx ≤ 0
0
16
(3.27)
In the sequel, the left hand side of (3.24) must be positive. This implies that there is no solution of (3.22) for c1 ̸= 0. Since it is easier to prove the assertion for the Neumann case, we left the detail to the readers. By Lemma 3.1 we can complete the proof of Theorem 1.3. Proof of Corollary 1.4: First consider the Neumann boundary case. Let ( ) ∫ −1 G(u, s) := D∆u + f u, −Du + s/τ − (1/τ − D)|Ω| udx , (3.28) Ω
which is a mapping G(·, ·) : C 2+α (Ω) × (s0 − δ, s0 + δ) −→ C α (Ω), where 0 < α < 1 (n ≥ 2), α = 0 (n = 1). We note that G is of class C 1 . By the assumption, we have G(u∗ , s0 ) = 0 and the invertibility of (∂G/∂u)(u∗ , s0 ). We can apply the implicit function theorem to obtain a solution u = u∗ (·; s) satisfying G(u∗ (·; s), s) = 0, |s − s0 | < η, where η is a small number. Then u∗ (·; s) is C 1 with respect to s near s = s0 . If ∥˜ u0 − u∗ (·; s0 )∥ is small enough, then |s − s0 | < η. By the assumption on L0 , we see that L has no zero eigenvalue if we consider it in the space ∫ ∫ {(ϕ, ψ) ∈ D(L) : ϕdx + τ ψdx = 0}. Ω
Ω
˜ 0 with By virtue of the conservation of the mass, given u ∫ 1 s= (˜ u0 + τ v˜0 )dx, |s − s0 | < η, |Ω| Ω ˜ 0 ) − u∗ (·; s). Then we may consider v(·, t) in the space we put v(·, t) := u(·, t; u ∫ Xs := {u ∈ X : (u + τ v)dx = 0}. Ω
Hence, Theorem 5.1.1 of [9] is applicable. We obtain the assertion for the Neumann case. ˜ 0 , we can restrict Next consider the periodic case. As seen above, given u ∗ ˜ 0 ) − u (·; s) in the space Xs if s ∈ Is . This implies that the v(·, t) := u(·, t; u zero eigenvalue of L corresponding to (w(·), cD − Dw(·)) is ruled out when the flow is considered in Xs . Since the map c 7→ u(· + c; s) creates a C 2 curve in Xs , we can apply the result of [9] (Chap. 5, Exc. 6) to obtain the desired assertion. This concludes the proof.
17
4
Bifurcation diagrams
In this section we study the bifurcation structure of stationary solutions to (1.3) which bifurcate from the uniform state under the assumption (1.28), that is, the Turing unstable case. It should be noticed here that solutions to the Neumann problem can be considered as those of the periodic boundary problem with the period double size of the original domain. More precisely, suppose (u(x, t), v(x, t)) is a solution to (1.3) in Ω = [0, L] with the Neumann boundary conditions (1.4), and we extend the functions u and v for x ∈ [0, 2L] as follows: { u(x, t), x ∈ [0, L], u(x, t) = { u(2L − x, t), x ∈ (L, 2L], (4.1) v(x, t), x ∈ [0, L], v(x, t) = v(2L − x, t), x ∈ (L, 2L]. It is easy to check that these extended functions solve (1.3) in Ω = [0, 2L] with the periodic boundary conditions of period 2L. On the contrary, suppose that a pair of 2L-periodic functions (u(x, t), v(x, t)) is a solution to (1.3) and that they have even symmetry: u(−x, t) = u(x, t) and v(−x, t) = v(x, t). Then it solves (1.3)- (1.4). Therefore, we only consider the bifurcation structure of the stationary solutions with the periodic boundary conditions and even symmetry throughout this section. Let (u, v) be the constant equilibrium solution, that is, f (u, v) = 0 and u + τ v = s. We study a local bifurcation about this constant equilibrium. By noticing that u and v are even and 2L-periodic, the Fourier expansions of them about the equilibrium can be described as ∑ u(x, t) − u = α0 + 2 αm cos mk0 x, m>0 ∑ (4.2) v(x, t) − v = β0 + 2 βm cos mk0 x. m>0
Here, k0 := π/L, and both αm and βm are real valued functions of t for m ≥ 0. In addition to (4.2) we can also express the Fourier expansions as follows: ∑ u(x, t) − u = αm eimk0 x , m∈Z ∑ (4.3) v(x, t) − v = βm eimk0 x , m∈Z
where, α−m = αm and β−m = βm .
18
The equation (1.3) now takes the following decomposition in terms of the Fourier modes: α˙m = −m2 k02 Dαm + fˆm (m ≥ 0), (4.4) τ β˙m = −m2 k02 βm − fˆm (m ≥ 0). Here, fˆm denotes the m-th Fourier coefficient of the function f (u, v). Notice that α0 + τ β0 is independent of time t. We look for the solution which satisfies α0 +∫τ β0 = 0 since it is equivalent to find solutions to (1.3) with the constraint s=
(u+τ v)dx/|Ω|. In sequel the equations we have to solve are the following: Ω
α˙m = −m2 k02 Dαm + fˆm (m ≥ 0), τ β˙m = −m2 k02 βm − fˆm (m > 0), τ β0 = −α0 .
(4.5)
By virtue of the constraint α0 + τ β0 = 0 we can remove the degeneracy of the stationary problem and calculate the continuation of the stationary solution for a given constant s. To determine a local bifurcation equation about the constant equilibrium, we take the Taylor expansion of f at (u, v): f (u, v) = f (u, v) +fu · (u − u) + fv · (v − v) fvv fuu + (u − u)2 + fuv · (u − u)(v − v) + (v − v)2 2 2 fuuu + (u − u)3 + · · · + O(4), 6
(4.6)
where each coefficient is a corresponding partial derivative of f evaluated at (u, v). Then the zero’th order equation of (4.5) can be described as α˙0 = ( fu · α0 + f)v · β0 + O(2) fv = fu − α0 + O(2). τ
(4.7)
We know that the equilibrium α0 = 0 is stable by the assumption (1.28). The linearized matrix for the non-zero mode m of (4.5) is given by fv −m2 k02 D + fu Am = fu m2 k02 + fv . − − τ τ Since traceAm = −(1 + τ D)m2 k02 /τ + fu − fv /τ < 0 from (1.28), the m’th mode is unstable if and only if det Am < 0. Now the expression m4 k04 D − (fu − Dfv )m2 k02 det Am = τ 19
√ shows that det Am can be negative for k := mk0 ∈ (0, √(fu − Dfv )/D) and a bifurcation to a nontrivial solution takes place at k = (fu − Dfv )/D. Here we used the Turing unstable assumption fu − Dfv > 0 in (1.28). Let us take k0 as a bifurcation parameter and search bifurcated branches numerically. In addition to the parameter k0 bifurcation diagrams for s are also exhibited in the successive subsections. We also show the local bifurcation analysis corresponding to the numerical results in the last subsection.
4.1
Numerical results for (1.6)
We introduce numerically obtained bifurcation diagrams of stationary solutions to (1.3) for the specific cases f : (1.6). We solve (4.5) numerically by cutting off the Fourier modes up to 64 modes. All the simulations are calculated by AUTO, a software package for bifurcation analysis ([3]). We first consider k0 as a bifurcation parameter and set the other constants as follows: s= 5 D = 0.01 a= 1 b = 0.1 τ= 1 We show a branch of stationary solutions for the first mode by taking k0 as the horizontal axis in Figure 1. It turns out that a stationary solution bifurcates from the uniform state in subcritical pitchfork, therefore, both of the nontrivial stationary solution and the trivial state can be stable simultaneously for an appropriate k0 . It should be remarked that any branch of higher mode solutions can be obtained by a simple scaling of the first mode since k0 = π/L. We also show a bifurcation diagram for the bifurcation parameter s and for k0 = 2 in Figure 1. Here, the branch of the uniform steady state forms a S-shaped curve and both ends of the branch of nontrivial stationary solutions connect to the S-shaped curve in subcritical pitchfork.
4.2
Numerical results for (1.7)
In similar to to the previous subsection we introduce the numerically obtained bifurcation diagrams of the stationary solutions to (1.3) for the specific case of f : (1.7). We first consider k0 as a bifurcation parameter and set the other constants as follows: D = 0.01 a1 = 2.5 a2 = 0.7 τ= 1 20
6.50 6.25 6.00 5.75 5.50 5.25 SN 5.00
BP4 BP3
BP2
BP1
4.75 BP5 4.50 0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
k
12.5
0
2.00 1.75 SN
10.0
1.50 1.25
7.5
1.00 5.0
0.75
SN
0.50
BP
2.5 BP
s 2.5
BP
0.25
0.0 0.0
SN
5.0
7.5
10.0
12.5
0.00 1.00
s 1.25
1.50
1.75
2.00
2.25
2.50
Figure 1: Bifurcation diagrams of the stationary solutions for (1.6). Horizontal axis denotes the bifurcation parameter k0 or s. Vertical axis denotes L2 norm of the solution.
Bifurcation diagrams for different values of s are shown in Figure 2. It turns out that the direction of pitchfork bifurcations from the uniform state depends on s. As a result both of the nontrivial stationary solution and the trivial state can be stable simultaneously by taking an appropriate k0 and s.
4.3
Numerical results for (1.8)
In similar to the previous subsections, we show the numerically obtained bifurcation diagrams of the stationary solutions to (1.3) for the specific case of f : (1.8). We first consider k0 as a bifurcation parameter and the other constants are given as follows: s = 10 D = 0.1 τ = 1 α = 0.1 a1 = 0.5 a2 = 2.2 We show branches of stationary solutions for several modes by taking k0 as the horizontal axis in Figure 3. Next, we show, in Figure 4, another continuation from the previous 1-mode stationary solution at k0 = 1 by taking s ∈ [0.1, 10] as the bifurcation parameter 21
5.5
7.5
5.0
7.0
4.5
6.5
4.0
6.0
3.5
5.5 BP6
3.0 2.5 0.0 12.0
BP4 BP3
BP2
5.0
BP1
k 0.5
1.0
1.5
2.0
2.5
3.0
0
SN
BP5 BP4 BP3
SN
BP2
BP1
k
4.5 0.0 15.0
0.5
1.0
1.5
2.0
2.5
0
SN
12.5
11.5
3.0
10.0
SN
11.0 7.5 10.5
BP1
SN 5.0
BP2 BP6
10.0
9.5 0.0
2.5
BP3
SN BP5
BP 0.5
BP4
k 1.0
1.5
2.0
2.5
3.0
0
s
0.0 0.0
2.5
5.0
7.5
10.0
12.5
15.0
Figure 2: Bifurcation diagrams of the stationary solutions for (1.7) for s = 3.0, 5.0, 10.0, respectively. The horizontal axis denotes the bifurcation parameter k0 and the vertical axis denotes L2 norm of the solution. We also show a bifurcation diagram for s in the bottom right.
35.
30.
25. BP6 20.
BP5 BP4 BP3
15.
BP2 BP1
10. 0.
k 1.
2.
3.
4.
5.
0
Figure 3: Bifurcation diagram of the stationary solutions for (1.8). The horizontal axis denotes the bifurcation parameter k0 and the vertical axis denotes L2 norm of the solution. The bifurcating branches for m = 1, 2, 3, 4, 5 modes from the right are observed.
instead of k0 . Then we find a Hopf bifurcation on this branch. In fact we can numerically observe periodic motions by solving (1.3) with small initial data.
22
17.5
35
15.0
30
28
’u29.txt’ ’u28.txt’ ’u27.txt’ ’u26.txt’ ’u25.txt’ ’u24.txt’ ’u23.txt’ ’u22.txt’ ’u21.txt’
29 25
12.5 27
20
10.0
15
26
7.5
10
5.0
24 21
2.5
22
5
25
0
23
s
0.0 0.0
2.5
5.0
7.5
10.0
12.5
-5 0
20
40
60
80
100
Figure 4: The left figure shows a branch of the stationary solutions for (1.8). The horizontal axis denotes the bifurcation parameter s and the vertical axis denotes L2 norm of the solution. The black square indicates a Hopf bifurcation point. The scaled profiles of the corresponding periodic solutions u on the branch are shown in the right figure.
4.4
Bifurcation analysis
As mentioned earlier, the bifurcation√ to the nontrivial solution from the uniform state takes place when k = mk0 = (fu − Dfv )/D. In this subsection let us focus on the bifurcation to the first branch, namely, m = 1. The purpose of this sub-section is to derive√the normal form of the local bifurcation to the first branch at k0 = k ∗ := (fu − Dfv )/D. The matrix A1 has 0 eigenvalue if and only if k0 = k ∗ . In order to calculate the critical mode explicitly we take a new coordinate (η, θ) instead of (α1 , β1 ) so that η corresponds to the critical mode. More precisely, when k0 = k ∗ the eigenvalues of A1 are 0 and fu traceA1 = Dfv − Dτ . The linearized matrix A1 for 1-mode can be diagonalized ( ) 0 0 −1 T A1 T = 0 traceA1 by
( T =
1 −fv fu −D Dτ
) .
The new coordinate (η, θ) for 1-mode is defined by ) ( ) ( ( ) η − fv θ η α1 =T = . fu β1 θ −Dη + Dτ θ We are interested in a small δ-neighborhood of 0 for the critical mode η. All the other modes are slave mode in the sense that those are of higher orders in the center manifold reduction. In sequel we can conclude α1 = η + O(δ 2 ), β1 = −Dη + O(δ 2 ). 23
(4.8)
By plugging the Fourier expansion into (4.6) we obtain the equation for 1-mode as follows: ( ) ( ) ( ) α˙1 α1 1 + O(δ 4 ), = A1 +N − τ1 β1 β˙1 where the quadratic and cubic nonlinear terms N of fˆ1 are given by ∑ fuu ∑ fvv ∑ αm1 αm2 + fuv αm1 βm2 + βm βm 2 m +m =1 2 m +m =1 1 2 m +m =1 1 2 1 2 1 2 ∑ ∑ fuuu fuuv + αm αm αm + αm αm βm 6 m +m +m =1 1 2 3 2 m +m +m =1 1 2 3 1 2 3 1 2 3 ∑ ∑ fuvv fvvv + αm βm βm + βm βm βm 2 m +m +m =1 1 2 3 6 m +m +m =1 1 2 3
N =
1
2
3
1
2
3
If m1 ̸= ±1 and m2 ̸= ±1 in the above quadratic terms, then the corresponding quadratic term such as αm1 αm2 is O(δ 4 ) since both αm1 and αm2 are slave modes. Therefore, from the above quadratic sums, O(δ 3 ) terms appear only when one of mj = ±1 (j = 1, 2). It turns out that the possible pairs (m1 , m2 ) are (1, 0), (0, 1) and (−1, 2), (2, −1). Thus we need to take into account the cooperate modes, that is, m = 0, 2. Moreover, if one of mj (j = 1, 2, 3) is a slave mode then the corresponding cubic term such as αm1 αm2 αm3 is O(δ 4 ). Therefore, from the above cubic sums, the O(δ 3 ) terms appear only when mj = ±1 for j = 1, 2, 3. Thus the possible triples (m1 , m2 , m3 ) are (1, 1, −1), (1, −1, 1), (−1, 1, 1). We then calculate the O(δ 3 ) terms of N as follows: N = fuu (α1 α0 + α−1 α2 ) + fuv (α1 β0 + α−1 β2 + α0 β1 + α2 β−1 ) + fvv (β1 β0 + β−1 β2 ) 3fuvv fvvv 3 fuuu 3 3fuuv 2 α1 + α1 β1 + α1 β12 + β + O(δ 4 ). + 2 2 2 2 1 Here we used the fact that α−1 = α1 and β−1 = β1 . By plugging (4.8) and the constraint β0 = −α0 /τ into the above we further obtain N = {fuu − Dfuv − (fuv − Dfvv )/τ }ηα0 +(fuu − Dfuv )ηα2 + (fuv − Dfvv )ηβ2 1 + (fuuu − 3fuuv D + 3fuvv D2 − fvvv D3 )η 3 + O(δ 4 ). 2 Now we write the equation for η as follows: η˙ = µη + T0 ηα0 + T1 ηα2 + T2 ηβ2 + T3 η 3 + O(δ 4 ),
(4.9)
where µ is the critical eigenvalue which satisfies µ = 0 at the critical point and
24
the constants are: fu − Dfv {fuu − Dfuv − (fuv − Dfvv )/τ }, fu − D2 τ fv fu − Dfv T1 = (fuu − Dfuv ), fu − D2 τ fv fu − Dfv T2 = (fuv − Dfvv ), fu − D2 τ fv fu − Dfv T3 = (fuuu − 3fuuv D + 3fuvv D2 − fvvv D3 ). 2(fu − D2 τ fv ) T0 =
(4.10)
Next, we obtain the equation for the cooperate modes similarly as follows. fv )α0 + 2N2 η 2 + O(δ 3 ), τ ) ( ) ( ) ( α˙2 α2 1 2 + O(δ 3 ), = A2 + N2 η β2 − τ1 β˙2 α˙0 = (fu −
(4.11) (4.12)
where N2 = (fuu − 2Dfuv + D2 fvv )/2. Now we find a suitable near-identity transformation for the critical mode η ξ = η − S0 ηα0 − S1 ηα2 − S2 ηβ2 ,
(4.13)
so that the equation for the critical mode can be written only by the critical mode. It is the so-called normal form transformation by which we obtain the reduced equation on the center manifold. Our task here is to determine the unknown constants S0 , S1 , S2 in order that the terms including slave modes vanish. By using the inverse transformation of the local diffeomorphism (4.13) η = ξ + S0 ξα0 + S1 ξα2 + S2 ξβ2 + O(ξ 4 ), we obtain the transformed equation for ξ: ξ˙ = µξ −2µS0 ξα0 − 2µS1 ξα2 − 2µS2 ξβ2 +T0 ξα0 + T1 ξα2 + T2 ξβ2 (
) α2 −S0 (fu − − ξ(S1 S2 )A2 β2 −N2 (2S0 + S1 − Sτ2 )ξ 3 + T3 ξ 3 + O(δ 4 ).
(4.14)
fv )ξα0 τ
Since it is sufficient to consider the critical case µ = 0, the unknown constants are obtained by solving the following equations: (fu − fτv )S0 = T0 , (S1 S2 )A2 = (T1 T2 ) 25
A simple calculation shows that det A2 =
12 (f τD u
− Dfv )2 and
( )−1 S2 fv DT2 − T1 2S0 + S1 − = 2 fu − T0 + . τ τ 3(fu − Dfv ) In consequence we have the cubic coefficient c1 of (4.14) explicitly as } { ( )−1 DT2 − T1 fv . T0 + c1 = T3 − N2 2 fu − τ 3(fu − Dfv )
(4.15)
The next theorem tells the local bifurcation analytically. Theorem 4.1 Suppose c1 is not zero at the constant equilibrium (u, v) then the dynamics on the center manifold of the equation (4.5) near the critical point k0 = k ∗ is given by ξ˙ = µξ + c1 ξ 3 + O(δ 4 ). By virtue of this theorem we can assert that the nontrivial stationary solution emerges in supercritical pitchfork bifurcation when c1 is negative. Remark 4.1 Not only k0 but also other constant such as s can be considered ∂µ ̸= 0. Moreover, we can see that the same type as a bifurcation parameter if ∂σ of bifurcation occurs as such a parameter varies. Remark 4.2 As D tends to zero, the limit of c1 gives a simpler criterion for the pitchfork bifurcation : lim c1 =
D→0
2 fuuu (τ fuu − fuv )fuu fuu − + 2 τ fu − fv 6fu
(4.16)
Let us determine the cubic coefficient c1 of the normal form in the case of (1.8): f (u, v) = a1 (u + v)[(αu + v)(u + v) − a2 ]. Suppose τ = 1 then u+v = s. Since f (u, v) = 0 means (αu+v)(u+v) = a2 , i.e., (αu + v) = a2 /s, we can calculate the partial derivatives of f at (u, v) explicitly. Then we have lim c1 = 3a1 α −
D→0
2a1 (4αs2 + 5a2 )(2αs2 + a2 ) . 3s2 (αs2 + a2 )
We calculate this value with the particular set of parameters given in the previous sub-section so that lim c1 = 3·0.5·0.1−
D→0
2 · 0.5(4 · 0.1 · 102 + 5 · 2.2)(2 · 0.1 · 102 + 2.2) = −0.1836 · · · , 3 · 102 (0.1 · 102 + 2.2) 26
therefore, we see that the pitchfork bifurcation observed in the case of (1.8) is supercritical for sufficiently small D. We can similarly evaluate c1 in the case of (1.7) as follows. lim c1 =
D→0
a1 a22 (11a22 s2 − 44a2 s + 35) 3(a2 s − 1)(a2 s + 1)5
It turns out that for s > 1/a2 ≈ 1.429 there are two critical values: s = s1 ≈ 1.565 and s = s2 ≈ 4.1493 where the sign of c1 changes. It is consistent to the results shown in Figure 2.
5
Appendix
Here we give a condition on f under which the semiflow S(t) defined by (1.10) is C k in X = H 1 (Ω) × H 1 (Ω). First let A in L2 (Ω) × L2 (Ω) be a closed operator defined by ( ) D 0 Au := − ∆u + u, 0 1 with the domain D(A) := {u ∈ H 2 (Ω) × H 2 (Ω) : ∂u/∂ν = 0 on ∂Ω} or D(A) := {u ∈ H 2 (0, L) × H 2 (0, L) : u(0) = u(L), ux (0) = ux (L)}. Then A is sectorial and it allows a fractional power Aα (0 ≤ α < 1). It is known that D(A1/2 ) = H 1 (Ω) × H 1 (Ω), or D(A1/2 ) = {u ∈ H 1 (0, L) × H 1 (0, L) : u(0) = u(L)}. By the existence theorem found in [9] we see that a solution map (u0 (·), v0 (·)) 7→ (u(·, t; u0 , v0 ), v(·, t; u0 , v0 )) is C k in H 1 (Ω) × H 1 (Ω) if fˆ : (u(·), v(·)) ∈ H 1 (Ω) × H 1 (Ω) 7→ f (u(·), v(·)) ∈ L2 (Ω) × L2 (Ω)
(5.1)
is C k . For the one-dimensional case, since C 0 [0, L] is imbedded in H 1 ([0, L]), fˆ of (5.1) is C k provided that f (u, v) is a C k function. As for the higher-dimensional case, we have the following lemma (given in Hale-Vegas[8]): 27
Lemma 5.1 Let Ω is a bounded domain of Rn (n ≥ 2) with smooth boundaries. Suppose that f (u, v) is C 1 and it has zero. If the partial derivatives fu and fv are uniformly bounded in R2 , then fˆ of (5.1) is C 1 . Proof.
For simplicity of notation we consider the case that f = f (u). Suppose sup |f ′ (u)| < ∞. u∈R
Then the map fˆ fˆ : H 1 (Ω) −→ L2 (Ω) makes sense since f (u) has zero. We prove that the Gataux derivative DG fˆ(u) (u = u(·) ∈ H 1 (Ω)) DG fˆ(u) : w(·) ∈ H 1 (Ω) 7→ f ′ (u(·))w(·) ∈ L2 (Ω) is continuous with respect to u ∈ H 1 (Ω). Then it turns that fˆ is Fr´eche continuously differentiable and the derivative coincides with DG fˆ. To prove the continuity of DG fˆ(u), let {uj } be any sequence in H 1 (Ω) which converges to u and estimate ∥DG fˆ(uj ) − DG fˆ(u)∥L(H 1 (Ω);L2 (Ω)) . (5.2) There exists a subsequence {uji } such that uji converges to u a.e. By the 2n Schwarz inequality and the imbedding H 1 (Ω) ⊂ L n−2 (Ω) (if n = 2, H 1 (Ω) ⊂ Lq (Ω) for q ≥ 2), we obtain ∥{DG fˆ(uji ) − DG fˆ(u)}(w)∥L2 ≤ ∥f ′ (uji (·)) − f ′ (u(·))∥Lp ∥w∥Lq ≤ C∥f ′ (uji (·)) − f ′ (u(·))∥Lp ∥w∥H 1 , where we take p = n, q = 2n/(n − 2) for n ≥ 3. Since |f ′ | is bounded, the Lebesgue convergence theorem tells ∫ ∫ ′ ′ p lim |f (uji (x)) − f (u(x))| dx = lim |f ′ (uji (x)) − f ′ (u(x))|p dx = 0. ji →∞
Ω
Ω ji →∞
Since {uj } is an arbitrary sequence, we obtain the continuity of DG fˆ(u). This concludes the proof.
Acknowledgements The authors would like to express their sincere acknowledgement to Professor Shuji Ishihara (University of Tokyo) for letting them know his works in a nice lecture. Their thanks also go to Professor Takashi Suzuki(Osaka University) and Doctor Souhei Tasaki (Osaka University) for the paper on a study of the phase field model. It turned to be very helpful for the present work. The research were supported in part by the Grant-in-Aid for Scientific Research (B) No.19340026 and (C) No.20540116, Japan Society for the Promotion of Science. 28
References [1] R. E. Baker, E. A. Gaffney and P. K. Maini, Partial differential equations for self-organization in cellular and developmental biology, Nonlinearity 21 (2008), R251-R290. [2] G. Caginalp, An analysis of a phase field model of a free boundary, Arch. Rational Mech. Anal. 92 (1986), 205-245. [3] Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A., Sandstede, B., Wang, X., AUTO 97: Continuation and bifurcation software for ordinary differential equations (with HomCont) (1997). [4] P. C. Fife, Mathematical Aspects of Reacting and Diffusing Systems, Springer-Verlag, Berlin-New York, 1979. [5] G. J. Fix, Phase filed methods for free boundary problems, Free Boundary Problems: Theory and Applications, Eds. A. Fasano, M. Primicero, Pitman, London, 1983, 580-589. [6] A. Gierer and H. Meinhardt, Biological pattern formation involving lateral inhibition, Some mathematical questions in biology, VI, Lectures on Math. in the Life Sciences 7, Amer. Math. Soc. Providence, R. I. (1974) 163-183. [7] J. K. Hale, Asymptotic behavior of dissipative systems, Math. Surveys and Monographs, 25, American Mathematical Society, Providence, R. I. 1988. [8] J. K. Hale and J. M. Vegas, A nonlinear parabolic equation with varying domain, Arch. Rational Mech. Anal. 86 (1984), 99-123. [9] D. Henry, Geometric Theory of Semilinear Parabolic Equations, SpringerVerlag, Berlin-New York, 1981. [10] S. Ishihara, M. Otsuji, and A. Mochizuki, Transient and steady state of mass-conserved reaction-diffusion systems, Physical Review E 75 (2007), 015203(R). [11] M. Otsuji, S. Ishihara, C. Co, K. Kaibuchi, A. Mochizuki, and S. Kuroda, A mass conserved reaction-diffusion system captures properties of cell polarity, PLoS Computational Biology 3 (2007), 1040-1054. [12] M. Mimura, M. Tabata and Y. Hosono, Multiple solutions of two-point boundary value problems of Neumann type with a small parameter, SIAM J. Math. Anal. 11 (1980), 613-631. [13] T. Miyasita, Global existence and exponential attractor of solutions of FixCaginalp equation, preprint. 29
[14] J. D. Murray, Mathematical Biology, Springer-Verlag, 1989. [15] Y. Nishiura and H. Fujii, Stability of singularly perturbed solutions to systems of reaction-diffusion equations, SIAM J. Math. Anal. 18 (1987), 1726-1770. [16] T. Suzuki, Mean Field Theories and Dual Variation, Atlantis Press, World Scientific, 2008. [17] T. Suzuki, and S. Tasaki, Stationary Fix-Caginalp equation with non-local term, to appear in Nonlinear Anal. [18] A. M. Turing, The chemical basis of morphogenesis, Phil. Trans. Roy. Soc. Lond. B237 (1952), 37-72.
30
