close

Вход

Забыли?

вход по аккаунту

?

j.jfranklin.2018.07.029

код для вставкиСкачать
Accepted Manuscript
Interaction between prey and mutually interfering predator in prey
reserve habitat: pattern formation and the Turing-Hopf bifurcation
Jai Prakash Tripathi, Syed Abbas, Gui-Quan Sun, Debaldev Jana,
Cui-Hua Wang
PII:
DOI:
Reference:
S0016-0032(18)30507-6
https://doi.org/10.1016/j.jfranklin.2018.07.029
FI 3576
To appear in:
Journal of the Franklin Institute
Received date:
Revised date:
Accepted date:
4 January 2018
6 June 2018
31 July 2018
Please cite this article as: Jai Prakash Tripathi, Syed Abbas, Gui-Quan Sun, Debaldev Jana,
Cui-Hua Wang, Interaction between prey and mutually interfering predator in prey reserve habitat: pattern formation and the Turing-Hopf bifurcation, Journal of the Franklin Institute (2018), doi:
https://doi.org/10.1016/j.jfranklin.2018.07.029
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service
to our customers we are providing this early version of the manuscript. The manuscript will undergo
copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please
note that during the production process errors may be discovered which could affect the content, and
all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Interaction between prey and mutually interfering predator in
prey reserve habitat: pattern formation and the Turing–Hopf
bifurcation
Jai Prakash Tripathi1 , Syed Abbas2 , Gui-Quan Sun3,† , Debaldev Jana4 ,
Cui-Hua Wang3
1 Department
of Mathematics, Central University of Rajasthan,
NH-8, Bandarsindri, Kishangarh-305817, Distt.-Ajmer, Rajasthan, India
2
School of Basic Sciences, Indian Institute of Technology Mandi,
Mandi-175001, Himachal Pradesh, India
Systems Research Center, Shanxi University
CR
IP
T
3 Complex
Taiyuan, Shanxi 030006, People’s Republic of China
4 Department
of Mathematics & SRM Research Institute,
SRM Institute of Science and Technology, Kattankulathur-603 203, Tamil Nadu, India
†
Corresponding Author. Email: sunguiquan@sxu.edu.cn
Abstract
ED
M
AN
US
In this paper, we propose a diffusive prey-predator system with mutually interfering predator (Crowley-Martin functional response) and prey reserve. In particular,
we develop and analyze both spatially homogeneous model based on ordinary differential equations and reaction-diffusion model. We mainly investigate the global
existence and boundedness of positive solution, stability properties of homogeneous
steady state, non-existence of non-constant positive steady state, conditions for
Turing instability and Hopf bifurcation of the diffusive system analytically. Conventional stability properties of the non-spatial counterpart of the system are also
studied. The analysis ensures that the prey reserve leaves stabilizing effect on the
stability of temporal system. The prey reserve and diffusive parameters leave parallel impact on the stability of the spatio-temporal system. Furthermore, we illustrate
the spatial patterns via numerical simulations, which show that the model dynamics
exhibits diffusion controlled pattern formation by different interesting patterns.
Introduction
CE
1
PT
Keywords: Crowley-Martin functional response; Prey reserve; Hopf bifurcation;
Global stability; Positive steady states; Turing pattern
AC
The interactions of individual organisms with each other and the environment are the
most important distinctive features in ecological systems. Predation is an important type
of interaction that affects the population dynamics of all species. As a result, one of the
main objectives of ecologists is to get insight into prey-predator relationship [1, 2, 3, 4].
The rate of prey consumption by an average predator (functional response) is one of
the significant characteristics of the prey-predator relationship. Functional response determines stability and bifurcation scenarios with diverse parameter range of ecological
systems. The real functional responses not only show the consistency of data fitting with
experimental results but also play an important role for accurate predictions of models. Functional responses can primarily be classified in two categories: (i) prey dependent
(linear and Holling type functional responses) [5, 6], (ii) predator dependent (ratio dependent, Beddington-DeAngelis, Crowley-Martin and Hassel-Varley) [7]–[13]. The intuitive
1
ACCEPTED MANUSCRIPT
= Y
CE
dY
dT
X
C(1 − m)Y
= X R(1 − ) −
,
K
A1 + B1 (1 − m)X + C1 Y + B1 C1 (1 − m)XY
PT
dX
dT
ED
M
AN
US
CR
IP
T
and experimental observations [14]–[16] indicate that the feeding rate of per unit consumer
decreases because of mutual interference among predators. The assumption that predator’s predation decreases due to high predator density (interference among the predator
individuals) even when prey density is high (presence of prey handling or searching of prey
by predator individual) was modelled mathematically by Crowley and Martin [17] (hereafter Crowley-Martin functional response and the CM model). For the detailed study
of Crowley-Martin functional response and its formulation, one can refer [18]–[22] and
references cited therein.
The experimental study on prey (Paramecium caudatum) and predator (Didinium
nasatum) shows that in a completely homogeneous environment, Didinium destructs all
prey and it collapses subsequently. Moreover, due to many reasons (such as, environmental pollution, over predation, unregulated harvesting, etc.) several species have became
extinct and many others are at the verge of extinction. This also signifies the persistence
(the survival of species for a longer time) of the system [23, 24]. Improving conditions
or adaptations [25] of natural habitat, imposing restriction on species harvesting, creating natural reserves, establishment of protected areas are the suitable environmental
strategies to save the species from getting extinct. More specifically, some fractions of the
prey population can partially be protected by using prey refuges (reserves). Environmental heterogeneity also enhances spatial refuges of prey. In a broader sense, any strategy
that decreases predation is refuge. Mite predator-prey interactions often exhibit spatial
refugia, which afford the prey some degree of protection from predation and reduce the
chance of extinction due to predation [26]. The term reserve zone/refuge has received
considerable attention in the dynamics of prey-predator model ([27]–[33]). In particular,
spatial prey refuge is one of the most relevant behavioural traits that affects the dynamics
of the prey-predator systems [34].
Motivated by the above discussion, in the proposed system, we introduce both predator
dependent functional response (Crowley-Martin) as well as prey cover (reserve). We
investigate the following prey-predator model:
−D+
(1)
EC(1 − m)X
,
A1 + B1 (1 − m)X + C1 Y + B1 C1 (1 − m)XY
AC
where X(t) and Y (t) represent the prey and predator populations at time t, respectively,
with X(0) > 0 and Y (0) > 0. Here, it is assumed that in the absence of predator population, prey population follows logistic growth rate while predators are subjected to natural
mortality. The system (1) is incorporated by a refuge which protects mX of prey, where
m ∈ [0, 1). Thus only (1 − m)X of the prey is available to predator. The parameters
C, E, B1 , C1 are also termed as saturating Crowley-Martin type response function parameters. The parameters E and C stand for the effect for capture rate and conversion factor
denoting the newly born predators for each captured prey and effect of handling time for
predators, respectively. The parameter C1 measures the magnitude of interference among
predators. The biological meaning of all other parameters can be found in [35]. All the
parameters, R, K, C, D, E, A1 , B1 , C1 in the system (1) assume only positive values and
will be considered as constants throughout our discussion. After non-dimensionalization
of the system (1) with the following set of variables and parameters:
2
ACCEPTED MANUSCRIPT
C
B1 K
C1
D
CEK
B1 C1 K
X
, y = Y, a =
,b =
,c =
,d = ,e =
and f =
,
K
RA1
A1
A1
R
RA1
A1
the system (1) takes the following form:
dx
a(1 − m)y
= x (1 − x) −
= xf1 (x, y),
dt
1 + b(1 − m)x + cy + f (1 − m)xy
(2)
dy
e(1 − m)x
= y −d+
= yf2 (x, y).
dt
1 + b(1 − m)x + cy + f (1 − m)xy
CR
IP
T
t = RT, x =
PT
ED
M
AN
US
Living status of population in their natural environment are very rarely spatially homogeneous [36]. Often this heterogeneity takes the form of a patchiness or other type of
structural complexity when areas with high population density alternate with areas where
the given species is virtually absent depending upon a particular environmental condition
[37]. For instance, an insect population living in corn fields in considerable densities may
be almost absent in the space between the fields simply because of unfavourable conditions
there. In spatially inhomogeneous case, the existence of a non-constant time-independent
positive solution (stationary pattern) ensures the dynamical richness of the system. Moreover, local interactions in space can give rise to large scale spatio-temporal patterns (e.g.,
spiral waves, spatio-temporal chaos (turbulence), stationary (Turing-type) patterns and
transitions between these modes). These phenomena are found in many organizational
levels of biotic systems, and their occurrence and properties are largely independent of
the precise interaction structure. Such type of spatio-temporal patterns have significant
consequences for fundamental biological processes. Pattern formation [38]–[43] study in
reaction-diffusion systems [44, 45] has been an informative and active research area since
the pioneering works of Alan Turing [46] in 1952. Turing patterns of various diffusive
prey-predator models have been studied extensively from several researchers [47, 48, 49].
For detailed study of pattern of diffusive prey-predator systems, one can refer [46, 50] and
references cited therein. Taking care of random movement of the individuals in two dimensional habitat, we incorporate diffusion terms into the growth equations for prey and
predator populations given by (2), to obtain the following reaction diffusion equations:
AC
CE
a(1 − m)y
∂u
= u (1 − u) −
+ d1 ∇2 u,
∂t
1 + b(1 − m)u + cv + f (1 − m)uv
(3)
∂v
e(1 − m)u
= v −d+
+ d2 ∇2 v,
∂t
1 + b(1 − m)u + cv + f (1 − m)uv
for (X, t) ∈ Ω × R, where Ω ⊂ R2 , and subjected to initial conditions
u(X, 0) = u0 (X) ≥ 0, v(X, 0) = v0 (X) ≥ 0, X ∈ Ω
(4)
and no-flux boundary condition
∂v
∂u
=
= 0, X ∈ ∂Ω, t > 0.
∂n
∂n
(5)
Here, we should note that the boundary of the region Ω is in C ∞ i.e., ∂Ω ∈ C ∞ The
boundary is smooth enough and there is no break or sharp edge in the boundary. In Eq.
3
ACCEPTED MANUSCRIPT
CR
IP
T
(3), the positive parameters d1 and d2 are the diffusion coefficients of prey and predator,
∂2
∂2
respectively, ∇2 =
+
is the usual two dimensional Laplacian operator in variable
∂x2 ∂y 2
X ≡ (x, y) ∈ Ω ⊂ R2 and Ω is a square bounded (simply connected) domain with smooth
∂
boundary ∂Ω. n is the outward unit normal vector of the boundary ∂Ω while
stands
∂n
for outward drawn normal derivative on the boundary ∂Ω. Thus the prey and predator
density at location ‘X’ and time ‘t’ is denoted by u(X, t) and v(X, t) respectively. The
zero flux boundary conditions (5) assure that there is no fluxes of populations through
the boundary, i.e., no external input is imposed from outside. In the present paper, we
will be concerned with the following two important questions:
i How the prey reserve parameter (m) and diffusion coefficients (d1 , d2 ) affect the
dynamical behaviour (co-existence, local and global stability and pattern structure)
of the systems (2) and 3).
AN
US
ii To explore the similarities/differences between sufficient parametric conditions for
co-existence of steady states, local and global stability of interior equilibrium of
temporal system and its spatio-temporal counterpart.
2.1
Analysis of the system (2)
Positivity and boundedness of the solutions
PT
2
ED
M
The rest of the paper is structured as follows. Permanence, persistence, local and
global stability of the equilibria and related dynamical behaviors of the system (2) are
considered in the next section. The dynamics of the associated spatio-temporal model
(3) is discussed in the Section 3. More precisely, in Section 3, we study the existence and
uniqueness of the global positive solution and nonexistence of nonconstant positive steady
states and Turing pattern formation in the system (3). In Section 4, we perform some
carefully designed numerical simulations to substantiate our analytical findings. Paper
ends with a section on discussion and the ecological significance of our analytical results.
Some future directions are also mentioned.
CE
n
Let R+ denotes the set of all non-negative real numbers and R+
= {x ∈ Rn : x =
(x1 , ..., xn ) where xi ∈ R+ , ∀ i = 1, 2, ...n}. Here we discuss the invariance of the boundary
2
and the interior of R+
along with conservation of overall energy or biomass flow.
AC
Lemma 2.1. The positive quadrant (Int(R2+ )) is invariant for the system (2).
Proof. The proof is straightforward and interested readers can refer to [51].
Lemma 2.2. All solutions of the system (2) that start in R2+ are confined to the region
1
D = {(x, y) ∈ R2+ : x(t) ≤ 1, 0 ≤ η(t) ≤ 1 + } as t → ∞ for all positive initial value
4d
a
2
(x(0), y(0)) ∈ R+ where η(t) = x + y(t).
e
Remark 2.3. The proof of this theorem follows similarly as in [52]. Geometrically, the
Lemma 2.2 assures that all the trajectories of the system (2) initiating from any point
in R2+ ultimately lie in the fixed bounded region defined by D. Hence, the flow/dynamical
system associated with system (2) is dissipative.
4
ACCEPTED MANUSCRIPT
2.2
Permanence
Theorem 2.4. Define M1 := 1, M2 := (e−bd)(1−m), m1 :=
Any positive solution {x(t), y(t)} of system (2)satisfies
m1 ≤ lim inf t→∞ x(t) ≤ lim supt→∞ x(t) ≤ M1 ,
m2 ≤ lim inf t→∞ y(t) ≤ lim supt→∞ y(t) ≤ M2 ,
o
n
d
c
.
whenever m < min 1 − , 1 −
a
(e − bd)m1
c−a(1−m)
, m2
c
:=
m1 (e−bd)(1−m)−d
.
(c+f m1 )d
Remark 2.5. For the detailed proof of the Theorem 2.4, the reader is referred to [53].
CR
IP
T
Using the definitions of lim sup and lim inf, one can claim that there exist four positive
constants ti (i = 1, 2, 3, 4) such that m1 ≤ x(t), ∀t ≥ t1 , x(t) ≤ M1 , ∀t ≥ t2 , m2 ≤ y(t),
∀t ≥ t3 , y(t) ≤ M2 , ∀t ≥ t4 . If we define m∗ = min(m1 , m2 ), M ∗ = max(M1 , M2 ) and
T = max(ti , 1 ≤ i ≤ 4, i ∈ N), then we have the following permanence result for (2).
AN
US
Theorem 2.6. Suppose that the time independent coefficients a, b, c, d, e, f are positive
and satisfy
n
o
c
d
m < min 1 − , 1 −
,
(6)
a
(e − bd)m1
then the model system (2) is permanent.
Theorem 2.7. If e < bd, then limt→∞ y(t) = 0.
2.3
ED
M
Remark 2.8. From the predator equation in (2), the Theorem 2.7 follows immediately.
One can easily observe that when f = m = 0 (absence of prey reserve and prey handling
at high prey density), the above discussion also remains valid. In this case, the system
(2) reduces to the system with Beddington-DeAngelis functional response by Cantrell and
Cosner [54] and Takeuchi et al. [55] with density dependence. The theorem 2.4 and 2.6
reduces to the corresponding results in [55].
Existence and stability around biomass equilibria
PT
The co-existence conditions of all trivial and non-trivial equilibria of the system (2) and
their dynamical behaviors are presented in the following table.
Feasibility
condition
(i)
E 0 (0, 0)
always
(ii)
E 1 (1, 0)
always
AC
CE
Equilibrium and
coordinate
(iii)
E ∗ (x∗ , y ∗ )
Jacobian matrix
and
eigenvalues
1
0
0 −d
λ1 = 1, λ2 = −d
a(1−m)
−1
− 1+b(1−m)
0
a(1 − m) >
(c + f (1 −
m)x∗ )(1
−
x∗ )
Stability
status
saddle
!
e(1−m)
1+b(1−m)
e(1−m)
−d + 1+b(1−m)
−d +
λ1 = −1, λ2 =
j11 j12
,
j21 j22
λ1 + λ2 = j11 + j22
λ1 λ2 = j11 j22 − j12 j21
stable
if m < 1 −
d
e−bd
stable if
j11 < 0
Here, x∗ is given by the root of the following cubic equation in z
αz 3 + βz 2 + γz + δ = 0,
(7)
2
2
where, α = ef (1 − m)
, β = e(1 − m) c − f (1 − m) + df (1 − m), γ = a(1 − m) (e −
bd) − c d + e(1 − m) , δ = −ad(1 − m).
5
ACCEPTED MANUSCRIPT
CR
IP
T
The constant term (δ) of the equation (7) is negative and the leading coefficient (α)
is always positive. Thus, the Descartes rule of sign ensures that the above cubic equation
(7) possesses at least one positive root, say x∗+ . Furthermore, for this value of x∗+ , the
(1 + b(1 − m)x∗+ )(1 − x∗+ )
∗
∗
=
.
will be given by y+
corresponding value of y+
a(1 − m) − (c + f (1 − m)x∗+ )(1 − x∗+ )
a(1 − m)x∗ y ∗ (b(1 − m) + f (1 − m)y ∗ )
,
Moreover, j11 = x∗ f1 x = −x∗ +
(1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2
−a(1 − m)x∗ (1 + b(1 − m)x∗ )
j12 = x∗ f1 y =
< 0,
(1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2
e(1 − m)y ∗ (1 + cy ∗ )
j21 = y ∗ f2 x =
> 0,
(1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2
e(1 − m)x∗ y ∗ (c + f (1 − m)x∗ )
j22 = y ∗ f2 y = −
< 0.
(1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2
The following theorem provides alternative conditions for local asymptotic stability of
E ∗.
AN
US
Theorem 2.9. If permanence condition (Theorem 2.6) and m < 1 +
hold, then the positive equilibrium E ∗ is locally asymptotically stable.
f − a + x∗ (bc − f )
ab(c − x∗ )
Global stability: We now discuss the global stability of interior equilibrium solution
E ∗ = (x∗ , y ∗ ) by constructing a suitable Liapunov functional.
M
Theorem 2.10. Let ξ = (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )(1 + b(1 − m)x +
cy + f (1 − m)xy). If permanence conditions given in Theorems (2.4 and 2.6) and m >
ay ∗
1−
hold then E ∗ = (x∗ , y ∗ ) is globally asymptotically stable.
(1 − x∗ )(b + aM2 )
ED
Proof. The system (2) can be re-written in the following form:
PT
dx
a(1 − m)y ∗
= x − (x − x∗ ) +
dt
1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗
a(1 − m)y
,
−
1 + b(1 − m)x + cy + f (1 − m)xy
(8)
dy
e(1 − m)x
e(1 − m)x
= y −
+
.
dt
1 + b(1 − m)y ∗ + cy ∗ + f (1 − m)x∗ y ∗ 1 + b(1 − m)y + cy + f (1 − m)xy
CE
∗
AC
Now, define V (x, y) : R2+ → R, such that.
V (x, y) = V1 (x) + V2 (y),
(9)
where, V1 (x) = x − x∗ − x∗ ln(x/x∗ ) , V2 (y) = η y − y ∗ − y ∗ ln(y/y ∗ ) . The constant η
is a positive constant and is defined in the discussion below. The time derivatives of V1
6
ACCEPTED MANUSCRIPT
and V2 along the solution of (8) are
CR
IP
T
dV1
dx
= (x − x∗ )
= (x − x∗ ) − (x − x∗ ) +
dt
xdt
a(1 − m)y ∗
a(1 − m)y
,
−
1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ 1 + b(1 − m)x + cy + f (1 − m)xy
dV2
dy
= η(y − y ∗ )
= η(y − y ∗ ) −
dt
ydt
e(1 − m)x∗
e(1 − m)x
.
+
1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ 1 + b(1 − m)x + cy + f (1 − m)xy
From the time derivative of (9) and a straightforward calculation, we obtain
dV
dt
a(1 − m)2 y ∗ (b + ay)
ηe(1 − m)x∗ (c + f (1 − m)x)
(x − x∗ )2 −
(y − y ∗ )2
ξ(x, y)
ξ(x, y)
ηe(1 − m)f (1 + cy ∗ ) − a(1 − m)(1 + a(1 − m)x∗ )
+
(x − x∗ )(y − y ∗ ).
ξ(x, y)
= −(x − x∗ )2 +
AN
US
With the help of Theorems (2.4) and (2.6) and choosing η in such a way that ηe(1 −
m)(1 + cy ∗ ) = a(1 − m)(1 + (1 − m)ax∗ ), we obtain
h
h
dV
b + aM2 ) i
(c + f (1 − m)m1 ) i
< − 1 − a(1 − m)2
(x − x∗ )2 − ηe(1 − m)x∗
(y − y ∗ )2 .(10)
dt
ξ(x, y)
ξ(x, y)
Using the value of ξ(x, y), coefficient of (x − x∗ )2 takes the following form:
M
a(1 − m)2 (b + aM2 )
(1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )(1 + b(1 − m)x + cy + f (1 − m)xy)
(1 − m)(b + aM2 )(1 − x∗ )
,
≤ −1 + (a1 + c1 L)
ay ∗
ED
− 1+
ay ∗
. Thus we conclude that if all the conditions
(1 − x∗ )(b + aM2 )
dV
of Theorem (2.10) are satisfied, then
< 0 along all the trajectories in R2+ except E ∗ .
dt
Therefore, E ∗ is globally asymptotically stable.
Spatiotemporal model and dynamics
AC
3
CE
PT
which is negative if m > 1−
In this section, to capture the dynamics and inhomogeneous distribution of predator and
prey, we discuss the diffusive model (3).
3.1
Global existence and boundedness of a positive solution of
(3)
Definition 3.1. A vector of the form α = (α1 , α2 , . . . , αn ) where each component αi is
a nonnegative integer, is called a multilinear of order |α| = α1 + α2 + . . . + αn . For any
∂ |α|
α
multiindex α, one can define, D U (x) =
.
∂xα1 1 , ∂xα2 2 , . . . ∂xαnn
7
ACCEPTED MANUSCRIPT
Definition 3.2. Sobolev Space: Let 1 ≤ p < ∞. The Sobolev space W k,p (U ) consists of
locally summable functions u : U → R such that for each multiindex α with |α| ≤ k, Dα u
exists in weak sense and belongs to Lp (U ). Here if p = 2 use notation H k (U ) = W k,2 (U )
(k = 0, 1, 2...).
The system of coupled partial differential equations (3) can be written as
∂w
= G(w) + M ∆w,
∂t
(11)
CR
IP
T
where w = (u, v)T and M = diag{d1 , d2 }, G(w) = (G1 (u, v), G2 (u, v)), G1 (u, v) = u (1 −
a(1−m)y
e(1−m)u
u) − 1+b(1−m)u+cv+f (1−m)uv , G2 (u, v) = v − d + 1+b(1−m)u+cv+f (1−m)uv .
Moreover, let H 2 (Ω) denotes Sobolev space of all square integrable functions such
that all the derivatives up to second order are also in L2 (Ω). Here the initial conditions
u(X, 0) = u0 (X), and v(X, 0) = v0 (X), belongs to Sobolev space H 2 (Ω) ∪ 0. We define
Aw = −M ∆w for w ∈ D(A), where A(w) = {w : w ∈ H 2 (Ω)×H 2 (Ω)}. Hence using this
notation, the system (11) can be rewritten in the following abstract differential equation:
AN
US
∂w
+ Aw = G(w).
(12)
∂t
One can easily establish the existence of a local solution for (12) and the non-negativity
of components u and v of the solution of initial-boundary value problem (3). The proof
of uniform boundedness of the solution of (3) suffices the global existence of the positive
solution. For this, we recall the following lemma:
M
Lemma 3.3. Let w(X, t) satisfies the following equation
PT
ED
∂w
− ∆w = w(1 − u) X ∈ Ω,
∂t
∂w
= 0 X ∈ ∂Ω, t > 0,
∂n
u(X, 0) = u0 (X) > 0, X ∈ Ω.
Then lim u(X, t) = 1 for each X ∈ Ω.
t→∞
AC
CE
Theorem 3.4. (Weak Persistence) Let Ω̄ = Ω ∪ ∂Ω. The non-negative solution (u, v)
satisfy
lim sup max u(., t) ≤ 1, lim sup max v(., t) ≤
t→∞
u∈Ω̄
t→∞
u∈Ω̄
(1 − m)(e − bd)
.
d
(13)
Proof. Using positivity of u and v, we obtain
∂u
= d1 ∆u + u(1 − u).
∂t
Applying Lemma 3.3, ∃ a T > 0 and an arbitrary small > 0 such that u(X, t) < 1 + in Ω × [T, ∞). Similarly, from second equation of (3), after using the bound of u, one can
also estimate the bound for v
lim sup max v(., t) ≤
t→∞
v∈Ω̄
e(1 − m)(1 + ) − bd(1 − m)
.
cd
8
ACCEPTED MANUSCRIPT
Moreover, the arbitrariness of gives that
lim sup max v(., t) ≤
v∈Ω̄
t→∞
(1 − m)(e − bd)
.
cd
The above discussion can also be summarized in the following theorem:
Theorem 3.5. The system (3) has a unique, nonnegative and bounded solution w =
(u, v)T such that
CR
IP
T
w ∈ C([0, ∞), X) ∩ C 1 ([0, ∞), X) ∩ C([0, ∞), D(A)).
Theorem 3.6. (Strong Persistence) The system (3) is strongly persistent i.e.
lim inf max u(., t) > 0 and lim inf max v(., t) > 0 if
t→∞
t→∞
u∈Ω̄
v∈Ω̄
(1 − m)(e − bd)(c − a(1 − m)) > d,
o
n
(1−m)(e−bd)(c−a(1−m))−d
.
where 0 = min c−a(1−m)
,
c
c+f (c−a(1−m))d
AN
US
c > a(1 − m),
Proof. Proof follows similarly as in Theorem 3.4.
3.2
M
Remark 3.7. Theorems 3.4 and 3.6 also ensure that conditions and bounds for the persistence of the temporal system (system (1) and Theorem 2.6) are same as results obtained
for the system with diffusion (system (3) Theorem 3.4).
Homogeneous steady state and its stability analysis
ED
One can easily establish the existence and uniqueness of a positive homogeneous steady
state U ∗ = (u∗ , v ∗ ) of the system (3). To investigate the local stability of U ∗ and the
effect of diffusion on the system (3), we linearize the system (3) around U ∗ . Substituting
T
PT
u(X, t), v(X, t)
= U ∗ + Υ(X, t) = U ∗ +
T
U (X, t), V (X, t)
,
(14)
AC
CE
in the system (3) and picking up all the terms which are linear in U and V, the system
(3) takes the following from:
∂Υ
= D∆Υ + J(U ∗ )Υ,
∂t
where D = diag(d1 , d2 ). More precisely,
∂U
= j11 U + j12 V + d1 ∆U,
∂t
∂V
= j21 U + j22 V + d1 ∆V,
∂t
(15)
where jij for i, j = 1, 2 are same as obtained for the model system without diffusion (refer
2
the section 2.3). Let λi ∞
i=1 denotes the set eigenvalues of −∆ on Ω ∈ R with zero flux
boundary condition such that 0 < λ1 < λ2 < · · · . Further, let χ(λi ) denote the eigenspace
corresponding to λi , and ψij , for j = 1, 2, 3, · · · , dimχ(λi ), be an orthonormal basis of
9
ACCEPTED MANUSCRIPT
∂v
∂u
χ(λi ). Let Hij = {cψij ; c ∈ R2 } and Y = {(u, v) ∈ C 1 (Ω̄) × C 1 (Ω̄) : ∂n
= ∂n
= 0}.
Ldimχ(λi )
L∞
Consider the following direct sum Xi = i=1
Hij and Y = i=1 Xi . Then the space
Xi is invariant under the operator P = D∆ + J(U ∗ ) for each i ≥ 1. Moreover, λ is an
eigenvalue of P if and only if it is also an eigenvalue of
−d1 λi + j11
j12
∗
∗
Ai = J(λi , U ) = −λi D + J(U ) =
.
j21
−d2 λi + j22
The eigenvalues of Ai are determined by
where
(16)
CR
IP
T
η 2 − Bη + C = 0,
B = −λi (d1 + d2 ) + j11 + j22 = −λi (d1 + d2 ) + trJ(U ∗ ),
C = d1 d1 λ2i − (d1 j22 + d2 j11 )λi + detJ(U ∗ ).
(17)
(18)
AN
US
It may be noted here that J(U ∗ ) is same as JE ∗ (refer the Jacobian matrix associated
with temporal system). From equations (16)-(18) and using Routh-Hurwitz criteria [56],
one can summarize the above discussion in the following theorem:
Theorem 3.8. (i) The positive equilibrium U ∗ of the system (3) is locally asymptotically if B > 0 and C > 0.
M
(ii) If j11 < 0 (refer the subsection 2.3), then the positive equilibrium is locally asymptotically stable in the presence as well as in the absence of diffusion.
ED
Remark 3.9. The sufficient condition (Section 2.3) for the asymptotic stability of positive equilibrium E ∗ of the temporal system (2) also implies the asymptotic stability of
homogeneous steady state solution U ∗ of the diffusive system (3).
In the next theorem, we discuss the global stability of positive constant steady-sate
U of (3). The existence of globally stable constant steady-state signifies two important
points: (i) Predator and prey species can maintain steady state for all future time under
certain parametric restriction. (ii) It also ensures the non-existence of spatial patterns.
PT
∗
AC
CE
Theorem 3.10. The positive homogeneous steady state U ∗ of the system (3) is globally
asymptotically stable if the corresponding equilibrium solution of (2) is globally asymptotically stable.
Proof. To prove the Theorem 3.10, we consider the following Lyapunov functional
Z
V (t) =
W1 (u(X, t), v(X, t))dX,
where
Ω
u
W1 (u, v) = u − u − u log ∗ + K
u
∗
∗
Z
v
v∗
η − v∗
dη.
η
Here K is a positive constant and is defined in discussion below. Differentiating W1 with
∗
v−v ∗
1
1
= u−u
, ∂W
= ve(1−m)
. The time derivative of V (t)
respect to u and v, we obtain ∂W
∂u
u
∂v
10
ACCEPTED MANUSCRIPT
CR
IP
T
along the trajectory of the system (3) gives
Z ∂W1 ∂u ∂W1 ∂v dV
=
+
dX
dt
∂u ∂t
∂v ∂t
Ω
Z i
u − u∗ h
a(1 − m)uv
=
d1 ∆u + u(1 − u) −
u
1 + b(1 − m)u + cv + f (1 − m)uv
Ω
i
∗ h
v−v
e(1 − m)uv
+
d2 ∆v − dv +
dX
e(1 − m)v
1 + b(1 − m)u + cv + f (1 − m)uv
Z v − v∗
u − u∗
d1 ∆u +
d2 ∆v
=
u
e(1 − m)v
Ω
Z
h
i
e(1 − m)v
∗
(u − u ) (1 − u) −
+
1 + b(1 − m)u + cv + f (1 − m)uv
Ω
i
∗ h
v−v
e(1 − m)u
+
−d+
dX
(19)
e(1 − m)
1 + b(1 − m)u + cv + f (1 − m)uv
= I1 (t) + I2 (t),
AN
US
where I1 (t) and I2 (t) denote first and second integral in (19), respectively. Using homogeneous boundary conditions, we obtain
Z u − u∗
v − v∗
I1 (t) =
d1 ∆u +
d2 ∆v dX
u
e(1 − m)v
Ω
Z d2 v ∗
d1 u∗
2
2
|∇u| +
|∇v| dX ≤ 0.
=−
(20)
u2
cdv 2
Ω
ED
M
Further, after some straightforward calculations and following the steps involved in the
proof of the Theorem 2.10, one can also show that I2 (t) < 0. This completes the proof of
Theorem 3.10.
Nonexistence of nonconstant positive steady state
CE
3.3
PT
Remark 3.11. The Theorem 3.10 also concludes that if all the conditions of the Theorem
2.10 are well satisfied, then the homogeneous steady state of the associated spatiotemporal
model (3) is globally asymptotically stable.
AC
For the non-existence of non-constant positive equilibrium of the system (3), we establish
the non-existence of non-constant positive equilibrium of the following second order elliptic
system:
−d1 ∆u = u (1 − u) −
−d2 ∆v = v − d +
a(1 − m)v
,
1 + b(1 − m)u + cv + f (1 − m)uv
e(1 − m)u
1 + b(1 − m)u + cv + f (1 − m)uv
with X ∈ Ω and zero-flux boundary condition
∂u
∂v
=
= 0, X ∈ ∂Ω.
∂n
∂n
11
,
(21)
ACCEPTED MANUSCRIPT
Before proving main theorem for non-existence of non-constant positive steady state,
we introduce two important lemmas and one theorem related to priori lower and upper
bounds for the positive solution of (21).
Lemma 3.12. (Harnack Inequality) [57] Let m(X) ∈ C(Ω). Further let ϕ(X) ∈ C 2 (Ω) ∩
C 1 (Ω) be a positive solution of ∆ϕ(X) + m(X)ϕ(X) = 0 in Ω subject to ∂ϕ
= 0 on ∂Ω.
∂ν
Then ∃ a positive constant M∗ = M∗ (km(X)k∞ ) such that max ϕ(X) ≤ M∗ min ϕ(X).
Ω
Ω
(a) If π ∈ C 2 (Ω) ∩ C 1 (Ω) satisfies
∆π(X) + ϕ(X, π(X) ≥ 0 in Ω,
and π(X0 ) = max π(X), then ϕ(X0 , π(X0 )) ≥ 0.
Ω
∂π
≤ 0 on ∂Ω,
∂ν
AN
US
(b) If π(X) ∈ C 2 (Ω) ∩ C 1 (Ω) satisfies
CR
IP
T
Lemma 3.13. [57] Let ϕ ∈ C(Ω × R).
∆π(X) + ϕ(X, π(X) ≤ 0 in Ω,
∂π
≥ 0 on ∂Ω,
∂ν
(22)
(23)
and π(X0 ) = min π(X), then ϕ(X0 , π(X0 )) ≤ 0.
Ω
M
Remark 3.14. The standard regularity theory of elliptic system [61] ensures that the
positive solution of (21) is contained in C 2 (Ω) × C 2 (Ω). Hence one can use the above
lemmas to find out some proper estimates (lower and upper bounds) for the solution of
(21).
(i) Any positive solution (u, v) of (21) satisfy
ED
Theorem 3.15.
max u(X) ≤ 1,
Ω
(1 − m)(e − bd)
,
cd
(24)
PT
Ω
max v(X) ≤
CE
(ii) Assume that d∗ be any fixed positive constant such that d1 , d2 ≥ d∗ and e(1 − m) >
d(1 + b(1 − m)). Then ∃ a positive constant Q = Q(a, b, d, d∗ , e, m) such that
min u(X) ≥ Q,
Ω
min v(X) ≥ Q.
(25)
Ω
AC
Proof. Both the inequalities of (24) are followed immediately by using Lemma 3.13 and
equation (21).
Furthermore, we denote
u(X1 ) = min u(X),
u(X2 ) = max u(X),
(26)
v(Y1 ) = min v(Y ),
v(Y2 ) = max v(Y ).
(27)
Ω
Ω
Ω
Ω
12
ACCEPTED MANUSCRIPT
Using the Lemma 3.13, equations (26) and (27) in (21), we obtain the following inequalities:
a(1 − m)v(X1 )
,
1 + b(1 − m)u(X1 ) + cv(X1 ) + f (1 − m)u(X1 )v(X1 )
e(1 − m) − d(1 + b(1 − m)) u(Y1 )
v(Y1 ) ≥
,
d(1 + f (1 − m)u(Y1 ))
e(1 − m) − d(1 + b(1 − m)) u(Y2 )
,
v(Y2 ) ≤
d(1 + f (1 − m)u(Y2 ))
1 − u(X1 ) ≤
(28)
(29)
(30)
AN
US
CR
IP
T
where e(1 − m) > d(1 + b(1 − m)). Equations (29) and (30) together with (26) and (27)
imply that
e(1 − m) − d(1 + b(1 − m)) u(Y1 )
v(Y1 ) ≥
d(1 + f (1 − m)u(Y1 ))
e(1 − m) − d(1 + b(1 − m)) u(X1 )
≥
,
(31)
d(1 + f (1 − m)u(X1 ))
e(1 − m) − d(1 + b(1 − m)) u(Y2 )
v(Y2 ) ≤
d(1 + f (1 − m)u(Y2 ))
e(1 − m) − d(1 + b(1 − m)) u(X2 )
≤
.
(32)
d(1 + f (1 − m)u(X2 ))
Using (27) and (32) in (28), we obtain
M
1 ≤ u(X1 ) + M ∗ u(X2 ),
(33)
e(1 − m) − d(1 + b(1 − m))
.
d
h
i
a(1−m)v
−1
Now we define m(X) = d1 1 − u − 1+a(1−m)u+cv+f (1−m)uv . Then applying the Harnack
inequality (Lemma 3.12) on ∆u(X) + m(X)u(X) = 0 (from first equation of (21)), we
have
(34)
max u(X) ≤ Q1 min u(X), i.e., u(X2 ) ≤ Q1 u(X1 ),
PT
ED
where M ∗ = a(1 − m)
X∈Ω
X∈Ω
AC
CE
where Q1 is a positive constant which depends on km(X)k∞ .
Combining (34) with (33), we have
min u(X) = u(X1 ) ≥ (1 + M ∗ Q1 )−1 := Q2 .
(35)
X∈Ω
From (31), we have
min v(y) = v(y1 ) ≥
y∈Ω
M ∗ (1 + M ∗ Q)−1
:= Q3 .
a(1 − m)(1 + f (1 − m)(1 + M ∗ Q)−1 )
(36)
Defining Q = min{Q2 , Q3 }, we have min u(X) = u(X1 ) ≥ Q, min v(Y ) = v(Y1 ) ≥ Q.
Theorem 3.16. Suppose that d2 ∗ > e(1−m)(1+b(1−m))−d
be a fixed positive constant. Then
λ1
one can obtain a positive constant d1 ∗ such that the model system (21) has no non-constant
positive solution provided d1 > d1 ∗ and d2 > d2 ∗ .
13
ACCEPTED MANUSCRIPT
Proof. Let us denote F1 (u, v) = u(1 − u) −
a(1−m)uv
,
1+b(1−m)u+cv+f (1−m)uv
e(1−m)uv
.
1+b(1−m)u+cv+f
(1−m)uv
R
1
ũ = Ω Ω u(X)dX, ṽ
F2 (u, v) = −dv +
Let (u, v) be a positive solution of (21). Moreover, we define
R
= Ω1 Ω v(X)dX, for any u, v ∈ L1 (Ω) and κ(u, v) = (1 + b(1 −
m)u + cv + f (1 − m)uv)(1 + b(1 − m)ũ + cṽ + f (1 − m)ũṽ). Integrating the first equation
of (21) after multiplying by (u − ũ), we obtain
Z
2
F1 (u, v)(u − ũ)dX = {F1 (u, v) − F1 (ũ, ṽ)}dX
Ω
ZΩ h
uv
=
(u − ũ) − (u2 − ũ2 ) − a(1 − m)
1 + b(1 − m)u + cv + f (1 − m)uv
Ω
i
ũṽ
dX
−
1 + b(1 − m)ũ + cṽ + f (1 − m)ũṽ
Z
h (v + cvṽ)(u − ũ)2
≤
(u − ũ)2 − a(1 − m)
κ(u, v)
Ω
(ũ + b(1 − m)uũ)(u − ũ)(v − ṽ) i
+
dX.
(37)
κ(u, v)
d1 |∇(u − ũ)| dX =
CR
IP
T
Ω
Z
Z
AN
US
Using the bounds of u and v from (25) and P1 = ab in (37), we have
Z
Z h
i
2
d1 |∇(u − ũ)| dX ≤
(u − ũ)2 + 2P1 |u − ũ||v − ṽ| dX.
Ω
(38)
Ω
Similarly integrating the second equation of (21) after multiplying by (v − ṽ), we obtain
Z
Z
F2 (u, v)(u − ũ)dX = {F2 (u, v) − F2 (ũ, ṽ)}dX
Ω
Z h
uv
=
− d(v − ṽ) + e(1 − m)
1 + b(1 − m)u + cv + f (1 − m)uv
Ω
i
ũṽ
−
dX.
(39)
1 + b(1 − m)ũ + cṽ + f (1 − m)ũṽ
d2 |∇(v − ṽ)| dX =
Ω
PT
ED
Ω
2
M
Z
2
AC
CE
Now by defining P2 = 1 + c(1 − m)2 (e−bd)
and using the bounds of u and v from (24) in
c2 d 2
(39), we obtain
Z
Z 2
d2 |∇(v − ṽ)| dX ≤
(e(1 − m)(1 + b(1 − m)) − d)(v − ṽ)2
Ω
Ω
+e(1 − m)P2 (u − ũ)(v − ṽ) dX.
(40)
Let us define R = P2 + 2P1 . Now from inequalities (38) and (40), we have
Z
Z
Z 2
2
{d1 |∇(u − ũ)| + d2 |∇(v − ṽ)| }dX ≤
(u − ũ)2 + (e(1 − m)(1 + b(1 − m))
Ω
Ω
2
−d)(v − ṽ) + 2R|u − ũ||v − ṽ| dX. (41)
14
ACCEPTED MANUSCRIPT
Now, using the Poincare inequality, the above inequality (41) yields
Z
Z
Z
R
2
2
d1 λ1 |∇(u − ũ)| + d2 λ1 |∇(v − ṽ)| dX ≤
1+
(u − ũ)2 + e(1 − m)(1
δ
Ω
Ω
+b(1 − m)) − d + δR (v − ṽ)2 dX,
(42)
CR
IP
T
for an arbitrary positive constant δ. As d2 λ1 > e(1 − m)(1 + b(1 − m)) − d, one can find a
sufficiently small δ ∗ such
that d2 λ1 ≥ e(1 − m)(1
+ b(1 − m)) − d + δ ∗ R. Moreover, if we
∗
define d1 = λ1 1 + δR∗ , then d1 λ1 − 1 + δR∗ > 0 (refer the condition of Theorem 3.16).
Hence one can conclude that u = ũ and v = ṽ.
3.4
AN
US
Remark 3.17. The Theorem 3.16 determines threshold values of diffusivity constants for
the non-existence of a non-constant positive solution of the system (21). The condition
d2 ∗ > e(1−m)(1+b(1−m))−d
also ensures that the prey reserve parameter (m) affects the nonλ1
existence of a non-constant positive solution which in turn implies that the prey refuge
may leave some effect for pattern formation.
Turing space and Hopf instabilities
As we know that any solution of the linear system (15) can be expanded into a Fourier
series. Hence
γ1 mn (t) sin k(km , kn )X,
m,n=0
∞
X
m,n=0
∞
X
∞
X
γ2 mn (t) sin k(km , kn )X,
∞
X
unm (X, t) =
ED
M
U (X, t) =
V (X, t) =
m,n=0
vnm (X, t) =
(43)
m,n=0
AC
CE
PT
where km and kn are corresponding wavenumbers. Now, replacing U and V from unm vnm
into (15), wee have
dγ1 nm
= (j11 − d1 k 2 )γ1 nm + j12 γ2 nm ,
dt
(44)
dγ2 nm
= j21 γ1 nm + (j22 − d2 k 2 )γ2 nm ,
dt
2
where k 2 = kn2 + km
. Equation (44) is a system of two linear differential equations and a
general solution of (44) has the form A exp(λ1 t) + B exp(λ2 t). Here the constants A and B
are determined with the help of initial conditions (4) while λ1 and λ2 are the eigenvalues
of the associated coefficient matrix
J=
j11 − d1 k 2
j12
j21
j22 − d2 k 2
15
.
ACCEPTED MANUSCRIPT
Here λ1 and λ2 are the solutions of the characteristic equation associated with the matrix
J
λ2 − tr(J) + det(J),
(45)
where tr(J) = tr(JE ∗ )−(d1 +d2 )k 2 and det(J) = det(D)k 4 +(−d2 j11 −d1 j22 )k 2 +det(JE ∗ ).
It is obvious that the diffusive instability holds, if tr(J) < 0 and det(J) < 0 and j11 < 0
are satisfied. Moreover, if j11 < 0 holds, then tr(JE ∗ ) < 0 and hence tr(J) < 0. The
condition det(J) < 0 is sufficient for the diffusive instability and which is equivalent to
the inequality
CR
IP
T
S(k 2 ) = det(D)k 4 + (−d2 j11 − d1 j22 )k 2 + det(JE ∗ ) < 0.
(46)
Here S is a quadratic function with respect to k 2 . Moreover the minimum value of S(k 2 )
is obtained at
1
2
d1 j22 + d2 j11 .
(47)
k 2 = kmin
=
2 det(D)
AN
US
2
Hence more stronger condition for diffusive instability becomes S(kmin
) < 0 i.e.,
1
d2 j11 + d1 j22 > 2(det(D)det(JE ∗ ) 2 .
(48)
As a consequence, we have the following theorem.
Theorem 3.18. The system (3) undergoes diffusion-induced instability if the condition
of Theorem 2.9, relations (46) and (48) hold, i.e.,
1
(49)
Numerical simulations
ED
4
M
d2 j11 + d1 j22 > 2(det(D)det(JE ∗ ) 2 .
PT
In order to show the analytical findings and stability results obtained in the previous
sections graphically, we numerically simulate the solutions of both the model systems
(diffusive and non-diffusive) with different sets of parameter values.
AC
CE
Case 1: The Non-diffusive System:
Here we execute two interesting aspects, namely coexistence and dynamics of the system
(2) with respect to two different significant parameters, viz., prey refuge (m) and maximum predation rate (a). The parameter values for the system (2) are listed in the Table
1. First fixing the parameter a = 20 and other parameters from Table 1, when we vary
prey refuge parameter m throughout its whole range [0, 1) (see bifurcation diagram with
respect to m: Fig.1(i)), then we observe that at the lower range of m coexistence equilibrium point (E ∗ ) of the system (2) is unstable focus and eventually system converges
to nearest stable limit cycle around E ∗ (Fig.1(iα): phase-plane diagram with m = 0.1).
At the intermediate range of m, system (2) converges to its coexistence equilibrium point
(E ∗ ) (Fig.1(iβ): phase-plane diagram with m = 0.6). Here, the system (2) changes its
unstable focus to stable focus status via Hopf-bifurcation [58, 59, 60] with respect to
parameter m. At the higher range of m, due to the lack of food (prey), predator goes
to extinction while prey goes to its highest density (carrying capacity) by gradually decreasing foraging efficiency of predator (Fig.1(iγ: phase-plane diagram with m = 0.97).
16
ACCEPTED MANUSCRIPT
CR
IP
T
Now coexistence equilibrium point (E ∗ ) becomes infeasible (ecological point of view) and
predator-free equilibrium point (E 1 ) becomes stable node. In this case, a transcritical
bifurcation occurs between E ∗ and E 1 [62]. On the other hand, keeping fixed the parameter m = 0.01 and other parameters from Table 1, if we vary maximum predation rate
parameter a throughout a large range [2, 20] (refer the bifurcation diagram with respect
to a: Fig.1(ii)), then we notice that at the lower range of a, system (2) converges to its
coexistence equilibrium point (E ∗ ) (Fig.1(iiα): phase-plane diagram with a = 4). Moreover, at higher range of a, coexistence equilibrium point (E ∗ ) of the system (2) becomes
an unstable focus and eventually system converges to nearest stable limit cycle around
E ∗ (Fig.1(iiβ): phase-plane diagram with a = 15). Here, system (2) changes its stable
focus to unstable focus status via Hopf-bifurcation with respect to parameter a.
Case 2: The Diffusive System:
M
AN
US
For the diffusive system, we will open up discussion from three different perspectives:
(1) suppose that there is no prey refuge (i.e., m = 0) and observe the effect of prey
diffusion on the pattern formation; (2) assume that there is some prey refuge (i.e., m 6= 0)
and examine the effect of prey diffusion on the pattern formation; (3) observation on the
system with different intensities of prey refuge along with other system parameters remain
fixed. In every situation, the relation between S(k 2 ) and k is shown and the corresponding
pattern structure are also demonstrated.
In the process of numerical simulation, we select the Neumann boundary condition
with a size of 200 × 200. Generally speaking, Neumann boundary is also called as zeroflux boundary, which implies that the domain is isolated or insulated from the external
environment. We assume that the space step is dx = 1 and the time step is dt = 0.005.
The initial distributions of both prey and predator are taken as follows:
v(x, y, 0) = v ∗ + ξ2 ,
ED
u(x, y, 0) = u∗ + ξ1 ,
(50)
CE
PT
where ξi (i = 1, 2) is a tiny random perturbation term. We run the simulations until
pattern reaches a stationary state or does not seem to change its main characteristics any
longer. Simultaneously, we fix d2 = 50 and other parameter values are given in Table
1 except the parameters a, m, d1 . Here, we just show the pattern formation of prey
population (u) in the xy-plane since predator population (v) exhibits extremely similar
distributions.
AC
Case 2.1: Without refuge and the effect of diffusion on pattern structure
In this subsection, we suppose that the prey have no any refuges (m = 0) and observe
the effect of prey diffusion on the pattern formation. From the equation (46), it is obvious
that the condition for the diffusion-induced instability is equivalent to S(k 2 ) < 0. According to the Fig. 2(i), it is clear that the interval of negativity of the function S(k 2 ) becomes
more smaller as d1 increases i.e., the corresponding diffusive system will be unstable only
when the prey diffusion coefficient is smaller and Turing pattern emerges under only this
circumstance. The pattern evolution with different prey diffusion coefficients are shown
in Figure 3. When the value of d1 is smaller, the center of region is occupied by shortstripe pattern and surrounded by hot-spot pattern with high density (Fig. 3(i)). If d1
takes a middle value, the density of center is uniform and the uniform density region is
17
ACCEPTED MANUSCRIPT
surrounded by fewer short-stripe pattern and the outermost part of the area is hot-spot
pattern with intermediate density (Fig. 3(ii)). When the value of d1 is high, the pattern
structure is same as the above situation, but the uniform density central region is more
larger and the hot-spot pattern becomes more sparse (Fig.3(iii)). It is obvious that the
diffusion coefficient of the prey reduces the range of pattern region, that is to say, the
diffusion coefficient have a stabilizing effect on the spatiotemporal system (3).
Case 2.2: With refuge and the effect of diffusion on pattern structure
PT
ED
M
AN
US
CR
IP
T
In this case, we consider a constant prey refuge (m = 0.1) and examine the effect
of prey diffusion on the pattern structure. The plot of S(k 2 ) as a function is shown in
Fig.2(ii), which is similar to Fig.2(i). The only difference is that both the interval of
negativity of the function S(k 2 ) and the corresponding functional value becomes more
smaller. Namely, when the prey refuge is taken into account, prey diffusion can still give
rise to Turing instability, however the prey refuge seemingly reduces the instability probability. The associated pattern formation is shown in Fig.4. When d1 is smaller, one can
observe that the density of central zone is uniform and surrounded by hot-spot pattern
(Fig.4(i)). If we increase d1 slightly, analogous-spiral pattern emerges and it is encompassed by cold-spot pattern (Fig.4(ii)). When d1 is higher, the cold-spot pattern fades
away and the region is occupied by analogous-spiral pattern (Fig.4(iii)). Moreover, the
central region appears the most safest region for the prey and the farther from the center,
the lower the security seems. From the perspective of the preys and predators, in order
to avoid being attacked by predators, the preys migrate towards the center. Moreover,
when the preys are protected at a certain percentage m, (1 − m) of preys is only available
to predator. The protective effect of population aggregation will lose a certain superiority
and the preys have to escape at a higher speed to avoid being eaten. This ensures that the
diffusion coefficient of preys will be higher because Turing pattern is induced by diffusion,
or more strictly speaking, it is caused by the value of dd12 . When the refuge is considered,
the value of dd21 becomes smaller, to some extent, which reduces the possibility of Turing
pattern. It is in complete agreement with the result of the Fig.2(ii).
Case 2.3: The effect of the refuge on pattern structure
AC
CE
In this situation, we fix the diffusion coefficient of prey at d1 = 0.01 and vary m from
0.01 to 0.1 and then we examine the effect of prey refuge on the pattern structure. The
plot of S(k 2 ) with different values of refuge m is shown in Fig. 2(iii). The Fig. 2(iii)
demonstrates that the interval of negativity of the function S(k 2 ) becomes more smaller
as m increases. Furthermore, if the value of m is greater than 0.14, the functional value
of S(k 2 ) will be positive and the Turing instability does not occur. Comparing Fig. 2(iii)
with Fig. 2(i), we can conclude that the prey refuge and diffusion have the parallel impact
on the diffusion-induced instability. The stationary-state patterns for different values of
m are shown in Fig. 5 and the pattern structure is almost consistent to Fig.3. This
perfectly illustrates that prey refuge and diffusion coefficient have the parallel impact
on pattern formation and structure. Thus one can summarize the above discussion by
saying that prey reserve and diffusive coefficients both leave an stabilizing effect on the
spatiotemporal system (3).
18
ACCEPTED MANUSCRIPT
5
Discussion
AC
CE
PT
ED
M
AN
US
CR
IP
T
Many efforts have been put to manage sub-population of the vulnerable cheetas (Acinonyx
jubatus) in small reserves. Like vulnerable cheetas there are many other species those
are either driven to extinction or at the verge of extinction, due to overexploitation,
over predation, indiscriminate harvesting, environmental pollution etc. There are several
biological examples where individual prey can decrease their predation risk by using a
refuge, which may be habitat that is spatially separate, a phase of prey life history that is
temporarily separate, a collective behaviour of prey such as aggregation into herds. This
motivates us for the creation of biospheres reserve zones/refuges and this may decrease
the interaction among species. While creating protected areas for a species several factors
should be taken into consideration like it’s capacity to house number of individuals of the
species to be protected, dynamics of the ecosystems supporting these species, economics
associated with the maintenance of the protected area, bio-economics of the ecosystem
and several others. The establishment of protected areas (reserve zones) is the potential
solution for optimal harvesting problems.
In this paper, we have investigated the dynamical behaviors of a prey-predator system (both for diffusive and non-diffusive) with Crowley-Martin functional response and
prey reserve. The prey refuge is also a significant factor as far as the size of habitat and
co-existence of species are concerned. The temporal system is thoroughly examined with
different types of bifurcations. Parametric conditions of co-existence and local asymptotic stability of each biomass equilibrium are presented systematically. We have derived
threshold values of prey reserve parameter (m) for the permanence and global stability
of the temporal system (2) (refer the Theorems 2.6 and 2.10). The parameter ranges of
existence of multiple bifurcations are identified. More precisely, bifurcation diagram with
respect to prey reserve (m) and maximum predation rate (a) ensures the occurrence of
transcritical bifurcation and Hopf-bifurcation (refer the Figure 1).
A rigorous investigation of the global dynamics and bifurcations of patterned solutions
of the spatio-temporal counterpart of the prey-predator system under zero-flux boundary
conditions is also presented. We have pointed out the similarity/difference between the
sufficient parametric conditions obtained for persistence, local and global stability of both
the temporal and its diffusive counterpart. The theoretical analysis reveals the fact that
global stability of the interior equilibrium solution (E ∗ ) of non-diffusive system (2) is
sufficient for the global asymptotic stability of positive homogeneous state (U ∗ ) of diffusive
system (3). The role of prey reserve in the non-existence of non-constant solution has
been investigated which suggests that prey reserve must leave some effect on pattern
formation (see the Theorem 3.16). We have carried out extensive numerical simulations
to explore the effect of prey reserve (m) and diffusion coefficients on pattern structure. The
obtained results ensure that the prey reserve and diffusion coefficient leave parallel impact
on pattern structure. In particular, the presence of prey reserve reduces the possibility of
Turing pattern.
In brief, the value of the present study can be summarized in three folds. First: it
analyzes the effects of prey refuge and maximum predation rate parameters on the stability
of the non-diffusive system. System shows diverse dynamical phenomena with respect
to these two significant system parameters and system switches its stability via Hopfbifurcation and changes stability with other feasible system equilibrium via transcritical
bifurcation. Second: the spatial system classifies diffusion-driven instability conditions
and gives corresponding theorems. Third: it illustrates spatial patterns close to the
19
ACCEPTED MANUSCRIPT
Acknowledgments
AC
6
CE
PT
ED
M
AN
US
CR
IP
T
onset of instability bifurcation via numerical simulations, indicating that the dynamics
of the diffusive system exhibits a prey refuge and diffusion controlled pattern formation.
Moreover, as our study shows that the effect of prey refuge for pattern formation in two
dimensional spatial domain has profound effect on predator-prey interaction. On one
side, this may enrich the spatial dynamics of the effect of prey refuge on the predatorprey model while on the other side it may provide important information to practical
ecologists about the impact of prey refuge on pattern formation. The stabilizing role of
prey reserve in the dynamics and stability of prey-predator system have been discussed
by few researchers [52, 53] but not much with Crowley-Martin functional response [54].
As the present study deals with combined effects of diffusion and prey refuge, it also
generalizes the existing works in [52, 53, 54]. The results obtained for global stability and
permanence in [52, 53], can easily be obtained from our study just by ignoring mutual
interference among predator species.
Future directions: It should be noted that pattern selection is an important method
to study reaction-diffusion equation, which has played a great role in chemistry, spatial
epidemiology, and ecosystems [63, 64, 65]. Predator-prey relationship is very crucial in
population and abundant pattern structures can easily be found in predator-prey systems,
such as spot, stripe, mixed, labyrinth, etc [65, 66]. In the biological systems, the pattern
formation has a very significant effect on the stability of community and ecosystem. For
instance, the spot pattern indicates that the prey (or predator) exhibits patchy distribution and population aggregation (the preys avoid predation and the predators add the
chance of catching preys). In this sense, for the pattern formation in the predator-prey
system, the two significant issues for further investigations are: (1) previous studies in
predator-prey system have focused on the local role of population, however, the nonlocal
effects of population have a great impact on themselves for the population movement.
As a result, the population model with nonlocal effect should be more reasonable [67].
The travelling wave solutions and entire solutions of the nonlocal models have received
some attention, but no one has ever done the research of pattern dynamics, which will
be a new perspective for population dynamics. (2) In order to maintain the ecosystem
stability, ecologists suggest several suitable measures, such as provision of prey refuge,
supply of additional food and so on. We wonder that what kind of control measures
are most reasonable or profitable. Maybe, we can combine the pattern formation with
control theory of the reaction-diffusion equation to find the optimal control strategy and
the corresponding optimal pattern structures.
The project is funded by the National Natural Science Foundation of China under Grant
11671241, 131 Talents of Shanxi University, Program for the Outstanding Innovative
Teams (OIT) of Higher Learning Institutions of Shanxi, and Natural Science Foundation
of Shanxi Province Grant No. 201601D021002.
References
[1] A.A. Berryman, The origin and evolution of predator-prey theory, Ecology 73(5)
(1992) 1530-1535.
20
ACCEPTED MANUSCRIPT
[2] C. Xu, Y. Shao, Bifurcations in a predator-prey model with discrete and distributed
time delay, Nonlinear Dynamics 67(3) (2012) 2207-2223.
[3] G.-Q. Sun, C.-H. Wang, L.-L. Chang, Y.-P. Wu, L. Li, Z. Jin, Effects of feedback
regulation on vegetation patterns in semi-arid environments, Appl. Math. Model. 61
(2018) 200-215.
[4] C. Xu, X. Tang, M. Liao, M., Stability and bifurcation analysis of a delayed predatorprey model of prey dispersal in two-patch environments, Appl. Math. Comput 216(10)
(2010) 2920-2936.
CR
IP
T
[5] C.S. Holling, Some characteristics of simple types of predation and parasitism,
Canada Entomology 91 (1959) 385-398.
[6] T.K. Kar, H. Matsuda, Global dynamics and controllability of a harvested preypredator system with Holling type III functional response, Nonlinear Analysis: Hybrid Systems 1 (2007) 59-67.
AN
US
[7] S. Liu, E. Beretta, A stage-structured predator-prey model of Beddington-DeAngelis
type , SIAM J. Appl. Math. 66(4) (2006) 1101-1129.
[8] J.R. Beddington, Mutual interference between parasites or predators and it’s effect
on searching efficiency, J. Anim. Ecol. 44(1) (1975) 331-340.
[9] S.B. Hsu, T.W. Hwang, Y. Kuang, Global analysis of the Michaelis-Menten-type
ratio-dependent predator-prey system, J. Math. Biol. 42 (2001) 489-506.
M
[10] C. Xu, P. Li, Oscillations for a delayed predator-prey model with Hassell-Varley-type
functional response, Comptes Rendus Biologies 338(4) (2015) 227-240.
ED
[11] J.P. Tripathi, S.S. Meghwani, M. Thakur, S. Abbas, A modified Leslie-Gower
predator-prey interaction model and parameter identifiability, Commun. Nonlinear
Sci. Numer. Simulat. 54 (2018) 331-346.
PT
[12] S. Ruan, D. Xiao, Global analysis in a predator-prey system with non-monotonic
functional response, SIAM J. Appl. Math. 61(4) (2001) 1445-1472.
AC
CE
[13] J.P. Tripathi, S. Tyagi, S. Abbas, Global analysis of a delayed density dependent
predator-prey model with Crowley-Martin functional response, Commun. Nonlin.
Sci. Numer. Simulat. 30 (2016) 45-69.
[14] M.P. Hassell, Mutual interference between searching insect parasites, J. Anim. Ecol.
40 (1971) 473-486.
[15] D. Jana, J.P. Tripathi, Impact of generalist type sexually reproductive top predator
interference on the dynamics of a food chain model, In. J. Dynam. Control (2016)
DOI: 10.1007/s40437-016-0255-9.
[16] B. Zimmermann, H.K. Sand, P. Wabakken, O. Liberg, H.P. Andreassen, Predatordependent functional response in wolves: from food limitation to surplus killing, J.
Anim. Ecol. 84 (2015) 102-112.
21
ACCEPTED MANUSCRIPT
[17] P.H. Crowley, E.K. Martin, Functional responses and interference within and between
year classes of a dragonfly population, J. N. Amer. Benth. Soc., 8 (1989) 211-221.
[18] G.T. Sklaski, J.F. Gilliam, Functional responses with predator interference: viable
alternative to Holling type II model, Ecology 82 (11) (2001) 3083-3092.
[19] R.D. Parshad, A. Basheer, D. Jana, J.P. Tripathi, Do prey handling predators really
matter: Subtle effects of a Crowley-Martin functional response, Chaos, Solitons and
Fractals 103 (2017) 410-421.
CR
IP
T
[20] R.K. Upadhyay, S.N. Raw, V. Rai., Dynamic complexities in a tri-trophic food chain
model with Holling type II and Crowley-Martin functional response, Nonlin. Anal.:
Model. Control 15(3) (2010) 361-375.
[21] J.P. Tripathi, Almost periodic solution and global attractivity for a density dependent predator-prey system with mutual interference and Crowley-Martin response
function, Differ. Equ. Dyn. Syst. (2016) DOI: 10.1007/s12591-016-0298-6.
AN
US
[22] H. Yin, X. Xiao, X. Wen, K. Liu, Pattern analysis of a modified Leslie-Gower
predator-prey model with Crowley-Martin functional response and diffusion, Comput. Math. Appl. 67 (2014) 1607-1621.
[23] V. Hutson, K. Schmitt, Permanence and the dynamics of biological systems, Math.
Biosci. 111 (1992) 1-71.
M
[24] S. Abbas, M. Banerjee, N. Hungerbuhler, Existence, uniqueness and stability analysis
of allelopathic stimulatory phytoplankton model, J. Math. Anal. Appl. 367 (2010)
249-259.
ED
[25] V. Krivan, D. Jana, Effects of animal dispersal on harvesting with protected areas,
J. Theor. Biol. 364 (2015) 131-138.
PT
[26] X. Guan, W. Wang, Y. Cai, Spatiotemporal dynamics of a Leslie-Gower predatorprey model incorporating a prey refuge, Nonlinear Anal. RWA 12 (2011) 2385-2395.
[27] J. Smith, Models in ecology, Cambridge University Press, Cambridge, 1974.
AC
CE
[28] D. Jana, R. Agrawal, R.K. Upadhyay, Dynamics of generalist predator in a stochastic
environment: effect of delayed growth and prey refuge, Appl. Math. Comput. 268
(2015) 1072-1094.
[29] D. Jana, S. Ray, Impact of physical and behavioral prey refuge on the stability and
bifurcation of Gause type Filippov prey-predator system. Model. Earth Syst. Environ.
2(1) (2016) 2-24.
[30] B. Dubey, A prey-predator model with a reserved area, Nonlinear Analysis: Modelling and Control 12(4) (2007) 479-494.
[31] B. Dubey, P. Chandra, P. Sinha, A model for fishery resource with reserve area,
Nonlinear Analysis. RWA. 4 (2003) 625-637.
22
ACCEPTED MANUSCRIPT
[32] Y. Huang, F. Chen, L. Zhong, Stability analysis of a predator-prey model with
Holling type III response function incorporating a prey refuge, Appl. Math. Comput.
182 (2006) 672-683.
[33] R. Yang, J. Wei, Stability and bifurcation analysis of a diffusive predator-prey system
in Holling type III with a prey refuge, Nonlinear Dynamics 79 (2015) 631-646.
[34] D. Jana, Chaotic dynamics of a discrete predator-prey system with prey refuge, Appl.
Math. Comput. 24 (2013) 848-865.
CR
IP
T
[35] D.L. DeAngelis, R.A. Goldstein, R.V. O’Neill, A model for tropic interaction, Ecological Society of America 56(4) (1975) 881-892.
[36] L.A.D. Rodrigues, D.C. Mistro, S. Petrovskii, Pattern formation, long-term transients, and the Turing-Hopf bifurcation in a space-and time-discrete predator-prey
system, Bull. Math. Biol. 73(8) (2011) 1812-1840.
AN
US
[37] R.S. Cantrell, C. Cosner, Spatial ecology via reaction-diffusion equations, Wiley,
West Sussex, 2003.
[38] M. Banarjee, S. Petrovskii, Self-organized spatial pattern and chaos in a ratio dependent predator-prey system, Theoretical Ecology 4 (2011) 37-53.
[39] G.-Q. Sun, A. Chakraborty, Q.X. Liu, Z. Jin, K.E. Anderson, B.L. Li, Influence of
time delay and nonlinear diffusion on herbivore outbreak, Commun. Nonlinear Sci.
Numer. Simulat. 19(5) (2014) 1507-1518.
M
[40] R.K. Upadhyay, V. Volpert, N.K. Thakur, Propagation of Turing patterns in plankton model, J. Biol. Dynam. 6 (2012) 524-538.
ED
[41] Y. Wang, J. Cao, G.-Q. Sun, J. Li, Effect of time delay on pattern dynamics in a
spatial epidemic model, Physica A 412 (2014) 137-148.
PT
[42] X-P. Yan, C-H Zhang, Stability and turing instability in a diffusive predator-prey
system with Beddington-DeAngelis functional response, Nonlinear Analysis: RWA.
20 (2014) 1-13.
CE
[43] X.C. Zhang, G.-Q., Sun, Z. Jin, Spatial dynamics in a predator-prey model with
Beddington-DeAngelis functional response, Physical Review E 85(2) (2012) 021924.
AC
[44] H. Merdan, S. Kayan, Hopf bifurcations in Lengyel–Epstein reaction–diffusion model
with discrete time delay, Nonlinear Dynam. 79(3) (2015) 1757-1770.
[45] W. Zuo, Y. Song, Stability and bifurcation analysis of a reaction–diffusion equation
with distributed delay, Nonlinear Dynamics 79(1) (2015) 437-454.
[46] A.M. Turing, The chemical basis of mokphogenesis, Philos. Trans. R. Soc. Lond. 237
(1952) 37-72.
[47] M. Banerjee, S. Abbas, Existence and Non-existence of spatial patterns in a ratiodependent predator-prey model, Ecol. Complex. 21 (2015) 199-214.
23
ACCEPTED MANUSCRIPT
[48] M. Banerjee, S. Banerjee, Turing instabilities and spatio-temporal chaos in ratiodependent Holling-Tanner model, Math. Biosci. 236 (2012) 64-76.
[49] G. Zhang, W. Wang, X. Wang, Coexistence states for a diffusive one-prey and twopredators with B-D functional response, Comput. Math. Appl. 387 (2012) 931-948.
[50] L. Xue, Pattern formation in a predator-prey model with spatial effect, Physica A
391 (2012) 5987-5996.
CR
IP
T
[51] J.P. Tripathi, S. Abbas, M. Thakur, Local and global stability analysis of two prey
one predator model with help, Commun. Nonlin. Sci. Numer. Simulat. 19 (2014)
3284-3297.
[52] L. Li, C.-H. Wang, S.-F. Wang, M.-T. Li, L. Yakob, B. Cazelles, Z. Jin, W.-Y. Zhang,
Hemorrhagic fever with renal syndrome in China: Mechanisms on two distinct annual
peaks and control measures, Internat. J. Biomath. 11 (2018) 1850030.
AN
US
[53] J.P. Tripathi, S. Abbas, M. Thakur, Dynamical analysis of a prey–predator model
with Beddington-DeAngelis type function response incorporating a prey refuge, Nonlinear Dynam. 80 (2015) 177-196.
[54] R.S. Cantrell, C. Cosner, On the dynamics of predator-prey models with the
Beddington-DeAngelis functional response, J. Math. Anal. Appl. 257 (2001) 206-222.
[55] H. Li, Y. Takeuchi, Dynamics of the density dependent predator-prey system with
Beddinton-DeAngelis functional response, J. Math. Anal. Appl. 374 (2011) 644-654.
M
[56] C. Xu, Y. Wu, Bifurcation and control of chaos in a chemical system, Appl. Math.
Model. 39(8) (2015) 2295-2310
ED
[57] Y. Lou, W-M. Ni, Diffusion self-diffusion and cross diffusion, J. Differ Equations 131
(1996) 79-131.
PT
[58] C. Xu, Stability and bifurcation analysis in a viral model with delay, Math. Meth.
Appl. Sci. 36(10) (2013) 1310-1320.
CE
[59] R. Curtu, Singular Hopf bifurcations and mixed-mode oscillations in a two-cell inhibitory neural network, Physica D 239(9) (2010) 504-514.
AC
[60] C. Xu, X. Tang, M. Liao, Frequency domain analysis for bifurcation in a simplified
tri-neuron BAM network model with two delays, Neural Networks 23(7) (2010) 872880.
[61] D. Gilbarg, N.S. Trudinger, Elliptic partial differential equations of second order,
Springer-Verlag, New York, 1996.
[62] D. Jana, N. Bairagi, Habitat complexity, dispersal and metapopulations: Macroscopic study of a predator–prey system, Ecol. Complex. 17 (2014) 131-139.
[63] Q. Ouyang, Patterns formation in reaction diffusion systems. Shanghai Sci-Tech.
Education Publishing House, Shanghai, 2000.
24
ACCEPTED MANUSCRIPT
[64] G.-Q. Sun, C.-H. Wang, Z.-Y. Wu, Pattern dynamics of a Gierer–Meinhardt model
with spatial effects, Nonlinear Dynamics 88 (2017) 1385-1396.
[65] G.-Q. Sun, Z.-Y. Wu, Z. Wang, Z. Jin, Influence of isolation degree of spatial patterns
on persistence of populations, Nonlinear Dynam. 83 (2016) 811-819.
[66] A. Yochelis, Y. Tintut, L.L. Demer, A. Garfinkel, The formation of labyrinths, spots
and stripe patterns in a biochemical approach to cardiovascular calcification. New J.
Phys. 10 (2008) 055002.
CR
IP
T
[67] X.-Q. Zhao, J. Borwein, P. Borwein, Dynamical systems in population biology,
Springer, New York, 2003.
AN
US
Parameters
Symbols
Prey saturation constant
b
Predation interference constant
c
Predator mortality rate
d
Conversion factor
e
Predator interference (in the presence of prey handling)
f
Default values
2.10
1.10
0.2
5
0.001
AC
CE
PT
ED
M
Table 1: System parameters and their default values chosen for simulation
25
ACCEPTED MANUSCRIPT
(i) a=20
1
0.9
0.9
(iα) m=0.1
y
0.8
0.8
0
0.7
0
0
0.8
0.7
(iβ) m=0.6
0
0.4
(iiβ) a=15
1
y
y
1
0.6
0.6
0.5
0
0
1
CR
IP
T
0
0.5
(iγ) m=0.97
1
y
Population
(iiα) a=4
1
y
1
(ii) m=0.01
1
x
0.4
0
x
0
1
0.4
x
x
0
1
0.3
0.2
0.2
y
AN
US
0.3
y
0.1
0
0.1
0
0.2
0.4
0.6
0.8
1
0
2
5
10
15
20
a
M
m
AC
CE
PT
ED
Figure 1: Bifurcation diagram of system (2) w.r.t m (Fig.(1)(i)) and three phase-plane diagrams w.r.t lower (Fig.(1)(iα)), intermediate (Fig.(1)(iβ)) and higher (Fig.(1)(iγ)) value
of m keeping a = 20. Bifurcation diagram of system (2) w.r.t a Fig.(1)(ii) and three
phase-plane diagrams w.r.t lower (Fig.(1)(iiα)) and higher (Fig.(1)(iiβ)) value of a keeping m = 0.01. Other parametric values are given in Table 1.
(i)
(ii)
(iii)
Figure 2: The graph of the function S(k 2 ) for parametric values given in Table 1 and for
(i) m = 0, a = 5, d2 = 50, (ii) m = 0.1, a = 5, d2 = 50 and (iii) a = 5, d1 = 0.01, d2 = 50.
26
ACCEPTED MANUSCRIPT
(i)
(ii)
(iii)
(ii)
(i)
AN
US
CR
IP
T
Figure 3: Snapshots of contour pictures of the time evolution of prey (u) in the xy-plane
with d2 = 50 at 300000 iterations where (i): d1 = 0.01; (ii): d1 = 0.02; (iii): d1 = 0.03
and other parametric values are given in Table 1.
(iii)
(i)
CE
PT
ED
M
Figure 4: Snapshots of contour pictures of the time evolution of prey (u) in the xy-plane
with m = 0.1, d2 = 50 at 300000 iterations where (i): d1 = 0.012; (ii): d1 = 0.015; (iii):
d1 = 0.018 and other parametric values are given in Table 1.
(ii)
(iii)
AC
Figure 5: Snapshots of contour pictures of the time evolution of prey (u) in the xy-plane
with d1 = 0.01 and d2 = 50 at 300000 iterations where (i): m = 0.01; (ii): m = 0.05; (iii):
m = 0.1 and other parametric values are given in Table 1.
27
Документ
Категория
Без категории
Просмотров
1
Размер файла
5 842 Кб
Теги
jfranklin, 2018, 029
1/--страниц
Пожаловаться на содержимое документа