Theory and applications of new fractional-order chaotic system under Caputo operator

Article History: Received 6 May 2021 Accepted 16 August 2021 Available 27 October 2021 This paper introduces the properties of a fractional-order chaotic system described by the Caputo derivative. The impact of the fractional-order derivative has been focused on. The phase portraits in different orders are obtained with the aids of the proposed numerical discretization, including the discretization of the Riemann-Liouville fractional integral. The stability analysis has been used to help us to delimit the chaotic region. In other words, the region where the order of the Caputo derivative involves and where the presented system in this paper is chaotic. The nature of the chaos has been established using the Lyapunov exponents in the fractional context. The schematic circuit of the proposed fractional-order chaotic system has been presented and simulated in via Mutltisim. The results obtained via Multisim simulation of the chaotic circuit are in good agreement with the results with Matlab simulations. That provided the fractional operators can be applied in realworlds applications as modeling electrical circuits. The presence of coexisting attractors for particular values of the parameters of the presented fractional-order chaotic model has been studied.


Introduction
The applications of fractional operators in modeling real-world problems continue to receive many impressions. In the literature, there exist many investigations, including fractional operators in mathematical physics [1][2][3], physics [4], science and engineering [4-6] and many other domains [7]. Modeling using fractional operators is nowadays preferred to modeling with integer-order derivative because it is mentioned in the literature that fractional operators consider memory effects and generalize the classical derivative. The Caputo derivative [1,8] receive more much interest because the classical initial condition used in physics can also be used. The Riemann-Liouville derivative is not physical tools because the initial conditions of the model problem should be in integral form. It is the main inconvenience in physics to have an initial condition in integral form. In the last ten years, there exist some interest in modeling chaotic systems using fractional-order derivatives. Theory of chaos in fractional context begin to receive attractions [9][10][11]. The reasons are, with the variation of the order of the used derivative, we can conserve or lost chaotic behaviors of the systems. Many new attractors born also with the fractional operators, see in [12][13][14][15]. Another reason is the definitions used to characterize the nature of the chaos. For example, in a fractional context, there exist hyperchaotic systems with one positive Lyapunov exponent.
Before the motivations and the works described in this paper, we give a brief literature review concerning modeling chaotic systems with and without fractional-order derivative. Recently, Sene et al. introduce several fractional chaotic systems via Caputo derivative and study their phase *Corresponding Author portraits, bifurcations diagrams, Lyapunov exponents, sensitivity to the initial conditions, and influence of the fractional orders in obtaining new types of attractors, see in the following investigations [11,14,[16][17][18]. In [13], In the fractional context, Pacheco et al. introduce and investigate fractional chaotic systems with hidden and selfexcited attractors. Akgul in [12] study chaotic oscillations using fractional operator memcapacitor. In this paper [12], the auteur draws the circuit schematic associated with the proposed fractional chaotic system and simulates it and proved the experimental results are in good agreement with the result obtained with the numerical scheme. Akgul and al. have also introduced many other chaotic systems where they use fractional operators; the readers can read the following investigations by the same author [19,20]. The fractionalorder chaotic systems with non-singular kernels have been recently introduced in the following papers [10,21]. In [22], Rajagopal et al. presented a multi-scroll chaotic system with sigmoid nonlinearity and have proposed the fractional context of their model; they notably investigate with Caputo fractional derivative and the predictorcorrector method of Adams-Bashforth-Moulton to obtain the phase portraits. For more investigations by Rajagopal in the fractional chaotic systems, look to [23]. For a fractional-order chaotic system with conformable derivative, the readers can refer to Perez works in the following paper [9]. In [24], Ruiz et al. have used the variable-order fractional derivative, which is a very novel operator to model fractional chaotic system. In [25], Danca stipulated the definition saying; we have hyperchaotic behaviors with two positive Lyapunov exponents; this definition seems not adequate in fractional context because there exist fractional hyperchaotic systems with one positive Lyapunov exponent. This conclusion addressed by Danca opens new problems in modeling a chaotic system with fractional operators. See more investigations in chaotic systems in [26].
In this paper, we propose modeling a chaotic system using the fractional operator. Motivated by the initial condition, we use the Caputo derivative in our modeling. In this paper, we will focus on the influence of the Caputo derivative in the dynamics of our considered chaotic system. The main question will be to determine the chaotic region and simulate the fractional chaotic system in Multisim using the circuit schematic of our proposed chaotic system. Note that the circuit schematic will include the capacitors and resistors, and other electrical tools. After modeling, we will study the influence generated by the variation of the model's parameters via the bifurcation diagrams on the dynamics of our chaotic system. We will confirm the chaotic region by calculating the Lyapunov exponents. It will be established our considered chaotic system is dissipative as well. The local stability of the equilibrium points of our system will be illustrated via the Matignon criterion.

On fractional operators in fractional calculus
We recall the fractional operators which we will utilize in this paper. There exist many types of fractional operators which can be classified into two families; the first family is the fractional operators with the singular kernels: Caputo derivative and Riemann-Liouville derivative, and the second family is composed of the fractional operators with nonsingular kernels: the fractional derivative with exponential kernel and the fractional derivative with Mittag-Leffer kernel. Many other operators exist in fractional calculus, but here we consider the operators with a great audience and englobing all the papers in the literature. Therefore the rest of the paper, we recall the Caputo derivative, the Riemann-Liouville derivative, the Riemann-Liouville integral, and the Laplace transform of the Caputo derivative. We begin by defining the fractional integral known as the name Riemann-Liouville fractional integral, which will be used to establish the numerical scheme utilized in this paper.
where we define the function Γ(...) as the Gamma Euler function and we set also the order α > 0.
The above derivative is generalized in the following definition established in the literature by Thabet and Fahd in [7].
with the orders defined as α, and ρ satisfying in particular the condition α, ρ > 0, the Gamma function is denoted by Γ(...), and for all t > 0.
The classical Riemann-Liouville derivative is obtained when the orders of the generalized fractional derivative satisfy the condition α = ρ = 1. We now define the Caputo derivative, which will be used in modeling the chaotic system studied in this paper. The choice is due to its initial condition, which is physic, and the memory effect generated by the Caputo derivative.
The generalization of the Caputo derivative is also introduced in the literature by Thabet and Jarad, and the following definition defines it.
The Caputo derivative also admits the Laplace transform and obeys many properties, which can be found in the literature of the fractional calculus. We advise the reader to take a look at the following investigations [1,8], in details.

Fractional-order attractor system
A new fractional chaotic system is addressed in this section. The fractional differential equations satisfied by this new chaotic system using the Caputo derivative is described as the following form The chaotic system is sensitive to the initial conditions. In the present work, the chaotic behaviors are obtained using the following initial conditions Later the impact of the initial conditions in the dynamics will be focused. The strange attractor is obtained with the values of the parameters of the model described by the system (5)-(7) set as the forms a = 9, b = 4, c = 1.5, d = 2.85, and e = 0.2. The values of the parameters of the chaotic systems also play a significant role in the attractors. In other words, the variation of the model's parameters can generate the removal of the chaotic behaviors. There exist many fractional operators in fractional calculus, here we opted to use the Caputo derivative due to the fact, it takes into account the memory effect, and the derivative takes into account also the initial conditions in form (8), contrary to the Riemann-Liouville derivative where the initial conditions should be in an integral form which is irrealistic in real-world problems. Another main motivation for the use of fractional operators is, in general, in the chaotic system, they generate new types of attractors according to their variations in specific intervals. The interval where the fractional chaotic system has chaotic behaviors will be determined with the Matignon criterion in the present paper. The next section will be the subject of the numerical approximation. The numerical algorithm will be used to depict the phase portraits of the present chaotic system. The last remark in this section is that when the order α = 1, then we obtain the following chaotic system defined by the equations The same initial conditions (8) are considered for the model (9)-(11) with integer-order derivative. Model (9)-(11) will not be considered in this section. Still, it is important to mention that all classical models with an integer-order derivative can be obtained when the fractional operators converge to 1..

Numerical procedure of the fractional chaotic system
In this section, we recall the numerical discretizations which we will use to depict the phase portraits of the present model. Numerical discretization exists in the literature; here, we describe the numerical discretization for our present model. We begin with the analytical solution of the model, which can be represented using the Riemann-Liouville integral; we have the following form Where v = (x, y, z). The functions ̟, κ and ω come from the primary model described in Eq.
(5)- (7). Therefore, we have the following forms Before beginning the discretization of the Riemann-Liouville integral at the point t n , we evaluate Eqs. (12)- (14) at this point t n , we have the following representations Let the step size h, and we define t n = nh, the Riemann-Liouville integral can be written as the following equations for the functions ̟,κ and ω with the parameters given as the following forms We consider the first-order interpolant polynomial of the functions ̟ (τ ), κ (τ ) and ω (τ ) represented as the following relationships [27] Plugin Eqs. (24)-(26) into the Eqs. (21)-(23), we obtain the following representations of the discretizations of the Riemann-Liouville integrals, we have that under the parameters set explicitly as the following representations and when n = 1, 2, ...., the parameters d of the discretization can be rewritten as the following equations and Finally, under the above representations, then the numerical scheme of the fractional differential equation (5)- (7) under consideration is represented by the following relationships The numerical discretization of the functions ̟, κ and ω are represented in the following equations The convergence and the stability are not detailed in this section; see in for more pieces of information [27]; here, we give pieces of information related to these two properties. We assume that x(t n ), y(t n ) and z(t n ) are the approximate numerical solutions of the fractional system represented in Eq. (5)-(7) and x n , y n and z n are the exact solutions, the residuals functions with Caputo derivative are given by the expressions We can notice the convergence of the numerical discretization presented in this section is obtained when the step-size h converges 0 [27]. The stability of our numerical scheme is obtained from the Lipchitz criterion of the functions u, v, and w. We finish this part by illustrating the numerical scheme described in this section. We consider for graphical representations the orders α = 0.86, α = 0.94 and α = 0.97. We begin by the order α = 0.86, we depict by the phase portraits in all We have the following Figures 1a-1b,2a-2b.
Let the order of the Caputo derivative given by α = 0.94, we draw by the phase portraits in all planes (x − y − z), (x − y), (x − z) and (y − z).
We have the Figures 3a-3b, 4a-4b. We continue with the order of the Caputo derivative defined by α = 0.97, we depict the phase portraits in all planes (x − y − z), (x − y), (x − z) and (y − z). We have the following Figures 5a,5b,6a,6b.

Bifurcation maps for detecting chaotic regions
In this section, we address the chaotic regions according to the variations of the parameters of the model. For the rest of the paper, we consider the order of the fractional derivative fixed to α = 0.97. We first consider the variation of the parameters a into the interval (8,12). We consider a small stepsize h = 0.0025 to see the structure of the curve. The bifurcation map related to the parameter a is represented in the Figure 7.
The chaotic region can be observed into the interval (8.8, 12 In the second part, we consider the parameter b. We suppose it varies into the interval (1, 4.5). The considered step-size for clarity is h = 0.0025. The bifurcation map associated to the small variation of the parameter b can be observed in Figure 10. We notice the system ends the chaotic region in the interval (4, 4.2) via period-doubling bifurcation. The region (4.2, 4.5) corresponds to our present fractional model to the periodic region. The chaotic region is observed into the interval (1, 4.2). For illustrations of the chaotic behaviors, we consider b = 3.8, and we calculate the Lyapunov exponents to confirm the chaotic behaviors when the parameter b takes this value. We consider the variation of the parameter c in a small interval (1, 2). The objective is to repeat the previous method to detect the chaotic region. The bifurcation map according to the variation of the parameter c is represented in Figure13.
In Figure 13 we observe the chaotic region is into the interval (1, 1.7). After this region, we get periodic behaviors for our fractional model. For the Figure 13 the considered step-size is h = 0.005. We can also notice the chaotic behaviors of the system in the chaotic region is not uniform; in the interval (1, 1.575) and (1.575, 1.7), the behaviors of the system are different. Let us consider the variation of the parameter d in a small interval (2, 3.5). The bifurcation map    versus the variation of the parameter d is represented in Figure 14.
The first remark is there exists a minor difference in the influence of parameter c and the parameter d. In our interval of variation, the chaotic region is observed into the interval (2, 2.9). The periodic behaviors are observed in the rest of the interval.
The main interest of this section is we can observe the chaotic region can be divided into two sections. The section (2, 2.87) and (2.87, 2.9) where the dynamics of the solutions are different as it can be observed in Figure 14. For illustrations of this difference we consider the phase portrait at a = 2.65 and a = 2.875.

Co-existing attractors for the fractional model
In this section, we illustrate the existence of a pair of attractors at two different initial conditions by considering the variation of the model's parameters (5)- (7). The considered initial conditions in this section are (1, 1, 1) and (−1, −1, −1). In the literature, the notion of coexisting attractors is present when at least two initial conditions are considered, and the parameters vary into a certain interval. According to the variation of the model's parameters, when the existence of a pair of attractors exists, it can also be observed via the bifurcation maps. Therefore this section will   be consecrated to observe the coexisting attractors via the bifurcations maps. We will also represent some phase portraits to support the findings of this section. We start with the variation of the parameter a. The coexisting bifurcation map according to the variation of the parameter a into the interval (8,11) with the initial condition (1, 1, 1) and (−1, −1, −1) are represented in Figure 19. We consider large intervals to observe clearly the coexisting attractor.
Via the bifurcation map, we can observe according to the variation of the parameter a, the system admits coexisting attractors into the interval (8.8, 9.5). Into this interval, the pair of attractors can not be confused contrary to the interval  In Figures 20a-20b, 21a-21b we are in the region where the coexisting attractors are clearly visible. We mean the region (8.8, 9.5). After this interval, the coexisting attractors are confused, as can be seen in the following figures at a = 10.
Similarly, we can provide the coexisting attractors for the other parameters of the model b and c. We finish this section by providing the coexisting attractors with the parameter b. Therefore we depict the coexisting bifurcation map according to the variation of the parameter b and taken at two different initial conditions (1, 1, 1) and (−1, −1, −1). The considered interval for the

Stability analysis of the equilibrium points
In this section, we trait the stability analysis of the equilibrium points of the considered model (5)- (7). We will use the Jacobian matrix of the system to give the values of the Lyapunov exponents. For the stability, we consider the Matignon criterion to prove the local stability of all the equilibrium points of the system (5)- (7). To obtain The eigenvalues of the matrix defined by J(A) are represented as the following forms λ 1 = 1.3233 + 5.3154i, λ 2 = 1.3233−5.3154i and λ 3 = −11.2966.  We can notice the last eigenvalue satisfies the Matignon criterion given by |arg (λ 3 )| = π > απ/2 for all orders of the Caputo derivative. The first and the second eigenvalues satisfy the condition |arg (λ 1,2 )| = 19π/45 > απ/2, for all orders satisfying that α ∈ (0, 0.84). Thus the point A is not stable when the order of the Caputo derivative is into the interval (0.84, 1). We We determine the eigenvalues of the matrix J(B) given as λ 1 = 1.3441 + 5.3418i, λ 2 = 1.3441 − 5.3418i and λ 3 = −11.3382. The next step consists to verify the Matignon criterion for all our eigenvalues. Note that the third eigenvalue respect the condition |arg (λ 3 )| = π > απ/2. We can notice the first and the last eigenvalues are conjugate and satisfies the condition |arg (λ 1,2 )| = 75.9π/180 > απ/2 for the order of the Caputo derivative obeying to α ∈ (0, 0.84). Thus the point B is not stable when the order describes the interval (0.84, 1). The third part consists to investigate the local stability of the point C. The Jacobian matrix at the point C is given by the following matrix The eigenvalues of the above matrix are summarized as the values λ 1 = −7.5, λ 2 = −4.0001 and λ 3 = 2.8501. The first and the second eigenvalues satisfy the matignon criterion because they are negative. The last eigenvalue obeys to |arg (λ 3 )| = 0 > απ/2 for all order of the Caputo derivative. Thus, the point C is not stable. The final conclusion about the stability analysis the equilibrium points of the fractional differential considered in this section; they are not stable into the interval (0.84, 1) which is probably the chaotic region of our considered system. The investigation about the chaotic region will be analyzed later for confirmation.

Lyapunov exponents for chaos characterization
In this section, we calculate the Lyapunov exponents associated with the different orders used in this paper to characterize the nature of the chaos. There are many types of behaviors, but we are interested in the chaotic, hyperchaotic, limit cycles behaviors, and convergence to an equilibrium point. For brief rappel, in the context of integer order derivative, the system's behavior is considered to describe limit cycle when all the Lyapunov exponents are negative, the same characterization in the context of the convergence of the solution and equilibrium point. The chaotic behavior is obtained when one of the Lyapunov exponents is negative. The hyperchaotic context follows when at least we have two positive Lyapunov exponents. In the context of fractional order derivative, the above cite characterizations are not adequate. Danca provided this assumption. He stipulates the determination of the Lyapunov exponents in fractional context is complex, and the number of positive Lyapunov exponents is not adequate to characterize hyperchaotic behaviors. There exist hyperchaotic behaviors with one positive Lyapunov exponents; see in [25]. The method of the determination use the Jacobian matrix, and the numerical scheme, the method of the determination can be found in [28], we have the following matrix For the Lyapunov exponents to validate the chaotic behaviors of the phase portraits depicted in Section 4, we calculate the Lyapunov exponents; we consider the first the order α = 0.86. We have the following values First of all, we can observe the sum of the Lyapunov exponents is negative, which means the considered system in this study is dissipative. We can also see our system has chaotic behaviors at the order α = 0.86 because one of the Lyapunov exponents is positive. We finish by its associated Kaplan Yorke dimension given by the following We have the same analysis with the order α = 0.94. We calculate the Lyapunov exponents. We have the following numbers L1 = 0.5720, L2 = 0, L3 = −11.2740.
We also have the chaotic behaviors at the considered order because one of the Lyapunov exponents is negative. The system stays dissipative in this order because the sum of the Lyapunov exponents is negative. For the order, α = 0.94, the Kaplan Yorke dimension obtained from the Lyapunov exponents is given by the following For the order α = 0.97, we repeat the same procedure of calculations. The values of the Lyapunov exponents are given by the numbers its associated dimension is given by With Eq. (51), we conclude our system is dissipative and has chaotic behaviors. We dress Table of Lyapunov exponents for the confirmation of the chaotic behaviors into the region (0.84, 1) obtained in the stability analysis. With Table ??, we can observe the fractional system (5)-(7) considered in the region (0.84, 1) has chaotic behaviors due to the fact there are one positive Lyapunov exponents at all orders in this interval. We can also observe the system is dissipative in this region. The Table previously presented confirm that our interval obtained in the stability investigations is our chaotic region.

Circuit realization of the uncommensurate fractional-order system
We finish this paper by realizing the circuit obtained by the fractional chaotic system considered in this section. To arrive at our objective, we draw the circuit from the system (5)- (7), and we simulate the model using Multisim. The interest of this procedure is to compare the phase portraits obtained with Matlab and the phase portraits obtained after simulation in Multisim. We will observe the results are in good agreement. This section plays an important role in the real applications of fractional calculus in engineering.
For the rest of this section, we consider the order α = 0.97. First of all, we recall the third approximation of fractional integrator, which we will use; it is given by the following relationship 1 s 0.98 ≈ 1.2234s 2 + 1463.28s + 4893.2 (s + 0.0106)(s + 3.7716)(s + 1341.4) .
(53) The above approximation will be compared to the transfer function associated with the fractional integrator given by the form Comparing Eq. (53) and Eq. (54), we get the following values for the Resistors and the Capacitors C a = 1.0656nF , C b = 8.5245n, C c = 7.596n, R a = 91.17M Ω,R b = 32.046kΩ, and R c = 101.12Ω. The second part consists of drawing the circuit associated with the fractional chaotic system (5)- (7). We first use the Capacitors and Resistors to rewrite Eqs. (5)- (7) we have the following system where, we set the following values for the Capacitors C 1 = C 2 = C 3 = 1nF , and the values of the resistors are R 1 = 87.719kΩ, R 2 = 25kΩ, R 3 = 33.33kΩ, R 4 = 6.25kΩ, R 5 = 100000kΩ, R 6 = 250kΩ, R 7 = 6.25kΩ and we set V = −1V . These values are obtained by calibrating Eq. (5)-(7) and the previous Eqs. (55)-(57). The phase portraits obtained after simulation in the Multisim, we have the following outputs.   For confirmation of the results in the oscilloscopes, we consider its corresponding fractional system and depict the portraits using the numerical scheme previously exposed in this paper. Thus, we consider the following system y ′ = e − ay + xz + cy, (59) D α c z = −bz + xy.
(61) We suppose the order α = 0.98 and we depict the portraits of system (58)-(60) in the following Figures 29a-29b,30a-30b. We can finally observe the results obtained after stimulation of the fractional chaotic (55)-(57) in Multisim are in good agreement with the results established for the model (58)-(60). This section also proves fractional calculus can be used in modeling electrical chaotic circuits.

Conclusion
In this paper, we have presented the bifurcation, Lyapunov exponents to characterize the chaos, the stability analysis of the equilibrium points, and the coexisting attractors for a class of fractional chaotic system described by the Caputo derivative. To arrives at our end, the fundamental tool for the effectiveness of this paper is the numerical procedure which is crucial for obtaining the phase portraits of the chaotic system. Our fractional chaotic system is simulated by Matlab; we have also proposed the circuit schematic for our considered fractional-order chaotic system. From the fact, we have not the materials to simulate the electrical chaotic circuit, we have simulated our model virtually by considering Mutlisim software. With the simulation in Multisim, we have proved the results in Matlab are in good agreement with the results simulated in the Multisim.