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

The Lorenz-Sten ﬂ o mathematical model describes a complex dynamical behavior related to atmospheric acoustic-gravity waves. In this study, qualitative analysis of the four-dimensional hyperchaotic Lorenz-Sten ﬂ o 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 di ﬀ erent 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-Sten ﬂ o system. Di ﬀ erent phase attractors of the system corresponding to di ﬀ erent fractional derivatives and parameters are presented to specify the dynamical nature of the system. The Lorenz-Sten ﬂ o system showed sensitivity to initial conditions. The master and slave systems showed a strong correlation among themselves, as veri ﬁ ed by graphs of time series solutions of the two systems.


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].
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 dynam-ical 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][21][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 oper-ators 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  Figure 1: Bifurcation diagram of (5) due to orders from 0.93 to 1.
2 Journal of Function Spaces 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][28][29]. The manuscript is organized as follows: Section 2 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 men-tioned 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.

Statement of the Problem
The Lorenz-Stenflo (LS) system of differential equations is given by [30,31].
The variable x 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.

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 function f : ½0,+∞Þ ⟶ ℝ is defined by Definition 2 (see [2,3]). The Riemann-Liouville fractional derivative for a continuous function f : ½0,+∞Þ ⟶ ℝ is defined by Definition 3 (see [2,3]). The Caputo fractional derivative for a continuous function f : ½0,+∞Þ ⟶ ℝ is defined by 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 where t (s) 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.

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 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 E 0 = ð0, 0, 0, 0Þ, and constant values of the parameters given by k = 0:7, r = 30, c = 4, and a = 10.  Journal of Function Spaces The matrix J at E 0 = ð0, 0, 0, 0Þ is given by : ð8Þ The eigenvalues of the matrix J 0 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 E 0 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].

Numerical Solution
This section presents the numerical method used to obtain the phase portraits and time series solutions of the fractionalorder 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 apply-ing the Reimann-Liouville fractional integral given by (2), we obtain the following system Applying the predictor-corrector method to (12), the integral terms are approximated and it becomes x and Moreover, the predictors are given as follows: 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 predictorcorrector 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.

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).

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. Based on the LE values in Table 1, 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.

Journal of Function Spaces
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 by dim ðLEÞ of system (5) for η ∈ ½0:857,1 can be calculated. For instance, for the fractional order η = 0:856, the Kaplan-Yorke dimension is From the above dim ðLEÞs, 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.

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 Figures 2-4 to verify the conclusion made in Figure 1.
It can be seen in Figures 2(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.

Bifurcation for Parameter r.
To get the bifurcation diagram of parameter r 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 y 0 = ½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.
It can be inferred in Figure 6 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.

Bifurcation for Parameter c.
For the bifurcation diagram of parameter c 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 y 0 = ½ 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.
It is observable in Figure 5 that system (5) shows hyperchaotic dynamics for parameter c in [3.5, 4.5]. Moreover, the phase portraits of system (5) for c = 3:6 and 4:4 shown in Figure 10 verify the conclusion drawn in Figure 9.

Bifurcation for Parameter k.
To obtain the bifurcation diagram due to parameter k 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 y 0 = ½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.
It is clear in Figure 11 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.

Bifurcation for Parameter a.
To obtain the bifurcation diagram due to parameter a of the hyperchaotic system (5), the fractional derivative, the initial condition, the time, and the values of the remaining parameters are the same as 18 Journal of Function Spaces those used in Section 5.5. The bifurcation diagram due to the variation of parameter a is illustrated in Figure 12. It is clear in Figure 12 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).

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 ðx 0 , 0:5,5, 0:5Þ . The different values of x 0 are shown in Figure 14.
It can be observed in Figure 14 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 x 0 of two trajectories that overlapped for the first 15.1 seconds are closer to each other than the value of x 0 of the third trajectory.

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 and the slave system be given by It can be shown that all the real parts of the eigenvalues of H are negative. Thus, the Matignon criterion is satisfied; H is a Hermitian matrix. Now, considering the coupling function (23), the slave system becomes 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, and a = 10 together with the fractional order η = 0:978: Moreover, the 19 Journal of Function Spaces initial value conditions used for the master and the slave system are y 0m = ½0:5,0:5,5, 0:5 and y 0s = ½5, 5, 5, 5, respectively.
It can be observed in Figure 15 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.

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.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no competing interests.

Authors' Contributions
Chernet Tuge Deressa did the actualization, methodology, formal analysis, validation, investigation, and initial draft. Sina Etemad did the actualization, validation, methodology, formal analysis, investigation, and initial draft. Mohammed K. A. Kaabar did the actualization, methodology, formal analysis, validation, investigation, initial draft, and supervision of the original draft and editing. Shahram Rezapour did the actualization, validation, methodology, formal analysis, investigation, and initial draft. All authors read and approved the final version.