# Bifurcation, Chaos, and Pattern Formation for the Discrete Predator-Prey Reaction-Diffusion Model

Authors: Lili Meng; Yutao Han; Zhiyi Lu; Guang Zhang
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

In this paper, a discrete predator-prey system with the periodic boundary conditions will be considered. First, we get the conditions for producing Turing instability of the discrete predator-prey system according to the linear stability analysis. Then, we show that the discrete model has the flip bifurcation and Turing bifurcation under the critical parameter values. Finally, a series of numerical simulations are carried out in the Turing instability region of the discrete predator-prey model; some new Turing patterns such as striped, bar, and horizontal bar are observed.


---

## Body

## 1. Introduction

Interaction between species and their natural environment is the main characteristic of ecological systems ([1]). Such interaction may occur over a wide range of spatial and temporal scales ([2]). Since the great works of Lotka (in 1925) and Volterra (in 1926) modeling predator-prey relations, interaction has been one of the central themes in mathematical ecology ([3–5]). In general, predator-prey models follow two principles: one is that population dynamics can be decomposed into birth and death processes; the other is the conservation of mass principle, stating that predators can grow only as a function of what they have eaten ([6]).Patterns are ever-present in the chemical and biological worlds; pattern formation is a fundamental problem in the study of far-from-equilibrium phenomena in spatially extended systems. Since Turing ([7]) first introduced his model of pattern formation, reaction-diffusion equations have been a primary means of predicting them. Similarly, structured systems of ordinary differential equations govern the spatiotemporal dynamics of ecological population models. A reaction-diffusion system exhibits diffusion-driven instability or Turing instability if the homogeneous steady state is stable to small perturbations in the absence of diffusion, but it is unstable to small spatial perturbations when diffusion is present. This approach allows us to understand and predict a variety of different phenomena, including the formation of structures that are similar to the patterns we observe in the living world ([8–10]).For the continuous predator-prey system, the mathematical problem is defined by(1)
u
′
t
=
u
t
ε
-
ε
u
t
K
-
α
v
t
+
d
∇
2
u
,
v
′
t
=
v
t
-
γ
+
β
u
t
-
δ
v
t
+
∇
2
v
.The parameters in the model (1) are summarized in the following list:u
(
t
): the quantities at time t of prey
v
(
t
): the quantities at time t of predator
ε: the intrinsic growth rate of the prey
α: the predation rate
γ: the intrinsic mortality of the predator
β: the conversation rate
K: the carrying capacity of the environment with respect to the prey
∇
2: ∂
2
/
∂
x
2
+
∂
2
/
∂
y
2
ε, α, γ, β, d, δ are all positive constants.For simplicity, we rewrite (1) in the following form:(2)
u
′
t
=
u
t
r
1
-
a
11
u
t
-
a
12
v
t
+
d
∇
2
u
,
v
′
t
=
v
t
-
r
2
+
a
21
u
t
-
a
22
v
t
+
∇
2
v
,where r
i
,
a
i
j
i
,
j
=
1,2
>
0. We let(3)
f
u
,
v
=
u
r
1
-
a
11
u
-
a
12
v
,
g
u
,
v
=
v
-
r
2
+
a
21
u
-
a
22
v
.The positive fixed point E
=
u
∗
,
v
∗ of (2) satisfies the system(4)
u
r
1
-
a
11
u
-
a
12
v
=
0
,
v
-
r
2
+
a
21
u
-
a
22
v
=
0
.Then (5)
E
=
r
1
a
22
+
r
2
a
12
a
11
a
22
+
a
12
a
21
,
r
1
a
21
-
r
2
a
11
a
11
a
22
+
a
12
a
21
, where r
1
/
r
2
>
a
11
/
a
21
. From (3), we have(6)
f
u
=
-
a
11
u
,
f
v
=
-
a
12
u
,
g
u
=
a
21
v
,
g
v
=
-
a
22
v
,at the fixed point E
u
∗
,
v
∗, f
u
+
g
v
<
0
,
f
u
g
v
-
f
v
g
u
>
0, but f
u
+
d
g
v
<
0; thus, we can conclude that such a simple continuous predator-prey system can not generate Turing instability. For the discrete time and space predator-prey system, research on the dynamics is not as common. This paper concerns the dynamical behaviors of the discrete predator-prey system.The paper is organized as follows. In Section2, the simple discrete form of (2) and the Turing instability analysis theoretically are studied. In Section 3, we analyze the flip bifurcation of (13). In Section 4, a series of numerical simulations and Lyapunov exponents are given to show the consistence with the theoretical analysis. Conclusions are drawn in Section 5.


## 2. Turing Instability Analysis for the Discrete L-V Predator-Prey Systems

Now, we discuss the Turing instability of the discrete predator-prey system. By Eular’s method, we have the discrete form of system (2):(7)
u
n
+
1
i
,
j
=
r
1
u
n
i
,
j
1
-
u
n
i
,
j
-
a
12
v
n
i
,
j
+
d
∇
2
u
n
i
,
j
,
v
n
+
1
i
,
j
=
r
2
v
n
i
,
j
1
+
a
21
u
n
i
,
j
-
v
n
i
,
j
+
∇
2
v
n
i
,
j
.The Laplacian diffusion parts are defined as (8)
∇
2
u
n
i
j
=
u
n
i
+
1
,
j
+
u
n
i
,
j
+
1
+
u
n
i
-
1
,
j
+
u
n
i
,
j
-
1
-
4
u
n
i
j
,and (9)
∇
2
v
n
i
j
=
v
n
i
+
1
,
j
+
v
n
i
,
j
+
1
+
v
n
i
-
1
,
j
+
v
n
i
,
j
-
1
-
4
v
n
i
j
.
u and v satisfy the periodic boundary conditions(10)
u
n
i
,
0
=
u
n
i
,
m
,
u
n
i
,
1
=
u
n
i
,
m
+
1
,
u
n
0
,
j
=
u
n
m
,
j
,
u
n
1
,
j
=
u
n
m
+
1
,
j
,and(11)
v
n
i
,
0
=
v
n
i
,
m
,
v
n
i
,
1
=
v
n
i
,
m
+
1
,
v
n
0
,
j
=
v
n
m
,
j
,
v
n
1
,
j
=
v
n
m
+
1
,
j
,where i
,
j
∈
1,2
,
⋯
,
m
=
[
1
,
m
] and n and m are a positive integer.In order to find the Turing instability region of the discrete predator-prey model, we first analyze the model with no spatial variation;u and v satisfy(12)
u
n
+
1
i
,
j
=
r
1
u
n
i
,
j
1
-
u
n
i
,
j
-
a
12
v
n
i
,
j
,
v
n
+
1
i
,
j
=
r
2
v
n
i
,
j
1
+
a
21
u
n
i
,
j
-
v
n
i
,
j
.For convenience, we make r
1
=
r
2
=
r, and a
12
=
a
21
=
a, so the model (12) is written as(13)
u
n
+
1
i
,
j
=
r
u
n
i
,
j
1
-
u
n
i
,
j
-
a
v
n
i
,
j
,
v
n
+
1
i
,
j
=
r
v
n
i
,
j
1
+
a
u
n
i
,
j
-
v
n
i
,
j
.There is a nonzero positive fixed point E
1
=
(
u
∗
,
v
∗
) of (13); that is, (14)
μ
∗
=
1
-
a
r
-
1
r
1
+
a
2
,
ν
∗
=
a
+
1
r
-
1
r
1
+
a
2
, where(15)
r
>
1
,
0
<
a
<
1
.The linearization of (13) about E
1 has the Jacobian matrix(16)
J
E
1
=
1
-
r
μ
∗
-
r
a
μ
∗
r
a
ν
∗
1
-
r
ν
∗
.The eigenvalues of (16) are (17)
λ
1
=
a
2
r
-
r
+
2
1
+
a
2
a
n
d
λ
2
=
2
-
r
. The linear stability of (13) is guaranteed if(18)
λ
1
<
1
,
λ
2
<
1
.Then, we consider that the diffusion parts are added. Now, let λ denote an eigenvalue of ∇
2 with the boundary condition (10); that is,(19)
∇
2
u
i
j
+
λ
u
i
j
=
0
.In view of ([11]), λ is in the form of(20)
λ
t
,
s
=
4
s
i
n
2
t
-
1
π
m
+
s
i
n
2
s
-
1
π
m
f
o
r
t
,
s
∈
1
,
m
.The linearized form of (7) is(21)
u
n
+
1
i
j
=
1
-
r
u
∗
u
n
i
j
-
r
a
u
∗
v
n
i
j
+
d
∇
2
u
n
i
j
,
v
n
+
1
i
j
=
r
a
v
∗
u
n
i
j
+
1
-
r
v
∗
v
n
i
j
+
∇
2
v
n
i
j
.Then, we, respectively, take the inner product of (21) with the corresponding eigenfunction X
t
s
i
j of the eigenvalue λ
t
,
s; then(22)
∑
i
,
j
=
1
m
X
t
s
i
j
u
n
+
1
i
j
=
1
-
r
u
∗
∑
i
,
j
=
1
m
X
t
s
i
j
u
n
i
j
-
r
a
u
∗
∑
i
,
j
=
1
m
X
t
s
i
j
v
n
i
j
+
d
∑
i
,
j
=
1
m
X
t
s
i
j
∇
2
u
n
i
j
,
∑
i
,
j
=
1
m
X
t
s
i
j
v
n
+
1
i
j
=
r
a
v
∗
∑
i
,
j
=
1
m
X
t
s
i
j
u
n
i
j
+
1
-
r
v
∗
∑
i
,
j
=
1
m
X
t
s
i
j
v
n
i
j
+
∑
i
,
j
=
1
m
X
t
s
i
j
∇
2
v
n
i
j
. Let (23)
U
n
=
∑
i
,
j
=
1
m
X
t
s
i
j
u
n
i
j and(24)
V
n
=
∑
i
,
j
=
1
m
X
t
s
i
j
v
n
i
j
.Then, we use the periodic boundary conditions (10) and (11), and the Abel transform; thus, we have that(25)
U
n
+
1
=
1
-
r
u
∗
U
n
-
r
a
u
∗
V
n
-
d
λ
t
,
s
U
n
,
V
n
+
1
=
r
a
v
∗
U
n
+
1
-
r
v
∗
V
n
-
λ
t
,
s
V
n
.If (
U
n
,
V
n
) is a solution of the system (25), then u
n
i
j
=
U
n
X
t
s
i
j
,
v
n
i
j
=
V
n
X
t
s
i
j is also clearly a solution of (21) with the periodic boundary conditions (10) and (11); thus, the unstable system (25) produces the problem (7), (10), and (11) is also unstable. The system (25) has the eigenvalue equation(26)
λ
2
+
d
+
1
λ
l
s
+
r
u
∗
+
v
∗
-
2
λ
+
h
λ
t
,
s
=
0
,where(27)
h
λ
t
,
s
=
d
λ
t
,
s
2
-
1
-
r
u
∗
+
d
1
-
r
v
∗
λ
t
,
s
+
1
-
r
u
∗
1
-
r
v
∗
+
r
2
a
2
u
∗
v
∗
. By calculating, the two eigenvalues are(28)
λ
±
t
,
s
,
r
=
1
2
p
t
,
s
,
r
±
1
2
p
t
,
s
,
r
2
-
4
q
t
,
s
,
r
, where(29)
p
t
,
s
,
r
=
-
d
+
1
λ
t
,
s
+
r
u
∗
+
v
∗
-
2
, and(30)
q
t
,
s
,
r
=
d
λ
t
,
s
2
-
1
-
r
u
∗
+
d
1
-
r
v
∗
λ
t
,
s
+
1
-
r
u
∗
1
-
r
v
∗
+
r
2
a
2
u
∗
v
∗
. On the basis of the two eigenvalues, we define(31)
Z
t
,
s
,
r
=
m
a
x
λ
+
t
,
s
,
λ
-
t
,
s
,
Z
m
r
=
m
a
x
t
=
1
,
s
=
1
m
Z
t
,
s
,
r
t
,
s
≠
1,1
,
Z
m
(
r
) represents the maximal value of absolute modulus of both eigenvalues λ
+ and λ
-. When Z
m
(
r
)
>
1, Turing instability occurs; when Z
m
(
r
)
<
1, the discrete system stabilizes at the homogeneous states. Thus, the threshold condition for the occurrence of Turing bifurcation requires Z
m
(
r
)
=
1. From this criterion, the critical value r
′ for Turing bifurcation can be described by the following cases.(1)p
(
t
,
s
,
r
)
2
>
4
q
(
t
,
s
,
r
) establishes a small neighborhood of r
=
r
′, and the critical value r
′ satisfies(32)
m
a
x
t
=
1
,
s
=
1
m
p
t
,
s
,
r
′
-
q
t
,
s
,
r
′
=
1
.(2) p
(
t
,
s
,
r
)
2
≤
4
q
(
t
,
s
,
r
) establishes a small neighborhood of r
=
r
′ and the critical value r
′ satisfies(33)
m
a
x
t
=
1
,
s
=
1
m
q
t
,
s
,
r
′
=
1
. Thus, the conditions of Turing instability of (7) are(34)
0
<
a
<
1
,
r
>
1
,
λ
1
<
1
,
λ
2
<
1
,
Z
m
r
>
1
.


## 3. Bifurcations and Center Manifolds

It is important to discuss the bifurcations and the center manifolds for the applications; for example, see ([12–16]). In this section, we study the flip bifurcation of (13) at the positive steady state E
1
=
μ
∗
,
ν
∗. When the flip bifurcation occurs, (
μ
∗
,
ν
∗
) loses its stability, and the discrete system (13) switches to a new behavior with period-2. At the flip bifurcation point, (
μ
∗
,
ν
∗
) is neither stable nor unstable. In this critical case, one of the two eigenvalues of J
(
(
μ
∗
,
ν
∗
) satisfies λ
1
=
-
1, λ
2
=
2
-
r
∗ when r
∗
=
a
2
+
3
/
1
-
a
2. Regarding r as the dependent bifurcation parameter, we have the following theorem.Theorem 1.
If the condition (15) is established, model (13) undergoes a flip bifurcation if α
1
≠
0
,
α
2
≠
0, and r
=
r
∗; furthermore, if α
2
>
0 is satisfied, then the bifurcated period-2 points are stable; if α
2
<
0, the bifurcated period-2 points are unstable.Proof.
Let(35)
ζ
n
=
μ
n
-
μ
∗
,
η
n
=
ν
n
-
ν
∗
,
δ
n
=
r
-
r
∗
,and parameter δ
n is a new and independent variable; the system (13) becomes(36)
ζ
n
+
1
=
a
11
ζ
n
+
a
12
η
n
+
a
13
δ
n
+
a
14
ζ
n
2
+
a
15
ζ
n
η
n
+
b
11
ζ
n
δ
n
+
b
12
η
n
δ
n
+
O
ζ
n
+
η
n
+
δ
n
3
,
η
n
+
1
=
a
21
ζ
n
+
a
22
η
n
+
a
23
δ
n
+
a
24
η
n
2
+
a
25
ζ
n
η
n
+
b
21
ζ
n
δ
n
+
b
22
η
n
δ
n
+
O
ζ
n
+
η
n
+
δ
n
3
,
δ
n
+
1
=
δ
n
,where(37)
a
11
=
1
-
r
∗
u
∗
,
a
12
=
-
a
r
∗
u
∗
,
a
13
=
-
u
∗
r
∗
,
a
14
=
-
r
∗
,
a
15
=
-
a
r
∗
,
a
21
=
a
r
∗
v
∗
,
a
22
=
1
-
r
∗
v
∗
,
a
23
=
v
∗
r
∗
,
a
24
=
-
r
∗
,
a
25
=
a
r
∗and(38)
b
11
=
1
-
2
u
∗
-
a
v
∗
,
b
12
=
-
a
u
∗
,
b
21
=
a
v
∗
,
b
22
=
1
+
a
u
∗
-
2
v
∗
.Let(39)
T
=
a
12
a
12
a
12
+
a
13
-
1
-
a
11
λ
2
-
a
11
1
-
a
11
0
0
1
-
a
11
; then,(40)
T
-
1
=
1
T
λ
2
-
a
11
1
-
a
11
-
a
12
1
-
a
11
a
12
1
-
a
11
-
a
12
+
a
13
λ
2
-
a
11
1
+
a
11
1
-
a
11
a
12
(
1
-
a
11
-
a
12
1
-
a
11
+
1
+
a
11
a
12
+
a
13
0
0
a
12
λ
2
-
a
11
+
1
+
a
11
a
12
. The transformation(41)
ζ
n
η
n
δ
n
=
T
u
n
v
n
δ
1
, changes (36) into(42)
u
n
+
1
v
n
+
1
δ
1
=
-
1
0
0
0
λ
2
0
0
0
1
u
n
v
n
δ
1
+
1
T
f
x
n
,
y
n
,
δ
n
g
x
n
,
y
n
,
δ
n
0
, where(43)
f
x
n
,
y
n
,
δ
n
=
λ
2
-
a
11
1
-
a
11
a
14
ζ
n
2
-
a
12
1
-
a
11
a
24
η
n
2
+
a
15
λ
2
-
a
11
1
-
a
11
-
a
25
a
12
1
-
a
11
ζ
n
η
n
+
b
11
λ
2
-
a
11
1
-
a
11
-
b
21
a
12
1
-
a
11
ζ
n
δ
n
+
b
12
λ
2
-
a
11
1
-
a
11
-
b
22
a
12
1
-
a
11
η
n
δ
n
+
O
x
n
+
+
y
n
+
δ
n
3
,
g
x
n
,
y
n
,
δ
n
=
a
14
1
-
a
11
2
ζ
n
2
+
a
24
a
12
1
-
a
11
η
n
2
+
a
15
1
-
a
11
2
+
a
25
a
12
1
-
a
11
ζ
n
η
n
+
b
11
1
-
a
11
2
+
b
21
a
12
1
-
a
11
ζ
n
δ
n
+
b
12
1
-
a
11
2
+
b
22
a
12
1
-
a
11
η
n
δ
n
+
O
x
n
+
+
y
n
+
δ
n
3
,
ζ
n
=
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
,
η
n
=
(
-
1
-
a
11
u
n
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
,
δ
n
=
1
-
a
11
δ
1
. Then,(44)
f
x
n
,
y
n
,
δ
n
=
λ
2
-
a
11
1
-
a
11
a
14
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
2
-
a
12
1
-
a
11
a
24
-
1
-
a
11
u
n
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
2
+
a
15
λ
2
-
a
11
1
-
a
11
-
a
25
a
12
1
-
a
11
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
-
1
-
a
11
u
n
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
+
b
11
λ
2
-
a
11
1
-
a
11
-
b
21
a
12
1
-
a
11
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
1
-
a
11
δ
1
+
b
12
λ
2
-
a
11
1
-
a
11
-
b
22
a
12
1
-
a
12
-
1
-
a
11
u
n
)
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
1
-
a
11
δ
1
,
g
x
n
,
y
n
,
δ
n
=
a
14
1
-
a
11
2
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
2
+
a
24
a
12
1
-
a
11
-
1
-
a
11
u
n
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
2
+
a
15
1
-
a
11
2
+
a
25
a
12
1
-
a
11
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
-
1
-
a
11
u
n
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
+
b
11
1
-
a
11
2
+
b
21
a
12
1
-
a
11
a
12
u
n
+
a
12
v
n
+
a
12
+
a
13
δ
1
1
-
a
11
δ
1
+
b
12
1
-
a
11
2
+
b
22
a
12
1
-
a
11
-
1
-
a
11
u
n
+
λ
2
-
a
11
v
n
+
1
-
a
11
δ
1
1
-
a
11
δ
1
.Applying the center manifold theorem, there exists a center manifold W
0
(
0,0
,
0
) of the model in a small neighborhood of δ
1
=
0, which can be represented as (45)
W
0
0,0
,
0
=
u
n
,
v
n
,
δ
1
∈
R
3
∣
v
n
=
h
u
n
,
δ
1
=
c
1
u
n
2
+
c
2
u
n
δ
1
+
c
3
δ
1
2
+
O
u
n
+
δ
1
3
)with(46)
v
n
+
1
=
λ
2
h
u
n
,
δ
1
+
g
x
n
,
y
n
,
δ
1
=
h
u
n
+
1
,
δ
1
=
h
-
u
n
+
f
x
n
,
y
n
,
δ
n
,
δ
1
.By calculating, we can obtain (47)
c
1
=
a
14
1
-
a
11
2
a
12
2
+
a
24
a
12
1
-
a
11
-
1
-
a
11
2
+
a
15
1
-
a
11
2
+
a
25
a
12
1
-
a
11
a
12
-
1
-
a
11
T
1
-
λ
2
,
c
2
=
2
a
14
1
-
a
11
2
a
12
a
12
+
a
13
+
2
a
24
a
12
1
-
a
11
2
-
1
-
a
11
T
1
-
λ
2
+
a
15
1
-
a
11
2
+
a
25
a
12
1
-
a
11
a
12
1
-
a
11
+
a
12
+
a
13
-
1
-
a
11
T
1
-
λ
2
+
b
11
1
-
a
11
2
+
b
21
a
12
1
-
a
11
a
12
1
-
a
11
T
1
-
λ
2
+
b
12
1
-
a
11
2
+
b
22
a
12
1
-
a
11
-
1
-
a
11
1
-
a
11
T
1
-
λ
2
,
c
3
=
0
. The model (36) restricted to the center manifold can be written as(48)
H
:
u
n
+
1
=
-
u
n
+
h
1
u
n
2
+
h
2
u
n
δ
1
+
h
3
u
n
2
δ
1
+
h
4
u
n
δ
1
2
+
h
5
u
n
3
+
O
u
n
+
δ
1
4
,
h
1
=
1
T
λ
2
-
a
11
1
-
a
11
a
14
a
12
2
-
a
12
a
24
1
-
a
11
-
1
-
a
11
2
+
a
15
λ
2
-
a
11
1
-
a
11
-
a
25
a
12
1
-
a
11
a
12
-
1
-
a
11
,
h
2
=
1
T
2
λ
2
-
a
11
1
-
a
11
a
14
a
12
a
12
+
a
13
-
2
a
12
1
-
a
11
a
24
1
-
a
11
-
1
-
a
11
+
a
15
1
-
a
11
λ
2
-
a
11
-
a
25
a
12
1
-
a
11
1
-
a
11
a
12
+
a
12
+
a
13
-
1
-
a
11
+
b
11
1
-
a
11
λ
2
-
a
11
-
b
21
a
12
1
-
a
11
a
12
1
-
a
11
,
h
3
=
1
T
c
2
2
λ
2
-
a
11
1
-
a
11
a
14
a
12
2
-
2
a
12
1
-
a
11
a
24
-
1
-
a
11
λ
2
-
a
11
+
a
15
1
-
a
11
λ
2
-
a
11
-
a
25
a
12
1
-
a
11
a
12
(
λ
2
-
a
11
+
a
12
-
1
-
a
11
+
1
T
c
1
2
λ
2
-
a
11
1
-
a
11
a
14
a
12
a
12
+
a
13
+
2
1
-
a
11
a
24
a
12
λ
2
-
a
11
1
-
a
11
+
a
15
λ
2
-
a
11
1
-
a
11
-
a
25
a
12
1
-
a
11
a
12
+
a
13
λ
2
-
a
11
+
a
12
1
-
a
11
+
b
11
1
-
a
11
λ
2
-
a
11
-
b
21
a
12
1
-
a
11
a
12
1
-
a
11
+
b
11
1
-
a
11
λ
2
-
a
11
-
b
22
a
12
1
-
a
11
λ
2
-
a
11
1
-
a
11
,
h
4
=
1
T
c
2
2
λ
2
-
a
11
1
-
a
11
a
14
a
12
a
12
+
a
13
+
2
1
-
a
11
a
24
a
12
λ
2
-
a
11
1
-
a
11
+
a
15
λ
2
-
a
11
1
-
a
11
-
a
25
a
12
1
-
a
11
a
12
+
a
13
λ
2
-
a
11
+
a
12
1
-
a
11
+
b
11
1
-
a
11
λ
2
-
a
11
-
b
21
a
12
1
-
a
11
a
12
1
-
a
11
+
b
11
1
-
a
11
λ
2
-
a
11
-
b
22
a
12
1
-
a
11
λ
2
-
a
11
1
-
a
11
,
h
5
=
1
T
c
1
2
λ
2
-
a
11
1
-
a
11
a
14
a
12
2
-
2
a
12
1
-
a
11
a
24
-
1
-
a
11
λ
2
-
a
11
+
a
15
1
-
a
11
λ
2
-
a
11
-
a
25
a
12
1
-
a
11
a
12
λ
2
-
a
11
+
a
12
-
1
-
a
11
. As stated by the flip bifurcation theorem in ([17]), the emergence of flip bifurcation for map (3) requires (49)
α
1
=
∂
2
H
∂
u
∂
δ
1
+
1
2
∂
H
∂
δ
1
∂
2
H
∂
u
2
0,0
=
h
2
≠
0
,
α
2
=
1
6
∂
3
H
∂
u
3
+
1
2
∂
2
H
∂
u
2
2
0,0
=
h
5
+
h
1
2
≠
0
.


## 4. Numerical Simulation

To illustrate the analytical results in the above sections and find new dynamics with different parameters, in this section, we provide some numerical evidence for the qualitative dynamic behavior of model (7).


### 4.1. Simulations about the Flip Bifurcation

In the following, we display the bifurcation diagrams and the Lyapunov exponents. Now, the fix point is(
u
∗
,
v
∗
)
=
(
0.37974684,0.88607595
), a
=
0.001, and r
∗ is considered a parameter with the range (
1.1
-
5
). From Figure 1(a) we see that equilibrium is stable for r
=
2.9, at this moment r
<
r
∗
=
3.000004, and the eigenvalues λ
1
=
0.9 and λ
2
=
0.8999, lose stability when r
=
r
∗. Furthermore, when r
>
r
∗, there is a period-doubling bifurcation. We also observe that there is a cascade of period doubling; moreover, a chaotic set emerges with increasing of r. At the same time, we show the maximum Lyapunov exponents, as shown in Figure 1(b); we find that the periodic windows are repeated. Quantitatively, we can determine the chaotic region; we find that when r
=
3.75 the maximum Lyapunov exponents are greater than zero. Model (7) may experience chaotic oscillating states.Figure 1
Flip bifurcation diagram of model (7) with r
∈
[
1.1
;
5
]; maximum Lyapunov exponents corresponding to flip bifurcation diagram.
(a)
Bifurcation diagram
(b)
The maximum Lyapunov exponent


### 4.2. Simulations of Related Spatial Turing Patterns

From the analysis of the second chapter, we know that, under conditions (15) and (18), the point (
u
∗
,
v
∗
) remains stable, and no Turing instability is induced. At the same time, under (15), (18), and Z
m
(
r
)
>
1, the point (
u
∗
,
v
∗
) becomes unstable from stable under diffusive effects, so Turing instability occurs. Briefly, we only display several patterns of the u and take t
=
5000. We fix a
=
0.001 and d
=
0.2, as shown in Figure 2. In the diagram of Figure 2, the abscissa is r, and the ordinate is the maximum eigenvalue of (25). When r
>
r
′, the critical value of Turing instability r
′
≈
1.4, the maximum eigenvalue of (25) is >
1; under conditions (15) and (18), we know r
<
3. At this time, some Turing patterns form. We find 0
<
r
′
<
r
∗; if 0
<
r
<
r
′, the stable homogeneous steady state keeps its stability; flip and Turing bifurcation cannot occur in this region. When r
′
<
r
<
3, the homogeneous steady state is not stable because of the diffusions. Flip bifurcation cannot appear this moment; only pure Turing instability appears in the homogeneous steady state. And when r
>
r
∗, both flip and Turing bifurcations arise. In the following, we give some new Turing patterns in the Turing instability region. The model (7) has three parameters: r, a, and d. To study the effects that parameters have on pattern formation, we assume that only one parameter is remaining fixed, and the others are changing. Then Figure 3 is obtained. The X-axis is time, and the Y-axis is the number of u in Figure 3. Firstly, in Figures 3(a), 3(b), and 3(c), we fix d
=
0.2207, and the other two parameters a and r are smaller. Comparing the pattern structure reveals the transition from diagonal striped to two horizontal stripes, and, finally, it becomes diagonal striped Turing patterns. In Figures 3(d), 3(e), and 3(f), when we fix d
=
0.2207, the other two parameters a and r are larger. By comparing the pattern structure, the two vertical striped patterns are found, as the parameters increase, and the two horizontal stripes arise; finally it becomes diagonal stripes. In the end, we fix r
=
1.27; the other two parameters a and d are increasing, and the two horizontal stripes first turn into diagonal stripes and, finally, turn into the vertical striped pattern; see Figures 3(g), 3(h), and 3(i).Figure 2
Turing bifurcation with critical pointr’ ≈ 1.4.Figure 3
Turing patterns.
(a)
d
=
0.2207, r
=
1.36, a
=
0.365
(b)
d
=
0.2207, r
=
1.35, a
=
0.246
(c)
d
=
0.2207, r
=
1.34, a
=
0.177
(d)
d
=
0.2207, r
=
1.29, a
=
0.332
(e)
d
=
0.2207, r
=
1.32, a
=
0.389
(f)
d
=
0.2207, r
=
1.35, a
=
0.446
(g)
r
=
1.27,, a
=
0.144, d
=
0.217
(h)
r
=
1.27, a
=
0.244, d
=
0.22
(i)
r
=
1.27, a
=
0.344, d
=
0.2234


## 4.1. Simulations about the Flip Bifurcation

In the following, we display the bifurcation diagrams and the Lyapunov exponents. Now, the fix point is(
u
∗
,
v
∗
)
=
(
0.37974684,0.88607595
), a
=
0.001, and r
∗ is considered a parameter with the range (
1.1
-
5
). From Figure 1(a) we see that equilibrium is stable for r
=
2.9, at this moment r
<
r
∗
=
3.000004, and the eigenvalues λ
1
=
0.9 and λ
2
=
0.8999, lose stability when r
=
r
∗. Furthermore, when r
>
r
∗, there is a period-doubling bifurcation. We also observe that there is a cascade of period doubling; moreover, a chaotic set emerges with increasing of r. At the same time, we show the maximum Lyapunov exponents, as shown in Figure 1(b); we find that the periodic windows are repeated. Quantitatively, we can determine the chaotic region; we find that when r
=
3.75 the maximum Lyapunov exponents are greater than zero. Model (7) may experience chaotic oscillating states.Figure 1
Flip bifurcation diagram of model (7) with r
∈
[
1.1
;
5
]; maximum Lyapunov exponents corresponding to flip bifurcation diagram.
(a)
Bifurcation diagram
(b)
The maximum Lyapunov exponent


## 4.2. Simulations of Related Spatial Turing Patterns

From the analysis of the second chapter, we know that, under conditions (15) and (18), the point (
u
∗
,
v
∗
) remains stable, and no Turing instability is induced. At the same time, under (15), (18), and Z
m
(
r
)
>
1, the point (
u
∗
,
v
∗
) becomes unstable from stable under diffusive effects, so Turing instability occurs. Briefly, we only display several patterns of the u and take t
=
5000. We fix a
=
0.001 and d
=
0.2, as shown in Figure 2. In the diagram of Figure 2, the abscissa is r, and the ordinate is the maximum eigenvalue of (25). When r
>
r
′, the critical value of Turing instability r
′
≈
1.4, the maximum eigenvalue of (25) is >
1; under conditions (15) and (18), we know r
<
3. At this time, some Turing patterns form. We find 0
<
r
′
<
r
∗; if 0
<
r
<
r
′, the stable homogeneous steady state keeps its stability; flip and Turing bifurcation cannot occur in this region. When r
′
<
r
<
3, the homogeneous steady state is not stable because of the diffusions. Flip bifurcation cannot appear this moment; only pure Turing instability appears in the homogeneous steady state. And when r
>
r
∗, both flip and Turing bifurcations arise. In the following, we give some new Turing patterns in the Turing instability region. The model (7) has three parameters: r, a, and d. To study the effects that parameters have on pattern formation, we assume that only one parameter is remaining fixed, and the others are changing. Then Figure 3 is obtained. The X-axis is time, and the Y-axis is the number of u in Figure 3. Firstly, in Figures 3(a), 3(b), and 3(c), we fix d
=
0.2207, and the other two parameters a and r are smaller. Comparing the pattern structure reveals the transition from diagonal striped to two horizontal stripes, and, finally, it becomes diagonal striped Turing patterns. In Figures 3(d), 3(e), and 3(f), when we fix d
=
0.2207, the other two parameters a and r are larger. By comparing the pattern structure, the two vertical striped patterns are found, as the parameters increase, and the two horizontal stripes arise; finally it becomes diagonal stripes. In the end, we fix r
=
1.27; the other two parameters a and d are increasing, and the two horizontal stripes first turn into diagonal stripes and, finally, turn into the vertical striped pattern; see Figures 3(g), 3(h), and 3(i).Figure 2
Turing bifurcation with critical pointr’ ≈ 1.4.Figure 3
Turing patterns.
(a)
d
=
0.2207, r
=
1.36, a
=
0.365
(b)
d
=
0.2207, r
=
1.35, a
=
0.246
(c)
d
=
0.2207, r
=
1.34, a
=
0.177
(d)
d
=
0.2207, r
=
1.29, a
=
0.332
(e)
d
=
0.2207, r
=
1.32, a
=
0.389
(f)
d
=
0.2207, r
=
1.35, a
=
0.446
(g)
r
=
1.27,, a
=
0.144, d
=
0.217
(h)
r
=
1.27, a
=
0.244, d
=
0.22
(i)
r
=
1.27, a
=
0.344, d
=
0.2234


## 5. Conclusions

This study demonstrates that space and time discrete predator-prey system can produce Turing instability, whereas continuous ones cannot. Through linear stability analysis, the conditions for producing Turing instability of the discrete predator-prey system are obtained. A series of numerical simulations are given and many new and interesting striped patterns are observed from the simulation. In addition, we show that the discrete model has the flip bifurcation and Turing bifurcation under the critical parameter values by the bifurcation theory and center manifold theorem. Finally, for the discrete time and space predator-prey model, a new instability mechanism is found, namely, flip-Turing instability, which is the basis of chaos. The bifurcation, chaos, and pattern formation provide us with a new and better understanding of dynamical complexity of the space and time discrete predator-prey system.



---

*Source: 9592878-2019-04-01.xml*