# Numerical Simulation of the Lorenz-Type Chaotic System Using Barycentric Lagrange Interpolation Collocation Method

Authors: Jun-Mei Li; Yu-Lan Wang; Wei Zhang
Journal: Advances in Mathematical Physics
 (2019)

---

## Abstract

Although some numerical methods of the Lorenz system have been announced, simple and efficient methods have always been the direction that scholars strive to pursue. Based on this problem, this paper introduces a novel numerical method to solve the Lorenz-type chaotic system which is based on barycentric Lagrange interpolation collocation method (BLICM). The system(
1
) is adopted as an example to elucidate the solution process. Numerical simulations are used to verify the effectiveness of the present method.


---

## Body

## 1. Introduction

In 1963, Edward Lorenz developed a simplified mathematical model for atmospheric convection. The model is a system of three ordinary differential equations now known as the Lorenz equations (see [1–3]):(1)
d
x
d
t
=
a
y
-
x
,
d
y
d
t
=
x
c
-
x
-
y
,
d
z
d
t
=
-
b
z
+
x
y
,
t
∈
0
,
Twith the initial conditions(2)
x
0
=
c
1
,
y
0
=
c
2
,
z
0
=
c
3
,wherex
(
t
) is proportional to the rate of convection, y
(
t
) to the horizontal temperature variation, and z
(
t
) to the vertical temperature variation. The constants a
,
c
,
b are system parameters proportional to the Prandtl number, Rayleigh number, and certain physical dimensions of the layer itself.As chaos theory progresses, many new Lorenz-type systems [4–6] have been proposed, specially Lorenz hyperchaotic systems [7–10]. The Lorenz system is widely used in electric circuits, chemical reactions, and forward osmosis. Although some numerical methods of the Lorenz system have been announced, simple and efficient methods have always been the direction that scholars strive to pursue.With the development of numerical analysis, there are some high-precision methods, such as variational iteration method [11–13], BLICM [14–22], and so on [23]. J.P. Berrut [24–26] introduced barycentric Lagrange interpolation, [27, 28] studied numerical stability of barycentric Lagrange interpolation, and [15, 16] give algorithm of BLICM. Some authors [14, 17–22] have used BLICM to solve all sorts of problems and show the BLICM is a high precision numerical method. This paper suggests the BLICM to solve the Lorenz system. The system (1) is adopted as an example to elucidate the solution process.


## 2. The Numerical Solution of the System (1)

First of all, we give initial functionx
0
(
t
)
,
y
0
(
t
)
,
z
0
(
t
) and construct following linear iterative format of system (1)(3)
d
x
n
d
t
+
a
x
n
t
-
a
y
n
t
=
0
,
d
y
n
d
t
-
c
x
n
t
+
y
n
t
=
x
n
-
1
t
z
n
-
1
t
,
d
z
n
d
t
+
b
z
n
t
=
x
n
-
1
t
y
n
-
1
t
.
n
=
1,2
,
3
,
…
,Next, we use BLICM to solve (3).Using the barycentric Lagrange interpolation functions [14–16, 24–26], we can get following:(4)
x
n
t
=
∑
j
=
1
M
ξ
j
t
x
n
t
j
,
y
n
t
=
∑
j
=
1
M
ξ
j
t
y
n
t
j
,
z
n
t
=
∑
j
=
1
M
ξ
j
t
z
n
t
j
.
(5)
x
n
′
t
=
∑
j
=
1
M
ξ
j
′
t
x
n
t
j
,
y
n
′
t
=
∑
j
=
1
M
ξ
j
′
t
y
n
t
j
,
z
n
′
t
=
∑
j
=
1
M
ξ
j
′
t
z
n
t
j
.
(6)
∑
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
.whereξ
j
(
t
)
=
ω
j
/
t
-
t
j
/
∑
k
=
1
M
(
ω
k
/
t
-
t
k
) is respectively barycentric interpolation primary function, ω
j
=
1
/
∏
i
=
1
,
j
≠
k
M
(
t
i
-
t
j
) is center of gravity Lagrange interpolation weight, and 0
≤
t
1
<
t
2
<
⋯
<
t
M
≤
T.Substitute formulae (4) and (5) into iterative format (3) and let t
=
t
i
,
(
i
=
1,2
,
⋯
M
). So, linear iterative format (3) can be written in following partitioned matrix form:(7)
D
+
a
I
-
a
I
0
-
c
I
D
+
I
0
0
0
D
+
b
I
x
n
y
n
z
n
=
0
x
n
-
1
z
n
-
1
x
n
-
1
y
n
-
1where I is M order unit matrix, D
=
(
ξ
j
′
(
t
i
)
)
i
,
j
=
1,2
,
⋯
M is M order matrix, and the vector(8)
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(9)
0
,
x
n
-
1
z
n
-
1
,
x
n
-
1
y
n
-
1
=
0
,
⋯
,
0
,
x
n
-
1
t
1
z
n
-
1
t
1
,
x
n
-
1
t
2
z
n
-
1
t
2
,
⋯
,
x
n
-
1
t
M
z
n
-
1
t
M
,
x
n
-
1
t
1
y
n
-
1
t
1
,
x
n
-
1
t
2
y
n
-
1
t
2
,
⋯
,
x
n
-
1
t
M
y
n
-
1
t
M
.The first line of (7) is replaced separately by the equation of initial conditions (6) in turn.So, we can get thatx
n
(
t
j
)
,
y
n
(
t
j
)
,
z
n
(
t
j
)
,
(
j
=
1,2
,
⋯
M
) are approximate solution of (1) and (2).


## 3. Numerical Experiment

In this section, some numerical examples are studied to find some new chaotic behaviors and verify the existing chaotic dynamic behaviors. In Experiments1–5, the accuracy of iteration control is ε
=
1
0
-
10, the initial iteration value x
0
=
y
0
=
z
0
=
0
;
x
1
=
y
1
=
z
1
=
T, and for parameters a
,
b
,
c, and m see Table 1.Table 1
Parameters used in the Experiments1–5.
Fig.
a
b
c
m
1
-1.5
5
1
2
-1.5
5
10
3
-1.5
5
100
4
7
1
5
5
7
2
5
6
7
3
5
8
10
15
1
9
10
15
2
10
10
15
3
12
10
8/3
28
13
0.5
-0.1
1.5
0.12Experiment 1.
We consider the model (1) with a
=
-
1.5
,
b
=
5 and the initial conditions x
(
0
)
=
0
,
y
(
0
)
=
1
,
z
(
0
)
=
0[4].We choose Chebyshev nodes, and the number of nodesM
=
40. Figure 1 is obtained by using the current method with c
=
1. Among them, (a) is the time series plot; (b) is the phase diagram of z; (c) is the three-dimensional space graph; (d) is the graph projected on (
x
,
z
)-plane; (e) is the graph projected on (
y
,
z
)-plane. Figures 2 and 3 are obtained by using the current method at c
=
10 and c
=
100, respectively. We can see that the fluctuation amplitude of x and y increases, while the fluctuation amplitude of z decreases with the increase of ρ. The corresponding graphs b
,
c
,
d, and e also have obvious changes.Figure 1
Lorenz system for Experiment1 at c
=
1: (a) time series plot; (b) phase diagram of z; (c) on the three-dimensional space; (d) projected on the (
x
,
z
)-plane; (e) projected on the (
y
,
z
)-plane.
(a)
(b)
(c)
(d)
(e)Figure 2
Lorenz system for Experiment1 at c
=
10: (a) time series plot; (b) phase diagram of z; (c) on the three-dimensional space; (d) projected on the (
x
,
z
)-plane; (e) projected on the (
y
,
z
)-plane.
(a)
(b)
(c)
(d)
(e)Figure 3
Lorenz system for Experiment1 at c
=
100: (a) time series plot; (b) phase diagram of z; (c) on the three-dimensional space; (d) projected on the (
x
,
z
)-plane; (e) projected on the (
y
,
z
)-plane.
(a)
(b)
(c)
(d)
(e)Experiment 2.
We consider the Lorenz-type system [6](10)
d
x
d
t
=
a
y
-
x
+
y
z
,
d
y
d
t
=
c
x
-
x
z
,
d
z
d
t
=
-
b
z
+
x
yWe choose Chebyshev nodes, the number of nodesM
=
40, and the parameters a
=
7
,
c
=
5 and the initial conditions x
(
0
)
=
0
,
y
(
0
)
=
1
,
z
(
0
)
=
0.Figure4 is strange attractors of chaotic system for Experiment 2 at b
=
1 by using the current method. (a) is the graph projected on (
x
,
y
)-plane; (b) is the graph projected on (
x
,
z
)-plane; (c) is the graph projected on (
y
,
z
)-plane. Figure 7 is time series plots of chaotic system for Experiment 2 at different parameter value b. We can see that the frequency of fluctuations of x
,
y, and z accelerates obviously with the increase of b. When b
=
3, the fluctuations of x
,
y, and z change obviously, and their fluctuations become smaller and smaller and finally stop at a certain value. Figures 5 and 6 are strange attractors of chaotic system for Experiment 2 by using the current method at b
=
2 and b
=
3 respectively.Figure 4
Strange attractors of chaotic system for Experiment2 at b
=
1: (a) (
x
,
y
)-plane; (b) (
x
,
z
)-plane; (c) (
y
,
z
)-plane.
(a)
(b)
(c)Figure 5
Strange attractors of chaotic system for Experiment2 at b
=
2: (a) (
x
,
y
)-plane; (b) (
x
,
z
)-plane; (c) (
y
,
z
)-plane.
(a)
(b)
(c)Figure 6
Strange attractors of chaotic system for Experiment2 at b
=
3: (a) (
x
,
y
)-plane; (b) (
x
,
z
)-plane; (c) (
y
,
z
)-plane.
(a)
(b)
(c)Figure 7
The time series plots of chaotic system for Experiment2: (a) b
=
1; (b) b
=
2; (c) b
=
3.
(a)
(b)
(c)Experiment 3.
We consider the 3D autonomous chaotic Lorenz-type system [7](11)
d
x
d
t
=
a
y
-
x
,
d
y
d
t
=
-
x
z
-
c
y
,
d
z
d
t
=
-
b
+
x
y
,We choose Chebyshev nodes and the number of nodesM
=
60 and the parameters a
=
10
,
b
=
15 and the initial conditions x
(
0
)
=
10
,
y
(
0
)
=
-
0.2
,
z
(
0
)
=
0.75.Figure8 is phase portraits of the 3D chaotic Lorenz type system for Experiment 3 at c
=
1 by using the current method. (a) is the graph projected on (
x
,
y
)-plane; (b) is the graph projected on (
x
,
z
)-plane; (c) is the graph projected on (
y
,
z
)-plane. Figure 11 is time series plots of the 3D chaotic Lorenz type system for Experiment 3 at different parameter value c. We can see that the fluctuation range of y changes obviously with the increase of c. When c
=
1, the fluctuation range of y is -
5 to 15, and when c
=
3, the fluctuation range of y is -
15 to 5. Figures 9 and 10 are phase portraits of the 3D chaotic Lorenz type system for Experiment 3 by using the current method at c
=
2 and c
=
3, respectively.Figure 8
Phase portraits of the 3D chaotic Lorenz type system for Experiment3 at c
=
1: (a) (
x
,
y
)-plane; (b) (
x
,
z
)-plane; (c) (
y
,
z
)-plane.
(a)
(b)
(c)Figure 9
Phase portraits of the 3D chaotic Lorenz type system for Experiment3 at c
=
2: (a) (
x
,
y
)-plane; (b) (
x
,
z
)-plane; (c) (
y
,
z
)-plane.
(a)
(b)
(c)Figure 10
Phase portraits of the 3D chaotic Lorenz type system for Experiment3 at c
=
3: (a) (
x
,
y
)-plane; (b) (
x
,
z
)-plane; (c) (
y
,
z
)-plane.
(a)
(b)
(c)Figure 11
The time series plots of 3D autonomous chaotic Lorenz-type system for Experiment3: (a) c
=
1; (b) c
=
2; (c) c
=
3.
(a)
(b)
(c)Experiment 4.
We consider the Lorenz system [5](12)
d
x
d
t
=
-
a
x
+
a
y
,
d
y
d
t
=
-
a
x
-
y
-
x
z
,
d
z
d
t
=
-
b
z
+
x
y
-
b
c
+
a
,where a
,
b, and c are real parameters, which satisfy the following initial conditions:(13)
x
0
=
1
,
y
0
=
2
,
z
0
=
1We choose Chebyshev nodes and the number of nodesM
=
40. Figure 12 is obtained by using the current method with the parameters a
=
10
,
b
=
8
/
3, and c
=
28. In Figure 12, (a) is the time series plot of x; (b) is the time series plot of z; (c) is the three-dimensional space graph; (d) is the graph projected on (
x
,
y
)-plane; (e) is the graph projected on (
x
,
z
)-plane; (f) is the graph projected on (
y
,
z
)-plane.Figure 12
Lorenz system for Experiment4 with a
=
10
,
b
=
8
/
3
,
c
=
28: (a) time series plot of x; (b) time series plot of z; (c) on the three-dimensional space; (d) projected on the (
x
,
y
)-plane; (e) projected on the (
x
,
z
)-plane; (f) projected on the (
y
,
z
)-plane.
(a)
(b)
(c)
(d)
(e)
(f)Experiment 5.
We consider the new chaotic system [8](14)
d
x
d
t
=
a
x
-
y
,
d
y
d
t
=
-
4
a
y
+
x
z
+
m
x
3
,
d
z
d
t
=
-
a
c
z
+
x
3
y
+
b
z
2
,where x
,
y, and z are state variables and a
,
b
,
c, and m are real parameters, which satisfy the following initial conditions:(15)
x
0
=
0.5
,
y
0
=
0
,
z
0
=
0We choose Chebyshev nodes and the number of nodesM
=
40. Figure 13 is obtained by using the current method with the parameters a
=
0.5
,
b
=
-
0.1, c
=
1.5, and m
=
0.12. In Figure 13, (a) is the time series plot; (b) is the three-dimensional space graph; (c) is the graph projected on (
x
,
y
)-plane; (d) is the graph projected on (
x
,
z
)-plane; (e) is the graph projected on (
y
,
z
)-plane.Figure 13
The new chaotic system for Experiment5 with a
=
0.5
,
b
=
-
0.1
,
c
=
1.5
,
m
=
0.12: (a) time series plot; (b) on the three-dimensional space; (c) projected on the (
x
,
y
)-plane; (d) projected on the (
x
,
z
)-plane; (e) projected on the (
y
,
z
)-plane.
(a)
(b)
(c)
(d)
(e)


## 4. Conclusions and Remarks

In this paper, the Lorenz System has solved by using BLICM. These numerical experiments illustrate that the numerical results of the present method are the same as the experimental results.All computations are performed by the MatlabR2017b software packages.



---

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