# A Novel Method to Determine the Local Stability of the n-Species Lotka-Volterra System with Multiple Delays

Authors: Xiao-Ping Chen; Hao Liu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The n-species Lotka-Volterra system with discrete delays is considered. The local asymptotic stability of positive equilibrium is investigated based on a contour integral method. The main purpose of this paper is to propose a new and general algorithm to study the local asymptotic stability of the positive equilibrium for then-dimensional Lotka-Volterra system. Some numerical experiments are carried out to show the effectiveness of the proposed method.


---

## Body

## 1. Introduction

In recent years, delay differential equations have been extensively used as models in population dynamics. In this paper, we study the local stability of the positive equilibrium of the following n-species Lotka-Volterra competition system with multiple delays:(1)
d
y
i
t
d
t
=
y
i
t
r
i
-
∑
j
=
1
n
a
i
j
y
j
t
-
τ
i
j
,
i
=
1
,
…
,
n
, where y
i
(
t
)  (
i
=
1
,
…
,
n
) denotes the population size of the ith species at time t and r
i, a
i
j, and τ
i
j  (
i
,
j
=
1
,
…
,
n
) are positive constants.As a rudimentary model on mathematical ecology, Lotka-Volterra system has received a lot of attention and has been researched extensively; see, for example, [1–15] and the references cited therein. A basic model is the two-species competitive system with multiple delays which takes the form of(2)
d
y
1
t
d
t
=
y
1
t
r
1
-
a
11
y
1
t
-
τ
11
-
a
12
y
2
t
-
τ
12
d
y
2
t
d
t
=
y
2
t
r
2
-
a
21
y
1
t
-
τ
21
-
a
22
y
2
t
-
τ
22
. Zhen and Ma [8] studied the local and global stability of the positive equilibrium of (2) by means of Lyapunov functions. There are also some other researches regarding this area; see, for example, [9–11]. Besides, there are some literatures on three-species models. Qin [12] first investigated the stability of positive equilibrium of a three-species Lotka-Volterra model with discrete delays. Ma et al. [13] studied the local stability of positive equilibrium of a three-species Lotka-Volterra model with discrete delays through constructing a proper Lyapunov function. For the general n-species models, Hofbauer and So [14] gave necessary and sufficient conditions for the global asymptotic stability of a positive equilibrium with the restriction τ
i
i
=
0  (
i
=
1
,
…
,
n
). Later on, Campbell [15] extended the results of [14] and studied the local and global stability of the trivial equilibrium of a functional differential equation used to model artificial neural networks with discrete delays. However, in these literatures, we emphasize that the usual approach to study the stability of the equilibrium for system (1) relies on the use of Lyapunov functionals. Generally speaking, it is not easy to construct a Lyapunov functional for an n-dimensional Lotka-Volterra model. Besides, a new Lyapunov functional is required for each model under consideration. As far as we know, there is no good all-purpose algorithm to study the stability of n-species Lotka-Volterra competition system with multiple delays at present. Therefore, it is necessary to find a general algorithm to study the local stability of the positive equilibrium of system (1) which owns more adaptability. In this paper, we propose a novel and general method to study the local asymptotic stability of positive equilibrium of system (1) by the contour integral method. To the best of our knowledge, this paper is the first to deal with the research of analyzing the local stability of Lotka-Volterra system by the contour integral method.The paper is organized as follows. In Section2, we recall an algorithm to find the rightmost eigenvalues of the nonlinear eigenvalue problem (NEP) in some disk and propose a novel algorithm for determining the local asymptotic stability of the positive equilibrium for the general n-species Lotka-Volterra system. Section 3 is devoted to some numerical experiments. Finally, some concluding remarks are made in Section 4.For convenience, we use the following notations:det
⁡
(
A
) and tr
⁡
(
A
) denote the trace and the determinant of a matrix A, respectively; A denotes ∞-norm of a matrix A; I denotes the identity matrix with suitable size; e
i is the ith column vector of I; (
·
)
T denotes the transpose of a vector or a matrix; Γ is the boundary of an open disk Ω and its radius is ρ; ι denotes the imaginary unit of a complex number.


## 2. Local Asymptotical Stability of the Positive Equilibrium

For convenience, in this paper, we impose the following hypothesis for the n-species Lotka-Volterra competition system (1):(
H
1
)
r
i
>
0 and r
i
-
∑
j
∈
J
i
(
a
i
j
/
a
j
j
)
r
j
‍
>
0,  i
=
1
,
…
,
n, where J
i
=
{
1
⩽
j
⩽
n
,
j
≠
i
}.It is well known (see [16]) that if condition (
H
1
) is satisfied, then system (1) has a unique positive equilibrium y
∗
=
(
y
1
∗
,
y
2
∗
,
…
,
y
n
∗
) and we can compute the unique positive equilibrium y
∗ by the following algebraic equation:(3)
∑
j
=
1
n
a
i
j
y
j
∗
=
r
i
,
i
=
1
,
…
,
n
.Let(4)
u
i
t
=
y
i
t
-
y
i
∗
,
i
=
1
,
…
,
n
; then system (1) reduces to the following system:(5)
d
u
i
t
d
t
=
u
i
t
+
y
i
∗
-
∑
j
=
1
n
a
i
j
u
j
t
-
τ
i
j
,
i
=
1
,
…
,
n
.According to expression (5), the variational system of (1) with respect to the positive equilibrium y
∗ is given by the linearized system of (5). That is,(6)
d
u
i
t
d
t
=
-
∑
j
=
1
n
a
i
j
u
j
t
-
τ
i
j
,
i
=
1
,
…
,
n
.Rewrite (6) in the matrix form as follows:(7)
z
˙
t
=
∑
i
,
j
=
1
n
A
i
j
z
t
-
τ
i
j
, where z
(
t
)
=
(
u
1
(
t
)
,
u
2
(
t
)
,
…
,
u
n
(
t
)
)
T and A
i
j
=
-
a
i
j
y
j
∗
e
i
e
j
T.It is easy to see that the characteristic matrix of system (7) is(8)
T
λ
=
λ
I
-
∑
i
,
j
=
1
n
A
i
j
e
-
λ
τ
i
j
.To establish an effective way to determine whether or not the positive equilibrium for system (1) is asymptotically stable, it suffices to study the asymptotic stability of system (7). It is well known that the local stability of trivial solution of system (7) depends on the locations on the complex plane of the characteristic roots of det
⁡
(
T
(
λ
)
)
=
0. And these characteristic roots of system (7) can be computed [17] as the eigenvalues of the following nonlinear eigenvalue problem:(9)
T
λ
x
=
0
, where λ
∈
C,  x
∈
C
n
(
x
≠
0
).As we know, if all the rightmost eigenvalues of NEP (9), which are the eigenvalues whose real parts are the biggest, have (strict) negative real parts, system (7) will be (asymptotically) stable; otherwise, it will be unstable; see, for example, [18].Next, we briefly recall the algorithm to compute the rightmost characteristic roots of system (7); see [18] for details.Assume thatm
(
m
≤
n
) is a positive integer; let f
(
λ
)
=
det
⁡
(
T
(
λ
)
) and the complex moments v
p be(10)
v
p
=
1
2
π
ι
∮
Γ
z
p
f
′
z
f
z
d
z
,
p
=
0,1
,
…
,
2
m
-
1
. By v
p
(
p
=
0,1
,
…
,
2
m
-
1
), two m
×
m Hankel matrices H
m and H
m
< are defined as follows:(11)
H
m
=
v
0
v
1
⋯
v
m
-
1
v
1
v
2
⋯
v
m
⋮
⋮
⋮
v
m
-
1
v
m
⋯
v
2
m
-
2
,
H
m
<
=
v
1
v
2
⋯
v
m
v
2
v
3
⋯
v
m
+
1
⋮
⋮
⋮
v
m
v
m
+
1
⋯
v
2
m
-
1
.Then, we consider that the above integration is evaluated via a trapezoidal rule on the circleΓ. Let N be the number of sample points on the circle Γ and ω
j
=
ρ
exp
⁡
(
2
π
j
ι
/
N
)
(
j
=
0,1
,
…
,
N
-
1
). From the trapezoidal rule, the complex moments v
p can be approximately calculated by the following formula:(12)
v
p
≈
v
^
p
=
1
N
∑
j
=
0
N
-
1
ω
j
p
+
1
f
′
ω
j
f
ω
j
. Based on the Trace-Theorem of Devidenko f
′
(
λ
)
/
f
(
λ
)
=
tr
⁡
(
T
-
1
(
λ
)
T
′
(
λ
)
)
,
v
^
p can be rewritten as(13)
v
^
p
=
1
N
∑
j
=
0
N
-
1
ω
j
p
+
1
tr
⁡
T
-
1
ω
j
T
′
ω
j
.Usingv
^
p
(
p
=
0,1
,
…
,
2
m
-
1
), two Hankel matrices H
^
m and H
^
m
< are formed as follows:(14)
H
m
≈
H
^
m
=
v
^
0
v
^
1
⋯
v
^
m
-
1
v
^
1
v
^
2
⋯
v
^
m
⋮
⋮
⋮
v
^
m
-
1
v
^
m
⋯
v
^
2
m
-
2
,
(15)
H
m
<
≈
H
^
m
<
=
v
^
1
v
^
2
⋯
v
^
m
v
^
2
v
^
3
⋯
v
^
m
+
1
⋮
⋮
⋮
v
^
m
v
^
m
+
1
⋯
v
^
2
m
-
1
.Then, the algorithm to compute the rightmost characteristic roots of system (7) is described as follows.Algorithm 1 (see [18]).
The following is an algorithm for determining the rightmost eigenvalues of NEP (9) in the disk Ω:(1)
Compute the radiusρ
=
∑
i
,
j
=
1
n
A
i
j of the disk Ω whose center is 0.
(2)
Choose the numberN of sample points on the circle Γ which is the boundary of the disk Ω.
(3)
Computeω
j
=
ρ
exp
⁡
(
2
π
j
ι
/
N
)
(
j
=
0,1
,
…
,
N
-
1
).
(4)
Computev
^
p
(
p
=
0,1
,
…
,
2
m
-
1
) by using (13).
(5)
Form twom
×
m Hankel matrices H
^
m and H
^
m
< by using (14) and (15), respectively.
(6)
Compute all the eigenvaluesλ
i
(
i
=
1,2
,
…
,
m
) of the GEP H
^
m
<
x
=
λ
H
^
m
x.
(7)
Compute the smallest singular value ofT
(
λ
i
)
(
i
=
1,2
,
…
,
m
) to determine whether λ
i is the eigenvalue of NEP (9).
(8)
Compute the rightmost eigenvalue of NEP (9) in the disk Ω.Remark 2.
Algorithm1 allows for using a larger positive integer m.Remark 3.
In actual computation,λ
i is considered as the eigenvalue of NEP (9) if the smallest singular value of T
(
λ
i
) is less than 1
0
-
10.Considering what was mentioned above and the stability theory, we can have the following stability criterion for system (1). That is, if all the rightmost eigenvalues of NEP (9) have strict negative real parts, then the positive equilibrium y
∗ of system (1) will be asymptotically stable; otherwise, it will be unstable.Therefore, we have the algorithm to determine the local asymptotic stability of the positive equilibrium for system (1) as follows.Algorithm 4.
The following is an algorithm for determining the local asymptotic stability of the positive equilibrium for system (1):(1)
Compute the unique positive equilibriumy
∗ of system (1) by solving the algebraic equation (3).
(2)
Run Algorithm1 to find the rightmost eigenvalue of NEP (9).
(3)
if the rightmost eigenvalue of NEP (9) has (strict) negative real partsthen
(4)
The positive equilibriumy
∗ of system (1) is asymptotically stable.
(5)
else
(6)
The positive equilibriumy
∗ of system (1) is not stable.
(7)
end if


## 3. Numerical Examples

In this section, we report some numerical examples to study the effectiveness of Algorithm4 for a set of systems (1). All computations are performed under MATLAB (version 2012b). λ
r
m denotes the rightmost eigenvalue obtained by Algorithm 1.Example 1 (see [10]).
We consider the following system:(16)
d
y
1
t
d
t
=
y
1
t
1
-
y
1
t
-
τ
11
-
0.5
y
2
t
-
τ
12
d
y
2
t
d
t
=
y
2
t
1
-
0.5
y
1
t
-
τ
21
-
y
2
t
-
τ
22
, which has a unique positive equilibrium (
2
/
3,2
/
3
).For the case when allτ
i
j
(
i
,
j
=
1,2
) are equal, we set τ
i
j
=
h. By the numerical simulation shown in Figure 1 in which h
=
1.3, the positive equilibrium of system (16) is asymptotically stable. And from Figure 2 with h
=
1.9, the positive equilibrium (
2
/
3,2
/
3
) loses its stability.Figure 1
Dynamic behavior of the system in Example1 with h
=
1.3.
(a)
(b)Figure 2
Dynamic behavior of the system in Example1 with h
=
1.9.
(a)
(b)By usingρ
=
∑
i
,
j
=
1
n
A
i
j, we obtain the radius ρ
=
1.667 of the disk. Table 1 lists the numerical results computed by Algorithm 4 with the parameters N
=
600, m
=
3, and h
=
1.3, from which we observe that the rightmost eigenvalue is approximately -
0.155934. So, the equilibrium of system (16) is locally asymptotically stable, which coincides with the result obtained in [10]. Besides, Table 1 also lists parts of the numerical results computed by Algorithm 4 with the parameters N
=
600, m
=
4, and h
=
1.9, from which we can find that the rightmost eigenvalues are approximately 0.020280
±
0
.
839447
ι. Therefore, the positive equilibrium of system (16) is unstable. The numerical results obtained by Algorithm 4 are in good agreement with the theoretical ones in [10].Table 1
Numerical results in Example1.
Case
h
λ
r
m
1
h
=
1.3
-
0.155934
2
h
=
1.9
0.020280
±
0.839447
ιExample 2 (see [19]).
In this example, we consider the following system:(17)
d
y
1
t
d
t
=
y
1
t
1
-
y
1
t
-
τ
11
-
0.8
y
2
t
-
τ
12
d
y
2
t
d
t
=
y
2
t
1
-
0.3
y
1
t
-
τ
21
-
y
2
t
-
τ
22
, which has a unique positive equilibrium 0.2632,0.9211.From Figure3, we see that positive equilibrium of system (17) is asymptotically stable when τ
11
=
τ
22
=
1.6, τ
12
=
0.6, and τ
21
=
1.3.Figure 3
Dynamic behavior of the system in Example2 with τ
11
=
τ
22
=
1.6, τ
12
=
0.6, and τ
21
=
1.3.
(a)
(b)We can get the radiusρ
=
1.0001 of the disk by using ρ
=
∑
i
,
j
=
1
n
A
i
j. The numerical results computed by Algorithm 4 are shown in Table 2 with the parameters N
=
300 and m
=
4. From Table 2, we observe that the rightmost eigenvalues of system (17) are approximately -
0.048663
±
1
.
005381
ι when τ
11
=
τ
22
=
1.6, τ
12
=
0.6, and τ
21
=
1.3. Hence, the equilibrium is asymptotically stable. Besides, when τ
11
=
τ
22
=
1.9, τ
12
=
0.6, and τ
21
=
1.3, the rightmost eigenvalues of system (17) are approximately 0.010671
±
0
.
885642
ι. Therefore, the positive equilibrium is unstable in this condition. The numerical results obtained by Algorithm 4 are in good agreement with the results of numerical simulation shown in Figures 3 and 4.Table 2
Numerical results in Example2.
Case
τ
i
j
(
i
,
j
=
1,2
)
λ
r
m
1
τ
11
=
τ
22
=
1.6, τ
12
=
0.6, τ
21
=
1.3
-
0.048663
±
1.005381
ι
2
τ
11
=
τ
22
=
1.9, τ
12
=
0.6, τ
21
=
1.3
0.010671
±
0.885642
ιFigure 4
Dynamic behavior of the system in Example2 with τ
11
=
τ
22
=
1.9, τ
12
=
0.6, and τ
21
=
1.3.
(a)
(b)Example 3 (see [20]).
We consider the following system:(18)
d
y
1
t
d
t
=
y
1
t
1
-
y
1
t
-
τ
11
-
0.6
y
2
t
-
τ
12
-
0.6
y
3
t
-
τ
13
d
y
2
t
d
t
=
y
2
t
1
-
0.6
y
1
t
-
τ
21
-
y
2
t
-
τ
22
-
0.6
y
3
t
-
τ
23
d
y
3
t
d
t
=
y
3
t
1
-
0.6
y
1
t
-
τ
31
-
0.6
y
2
t
-
τ
32
-
y
3
t
-
τ
33
, which has a unique positive equilibrium (
5
/
11,5
/
11,5
/
11
).We obtain the radiusρ
=
3 of the disk according to ρ
=
∑
i
,
j
=
1
n
A
i
j. Table 3 lists the numerical results computed by Algorithm 4 with the parameters N
=
300, m
=
6, τ
11
=
3.4, and τ
12
=
τ
13
=
τ
21
=
τ
22
=
τ
23
=
τ
31
=
τ
32
=
τ
33
=
0, from which we observe that the rightmost eigenvalues are approximately -
0.006065
±
0
.
349566
ι. Therefore, the unique positive equilibrium of system (18) is stable, which corresponds to the theoretical ones in [20]. Besides, Algorithm 4 can also determine the local asymptotic stability of the positive equilibrium for system (18) with τ
12
=
τ
13
=
τ
21
=
τ
22
=
τ
23
=
τ
31
=
τ
32
=
τ
33
=
2 (see Table 3). In this condition, the rightmost eigenvalues of system (18) are approximately 0.068341
±
0
.
788396
ι. Therefore, the positive equilibrium for system (18) is not stable. The numerical results obtained by Algorithm 4 are in good agreement with the result of numerical simulation shown in Figure 5.Table 3
Numerical results in Example3.
Case
τ
i
j
(
i
,
j
=
1,2
,
3
)
λ
r
m
1
τ
11
=
3.4, τ
12
=
τ
13
=
τ
21
=
τ
22 = τ
23
=
τ
31
=
τ
32 = τ
33
=
0
-
0.006065
±
0.349566
ι
2
τ
11
=
3.4, τ
12
=
τ
13
=
τ
21 = τ
22
=
τ
23
=
τ
31 = τ
32
=
τ
33
=
2
0.068341
±
0.788396
ιFigure 5
Dynamic behavior of the system in Example3 with τ
11
=
3.4 and τ
12
=
τ
13
=
τ
21
=
τ
22
=
τ
23
=
τ
31
=
τ
32
=
τ
33
=
2.
(a)
(b)
(c)


## 4. Conclusion

In this paper, we have proposed a novel method based on the contour integral method to determine the local stability of positive equilibrium for system (1). Our numerical results show that the proposed method is extremely convenient and useful for the stability analysis of the n-species Lotka-Volterra system with discrete delays.



---

*Source: 7920482-2016-12-13.xml*