# The Dynamics of a Diffusive Nutrient-Algae Model Based upon the Sanyang Wetland

Authors: Yi Wang; Min Zhao; Chuanjun Dai; Yunli Deng
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The stability and spatiotemporal dynamics of a diffusive nutrient-algae model are investigated mathematically and numerically. Mathematical theoretical studies have considered the positivity and boundedness of the solution and the existence, local stability, and global stability of equilibria. Turing instability has also been studied. Furthermore, a series of numerical simulations was performed and a complex Turing pattern found. These results indicate that the nutrient input rate has an important influence on the density and spatial distribution of algae populations. This may help us to obtain a better understanding of the interactions of nutrient and algae and to investigate plankton dynamics in aquatic ecosystems.


---

## Body

## 1. Introduction

Diffusive processes are often used to represent the formation of spatial patterns in biological systems [1]. Pattern formation in nonlinear complex systems is one of the central problems of the natural, social, and technological sciences. The occurrence of multiple steady states and transitions from one to another after critical fluctuations, the phenomena of excitability, oscillations, and waves, and the emergence of macroscopic order from microscopic interactions in various nonlinear nonequilibrium systems in nature and society have been the subject of many theoretical and experimental studies [2, 3]. Plankton plays an important role in ocean ecology and climate because of their participation in the global carbon cycle at the base of the food chain [4]. Now, using the reaction-diffusion equation and patterns to study the spatiotemporal dynamics of the planktonic ecosystem has aroused the interest of many researchers.In recent years, the relationships among nutrients, phytoplankton, zooplankton, and fish have become the focus of researchers [5–8], including the control of phytoplankton [9, 10], and dynamic analysis of the relationships [11, 12]. Numerous theoretical ecologists have built many models to reveal the inner relationships among these populations and to investigate their dynamics. Drago et al. [13] presented a three-dimensional numerical model to analyze the dispersion of suspended solids and conservative pollutants released into ambient water and their effect on trophic behavior. Luo [14] derived and analyzed a mathematical model for interactions between phytoplankton and zooplankton in a periodic environment, where the growth rate and the intrinsic carrying capacity of phytoplankton are changing with respect to time and nutrient concentration. James et al. [15] built a model of the evolution of phytoplankton, zooplankton, and fish to investigate the mechanisms that trigger plankton blooms. Alvarez-Vázquez et al. [16] presented a mathematical model involving nutrients, phytoplankton, zooplankton, organic detritus, and dissolved oxygen to simulate eutrophication processes in aquatic media.The Sanyang wetland is located in Wenzhou, in a subtropical region. It covers an area of 13 square kilometers and contains many streams. The ratio of water area to land in the Sanyang wetland is as high as 1.1 : 1. With rapid economic development in Wenzhou, most rivers and wetlands are in a state of eutrophication. Because of eutrophication, nuisance cyanobacterial blooms have taken place several times in the Sanyang wetland in recent years.Recently, more and more researchers have investigated the dynamics of the planktonic ecosystem by means of spatial patterns [17–24]. Self-organizing spatial patterns appeal to the theoretical biologist because of the pioneering work of Turing [25]. Many self-organizing spatial patterns have been found in their research, including patchiness, bands, and spiral waves. Wang et al. [26] used a reaction-diffusion-advection model of algae and mussels to demonstrate that young mussel beds on soft sediments can display large-scale regular spatial patterns and banded patterns. Their work is significant to study spatial patterns. Serizawa et al. [27] presented a minimal nutrient-phytoplankton model that can exhibit various types of spatial patterns, including patchiness. Although the model is simple, their work is still of interest. Liu et al. [28] studied a spatially extended nutrient-phytoplankton-zooplankton-fish reaction-diffusion system and found spiral waves and spatial chaos patterns. Their work is important to plankton ecological system. van de Koppel et al. [29] proposed a simple spatial model of algae and mussels and found simulated spatial patterns which are consistent with previous observations. Their works vigorously promote the study of emergent spatial patterning at larger spatial scales.The nutrients in the Sanyang wetland comes mainly from sewage output by nearby residents and wastewater discharges from industries and businesses. Therefore, the input rate of nutrients flowing into the Sanyang wetland can be assumed constant. Many researchers have usedb
x to describe natural nutrient removal [27, 30, 31]. Therefore, the term b
u has been used here to describe the natural removal of nutrients (mainly in the form of sedimentation) in the Sanyang wetland because of various complex long-term factors. However, this term has another important meaning in the mathematical model. Without this factor, the nutrient concentration will increase without limit when the algae are completely extinct. Wenzhou belongs to a subtropical monsoon climate zone and is perennially windy. Therefore, the effects of wind and waves on nutrients in the Sanyang wetland must be considered. Because of wind and wave action, a positive nutrient input is released into the water from bottom sediments. Ecologists usually approximate this input by the following sigmoid function [32]:(1)
f
(
x
)
=
x
2
1
+
x
2
.The HollingI functional response function is suitable for algae, cells, and lower organisms. Many studies have shown that algae density in the Sanyang wetland is positively correlated with nutrient concentration within certain limits. Hence, in this paper, the Holling I response function has been used to describe absorption of nutrients by algae. In the Sanyang wetland, nutrients are limited, and density restrictions exist in the algae population. Therefore, it is suitable to use a logistic growth function to describe algae population growth in the Sanyang wetland.Based on this discussion, the following dynamic reaction-diffusion model has been proposed for nutrients and algae in the Sanyang wetland:(2)
∂
u
∂
t
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
+
d
1
Δ
u
,
∂
v
∂
t
=
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
+
d
2
Δ
v
,where u and v denote the nutrient concentration and the algae density. Research has shown that phosphorus is limited and nitrogen is abundant in this wetland. Therefore, the main factor affecting algae population growth is phosphorus, and, in this paper, the nutrient has been assumed to be phosphorus. Let I be the input rate of nutrients flowing into the water, a the maximum growth rate of the algae population, e the efficiency of conversion, r the intrinsic growth rate, k the carrying capacity of the algae population, m the death rate of the algae population, and d
1 and d
2 the diffusion coefficients of nutrients and algae, respectively. Δ
=
∂
2
/
∂
x
2
+
∂
2
/
∂
y
2 is the usual Laplacian operator in two-dimensional space, and I, b, a, e, r, k, and m are positive constants.Model (2) will be studied under the following nonzero initial conditions,(3)
u
(
x
,
y
,
0
)
>
0
,
v
(
x
,
y
,
0
)
>
0
,
(
x
,
y
)
∈
Ω
=
0
,
L
x
×
0
,
L
y
,and the following zero-flux boundary conditions,(4)
∂
u
∂
n
=
∂
v
∂
n
=
0
,
(
x
,
y
)
∈
∂
Ω
,where L
x and L
y are the size of the system in the x- and y-directions and n is the outward unit normal vector of the boundary ∂
Ω, which is assumed to be smooth.For the reaction-diffusion nutrient-algae system (2), the reduced system is an ordinary differential equation of the form(5)
d
u
d
t
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
,
d
v
d
t
=
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
.Next, the ODE system (5) will be analyzed.


## 2. Dynamics Analysis of the Reduced ODE System

### 2.1. Positivity and Boundedness of the Solutions

Theorem 1.
Suppose thata
e
L
+
r
-
m
>
0
(
L
=
max
⁡
u
(
0
)
,
(
I
+
1
)
/
b
) holds; then all the solutions of system (5) with initial conditions are positive and bounded for all t
≥
0.Proof.
From the first equation of system (5), (6)
d
u
d
t
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
≥
-
b
u
+
u
2
1
+
u
2
-
a
u
v
.It follows that d
u
/
d
t
≥
(
-
b
+
u
/
(
1
+
u
2
)
-
a
v
)
u.
Hence,u
(
t
)
≥
u
(
0
)
exp
⁡
∫
0
t
[
-
b
+
u
(
s
)
/
(
1
+
(
u
s
)
2
)
-
a
v
(
s
)
]
d
s
>
0.
From the second equation of system (5), it can be determined that(7)
v
(
t
)
=
v
(
0
)
exp
⁡
∫
0
t
[
r
(
1
-
v
(
s
)
k
)
+
a
e
u
(
s
)
-
m
]
d
s
>
0
,
d
u
d
t
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
≤
I
-
b
u
+
u
2
1
+
u
2
≤
I
+
1
-
b
u
,leading to d
u
/
d
t
≤
I
+
1
-
b
u.
Then,(8)
u
(
t
)
≤
I
+
1
b
(
u
(
0
)
-
I
+
1
b
)
e
-
b
t
≤
max
⁡
u
(
0
)
,
I
+
1
b
=
L
,
d
v
d
t
=
v
[
r
(
1
-
v
k
)
+
a
e
u
-
m
]
≤
v
[
r
(
1
-
v
k
)
+
a
e
L
-
m
]
=
v
(
a
e
L
+
r
-
m
-
r
k
v
)
.If a
e
L
+
r
-
m
>
0, then, for ε
>
0, there exists T
(
>
0
) such that, for any t
>
T, (9)
v
(
t
)
≤
k
(
a
e
L
+
r
-
m
)
r
+
ε
.Moreover, v
(
t
)
=
v
(
0
)
exp
⁡
∫
0
t
[
r
(
1
-
v
(
s
)
/
k
)
+
a
e
u
(
s
)
-
m
]
d
s is continuous when t
∈
[
0
,
T
].
This means that there existsM
(
>
0
) such that v
(
t
)
≤
M
t
∈
[
0
,
T
].
Therefore,v
(
t
)
≤
max
⁡
k
(
a
e
L
+
r
-
m
)
/
r
+
ε
,
M when t
≥
0, and the proof is complete.


### 2.2. Existence of Equilibria

Theorem 2.
System (5) has boundary equilibrium E
1
=
(
u
-
,
0
) (extinction of algae population).Proof.
IfE
1
=
(
u
-
,
0
) is the boundary equilibrium of (5), then(10)
I
-
b
u
-
+
u
-
2
1
+
u
-
2
=
0
that
is
,
b
u
-
3
-
(
I
+
1
)
u
-
2
+
b
u
-
-
I
=
0
.Assume the following function:(11)
f
(
u
)
=
b
u
3
-
(
I
+
1
)
u
2
+
b
u
-
I
.Obviously,(12)
f
(
0
)
=
-
I
,
lim
u
→
+
∞
f
(
u
)
=
+
∞
.
It is easy to establish that the curve of functionf
(
u
) intersects the u-axis, and, therefore, the existence of the boundary equilibrium is guaranteed, which completes the proof.Theorem 3.
System (5) has a unique interior equilibrium E
∗
=
(
u
∗
,
v
∗
) (coexistence of nutrients and algae) if a
2
k
e
-
r
I
-
r
>
0 and a
k
(
r
-
m
)
+
b
r
>
0.Proof.
IfE
∗
=
(
u
∗
,
v
∗
) is the interior equilibrium of (5), then(13)
I
-
b
u
∗
+
u
∗
2
1
+
u
∗
2
-
a
u
∗
v
∗
=
0
,
r
v
∗
(
1
-
v
∗
k
)
+
a
e
u
∗
v
∗
-
m
v
∗
=
0
.Consequently, u
∗ is the positive root of the fourth-degree equation:(14)
a
2
k
e
u
∗
4
+
a
k
(
r
-
m
)
+
b
r
u
∗
3
+
(
a
2
k
e
-
r
I
-
r
)
u
∗
2
+
a
k
(
r
-
m
)
+
b
r
u
∗
-
r
I
=
0
,
v
∗
=
k
r
(
a
e
u
∗
+
r
-
m
)
.Consider the following function:(15)
h
x
=
a
2
k
e
x
4
+
a
k
(
r
-
m
)
+
b
r
x
3
+
(
a
2
k
e
-
r
I
-
r
)
x
2
+
a
k
(
r
-
m
)
+
b
r
x
-
r
I
.It is easy to establish that(16)
h
′
x
=
4
a
2
k
e
x
3
+
3
a
k
(
r
-
m
)
+
b
r
x
2
+
2
(
a
2
k
e
-
r
I
-
r
)
x
+
a
k
(
r
-
m
)
+
b
r
=
[
4
a
2
k
e
x
2
+
2
(
a
2
k
e
-
r
I
-
r
)
]
x
+
a
k
r
-
m
+
b
r
3
x
2
+
1
.It is also easy to determine that h
′
(
x
)
>
0
(
x
>
0
) if a
2
k
e
-
r
I
-
r
>
0 and a
k
(
r
-
m
)
+
b
r
>
0.
Obviously,(17)
h
(
0
)
=
-
r
I
,
lim
x
→
+
∞
h
(
x
)
=
+
∞
.
It can easily be shown that the curve of functionh
(
x
) intersects the x-axis only once. Hence, the existence and uniqueness of the interior equilibrium is guaranteed, and this completes the proof.


### 2.3. Local Stability of Equilibria

Theorem 4.
The boundary equilibriumE
1
=
(
u
-
,
0
) of system (5) is locally asymptotically stable if (18)
-
b
+
2
u
-
1
+
u
-
2
2
<
0
,
a
e
u
-
+
r
-
m
<
0
.The Jacobian matrix of system (5) is (19)
J
=
-
b
-
a
v
+
2
u
1
+
u
2
2
-
a
u
a
e
v
r
-
2
r
k
v
+
a
e
u
-
m
.The Jacobian matrix of system (5) at E
1
=
(
u
-
,
0
) is(20)
J
(
E
1
)
=
-
b
+
2
u
-
1
+
u
-
2
2
-
a
u
-
0
a
e
u
-
+
r
-
m
.It is clear that the Jacobian matrix J
(
E
1
) has two eigenvalues: (21)
λ
1
=
-
b
+
2
u
-
1
+
u
-
2
2
,
λ
2
=
a
e
u
-
+
r
-
m
.Hence, the equilibrium E
1
=
(
u
-
,
0
) is locally asymptotically stable if (22)
-
b
+
2
u
-
1
+
u
-
2
2
<
0
,
a
e
u
-
+
r
-
m
<
0
.Theorem 5.
The interior equilibriumE
∗
=
(
u
∗
,
v
∗
) of system (5) is locally asymptotically stable if -
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
<
0.Proof.
The Jacobian matrix of system (5) at E
∗
=
(
u
∗
,
v
∗
) is(23)
J
(
E
∗
)
=
J
11
J
12
J
21
J
22
=
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
-
a
u
∗
a
e
v
∗
-
r
k
v
∗
.Moreover, (24)
det
⁡
(
λ
E
-
J
(
E
∗
)
)
=
λ
2
-
tr
⁡
(
J
(
E
∗
)
)
λ
+
det
⁡
(
J
(
E
∗
)
)
,where tr
⁡
(
J
E
∗
)
=
a
11
-
(
r
/
k
)
v
∗, det
⁡
(
J
E
∗
)
=
(
a
2
e
u
∗
-
a
11
(
r
/
k
)
)
v
∗, and a
11
=
-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2.
It can be easily determined thattr
⁡
(
J
(
E
∗
)
)
<
0 and det
⁡
(
J
(
E
∗
)
)
>
0 if a
11
<
0. This means that the two eigenvalues of the Jacobian matrix J
(
E
∗
) have a negative real part. Therefore, the interior equilibrium E
∗
=
(
u
∗
,
v
∗
) is locally asymptotically stable.


### 2.4. Global Stability of Equilibria

Theorem 6.
Assuming that(25)
a
e
u
-
+
r
-
m
<
0
,
-
b
+
u
-
2
+
1
+
u
-
2
(
u
-
2
+
1
)
<
0
,then the boundary equilibrium E
1
=
(
u
-
,
0
) of system (5) is globally asymptotically stable.Proof.
Obviously,u
- satisfies the equation I
-
b
u
-
+
u
-
2
/
(
1
+
u
-
2
)
=
0.
To prove the above statement, let us consider the following Lyapunov function:(26)
V
(
u
,
v
)
=
∫
u
¯
u
e
-
u
-
e
d
e
+
1
e
∫
0
v
d
f
.
It is easy to show thatV
(
u
,
v
)
≥
0 for all t
≥
0. Differentiating V
(
u
,
v
) along the solutions of system (5), (27)
d
V
(
u
,
v
)
d
t
=
u
-
u
-
u
d
u
d
t
+
1
e
d
v
d
t
=
u
-
u
-
u
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
)
+
1
e
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
]
=
u
-
u
-
u
[
-
b
(
u
-
u
-
)
+
u
2
1
+
u
2
-
u
-
2
1
+
u
-
2
-
a
u
v
]
+
1
e
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
]
=
u
-
u
-
2
u
[
-
b
+
1
1
+
u
-
2
u
+
u
-
1
+
u
2
]
+
1
e
[
-
r
k
v
2
+
(
a
e
u
-
+
r
-
m
)
v
]
.Assuming the following function,(28)
h
(
u
)
=
1
1
+
u
-
2
u
+
u
-
1
+
u
2
,then, h
′
(
u
)
=
(
1
/
1
+
u
-
2
)
(
(
1
-
u
2
-
2
u
-
u
)
/
(
1
+
u
2
)
2
).
Leth
′
(
u
)
=
0; then, u
=
-
u
-
±
1
+
u
-
2.
Obviously,lim
u
→
0
⁡
h
(
u
)
=
u
-
/
(
1
+
u
-
2
), lim
u
→
+
∞
⁡
h
(
u
)
=
0.
Therefore,h
(
u
)
≤
h
(
-
u
-
+
1
+
u
-
2
)
=
(
u
-
2
+
1
+
u
-
)
/
2
(
u
-
2
+
1
).
Then(29)
d
V
(
u
,
v
)
d
t
≤
u
-
u
-
2
u
[
-
b
+
u
-
2
+
1
+
u
-
2
(
u
-
2
+
1
)
]
+
1
e
[
-
r
k
v
2
+
(
a
e
u
-
+
r
-
m
)
v
]
.Obviously, if a
e
u
-
+
r
-
m
<
0 and -
b
+
(
u
-
2
+
1
+
u
-
)
/
2
(
u
-
2
+
1
)
<
0, then d
V
(
u
,
v
)
/
d
t
<
0 strictly for all u
,
v
>
0 except the boundary equilibrium (
u
-
,
0
), where d
V
(
u
,
v
)
/
d
t
=
0.
Therefore,V
(
u
,
v
) satisfies Lyapunov’s asymptotic stability theorem, and the boundary equilibrium (
u
-
,
0
) of system (5) is globally asymptotically stable. This completes the proof.Remark 7.
It is easy to verify that if condition (25) holds, then condition (18) is true.Theorem 8.
The positive equilibriumE
∗
=
(
u
∗
,
v
∗
) of system (5) is globally asymptotically stable if -
4
I
(
I
+
1
)
u
∗
4
+
(
1
-
4
I
-
8
I
2
)
u
∗
2
-
4
I
2
<
0.Proof.
Obviously,(
u
∗
,
v
∗
) satisfies the equations(30)
I
-
b
u
∗
+
u
∗
2
1
+
u
∗
2
-
a
u
∗
v
∗
=
0
,
r
(
1
-
v
∗
k
)
+
a
e
u
∗
-
m
=
0
.Now let us construct the following Lyapunov function:(31)
V
(
u
,
v
)
=
u
-
u
∗
-
u
∗
ln
⁡
u
u
∗
+
1
e
(
v
-
v
∗
-
v
∗
ln
⁡
v
v
∗
)
.Obviously, V
(
u
,
v
) is continuous for all u
,
v
>
0.
By simple computation,(32)
∂
V
∂
u
=
1
-
u
∗
u
,
∂
V
∂
v
=
1
e
(
1
-
v
∗
v
)
.This shows that the positive equilibrium (
u
∗
,
v
∗
) is the only extremum of the function V
(
u
,
v
) in the first quadrant. Obviously, (33)
lim
u
→
0
V
(
u
,
v
)
=
lim
v
→
0
V
(
u
,
v
)
=
lim
u
→
+
∞
V
(
u
,
v
)
=
lim
v
→
+
∞
V
(
u
,
v
)
=
+
∞
.Equations  (32) and (33) show that the positive equilibrium (
u
∗
,
v
∗
) is the global minimum in the first quadrant.
Therefore,V
(
u
,
v
)
≥
V
(
u
∗
,
v
∗
)
=
0 for all u
,
v
>
0.
DifferentiatingV
(
u
,
v
) along the solutions of model (5), (34)
d
V
(
u
,
v
)
d
t
=
(
1
-
u
∗
u
)
d
u
d
t
+
1
e
(
1
-
v
∗
v
)
d
v
d
t
=
u
-
u
∗
u
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
)
+
1
e
(
v
-
v
∗
)
[
r
(
1
-
v
k
)
+
a
e
u
-
m
]
=
-
b
u
u
-
u
∗
2
-
a
v
∗
u
u
-
u
∗
2
-
r
k
e
v
-
v
∗
2
+
(
u
+
u
∗
)
u
-
u
∗
2
u
(
1
+
u
2
)
(
1
+
u
∗
2
)
=
[
u
+
u
∗
u
1
+
u
2
1
+
u
∗
2
-
b
+
a
v
∗
u
]
·
u
-
u
∗
2
-
r
k
e
v
-
v
∗
2
=
H
1
-
H
2
,
H
1
=
[
u
+
u
∗
u
(
1
+
u
2
)
(
1
+
u
∗
2
)
-
b
+
a
v
∗
u
]
u
-
u
∗
2
=
-
I
+
1
u
∗
2
+
I
u
2
+
u
∗
u
-
I
1
+
u
∗
2
u
∗
u
1
+
u
2
1
+
u
∗
2
·
u
-
u
∗
2
≤
-
4
I
I
+
1
u
∗
4
+
1
-
4
I
-
8
I
2
u
∗
2
-
4
I
2
u
∗
u
1
+
u
2
1
+
u
∗
2
4
I
+
1
u
∗
2
+
I
·
u
-
u
∗
2
.
Obviously,H
2
≥
0. Furthermore, if -
4
I
(
I
+
1
)
u
∗
4
+
(
1
-
4
I
-
8
I
2
)
u
∗
2
-
4
I
2
<
0, H
1
≤
0, and, therefore, d
V
(
u
,
v
)
/
d
t
=
H
1
-
H
2
<
0 strictly for all u
,
v
>
0 except the positive equilibrium (
u
∗
,
v
∗
), where d
V
(
u
,
v
)
/
d
t
=
0.
Therefore,V
(
u
,
v
) satisfies Lyapunov’s asymptotic stability theorem, and the positive equilibrium (
u
∗
,
v
∗
) of system (5) is globally asymptotically stable. This completes the proof.


## 2.1. Positivity and Boundedness of the Solutions

Theorem 1.
Suppose thata
e
L
+
r
-
m
>
0
(
L
=
max
⁡
u
(
0
)
,
(
I
+
1
)
/
b
) holds; then all the solutions of system (5) with initial conditions are positive and bounded for all t
≥
0.Proof.
From the first equation of system (5), (6)
d
u
d
t
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
≥
-
b
u
+
u
2
1
+
u
2
-
a
u
v
.It follows that d
u
/
d
t
≥
(
-
b
+
u
/
(
1
+
u
2
)
-
a
v
)
u.
Hence,u
(
t
)
≥
u
(
0
)
exp
⁡
∫
0
t
[
-
b
+
u
(
s
)
/
(
1
+
(
u
s
)
2
)
-
a
v
(
s
)
]
d
s
>
0.
From the second equation of system (5), it can be determined that(7)
v
(
t
)
=
v
(
0
)
exp
⁡
∫
0
t
[
r
(
1
-
v
(
s
)
k
)
+
a
e
u
(
s
)
-
m
]
d
s
>
0
,
d
u
d
t
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
≤
I
-
b
u
+
u
2
1
+
u
2
≤
I
+
1
-
b
u
,leading to d
u
/
d
t
≤
I
+
1
-
b
u.
Then,(8)
u
(
t
)
≤
I
+
1
b
(
u
(
0
)
-
I
+
1
b
)
e
-
b
t
≤
max
⁡
u
(
0
)
,
I
+
1
b
=
L
,
d
v
d
t
=
v
[
r
(
1
-
v
k
)
+
a
e
u
-
m
]
≤
v
[
r
(
1
-
v
k
)
+
a
e
L
-
m
]
=
v
(
a
e
L
+
r
-
m
-
r
k
v
)
.If a
e
L
+
r
-
m
>
0, then, for ε
>
0, there exists T
(
>
0
) such that, for any t
>
T, (9)
v
(
t
)
≤
k
(
a
e
L
+
r
-
m
)
r
+
ε
.Moreover, v
(
t
)
=
v
(
0
)
exp
⁡
∫
0
t
[
r
(
1
-
v
(
s
)
/
k
)
+
a
e
u
(
s
)
-
m
]
d
s is continuous when t
∈
[
0
,
T
].
This means that there existsM
(
>
0
) such that v
(
t
)
≤
M
t
∈
[
0
,
T
].
Therefore,v
(
t
)
≤
max
⁡
k
(
a
e
L
+
r
-
m
)
/
r
+
ε
,
M when t
≥
0, and the proof is complete.


## 2.2. Existence of Equilibria

Theorem 2.
System (5) has boundary equilibrium E
1
=
(
u
-
,
0
) (extinction of algae population).Proof.
IfE
1
=
(
u
-
,
0
) is the boundary equilibrium of (5), then(10)
I
-
b
u
-
+
u
-
2
1
+
u
-
2
=
0
that
is
,
b
u
-
3
-
(
I
+
1
)
u
-
2
+
b
u
-
-
I
=
0
.Assume the following function:(11)
f
(
u
)
=
b
u
3
-
(
I
+
1
)
u
2
+
b
u
-
I
.Obviously,(12)
f
(
0
)
=
-
I
,
lim
u
→
+
∞
f
(
u
)
=
+
∞
.
It is easy to establish that the curve of functionf
(
u
) intersects the u-axis, and, therefore, the existence of the boundary equilibrium is guaranteed, which completes the proof.Theorem 3.
System (5) has a unique interior equilibrium E
∗
=
(
u
∗
,
v
∗
) (coexistence of nutrients and algae) if a
2
k
e
-
r
I
-
r
>
0 and a
k
(
r
-
m
)
+
b
r
>
0.Proof.
IfE
∗
=
(
u
∗
,
v
∗
) is the interior equilibrium of (5), then(13)
I
-
b
u
∗
+
u
∗
2
1
+
u
∗
2
-
a
u
∗
v
∗
=
0
,
r
v
∗
(
1
-
v
∗
k
)
+
a
e
u
∗
v
∗
-
m
v
∗
=
0
.Consequently, u
∗ is the positive root of the fourth-degree equation:(14)
a
2
k
e
u
∗
4
+
a
k
(
r
-
m
)
+
b
r
u
∗
3
+
(
a
2
k
e
-
r
I
-
r
)
u
∗
2
+
a
k
(
r
-
m
)
+
b
r
u
∗
-
r
I
=
0
,
v
∗
=
k
r
(
a
e
u
∗
+
r
-
m
)
.Consider the following function:(15)
h
x
=
a
2
k
e
x
4
+
a
k
(
r
-
m
)
+
b
r
x
3
+
(
a
2
k
e
-
r
I
-
r
)
x
2
+
a
k
(
r
-
m
)
+
b
r
x
-
r
I
.It is easy to establish that(16)
h
′
x
=
4
a
2
k
e
x
3
+
3
a
k
(
r
-
m
)
+
b
r
x
2
+
2
(
a
2
k
e
-
r
I
-
r
)
x
+
a
k
(
r
-
m
)
+
b
r
=
[
4
a
2
k
e
x
2
+
2
(
a
2
k
e
-
r
I
-
r
)
]
x
+
a
k
r
-
m
+
b
r
3
x
2
+
1
.It is also easy to determine that h
′
(
x
)
>
0
(
x
>
0
) if a
2
k
e
-
r
I
-
r
>
0 and a
k
(
r
-
m
)
+
b
r
>
0.
Obviously,(17)
h
(
0
)
=
-
r
I
,
lim
x
→
+
∞
h
(
x
)
=
+
∞
.
It can easily be shown that the curve of functionh
(
x
) intersects the x-axis only once. Hence, the existence and uniqueness of the interior equilibrium is guaranteed, and this completes the proof.


## 2.3. Local Stability of Equilibria

Theorem 4.
The boundary equilibriumE
1
=
(
u
-
,
0
) of system (5) is locally asymptotically stable if (18)
-
b
+
2
u
-
1
+
u
-
2
2
<
0
,
a
e
u
-
+
r
-
m
<
0
.The Jacobian matrix of system (5) is (19)
J
=
-
b
-
a
v
+
2
u
1
+
u
2
2
-
a
u
a
e
v
r
-
2
r
k
v
+
a
e
u
-
m
.The Jacobian matrix of system (5) at E
1
=
(
u
-
,
0
) is(20)
J
(
E
1
)
=
-
b
+
2
u
-
1
+
u
-
2
2
-
a
u
-
0
a
e
u
-
+
r
-
m
.It is clear that the Jacobian matrix J
(
E
1
) has two eigenvalues: (21)
λ
1
=
-
b
+
2
u
-
1
+
u
-
2
2
,
λ
2
=
a
e
u
-
+
r
-
m
.Hence, the equilibrium E
1
=
(
u
-
,
0
) is locally asymptotically stable if (22)
-
b
+
2
u
-
1
+
u
-
2
2
<
0
,
a
e
u
-
+
r
-
m
<
0
.Theorem 5.
The interior equilibriumE
∗
=
(
u
∗
,
v
∗
) of system (5) is locally asymptotically stable if -
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
<
0.Proof.
The Jacobian matrix of system (5) at E
∗
=
(
u
∗
,
v
∗
) is(23)
J
(
E
∗
)
=
J
11
J
12
J
21
J
22
=
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
-
a
u
∗
a
e
v
∗
-
r
k
v
∗
.Moreover, (24)
det
⁡
(
λ
E
-
J
(
E
∗
)
)
=
λ
2
-
tr
⁡
(
J
(
E
∗
)
)
λ
+
det
⁡
(
J
(
E
∗
)
)
,where tr
⁡
(
J
E
∗
)
=
a
11
-
(
r
/
k
)
v
∗, det
⁡
(
J
E
∗
)
=
(
a
2
e
u
∗
-
a
11
(
r
/
k
)
)
v
∗, and a
11
=
-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2.
It can be easily determined thattr
⁡
(
J
(
E
∗
)
)
<
0 and det
⁡
(
J
(
E
∗
)
)
>
0 if a
11
<
0. This means that the two eigenvalues of the Jacobian matrix J
(
E
∗
) have a negative real part. Therefore, the interior equilibrium E
∗
=
(
u
∗
,
v
∗
) is locally asymptotically stable.


## 2.4. Global Stability of Equilibria

Theorem 6.
Assuming that(25)
a
e
u
-
+
r
-
m
<
0
,
-
b
+
u
-
2
+
1
+
u
-
2
(
u
-
2
+
1
)
<
0
,then the boundary equilibrium E
1
=
(
u
-
,
0
) of system (5) is globally asymptotically stable.Proof.
Obviously,u
- satisfies the equation I
-
b
u
-
+
u
-
2
/
(
1
+
u
-
2
)
=
0.
To prove the above statement, let us consider the following Lyapunov function:(26)
V
(
u
,
v
)
=
∫
u
¯
u
e
-
u
-
e
d
e
+
1
e
∫
0
v
d
f
.
It is easy to show thatV
(
u
,
v
)
≥
0 for all t
≥
0. Differentiating V
(
u
,
v
) along the solutions of system (5), (27)
d
V
(
u
,
v
)
d
t
=
u
-
u
-
u
d
u
d
t
+
1
e
d
v
d
t
=
u
-
u
-
u
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
)
+
1
e
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
]
=
u
-
u
-
u
[
-
b
(
u
-
u
-
)
+
u
2
1
+
u
2
-
u
-
2
1
+
u
-
2
-
a
u
v
]
+
1
e
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
]
=
u
-
u
-
2
u
[
-
b
+
1
1
+
u
-
2
u
+
u
-
1
+
u
2
]
+
1
e
[
-
r
k
v
2
+
(
a
e
u
-
+
r
-
m
)
v
]
.Assuming the following function,(28)
h
(
u
)
=
1
1
+
u
-
2
u
+
u
-
1
+
u
2
,then, h
′
(
u
)
=
(
1
/
1
+
u
-
2
)
(
(
1
-
u
2
-
2
u
-
u
)
/
(
1
+
u
2
)
2
).
Leth
′
(
u
)
=
0; then, u
=
-
u
-
±
1
+
u
-
2.
Obviously,lim
u
→
0
⁡
h
(
u
)
=
u
-
/
(
1
+
u
-
2
), lim
u
→
+
∞
⁡
h
(
u
)
=
0.
Therefore,h
(
u
)
≤
h
(
-
u
-
+
1
+
u
-
2
)
=
(
u
-
2
+
1
+
u
-
)
/
2
(
u
-
2
+
1
).
Then(29)
d
V
(
u
,
v
)
d
t
≤
u
-
u
-
2
u
[
-
b
+
u
-
2
+
1
+
u
-
2
(
u
-
2
+
1
)
]
+
1
e
[
-
r
k
v
2
+
(
a
e
u
-
+
r
-
m
)
v
]
.Obviously, if a
e
u
-
+
r
-
m
<
0 and -
b
+
(
u
-
2
+
1
+
u
-
)
/
2
(
u
-
2
+
1
)
<
0, then d
V
(
u
,
v
)
/
d
t
<
0 strictly for all u
,
v
>
0 except the boundary equilibrium (
u
-
,
0
), where d
V
(
u
,
v
)
/
d
t
=
0.
Therefore,V
(
u
,
v
) satisfies Lyapunov’s asymptotic stability theorem, and the boundary equilibrium (
u
-
,
0
) of system (5) is globally asymptotically stable. This completes the proof.Remark 7.
It is easy to verify that if condition (25) holds, then condition (18) is true.Theorem 8.
The positive equilibriumE
∗
=
(
u
∗
,
v
∗
) of system (5) is globally asymptotically stable if -
4
I
(
I
+
1
)
u
∗
4
+
(
1
-
4
I
-
8
I
2
)
u
∗
2
-
4
I
2
<
0.Proof.
Obviously,(
u
∗
,
v
∗
) satisfies the equations(30)
I
-
b
u
∗
+
u
∗
2
1
+
u
∗
2
-
a
u
∗
v
∗
=
0
,
r
(
1
-
v
∗
k
)
+
a
e
u
∗
-
m
=
0
.Now let us construct the following Lyapunov function:(31)
V
(
u
,
v
)
=
u
-
u
∗
-
u
∗
ln
⁡
u
u
∗
+
1
e
(
v
-
v
∗
-
v
∗
ln
⁡
v
v
∗
)
.Obviously, V
(
u
,
v
) is continuous for all u
,
v
>
0.
By simple computation,(32)
∂
V
∂
u
=
1
-
u
∗
u
,
∂
V
∂
v
=
1
e
(
1
-
v
∗
v
)
.This shows that the positive equilibrium (
u
∗
,
v
∗
) is the only extremum of the function V
(
u
,
v
) in the first quadrant. Obviously, (33)
lim
u
→
0
V
(
u
,
v
)
=
lim
v
→
0
V
(
u
,
v
)
=
lim
u
→
+
∞
V
(
u
,
v
)
=
lim
v
→
+
∞
V
(
u
,
v
)
=
+
∞
.Equations  (32) and (33) show that the positive equilibrium (
u
∗
,
v
∗
) is the global minimum in the first quadrant.
Therefore,V
(
u
,
v
)
≥
V
(
u
∗
,
v
∗
)
=
0 for all u
,
v
>
0.
DifferentiatingV
(
u
,
v
) along the solutions of model (5), (34)
d
V
(
u
,
v
)
d
t
=
(
1
-
u
∗
u
)
d
u
d
t
+
1
e
(
1
-
v
∗
v
)
d
v
d
t
=
u
-
u
∗
u
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
)
+
1
e
(
v
-
v
∗
)
[
r
(
1
-
v
k
)
+
a
e
u
-
m
]
=
-
b
u
u
-
u
∗
2
-
a
v
∗
u
u
-
u
∗
2
-
r
k
e
v
-
v
∗
2
+
(
u
+
u
∗
)
u
-
u
∗
2
u
(
1
+
u
2
)
(
1
+
u
∗
2
)
=
[
u
+
u
∗
u
1
+
u
2
1
+
u
∗
2
-
b
+
a
v
∗
u
]
·
u
-
u
∗
2
-
r
k
e
v
-
v
∗
2
=
H
1
-
H
2
,
H
1
=
[
u
+
u
∗
u
(
1
+
u
2
)
(
1
+
u
∗
2
)
-
b
+
a
v
∗
u
]
u
-
u
∗
2
=
-
I
+
1
u
∗
2
+
I
u
2
+
u
∗
u
-
I
1
+
u
∗
2
u
∗
u
1
+
u
2
1
+
u
∗
2
·
u
-
u
∗
2
≤
-
4
I
I
+
1
u
∗
4
+
1
-
4
I
-
8
I
2
u
∗
2
-
4
I
2
u
∗
u
1
+
u
2
1
+
u
∗
2
4
I
+
1
u
∗
2
+
I
·
u
-
u
∗
2
.
Obviously,H
2
≥
0. Furthermore, if -
4
I
(
I
+
1
)
u
∗
4
+
(
1
-
4
I
-
8
I
2
)
u
∗
2
-
4
I
2
<
0, H
1
≤
0, and, therefore, d
V
(
u
,
v
)
/
d
t
=
H
1
-
H
2
<
0 strictly for all u
,
v
>
0 except the positive equilibrium (
u
∗
,
v
∗
), where d
V
(
u
,
v
)
/
d
t
=
0.
Therefore,V
(
u
,
v
) satisfies Lyapunov’s asymptotic stability theorem, and the positive equilibrium (
u
∗
,
v
∗
) of system (5) is globally asymptotically stable. This completes the proof.


## 3. Stability Analysis in the Presence of Diffusion and Turing Instability

### 3.1. Stability Analysis of Equilibria in the Presence of Diffusion

Theorem 9.
Ifa
e
u
-
+
r
-
m
>
0, the steady state E
1
=
(
u
-
,
0
) of the diffusive system (2) is unstable.Proof.
To investigate the stability of the equilibriumE
1, let us consider the corresponding eigenvalue problem of the linearized operator around E
1.
The linearization of system (2) at steady state E
1
=
(
u
-
,
0
) can be expressed as (35)
∂
Z
∂
t
=
(
D
Δ
+
J
(
E
1
)
)
Z
,where Z
=
(
u
,
v
)
T, D
=
diag
(
d
1
,
d
2
), and (36)
J
(
E
1
)
=
-
b
+
2
u
-
1
+
u
-
2
2
-
a
u
-
0
a
e
u
-
+
r
-
m
.From the above, the linearized result of system (2) around E
1 is(37)
∂
u
∂
t
=
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
¯
2
2
]
u
-
a
u
-
v
,
∂
v
∂
t
=
d
2
Δ
v
+
(
a
e
u
-
+
r
-
m
)
v
,
∂
u
∂
n
∂
Ω
=
∂
v
∂
n
∂
Ω
=
0
.The corresponding characteristic equation of (37) can be obtained as(38)
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
-
2
2
]
u
-
a
u
-
v
=
λ
u
,
d
2
Δ
v
+
(
a
e
u
-
+
r
-
m
)
v
=
λ
v
,
∂
u
∂
n
∂
Ω
=
∂
v
∂
n
∂
Ω
=
0
,where λ is an eigenvalue of (38) and the corresponding eigenvector is (
u
,
v
). If v
≠
0, then λ is an eigenvalue of the operator d
2
Δ
+
(
a
e
u
-
+
r
-
m
) with a homogeneous Neumann boundary condition. Therefore, λ must be real. In the same way, λ is also real if u
≠
0. Hence, all eigenvalues of (38) are real. Let λ
max
⁡ be the largest eigenvalue. Consider the principal eigenvalue λ
^ of the following equation:(39)
d
2
Δ
v
+
(
a
e
u
-
+
r
-
m
)
v
=
λ
v
,
∂
v
∂
n
∂
Ω
=
0
.
From the above equation,λ
^
>
0 and the associated eigenvector v
^
>
0 if a
e
u
-
+
r
-
m
>
0. It is claimed here that λ
^ is also an eigenvalue of (38). In fact, u
^
>
0 is taken to be a solution of(40)
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
-
2
2
]
u
-
a
u
-
v
=
λ
^
u
,
∂
v
∂
n
∂
Ω
=
0
;then (
u
,
v
)
=
(
u
^
,
v
^
) satisfies (38) with λ
=
λ
^. Therefore, λ
^
>
0 is an eigenvalue of (38). Hence, λ
max
⁡
≥
λ
^
>
0, and the steady state E
1 of the diffusive system (2) is unstable. This completes the proof.Theorem 10.
If-
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0, the diffusive system (2) is stable at E
1.Proof.
In the case-
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0, let (
u
^
,
v
^
) be the principal eigenvector of (38) corresponding to the largest eigenvalue λ
max
⁡.
Ifv
^
≠
0, then λ
max
⁡ is also an eigenvalue of (39). Therefore, λ
max
⁡
<
0 can be obtained if a
e
u
-
+
r
-
m
<
0 holds.
Ifv
^
=
0, then u
^
≠
0. Hence, λ
max
⁡ is an eigenvalue of (41)
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
-
2
2
]
u
=
λ
u
,
∂
u
∂
n
∂
Ω
=
0
.
Obviously, the largest eigenvalue of (41) is -
b
+
2
u
-
/
(
1
+
u
-
2
)
2. Therefore, λ
max
⁡
<
0 if -
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0. Therefore, the diffusive system (2) is stable at E
1 if -
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0. This completes the proof.To prove Theorem12, let us introduce the following lemma [33].Lemma 11.
Consider the following equation:(42)
∂
u
∂
t
=
d
1
Δ
u
+
F
1
(
u
,
w
)
,
t
>
0
,
x
∈
Ω
∂
w
∂
t
=
d
2
Δ
w
+
F
2
(
u
,
w
)
,
t
>
0
,
x
∈
Ω
∂
u
∂
n
=
∂
w
∂
n
=
0
,
t
>
0
,
x
∈
∂
Ω
u
0
,
x
≥
0
,
w
0
,
x
≥
0
,
x
∈
Ω
.Let (
u
~
,
w
~
) be a steady state of (42); that is, F
1
(
u
~
,
w
~
)
=
F
2
(
u
~
,
w
~
)
=
0.
Ifθ
1
-
θ
4
<
0, θ
2
≥
0, θ
3
≥
0, θ
2
θ
3
-
θ
1
θ
4
>
0, where(43)
θ
1
=
∂
F
1
∂
u
u
=
u
~
,
w
=
w
~
θ
2
=
-
∂
F
1
∂
w
u
=
u
~
,
w
=
w
~
,
θ
3
=
∂
F
2
∂
u
u
=
u
~
,
w
=
w
~
θ
4
=
-
∂
F
2
∂
w
u
=
u
~
,
w
=
w
~
,then (
u
~
,
w
~
) is uniformly asymptotically stable. Furthermore, if θ
2
θ
3
-
θ
1
θ
4
<
0, then (
u
~
,
w
~
) is unstable.Theorem 12.
The following conclusions hold.(
1
)
If(44)
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
<
0
,
then the positive constant solutionE
∗
=
(
u
∗
,
v
∗
) of system (2) is uniformly asymptotically stable.
(
2
)
If(45)
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
>
k
r
a
2
e
u
∗
,
then the positive constant solutionE
∗
=
(
u
∗
,
v
∗
) of system (2) is unstable.Proof.
Consider(46)
F
1
u
,
v
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
,
F
2
(
u
,
v
)
=
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
,
θ
1
=
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
,
θ
2
=
a
u
∗
,
θ
3
=
a
e
v
∗
,
θ
4
=
r
k
v
∗
.Obviously, θ
1
-
θ
4
<
0, θ
2
>
0, θ
3
>
0, θ
2
θ
3
-
θ
1
θ
4
>
0 if -
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
<
0.
Applying Lemma11, the first conclusion of Theorem 12 follows.
On the other hand, if-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
>
(
k
/
r
)
a
2
e
u
∗, it can easily be established that θ
2
θ
3
-
θ
1
θ
4
=
a
2
e
u
∗
v
∗
-
(
r
/
k
)
v
∗
[
-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
]
<
0. Therefore, E
∗ of system (2) is unstable by Lemma 11. This completes the proof.Theorem 13.
If(47)
-
b
-
a
v
∗
+
1
1
+
u
∗
2
1
1
+
1
+
u
∗
2
-
u
∗
2
<
0
,then the positive equilibrium E
∗
=
(
u
∗
,
v
∗
) of system (2) is globally asymptotically stable.Proof.
To prove Theorem13, it is necessary to construct a Lyapunov function. Define (48)
H
(
u
,
v
)
=
∫
u
∗
u
ξ
-
u
∗
ξ
d
ξ
+
1
e
∫
v
∗
v
η
-
v
∗
η
d
η
,
H
1
(
t
)
=
∬
Ω
H
(
u
,
v
)
d
Ω
.
Obviously,H
(
u
,
v
)
≥
0, H
(
u
,
v
)
=
0 if and only if (
u
,
v
)
=
(
u
∗
,
v
∗
). Calculating the derivative of H
1
(
t
) along the solutions of model (2), (49)
d
H
1
(
t
)
d
t
=
∬
Ω
(
∂
H
∂
u
∂
u
∂
t
+
∂
H
∂
v
∂
v
∂
t
)
d
Ω
=
∬
Ω
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
+
d
1
Δ
u
)
∂
H
∂
u
m
m
m
m
c
+
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
+
d
2
Δ
v
]
m
m
m
m
c
u
2
1
+
u
2
·
∂
H
∂
v
d
Ω
=
∬
Ω
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
)
∂
H
∂
u
m
m
m
m
c
u
2
1
+
u
2
+
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
]
∂
H
∂
v
d
Ω
+
∬
Ω
(
d
1
Δ
u
∂
H
∂
u
+
d
2
Δ
v
∂
H
∂
v
)
d
Ω
=
∬
Ω
u
-
u
∗
2
u
u
+
u
∗
1
+
u
2
1
+
u
∗
2
-
b
-
a
v
∗
d
Ω
-
r
k
e
∬
Ω
v
-
v
∗
2
d
Ω
-
d
1
∬
Ω
u
∗
u
2
∂
u
∂
x
2
+
∂
u
∂
y
2
d
Ω
-
d
2
e
∬
Ω
v
∗
v
2
∂
v
∂
x
2
+
∂
v
∂
y
2
d
Ω
=
I
1
-
I
2
-
I
3
-
I
4
.Obviously, I
2
≥
0, I
3
≥
0, and I
4
≥
0.
Consider the following function:(50)
f
(
x
)
=
x
+
u
∗
1
+
x
2
(
x
≥
0
)
.Then, f
′
(
x
)
=
(
-
x
2
-
2
u
∗
x
+
1
)
/
(
1
+
x
2
)
2. Because lim
x
→
0
⁡
f
(
x
)
=
u
∗
,
lim
x
→
+
∞
⁡
f
(
x
)
=
0, then(51)
f
(
x
)
≤
f
(
-
u
∗
+
u
∗
2
+
1
)
=
u
∗
2
+
1
1
+
-
u
∗
+
u
∗
2
+
1
2
,
I
1
=
∬
Ω
u
-
u
∗
2
u
u
+
u
∗
(
1
+
u
2
)
(
1
+
u
∗
2
)
-
b
-
a
v
∗
d
Ω
≤
∬
Ω
u
-
u
∗
2
u
-
b
-
a
v
∗
+
1
1
+
u
∗
2
1
1
+
u
∗
2
+
1
-
u
∗
2
m
m
m
m
m
m
m
m
m
c
m
·
1
1
+
u
∗
2
+
1
-
u
∗
2
d
Ω
.
Therefore, if-
b
-
a
v
∗
+
(
1
/
1
+
u
∗
2
)
(
1
/
(
1
+
(
1
+
u
∗
2
-
u
∗
)
2
)
)
<
0, I
1
<
0, which is equivalent to d
H
1
(
t
)
/
d
t
<
0.
Thus,H
1
(
t
) satisfies Lyapunov’s asymptotic stability theorem, and the positive equilibrium (
u
∗
,
v
∗
) of system (2) is globally asymptotically stable. This completes the proof.Remark 14.
It is easy to verify that if condition (47) holds, then condition (44) is true.


### 3.2. Turing Instability

To analyze the spatial system and how a small heterogeneous perturbation of a homogeneous steady state develops over a long time period, small space- and time-dependent perturbations forE
∗ of system (2) will be considered:(52)
u
=
u
∗
+
ε
exp
⁡
(
K
z
i
+
λ
t
)
,
v
=
v
∗
+
τ
exp
⁡
K
z
i
+
λ
t
,where ε, τ are small enough and K is the wave number. Substituting (52) into (2) yields its characteristic matrix:(53)
A
=
J
11
-
d
1
K
2
J
12
J
21
J
22
-
d
2
K
2
.Correspondingly, λ
1 and λ
2 are the roots of the following characteristic equation:(54)
λ
2
-
t
r
K
λ
+
Δ
K
=
0
,
t
r
K
=
J
11
+
J
22
-
(
d
1
+
d
2
)
K
2, Δ
K
=
d
1
d
2
K
4
-
(
d
1
J
22
+
d
2
J
11
)
K
2
+
J
11
J
22
-
J
12
J
21.The roots of (54) yield the dispersion relation:(55)
λ
1,2
(
K
)
=
t
r
K
±
t
r
K
2
-
4
Δ
K
2
.As is well known, Turing instability means that the stable equilibrium point is driven to become unstable by the local dynamics and diffusion of species. The stability conditions forE
∗of system (5) are t
r
0
=
J
11
+
J
22
<
0 and Δ
0
=
J
11
J
22
-
J
12
J
21
>
0. It is clear that t
r
K
<
t
r
0
<
0. Therefore, the stability of equilibrium E
∗ of system (2) changes with the sign of Δ
K. It is easy to establish that Δ
K
<
0 for K
1
2
<
K
2
<
K
2
2, where(56)
K
1,2
2
=
d
1
J
22
+
d
2
J
11
d
1
J
22
+
d
2
J
11
2
-
4
d
1
d
2
(
J
11
J
22
-
J
12
J
21
)
m
m
c
∓
d
1
J
22
+
d
2
J
11
2
-
4
d
1
d
2
(
J
11
J
22
-
J
12
J
21
)
·
2
d
1
d
2
-
1
.IfK
1,2
2 has positive values, then the range of instability for E
∗ can be obtained, which is called the Turing space. To illustrate the Turing space, the dispersion relations corresponding to several values of the bifurcation parameter I are plotted in Figure 1. The other parameters come from the actual monitoring data from the Sanyang wetland from 2012 to 2013 and from the literature:(57)
a
=
0.6
,
b
=
0.015
,
e
=
0.28
,
r
=
0.533
,
k
=
1
,
m
=
0.24
,
d
1
=
0.1
,
d
2
=
1.4
.Figure 1
Variation of the dispersion relation of model (2) around E
∗. The yellow line corresponds to I
=
0.0055, the green to I
∗
=
0.006195, and the red to I
=
0.008.In Figure1, the green line corresponds to the critical Turing value, I
∗
=
0.006195. When I
=
0.0055
<
I
∗ (the yellow line in Figure 1), Turing instability occurs, whereas when I
=
0.008
>
I
∗ (the red line in Figure 1), the Turing instability decays.


## 3.1. Stability Analysis of Equilibria in the Presence of Diffusion

Theorem 9.
Ifa
e
u
-
+
r
-
m
>
0, the steady state E
1
=
(
u
-
,
0
) of the diffusive system (2) is unstable.Proof.
To investigate the stability of the equilibriumE
1, let us consider the corresponding eigenvalue problem of the linearized operator around E
1.
The linearization of system (2) at steady state E
1
=
(
u
-
,
0
) can be expressed as (35)
∂
Z
∂
t
=
(
D
Δ
+
J
(
E
1
)
)
Z
,where Z
=
(
u
,
v
)
T, D
=
diag
(
d
1
,
d
2
), and (36)
J
(
E
1
)
=
-
b
+
2
u
-
1
+
u
-
2
2
-
a
u
-
0
a
e
u
-
+
r
-
m
.From the above, the linearized result of system (2) around E
1 is(37)
∂
u
∂
t
=
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
¯
2
2
]
u
-
a
u
-
v
,
∂
v
∂
t
=
d
2
Δ
v
+
(
a
e
u
-
+
r
-
m
)
v
,
∂
u
∂
n
∂
Ω
=
∂
v
∂
n
∂
Ω
=
0
.The corresponding characteristic equation of (37) can be obtained as(38)
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
-
2
2
]
u
-
a
u
-
v
=
λ
u
,
d
2
Δ
v
+
(
a
e
u
-
+
r
-
m
)
v
=
λ
v
,
∂
u
∂
n
∂
Ω
=
∂
v
∂
n
∂
Ω
=
0
,where λ is an eigenvalue of (38) and the corresponding eigenvector is (
u
,
v
). If v
≠
0, then λ is an eigenvalue of the operator d
2
Δ
+
(
a
e
u
-
+
r
-
m
) with a homogeneous Neumann boundary condition. Therefore, λ must be real. In the same way, λ is also real if u
≠
0. Hence, all eigenvalues of (38) are real. Let λ
max
⁡ be the largest eigenvalue. Consider the principal eigenvalue λ
^ of the following equation:(39)
d
2
Δ
v
+
(
a
e
u
-
+
r
-
m
)
v
=
λ
v
,
∂
v
∂
n
∂
Ω
=
0
.
From the above equation,λ
^
>
0 and the associated eigenvector v
^
>
0 if a
e
u
-
+
r
-
m
>
0. It is claimed here that λ
^ is also an eigenvalue of (38). In fact, u
^
>
0 is taken to be a solution of(40)
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
-
2
2
]
u
-
a
u
-
v
=
λ
^
u
,
∂
v
∂
n
∂
Ω
=
0
;then (
u
,
v
)
=
(
u
^
,
v
^
) satisfies (38) with λ
=
λ
^. Therefore, λ
^
>
0 is an eigenvalue of (38). Hence, λ
max
⁡
≥
λ
^
>
0, and the steady state E
1 of the diffusive system (2) is unstable. This completes the proof.Theorem 10.
If-
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0, the diffusive system (2) is stable at E
1.Proof.
In the case-
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0, let (
u
^
,
v
^
) be the principal eigenvector of (38) corresponding to the largest eigenvalue λ
max
⁡.
Ifv
^
≠
0, then λ
max
⁡ is also an eigenvalue of (39). Therefore, λ
max
⁡
<
0 can be obtained if a
e
u
-
+
r
-
m
<
0 holds.
Ifv
^
=
0, then u
^
≠
0. Hence, λ
max
⁡ is an eigenvalue of (41)
d
1
Δ
u
+
[
-
b
+
2
u
-
1
+
u
-
2
2
]
u
=
λ
u
,
∂
u
∂
n
∂
Ω
=
0
.
Obviously, the largest eigenvalue of (41) is -
b
+
2
u
-
/
(
1
+
u
-
2
)
2. Therefore, λ
max
⁡
<
0 if -
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0. Therefore, the diffusive system (2) is stable at E
1 if -
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0. This completes the proof.To prove Theorem12, let us introduce the following lemma [33].Lemma 11.
Consider the following equation:(42)
∂
u
∂
t
=
d
1
Δ
u
+
F
1
(
u
,
w
)
,
t
>
0
,
x
∈
Ω
∂
w
∂
t
=
d
2
Δ
w
+
F
2
(
u
,
w
)
,
t
>
0
,
x
∈
Ω
∂
u
∂
n
=
∂
w
∂
n
=
0
,
t
>
0
,
x
∈
∂
Ω
u
0
,
x
≥
0
,
w
0
,
x
≥
0
,
x
∈
Ω
.Let (
u
~
,
w
~
) be a steady state of (42); that is, F
1
(
u
~
,
w
~
)
=
F
2
(
u
~
,
w
~
)
=
0.
Ifθ
1
-
θ
4
<
0, θ
2
≥
0, θ
3
≥
0, θ
2
θ
3
-
θ
1
θ
4
>
0, where(43)
θ
1
=
∂
F
1
∂
u
u
=
u
~
,
w
=
w
~
θ
2
=
-
∂
F
1
∂
w
u
=
u
~
,
w
=
w
~
,
θ
3
=
∂
F
2
∂
u
u
=
u
~
,
w
=
w
~
θ
4
=
-
∂
F
2
∂
w
u
=
u
~
,
w
=
w
~
,then (
u
~
,
w
~
) is uniformly asymptotically stable. Furthermore, if θ
2
θ
3
-
θ
1
θ
4
<
0, then (
u
~
,
w
~
) is unstable.Theorem 12.
The following conclusions hold.(
1
)
If(44)
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
<
0
,
then the positive constant solutionE
∗
=
(
u
∗
,
v
∗
) of system (2) is uniformly asymptotically stable.
(
2
)
If(45)
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
>
k
r
a
2
e
u
∗
,
then the positive constant solutionE
∗
=
(
u
∗
,
v
∗
) of system (2) is unstable.Proof.
Consider(46)
F
1
u
,
v
=
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
,
F
2
(
u
,
v
)
=
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
,
θ
1
=
-
b
-
a
v
∗
+
2
u
∗
1
+
u
∗
2
2
,
θ
2
=
a
u
∗
,
θ
3
=
a
e
v
∗
,
θ
4
=
r
k
v
∗
.Obviously, θ
1
-
θ
4
<
0, θ
2
>
0, θ
3
>
0, θ
2
θ
3
-
θ
1
θ
4
>
0 if -
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
<
0.
Applying Lemma11, the first conclusion of Theorem 12 follows.
On the other hand, if-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
>
(
k
/
r
)
a
2
e
u
∗, it can easily be established that θ
2
θ
3
-
θ
1
θ
4
=
a
2
e
u
∗
v
∗
-
(
r
/
k
)
v
∗
[
-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
]
<
0. Therefore, E
∗ of system (2) is unstable by Lemma 11. This completes the proof.Theorem 13.
If(47)
-
b
-
a
v
∗
+
1
1
+
u
∗
2
1
1
+
1
+
u
∗
2
-
u
∗
2
<
0
,then the positive equilibrium E
∗
=
(
u
∗
,
v
∗
) of system (2) is globally asymptotically stable.Proof.
To prove Theorem13, it is necessary to construct a Lyapunov function. Define (48)
H
(
u
,
v
)
=
∫
u
∗
u
ξ
-
u
∗
ξ
d
ξ
+
1
e
∫
v
∗
v
η
-
v
∗
η
d
η
,
H
1
(
t
)
=
∬
Ω
H
(
u
,
v
)
d
Ω
.
Obviously,H
(
u
,
v
)
≥
0, H
(
u
,
v
)
=
0 if and only if (
u
,
v
)
=
(
u
∗
,
v
∗
). Calculating the derivative of H
1
(
t
) along the solutions of model (2), (49)
d
H
1
(
t
)
d
t
=
∬
Ω
(
∂
H
∂
u
∂
u
∂
t
+
∂
H
∂
v
∂
v
∂
t
)
d
Ω
=
∬
Ω
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
+
d
1
Δ
u
)
∂
H
∂
u
m
m
m
m
c
+
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
+
d
2
Δ
v
]
m
m
m
m
c
u
2
1
+
u
2
·
∂
H
∂
v
d
Ω
=
∬
Ω
(
I
-
b
u
+
u
2
1
+
u
2
-
a
u
v
)
∂
H
∂
u
m
m
m
m
c
u
2
1
+
u
2
+
[
r
v
(
1
-
v
k
)
+
a
e
u
v
-
m
v
]
∂
H
∂
v
d
Ω
+
∬
Ω
(
d
1
Δ
u
∂
H
∂
u
+
d
2
Δ
v
∂
H
∂
v
)
d
Ω
=
∬
Ω
u
-
u
∗
2
u
u
+
u
∗
1
+
u
2
1
+
u
∗
2
-
b
-
a
v
∗
d
Ω
-
r
k
e
∬
Ω
v
-
v
∗
2
d
Ω
-
d
1
∬
Ω
u
∗
u
2
∂
u
∂
x
2
+
∂
u
∂
y
2
d
Ω
-
d
2
e
∬
Ω
v
∗
v
2
∂
v
∂
x
2
+
∂
v
∂
y
2
d
Ω
=
I
1
-
I
2
-
I
3
-
I
4
.Obviously, I
2
≥
0, I
3
≥
0, and I
4
≥
0.
Consider the following function:(50)
f
(
x
)
=
x
+
u
∗
1
+
x
2
(
x
≥
0
)
.Then, f
′
(
x
)
=
(
-
x
2
-
2
u
∗
x
+
1
)
/
(
1
+
x
2
)
2. Because lim
x
→
0
⁡
f
(
x
)
=
u
∗
,
lim
x
→
+
∞
⁡
f
(
x
)
=
0, then(51)
f
(
x
)
≤
f
(
-
u
∗
+
u
∗
2
+
1
)
=
u
∗
2
+
1
1
+
-
u
∗
+
u
∗
2
+
1
2
,
I
1
=
∬
Ω
u
-
u
∗
2
u
u
+
u
∗
(
1
+
u
2
)
(
1
+
u
∗
2
)
-
b
-
a
v
∗
d
Ω
≤
∬
Ω
u
-
u
∗
2
u
-
b
-
a
v
∗
+
1
1
+
u
∗
2
1
1
+
u
∗
2
+
1
-
u
∗
2
m
m
m
m
m
m
m
m
m
c
m
·
1
1
+
u
∗
2
+
1
-
u
∗
2
d
Ω
.
Therefore, if-
b
-
a
v
∗
+
(
1
/
1
+
u
∗
2
)
(
1
/
(
1
+
(
1
+
u
∗
2
-
u
∗
)
2
)
)
<
0, I
1
<
0, which is equivalent to d
H
1
(
t
)
/
d
t
<
0.
Thus,H
1
(
t
) satisfies Lyapunov’s asymptotic stability theorem, and the positive equilibrium (
u
∗
,
v
∗
) of system (2) is globally asymptotically stable. This completes the proof.Remark 14.
It is easy to verify that if condition (47) holds, then condition (44) is true.


## 3.2. Turing Instability

To analyze the spatial system and how a small heterogeneous perturbation of a homogeneous steady state develops over a long time period, small space- and time-dependent perturbations forE
∗ of system (2) will be considered:(52)
u
=
u
∗
+
ε
exp
⁡
(
K
z
i
+
λ
t
)
,
v
=
v
∗
+
τ
exp
⁡
K
z
i
+
λ
t
,where ε, τ are small enough and K is the wave number. Substituting (52) into (2) yields its characteristic matrix:(53)
A
=
J
11
-
d
1
K
2
J
12
J
21
J
22
-
d
2
K
2
.Correspondingly, λ
1 and λ
2 are the roots of the following characteristic equation:(54)
λ
2
-
t
r
K
λ
+
Δ
K
=
0
,
t
r
K
=
J
11
+
J
22
-
(
d
1
+
d
2
)
K
2, Δ
K
=
d
1
d
2
K
4
-
(
d
1
J
22
+
d
2
J
11
)
K
2
+
J
11
J
22
-
J
12
J
21.The roots of (54) yield the dispersion relation:(55)
λ
1,2
(
K
)
=
t
r
K
±
t
r
K
2
-
4
Δ
K
2
.As is well known, Turing instability means that the stable equilibrium point is driven to become unstable by the local dynamics and diffusion of species. The stability conditions forE
∗of system (5) are t
r
0
=
J
11
+
J
22
<
0 and Δ
0
=
J
11
J
22
-
J
12
J
21
>
0. It is clear that t
r
K
<
t
r
0
<
0. Therefore, the stability of equilibrium E
∗ of system (2) changes with the sign of Δ
K. It is easy to establish that Δ
K
<
0 for K
1
2
<
K
2
<
K
2
2, where(56)
K
1,2
2
=
d
1
J
22
+
d
2
J
11
d
1
J
22
+
d
2
J
11
2
-
4
d
1
d
2
(
J
11
J
22
-
J
12
J
21
)
m
m
c
∓
d
1
J
22
+
d
2
J
11
2
-
4
d
1
d
2
(
J
11
J
22
-
J
12
J
21
)
·
2
d
1
d
2
-
1
.IfK
1,2
2 has positive values, then the range of instability for E
∗ can be obtained, which is called the Turing space. To illustrate the Turing space, the dispersion relations corresponding to several values of the bifurcation parameter I are plotted in Figure 1. The other parameters come from the actual monitoring data from the Sanyang wetland from 2012 to 2013 and from the literature:(57)
a
=
0.6
,
b
=
0.015
,
e
=
0.28
,
r
=
0.533
,
k
=
1
,
m
=
0.24
,
d
1
=
0.1
,
d
2
=
1.4
.Figure 1
Variation of the dispersion relation of model (2) around E
∗. The yellow line corresponds to I
=
0.0055, the green to I
∗
=
0.006195, and the red to I
=
0.008.In Figure1, the green line corresponds to the critical Turing value, I
∗
=
0.006195. When I
=
0.0055
<
I
∗ (the yellow line in Figure 1), Turing instability occurs, whereas when I
=
0.008
>
I
∗ (the red line in Figure 1), the Turing instability decays.


## 4. Numerical Simulations

On the basis of actual monitoring data from the Sanyang wetland from 2012 to 2013, the following parameter set is considered:I
=
0.008, initial value u
(
0
)
,
v
(
0
)
=
0.1,0.1, t
=
40, and other parameter values as given in (57). To verify the feasibility and correctness of the theoretical results, numerical simulations were performed. It is obvious from Figures 2 and 3 that the interior equilibrium E
∗
=
(
0.02462459569,0.5574801727
) of system (5) is locally asymptotically stable.Figure 2
Time-sequence diagram of nutrients and algae. The red line corresponds to nutrients and the blue line to algae.Figure 3
Phase trajectories of nutrients and algae beginning at different initial levels.To investigate more thoroughly how the input rateI of nutrients flowing into the water affects the spatiotemporal dynamics of system (2), spatial distribution diagrams of the system are generated for various values of I. All these numerical simulations employ the zero-flux boundary conditions and a discrete two-dimensional domain with L
x
=
L
y
=
200. The spatial step size is h
=
1
/
3, and the temporal step size is τ
=
0.01. The initial value of system (2) is placed at the interior equilibrium E
∗
=
(
u
∗
,
v
∗
), and the initial perturbation is 0.0005 space units per time unit. As is well known, in such numerical simulations, the spatial distributions of predator and prey are always of the same type. Therefore, it is necessary to analyze the spatial distribution of only one of these. Here, the spatial distribution of the algae population is considered. Some snapshots have been extracted and are shown in red (blue) corresponding to the high (low) value of algae density v when the value of I increases from 0.005 to 0.006.Figure4 shows the process of pattern formation in system (2) with I
=
0.0055 and other parameters fixed as in (57).Figure 4
Spatial patterns obtained with model (2) for I
=
0.0055. Other parameters are fixed as in (57). Number of iterations: (a) 0; (b) 40000; (c) 55000; (d) 60000; (e) 80000; (f) 300000.
(a)
(b)
(c)
(d)
(e)
(f)From Figure4, it is clear that the algae population pattern takes a long time to settle down, beginning with a uniform state (
u
∗
,
v
∗
)
=
(
0.7722029183,0.7931146159
) (Figure 4(a)). After 40000 iterations, the spatial distribution of the algae population consists mainly of a number of spots and interconnected stripes (Figure 4(b)). The remaining four images consist of blue/red spots on a red/blue background. These are referred to here as “hot spots” and “cold spots” [34]. The hot spots are isolated zones with high algae density and low nutrient concentration, whereas the cold spots represent the opposite case. Figures 4(c) and 4(d) show some isolated “cold spots,” whereas Figures 4(e) and 4(f) contain some isolated “hot spots” and stripes.Figure5 shows two patterns obtained from system (2) at 300000 iterations. When I
=
0.005, the homogeneous state is (
u
∗
,
v
∗
)
=
(
0.7498330425,0.7860636982
), and the spot pattern is made up of “hot spots” around the edge (Figure 5(a)). However, with I
=
0.006, (
u
∗
,
v
∗
)
=
(
0.7860324719,0.7974736497
), and the spot-stripe pattern consists of a certain number of “hot spots” and stripes (Figure 5(b)).Figure 5
Hot-spot patterns (a) and hot spots-stripe patterns (b) obtained with model (2) for (a) I
=
0.005; (b) I
=
0.006, at 300000 iterations. Other parameters are fixed as in (57).
(a)
(b)By comparing Figures5(a), 4(f), and 5(b), it is apparent that the algae population density increases with I within certain limits. Moreover, as the value of I increases, the spatial distribution of the algae population becomes relatively intensive and the distribution relatively wide.


## 5. Discussion and Conclusions

In this paper, a reaction-diffusion nutrient-algae model has been used to investigate the interaction between nutrients and algae mathematically and numerically. Mathematical theoretical work has examined the positivity and boundedness of solutions and the existence and stability of equilibria. Local and global stability analyses of equilibria in the presence of diffusion have also been performed. Mathematical analysis indicated that all solutions of the reduced ODE system are positive and bounded under some certain conditions and that the ODE system always has boundary equilibrium. If-
b
+
2
u
-
/
(
1
+
u
-
2
)
2
<
0 and a
e
u
-
+
r
-
m
<
0, then the nonspatial system and the diffusive system are stable at E
1. When -
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
<
0, the interior equilibrium E
∗
=
(
u
∗
,
v
∗
) is always stable, regardless of whether the system is an ODE or a PDE system. Theorem 12 gives the stability condition of the positive constant solution E
∗
=
(
u
∗
,
v
∗
) of system (2). However, when 0
<
-
b
-
a
v
∗
+
2
u
∗
/
(
1
+
u
∗
2
)
2
<
(
k
/
r
)
a
2
e
u
∗, the stability of E
∗ may vary. In this paper, only the positive constant solution E
∗ of system (2) is discussed. Other questions remain to be answered: whether system (2) has other nonnegative constant solutions or nonconstant positive steady states; if there are positive solutions, what is their exact multiplicity and how stable are any positive solutions.Numerical simulations have indicated that the input rateI of nutrients flowing into the water has an important influence on the density and spatial distribution of the algae population. The spatial distribution of the algae population becomes relatively intensive and the distribution relatively wide as I increases. These results may help to provide a better understanding of the interactions of nutrients and algae and the variations in the spatial distribution of algae over time. This research is also expected to contribute to exploring the mechanisms of eutrophication.



---

*Source: 251937-2015-10-25.xml*