Entropy/Energy stable schemes for evolutionary dispersal models
This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/authorsrights
Author's personal copy Journal of Computational Physics 256 (2014) 656–677
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Entropy/energy stable schemes for evolutionary dispersal models Hailiang Liu, Hui Yu Iowa State University, Mathematics Department, Ames, IA 50011, United States
a r t i c l e
i n f o
Article history: Received 4 March 2013 Received in revised form 13 August 2013 Accepted 17 August 2013 Available online 17 September 2013 Keywords: Dispersal Diffusion–advection–reaction equation Entropy Energy Positivity preserving
a b s t r a c t In this paper we propose some entropy/energy stable finite difference schemes for the reaction–diffusion–advection equation arising in the evolution of biased dispersal of population dynamics. The peculiar feature of these active dispersal models is that the transient solution converges to the stable steady state when time goes to infinity. For the numerical method to capture the long-time pattern of persistence or extinction, we use the relative entropy when the resource potential is logarithmic, and explore the usual energy for other resource potentials. The present schemes are shown to satisfy three important properties of the continuous model for the population density: (i) positivity preserving, (ii) equilibrium preserving, and (iii) entropy or energy satisfying. These ensure that our schemes provide a satisfying long-time behavior, thus revealing the desired dispersal pattern. Moreover, we present several numerical results which confirm the second-order accuracy for various resource potentials and underline the efficiency to preserve the large time asymptotic. © 2013 Elsevier Inc. All rights reserved.
1. Introduction This work is concerned with the numerical approximation of a class of reaction–diffusion–advection equations arising in the evolution of biased dispersal of population dynamics, with emphasis on exploring the entropy/energy structure of the dispersal model so that the resulting methods provide a satisfying long-time behavior. 1.1. Mathematical formulations Reaction–diffusion equations have been widely used to model the biological problems [35,36]. One of the well-known examples is the logistic reaction–diffusion model for the population growth with random dispersal,
∂t u = u + λu (m − u ) in Ω × (0, ∞),
(1.1)
where the population inhabiting a bounded domain Ω ⊂ Rd has density u (x, t ) at location x and time t, and local growth rate m. The parameter λ > 0 is the reverse of the dispersal rate. If the environment is heterogeneous, i.e., m(x) is not a constant, then the population may have a tendency to move toward resources in addition to the random movement. The model may be upgraded to the following form
∂t u = ∇ · (∇ u + u ∇ P ) + λu (m − u ) in Ω × (0, ∞), E-mail addresses: [email protected] (H. Liu), [email protected] (H. Yu). 0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.08.032
(1.2)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
657
where P = P (m), which we call resource potential, reflects the movement tendency of the population. The time evolution is subjected to both the initial density u (x, 0) = u 0 (x) and the zero flux boundary condition
(∇ u + u ∇ P ) · ν = 0 on ∂Ω × (0, ∞),
(1.3)
where ν is the outward normal vector on the boundary ∂Ω which is assumed to be smooth. Several dispersal strategies have been studied in literature, for example P = −αm in [2,5,9,13,15], and the obtained results may apply to a more general reaction than
F (u , m) = u (m − u ), as long as it satisfies F (0, m) = 0 and F (+∞, m) < 0. With the logistical reaction term or its variation, the biologically relevant solutions to (1.2) are nonnegative and ultimately bounded. It follows from the regularity theory of parabolic partial differential equations that in the state space bounded orbits are pre-compact, hence the semi-flow will have a compact attractor, and the ω -limit set of any initial state u 0 (x) will be a compact invariant set. In other words, as time evolves the solution of (1.2) is expected to approach some stable patterns, independent of the choice of initial density; see, e.g., [7]. The main result reviewed in Section 2 may be stated as follows.
¯ is positive somewhere in Ω and P is smooth in m. There exists a unique λ∗ > 0 and a positive Theorem 1.1. Suppose that m ∈ C 2 (Ω) equilibrium solution to (1.2) if and only if
me − P dx < 0.
Ω
Moreover, λ∗ = 0 if and only if Ω me − P dx 0. (1) If 0 < λ λ∗ , all nonnegative solutions of (1.2) decay toward zero as t → ∞. (2) If λ > λ∗ , the positive equilibrium is globally attractive among nonzero nonnegative solutions. Of special interest is the ideal free dispersal strategy determined by P such that at an equilibrium u eq (x) [8,10,11], both
∇ · (∇ u eq + u eq ∇ P ) = 0 and m − u eq = 0 in Ω. This will hold if P = − ln m (unique up to some constant). In such a setting, the species u can perfectly match the environmental resource. The role of ideal free distributions has been well recognized in literature. Moreover, the quantity
E [u ] =
u log
u m
+ m − u dx
(1.4)
Ω
is nonincreasing in time. This corresponds to some physical entropy relative to the equilibrium state u eq (x) = m(x), called relative entropy. For general resource potential P , the quantity
V [u ] =
1 − P P 2 e ∇ ue − λe P 2
u
F (ξ, m) dξ dx
(1.5)
0
Ω
is nonincreasing in time. Moreover, Eq. (1.2) can be rewritten as a gradient flow,
∂t u = −e − P
δV , δu
(1.6)
where δδVu denotes the standard variational derivative of V with respect to u, and the functional V often corresponds to some free energy of the underlying physics, and thus is called the energy. In the context of two competing species, the evolutionary dynamics of conditional dispersal becomes much more complex, see [13,22]. A typical model may be described as
⎧ ut = α ∇ · (∇ u + u ∇ P ) + u (m − u − v ) ⎪ ⎪ ⎨ v = β∇ · (∇ v + v ∇ Q ) + v (m − u − v ) t ⎪ u + u ∂ ν P = ∂ν v + v ∂ ν Q = 0 ∂ ⎪ ⎩ ν u (x, 0) = u 0 (x), v (x, 0) = v 0 (x)
in Ω, in Ω, on ∂Ω,
(1.7)
in Ω.
The two nonnegative constants α and β represent the dispersal rates of two species, respectively. In this paper we shall focus on the case when α = β , and we can envision that system (1.7) models two competing species that are identical except their different dispersal strategies P and Q . An interesting question is whether there is any strategy in system (1.7)
Author's personal copy 658
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
that is globally evolutionarily stable, in the sense that the semi-trivial steady state (m, 0) is always locally stable. This is the case if the dispersal strategy is ideal free at (m, 0), i.e., P (m) = − ln m + C and Q (m) = P (m); see [11,24]. Moreover, the following functional
E [u , v ] =
u dx + Ω
v dx − Ω
m(x) ln u dx
(1.8)
Ω
is nonincreasing in time. For other results on dispersal strategies in heterogeneous landscapes, we refer to [13,27,28] for the study of aggregation profiles in terms of the resource distribution, and [16,26] for nonlocal dispersal strategies. The aim of this paper is to give reliable numerical schemes for (1.2) and (1.7) from the perspective of providing a satisfying long-time behavior. A key fact is that they both admit certain entropy/energy structure, and we demand our numerical schemes to satisfy the entropy/energy decreasing property in the discrete setting. In addition, both positivity and steady states for (1.1) are required to be preserved as well. These three requirements together are important for system (1.1), yet adding levels of difficulty to the design of a numerical method with high accuracy. As a preliminary attempt, only some simple time stepping is discussed in the present paper. In recent years, energy preservation or dissipation numerical schemes have drawn much interest and been extensively studied for several PDEs, since they are more likely to give a better performance in the long-time simulation; see e.g. [1, 12,23,30,33,34]. To the best of our knowledge, however, no entropy stable scheme for (1.2) or (1.7) has yet appeared in the literature. Though, due to the physical and mathematical significance of advection–diffusion–reaction equations, numerical methods have been developed and analyzed through decades; see, e.g., [3,4,29,31] and the recent monograph [25]. If the model equation has a nice property such as energy conservation/dissipation, then design of an energy satisfying numerical method is often desired. The Crank–Nicolson time discretization has been proven successful in many cases. For example, in [14], eight finite difference schemes are compared for solving the nonlinear Schrödinger equation, and two of three conservative schemes are of the Crank–Nicolson type. Efforts have also been devoted to the development of procedures to mimic the continuous energy property, see e.g., [1,19,23,32]. For example, the average vector field method is developed in [1] for some Hamiltonian PDEs and dissipative PDEs. In an earlier work by Furihata [23], a procedure was given to derive energy satisfying (conservation or dissipation) finite difference schemes for the class of ∂t u = ∂xα ( δδGu ), about which the
d G = 0 or 0. For higher order nonlinear dissipative discrete calculus is used to mimic the continuous energy property dt PDEs such as the Cahn–Hilliard equation, efficient time stepping techniques which can inherit the unconditional energy stability are particularly in demand. Several methods such as the convex splitting by Eyre [20] and the stabilized time discretization which allow large time steps have been explored in literature, see e.g., [17,18] and references therein.
1.2. Main steps of this paper In Section 2 mathematical theory and properties regarding the target equation (1.2) are presented; we particularly highlight the energy/entropy structure, as well as the large time pattern formation in terms of the parameter λ. The key idea is to reformulate (1.2) and (1.7) into the nonlogarithmic Landau form in terms of g = ue P and h = ve Q , following the strategy in [30] by the authors for the Fokker–Planck equation of FENE dumbbell polymers. Section 3 is devoted to the case with the logarithmic potential and the proposed scheme is presented. The schemes at both semidiscrete and fully discrete levels are shown to satisfy three important properties of the continuous model for the population density: (i) positivity preserving, (ii) equilibrium preserving, and (iii) entropy satisfying. The schemes presented in Section 4 apply to the case with general resource potential and is extensible to handle more general reaction terms. In that section we first present a semidiscrete scheme with discretization only in time. This scheme is proved to be uniquely solvable provided the time step satisfies a bound depending on the growth rate m and the numerical solution at the previous time step. The scheme preserves the equilibrium solution with numerical energy being nonincreasing in time for arbitrary time step. In addition, the scheme can be expressed as a second-order prediction–correction method, and the prediction is shown to be in the same range of the numerical solution at the previous time step. As a consequence, we prove that the fully discrete scheme satisfies both desired properties: energy satisfying and equilibrium preserving. We also propose an entropy stable scheme in Section 5 to simulate the ideal free distribution for two competing species. Section 6 is devoted to numerical tests of proposed schemes. Finally some concluding remarks are presented in Section 7. 2. Model reformulation and mathematical properties In order to reveal the dynamic picture of (1.2), following [30], we rewrite (1.2) in terms of g = u / M with M = e − P :
M ∂t g + P g = 0 in Ω,
(2.1)
where the operator P is defined as
P g = −∇ · ( M ∇ g ) − λ M g (m − M g ).
(2.2)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
The initial condition is g 0 (x) = M g (m − M g ).
u 0 (x) M (x)
659
and on the boundary ∂Ω we have M g ν = 0. The reaction term becomes F ( g , m) =
Remark 2.1. M is allowed to be zero on boundary, which corresponds to the case that P is singular in the original equation. This is an important case in biological dispersal with an inapproachable boundary. Our reformulation has advantage to convert this singularity to certain degeneracy in term M ∂t g, which can be easily treated with some implicit time stepping, as done in [30] for the FENE polymer model. 2.1. Entropy for log potential P (m) = − ln m In such case M = e − P = m and (2.1) reduces to
m∂t g = ∇ · (m∇ g ) + λm2 g (1 − g )
in Ω.
(2.3)
The logistic reaction term ensures that any solution remains nonnegative and bounded, i.e.,
0 g (x, t ) max 1, g 0 ∞
∀(x, t ) ∈ Ω × [0, ∞).
As t → ∞, g will converge to the stable equilibrium solution g˜ ≡ 1, independent of the choice of the initial data. This motivates us to explore the entropy which might have some biological significance. The entropy defined in (1.4) becomes
E [g] =
m g log g + (1 − g ) dx. Ω
As a convex functional of g, it defines the entropy structure of (2.3):
m∂t g = ∇ · mg ∇ m−1
δE δg
+ λ F ( g , m),
where the variational derivative δδ Eg = m log g. The reaction term F ( g , m) = m2 g (1 − g ) satisfies
F·
δE 0 δg
for nonnegative g. Using integration by parts, we have
d dt
E =−
2 δE δ E mg ∇ m−1 dx + λ m −1 F · dx 0, δg δg
Ω
Ω
where the zero flux boundary condition has been taken into account. We shall propose a numerical method to satisfy this entropy property. 2.2. Energy for all other potential In general, the evolution may lead to extinction or persistence of the species, and the concept of entropy is no longer adequate to control the reaction term. We shall follow the standard gradient flow idea [7] by using the functional defined in (1.5),
V [g] =
1 2
M |∇ g | − G x, g (x, t ) dx, 2
(2.4)
Ω
where G is chosen as
g G (x, g ) = λ
F (s, m) ds =
λ 2
m(x) M (x) g 2 −
λ 3
M 2 (x) g 3 ,
0
so that (2.1) can be rewritten as
M ∂t g = ∇ · ( M ∇ g ) + ∂ g G . The above energy also defines the gradient flow structure of (2.5):
(2.5)
Author's personal copy 660
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
M ∂t g = −
δV , δg
where the zero flux has been used in the variational derivative. A formal calculation gives
d dt
δV ∂t g dx = − δg
V = Ω
M (∂t g )2 dx 0.
(2.6)
Ω
We shall propose a numerical method to satisfy this energy dissipation inequality in the discrete setting. 2.3. Large time pattern formation The nonnegative equilibrium solution g˜ (x), if it exists, solves
P g˜ = 0 in Ω
M g˜ ν |∂Ω = 0.
and
Let g = g˜ + ε w in (2.1) and then send
(2.7)
ε → 0 to obtain the linearized problem
M ∂t w = ∇ · ( M ∇ w ) + λ M w (m − 2M g˜ )
in Ω
M w ν |∂Ω = 0.
and
(2.8)
It is known that if g˜ is linearly stable, then g˜ is stable; for a proof see, e.g., [38, Theorem 11.22]. The linear stability of g˜ is determined by the spectrum of the linearized operator of P . For g˜ = 0, it boils down to the eigenvalue problem:
∇ · ( M ∇φ) + λmM φ = σ M φ in Ω and M φν |∂Ω = 0. (2.9) When σ < 0, 0 is stable; when σ > 0, 0 is unstable. Hence σ = 0 is a threshold value, corresponding to a special parameter λ∗ which satisfies
∇ · ( M ∇ψ) + λ∗mM ψ = 0 in Ω and M ψν |∂Ω = 0.
(2.10)
In what follows we shall denote the weighted L 2 norm by v M with
v 2M =
M v 2 dx. Ω
The basic results for (2.1), (2.7)–(2.10), leading to Theorem 1.1, are presented in Appendix A. These results are known for some specific choices of P such as P = −αm [2,7], yet a complete account for general resource potentials seems unavailable in literature. For the sake of completeness, we present detailed results for (2.1) together with the key arguments in Appendix A. 3. Entropy stable method for one-species model 3.1. Semidiscrete scheme in space We present only the one-dimensional case. Generalizations to the multidimensional space are straightforward for tensor product grids and the results remain valid without modifications. Given an integer N, we partition the domain Ω := [a, b] by using a uniform mesh x, where
b−a
x =
N
and
xj = a +
j−
1 2
x.
Note that x 1 = a and x N + 1 = b. Let g j (t ) be the approximation of g (x j , t ). We propose the following semidiscrete scheme: 2
2
for 1 j N,
mj
d dt
g j = D − (m j + 1 D + g j ) + λm2j g j (1 − g j ) 2
and
g j (0) = g 0 (x j ),
(3.1)
where
m j = m(x j ),
m j + 1 = m(x j + 1 ), 2
2
and the two difference operators are
D−u j =
u j − u j −1
x
,
D+u j =
u j +1 − u j
x
.
In order to incorporate the zero flux boundary condition, we force
m 1 D + g 0 = m N + 1 D + g N = 0. 2
2
(3.2)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
661
Theorem 3.1. The semidiscrete scheme (3.1) possesses three properties: (1) If g j (0) 0 ∀ j, then
0 g j (t ) K
∀ j and ∀t > 0,
where K = max{1, max1 j N g j (0)}. (2) The semidiscrete entropy
E h (t ) =
N
m j g j log g j + (1 − g j ) x
j =1
is nonincreasing as time evolves. (3) It preserves the equilibrium solution g˜ ≡ 1; i.e., if g j (0) = 1, then
g j (t ) = 1,
1 j N ∀t > 0.
Proof. (1) Define
Σ = g : g j ∈ [0, K ], 1 j N ⊂ R N . is an equilibrium solution of the ODE system The claimed bounds follow if we show that Σ is an invariant region. Since 0 (3.1), it suffices to prove that m where
d g dt
· ν :=
N
m i νi
j =1
dg i
< 0 ∀ g ∈ ∂Σ \ {0 },
dt
ν is the outward normal vector at g . For each g , we define two sets of indices s0 and s K such that s0 = {1 j N ; g j = 0},
s K = {1 j N ; g j = K }.
Note that s0 and s K cannot be empty at the same time. There are two cases: (a) s0 = {1, . . . , N } and s0 = ∅. Then there exists j such that
g j = 0, and
νj =
g j +1 > 0 or
g j +1 = 0,
g j > 0,
−α j < 0, if j ∈ s0 , 0, if j ∈ / s0 .
Therefore,
m
d dt
g · ν = −
m j + 1 g j +1 + m j − 1 g j −1 2
2
(x)2
j ∈s0
(b) s0 = ∅, then s K = ∅. The normal vector
νj =
α j < 0.
ν is in the following form
α j > 0, if j ∈ s K , if j ∈ / sK .
0,
It follows that
m
d dt
g · ν =
1
(x)2
j ∈s K
m j + 1 ( g j +1 − K ) + m j − 1 ( g j −1 − K ) + λm2j K (1 − K ) 2
2
α j < 0.
So Σ is an invariant region. In other words, the scheme preserves the positivity of the solution and g j K ∀ j. (2) The time derivative of the entropy E h (t ) is given by
dE h (t ) dt
=
N
mj
j =1
=
N j =1
d dt
g j log g j x
D − (m j + 1 D + g j ) + λm2j g j (1 − g j ) log g j x 2
Author's personal copy 662
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
=−
N −1 j =1
m j + 1 D + g j (log g j +1 − log g j ) − λ 2
N
m2j g j ( g j − 1) log g j x
j =1
0, where we used the zero flux boundary condition and ( X − Y ) log YX 0 for X , Y > 0. (3) If g j (0) = 1 for 1 j N, then E h (0) = 0. So E h (t ) 0 by (2). However E h (t ) 0 ∀t 0. Hence E h (t ) = 0, implying that g j must be 1 for 1 j N. 2 3.2. Fully discrete scheme To preserve both positivity and entropy property at fully discrete level, we start with a simple time discretization:
mj
g nj +1 − g nj
t
= D − m j + 1 D + g nj +1 + λm2j g nj 1 − g nj +1 for 1 j N ,
(3.3)
2
subject to initial data
g 0j = g 0 (x j ),
1 j N,
and the zero flux boundary condition
m 1 D + g 0n+1 = m N + 1 D + g nN+1 = 0. 2
2
The scheme enjoys following properties. Theorem 3.2. The fully discrete scheme (3.3) is unconditionally stable in the sense that: (1) If g 0j 0 ∀ j, then for any n 1,
0 g nj K ,
1 j N,
where K = max{1, max1 j N g 0j }. (2) The fully discrete entropy
En =
N
m j g nj log g nj + 1 − g nj x
(3.4)
j =1
is nonincreasing, i.e., E n+1 E n ∀n 0. (3) It preserves the equilibrium solution g˜ ≡ 1; i.e. if g 0j = 1, then g nj = 1 ∀n 1. Proof. (1) Rearranging (3.3), we have
+1 +1 −tm j − 1 g nj − + t (m j − 1 + m j + 1 ) + (x)2m j 1 + λtm j g nj g nj +1 − tm j + 1 g nj + 1 1 2
2
2
2
= (x)2m j (1 + t λm j ) g nj .
(3.5)
This is a tridiagonal linear system with the coefficient matrix being diagonally dominant. So { g nj +1 } is nonnegative, as long as { g nj } is nonnegative. It also implies the existence and uniqueness of the numerical solution.
As to the upper bound, it suffices to show that if g nj K ∀ j, then g nj +1 K ∀ j. Let g ni +1 = max1 j N g nj +1 . Then the i-th row of the linear system gives
(x)2mi (1 + t λmi ) g ni (x)2mi 1 + t λmi g ni g in+1 . Therefore
g in+1 − g ni t λmi g ni 1 − g in+1 . This implies that g ni +1 max{1, g ni }. By the definition of K , we always have
g nj +1 K ,
1 j N.
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
663
(2) A simple calculation gives
E n +1 − E n =
N
j =1
=
N
m j g nj +1 log g nj +1 − g nj log g nj − g nj +1 − g nj x
mj
g nj +1 − g nj log g nj +1 + g nj log
j =1
g nj +1 g nj
− g nj +1 − g nj x.
Because m j > 0, g nj 0, and log X X − 1 for X > 0, we have
E
n +1
n
−E
N
mj
g nj +1
g nj
−
log g nj +1
+
g nj
j =1
=
N
g nj +1 g nj
−1 −
g nj +1
−
g nj
x
m j g nj +1 − g nj log g nj +1 x
j =1
= −t
N −1 j =1
m j+ 1
2
+1 g nj + 1
−
g nj +1
log
+1 g nj + 1
g nj +1
− t xλ
N
m2j g nj g nj +1 − 1 log g nj +1
j =1
0. Here we have used the facts that m 1 D + g 0n+1 = m N + 1 D + g nN+1 = 0 and ( X − Y ) log g nj
(3) Suppose N
2
2
X Y
0 for X , Y > 0.
≡ 1. We sum (3.5) over j and end up with
m j (1 + t λm j ) g nj +1 − 1 = 0.
j =1
From (1) we know that 0 g nj +1 1. Therefore g nj +1 = 1 for all j’s. So the scheme preserves the equilibrium solution.
2
Remark 3.1. The scheme (3.3) is easy to implement, yet it is only first order in time. One may easily construct a second order semi-implicit scheme, for example: find g nj +1 = 2g ∗j − g nj by solving
2m j
g ∗j − g nj
t
2 = D − m j + 1 D + g ∗j + λm2j g ∗j 1 − 2g nj + g nj 2
(3.6)
for 1 j N. However, this scheme is no longer unconditionally entropy stable. In practice, one might still prefer to use this scheme due to its second order accuracy in time. 4. Energy stable method for one-species system 4.1. Semidiscrete scheme Let g n (x) be an approximation of g (x, tn ) and G n = G (x, g n ). The initial data is given by g 0 (x) = g 0 (x). We first discretize (2.5) in time to obtain
M
g n +1 − g n
t
G n +1 − G n = ∇ · M ∇ g ∗ + n +1 g − gn
with
g∗ =
g n +1 + g n 2
.
Notice that
G n +1 − G n g n +1
−
gn
λ 2 = λM g ∗ m − M g ∗ − M 2 g ∗ − gn . 3
So (4.1) can be rewritten as a prediction–correction scheme:
(4.1)
Author's personal copy 664
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
(1) Given g n , compute the intermediate solution g ∗ which solves
⎧ ⎨ ⎩
2M
g ∗ − gn
λ
+ P g∗ +
t M g ν∗ = 0
3
M 2 g ∗ − gn
2
= 0 in Ω,
(4.2)
on ∂Ω.
(2) Obtain g n+1 from g n and g ∗ :
g n+1 = 2g ∗ − g n .
(4.3)
Theorem 4.1. Given g n 0 in Ω , if
t
0, hnj 0, 1 j N, then
g nj +1 > 0
and hnj +1 0
1 j N ∀n 0,
provided t < m1 . ∞ (2) The fully discrete entropy
E n = x
N
m j g nj + W j hnj − m j ln m j g nj
j =1
is nonincreasing unconditionally, i.e., E n+1 E n .
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
669
n+1 } depends on t continuously. If we increase t continuously from Proof. (1) Notice that the numerical solution { g n+1 , h 1 n + 1 zero up to m , no component of g can become negative without passing through zero since g n > 0. Assume for some ∞
t, we have g ni +1 = min1 j g nj +1 = 0. Then Eq. (5.4a) gives
−
g ni
t
0,
which contradicts with g ni > 0. So we must have g ni +1 > 0, hence g nj +1 > 0 ∀ j.
n+1 . Assume for some t, we have hn+1 = min1 j hn+1 = −ε , where We apply the similar argument to h i j small number. Then Eq. (5.4b) gives
ε is a sufficiently
hn − mi + mi g in+1 − W i ε ε − i 0, t t 1
1
t
− mi −mi g in+1 + W i ε 0,
since g ni +1 > 0, W i is bounded and
ε is sufficiently small. This contradicts with the choice of t. Therefore, we must have
hni +1 0. (2) For the fully discrete entropy E n , we have
E
n +1
n
−E −
N
n +1
j =1
+1 g nj +1 g nj + 1
α t M j+ 12 ( g j+1 − g j ) x
n 2
− t x
N
F nj +1
2
0.
2
j =1
6. Numerical tests We provide numerical results to demonstrate the accuracy of the schemes and the capacity to capture solution features for large times. Denote the initial data by u 0 (x), the numerical solution at (x j , tn ) by unj = M (x j ) g nj , and the exact solution by u (x j , tn ). Define the L ∞ error as
max unj − u (x j , tn ),
1 j N
and the L 1 error as N n u − u (x j , tn )h. j
j =1
When the exact solution is not available, we replace u (x j , tn ) by a reference solution to compute the errors. 6.1. Tests when P (m) = − ln m Consider the problem with Ω = [0, 1] and
m(x) = 10x2 (x − 1)2 (x − 0.4005)2 . Note that m(0) = m(1) = 0, so Eq. (2.3) becomes degenerate on ∂Ω . We first illustrate the accuracy of scheme (3.3) with the continuous initial data:
u 0 (x) =
sin(100x) + 1 2
(6.1)
.
Table 1 shows that the scheme is of second order in space. Here the reference solution is taken as the numerical solution with N = 1280. The entropy decreasing property is shown in Table 2, where the relative entropy (3.4) is decreasing to 0 as t increases, indicating the numerical convergence to the equilibrium solution m(x) as time becomes large. We also test the numerical convergence to the equilibrium solution m(x) by two numerical examples shown in Fig. 1. The initial data are respectively given by random data and the δ -like function
u 0 (x) =
10[cos(20π (x − 0.8)) + 1]
if |x − 0.8| 0.05;
0
elsewhere.
Author's personal copy 670
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
Table 1 Error and order of accuracy of the numerical solution of (3.3) with the initial data (6.1) on a uniform mesh of N cells, the final time t = 1.0. u 0 (x)
sin(100x)+1 2
N
L ∞ error
Order
L 1 error
Order
20 40 80 160
2.352E−03 2.601E−04 5.763E−05 1.524E−05
– 3.176 2.174 1.919
1.175E−03 1.305E−04 2.821E−05 6.800E−06
– 3.170 2.210 2.053
Table 2 The entropy of the numerical solutions on a uniform mesh with N = 80. t
0
5
10
50
100
200
400
E (t )
1.770E+00
1.534E−01
5.937E−02
2.702E−03
3.719E−04
1.905E−05
9.530E−08
Fig. 1. For the figures in the first row, u 0 (x) is given by the random data; for the second row, u 0 (x) is given by the δ -like function. N = 80.
6.2. Tests when P (m) = −αm Here
α 0 is a parameter. We consider the growth rate
m(x) = 10x2 (x − 1)2 (x − 0.4005)2 − 0.02, which is positive for some x ∈ Ω , and set up
m(x) dx ≈ −0.0048 < 0,
Ω
and
α = 40 and λ = 0.1 such that m(x)e − P dx ≈ 0.0045 > 0.
Ω
In such a setting we have λ∗ = 0: the solution will converge to the positive equilibrium other than m(x), regardless of the size of λ as well as the initial density. We first illustrate the accuracy of scheme (4.6)–(4.7) with Newton’s iteration. The initial data is taken as (6.1), again using the numerical solution with N = 1280 as the reference solution, the results in Table 3 show that the scheme is of the second order in space. We display in Table 4 the discrete energy V at different times, using the same initial data (6.1). Fig. 2 shows that solutions with different initial data converge to the same positive equilibrium solution, as predicted by the theory. Moreover, we observe that the equilibrium solution concentrates at the local positive maximum of m(x). This is consistent with the theoretical result obtained in [5, Theorem 4].
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
671
Table 3 Error and order of accuracy of the numerical solution of scheme (4.6)–(4.7) with initial data (6.1) on a uniform mesh of N cells, the final time t = 1.0. u 0 (x)
sin(100x)+1 2
N
L ∞ error
Order
L 1 error
Order
10 20 40 80 160
4.122E−03 2.010E−03 2.024E−04 4.354E−05 1.040E−05
– 1.036 3.312 2.217 2.066
6.491E−02 1.762E−02 4.156E−03 1.033E−03 2.552E−04
– 1.881 2.084 2.008 2.018
Table 4 The energy of the numerical solutions on a uniform mesh of N = 80 cells. tn
0
5
10
50
100
200
Vn
7.387E+02
4.405E−04
7.523E−05
6.435E−07
4.125E−08
−1.522E−08
Fig. 2. For the figures in the first row, u 0 (x) is given by the random data; for the second row, u 0 (x) is given by the δ -like function. N = 80.
6.3. Two-dimensional tests with P (m) = −αm2 We consider the resource
2 2 m(x) = 10 0.25 − (x1 − 0.5)2 − (x2 − 0.5)2 e (x1 −0.5) −(x2 −0.5) − 0.3 (see Fig. 3), and set
α = 4 so that
m(x) dx ≈ −0.1198 < 0 and Ω
m(x)e − P dx ≈ −0.1590 < 0.
Ω
In such a setting we have λ∗ > 0: the solution will converge to the positive equilibrium if λ > λ∗ , and decay to zero if λ λ∗ , regardless of the choice of initial density. We test scheme (4.6)–(4.7) with two initial data: (i) the cos-like function u 0 (x) = cos(3π |x|2 ) + 1; (ii) the δ -like function concentrating at (x0i , x0 j ) for i , j = 1, 2 with x01 = 0.2 and x02 = 0.8.
u 0 (x) =
25[cos(10π (x1 − x0i )) + 1][cos(10π (x2 − x0 j )) + 1]
if |xi − x0 j | 0.1 and i , j = 1, 2;
0
elsewhere.
Author's personal copy 672
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
2 2 Fig. 3. Plot of m(x) = 10[0.25 − (x1 − 0.5)2 − (x2 − 0.5)2 ]e (x1 −0.5) −(x2 −0.5) − 0.3.
Fig. 4. Figures in the first row, u 0 (x) is given by the cos-like function; in the second row, u 0 (x) is given by the δ -like function. N = 20 × 20, λ = 1 < λ∗ .
From Fig. 4, we observe that when t = 100, the numerical solutions is of order O(10−7 ), approaching zero. In other words, the species will extinct eventually, consistent with the theoretical result stated in Theorem 1.1. Approach toward the positive equilibrium is shown with two examples in Fig. 5, corresponding to the persistence of the species as predicted by Theorem 1.1. Moreover, the energy V ≈ −0.13439 when t = 0.5 and 100 for both of the initial data. These tests indicate that the steady state is asymptotically approached in time and well preserved. 6.4. Two-dimensional test for the two-species system We consider the two-species system with resource potentials P (m) = − ln m and Q (m) = −m, where the resource m(x) is given by
m(x) = 0.5 sin(π x1 ) sin(π x2 ) + 0.6, (See Fig. 6.) We set
x ∈ Ω = [0, 2.5]2 .
α = 5, β = 1, and take initial value as
u 0 (x) = 1 − 2v 0 (x),
v 0 (x) =
1 4
sin2 (4π x1 ) sin2 (4π x2 ) if 1 x1 , x2 1.5,
0
elsewhere.
Visually, we observe the approach toward the equilibrium solution (m, 0) through the numerical result shown in Fig. 7. Quantitatively, Table 5 shows that the entropy is decreasing as time evolves and indicates the asymptotic approach toward the equilibrium state, while we note that E ≈ 5.34191 when the solution approaches the steady state. 7. Concluding remarks We have developed finite difference schemes for a class of reaction–diffusion–advection equations arising in the evolution of biased dispersal of population dynamics, including the single species model and two competing species system. The schemes are shown to preserve both positivity and equilibrium, satisfying either entropy or energy dissipation structure.
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
673
Fig. 5. Figures in the first row, u 0 (x) is given by the cos-like function; in the second row, u 0 (x) is given by the δ -like function. N = 20 × 20, λ = 1000 > λ∗ .
Fig. 6. Plot of m(x) = 0.5 sin(π x1 ) sin(π x2 ) + 0.6. Table 5 The discrete entropy on a uniform mesh with N = 40 × 40. tn
0
0.1
0.3
0.5
100
En
6.26436
5.80962
5.66577
5.57095
5.35172
Fig. 7. The first row: u component; the second row: v component. N = 40 × 40.
Author's personal copy 674
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
These ensure that the numerical solution provides a satisfying long-time behavior, thus revealing the desired dispersal pattern. Numerical examples are shown to illustrate the second-order accuracy for various resource potentials and underline the efficiency to preserve the large-time asymptotic. We also proved that the unique existence of the numerical solutions for the proposed schemes. Our future work would develop higher order methods based on the proposed semidiscrete formulations. For dispersal models with multi-species, the dynamics become much more complex, further numerical investigation into multi-species dispersal models deserves serious consideration. Acknowledgements This research was partially supported by the National Science Foundation under grant DMS-1312636 and the Mathematical Research Network KI-net grant. Appendix A Theorem A.1. Problem (2.10) has a unique λ∗ > 0 characterized by a positive eigenfunction if and only if
mM dx < 0. Ω
Moreover, λ∗ = 0 if and only if Ω mM dx 0. Proof. The existence of (λ∗ , ψ) with ψ > 0 in Ω and M = 1 is given in [6]. In a similar manner, one can prove the existence of (λ∗ , ψ) with ψ > 0 and nonconstant M > 0. Also, λ∗ can be expressed by the variational formulation
λ∗ = inf
2 Ω mM v dx
v∈S
where
∇ v 2M
,
S = v ∈ H 1 , M v ν = 0 on ∂Ω and
mM v 2 dx > 0 .
Ω
In fact, (2.10) yields
∇ψ 2M 0. 2 Ω mM ψ dx
λ∗ =
(A.1)
Divide both sides of (2.10) by ψ , and integrate over Ω to obtain
λ∗
mM dx = −
Ω
M |∇ψ|2
ψ2
dx 0.
(A.2)
Ω
For Ω mM dx < 0, (A.2) implies that λ∗ > 0. For Ω mM dx = 0, (A.2) implies that ∇ψ 2M = 0. Therefore λ∗ = 0 by (A.1). For Ω mM dx > 0, we must have λ∗ = 0 by (A.2). 2 Let (σ ∗ , φ) denote the pair of the principle eigenvalue and eigenfunction of (2.9), we thus have φ > 0 in Ω , see [38, Theorem 11.10]. The following result relates the sign of σ ∗ to the relative size of λ in terms of λ∗ . Theorem A.2. sign(σ ∗ ) = sign(λ − λ∗ ). Proof. Rewrite (2.9) as
∇ · ( M ∇φ) + λ m −
σ∗ λ
M φ = 0,
∗ = m − σλ . Then λ is determined by Let m
∇ v 2M , 2 Ω m M v dx
λ = inf v ∈ S˜
M φν |∂Ω = 0.
(A.3)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
˜ If where S˜ is a modification of S with m replaced by m. < m and eigenvalue. If σ ∗ > 0, then m
675
σ ∗ = 0, then φ = ψ and λ = λ∗ by the simpleness of the principle
M φ dx < m 2
0< Ω
mM φ 2 dx. Ω
Therefore
∇φ 2M ∇φ 2M > λ∗ . 2 dx 2 dx m M φ mM φ Ω Ω
λ= If
> m and σ ∗ < 0, then m 2 M ψ dx > mM ψ 2 dx > 0. m Ω
Ω
Hence,
∇ψ 2M ∇ψ 2M < = λ∗ . 2 dx 2 dx m M ψ mM ψ Ω Ω
λ
2
The result about the equilibrium solution g˜ is the following.
¯ is positive somewhere in Ω and P is bounded. Theorem A.3. Suppose that m ∈ C 2 (Ω) (1) If 0 < λ λ∗ , (2.7) only has a zero solution. (2) If λ > λ∗ , (2.7) has a zero solution and a unique positive solution. Proof. (1) For 0 < λ λ∗ , multiply (2.7) by g˜ to get
mM g˜ 2 dx =
λ Ω
M 2 g˜ 2 dx 0.
M |∇ g˜ |2 dx + λ Ω
Ω
By the definition of λ∗ , we have
λ
mM g˜ 2 dx λ∗
Ω
Ω
M g˜ dx λ − λ 2 2
λ Ω
M 2 g˜ 2 dx. Ω
It follows that
mM g˜ 2 dx + λ
∗
mM g˜ 2 dx 0. Ω
Therefore g˜ ≡ 0 on Ω . m(x) (2) One can verify that any constant g¯ greater than M (x) ∞ is a super-solution of (2.7). We normalize the positive eigenfunction φ such that φ ∞ = 1. Then g = ε φ with
0 < ε min
σ∗ λ M ∞
, g¯
is a sub-solution of (2.7) because
P g = ε M φ ε λ M φ − σ ∗ ε M φ ε λ M − σ ∗ 0 in Ω Moreover,
0 < g g¯ and
and
M g |∂Ω = 0. ν
∂ λ F ( g , m) = −λ M 2 < 0. ∂g g
So the conditions of [37, Theorem 3.4 in Chapter 3] are satisfied, and it ensures that there exists a unique solution g˜ of (2.7) such that g g˜ g¯ . The above claimed uniqueness is actually restricted to the situation that any two solutions, if g˜ 1 g˜ 2 , must coincide. For the general case, we define the minimum of two arbitrary positive solutions as
Author's personal copy 676
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
g ∗ (x) = min g˜ 1 (x), g˜ 2 (x) . A slight modification of the argument in [21, Lemma 1.5] when applied to (2.7) (care is needed for treating the Neumann boundary condition) enables us to show that g ∗ is a weak super-solution. Subtracting the weak formulations of g ∗ and g˜ i against the test function taken as g˜ i and g ∗ respectively, we have
λ
M 2 g ∗ g˜ i g˜ i − g ∗ dx 0 for i = 1, 2,
Ω
which implies that g ∗ = g˜ i in Ω , hence g˜ 1 = g˜ 2 .
2
The above results when combined with the nonincreasing property of the energy V defined in (2.4) enable us to predict persistence or extinction of the evolving species; see [7] for more background details. Theorem A.4. Under the same assumptions as in Theorem A.3, we have: (1) If 0 < λ λ∗ , all nonnegative solutions of (2.1) decay toward zero as t → ∞. (2) If λ > λ∗ , the positive equilibrium is globally attractive among nonzero nonnegative solutions. Proof. We first prove the convergence to the steady state. The discussion on V in Section 2.2 shows that
d dt
M (∂t g )2 dx 0.
V =−
(A.4)
Ω
Thus V is nonincreasing along the orbit starting at g 0 (x). The boundedness and precompactness of the orbits imply that V [ g ] is bounded from below, so there exists V eq = limt →∞ V [ g ]. If w is in the ω -limit set of g 0 (x), then V [ w ] = V eq . Let w ∗ (x, t ) be the orbit starting at some point w in the ω -limit set of g 0 (x). Since that set is invariant, the entire orbit w ∗ (x, t ) must belong to it, so V [ w ∗ ] = V eq . In view of (A.4), this is possible if and only if ∂t w ∗ = 0. Therefore, w must be an equilibrium of (2.1). Since w is an arbitrary element of the ω -limit set of an arbitrary initial point g 0 , the ω -limit set for the semi-flow generated by (2.1) must consist of all the equilibria, which shows the global attractiveness of stable equilibrium solutions. We now distinguish two cases in terms of λ. The claimed results follow from the fact that zero equilibrium solution is stable for 0 < λ λ∗ and unstable for λ > λ∗ , as inferred from Theorems A.2 and A.3. 2 References [1] E. Celledoni, V. Grimm, R.I. McLachlan, D.I. McLaren, D. O’Neale, B. Owren, G.R.W. Quispel, Preserving energy resp. dissipation in numerical PDEs using the “Average Vector Field” method, arXiv:1202.4555, 2012. [2] Fethi Belgacem, Chris Cosner, The effects of dispersal along environmental gradients on the dynamics of populations in heterogeneous environments, Can. Appl. Math. Q. 3 (4) (1995) 379–397. [3] R. Bermejo, J. Carpio, An adaptive finite element semi-Lagrangian implicit–explicit Runge–Kutta–Chebyshev method for convection dominated reaction– diffusion problems, Appl. Numer. Math. 58 (1) (January 2008) 16–39. [4] Rodolfo Bermejo, Mofdi El Amrani, A finite element semi-Lagrangian explicit Runge–Kutta–Chebyshev method for convection dominated reaction– diffusion problems, J. Comput. Appl. Math. 154 (1) (May 2003) 27–61. [5] Andriy Bezuglyy, Yuan Lou, Reaction–diffusion models with large advection coefficients, Appl. Anal. 89 (7) (2010) 983–1004. [6] K.J. Brown, S.S. Lin, On the existence of positive eigenfunctions for an eigenvalue problem with indefinite weight function, J. Math. Anal. Appl. 75 (1) (1980) 112–120. [7] Robert Stephen Cantrell, Chris Cosner, Spatial Ecology via Reaction–Diffusion Equations, John Wiley & Sons, Ltd., 2004. [8] Robert Stephen Cantrell, Chris Cosner, Donald L. Deangelis, Victor Padron, The ideal free distribution as an evolutionarily stable strategy, J. Biol. Dyn. 1 (3) (2007) 249–271. [9] Robert Stephen Cantrell, Chris Cosner, Yuan Lou, Movement toward better environments and the evolution of rapid diffusion, Math. Biosci. 204 (2) (2006) 199–214. [10] Robert Stephen Cantrell, Chris Cosner, Yuan Lou, Approximating the ideal free distribution via reaction–diffusion–advection equations, J. Differ. Equ. 245 (12) (2008) 3687–3703. [11] Robert Stephen Cantrell, Chris Cosner, Yuan Lou, Evolution of dispersal and the ideal free distribution, Math. Biosci. Eng. 7 (1) (2010) 17–36. [12] J.A. Carrillo, M.P. Gualdani, A. Jüngel, Convergence of an entropic semi-discretization for nonlinear Fokker–Planck equations in Rd , Publ. Mat. 52 (2) (2008) 413–433. [13] Xinfu Chen, Richar Hambrock, Yuan Lou, Evolution of conditional dispersal: a reaction–diffusion–advection model, J. Math. Biol. 57 (3) (2008) 361–386. [14] Q. Chang, E. Jia, W. Sun, Difference schemes for solving the generalized nonlinear Schrödinger equation, J. Comput. Phys. 148 (1999) 397–415. [15] Chris Cosner, Yuan Lou, Does movement toward better environments always benefit a population?, J. Math. Anal. Appl. 277 (2) (2003) 489–503. [16] Jérôme Coville, Juan Dávila, Salomé Martínez, Nonlocal anisotropic dispersal with monostable nonlinearity, J. Differ. Equ. 244 (12) (2008) 3080–3118. [17] Feng Chen, Jie Shen, Efficient energy stable schemes with spectral discretization in space for anisotropic Cahn–Hilliard systems, Commun. Comput. Phys. 13 (2013) 1189–1208. [18] Craig Collins, Jie Shen, Steven M. Wise, An efficient, energy stable scheme for the Cahn–Hilliard–Brinkman system, Commun. Comput. Phys. 13 (2013) 929–957. [19] M. Dahlby, B. Owren, A general framework for deriving integral-preserving numerical methods for PDEs, SIAM J. Sci. Comput. 33 (2011) 2318–2340.
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
677
[20] David J. Eyre, Unconditionally gradient stable time marching the Cahn–Hilliard equation, in: Computational and Mathematical Models of Micro Structural Evolution, San Francisco, CA, 1998, in: Mater. Res. Soc. Symp. Proc., vol. 529, MRS, Warrendale, PA, 1998, pp. 39–46. [21] E.N. Dancer, G. Sweers, On the existence of a maximal weak solution for a semilinear elliptic equation, Differ. Integral Equ. 2 (4) (1989) 533–540. [22] Jack Dockery, Vivian Hutson, Konstantin Mischaikow, Mark Pernarowski, The evolution of slow dispersal rates: a reaction diffusion model, J. Math. Biol. 37 (1) (1998) 61–83. [23] Daisuke Furihata, Finite difference schemes for ∂ u /∂ t = (∂/∂ x)α δ G /δ u that inherit energy conservation or dissipation property, J. Comput. Phys. 156 (1) (1999) 181–205. [24] Richard Gejji, Yuan Lou, Daniel Munther, Justin Peyton, Evolutionary convergence to ideal free dispersal strategies and coexistence, Bull. Math. Biol. 74 (2) (2011) 1–43. [25] W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion–Reaction Equations, first edition, Springer Ser. Comput. Math., vol. 33, Springer, 2003. [26] Chiu-Yen Kao, Yuan Lou, Wenxian Shen, Random dispersal vs. non-local dispersal, Discrete Contin. Dyn. Syst. 26 (2) (2010) 551–596. [27] King-Yeung Lam, Concentration phenomena of a semilinear elliptic equation with large advection in an ecological model, J. Differ. Equ. 250 (1) (2011) 161–181. [28] King-Yeung Lam, Wei-Ming Ni, Limiting profiles of semilinear elliptic equations with large advection in population dynamics, Discrete Contin. Dyn. Syst. 28 (3) (2010) 1051–1067. [29] D. Lanser, J.G. Verwer, Analysis of operator splitting for advection–diffusion–reaction problems from air pollution modeling, J. Comput. Appl. Math. 111 (1–2) (1999) 201–216. [30] H. Liu, H. Yu, An entropy satisfying conservative method for the Fokker–Planck equation of the finitely extensible nonlinear elastic dumbbell model, SIAM J. Numer. Anal. 50 (3) (2012) 1207–1239. [31] Christian Lubich, Alexander Ostermann, Runge–Kutta approximation of quasi-linear parabolic equations, Math. Comput. 64 (210) (1995) 601–627. [32] S. Li, L. Vu-Quoc, Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein–Gordon equation, SIAM J. Numer. Anal. 32 (1995) 1839–1875. [33] Takayasu Matsuo, Dissipative/conservative Galerkin method using discrete partial derivatives for nonlinear evolution equations, J. Comput. Appl. Math. 218 (2) (2008) 506–521. [34] Takayasu Matsuo, Hisashi Yamaguchi, An energy-conserving Galerkin scheme for a class of nonlinear dispersive equations, J. Comput. Phys. 228 (12) (2009) 4346–4358. [35] J.D. Murray, Spatial models and biomedical applications, in: Mathematical Biology. II, third edition, in: Interdiscip. Appl. Math., vol. 18, Springer-Verlag, New York, 2003. [36] Akira Okubo, Diffusion and Ecological Problems: Mathematical Models, Biomathematics, vol. 10, Springer-Verlag, Berlin, 1980, an extended version of the Japanese edition. Ecology and Diffusion, translated by G.N. Parker. [37] Chia-Ven Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, 1992. [38] Joel Smoller, Shock Waves and Reaction–Diffusion Equations, Grundlehren Math. Wiss. (Fundam. Princ. Math. Sci.), vol. 258, Springer-Verlag, New York, 1983.
Author's personal copy Journal of Computational Physics 256 (2014) 656–677
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Entropy/energy stable schemes for evolutionary dispersal models Hailiang Liu, Hui Yu Iowa State University, Mathematics Department, Ames, IA 50011, United States
a r t i c l e
i n f o
Article history: Received 4 March 2013 Received in revised form 13 August 2013 Accepted 17 August 2013 Available online 17 September 2013 Keywords: Dispersal Diffusion–advection–reaction equation Entropy Energy Positivity preserving
a b s t r a c t In this paper we propose some entropy/energy stable finite difference schemes for the reaction–diffusion–advection equation arising in the evolution of biased dispersal of population dynamics. The peculiar feature of these active dispersal models is that the transient solution converges to the stable steady state when time goes to infinity. For the numerical method to capture the long-time pattern of persistence or extinction, we use the relative entropy when the resource potential is logarithmic, and explore the usual energy for other resource potentials. The present schemes are shown to satisfy three important properties of the continuous model for the population density: (i) positivity preserving, (ii) equilibrium preserving, and (iii) entropy or energy satisfying. These ensure that our schemes provide a satisfying long-time behavior, thus revealing the desired dispersal pattern. Moreover, we present several numerical results which confirm the second-order accuracy for various resource potentials and underline the efficiency to preserve the large time asymptotic. © 2013 Elsevier Inc. All rights reserved.
1. Introduction This work is concerned with the numerical approximation of a class of reaction–diffusion–advection equations arising in the evolution of biased dispersal of population dynamics, with emphasis on exploring the entropy/energy structure of the dispersal model so that the resulting methods provide a satisfying long-time behavior. 1.1. Mathematical formulations Reaction–diffusion equations have been widely used to model the biological problems [35,36]. One of the well-known examples is the logistic reaction–diffusion model for the population growth with random dispersal,
∂t u = u + λu (m − u ) in Ω × (0, ∞),
(1.1)
where the population inhabiting a bounded domain Ω ⊂ Rd has density u (x, t ) at location x and time t, and local growth rate m. The parameter λ > 0 is the reverse of the dispersal rate. If the environment is heterogeneous, i.e., m(x) is not a constant, then the population may have a tendency to move toward resources in addition to the random movement. The model may be upgraded to the following form
∂t u = ∇ · (∇ u + u ∇ P ) + λu (m − u ) in Ω × (0, ∞), E-mail addresses: [email protected] (H. Liu), [email protected] (H. Yu). 0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.08.032
(1.2)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
657
where P = P (m), which we call resource potential, reflects the movement tendency of the population. The time evolution is subjected to both the initial density u (x, 0) = u 0 (x) and the zero flux boundary condition
(∇ u + u ∇ P ) · ν = 0 on ∂Ω × (0, ∞),
(1.3)
where ν is the outward normal vector on the boundary ∂Ω which is assumed to be smooth. Several dispersal strategies have been studied in literature, for example P = −αm in [2,5,9,13,15], and the obtained results may apply to a more general reaction than
F (u , m) = u (m − u ), as long as it satisfies F (0, m) = 0 and F (+∞, m) < 0. With the logistical reaction term or its variation, the biologically relevant solutions to (1.2) are nonnegative and ultimately bounded. It follows from the regularity theory of parabolic partial differential equations that in the state space bounded orbits are pre-compact, hence the semi-flow will have a compact attractor, and the ω -limit set of any initial state u 0 (x) will be a compact invariant set. In other words, as time evolves the solution of (1.2) is expected to approach some stable patterns, independent of the choice of initial density; see, e.g., [7]. The main result reviewed in Section 2 may be stated as follows.
¯ is positive somewhere in Ω and P is smooth in m. There exists a unique λ∗ > 0 and a positive Theorem 1.1. Suppose that m ∈ C 2 (Ω) equilibrium solution to (1.2) if and only if
me − P dx < 0.
Ω
Moreover, λ∗ = 0 if and only if Ω me − P dx 0. (1) If 0 < λ λ∗ , all nonnegative solutions of (1.2) decay toward zero as t → ∞. (2) If λ > λ∗ , the positive equilibrium is globally attractive among nonzero nonnegative solutions. Of special interest is the ideal free dispersal strategy determined by P such that at an equilibrium u eq (x) [8,10,11], both
∇ · (∇ u eq + u eq ∇ P ) = 0 and m − u eq = 0 in Ω. This will hold if P = − ln m (unique up to some constant). In such a setting, the species u can perfectly match the environmental resource. The role of ideal free distributions has been well recognized in literature. Moreover, the quantity
E [u ] =
u log
u m
+ m − u dx
(1.4)
Ω
is nonincreasing in time. This corresponds to some physical entropy relative to the equilibrium state u eq (x) = m(x), called relative entropy. For general resource potential P , the quantity
V [u ] =
1 − P P 2 e ∇ ue − λe P 2
u
F (ξ, m) dξ dx
(1.5)
0
Ω
is nonincreasing in time. Moreover, Eq. (1.2) can be rewritten as a gradient flow,
∂t u = −e − P
δV , δu
(1.6)
where δδVu denotes the standard variational derivative of V with respect to u, and the functional V often corresponds to some free energy of the underlying physics, and thus is called the energy. In the context of two competing species, the evolutionary dynamics of conditional dispersal becomes much more complex, see [13,22]. A typical model may be described as
⎧ ut = α ∇ · (∇ u + u ∇ P ) + u (m − u − v ) ⎪ ⎪ ⎨ v = β∇ · (∇ v + v ∇ Q ) + v (m − u − v ) t ⎪ u + u ∂ ν P = ∂ν v + v ∂ ν Q = 0 ∂ ⎪ ⎩ ν u (x, 0) = u 0 (x), v (x, 0) = v 0 (x)
in Ω, in Ω, on ∂Ω,
(1.7)
in Ω.
The two nonnegative constants α and β represent the dispersal rates of two species, respectively. In this paper we shall focus on the case when α = β , and we can envision that system (1.7) models two competing species that are identical except their different dispersal strategies P and Q . An interesting question is whether there is any strategy in system (1.7)
Author's personal copy 658
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
that is globally evolutionarily stable, in the sense that the semi-trivial steady state (m, 0) is always locally stable. This is the case if the dispersal strategy is ideal free at (m, 0), i.e., P (m) = − ln m + C and Q (m) = P (m); see [11,24]. Moreover, the following functional
E [u , v ] =
u dx + Ω
v dx − Ω
m(x) ln u dx
(1.8)
Ω
is nonincreasing in time. For other results on dispersal strategies in heterogeneous landscapes, we refer to [13,27,28] for the study of aggregation profiles in terms of the resource distribution, and [16,26] for nonlocal dispersal strategies. The aim of this paper is to give reliable numerical schemes for (1.2) and (1.7) from the perspective of providing a satisfying long-time behavior. A key fact is that they both admit certain entropy/energy structure, and we demand our numerical schemes to satisfy the entropy/energy decreasing property in the discrete setting. In addition, both positivity and steady states for (1.1) are required to be preserved as well. These three requirements together are important for system (1.1), yet adding levels of difficulty to the design of a numerical method with high accuracy. As a preliminary attempt, only some simple time stepping is discussed in the present paper. In recent years, energy preservation or dissipation numerical schemes have drawn much interest and been extensively studied for several PDEs, since they are more likely to give a better performance in the long-time simulation; see e.g. [1, 12,23,30,33,34]. To the best of our knowledge, however, no entropy stable scheme for (1.2) or (1.7) has yet appeared in the literature. Though, due to the physical and mathematical significance of advection–diffusion–reaction equations, numerical methods have been developed and analyzed through decades; see, e.g., [3,4,29,31] and the recent monograph [25]. If the model equation has a nice property such as energy conservation/dissipation, then design of an energy satisfying numerical method is often desired. The Crank–Nicolson time discretization has been proven successful in many cases. For example, in [14], eight finite difference schemes are compared for solving the nonlinear Schrödinger equation, and two of three conservative schemes are of the Crank–Nicolson type. Efforts have also been devoted to the development of procedures to mimic the continuous energy property, see e.g., [1,19,23,32]. For example, the average vector field method is developed in [1] for some Hamiltonian PDEs and dissipative PDEs. In an earlier work by Furihata [23], a procedure was given to derive energy satisfying (conservation or dissipation) finite difference schemes for the class of ∂t u = ∂xα ( δδGu ), about which the
d G = 0 or 0. For higher order nonlinear dissipative discrete calculus is used to mimic the continuous energy property dt PDEs such as the Cahn–Hilliard equation, efficient time stepping techniques which can inherit the unconditional energy stability are particularly in demand. Several methods such as the convex splitting by Eyre [20] and the stabilized time discretization which allow large time steps have been explored in literature, see e.g., [17,18] and references therein.
1.2. Main steps of this paper In Section 2 mathematical theory and properties regarding the target equation (1.2) are presented; we particularly highlight the energy/entropy structure, as well as the large time pattern formation in terms of the parameter λ. The key idea is to reformulate (1.2) and (1.7) into the nonlogarithmic Landau form in terms of g = ue P and h = ve Q , following the strategy in [30] by the authors for the Fokker–Planck equation of FENE dumbbell polymers. Section 3 is devoted to the case with the logarithmic potential and the proposed scheme is presented. The schemes at both semidiscrete and fully discrete levels are shown to satisfy three important properties of the continuous model for the population density: (i) positivity preserving, (ii) equilibrium preserving, and (iii) entropy satisfying. The schemes presented in Section 4 apply to the case with general resource potential and is extensible to handle more general reaction terms. In that section we first present a semidiscrete scheme with discretization only in time. This scheme is proved to be uniquely solvable provided the time step satisfies a bound depending on the growth rate m and the numerical solution at the previous time step. The scheme preserves the equilibrium solution with numerical energy being nonincreasing in time for arbitrary time step. In addition, the scheme can be expressed as a second-order prediction–correction method, and the prediction is shown to be in the same range of the numerical solution at the previous time step. As a consequence, we prove that the fully discrete scheme satisfies both desired properties: energy satisfying and equilibrium preserving. We also propose an entropy stable scheme in Section 5 to simulate the ideal free distribution for two competing species. Section 6 is devoted to numerical tests of proposed schemes. Finally some concluding remarks are presented in Section 7. 2. Model reformulation and mathematical properties In order to reveal the dynamic picture of (1.2), following [30], we rewrite (1.2) in terms of g = u / M with M = e − P :
M ∂t g + P g = 0 in Ω,
(2.1)
where the operator P is defined as
P g = −∇ · ( M ∇ g ) − λ M g (m − M g ).
(2.2)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
The initial condition is g 0 (x) = M g (m − M g ).
u 0 (x) M (x)
659
and on the boundary ∂Ω we have M g ν = 0. The reaction term becomes F ( g , m) =
Remark 2.1. M is allowed to be zero on boundary, which corresponds to the case that P is singular in the original equation. This is an important case in biological dispersal with an inapproachable boundary. Our reformulation has advantage to convert this singularity to certain degeneracy in term M ∂t g, which can be easily treated with some implicit time stepping, as done in [30] for the FENE polymer model. 2.1. Entropy for log potential P (m) = − ln m In such case M = e − P = m and (2.1) reduces to
m∂t g = ∇ · (m∇ g ) + λm2 g (1 − g )
in Ω.
(2.3)
The logistic reaction term ensures that any solution remains nonnegative and bounded, i.e.,
0 g (x, t ) max 1, g 0 ∞
∀(x, t ) ∈ Ω × [0, ∞).
As t → ∞, g will converge to the stable equilibrium solution g˜ ≡ 1, independent of the choice of the initial data. This motivates us to explore the entropy which might have some biological significance. The entropy defined in (1.4) becomes
E [g] =
m g log g + (1 − g ) dx. Ω
As a convex functional of g, it defines the entropy structure of (2.3):
m∂t g = ∇ · mg ∇ m−1
δE δg
+ λ F ( g , m),
where the variational derivative δδ Eg = m log g. The reaction term F ( g , m) = m2 g (1 − g ) satisfies
F·
δE 0 δg
for nonnegative g. Using integration by parts, we have
d dt
E =−
2 δE δ E mg ∇ m−1 dx + λ m −1 F · dx 0, δg δg
Ω
Ω
where the zero flux boundary condition has been taken into account. We shall propose a numerical method to satisfy this entropy property. 2.2. Energy for all other potential In general, the evolution may lead to extinction or persistence of the species, and the concept of entropy is no longer adequate to control the reaction term. We shall follow the standard gradient flow idea [7] by using the functional defined in (1.5),
V [g] =
1 2
M |∇ g | − G x, g (x, t ) dx, 2
(2.4)
Ω
where G is chosen as
g G (x, g ) = λ
F (s, m) ds =
λ 2
m(x) M (x) g 2 −
λ 3
M 2 (x) g 3 ,
0
so that (2.1) can be rewritten as
M ∂t g = ∇ · ( M ∇ g ) + ∂ g G . The above energy also defines the gradient flow structure of (2.5):
(2.5)
Author's personal copy 660
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
M ∂t g = −
δV , δg
where the zero flux has been used in the variational derivative. A formal calculation gives
d dt
δV ∂t g dx = − δg
V = Ω
M (∂t g )2 dx 0.
(2.6)
Ω
We shall propose a numerical method to satisfy this energy dissipation inequality in the discrete setting. 2.3. Large time pattern formation The nonnegative equilibrium solution g˜ (x), if it exists, solves
P g˜ = 0 in Ω
M g˜ ν |∂Ω = 0.
and
Let g = g˜ + ε w in (2.1) and then send
(2.7)
ε → 0 to obtain the linearized problem
M ∂t w = ∇ · ( M ∇ w ) + λ M w (m − 2M g˜ )
in Ω
M w ν |∂Ω = 0.
and
(2.8)
It is known that if g˜ is linearly stable, then g˜ is stable; for a proof see, e.g., [38, Theorem 11.22]. The linear stability of g˜ is determined by the spectrum of the linearized operator of P . For g˜ = 0, it boils down to the eigenvalue problem:
∇ · ( M ∇φ) + λmM φ = σ M φ in Ω and M φν |∂Ω = 0. (2.9) When σ < 0, 0 is stable; when σ > 0, 0 is unstable. Hence σ = 0 is a threshold value, corresponding to a special parameter λ∗ which satisfies
∇ · ( M ∇ψ) + λ∗mM ψ = 0 in Ω and M ψν |∂Ω = 0.
(2.10)
In what follows we shall denote the weighted L 2 norm by v M with
v 2M =
M v 2 dx. Ω
The basic results for (2.1), (2.7)–(2.10), leading to Theorem 1.1, are presented in Appendix A. These results are known for some specific choices of P such as P = −αm [2,7], yet a complete account for general resource potentials seems unavailable in literature. For the sake of completeness, we present detailed results for (2.1) together with the key arguments in Appendix A. 3. Entropy stable method for one-species model 3.1. Semidiscrete scheme in space We present only the one-dimensional case. Generalizations to the multidimensional space are straightforward for tensor product grids and the results remain valid without modifications. Given an integer N, we partition the domain Ω := [a, b] by using a uniform mesh x, where
b−a
x =
N
and
xj = a +
j−
1 2
x.
Note that x 1 = a and x N + 1 = b. Let g j (t ) be the approximation of g (x j , t ). We propose the following semidiscrete scheme: 2
2
for 1 j N,
mj
d dt
g j = D − (m j + 1 D + g j ) + λm2j g j (1 − g j ) 2
and
g j (0) = g 0 (x j ),
(3.1)
where
m j = m(x j ),
m j + 1 = m(x j + 1 ), 2
2
and the two difference operators are
D−u j =
u j − u j −1
x
,
D+u j =
u j +1 − u j
x
.
In order to incorporate the zero flux boundary condition, we force
m 1 D + g 0 = m N + 1 D + g N = 0. 2
2
(3.2)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
661
Theorem 3.1. The semidiscrete scheme (3.1) possesses three properties: (1) If g j (0) 0 ∀ j, then
0 g j (t ) K
∀ j and ∀t > 0,
where K = max{1, max1 j N g j (0)}. (2) The semidiscrete entropy
E h (t ) =
N
m j g j log g j + (1 − g j ) x
j =1
is nonincreasing as time evolves. (3) It preserves the equilibrium solution g˜ ≡ 1; i.e., if g j (0) = 1, then
g j (t ) = 1,
1 j N ∀t > 0.
Proof. (1) Define
Σ = g : g j ∈ [0, K ], 1 j N ⊂ R N . is an equilibrium solution of the ODE system The claimed bounds follow if we show that Σ is an invariant region. Since 0 (3.1), it suffices to prove that m where
d g dt
· ν :=
N
m i νi
j =1
dg i
< 0 ∀ g ∈ ∂Σ \ {0 },
dt
ν is the outward normal vector at g . For each g , we define two sets of indices s0 and s K such that s0 = {1 j N ; g j = 0},
s K = {1 j N ; g j = K }.
Note that s0 and s K cannot be empty at the same time. There are two cases: (a) s0 = {1, . . . , N } and s0 = ∅. Then there exists j such that
g j = 0, and
νj =
g j +1 > 0 or
g j +1 = 0,
g j > 0,
−α j < 0, if j ∈ s0 , 0, if j ∈ / s0 .
Therefore,
m
d dt
g · ν = −
m j + 1 g j +1 + m j − 1 g j −1 2
2
(x)2
j ∈s0
(b) s0 = ∅, then s K = ∅. The normal vector
νj =
α j < 0.
ν is in the following form
α j > 0, if j ∈ s K , if j ∈ / sK .
0,
It follows that
m
d dt
g · ν =
1
(x)2
j ∈s K
m j + 1 ( g j +1 − K ) + m j − 1 ( g j −1 − K ) + λm2j K (1 − K ) 2
2
α j < 0.
So Σ is an invariant region. In other words, the scheme preserves the positivity of the solution and g j K ∀ j. (2) The time derivative of the entropy E h (t ) is given by
dE h (t ) dt
=
N
mj
j =1
=
N j =1
d dt
g j log g j x
D − (m j + 1 D + g j ) + λm2j g j (1 − g j ) log g j x 2
Author's personal copy 662
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
=−
N −1 j =1
m j + 1 D + g j (log g j +1 − log g j ) − λ 2
N
m2j g j ( g j − 1) log g j x
j =1
0, where we used the zero flux boundary condition and ( X − Y ) log YX 0 for X , Y > 0. (3) If g j (0) = 1 for 1 j N, then E h (0) = 0. So E h (t ) 0 by (2). However E h (t ) 0 ∀t 0. Hence E h (t ) = 0, implying that g j must be 1 for 1 j N. 2 3.2. Fully discrete scheme To preserve both positivity and entropy property at fully discrete level, we start with a simple time discretization:
mj
g nj +1 − g nj
t
= D − m j + 1 D + g nj +1 + λm2j g nj 1 − g nj +1 for 1 j N ,
(3.3)
2
subject to initial data
g 0j = g 0 (x j ),
1 j N,
and the zero flux boundary condition
m 1 D + g 0n+1 = m N + 1 D + g nN+1 = 0. 2
2
The scheme enjoys following properties. Theorem 3.2. The fully discrete scheme (3.3) is unconditionally stable in the sense that: (1) If g 0j 0 ∀ j, then for any n 1,
0 g nj K ,
1 j N,
where K = max{1, max1 j N g 0j }. (2) The fully discrete entropy
En =
N
m j g nj log g nj + 1 − g nj x
(3.4)
j =1
is nonincreasing, i.e., E n+1 E n ∀n 0. (3) It preserves the equilibrium solution g˜ ≡ 1; i.e. if g 0j = 1, then g nj = 1 ∀n 1. Proof. (1) Rearranging (3.3), we have
+1 +1 −tm j − 1 g nj − + t (m j − 1 + m j + 1 ) + (x)2m j 1 + λtm j g nj g nj +1 − tm j + 1 g nj + 1 1 2
2
2
2
= (x)2m j (1 + t λm j ) g nj .
(3.5)
This is a tridiagonal linear system with the coefficient matrix being diagonally dominant. So { g nj +1 } is nonnegative, as long as { g nj } is nonnegative. It also implies the existence and uniqueness of the numerical solution.
As to the upper bound, it suffices to show that if g nj K ∀ j, then g nj +1 K ∀ j. Let g ni +1 = max1 j N g nj +1 . Then the i-th row of the linear system gives
(x)2mi (1 + t λmi ) g ni (x)2mi 1 + t λmi g ni g in+1 . Therefore
g in+1 − g ni t λmi g ni 1 − g in+1 . This implies that g ni +1 max{1, g ni }. By the definition of K , we always have
g nj +1 K ,
1 j N.
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
663
(2) A simple calculation gives
E n +1 − E n =
N
j =1
=
N
m j g nj +1 log g nj +1 − g nj log g nj − g nj +1 − g nj x
mj
g nj +1 − g nj log g nj +1 + g nj log
j =1
g nj +1 g nj
− g nj +1 − g nj x.
Because m j > 0, g nj 0, and log X X − 1 for X > 0, we have
E
n +1
n
−E
N
mj
g nj +1
g nj
−
log g nj +1
+
g nj
j =1
=
N
g nj +1 g nj
−1 −
g nj +1
−
g nj
x
m j g nj +1 − g nj log g nj +1 x
j =1
= −t
N −1 j =1
m j+ 1
2
+1 g nj + 1
−
g nj +1
log
+1 g nj + 1
g nj +1
− t xλ
N
m2j g nj g nj +1 − 1 log g nj +1
j =1
0. Here we have used the facts that m 1 D + g 0n+1 = m N + 1 D + g nN+1 = 0 and ( X − Y ) log g nj
(3) Suppose N
2
2
X Y
0 for X , Y > 0.
≡ 1. We sum (3.5) over j and end up with
m j (1 + t λm j ) g nj +1 − 1 = 0.
j =1
From (1) we know that 0 g nj +1 1. Therefore g nj +1 = 1 for all j’s. So the scheme preserves the equilibrium solution.
2
Remark 3.1. The scheme (3.3) is easy to implement, yet it is only first order in time. One may easily construct a second order semi-implicit scheme, for example: find g nj +1 = 2g ∗j − g nj by solving
2m j
g ∗j − g nj
t
2 = D − m j + 1 D + g ∗j + λm2j g ∗j 1 − 2g nj + g nj 2
(3.6)
for 1 j N. However, this scheme is no longer unconditionally entropy stable. In practice, one might still prefer to use this scheme due to its second order accuracy in time. 4. Energy stable method for one-species system 4.1. Semidiscrete scheme Let g n (x) be an approximation of g (x, tn ) and G n = G (x, g n ). The initial data is given by g 0 (x) = g 0 (x). We first discretize (2.5) in time to obtain
M
g n +1 − g n
t
G n +1 − G n = ∇ · M ∇ g ∗ + n +1 g − gn
with
g∗ =
g n +1 + g n 2
.
Notice that
G n +1 − G n g n +1
−
gn
λ 2 = λM g ∗ m − M g ∗ − M 2 g ∗ − gn . 3
So (4.1) can be rewritten as a prediction–correction scheme:
(4.1)
Author's personal copy 664
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
(1) Given g n , compute the intermediate solution g ∗ which solves
⎧ ⎨ ⎩
2M
g ∗ − gn
λ
+ P g∗ +
t M g ν∗ = 0
3
M 2 g ∗ − gn
2
= 0 in Ω,
(4.2)
on ∂Ω.
(2) Obtain g n+1 from g n and g ∗ :
g n+1 = 2g ∗ − g n .
(4.3)
Theorem 4.1. Given g n 0 in Ω , if
t
0, hnj 0, 1 j N, then
g nj +1 > 0
and hnj +1 0
1 j N ∀n 0,
provided t < m1 . ∞ (2) The fully discrete entropy
E n = x
N
m j g nj + W j hnj − m j ln m j g nj
j =1
is nonincreasing unconditionally, i.e., E n+1 E n .
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
669
n+1 } depends on t continuously. If we increase t continuously from Proof. (1) Notice that the numerical solution { g n+1 , h 1 n + 1 zero up to m , no component of g can become negative without passing through zero since g n > 0. Assume for some ∞
t, we have g ni +1 = min1 j g nj +1 = 0. Then Eq. (5.4a) gives
−
g ni
t
0,
which contradicts with g ni > 0. So we must have g ni +1 > 0, hence g nj +1 > 0 ∀ j.
n+1 . Assume for some t, we have hn+1 = min1 j hn+1 = −ε , where We apply the similar argument to h i j small number. Then Eq. (5.4b) gives
ε is a sufficiently
hn − mi + mi g in+1 − W i ε ε − i 0, t t 1
1
t
− mi −mi g in+1 + W i ε 0,
since g ni +1 > 0, W i is bounded and
ε is sufficiently small. This contradicts with the choice of t. Therefore, we must have
hni +1 0. (2) For the fully discrete entropy E n , we have
E
n +1
n
−E −
N
n +1
j =1
+1 g nj +1 g nj + 1
α t M j+ 12 ( g j+1 − g j ) x
n 2
− t x
N
F nj +1
2
0.
2
j =1
6. Numerical tests We provide numerical results to demonstrate the accuracy of the schemes and the capacity to capture solution features for large times. Denote the initial data by u 0 (x), the numerical solution at (x j , tn ) by unj = M (x j ) g nj , and the exact solution by u (x j , tn ). Define the L ∞ error as
max unj − u (x j , tn ),
1 j N
and the L 1 error as N n u − u (x j , tn )h. j
j =1
When the exact solution is not available, we replace u (x j , tn ) by a reference solution to compute the errors. 6.1. Tests when P (m) = − ln m Consider the problem with Ω = [0, 1] and
m(x) = 10x2 (x − 1)2 (x − 0.4005)2 . Note that m(0) = m(1) = 0, so Eq. (2.3) becomes degenerate on ∂Ω . We first illustrate the accuracy of scheme (3.3) with the continuous initial data:
u 0 (x) =
sin(100x) + 1 2
(6.1)
.
Table 1 shows that the scheme is of second order in space. Here the reference solution is taken as the numerical solution with N = 1280. The entropy decreasing property is shown in Table 2, where the relative entropy (3.4) is decreasing to 0 as t increases, indicating the numerical convergence to the equilibrium solution m(x) as time becomes large. We also test the numerical convergence to the equilibrium solution m(x) by two numerical examples shown in Fig. 1. The initial data are respectively given by random data and the δ -like function
u 0 (x) =
10[cos(20π (x − 0.8)) + 1]
if |x − 0.8| 0.05;
0
elsewhere.
Author's personal copy 670
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
Table 1 Error and order of accuracy of the numerical solution of (3.3) with the initial data (6.1) on a uniform mesh of N cells, the final time t = 1.0. u 0 (x)
sin(100x)+1 2
N
L ∞ error
Order
L 1 error
Order
20 40 80 160
2.352E−03 2.601E−04 5.763E−05 1.524E−05
– 3.176 2.174 1.919
1.175E−03 1.305E−04 2.821E−05 6.800E−06
– 3.170 2.210 2.053
Table 2 The entropy of the numerical solutions on a uniform mesh with N = 80. t
0
5
10
50
100
200
400
E (t )
1.770E+00
1.534E−01
5.937E−02
2.702E−03
3.719E−04
1.905E−05
9.530E−08
Fig. 1. For the figures in the first row, u 0 (x) is given by the random data; for the second row, u 0 (x) is given by the δ -like function. N = 80.
6.2. Tests when P (m) = −αm Here
α 0 is a parameter. We consider the growth rate
m(x) = 10x2 (x − 1)2 (x − 0.4005)2 − 0.02, which is positive for some x ∈ Ω , and set up
m(x) dx ≈ −0.0048 < 0,
Ω
and
α = 40 and λ = 0.1 such that m(x)e − P dx ≈ 0.0045 > 0.
Ω
In such a setting we have λ∗ = 0: the solution will converge to the positive equilibrium other than m(x), regardless of the size of λ as well as the initial density. We first illustrate the accuracy of scheme (4.6)–(4.7) with Newton’s iteration. The initial data is taken as (6.1), again using the numerical solution with N = 1280 as the reference solution, the results in Table 3 show that the scheme is of the second order in space. We display in Table 4 the discrete energy V at different times, using the same initial data (6.1). Fig. 2 shows that solutions with different initial data converge to the same positive equilibrium solution, as predicted by the theory. Moreover, we observe that the equilibrium solution concentrates at the local positive maximum of m(x). This is consistent with the theoretical result obtained in [5, Theorem 4].
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
671
Table 3 Error and order of accuracy of the numerical solution of scheme (4.6)–(4.7) with initial data (6.1) on a uniform mesh of N cells, the final time t = 1.0. u 0 (x)
sin(100x)+1 2
N
L ∞ error
Order
L 1 error
Order
10 20 40 80 160
4.122E−03 2.010E−03 2.024E−04 4.354E−05 1.040E−05
– 1.036 3.312 2.217 2.066
6.491E−02 1.762E−02 4.156E−03 1.033E−03 2.552E−04
– 1.881 2.084 2.008 2.018
Table 4 The energy of the numerical solutions on a uniform mesh of N = 80 cells. tn
0
5
10
50
100
200
Vn
7.387E+02
4.405E−04
7.523E−05
6.435E−07
4.125E−08
−1.522E−08
Fig. 2. For the figures in the first row, u 0 (x) is given by the random data; for the second row, u 0 (x) is given by the δ -like function. N = 80.
6.3. Two-dimensional tests with P (m) = −αm2 We consider the resource
2 2 m(x) = 10 0.25 − (x1 − 0.5)2 − (x2 − 0.5)2 e (x1 −0.5) −(x2 −0.5) − 0.3 (see Fig. 3), and set
α = 4 so that
m(x) dx ≈ −0.1198 < 0 and Ω
m(x)e − P dx ≈ −0.1590 < 0.
Ω
In such a setting we have λ∗ > 0: the solution will converge to the positive equilibrium if λ > λ∗ , and decay to zero if λ λ∗ , regardless of the choice of initial density. We test scheme (4.6)–(4.7) with two initial data: (i) the cos-like function u 0 (x) = cos(3π |x|2 ) + 1; (ii) the δ -like function concentrating at (x0i , x0 j ) for i , j = 1, 2 with x01 = 0.2 and x02 = 0.8.
u 0 (x) =
25[cos(10π (x1 − x0i )) + 1][cos(10π (x2 − x0 j )) + 1]
if |xi − x0 j | 0.1 and i , j = 1, 2;
0
elsewhere.
Author's personal copy 672
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
2 2 Fig. 3. Plot of m(x) = 10[0.25 − (x1 − 0.5)2 − (x2 − 0.5)2 ]e (x1 −0.5) −(x2 −0.5) − 0.3.
Fig. 4. Figures in the first row, u 0 (x) is given by the cos-like function; in the second row, u 0 (x) is given by the δ -like function. N = 20 × 20, λ = 1 < λ∗ .
From Fig. 4, we observe that when t = 100, the numerical solutions is of order O(10−7 ), approaching zero. In other words, the species will extinct eventually, consistent with the theoretical result stated in Theorem 1.1. Approach toward the positive equilibrium is shown with two examples in Fig. 5, corresponding to the persistence of the species as predicted by Theorem 1.1. Moreover, the energy V ≈ −0.13439 when t = 0.5 and 100 for both of the initial data. These tests indicate that the steady state is asymptotically approached in time and well preserved. 6.4. Two-dimensional test for the two-species system We consider the two-species system with resource potentials P (m) = − ln m and Q (m) = −m, where the resource m(x) is given by
m(x) = 0.5 sin(π x1 ) sin(π x2 ) + 0.6, (See Fig. 6.) We set
x ∈ Ω = [0, 2.5]2 .
α = 5, β = 1, and take initial value as
u 0 (x) = 1 − 2v 0 (x),
v 0 (x) =
1 4
sin2 (4π x1 ) sin2 (4π x2 ) if 1 x1 , x2 1.5,
0
elsewhere.
Visually, we observe the approach toward the equilibrium solution (m, 0) through the numerical result shown in Fig. 7. Quantitatively, Table 5 shows that the entropy is decreasing as time evolves and indicates the asymptotic approach toward the equilibrium state, while we note that E ≈ 5.34191 when the solution approaches the steady state. 7. Concluding remarks We have developed finite difference schemes for a class of reaction–diffusion–advection equations arising in the evolution of biased dispersal of population dynamics, including the single species model and two competing species system. The schemes are shown to preserve both positivity and equilibrium, satisfying either entropy or energy dissipation structure.
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
673
Fig. 5. Figures in the first row, u 0 (x) is given by the cos-like function; in the second row, u 0 (x) is given by the δ -like function. N = 20 × 20, λ = 1000 > λ∗ .
Fig. 6. Plot of m(x) = 0.5 sin(π x1 ) sin(π x2 ) + 0.6. Table 5 The discrete entropy on a uniform mesh with N = 40 × 40. tn
0
0.1
0.3
0.5
100
En
6.26436
5.80962
5.66577
5.57095
5.35172
Fig. 7. The first row: u component; the second row: v component. N = 40 × 40.
Author's personal copy 674
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
These ensure that the numerical solution provides a satisfying long-time behavior, thus revealing the desired dispersal pattern. Numerical examples are shown to illustrate the second-order accuracy for various resource potentials and underline the efficiency to preserve the large-time asymptotic. We also proved that the unique existence of the numerical solutions for the proposed schemes. Our future work would develop higher order methods based on the proposed semidiscrete formulations. For dispersal models with multi-species, the dynamics become much more complex, further numerical investigation into multi-species dispersal models deserves serious consideration. Acknowledgements This research was partially supported by the National Science Foundation under grant DMS-1312636 and the Mathematical Research Network KI-net grant. Appendix A Theorem A.1. Problem (2.10) has a unique λ∗ > 0 characterized by a positive eigenfunction if and only if
mM dx < 0. Ω
Moreover, λ∗ = 0 if and only if Ω mM dx 0. Proof. The existence of (λ∗ , ψ) with ψ > 0 in Ω and M = 1 is given in [6]. In a similar manner, one can prove the existence of (λ∗ , ψ) with ψ > 0 and nonconstant M > 0. Also, λ∗ can be expressed by the variational formulation
λ∗ = inf
2 Ω mM v dx
v∈S
where
∇ v 2M
,
S = v ∈ H 1 , M v ν = 0 on ∂Ω and
mM v 2 dx > 0 .
Ω
In fact, (2.10) yields
∇ψ 2M 0. 2 Ω mM ψ dx
λ∗ =
(A.1)
Divide both sides of (2.10) by ψ , and integrate over Ω to obtain
λ∗
mM dx = −
Ω
M |∇ψ|2
ψ2
dx 0.
(A.2)
Ω
For Ω mM dx < 0, (A.2) implies that λ∗ > 0. For Ω mM dx = 0, (A.2) implies that ∇ψ 2M = 0. Therefore λ∗ = 0 by (A.1). For Ω mM dx > 0, we must have λ∗ = 0 by (A.2). 2 Let (σ ∗ , φ) denote the pair of the principle eigenvalue and eigenfunction of (2.9), we thus have φ > 0 in Ω , see [38, Theorem 11.10]. The following result relates the sign of σ ∗ to the relative size of λ in terms of λ∗ . Theorem A.2. sign(σ ∗ ) = sign(λ − λ∗ ). Proof. Rewrite (2.9) as
∇ · ( M ∇φ) + λ m −
σ∗ λ
M φ = 0,
∗ = m − σλ . Then λ is determined by Let m
∇ v 2M , 2 Ω m M v dx
λ = inf v ∈ S˜
M φν |∂Ω = 0.
(A.3)
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
˜ If where S˜ is a modification of S with m replaced by m. < m and eigenvalue. If σ ∗ > 0, then m
675
σ ∗ = 0, then φ = ψ and λ = λ∗ by the simpleness of the principle
M φ dx < m 2
0< Ω
mM φ 2 dx. Ω
Therefore
∇φ 2M ∇φ 2M > λ∗ . 2 dx 2 dx m M φ mM φ Ω Ω
λ= If
> m and σ ∗ < 0, then m 2 M ψ dx > mM ψ 2 dx > 0. m Ω
Ω
Hence,
∇ψ 2M ∇ψ 2M < = λ∗ . 2 dx 2 dx m M ψ mM ψ Ω Ω
λ
2
The result about the equilibrium solution g˜ is the following.
¯ is positive somewhere in Ω and P is bounded. Theorem A.3. Suppose that m ∈ C 2 (Ω) (1) If 0 < λ λ∗ , (2.7) only has a zero solution. (2) If λ > λ∗ , (2.7) has a zero solution and a unique positive solution. Proof. (1) For 0 < λ λ∗ , multiply (2.7) by g˜ to get
mM g˜ 2 dx =
λ Ω
M 2 g˜ 2 dx 0.
M |∇ g˜ |2 dx + λ Ω
Ω
By the definition of λ∗ , we have
λ
mM g˜ 2 dx λ∗
Ω
Ω
M g˜ dx λ − λ 2 2
λ Ω
M 2 g˜ 2 dx. Ω
It follows that
mM g˜ 2 dx + λ
∗
mM g˜ 2 dx 0. Ω
Therefore g˜ ≡ 0 on Ω . m(x) (2) One can verify that any constant g¯ greater than M (x) ∞ is a super-solution of (2.7). We normalize the positive eigenfunction φ such that φ ∞ = 1. Then g = ε φ with
0 < ε min
σ∗ λ M ∞
, g¯
is a sub-solution of (2.7) because
P g = ε M φ ε λ M φ − σ ∗ ε M φ ε λ M − σ ∗ 0 in Ω Moreover,
0 < g g¯ and
and
M g |∂Ω = 0. ν
∂ λ F ( g , m) = −λ M 2 < 0. ∂g g
So the conditions of [37, Theorem 3.4 in Chapter 3] are satisfied, and it ensures that there exists a unique solution g˜ of (2.7) such that g g˜ g¯ . The above claimed uniqueness is actually restricted to the situation that any two solutions, if g˜ 1 g˜ 2 , must coincide. For the general case, we define the minimum of two arbitrary positive solutions as
Author's personal copy 676
H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
g ∗ (x) = min g˜ 1 (x), g˜ 2 (x) . A slight modification of the argument in [21, Lemma 1.5] when applied to (2.7) (care is needed for treating the Neumann boundary condition) enables us to show that g ∗ is a weak super-solution. Subtracting the weak formulations of g ∗ and g˜ i against the test function taken as g˜ i and g ∗ respectively, we have
λ
M 2 g ∗ g˜ i g˜ i − g ∗ dx 0 for i = 1, 2,
Ω
which implies that g ∗ = g˜ i in Ω , hence g˜ 1 = g˜ 2 .
2
The above results when combined with the nonincreasing property of the energy V defined in (2.4) enable us to predict persistence or extinction of the evolving species; see [7] for more background details. Theorem A.4. Under the same assumptions as in Theorem A.3, we have: (1) If 0 < λ λ∗ , all nonnegative solutions of (2.1) decay toward zero as t → ∞. (2) If λ > λ∗ , the positive equilibrium is globally attractive among nonzero nonnegative solutions. Proof. We first prove the convergence to the steady state. The discussion on V in Section 2.2 shows that
d dt
M (∂t g )2 dx 0.
V =−
(A.4)
Ω
Thus V is nonincreasing along the orbit starting at g 0 (x). The boundedness and precompactness of the orbits imply that V [ g ] is bounded from below, so there exists V eq = limt →∞ V [ g ]. If w is in the ω -limit set of g 0 (x), then V [ w ] = V eq . Let w ∗ (x, t ) be the orbit starting at some point w in the ω -limit set of g 0 (x). Since that set is invariant, the entire orbit w ∗ (x, t ) must belong to it, so V [ w ∗ ] = V eq . In view of (A.4), this is possible if and only if ∂t w ∗ = 0. Therefore, w must be an equilibrium of (2.1). Since w is an arbitrary element of the ω -limit set of an arbitrary initial point g 0 , the ω -limit set for the semi-flow generated by (2.1) must consist of all the equilibria, which shows the global attractiveness of stable equilibrium solutions. We now distinguish two cases in terms of λ. The claimed results follow from the fact that zero equilibrium solution is stable for 0 < λ λ∗ and unstable for λ > λ∗ , as inferred from Theorems A.2 and A.3. 2 References [1] E. Celledoni, V. Grimm, R.I. McLachlan, D.I. McLaren, D. O’Neale, B. Owren, G.R.W. Quispel, Preserving energy resp. dissipation in numerical PDEs using the “Average Vector Field” method, arXiv:1202.4555, 2012. [2] Fethi Belgacem, Chris Cosner, The effects of dispersal along environmental gradients on the dynamics of populations in heterogeneous environments, Can. Appl. Math. Q. 3 (4) (1995) 379–397. [3] R. Bermejo, J. Carpio, An adaptive finite element semi-Lagrangian implicit–explicit Runge–Kutta–Chebyshev method for convection dominated reaction– diffusion problems, Appl. Numer. Math. 58 (1) (January 2008) 16–39. [4] Rodolfo Bermejo, Mofdi El Amrani, A finite element semi-Lagrangian explicit Runge–Kutta–Chebyshev method for convection dominated reaction– diffusion problems, J. Comput. Appl. Math. 154 (1) (May 2003) 27–61. [5] Andriy Bezuglyy, Yuan Lou, Reaction–diffusion models with large advection coefficients, Appl. Anal. 89 (7) (2010) 983–1004. [6] K.J. Brown, S.S. Lin, On the existence of positive eigenfunctions for an eigenvalue problem with indefinite weight function, J. Math. Anal. Appl. 75 (1) (1980) 112–120. [7] Robert Stephen Cantrell, Chris Cosner, Spatial Ecology via Reaction–Diffusion Equations, John Wiley & Sons, Ltd., 2004. [8] Robert Stephen Cantrell, Chris Cosner, Donald L. Deangelis, Victor Padron, The ideal free distribution as an evolutionarily stable strategy, J. Biol. Dyn. 1 (3) (2007) 249–271. [9] Robert Stephen Cantrell, Chris Cosner, Yuan Lou, Movement toward better environments and the evolution of rapid diffusion, Math. Biosci. 204 (2) (2006) 199–214. [10] Robert Stephen Cantrell, Chris Cosner, Yuan Lou, Approximating the ideal free distribution via reaction–diffusion–advection equations, J. Differ. Equ. 245 (12) (2008) 3687–3703. [11] Robert Stephen Cantrell, Chris Cosner, Yuan Lou, Evolution of dispersal and the ideal free distribution, Math. Biosci. Eng. 7 (1) (2010) 17–36. [12] J.A. Carrillo, M.P. Gualdani, A. Jüngel, Convergence of an entropic semi-discretization for nonlinear Fokker–Planck equations in Rd , Publ. Mat. 52 (2) (2008) 413–433. [13] Xinfu Chen, Richar Hambrock, Yuan Lou, Evolution of conditional dispersal: a reaction–diffusion–advection model, J. Math. Biol. 57 (3) (2008) 361–386. [14] Q. Chang, E. Jia, W. Sun, Difference schemes for solving the generalized nonlinear Schrödinger equation, J. Comput. Phys. 148 (1999) 397–415. [15] Chris Cosner, Yuan Lou, Does movement toward better environments always benefit a population?, J. Math. Anal. Appl. 277 (2) (2003) 489–503. [16] Jérôme Coville, Juan Dávila, Salomé Martínez, Nonlocal anisotropic dispersal with monostable nonlinearity, J. Differ. Equ. 244 (12) (2008) 3080–3118. [17] Feng Chen, Jie Shen, Efficient energy stable schemes with spectral discretization in space for anisotropic Cahn–Hilliard systems, Commun. Comput. Phys. 13 (2013) 1189–1208. [18] Craig Collins, Jie Shen, Steven M. Wise, An efficient, energy stable scheme for the Cahn–Hilliard–Brinkman system, Commun. Comput. Phys. 13 (2013) 929–957. [19] M. Dahlby, B. Owren, A general framework for deriving integral-preserving numerical methods for PDEs, SIAM J. Sci. Comput. 33 (2011) 2318–2340.
Author's personal copy H. Liu, H. Yu / Journal of Computational Physics 256 (2014) 656–677
677
[20] David J. Eyre, Unconditionally gradient stable time marching the Cahn–Hilliard equation, in: Computational and Mathematical Models of Micro Structural Evolution, San Francisco, CA, 1998, in: Mater. Res. Soc. Symp. Proc., vol. 529, MRS, Warrendale, PA, 1998, pp. 39–46. [21] E.N. Dancer, G. Sweers, On the existence of a maximal weak solution for a semilinear elliptic equation, Differ. Integral Equ. 2 (4) (1989) 533–540. [22] Jack Dockery, Vivian Hutson, Konstantin Mischaikow, Mark Pernarowski, The evolution of slow dispersal rates: a reaction diffusion model, J. Math. Biol. 37 (1) (1998) 61–83. [23] Daisuke Furihata, Finite difference schemes for ∂ u /∂ t = (∂/∂ x)α δ G /δ u that inherit energy conservation or dissipation property, J. Comput. Phys. 156 (1) (1999) 181–205. [24] Richard Gejji, Yuan Lou, Daniel Munther, Justin Peyton, Evolutionary convergence to ideal free dispersal strategies and coexistence, Bull. Math. Biol. 74 (2) (2011) 1–43. [25] W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion–Reaction Equations, first edition, Springer Ser. Comput. Math., vol. 33, Springer, 2003. [26] Chiu-Yen Kao, Yuan Lou, Wenxian Shen, Random dispersal vs. non-local dispersal, Discrete Contin. Dyn. Syst. 26 (2) (2010) 551–596. [27] King-Yeung Lam, Concentration phenomena of a semilinear elliptic equation with large advection in an ecological model, J. Differ. Equ. 250 (1) (2011) 161–181. [28] King-Yeung Lam, Wei-Ming Ni, Limiting profiles of semilinear elliptic equations with large advection in population dynamics, Discrete Contin. Dyn. Syst. 28 (3) (2010) 1051–1067. [29] D. Lanser, J.G. Verwer, Analysis of operator splitting for advection–diffusion–reaction problems from air pollution modeling, J. Comput. Appl. Math. 111 (1–2) (1999) 201–216. [30] H. Liu, H. Yu, An entropy satisfying conservative method for the Fokker–Planck equation of the finitely extensible nonlinear elastic dumbbell model, SIAM J. Numer. Anal. 50 (3) (2012) 1207–1239. [31] Christian Lubich, Alexander Ostermann, Runge–Kutta approximation of quasi-linear parabolic equations, Math. Comput. 64 (210) (1995) 601–627. [32] S. Li, L. Vu-Quoc, Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein–Gordon equation, SIAM J. Numer. Anal. 32 (1995) 1839–1875. [33] Takayasu Matsuo, Dissipative/conservative Galerkin method using discrete partial derivatives for nonlinear evolution equations, J. Comput. Appl. Math. 218 (2) (2008) 506–521. [34] Takayasu Matsuo, Hisashi Yamaguchi, An energy-conserving Galerkin scheme for a class of nonlinear dispersive equations, J. Comput. Phys. 228 (12) (2009) 4346–4358. [35] J.D. Murray, Spatial models and biomedical applications, in: Mathematical Biology. II, third edition, in: Interdiscip. Appl. Math., vol. 18, Springer-Verlag, New York, 2003. [36] Akira Okubo, Diffusion and Ecological Problems: Mathematical Models, Biomathematics, vol. 10, Springer-Verlag, Berlin, 1980, an extended version of the Japanese edition. Ecology and Diffusion, translated by G.N. Parker. [37] Chia-Ven Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, 1992. [38] Joel Smoller, Shock Waves and Reaction–Diffusion Equations, Grundlehren Math. Wiss. (Fundam. Princ. Math. Sci.), vol. 258, Springer-Verlag, New York, 1983.
