# Numerical Simulation of a Class of Three-Dimensional Kolmogorov Model with Chaotic Dynamic Behavior by Using Barycentric Interpolation Collocation Method

Authors: Mingjing Du; Junmei Li; Yulan Wang; Wei Zhang
Journal: Complexity
 (2019)

---

## Abstract

This paper numerically simulates three-dimensional Kolmogorov model with chaotic dynamic behavior by barycentric Lagrange interpolation collocation method. Some numerical examples are studied for finding some new chaotic behaviors and demonstrating some existing chaotic dynamic behaviors of the Kolmogorov model. Results obtained by the present method indicate that the method has merits of small operations and good numerical stability.


---

## Body

## 1. Introduction

J.P. Berrut [1, 2] introduced barycentric Lagrange interpolation and studied its numerical stability and convergence. Barycentric Lagrange interpolation is unconditionally stable at the Chebyshev points. S.P. Li, Z.Q. Wang [3, 4] gave some algorithms of barycentric Lagrange interpolation collocation method (BLICM). Some authors [3, 5–8] solved all sorts of equations and showed the BLICM has merits of small operations and high precision (see [3, 4, 9]). This paper numerically simulates some three-dimensional Kolmogorov models with chaotic dynamic behavior. The purpose of this paper is to find some new chaotic behaviors and verify the existing chaotic dynamic behaviors by the BLICM.Three-dimensional Kolmogorov models comprise a significant class of ecological models that are used widely in ecology to represent the dynamic behavior of prey and predators, which are expressed in the following form:(1)
d
x
d
t
=
x
f
1
x
,
y
,
z
,
d
y
d
t
=
y
f
2
x
,
y
,
z
,
d
z
d
t
=
z
f
3
x
,
y
,
z
,
x
,
y
,
z
≥
0
,
0
≤
t
≤
T
,where x
,
y
,
z represents the population density of the species and f
i
(
x
,
y
,
z
) represents the per capita growth rate of the x
,
y
,
z species.In ecology, the most frequently used model is the Lotka-Volterra system; that is, each per capita growth functionf
i is affine and chosen as the logistic growth. In this circumstance, model (1) reads as(2)
d
x
d
t
=
x
r
1
-
a
11
x
-
a
12
y
-
a
13
z
,
d
y
d
t
=
y
r
2
-
a
21
x
-
a
22
y
-
a
23
z
,
d
z
d
t
=
z
r
3
-
a
31
x
-
a
32
y
-
a
33
z
.
x
,
y
,
z
≥
0
,
0
≤
t
≤
T
,Model (2) is a totally competitive system if all parameters r
i
,
a
i
j are positive.There are several famous functional responses in the Kolmogorov model which are referred to as Holling type I, type II, type III, type IV, Monod-Haldane type, Hassel-Verley type, Beddington-DeAngelis type functional response, etc.In this paper, we consider the following three-dimensional Kolmogorov model with functional response:(3)
d
x
d
t
=
x
r
1
-
a
11
x
-
a
12
y
-
a
13
z
+
h
1
x
,
y
,
z
,
d
y
d
t
=
y
r
2
-
a
21
x
-
a
22
y
-
a
23
z
+
h
2
x
,
y
,
z
,
d
z
d
t
=
z
r
3
-
a
31
x
-
a
32
y
-
a
33
z
+
h
3
x
,
y
,
z
,
x
,
y
,
z
≥
0
,
0
≤
t
≤
T
,with the following initial condition:(4)
x
0
=
c
1
,
y
0
=
c
2
,
z
0
=
c
3
,where x
,
y
,
z represent the population density of the species and h
i
(
x
,
y
,
z
) is the known functional response. r
i
,
c
i
,
a
i
j are known constants. x
,
y
,
z are unknown functions of time t.Many researchers [10–20] studied the dynamics of three dimensional Kolmogorov model with different type of functional response in theory. They found some chaotic dynamics [13, 18–23] of the three-dimensional Kolmogorov model. Chaos and hyperchaos exist in many natural processes and are one of the main contents of nonlinear science research. Although many kinds of numerical methods of the Kolmogorov model have been announced, simple and efficient methods have always been the direction that scholars strive to pursue. This paper suggests the BLICM to solve the three-dimensional Kolmogorov model. Model (3) is adopted as an example to elucidate the solution process.


## 2. Barycentric Lagrange Interpolation Collocation Method

First of all, we discrete computational interval[
0
,
T
] by Chebyshev points into 0
=
t
1
<
t
2
<
⋯
<
t
M
=
T and construct following linear iterative format of model (3):(5)
d
x
n
d
t
=
r
1
x
n
-
a
11
x
n
+
a
12
y
n
+
a
13
z
n
x
n
-
1
+
h
1
x
n
-
1
,
y
n
-
1
,
z
n
-
1
,
d
y
n
d
t
=
r
2
y
n
-
a
21
x
n
+
a
22
y
n
+
a
23
z
n
y
n
-
1
+
h
2
x
n
-
1
,
y
n
-
1
,
z
n
-
1
,
d
z
n
d
t
=
r
3
z
n
-
a
31
x
n
+
a
32
y
n
+
a
33
z
n
z
n
-
1
+
h
3
x
n
-
1
,
y
n
-
1
,
z
n
-
1
.The format (5) is convergent; then x
n
(
t
)
→
x
(
t
), y
n
(
t
)
→
y
(
t
), z
n
(
t
)
→
z
(
t
)
,
(
n
→
∞
).Next, we transform format (5) into the following linear algebraic equations.(6)
D
-
r
1
I
+
a
11
d
i
a
g
x
n
-
1
a
12
d
i
a
g
x
n
-
1
a
13
d
i
a
g
x
n
-
1
a
21
d
i
a
g
y
n
-
1
D
-
r
2
I
+
a
22
d
i
a
g
y
n
-
1
a
23
d
i
a
g
y
n
-
1
a
31
d
i
a
g
z
n
-
1
a
32
d
i
a
g
z
n
-
1
D
-
r
3
I
+
a
33
d
i
a
g
z
n
-
1
x
n
y
n
z
n
=
h
1
x
n
-
1
,
y
n
-
1
,
z
n
-
1
h
2
x
n
-
1
,
y
n
-
1
,
z
n
-
1
h
3
x
n
-
1
,
y
n
-
1
,
z
n
-
1where D
=
(
ξ
j
′
(
t
i
)
)
i
,
j
=
1,2
,
⋯
M is M order matrix. I is M order unit matrix and d
i
a
g is a symbol of diagonal matrix composed of vectors. ξ
j
(
t
)
=
ω
j
/
t
-
t
j
/
∑
k
=
1
M
ω
k
/
t
-
t
k is, respectively, barycentric interpolation primary function and ω
j
=
1
/
∏
i
=
1
,
j
≠
k
M
(
t
i
-
t
j
) is center of gravity Lagrange interpolation weight. The vector (7)
x
n
,
y
n
,
z
n
=
x
n
t
1
,
x
n
t
2
,
…
,
x
n
t
M
,
y
n
t
1
,
y
n
t
2
,
…
,
y
n
t
M
,
z
n
t
1
,
z
n
t
2
,
…
,
z
n
t
M
.The vector(8)
h
1
x
n
-
1
,
y
n
-
1
,
z
n
-
1
,
h
2
x
n
-
1
,
y
n
-
1
,
z
n
-
1
,
h
3
x
n
-
1
,
y
n
-
1
,
z
n
-
1
=
h
1
x
n
-
1
t
1
,
y
n
-
1
t
1
,
z
n
-
1
t
1
,
…
,
h
1
x
n
-
1
t
M
,
y
n
-
1
t
M
,
z
n
-
1
t
M
,
h
2
x
n
-
1
t
1
,
y
n
-
1
t
1
,
z
n
-
1
t
1
,
…
,
h
2
x
n
-
1
t
M
,
y
n
-
1
t
M
,
z
n
-
1
t
M
,
h
3
x
n
-
1
t
1
,
y
n
-
1
t
1
,
z
n
-
1
t
1
,
…
,
h
3
x
n
-
1
t
M
,
y
n
-
1
t
M
,
z
n
-
1
t
M
.The last and the first line of equations (6) are replaced separately by the equation of the initial condition (9) in turn.(9)
∑
i
=
1
M
ξ
i
0
x
n
t
i
=
c
1
,
∑
i
=
1
M
ξ
i
0
y
n
t
i
=
c
2
,
∑
i
=
1
M
ξ
i
0
z
n
t
i
=
c
3
.So, we can get a numerical solution of (3) and (4).


## 3. Numerical Experiments

In this section, some numerical examples are studied to demonstrate the accuracy of the present method. The examples are computed using MatlabR2017a. In numerical experiments, the number of nodesM
=
40. The accuracy of iteration control ε
=
1
0
-
10 and the initial iteration value x
0
=
y
0
=
z
0
=
0
;
x
1
=
y
1
=
z
1
=
t.Experiment 1.
We consider the following three-species food chain model [18]:(10)
x
˙
=
r
x
1
-
x
50
-
0.1
x
y
1
+
0.1
x
,
y
˙
=
-
y
+
0.05
x
y
-
0.1
y
z
1
+
0.2
y
,
z
˙
=
-
z
+
0.03
y
z
,where r is real parameters, which satisfy the initial condition x
(
0
)
=
1
,
y
(
0
)
=
2
,
z
(
0
)
=
1
. Results of Experiment 1 are given in Figures 1-2.
Figure1 is obtained by using the current method with r
=
1. Among them, (
a
1
) is the time series plot; (
b
1
) is the phase diagram of x; (
c
1
) is the phase diagram of y; (
d
1
) is the graph projected on (
x
,
y
)-plane; (
e
1
) is the three-dimensional space graph. Figure 2 is obtained by using the current method with r
=
4.2. (
a
2
) is the time series plot; (
b
2
) is the phase diagram of y; (
c
2
) is the three-dimensional space graph; (
d
2
) is the graph projected on (
x
,
y
)-plane; (
e
2
) is the the graph projected on (
x
,
z
)-plane; (
f
) is the graph projected on (
y
,
z
)-plane.
Figure1 gets some new chaotic behaviors. Figure 2 verifies the existing chaotic dynamic behaviors [18]. Our study suggests that model (10) will go chaotic when the rate of the self-reproduction of the prey is large. Our numerical results are in good agreement with the theory [18].Figure 1
The three species food chain model for Experiment1 with r
=
1: (
a
1
) time series plot; (
b
1
) phase diagram of x; (
c
1
) phase diagram of y; (
d
1
) projected on the (
x
,
y
)-plane; (
e
1
) on the three-dimensional space.Figure 2
The three-species food chain model for Experiment1 with r
=
4.2: (
a
2
) time series plot; (
b
2
) phase diagram of y; (
c
2
) on the three-dimensional space; (
d
2
) projected on the (
x
,
y
)-plane; (
e
2
) projected on the (
x
,
z
)-plane; (
f
) projected on the (
y
,
z
)-plane.Experiment 2.
We consider the following tritrophic food chain model [19]:(11)
x
˙
=
x
1
-
x
-
y
1
+
1.4
x
-
z
1
+
5
x
+
8
y
,
y
˙
=
y
-
0.16
+
x
1
+
1.4
x
-
0.1
z
1
+
5
x
+
8
y
,
z
˙
=
0.1
z
2
-
0.5
z
2
1
+
8
x
+
8
y
,
In comparison with [19], we get better numerical results; results of Experiment 2 are given in Figures 3-4. Our numerical results are in good agreement with the theory [19].
Figures3 and 4 are obtained by using the current method with the initial conditions (
1
): x
(
0
)
=
1
,
y
(
0
)
=
1
,
z
(
0
)
=
1 and the initial conditions (
2
): x
(
0
)
=
1
,
y
(
0
)
=
0
,
z
(
0
)
=
1, respectively. (
a
1
) and (
a
2
) are the time series plot of x
,
y
,
z; (
b
1
) and (
b
2
) are the phase diagram of z; (
c
1
) and (
c
2
) are the three-dimensional space graph; (
d
1
) and (
d
2
) are the graph projected on (
x
,
y
)-plane; (
e
1
) and (
e
2
) are the graph projected on (
x
,
z
)-plane; (
f
1
) and (
f
2
) are the graph projected on (
y
,
z
)-plane. It can be seen that the change of initial conditions leads to significant changes in the time series diagrams of x
,
y, and z.Figure 3
A tritrophic food web system for Experiment2 with the initial condition (
1
): (
a
1
) time series plot of x
,
y
,
z; (
b
1
) phase diagram of z; (
c
1
) on the three-dimensional space; (
d
1
) projected on the (
x
,
y
)-plane; (
e
1
) projected on the (
x
,
z
)-plane; (
f
1
) projected on the (
y
,
z
)-plane.Figure 4
A tritrophic food web system for Experiment2 with the initial condition (
2
): (
a
2
) time series plot of x
,
y
,
z; (
b
2
) phase diagram of z; (
c
2
) on the three-dimensional space; (
d
2
) projected on the (
x
,
y
)-plane; (
e
2
) projected on the (
x
,
z
)-plane; (
f
2
) projected on the (
y
,
z
)-plane.Experiment 3.
We consider the following tritrophic food chain model [20]:(12)
x
˙
=
x
1
-
x
30
-
w
y
x
2
/
a
+
x
+
60
,
y
˙
=
y
-
0.05
+
0.6
x
x
2
/
a
+
x
+
60
-
0.2
z
y
+
50
+
0.21
z
,
z
˙
=
z
-
0.06
+
0.12
y
y
+
50
+
0.21
z
,
wherea is real parameters, which satisfy the initial condition x
(
0
)
=
5
,
y
(
0
)
=
4
,
z
(
0
)
=
2
. Our numerical simulations and experimental results of [20] are confirmed. Results of Experiment 3 are given in Figure 5.
Figure5 is obtained by using the current method. (
a
) is time series plot and (
d
) is the three-dimensional space graph with a
=
7; (
b
) is time series plot and (
e
) is the three-dimensional space graph with a
=
9; (
c
) is time series plot and (
f
) is the three-dimensional space graph with a
=
12. We can see that, from the figure, the fluctuation range of x
,
y
,
z at t
=
900 is very large when a
=
7; the fluctuation range of x
,
y
,
z at t
=
900 is obviously reduced when a
=
9; x
,
y
,
z at t
=
900 almost all tend to be stationary when a
=
12.Figure 5
Time series plots and the three-dimensional space graphs of the tritrophic food chain model for Experiment3 with different parameters: (
a
) and (
d
)
j
=
7; (
b
) and (
e
)
j
=
9; (
c
) and (
f
)
j
=
12.Experiment 4.
We consider the following model with Beddington-DeAngelis-type functional response [21]:(13)
x
˙
=
r
x
1
-
x
k
-
a
x
z
1
+
b
x
+
c
z
,
y
˙
=
k
1
a
x
z
1
+
b
x
+
c
z
-
α
y
-
d
1
y
,
z
˙
=
α
y
-
d
2
z
,where r
,
k
,
a
,
b
,
c
,
α
,
d
1
,
d
2, and k
1 are real parameters, which satisfy the initial condition x
(
0
)
=
1.6
,
y
(
0
)
=
0.5
,
z
(
0
)
=
0.5
. We choose the parameters r
=
1.5
,
k
=
1.6
,
a
=
1.4
,
b
=
1
,
c
=
0.1
,
α
=
0.09
,
d
1
=
0.5
,
d
2
=
0.05. Results of Experiment 4 are given in Figures 6-7.
Figure6 is time series plots and the three-dimensional space graphs of a predator-prey model with Beddington-DeAngelis type functional response for Experiment 4 with different parameters obtained by using the current method. (
a
1
) is time series plot and (
b
1
) is the three-dimensional space graph with k
1
=
0.958; (
a
2
) is time series plot and (
b
2
) is the three-dimensional space graph with k
1
=
0.989; (
a
3
) is time series plot and (
b
3
) is the three-dimensional space graph with k
1
=
1.089. Figure 7 is obtained by using the current method with k
1
=
0.958 and k
1
=
0.989, respectively. (
c
1
) and (
d
1
) are the graph projected on (
x
,
y
)-plane; (
c
2
) and (
d
2
) are the graph projected on (
x
,
z
)-plane; (
c
3
) and (
d
3
) are the graph projected on (
y
,
z
)-plane.Figure 6
Time series plots and the three-dimensional space graphs of Experiment4 with different parameters: (
a
1
) and (
b
1
)
k
1
=
0.958; (
a
2
) and (
b
2
)
k
1
=
0.989; (
a
3
) and (
b
3
)
k
1
=
1.089.Figure 7
The results of Experiment4: (
c
i
) with k
1
=
0.958 and (
d
i
) with k
1
=
0.989
(
i
=
1,2
,
3
), (
c
1
) and (
d
1
) projected on the (
x
,
y
)-plane; (
c
2
) and (
d
2
) projected on the (
x
,
z
)-plane; (
c
3
) and (
d
3
) projected on the (
y
,
z
)-plane.Experiment 5.
We consider the following model with Monod-Haldane type response function [22]:(14)
x
˙
=
ρ
1
x
1
-
x
k
1
-
β
x
z
ϕ
+
γ
x
2
,
y
˙
=
c
β
x
z
ϕ
+
γ
x
2
-
α
+
δ
1
y
,
z
˙
=
α
y
-
δ
2
z
,where ρ
1
,
k
1
,
β
,
γ
,
c
,
ϕ
,
α
,
δ
1, and δ
2 are real parameters, which satisfy the initial condition x
(
0
)
=
1
,
y
(
0
)
=
1
,
z
(
0
)
=
1
. We choose the parameter ρ
1
=
1.6
,
k
1
=
1.72
,
β
=
1.4
,
γ
=
0.01
,
α
=
0.09
,
ϕ
=
0.7
,
δ
1
=
0.5
,
δ
2
=
0.05.
Our numerical simulation results are in good agreement with the theory [22]. We give some new chaotic behaviors. Results of Experiment 5 are given in Figures 8-9.
Figure8 is time series plots and the three-dimensional space graphs of a predator-prey system with Monod-Haldane type response function for Experiment 5 with different parameters obtained by using the current method. (
a
1
) is time series plot and (
b
1
) is the three-dimensional space graph with c
=
1.4; (
a
2
) is time series plot and (
b
2
) is the three-dimensional space graph with c
=
1.575; (
a
3
) is time series plot and (
b
3
) is the three-dimensional space graph with c
=
1.62. Figure 9 is obtained by using the current method with c
=
1.4 and c
=
1.575, respectively. (
c
1
) and (
d
1
) are the graph projected on (
x
,
y
)-plane; (
c
2
) and (
d
2
) are the graph projected on (
x
,
z
)-plane; (
c
3
) and (
d
3
) are the graph projected on (
y
,
z
)-plane.Figure 8
Time series plots and the three-dimensional space graphs of Experiment5 with different parameters: (
a
1
) and (
b
1
)
c
=
1.4; (
a
2
) and (
b
2
)
c
=
1.575; (
a
3
) and (
b
3
)
c
=
1.62.Figure 9
The results of Experiment5: (
c
i
) with c
=
1.4 and (
d
i
) with c
=
1.575
(
i
=
1,2
,
3
), (
c
1
) and (
d
1
) projected on the (
x
,
y
)-plane; (
c
2
) and (
d
2
) projected on the (
x
,
z
)-plane; (
c
3
) and (
d
3
) projected on the (
y
,
z
)-plane.Experiment 6.
We consider the following model with a Watt-type functional response [23]:(15)
x
˙
=
x
b
1
t
-
d
1
t
x
-
a
1
t
z
1
-
e
x
p
-
c
x
z
m
+
D
1
t
y
-
x
,
y
˙
=
y
b
2
t
-
d
2
t
y
+
D
2
t
x
-
y
,
z
˙
=
z
b
3
t
-
d
3
t
z
+
a
2
t
z
1
-
e
x
p
-
c
x
z
m
,
Our numerical simulation results are in good agreement with the theory [23]. Compared with [23], we get better numerical results and give some new chaotic behaviors. Results of Experiment 6 are given in Figures 10–15. We choose the parameter; see Tables 1-2.
Figures10–12 are obtained by using the current method. (
a
1
), (
a
2
), and (
a
3
) are the time series plot of x
,
y
,
z; (
b
1
), (
b
2
), and (
b
3
) are the phase diagram of x; (
c
1
), (
c
2
), and (
c
3
) are the phase diagram of y; (
d
1
), (
d
2
), and (
d
3
) are the graph projected on (
x
,
y
)-plane; (
e
1
), (
e
2
), and (
e
3
) are the three-dimensional space graph. Figures 13–15 are obtained by using the current method. (
a
4
), (
a
5
), and (
a
6
) are the time series plot of x
,
y
,
z; (
b
4
), (
b
5
), and (
b
6
) are the phase diagram of x; (
c
4
), (
c
5
), and (
c
6
) are the phase diagram of y; (
d
4
), (
d
5
), and (
d
6
) are the graph projected on (
x
,
y
)-plane; (
e
4
), (
e
5
), and (
e
6
) are the three-dimensional space graph.Table 1
Parameters used in Experiment6.
Figure
b
1
(
t
)
b
2
(
t
)
b
3
(
t
)
d
1
(
t
)
d
2
(
t
)
d
3
(
t
)
a
1
(
t
)
a
2
(
t
)
10, 11, 12
0.85
+
s
i
n
(
2
π
t
)
0.75
+
s
i
n
(
2
π
t
)
0.8
+
s
i
n
(
2
π
t
)
0.35
+
c
o
s
(
2
π
t
)
0.5
0.4
0.6
0.45
+
s
i
n
(
2
π
t
)
13, 14, 15
0.45
+
s
i
n
(
2
π
t
)
0.75
+
s
i
n
(
2
π
t
)
0.8
+
s
i
n
(
2
π
t
)
0.5
+
c
o
s
(
2
π
t
)
0.35
0.4
0.26
0.3
+
s
i
n
(
2
π
t
)Table 2
Parameters used in Experiment6.
Figure
c
m
D
1
(
t
)
D
2
(
t
)
c
1
c
2
c
3
10
1
0.5
0.2
+
c
o
s
(
2
π
t
)
0.3
+
s
i
n
(
2
π
t
)
0
1
1
11
1
0.5
0.2
+
c
o
s
(
2
π
t
)
0.3
+
s
i
n
(
2
π
t
)
1
0
1
12
1
0.5
0.2
+
c
o
s
(
2
π
t
)
0.3
+
s
i
n
(
2
π
t
)
1
1
1
13
1
0.35
0.4
+
c
o
s
(
2
π
t
)
0.1
+
s
i
n
(
2
π
t
)
0
1
1
14
1
0.35
0.4
+
c
o
s
(
2
π
t
)
0.1
+
s
i
n
(
2
π
t
)
1
0
1
15
1
0.35
0.4
+
c
o
s
(
2
π
t
)
0.1
+
s
i
n
(
2
π
t
)
1
1
1Figure 10
The results of Experiment6: (
a
1
)  time series plot of x
,
y
,
z; (
b
1
) phase diagram of x; (
c
1
) phase diagram of y; (
d
1
) projected on the (
x
,
y
)-plane; (
e
1
) on the three-dimensional space.Figure 11
The results of Experiment6: (
a
2
) time series plot of x
,
y
,
z; (
b
2
) phase diagram of x; (
c
2
) phase diagram of y; (
d
2
) projected on the (
x
,
y
)-plane; (
e
2
) on the three-dimensional space.Figure 12
The results of Experiment6: (
a
3
) time series plot of x
,
y
,
z; (
b
3
) phase diagram of x; (
c
3
) phase diagram of y; (
d
3
) projected on the (
x
,
y
)-plane; (
e
3
) on the three-dimensional space.Figure 13
The results of Experiment6: (
a
4
) time series plot of x
,
y
,
z; (
b
4
) phase diagram of x; (
c
4
) phase diagram of y; (
d
4
) projected on the (
x
,
y
)-plane; (
e
4
) on the three-dimensional space.Figure 14
The results of Experiment6: (
a
5
) time series plot of x
,
y
,
z; (
b
5
) phase diagram of x; (
c
5
) phase diagram of y; (
d
5
) projected on the (
x
,
y
)-plane; (
e
5
) on the three-dimensional space.Figure 15
The results of Experiment6: (
a
6
) time series plot of x
,
y
,
z; (
b
6
) phase diagram of x; c
6 phase diagram of y; (
d
6
) projected on the (
x
,
y
)-plane; (
e
6
) on the three-dimensional space.


## 4. Conclusions and Remarks

In this paper, the three dimensional Kolmogorov model was solved by using the barycentric Lagrange interpolation collocation method. These numerical experiments illustrate that numerical results of the present method are the same as experimental results.All computations are performed by the MatlabR2007b software packages.



---

*Source: 3426974-2019-04-09.xml*