Numerical Simulation of the Fractional-Order Lorenz Chaotic Systems with Caputo Fractional Derivative

1  Introduction

In 1963, Edward Lorenz discovered a mathematical model for atmospheric convection. This model is also known as Lorenz chaotic system. The Lorenz system is widely used in electric circuits, forward osmosis and chemical reactions. In recent years, people have studied the chaotic behavior in the fractional dynamic system and found that the fractional dynamic system has unique properties that the integer dynamic system does not have. Therefore, the numerical simulation of fractional chaotic system is very important. In this paper, we simulate the fractional-order Lorenz chaotic dynamical systems [1–6] is as

{Dtα1x(t)=σ(y(t)−x(t)),Dtα2y(t)=ρx(t)−y(t)−x(t)z(t),t∈[0,T],Dtα3z(t)=x(t)y(t)−βz(t),(1)

with the initial conditions

x(0)=c1,y(0)=c2,z(0)=c3.(2)

The solutions of the fractional-order Lorenz chaotic dynamical systems are very hard to obtain analytically, and researchers, therefore, rely on numerical methods to provide an approximate solution that could be used for prediction. In the last decades, several numerical methods have been proposed. In [1], based on the qualitative theory, the authors investigated the existence and uniqueness of solutions for a class of fractional-order Lorenz chaotic systems (1). In [2], compared the dynamical regimes of fractional-order systems with dynamical regimes of the corresponding standard systems. In [3], Complex dynamics with interesting characteristics were presented by means of phase portraits, the largest Lyapunov exponent and bifurcation diagrams. In [4–6], the authors gave a dynamic analysis of a fractional-order Lorenz chaotic system. Although some numerical and analytical methods of the FDEs have been announced, such as spectral method [7–11], reproducing kernel method [12–19], homotopy perturbation method [20–23], high-precision numerical approach [24–27], and so on numerical and analytical methods [28–36]. These researchers all say their own approach can accurately simulate chaotic systems. In fact, since chaotic systems have no exact solution, researchers do not know which method is more accurate. For the numerical simulation of chaotic systems, it is necessary to use numerical methods to study the long time properties of solutions of the fractional order chaotic systems. This paper introduces a high-precision numerical method [24–27] for solving system (1). Some complex dynamic behavior of the fractional-order Lorenz chaotic systems are discovered by using the present numerical approach. We observe some novel dynamic behavior in numerical simulations which are unlike any that have been previously discovered in numerical simulations or theoretical studies. The simulation results of numerical experiments demonstrate the effectiveness of the present method.

Fractal and fractional calculus [37–48] have been widely concerned. In the last three decades, there have existed many inequivalent definitions [49–51] of fractional derivatives. The most famous of these definitions that have been widely popularized in the world of fractional calculus is Riemann-Liouville fractional definition, Gru¨nwald-Letnikov fractional derivative (GLFD) and Caputo fractional derivative definition.

Definition 1.1. Riemann-Liouville fractional derivative of order α of a function y(t) about t on the interval (t0,t) is defined as

t0RLDtαy(t)={1Γ(n−α)dndtn∫t0ty(τ)(t−τ)α−n+1dτ,0≤n−1<α<n.dny(t)dtn,α=n∈N.(3)

Definition 1.2. Caputo fractional derivative of order α of a function y(t)∈Cn[t0,t] about t on the interval (t0,t) is defined as

t0CDtαy(t)={1Γ(n−α)∫t0ty(n)(τ)(t−τ)α−n+1dτ,0≤n−1<α<n.dny(t)dtn,α=n∈N.(4)

Definition 1.3. The Gru¨nwald-Letnikov fractional derivative of α-order on the interval (t0,t) is defined as

t0GLDtαy(t)=limh→01hα∑j=0[(t−t0)/h]⁡(−1)j(jα)y(t−jh),(5)

where

(−1)j(jα)=(−1)jΓ(α+1)Γ(j+1)Γ(α−j+1).

Let

ωj(α)=(−1)j(jα)=(−1)jΓ(α+1)Γ(j+1)Γ(α−j+1),j=0,1,…,n.(6)

Therefore, the α-order Gru¨nwald-Letnikov derivative in Eq. (4) is transformed into the following form:

t0GLDtαy(t)=limh→01hα∑j=0[(t−t0)/h]⁡ωj(α)y(t−jh).(7)

Using Newton’s binomial theorem, we know

(1−z)n=∑j=0n⁡(−1)j(jn)zj=∑j=0n⁡ωj(n)zj.(8)

We can proof the following form:

∑j=0∞⁡ωj(α)zj=(1−z)α.(9)

(1−z)α is called 1-order generating function of α-order Gru¨nwald-Letnikov derivative on the interval (t0,t).

t0GLDtαy(t)=1hα∑i=0[(t−t0)/h]⁡ωjy(t−ih)+o(h).(10)

Theorem 1.4. If f(t) is continuous and differentiable of order n−1 on the interval (t0,t), and f(n)(t) is integrable, then Gru¨nwald-Letnikov derivative can be written in the following integral form:

t0GLDtαy(t)=∑j=0n−1⁡y(j)(t0)Γ(1+j−α)(t−t0)j−α+1Γ(n−α)∫0ty(n)(τ)(t−τ)α−n+1dτ,0≤n−1<α<n.(11)

and, the Gru¨nwald-Letnikov derivative, Riemann-Liouville derivative and Caputo derivative have the following relationship:

t0GLDtαy(t)=t0RLDtαy(t),t0CDtαy(t)=t0RLDtαy(t)−∑j=0n−1⁡y(j)(t0)Γ(1+j−α)(t−t0)j−α.(12)

Proof 1. For the proof, please refer to [49].

From Theorem (1.4), it follows that, if y(j)(t0)=0,j=0,1,…,n−1, then t0GLDtαy(t)=t0RLDtαy(t)=t0CDtαy(t). If exists j, such that y(j)(t0)≠0,j=0,1,…,n−1, we use Taylor form, let y¯(t)=y(t)−∑j=0n−1y(j)(t0)Γ(1+j−α)(t−t0)j−α, then t0GLDtαy¯(t)=t0RLDtαy¯(t)=t0CDtαy¯(t). For convenience of expression, we denote

t0GLDtαy(t)=t0RLDtαy(t)=t0CDtαy(t)=Dtαy(t).

2  Numerical Approach

In [24–27], in order to obtain higher precision, author given the construction method of generating function for arbitrary p, and then give a recursive method of fractional derivative and integral based on the generating function.

Definition 2.1. A p-order polynomial function is defined as

gp(z)=∑k=1p⁡1k(1−z)k.(13)

Theorem 2.2. The p order polynomial function gp(z) could be written as

gp(z)=∑k=0p⁡gkzk,(14)

where gk is the solutions of the following equation [24–27]:

(111⋯1123⋯p+112232⋯(p+1)2⋮⋮⋮⋯⋮12p3p⋯(p+1)p)(g0g1g2⋮gp)=−(012⋮p).(15)

Proof 2. It can be seen from Eqs. (13) and (14)

∑k=0p⁡gkzk=∑k=1p⁡1k(1−z)k.(16)

Substituting z=1 into Eq. (16), there is

∑k=0p⁡gk=0.(17)

Multiply both ends of Eq. (16) by z and find the first derivative of z, then

∑k=0p⁡(k+1)gkzk=∑k=1p⁡1k(1−z)k−z∑k=1p⁡(1−z)k−1.(18)

Then substitute z=1 into Eq. (18), then there is

∑k=0p⁡(k+1)gk=−1.

Multiply both ends of Eq. (16) by z, and then find the first derivative of z

∑k=0p⁡(k+1)2gkzk=∑k=1p⁡1k(1−z)k−3z∑k=1p⁡(1−z)k−1+z2∑k=2p⁡1k−1(1−z)k−2.(19)

Substituting z=1 into Eq. (19), we can derive

∑k=0p⁡(k+1)2gk=−2.

Repeating the above process, the following equation can be established:

{g0+g1+g2+⋯+gp=0,g0+2g1+3g2+⋯+(p+1)gp=−1,g0+22g1+32g2+⋯+(p+1)2gp=−2,⋮g0+2pg1+3pg2+⋯+(p+1)pgp=−p.

The matrix form of the equation is formula (15), which is proved by the theorem.

Definition 2.3. Ihe p-order generating function gpα(z) of α-order Gru¨nwald-Letnikov derivative is defined as

gpα(z)=(g0+g1z+⋯+gpzp)α.(20)

Theorem 2.4. If the generating function gpα(z) can be written as

gpα(z)=∑k=0∞⁡ωk(α,p)zk.(21)

where the subsequent coefficient ω(α,p)k can be recursively calculated by the following formula [24–27]:

{ωk(α,p)=g0α,k=0,ωk(α,p)=−1g0∑i=1p⁡gi(1−i1+αk)ωk−i(α,p),k=1,2,…,ωk(α,p)=0,k<0.(22)

Proof 3. Let z=0 in Eq. (20), then it can be proved ω0(α,p)=g0α. Rewritten Eq. (20), it can be seen that

(g0+g1z+⋯+gpzp)α=∑k=−∞∞⁡ωk(α,p)zk,(23)

where, if k<0, then ωk(α,p)=0. Find the first derivative of z on both sides of Eq. (23), it can be concluded that

α(g1+2g2z+⋯+pgpzp−1)(g0+g1z+⋯+gpzp)α−1=∑k=−∞∞⁡kωk(α,p)zk−1.(24)

Both sides of Eq. (24) are multiplied by (g0+g2z+⋯+gpzp) at the same time, it can be seen

α(g1+⋯+pgpzp−1)(g0+⋯+gpzp)α=(g0+⋯+gpzp)∑k=−∞∞⁡kωk(α,p)zk−1.(25)

Substituting Eq. (23) into Eq. (25), we can get

α(g1+⋯+pgpzp−1)∑k=−∞∞⁡ωk(α,p)zk=(g0+⋯+gpzp)∑k=−∞∞⁡kωk(α,p)zk−1.(26)

Using the zdωk(α,p)=ωk−d(α,p) property, the left end of Eq. (26) can be written as

∑k=−∞∞⁡α(g1ωk(α,p)+2g2ωk−1(α,p)+⋯+pgpωk−p+1(α,p))zk,(27)

then its right end can be written

∑k=−∞∞⁡[(k+1)g0ωk+1(α,p)+kg1ωk(α,p)+⋯+(k−p+1)gpωk−p+1(α,p)]zk.(28)

By comparing the same square coefficient of z in Eq. (27), we can get

α(g1ωk(α,p)+⋯+pgpωk−p+1(α,p))=(k+1)g0ωk+1(α,p)+⋯+(k−p+1)gpωk−p+1(α,p).(29)

Move Eq. (29) back one step, let k=k−1, then the equation can become

g0kωk(α,p)+g1(k−1−α)ωk−1(α,p)+⋯+gp(k−p−pα)ωk−p(α,p)=0.(30)

If k≠0, it is thus clear that

ωk(α,p)=−1g0[g1(1−1+αk)ωk−1(α,p)+g2(1−21+αk)ωk−2(α,p)+⋯+gp(1−p1+αk)ωk−p(α,p)](31)

where, when k<0, ωk(α,p)=0.The above formula is recursive, so the theorem is proved.

Corollary 2.5. The p-order generating function gpα(z) of α-order Gru¨nwald-Letnikov derivative could be written as [24–27]

gpα(z)=∑j=0∞⁡ωj(α,p)zj.(32)

where

{ωk(α,p)=g0,k=0,ωk(α,p)=−1g0∑i=0k−1⁡(1−i1+αk)giωk−i(α,p),k=1,2,….,p−1,ωk(α,p)=−1g0∑i=0p⁡(1−i1+αk)giωk−i(α,p),k=p,p+1,p+2,….,ωk(α,p)=0,k<0,(33)

The α-order Gru¨nwald-Letnikov fractional derivative with p-order generating function gpα(z) is given as

t0GLDtαy(t)=limh→01hα∑i=0[(t−t0)/h]⁡ωi(α,p)y(t−ih).(34)

Applying (34), an approximate computation scheme of the α-order Caputo fractional derivative with p-order generating function gpα(z) is given as

Dtαy(t)≃ym=1hα∑j=0m⁡ωj(α,p)y(t−jh)=1hαy(tk)−1hα∑i=1m⁡ωj(α,p)y(tk−j).(35)

and limm→∞ym=Dtαy(t).

So, a high-precision numerical approach of the fractional-order Lorenz chaotic systems (1) is given by

{x(tk)=σhα1(y(tk−1)−x(tk))+x(0)−∑j=1m⁡ωj(α1,p)(x(tk−j)−x(0)),y(tk)=hα2(ρx(tk)−y(tk−1)−x(tk)z(tk−1))+y(0)−∑j=1m⁡ωj(α2,p)(y(tk−j)−y(0)),z(tk)=hα3(x(tk)y(tk)−βz(tk−1))+z(0)−∑j=1m⁡ωj(α3,p)(z(tk−j)−z(0)).(36)

3  Numerical Experiment

In this section, some numerical examples are studied. Some novel chaotic behaviors are shown. We consider the systems (1) with the initial conditions x(0)=−15,y(0)=−15,z(0)=20. h=0.01,p=20. The complex dynamic behaviors of the systems (1) are shown in Figs. 1–21. Fig. 1 shows time series plots of the systems (1) with parameters [σ,ρ,β]=[10,28,8/3] at different α. Figs. 2 and 3 show time series plots of the systems (1) with different parameters [σ,ρ,β] and different fractional derivative α. Figs. 4–19 projected on the (x,y), (x,z), (y,z)-plane and show phase diagram of x, y, z of the systems (1) with different parameters [σ,ρ,β] and different fractional derivative α. Figs. 20 and 21 show chaotic attractor of fractional-order Lorenz systems (1) with different parameters [σ,ρ,β] and different fractional derivative α. The rich chaotic attractor of fractional-order Lorenz systems (1) is shown in Figs. 20 and 21. The simulation results of integer order are in good agreement with those of other methods. In this paper, many novel chaotic attractors for fractional systems are obtained.

Figure 1: Time series plots of the systems (1) at [σ,ρ,β]=[10,28,8/3] with parameters α

Figure 2: Time series plots of the systems (1) with different parameters

Figure 3: Numerical results for the systems (1) at α=[1.16,1.18,1.15],[σ,ρ,β]=[35,20,15], T=100

Figure 4: Numerical results for the systems (1) at α=[0.9,1.29,1.28],[σ,ρ,β]=[8,12,10], T=100

Figure 5: Numerical results for the systems (1) at α=[1.29,1.2,1.1],[σ,ρ,β]=[8,12,10], T=100

Figure 6: Numerical results for the systems (1) at α=[0.95,1.02,1.03],[σ,ρ,β]=[10,28,8/3], T=100

Figure 7: Numerical results for the systems (1) at α=[1.1,1.15,1.18],[σ,ρ,β]=[10,28,8/3], T=100

Figure 8: Numerical results for the systems (1) at α=[0.9,1.02,1.27],[σ,ρ,β]=[5,12,10], T=100

Figure 9: Numerical results for the systems (1) at α=[0.9,1.18,1.38],[σ,ρ,β]=[3,10,11], T=100

Figure 10: Numerical results for the systems (1) at α=[1,1.18,1.38],[σ,ρ,β]=[1,8,15], T=100

Figure 11: Numerical results for the systems (1) at α=[1.22,2.21,1.19],[σ,ρ,β]=[20,10,9]

Figure 12: Numerical results for the systems (1) at α=[1,1.12,1.18],[σ,ρ,β]=[4,9,9], T=100

Figure 13: Numerical results for the systems (1) at α=[1,1.12,1.18],[σ,ρ,β]=[6,8,9], T=100

Figure 14: Numerical results for the systems (1) at α=[1,0.98,1.15],[σ,ρ,β]=[15,20,10], T=100

Figure 15: Numerical results for the systems (1) at α=[1.12,1.26,1.28],[σ,ρ,β]=[30,30,10], T=100

Figure 16: Numerical results for (1) at α=[1.35,1.25,1.2],[c1,c2,c3]=[−18,−25,20],[σ,ρ,β]=[35,15,9], T=100

Figure 17: Numerical results for the systems (1) at α=[0.92,0.99,1.01],[σ,ρ,β]=[20,30,8/3], T=100

Figure 18: Numerical results for the systems (1) at α=[0.92,1,1],[σ,ρ,β]=[5,30,8/3], T=100

Figure 19: Numerical results for the systems (1) at α=[1,1,1],[σ,ρ,β]=[10,28,8/3], T=100

Figure 20: Chaotic attractor of fractional-order Lorenz systems (1) with different parameters

Figure 21: Chaotic attractor of fractional-order Lorenz systems (1) with different parameters

4  Conclusions and Remarks

In this paper, some complex dynamic behavior of fractional-order Lorenz chaotic systems are shown by using the present method. We observe many novel dynamic behaviors in numerical experiments which are unlike any that have been previously discovered in numerical experiments or theoretical studies. We investigate the influence of α1, α2, α3 on the numerical solution of fractional-order Lorenz chaotic systems. The simulation results of integer order are in good agreement with those of other methods. The results presented in this paper suggested that the present numerical method is also readily applicable to a more chaotic system.

All computations are performed by the MatlabR2017b software.

Acknowledgement: The authors would like to express their thanks to the unknown referees for their careful reading and helpful comments.

Funding Statement: This paper is supported by the Natural Science Foundation of Inner Mongolia [2021MS01009] and Jining Normal University [JSJY2021040, Jsbsjj1704, jsky202145].

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.