# Qualitative Analysis of a Hyperchaotic Lorenz-Stenflo Mathematical Model via the Caputo Fractional Operator

Authors: Chernet Tuge Deressa; Sina Etemad; Mohammed K. A. Kaabar; Shahram Rezapour
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

The Lorenz-Stenflo mathematical model describes a complex dynamical behavior related to atmospheric acoustic-gravity waves. In this study, qualitative analysis of the four-dimensional hyperchaotic Lorenz-Stenflo system via the Caputo fractional derivative is implemented. By using the Matignon stability criterion, the local stability analysis of the system showed that all the equilibrium points of the system are locally unstable. Calculation of the Lyapunov exponents along with the relevant bifurcation diagrams with respect to different fractional orders exposed the hyperchaotic dynamical behavior for the system. Bifurcation diagrams for all the four parameters in the system also showed the hyperchaotic nature of the Lorenz-Stenflo system. Different phase attractors of the system corresponding to different fractional derivatives and parameters are presented to specify the dynamical nature of the system. The Lorenz-Stenflo system showed sensitivity to initial conditions. The master and slave systems showed a strong correlation among themselves, as verified by graphs of time series solutions of the two systems.


---

## Body

## 1. Introduction

In recent years, fractional operators have been applied to develop mathematical models for which we can investigate different dynamical systems in some areas such as mathematical biology, epidemiology, and engineering [1].Different researchers introduced several concepts of fractional derivatives. Some of them are the Caputo fractional-order derivative [2–4], the Riemann Liouville fractional-order derivative [4], the Caputo-Fabrizio fractional-order derivative [5], and the Atangana-Baleanu fractional derivatives [6].In the mathematical models with integer derivatives, the derivative orders give the instantaneous rate of change of the function. On the other hand, in the case of mathematical models with fractional derivatives, the parameters of dynamical systems represent the memory index of variation of the function [7].Thus, the advantage of using the memory index of a dynamical system with fractional derivative made fractional derivatives advantageously applied in the fields of chaotic dynamics, epidemiological modeling, and several other fields. In particular, one can see these applications in the modeling of COVID-19 based on real information from Pakistan [8], designing the SEIR model of COVID-19 [9], the modeling of the memristor-based hyperchaotic circuit via nonsingular operator [10, 11], the thermostat model via Bernstein polynomials [12], the modeling of coronavirus by the Caputo operator [13], the optimal control of nonsingular tumor-immune surveillance [14], the SEIRA model [15], designing a bank data with fractal-fractional operators [16], the nonlinear modeling of the Navier system [17], analyzing an Atangana–Baleanu SEAIR model [18], compartmental disease modeling [19], and the references therein.Lorenz developed the famous three-dimensional mathematical model to study the dynamics and behaviors of the atmosphere in 1963 [20–22]. Later on, in 1996, Lennart Stenflo modified the Lorenz system, adding a fourth parameter to consider the evolution of finite-amplitude acoustic-gravity waves in a rotating atmosphere, and developed a new four-dimensional Lorenz-Stenflo mathematical model. Since then, many studies have been conducted on the complex dynamical behaviors of the system by applying both integer and fractional derivatives. For further information, see [21] and the references therein.Some of the studies conducted on the complex dynamics of Lorenz-Stenflo systems are reviewed as follows: Zhang et al. [22] investigated the qualitative properties of the higher integer-order Lorenz-Stenflo Chaotic system which appeared in mathematical physics. The authors proved that the higher-order Lorenz-Stenflo system is globally stable. A globally attractive set of Lorenz-Stenflo systems indicating the evolution of a finite amplitude of acoustic gravity waves contained in the rotating atmosphere is investigated by Zhang et al. [20]. Adaptive control synchronization and circuit implementation of the ordinary Lorenz-Stenflo system with an integer order was established by Yang and Wu [23]. An analysis on the dynamics of the Lorenz-Stenflo system via fractional operators along with sets of parameters is done through the spectra of the Lyapunov exponent type, bifurcation graphs, and 0-1 test.The complex dynamics of Lorenz-Stenflo dynamical systems are analyzed with Lyapunov exponents, bifurcation diagrams, and the 0-1 test by using the Adomian decomposition numerical scheme and fractional-order representation of the model in the sense of the Riemann-Liouville integral [24]. A robust chaos suppression control is designed and applied via sliding mode control to the integer-order Lorenz-Stenflo system constrained to uncertainties and nonlinearities [25].There is only limited literature devoted to studying the complex dynamics of the Lorenz-Stenflo systems by using integer-order derivatives and fractional-order derivatives. Moreover, to the best of the authors’ knowledge, no study is conducted on the qualitative analysis of the Lorenz-Stenflo mathematical model in the sense of the Caputo fractional derivative and by using a numerical scheme developed by Garrappa [26], which is the main focus of this study.Among the several notions of fractional derivatives presented above, the Caputo fractional derivative gives the opportunity of including the classical initial conditions in a mathematical model and the Caputo fractional derivative of a constant is zero, which is not the case for instance in the Riemann-Liouville fractional derivative. Moreover, it has a MATLAB code that can obtain phase portraits and time series solutions of chaotic or hyperchaotic systems using a numerical scheme named the predictor-corrector method reported by Garrappa and applied in several studies on chaotic systems [27–29].The manuscript is organized as follows: Section2 involves the statement of the problem. In this section, the fractional representation of the Lorenz-Stenflo system (LS) is designed via the Caputo derivatives after recalling some critical preliminary definitions of fractional operators. Section 3 is devoted to analyzing the local stability of the mentioned LS system. In Section 4, the numerical scheme for the fractional-order LS system is developed based on a research reported by Garrappa [26]. The broader part of the manuscript is devoted to Lyapunov exponents, bifurcation, and chaos of the LS system in Section 5, followed by sensitivity analysis to initial conditions in Section 6. In Section 7, the development of the master and slave system to create a strong relationship between the systems using a coupling function is considered. Finally, a conclusion and a reference list are provided.


## 2. Statement of the Problem

The Lorenz-Stenflo (LS) system of differential equations is given by [30, 31].
(1)ẋ=ay−x+cw,ẏ=rx−y−xz,ż=−kz+xy,ẇ=−x−aw.The variablex is the intensity of motion of the fluid; y and z denote the horizontal and vertical direction temperature variations of the atmosphere. Parameters a, r, and k are all positive real numbers and are the Prandtl number, c is the rotation number, r is the Rayleigh number, and k is a geometric parameter. Parameters a and k depend on the material and geometrical properties of the layer of the fluid. Parameter r is proportional to the difference in temperature.We used fractional operators to find the hidden properties of the nature of solution of the LS system that are not observable via integer-order derivatives.


### 2.1. The Formulation of the LS System via the Caputo Fractional Derivative

Here, we shortly give several definitions of fractional operators pertinent to our study.Definition 1 (see [2, 3]).
The Riemann-Liouville fractional integral for a continuous functionf:0,+∞⟶ℝ is defined by
(2)I0RLtηft=1Γη∫0tt−τη−1fτdτ,η∈0,1,t>0.Definition 2 (see [2, 3]).
The Riemann-Liouville fractional derivative for a continuous functionf:0,+∞⟶ℝ is defined by
(3)D0RLtηft=1Γ1−ηddt∫0tt−τ−ηfτdτ,η∈0,1,t>0.Definition 3 (see [2, 3]).
The Caputo fractional derivative for a continuous functionf:0,+∞⟶ℝ is defined by
(4)D0Ctηf=1Γ1−η∫0tt−τ−ηddτfτdτ,η∈0,1,t>0.This section develops the fractional order representation of the LS system (1). The Caputo fractional derivative is used because it can incorporate customary initial conditions in the model, unlike the Riemann-Liouville fractional derivative. Furthermore, the Caputo fractional derivative of a constant is zero, which is not the case in the Riemann-Liouville fractional derivative. Accordingly, the Caputo fractional representation of the LS system (1) is given by
(5)D0Ctηx=L1x,y,z,w,D0Ctηy=L2x,y,z,w,D0Ctηz=L3x,y,z,w,D0Ctηw=L4x,y,z,w,where
(6)L1t,x,y,z,w=ay−x+cw,L2t,x,y,z,w=rx−y−xz,L3t,x,y,z,w=−kz+xy,L4t,x,y,z,w=−x−aw,with the initial value conditions x0=0.5,y0=0.5,z0=5, and w0=0.5.The fractional-order derivativeη∈0,1 and the parameter values are a=10,r=30,k=0.7, and c=4.The specific objectives of this study are to analyze the local stability of system (5), to investigate the chaotic behavior of system (5) via Lyapunov exponents, and to plot the bifurcation diagrams, attractors, and time series trajectories of the system by variation of fractional orders and four parameters of system (5). In addition, sensitivity to initial conditions and synchronization of the LS fractional system are considered.


## 2.1. The Formulation of the LS System via the Caputo Fractional Derivative

Here, we shortly give several definitions of fractional operators pertinent to our study.Definition 1 (see [2, 3]).
The Riemann-Liouville fractional integral for a continuous functionf:0,+∞⟶ℝ is defined by
(2)I0RLtηft=1Γη∫0tt−τη−1fτdτ,η∈0,1,t>0.Definition 2 (see [2, 3]).
The Riemann-Liouville fractional derivative for a continuous functionf:0,+∞⟶ℝ is defined by
(3)D0RLtηft=1Γ1−ηddt∫0tt−τ−ηfτdτ,η∈0,1,t>0.Definition 3 (see [2, 3]).
The Caputo fractional derivative for a continuous functionf:0,+∞⟶ℝ is defined by
(4)D0Ctηf=1Γ1−η∫0tt−τ−ηddτfτdτ,η∈0,1,t>0.This section develops the fractional order representation of the LS system (1). The Caputo fractional derivative is used because it can incorporate customary initial conditions in the model, unlike the Riemann-Liouville fractional derivative. Furthermore, the Caputo fractional derivative of a constant is zero, which is not the case in the Riemann-Liouville fractional derivative. Accordingly, the Caputo fractional representation of the LS system (1) is given by
(5)D0Ctηx=L1x,y,z,w,D0Ctηy=L2x,y,z,w,D0Ctηz=L3x,y,z,w,D0Ctηw=L4x,y,z,w,where
(6)L1t,x,y,z,w=ay−x+cw,L2t,x,y,z,w=rx−y−xz,L3t,x,y,z,w=−kz+xy,L4t,x,y,z,w=−x−aw,with the initial value conditions x0=0.5,y0=0.5,z0=5, and w0=0.5.The fractional-order derivativeη∈0,1 and the parameter values are a=10,r=30,k=0.7, and c=4.The specific objectives of this study are to analyze the local stability of system (5), to investigate the chaotic behavior of system (5) via Lyapunov exponents, and to plot the bifurcation diagrams, attractors, and time series trajectories of the system by variation of fractional orders and four parameters of system (5). In addition, sensitivity to initial conditions and synchronization of the LS fractional system are considered.


## 3. Stability Analysis of the LS System (5)

In this section of the manuscript, analysis of the local stability of the mentioned system (5) is conducted. There are several methods for performing local stability analysis of systems. Some of these methods are the Laplace transform techniques and the Matignon criterion. For the local stability analysis of the LS system (5), we used the Matignon method, because it is most commonly used in literature [30, 31]. The Matignon condition for fractional stability is given by
(7)∣argλJ∣>ηπ2,where J is the Jacobian matrix of the system, λJ denotes the class of all eigenvalues of J, and η∈0,1 is the order of the LS system (5). The LS system (5) is said to be locally asymptotically stable whenever all the eigenvalues of J satisfy the Matignon criterion. Firstly, we considered the equilibrium points of (5), given by E0=0,0,0,0, and constant values of the parameters given by k=0.7, r=30, c=4, and a=10.The matrixJ at E0=0,0,0,0 is given by
(8)J0=−aa0cr−10000−k0−100−a.The eigenvalues of the matrixJ0 are λ1=12.3287,λ2=−23.2066,λ3=−10.1221, and λ4=−0.7000, and the corresponding arguments of the eigenvalues are 0 for λ1 and π for the remaining eigenvalues. It is easy to conclude that E0 is locally unstable, because the absolute value of the arguments does not satisfy the Matignon criteria (7). The presence of an eigenvalue with a positive real part, λ1=12.3287, is necessary for the SL system to exhibit a double-scroll attractor [31].The general equilibrium points of the system in terms of the parameters are given by(9)∓aθ,∓a2+caθ,a2r−a2+ca2,±θ,θ=ka2r−a2+ca2a2+c.Substituting the parameter valuesk=0.7, r=30, c=4, and a=10, the equilibrium points are given by E1=−4.4150,−4.5916,28.9600,+0.4415 and E2=+4.4150,+4.5916,28.9600,−0.4415.The Jacobian matrix corresponding to the equilibrium pointE1=−4.4150,−4.5916,28.9600,+0.4415 is given by
(10)J1=−1010041.04−14.4150−4.5916−4.415−0.70−100−10.The eigenvalues ofJ1 are λ1,2=0.2026±5.7129i and λ3,4=−11.0526±1.3824i. The arguments of the eigenvalues are ±1.5353forλ1,2 and ±3.0172forλ3,4. It can then be inferred that the Matignon criterion of local stability is satisfied for λ3,4, since ∣±3.0172∣>ηπ/2, ∀η∈0,1. On the other hand, the Matignon criterion is not satisfied for λ1,2 since ∣±1.5353∣>ηπ/2 only for η<0.9774. That is, the equilibrium point E1is unstable for η≥0.9774.The Jacobian matrix corresponding to the equilibrium pointE2=+4.4150,+4.5916,28.9600,−0.4415 is given by
(11)J2=−1010041.04−1−4.41504.59164.415−0.70−100−10.The eigenvalues and their corresponding arguments of Jacobian matrixJ2 are identical to those of Jacobian matrix J1. Hence, the conclusion in relation to the local stability of equilibrium point E2 is the same as that in equilibrium point E1.In summary, the equilibria of the fractional-order system (5) are unstable for the fractional order of η∈0.9774,1,k=0.7,r=30,c=4, and a=10.


## 4. Numerical Solution

This section presents the numerical method used to obtain the phase portraits and time series solutions of the fractional-order LS system (5). To solve fractional-order systems, there are several numerical and analytical methods such as the predictor-corrector method, Homotopy method, and Adomian decomposition method. In this study, we use the predicator-corrector method, because it has a MATLAB code that can be used to obtain phase portraits and time series solutions of chaotic or hyperchaotic systems.The numerical method in [26] is applied to our system as follows: by starting from the first equation of (5) and applying the Reimann-Liouville fractional integral given by (2), we obtain the following system
(12)xt=x0+I0RLtηL1t,x,y,z,w,yt=y0+I0RLtηL2t,x,y,z,w,zt=z0+I0RLtηL3t,x,y,z,w,wt=w0+I0RLtηL4t,x,y,z,w.Applying the predictor-corrector method to (12), the integral terms are approximated and it becomes
(13)xt=x0+hηκmηL10+∑j=1m−1ψm−jηL1tj,xj,yj,zj,wj+ψ0ηL1tj,xmp,ymp,zmp,wmp,yt=y0+hηκmηL20+∑j=1m−1ψm−jηL2tj,xj,yj,zj,wj+ψ0ηL2tj,xmp,ymp,zmp,wmp,zt=z0+hηκmηL30+∑j=1m−1ψn−jηL3tj,xj,yj,zj,wj+ψ0ηL3tj,xmp,ymp,zmp,wmp,wt=w0+hηκmηL40+∑j=1m−1ψm−jηL4tj,xj,yj,zj,wj+ψ0ηL4tj,xmp,ymp,zmp,wmp,whereL1,L2,L3, and L4 are defined by (5)
(14)Li0=Lix0,y0,z0,w0,i=1,2,3,and
(15)κmη=m−1η−mηm−η−1Γ2+η,ψmη=m−1η+1−2mη+1+mηm+1Γ2+η,m=1,2,3,⋯,ψ0η=1Γ2+η.Moreover, the predictors are given as follows:(16)xptm=x0+hη∑j=1m−1ψm−j−1ηL1tj,xj,yj,zj,wj,yptm=y0+hη∑j=1m−1ψm−j−1ηL2tj,xj,yj,zj,wj,zptm=z0+hη∑j=1m−1ψm−j−1ηL3tj,xj,yj,zj,wj,wptm=w0+hη∑j=1m−1ψm−j−1ηL4tj,xj,yj,zj,wj.By Garrappa [26], the discretization method applied in this study is stable and convergent. Thus, the advantages of this discretization method are that it is both stable and convergent, and of course, it has a MATLAB code that can be used for plotting the attractors of the LS system (5). Several other advantages of the method are described in [27].As mentioned above, in this research, the predictor-corrector technique is utilized to obtain the attractors and time series trajectories of the LS system (5). Moreover, the Lyapunov exponents are obtained using an algorithm by Danca et al. [32]. The code is developed to determine all Lyapunov exponents of a class fractional-order system modeled by the Caputo derivative. The predictor-corrector Adams-Bashforth-Moulton numerical method is the underlying numerical method used in this code.


## 5. Lyapunov Exponents, Bifurcation Diagrams, and Hyperchaotic Behavior of the LS System

This section is devoted to investigating the chaotic or hyperchaotic nature of the LS system (5). The magnitude of the chaos is quantified via the Lyapunov exponents (LE). Furthermore, bifurcation diagrams caused by the variation of the parameters of (5) are portrayed using the values of the parameters and the initial value condition presented in (5).


### 5.1. Lyapunov Exponents for the Fractional Orderη

Based on the Danca algorithm mentioned above, some of the Lyapunov exponents corresponding to different fractional orders of the LS system (5) are shown in Table 1. The simulation is made to run for 300 s.Table 1
LEs for different fractional ordersη.
η0.8560.8570.8660.8680.8880.90.9560.980.9880.9991LE114.177314.121413.628213.521112.496911.92179.58118.73028.46448.11278.0815LE21.08241.07881.04611.0390.97040.93130.76820.7070.68770.6620.6597LE3−5.5673−5.553−5.4255−5.3973−5.1187−4.9549−4.2295−3.9409−3.8478−3.7225−3.7113LE4−9.6713−9.6712−9.6592−9.654−9.5574−9.4648−8.7868−8.4139−8.2829−8.0987−8.0818Sum0.0211−0.024−0.4104−0.4912−1.2088−1.5667−2.667−2.9176−2.9786−3.0465−3.0519Based on the LE values in Table1, it is followed that the dynamical system (5) possibly exhibits hyperchaotic behavior for η∈0.857,1, because in each column of Table 1, there are two positive LEs. The positive Lyapunov exponent ensures the sensitive dependence on the choice of initial conditions (local instability of the system in the state space), which is only one property of a chaotic system. Crudely, for a given dynamical system to be chaotic, it must have the following properties: sensitivity to initial values, topological transitivity, and also dense periodic orbits.System (5) is dissipative for η∈0.857,1 as the sum of the LEs in every column of Table 1 is negative. The dissipativity of the system ensures the existence of an attractor. From the experimental result in Table 1, it is observed that for values of η∈0,0.857, the system is not dissipative. For instance, the sum of LEs for the fractional order η=0.856 is 0.0211. Therefore, to get a dissipative hyperchaotic system in this study, we considered the fractional order η∈0.857,1.The Kaplan-Yorke dimension denoted bydimLE of system (5) for η∈0.857,1 can be calculated. For instance, for the fractional order η=0.856, the Kaplan-Yorke dimension is
(17)dimLE=3+14.1214+1.0788−5.5530−9.6712=3.9975.The Kaplan-Yorke dimension corresponding toη=0.888 is
(18)dimLE=3+12.4969+0.9704−5.1187−9.5574=3.8735.The Kaplan-Yorke dimension corresponding toη=0.988 is
(19)dimLE=3+8.4644+0.6620−3.7225−8.2829=3.6524.The Kaplan-Yorke dimension corresponding toη=1 is
(20)dimLE=3+8.0815+0.6597−3.7113−8.0818=3.6224.From the abovedimLEs, it can be said that the Kaplan-Yorke dimension is decreased from 3.9975 to 3.6224, as the fractional order is increased from 0.857 to 1. Note that the Kaplan-York dimension determines the upper bound of the Hausdorff dimensions. In this study, the Hausdorff dimensions of the attractors corresponding to different fractional orders are nonintegers and system (5) has strange attractors with fractal structures. Moreover, in Table 1 and the Kaplan-Yorke dimensions calculated above, one can find that the significance of the hyperchaotic attractors decreases as the fractional order increases from 0.857 to 1. The loss of significance is described by decreasing the sum of the positive LEs in Table 1 and also by decreasing the dimensions as the fractional order increases from 0.857 to 1. The loss of significance is also observable from the bifurcation diagram of the fractional orders shown in Figure 1 in the next section.Figure 1
Bifurcation diagram of (5) due to orders from 0.93 to 1.


### 5.2. Bifurcation for the Fractional Orderη

The bifurcation diagram of the fractional orderη is obtained by varying its value in the interval 0.857,1.00 with a time step of h=0.001. The bifurcation diagram due to the variations of the order is illustrated in Figure 1. According to this figure, for values of η∈0.933,0.953, the given system exhibits oscillation with stability and the system excites a chaotic behavior for fractional values η≥0.953.The simulation results corresponding to different fractional orders are portrayed in Figures2–4 to verify the conclusion made in Figure 1.Figure 2
Phase portraits of (5) for η=0.950 projected on three planes: (a) x-y, (b) y-w, and (c) y-z.
(a)(b)(c)Figure 3
Time series trajectories of the LS system (5) for the fractional order η=0.950.Figure 4
Attractors of the hyperchaotic LS system (5) projected on the planes: (a) x-y, (b) x-z, (c) x-w, (d) y-z, and (e) y-w.
(a)(b)(c)(d)(e)It can be seen in Figures2(a)–2(c) that the system exhibits oscillation with stability for the fractional order η=0.95∈0.933,0.953, which is in agreement with the conclusion in Figure 1. The time series solution of system (5) is shown for η=0.95 in Figure 3, where the solution trajectories have the property of oscillation and stability as claimed from the bifurcation diagram shown in Figure 1. Moreover, Figures 4 and 5 characterize that system (5) exhibits a hyperchaotic behavior at η=0.974 confirming what is predicted by the bifurcation diagram in Figure 1.Figure 5
Time series trajectories of system (5) for the fractional order η=0.974.


### 5.3. Bifurcation for Parameterr

To get the bifurcation diagram of parameterr of the hyperchaotic system (5), the fractional order used is η=0.97 and k=0.7, c=4, a=10, and r∈29.5,30.5 with a time step of h=0.001 and the initial value condition is y0=0.5,0.5,5,0.5. The simulation is made to run for 100 s. The bifurcation graph due to the variation of r is illustrated in Figure 6.Figure 6
Bifurcation diagram forr∈29.5,30.5 and η=0.974.It can be inferred in Figure6 that system (5) exhibits a hyperchaotic behavior for the parameters considered in the simulation for r∈29.5,30.5. The phase portraits of system (5) shown in Figures 7 and 8 for r=29.6 and r=30.4 verify the conclusion from the bifurcation diagram of Figure 6.Figure 7
Attractors of the LS system (5) for r=29.6 projected on the planes: (a) z-y, (b) y-w, and (c) z-w.
(a)(b)(c)Figure 8
Attractors of the LS system (5) for r=30.4 projected on the planes: (a) x-y, (b) x-z, and (c) x-w.
(a)(b)(c)


### 5.4. Bifurcation for Parameterc

For the bifurcation diagram of parameterc of the hyperchaotic system (5), the fractional order used is η=0.974, k=0.7, and a=10 and the bifurcation parameter is made to vary in the interval [3.5, 4.5] with a time step of h=0.001, with initial value condition y0=0.5;0.5;5;0.5, and the simulation is made to run for 100 s. The bifurcation diagram due to the variation of parameter c is illustrated in Figure 9.Figure 9
Bifurcation diagram forc∈3.5,4.5 and η=0.974.It is observable in Figure5 that system (5) shows hyperchaotic dynamics for parameter c in [3.5, 4.5]. Moreover, the phase portraits of system (5) for c=3.6and4.4 shown in Figure 10 verify the conclusion drawn in Figure 9.Figure 10
Attractors of the LS system (5) for c=3.6 and 4.4 projected on the planes: (a) y-w, (b) x-z, and (c) x-w.
(a)(b)


### 5.5. Bifurcation for Parameterk

To obtain the bifurcation diagram due to parameterk of the hyperchaotic system (5), the fractional order used is η=0.974. The values of the remaining parameters are k=0.7 and a=10, and the bifurcation parameter k is made to vary in the interval [0.5, 1] with a time step of h=0.001. The initial value condition is y0=0.5;0.5;5;0.5 and the simulation is made to run for 100 s. The bifurcation graph due to the variation of k is illustrated in Figure 11.Figure 11
Bifurcation diagram of the LS system (5) for parameter k in the interval [0.5, 1].It is clear in Figure11 that system (5) shows hyperchaotic dynamics for parameter k of the interval [0.5, 1], the given values of the parameters, and fractional orders. Moreover, the attractors of system (5) for values of parameter k in [0.5, 1] are approximately similar to the figures shown in Figure 10.


### 5.6. Bifurcation for Parametera

To obtain the bifurcation diagram due to parametera of the hyperchaotic system (5), the fractional derivative, the initial condition, the time, and the values of the remaining parameters are the same as those used in Section 5.5. The bifurcation diagram due to the variation of parameter a is illustrated in Figure 12.Figure 12
Bifurcation diagram of the LS system (5), for the values of a in [0.5,1].It is clear in Figure12 that system (5) shows the hyperchaotic dynamics for parameter a in [9.5,10.5], the given parameters, and the fractional orders. Moreover, the attractors of system (5) for different values of a are shown in Figure 13, verifying the hyperchaotic nature of system (5).Figure 13
Attractors of the LS system (5) for different values of a in [0.5, 1] projected on: (a) x-y-z space, (b) x-w-z space, (c) z-w plane, and (d) y-w plane.
(a)(b)(c)(d)


## 5.1. Lyapunov Exponents for the Fractional Orderη

Based on the Danca algorithm mentioned above, some of the Lyapunov exponents corresponding to different fractional orders of the LS system (5) are shown in Table 1. The simulation is made to run for 300 s.Table 1
LEs for different fractional ordersη.
η0.8560.8570.8660.8680.8880.90.9560.980.9880.9991LE114.177314.121413.628213.521112.496911.92179.58118.73028.46448.11278.0815LE21.08241.07881.04611.0390.97040.93130.76820.7070.68770.6620.6597LE3−5.5673−5.553−5.4255−5.3973−5.1187−4.9549−4.2295−3.9409−3.8478−3.7225−3.7113LE4−9.6713−9.6712−9.6592−9.654−9.5574−9.4648−8.7868−8.4139−8.2829−8.0987−8.0818Sum0.0211−0.024−0.4104−0.4912−1.2088−1.5667−2.667−2.9176−2.9786−3.0465−3.0519Based on the LE values in Table1, it is followed that the dynamical system (5) possibly exhibits hyperchaotic behavior for η∈0.857,1, because in each column of Table 1, there are two positive LEs. The positive Lyapunov exponent ensures the sensitive dependence on the choice of initial conditions (local instability of the system in the state space), which is only one property of a chaotic system. Crudely, for a given dynamical system to be chaotic, it must have the following properties: sensitivity to initial values, topological transitivity, and also dense periodic orbits.System (5) is dissipative for η∈0.857,1 as the sum of the LEs in every column of Table 1 is negative. The dissipativity of the system ensures the existence of an attractor. From the experimental result in Table 1, it is observed that for values of η∈0,0.857, the system is not dissipative. For instance, the sum of LEs for the fractional order η=0.856 is 0.0211. Therefore, to get a dissipative hyperchaotic system in this study, we considered the fractional order η∈0.857,1.The Kaplan-Yorke dimension denoted bydimLE of system (5) for η∈0.857,1 can be calculated. For instance, for the fractional order η=0.856, the Kaplan-Yorke dimension is
(17)dimLE=3+14.1214+1.0788−5.5530−9.6712=3.9975.The Kaplan-Yorke dimension corresponding toη=0.888 is
(18)dimLE=3+12.4969+0.9704−5.1187−9.5574=3.8735.The Kaplan-Yorke dimension corresponding toη=0.988 is
(19)dimLE=3+8.4644+0.6620−3.7225−8.2829=3.6524.The Kaplan-Yorke dimension corresponding toη=1 is
(20)dimLE=3+8.0815+0.6597−3.7113−8.0818=3.6224.From the abovedimLEs, it can be said that the Kaplan-Yorke dimension is decreased from 3.9975 to 3.6224, as the fractional order is increased from 0.857 to 1. Note that the Kaplan-York dimension determines the upper bound of the Hausdorff dimensions. In this study, the Hausdorff dimensions of the attractors corresponding to different fractional orders are nonintegers and system (5) has strange attractors with fractal structures. Moreover, in Table 1 and the Kaplan-Yorke dimensions calculated above, one can find that the significance of the hyperchaotic attractors decreases as the fractional order increases from 0.857 to 1. The loss of significance is described by decreasing the sum of the positive LEs in Table 1 and also by decreasing the dimensions as the fractional order increases from 0.857 to 1. The loss of significance is also observable from the bifurcation diagram of the fractional orders shown in Figure 1 in the next section.Figure 1
Bifurcation diagram of (5) due to orders from 0.93 to 1.


## 5.2. Bifurcation for the Fractional Orderη

The bifurcation diagram of the fractional orderη is obtained by varying its value in the interval 0.857,1.00 with a time step of h=0.001. The bifurcation diagram due to the variations of the order is illustrated in Figure 1. According to this figure, for values of η∈0.933,0.953, the given system exhibits oscillation with stability and the system excites a chaotic behavior for fractional values η≥0.953.The simulation results corresponding to different fractional orders are portrayed in Figures2–4 to verify the conclusion made in Figure 1.Figure 2
Phase portraits of (5) for η=0.950 projected on three planes: (a) x-y, (b) y-w, and (c) y-z.
(a)(b)(c)Figure 3
Time series trajectories of the LS system (5) for the fractional order η=0.950.Figure 4
Attractors of the hyperchaotic LS system (5) projected on the planes: (a) x-y, (b) x-z, (c) x-w, (d) y-z, and (e) y-w.
(a)(b)(c)(d)(e)It can be seen in Figures2(a)–2(c) that the system exhibits oscillation with stability for the fractional order η=0.95∈0.933,0.953, which is in agreement with the conclusion in Figure 1. The time series solution of system (5) is shown for η=0.95 in Figure 3, where the solution trajectories have the property of oscillation and stability as claimed from the bifurcation diagram shown in Figure 1. Moreover, Figures 4 and 5 characterize that system (5) exhibits a hyperchaotic behavior at η=0.974 confirming what is predicted by the bifurcation diagram in Figure 1.Figure 5
Time series trajectories of system (5) for the fractional order η=0.974.


## 5.3. Bifurcation for Parameterr

To get the bifurcation diagram of parameterr of the hyperchaotic system (5), the fractional order used is η=0.97 and k=0.7, c=4, a=10, and r∈29.5,30.5 with a time step of h=0.001 and the initial value condition is y0=0.5,0.5,5,0.5. The simulation is made to run for 100 s. The bifurcation graph due to the variation of r is illustrated in Figure 6.Figure 6
Bifurcation diagram forr∈29.5,30.5 and η=0.974.It can be inferred in Figure6 that system (5) exhibits a hyperchaotic behavior for the parameters considered in the simulation for r∈29.5,30.5. The phase portraits of system (5) shown in Figures 7 and 8 for r=29.6 and r=30.4 verify the conclusion from the bifurcation diagram of Figure 6.Figure 7
Attractors of the LS system (5) for r=29.6 projected on the planes: (a) z-y, (b) y-w, and (c) z-w.
(a)(b)(c)Figure 8
Attractors of the LS system (5) for r=30.4 projected on the planes: (a) x-y, (b) x-z, and (c) x-w.
(a)(b)(c)


## 5.4. Bifurcation for Parameterc

For the bifurcation diagram of parameterc of the hyperchaotic system (5), the fractional order used is η=0.974, k=0.7, and a=10 and the bifurcation parameter is made to vary in the interval [3.5, 4.5] with a time step of h=0.001, with initial value condition y0=0.5;0.5;5;0.5, and the simulation is made to run for 100 s. The bifurcation diagram due to the variation of parameter c is illustrated in Figure 9.Figure 9
Bifurcation diagram forc∈3.5,4.5 and η=0.974.It is observable in Figure5 that system (5) shows hyperchaotic dynamics for parameter c in [3.5, 4.5]. Moreover, the phase portraits of system (5) for c=3.6and4.4 shown in Figure 10 verify the conclusion drawn in Figure 9.Figure 10
Attractors of the LS system (5) for c=3.6 and 4.4 projected on the planes: (a) y-w, (b) x-z, and (c) x-w.
(a)(b)


## 5.5. Bifurcation for Parameterk

To obtain the bifurcation diagram due to parameterk of the hyperchaotic system (5), the fractional order used is η=0.974. The values of the remaining parameters are k=0.7 and a=10, and the bifurcation parameter k is made to vary in the interval [0.5, 1] with a time step of h=0.001. The initial value condition is y0=0.5;0.5;5;0.5 and the simulation is made to run for 100 s. The bifurcation graph due to the variation of k is illustrated in Figure 11.Figure 11
Bifurcation diagram of the LS system (5) for parameter k in the interval [0.5, 1].It is clear in Figure11 that system (5) shows hyperchaotic dynamics for parameter k of the interval [0.5, 1], the given values of the parameters, and fractional orders. Moreover, the attractors of system (5) for values of parameter k in [0.5, 1] are approximately similar to the figures shown in Figure 10.


## 5.6. Bifurcation for Parametera

To obtain the bifurcation diagram due to parametera of the hyperchaotic system (5), the fractional derivative, the initial condition, the time, and the values of the remaining parameters are the same as those used in Section 5.5. The bifurcation diagram due to the variation of parameter a is illustrated in Figure 12.Figure 12
Bifurcation diagram of the LS system (5), for the values of a in [0.5,1].It is clear in Figure12 that system (5) shows the hyperchaotic dynamics for parameter a in [9.5,10.5], the given parameters, and the fractional orders. Moreover, the attractors of system (5) for different values of a are shown in Figure 13, verifying the hyperchaotic nature of system (5).Figure 13
Attractors of the LS system (5) for different values of a in [0.5, 1] projected on: (a) x-y-z space, (b) x-w-z space, (c) z-w plane, and (d) y-w plane.
(a)(b)(c)(d)


## 6. Sensitivity to Initial Conditions

This section examines the effect of different initial values on the dynamics of the hyperchaotic system (5). Since sensitivity to initial conditions is one of the important properties of chaotic or hyperchaotic systems, it is necessary to examine this property for the LS system (5) using different initial conditions. In this study, to examine sensitivity to initial conditions, we used the parameter values indicated in (5) and η=0.974. The initial condition considered is x0,0.5,5,0.5. The different values of x0 are shown in Figure 14.Figure 14
The time series solution of LS (5) corresponding to different initial conditions due to different values of x0.
(a)(b)It can be observed in Figure14 that all the trajectories coincided for about the first 5.1 seconds and they were divided into two for about 10.1 seconds. At about 15.1 seconds, they divided into three, diverging from each other. The values of x0of two trajectories that overlapped for the first 15.1 seconds are closer to each other than the value of x0of the third trajectory.


## 7. Synchronization of the Lorenz-Stenflo System

In this part of the study, we developed a master-slave system of the dynamic system (5). Firstly, two identical copies of system (5) are associated via a coupling function. Then, we performed a simulation of the coupled systems using different initial conditions. Finally, the error dynamics showed asymptotic stability and the phase portraits of the master and slave system showed a strong correlation.Let the master hyperchaotic model be given by(21)D0Ctqxt=ay−x+cw,D0Ctqyt=rx−y−xz,D0Ctqzt=−kz+xy,D0Ctqwt=−x−aw,and the slave system be given by
(22)D0Ctqxst=ays−xs+cws,D0Ctqyst=rxs−ys−xszs,D0Ctqzst=−kzs+xsys,D0Ctqwst=−xs−aws,+H−∂F∂XXs−X,where the coupling function is H−∂F/∂XXs−X,∂F/∂X is the Jacobian matrix of the system, H is a Hermitian matrix of appropriate size, X=x,y,z,w, and Xs=xs,ys,zs,ws. Moreover, the coupling function is given by
(23)H−∂F∂XXs−X=000−5.5z−32−2x0−y−x−2.60−0.5007xs−xys−yzs−zws−w,where a Hermitian matrixH is chosen to be
(24)H=−10100−1.5−2−30000−3.30−1.500−3.It can be shown that all the real parts of the eigenvalues ofH are negative. Thus, the Matignon criterion is satisfied; H is a Hermitian matrix.Now, considering the coupling function (23), the slave system becomes
(25)D0Ctqxst=ays−xs+cws−5.5ws−w,D0Ctqyst=rxs−ys−xszs+−32+zxs−x−2ys−y+xzs−z,D0Ctqzst=−kzs+xsys−yxs−x−xys−y−2.6zs−z,D0Ctqwst=−xs−aws−0.5xs−x+7ws−w.Defining the error as(26)e1=xs−x,e2=ys−y,e3=zs−z,e4=ws−w,the error dynamics is shown in
(27)D0Ctqe1D0Ctqe2D0Ctqe3D0Ctqe4=He=−10100−1.5−2−30000−3.30−1.500−3e1e2e3e4.To verify if the master and slave systems in equations (21) and (25) are correlated, the time series trajectories of the two systems, including the graphs of the dynamics of the error, are illustrated in Figures 15–19. The values of the parameters used are k=0.7,r=30,c=4,anda=10 together with the fractional order η=0.978. Moreover, the initial value conditions used for the master and the slave system are y0m=0.5,0.5,5,0.5 and y0s=5,5,5,5, respectively.Figure 15
Graphs of the error dynamics (27).Figure 16
Time series trajectories ofxt and xst.Figure 17
Time series orbits ofyt and yst.Figure 18
Time series orbits ofzt and zst.Figure 19
Time series orbits ofwt and wst.It can be observed in Figure15 that there is a robust correlation between the master and slave systems because the error graphs converge to zero approximately in the first 2.5 seconds, as shown in Figure 15.


## 8. Conclusion

This research study implemented a qualitative analysis on the four-dimensional Lorenz-Stenflo mathematical model in the sense of the Caputo operator. The Matignon local stability criterion showed that the equilibrium points of the LS system are locally unstable, implying that the LS system led to hyperchaotic dynamical behavior. The scheme developed by Garrappa approximated the numerical solution, and the corresponding MATLAB code was used to obtain all the figures in the study. The authors believe that the qualitative analysis made in this manuscript, the displayed figures, and the synchronization systems developed have revealed complex dynamical behaviors of the Lorenz-Stenflo system that are not obtained earlier by other studies on the system. It is also worth mentioning that this study have been done by including other concepts of fractional derivatives and comparing the results to the results obtained due to the Caputo fractional derivative, to be treated in the future by the authors or interested researchers in the area.



---

*Source: 4975104-2022-01-27.xml*