# On a Memristor-Based Hyperchaotic Circuit in the Context of Nonlocal and Nonsingular Kernel Fractional Operator

Authors: Shahram Rezapour; Chernet Tuge Deressa; Sina Etemad
Journal: Journal of Mathematics
 (2021)

---

## Abstract

Memristor is a nonlinear and memory element that has a future of replacing resistors for nonlinear circuit computation. It exhibits complex properties such as chaos and hyperchaos. A five-dimensional memristor-based circuit in the context of a nonlocal and nonsingular fractional derivative is considered for analysis. The Banach fixed point theorem and contraction principle are utilized to verify the existence and uniqueness of the solution of the five-dimensional system. A numerical method developed by Toufik and Atangana is used to get approximate solutions of the system. Local stability analysis is examined using the Matignon fractional-order stability criteria, and it is shown that the trivial equilibrium point is unstable. The Lyapunov exponents for different fractional orders exposed that the nature of the five-dimensional fractional-order system is hyperchaotic. Bifurcation diagrams are obtained by varying the fractional order and two of the parameters in the model. It is shown using phase-space portraits and time-series orbit figures that the system is sensitive to derivative order change, parameter change, and small initial condition change. Master-slave synchronization of the hyperchaotic system was established, the error analysis was made, and the simulation results of the synchronized systems revealed a strong correlation among themselves.


---

## Body

## 1. Introduction

In the last decade, fractional differential equations started gaining much attention in modeling several real-world problems in different areas including in mathematical epidemiology, physics, engineering, and many others. Fractional-order operators are either with singular kernels such as the Caputo derivative and the Reimann–Liouville fractional derivatives or with nonsingular kernels such as the Caputo–Fabrizio and Atangana–Baleanu fractional derivatives [1–4].One of the differences between integer- and fractional-order derivatives is that the integer-order derivative describes local properties of a certain dynamic system, whereas the fractional-order derivative representation of a dynamic system involves the whole space of the process [5]. That is, applying fractional derivative orders in modeling real-world problems is essential for describing the hereditary specifications and effectiveness of the memory as essential aspects of different mechanisms in the problem [6, 7].The recent increase in the study of different dynamical systems using fractional-order derivatives is attributed to the fact that most of the dynamic systems associated with complex systems are found to be nonlocal involving long memory in time and intrinsically fractional derivative operators can describe such systems more accurately than the integer derivatives [8]. In other words, important features of many physical systems are best described or exposed by using fractional-order operators.There are different studies conducted using fractional derivatives for diverse types of systems pertinenet to different problems. For instance, stability analysis of fractional differential equations with unknown parameters was considered using the D-decomposition method by Koksal [9]. Time- and frequency-domain response of the RLC circuit with fractional-order derivative was investigated in [10]. Telegraph equation used for the power transmission line with a nonlocal boundary value problem was deliberated with a different type of finite difference numerical schemes in [11].Chaos theory has attracted many researchers and has applications in the fields of encryption and secure of communication [12], modeling financial systems, representing circuit diagrams, and many others [13, 14].Nowadays, there are many pieces of literature dedicated to the analysis of chaotic systems using fractional derivative operators. This includes chaotic systems of Chua’s electrical circuits and memristor-based circuit systems. Some of the literature are reviewed as follows.Sene applied the Caputo fractional derivative operator in detecting the chaotic behavior of different 3D and 4D chaotic systems. He used Lyapunov exponent and bifurcation diagrams to identify the nature of the chaos and impact of parameter variation for the different chaotic models investigated [15–18]. Different chaotic systems including Chua’s electric circuit and several other chaotic systems were analyzed using fractional derivative order mathematical models by Petráš [19]. Bifurcation and chaotic behaviors in fractional-order simplified Lorenz system using Adams–Bashforth–Moulton predictor-corrector scheme were considered in [20]. Atangana–Baleanu fractional derivative operators were used for modeling and analysis of different chaotic and hyperchaotic systems, and solutions were approximated using a two-step Adams–Bashforth numerical scheme in [21].It was in 1971 that circuit theorist Chua proposed memristor as a missing two-terminal nonlinear electrical component. The three basic components of a circuit are resistor, capacitor, and inductor. Memristor known for its memory effect and nonlinear characteristics is the fourth circuit component. Memristor relates magnetic flux and electric charge linkage in which case it is called a charge-controlled electric modelϕ=ϕq or it models a relationship between charges and flux q=qϕ in which case it is called a flux-controlled model [12, 22].At present, there are several studies conducted on memristor-based chaotic circuits using both integer and fractional-order derivatives. In [23], a conformal fractional-order simplest memristor-based chaotic circuit was investigated based on conformable Adomian decomposition method, Lyapunov exponent, bifurcation diagram, and Poincare’ sections. Buscarino et al. introduced a chaotic circuit based on a realistic model of the HP memristor, and numerical results showed a generation of chaotic attractors [22]. A novel 5D chaotic system with flux-controlled memristor and integer-order derivative, extracted from Wang’s 4D hyperchaotic system, was proposed by Wang et al. [24]. A memristor-based chaotic circuit modified by replacing the nonlinear resistor in Chua’s circuit with a flux-controlled memristor was analyzed in [12].This research aims to study the memory effect properties and detection of chaos in a five-dimensional memristor-based system. Accordingly, an integer model memristor-based circuit is represented by Atangana–Baleanu fractional derivatives in the Caputo sense (ABC), and the existence and uniqueness of solution of the ABC fractional model are analyzed based on Banach fixed point theorem for contraction principle. From the different concepts of fractional-order operators introduced above, the purpose of choosing the ABC fractional derivative is due to the fact that it possesses nonlocal kernels and of course it allows the inclusion of traditional initial conditions in the formulation of a mathematical model. Numerical approximation of the ABC fractional model is made using the newly developed numerical approximation for fractional derivative developed by Toufik and Atangana in [25]. Local stability of the fractional model representation is accomplished using the Matignon stability criterion. The existence and nature of chaos in the fractional model are investigated using Lyapunov exponents. A bifurcation diagram for different fractional derivatives and parameter variation is given. Several phase portraits are depicted as a verification for the impact of different parameter values and different fractional derivative orders. The impact of initial conditions on the solution trajectory of the system is also investigated using simulation of the trajectories of the system for different initial conditions. Lastly, master-slave synchronization of the system is performed accompanied by the corresponding error analysis and simulation of several cases of the synchronized system. All the phase portraits and solution trajectories in this work were obtained from the numerical scheme of Toufik and Atangana adapted for the memristor chaotic model considered in this work. A computing software application called Matlab 2019a is used for the simulation of different results.A Matlab code for Lyapunov exponents and bifurcation diagram of fractional-order systems called Danca algorithm [26] is used to quantify the chaos by calculating Lyapunov exponents and obtaining bifurcation diagrams for different fractional orders and different parameter values of the model. Some of the evidence for the originality of this work includes the application of Atangana–Baleanu fractional operator to the memristor-based system considered in this study, application of the newly developed numerical approximation by Toufik and Atangana for the fractional-order systems, obtaining the phase portraits of the system from the numerical scheme, and performing synchronization of the five-dimensional system using the numerical approximation.The remaining part of this paper is arranged as follows. In Section2, the memristor-based circuit model considered in this study is described in the context of integer derivative. Section 3 is devoted to the fractional derivative representation of the memristor-based systems following some recapping of preliminary concepts and definitions of Atangana–Baleanu fractional derivatives. The existence and uniqueness of the solution for the fractional derivative representation of the model are portrayed in Section 4 of the paper. The numerical scheme applied to get all the phase portraits and approximate time-series solution of the memristor-based system are developed in Section 5. Section 6 is concerned with the local stability analysis of the fractional model followed by Lyapunov exponents, bifurcation diagrams with different fractional orders, and parameter variation in Section 7. Investigation of the impact of small change in the initial conditions on the dynamics of the system is considered in Section 8 followed by synchronization of the hyperchaotic system in Section 9. Finally, in Section 10, conclusions are given followed by list of references.


## 2. Mathematical Model Description of Memristor-Based Circuit

Wang et al. [27] developed a new four-dimensional hyperchaotic circuit system using an improved modularized design and proposed implementation of the circuit. The hyperchaotic circuit system mentioned in [24] is described under an integer-order system of differential equations. Later on, Wang et al. [24] proposed a five-dimensional flux-controlled memristor-based circuit system and demonstrated that the system exhibits hyperchaotic character under an integer-order system of differential equations. The five-dimensional nondimensionalized flux-controlled memristor-based circuit developed by Wang et al. is described using integer-order differential as given in the following equation:(1)dxdt=β1y−x+4yz−0.02xtWu,dydt=−x+16y−xz+w,dzdt=−β2z+xy−xu−yw,dwdt=−10y+0.15xz−0.3zu,dudt=−x,where the memductance used in this study is(2)Wu=a+3bu2,a>0,b>0.In this study, motivated by the contribution of Wang et al. [24] in developing the five-dimensional memristor-based circuit, we analyzed system (1) under ABC fractional derivatives.


## 3. Fractional Derivative Representation of the Model

In this section, we recall the definitions and basic properties of nonsingular and nonlocal Atangana–Baleanu fractional derivative of Caputo type (ABC) used for our analysis of the circuit.Definition 1.
Letf∈C1a,b,a<b, and let q∈0,1. The AB fractional derivative of order q is defined as [7, 28, 29](3)DaABCtqft=Fq1−q∫atdfdkEq−q1−qt−kqdk,where Fq=1−q−q/Γq and the Mittag-Leffler function is(4)Eqz=∑β=0∞zβΓ1+qβ,q,z∈ℂ,ℜq>0.Definition 2.
The Atangana–Baleanu (AB) fractional integral of the functionf∈C1a,b,a<b, is given by [7, 28, 29](5)IaABtqft=1−qFqft+qFqΓq∫atfkt−kq−1dk.Lemma 1.
The AB fractional derivative and AB fractional integral off∈C1a,b,a<b, fulfill [30](6)IaABtqDaABCtqft=ft−fa.Lemma 2.
Forf,g∈C1a,b,b>a, the AB fractional derivative in the Caputo sense satisfies the Lipschitz condition [29]:(7)DaABtqft−DaABtqgt≤Λft−gt.Now we continue with the reformulation of (1) in terms of ABC fractional derivatives.Lemma 3.
LetGt∈H1a,b. Then, D0ABCtqft=Gt,f0=G0, admits one solution given by(8)ft=G0+1−qFqGt+qFqΓq∫0tt−kq−1Gkdk.We can describe the result in Lemma3 in the form of the Banach fixed theorem point as follows. We begin by defining a Banach space(9)X=x∈CJ,ℝ,J=0,1,with a norm defined as xx=Supt∈Jxt.Now define the operatorsM1,M2,M3,M4,M5:X⟶X by(10)M1xt=x0+1−qFqβ1yt−xt+4ytzt−0.02xtWut+qFqΓq∫0tt−kq−1β1y−x+4yz−0.02xtWukdk,(11)M2yt=y0+1−qFq−xt+16yt−xtzt+wt+qFqΓq∫0tt−kα−1−xk+16yk−xkzk+wkdk,(12)M3zt=z0+1−qFq−β2zt+xtyt−xtut−ytwt+qFqΓq∫0tt−kq−1−β2zk+xkyk−xkuk−ykwkdk,(13)M4wt=w0+1−qFq−10yt+0.15xtzt−0.3ztut+qFqΓq∫0tt−kq−1−10yk+0.15xkzk−0.3zkukdk,(14)M5ut=z0+1−qFq−xt+qFqΓq∫0tt−kq−1−xkdk.We are now ready to describe the dynamic equation for the memristor-based chaotic circuit given in (1) using ABC fractional derivative:(15)D0ABCtqxt=β1y−x+4yz−0.02xtWu,D0ABCtqyt=−x+16y−xz+w,D0ABCtqzt=−β2z+xy−xu−yw,D0ABCtqwt=−10y+0.15xz−0.3zu,D0ABCtqut=−x.The initial condition is given asx0,y0,z0,w0,u0=x0,y0,z0,w0,u0 and Wu=a+3bu2.


## 4. Existence Theory on the Model

Here, the existence and uniqueness of the solution for the ABC model given in (15) are shown using Banach fixed point theorem for contraction mapping. The following two theorems are worth recalling before proceeding further.Theorem 1.
LetΠ be any nonempty closed subset of a Banach space X. Then, any contraction M:Π⟶Π has a unique fixed point [31, 32].Theorem 2.
Assume thatx,y,z,w,u are continuous functions satisfying the following conditions:(C1) There exist constantsk1,k2,k3,k4,k5>0 such that xt<k1,yt<k2,zt<k3,wt<k4,ut<k5.(C2)Γq1−q+1β1+0.02W<FqΓq, 161−qΓq+16<FqΓq, and Γq1−q+1β2<FqΓq; then, the ABC fractional derivative system given by (15) has a unique solution in the region X.Proof.
Let us show that the operatorM1 defined in (10) is well defined in the sense that M1xt∈ℑr and D0ABCtqM1xt is continuous on J=0,1, where(16)ℑr=B∈X,B≤r,r≥0,for β1k2+k1+4k2k3+0.02k1a+3bk52x0XFqΓq+Γq1−q+1/FqΓq<r.
Now for anyx∈ℑr and from (10), we have(17)M1x=x0X+1−qFqsupt∈Jβ1yt−xt+4ytzt−0.02xtWut+qFqΓqsupt∈J∫0tt−kq−1β1y−x+4yz−0.02xtWukdk≤x0+β1k2−k1+4k2k3−0.02k1a+3bk521−qFq+qFqΓqsupt∈J∫0tt−kq−1dk≤x0FqΓq+Γq1−q+1FqΓqβ1k2+k1+4k2k3+0.02k1a+3bk52<r.
Consequently, we haveM1xt∈ℑr.
To show continuous differentiability onJ=0,1 we proceed from D0ABCtqM1xt=D0ABCtqI0ABtqβ1y−x+4yz−0.02xtWuβ1y−x+4yz−0.02xtWu which is continuous on J=0,1, and then we conclude that D0ABCtqM1xt is continuous on J; as a result, M1ℑr⊂ℑr.
To show that the operatorM1 has a fixed point based on Theorem 1, it is enough to show that M1 is a contraction mapping. Indeed, let x1,x2∈X,t∈J, and it is not difficult to show that(18)M1x1−M1x2X≤Γq1−q+1FqΓqβ1+0.02WuXx2−x1X≤Hx2−x1X,where H=Γq1−q+1β1+0.02WuX/FqΓq. Since H<1 by hypothesis of Theorem 2, we conclude that M1 is a contraction mapping.
To show that the remaining operatorsM2,M3,M4, and M4 are contraction mappings, we proceed in the same manner.
Let us show that the operatorM2 defined in (14) is well defined in the sense that M2yt∈ℑ′ and D0ABCtqM2yt is continuous on J=0,1, where(19)ℑr′=B∈X,B≤r′,r′≥0,for r′>k1+16k2+k1k3+k4y0XFqΓq+Γq1−q+1/FqΓq.
Now for anyy∈ℑr′ and from (11), we have(20)M2ytx=y0X+1−qFqsupt∈J−xt+16yt−xtzt+wt+qFqΓqsupt∈J∫0tt−kq−1−xk+16yk−xkzk+wkdk≤k1+16k2+k1k3+k4y0XFqΓq+Γq1−q+1FqΓq≤r′.
Consequently, we haveM2yt∈ℑr′. To show continuous differentiability on J=0,1, we proceed from D0ABCtqM2yt=D0ABCtqI0ABtq−x+16y−xz+w=−x+16y−xz+w which is continuous on J=0,1, and we conclude that D0ABCtqM2yt is continuous on J; as a result, D0ABCtqM2yt.
To show that the operatorM2 has a fixed point, we apply Theorem 1. Following the theorem, it is sufficient to show that M2 is a contraction mapping. Indeed, let y1,y2∈X,t∈J. Then,(21)M2y2−M2y1X≤1−qFq16y2−y1+q16y2−y1FqΓqsupt∈J∫0tt−kq−1dk≤161−qΓq+16FqΓqy2−y1≤H′y2−y1X,where H′=161−qΓq+16/FqΓq. Since H′<1, we conclude that M2 is a contraction mapping.
Let us show that the operatorM3 defined in (12) is well defined in the sense that D0ABCtqM2yt=D0ABCtqI0ABtq−x+16y−xz+w and D0ABCtqM3zt is continuous on J=0,1, where(22)ℑr″″=u∈X,u≤r″,r″≥0,for r″>β2k3+k1k2+k1k5+k2k4z0FqΓq+Γq1−q+1/FqΓq.
Now for anyz∈ℑr″″ and from (12), we have(23)M3ztX=z0X+1−qFqsupt∈J−β2zt+xtyt−xtut−ytwt+qFqΓqsupt∈J∫0tt−kq−1−β2zk+xkyk−xkuk−ykwkdk≤z0X+1−qFq−β2k3+k1k2−k1k5−k2k4+qFqΓq−β2k3+k1k2−k1k5−k2k4supt∈J∫0tt−kq−1dk≤β2k3+k1k2+k1k5+k2k4z0XFqΓq+Γq1−q+1FqΓq<r″.
Consequently, we haveM3zt∈ℑr″″. To show continuous differentiability on J=0,1, we proceed from D0ABCtqM3zt=D0ABCtqI0ABtq−β2z+xy−xu−yw=−β2z+xy−xu−yw which is continuous on J=0,1, and we conclude that D0ABCtqM3zt is continuous on J; as a result, M3ℑr″″⊂ℑr″″.
To show that the operatorM3 has a fixed point, we apply Theorem 1. Following the theorem, it is sufficient to show that M3 is a contraction mapping. Indeed, let z1,z2∈X,t∈J. Then,(24)M3z1−M3z2X≤1−qFqβ2z2−z1X+qFqΓqβ2z2−z1Xsupt∈J∫0tt−kq−1dk≤Γq1−q+1β2FqΓqz2−z1X≤H″z2−z1X,where H″=Γq1−q+1β2/FqΓq. Since by hypothesis H″<1, we conclude that M3 is a contraction mapping.
We have then proved that the operatorsM1,M2, and M3 are well defined and are contraction mappings. The case for the operators M4 and M5 follows immediately. Hence, by the Banach fixed point theorem, system (17) has a unique solution in X.


## 5. Numerical Solutions

In this part, the numerical scheme applied to get the phase portrait of the dynamic system (15) is introduced. In the context of chaotic or hyperchaotic fractional differential equations, the use of analytical methods such as the Sudumu transform method, the Laplace transform method, the homotopy analysis method, and the homotopy perturbation method cannot easily be applied because of the nonlinearities of the system [13]. This leads to the need for using numerical methods to approximate the solutions of systems of fractional differential equations. Some of the numerical methods that can be applied for this case include Adams–Bashforth and Toufik–Atangana numerical schemes [8]. Both of these methods are based on Lagrange interpolation polynomials. In this study, the newly developed numerical approximation for fractional derivatives by Toufik and Atangana is employed. The numerical scheme is particularly developed for approximation of Atangana–Baleanu fractional derivative considered in this study, and it is proved to be convergent, stable, and consistent [8].For convenience, let us write (15) in the following form:(25)D0ABCtqx=N1t,x,y,z,w,u,D0ABCtqy=N2t,x,y,z,w,u,D0ABCtqz=N3t,x,y,z,w,u,D0ABCtqw=N4t,x,y,z,w,u,D0ABCtqu=N5t,x,y,z,w,u,where(26)N1t,x,y,z,w,u=β1y−x+4yz−0.02xtWu,N2t,x,y,z,w,u=−x+16y−xz+w,N3t,x,y,z,w,u=−β2z+xy−xu−yw,N4t,x,y,z,w,u=−10y+0.15xz−0.3zu,N5t,x,y,z,w,u=−x.Now from Lemma3 and the first equation of (25), we have(27)D0ABCtqxt=N1t,xt,yt,zt,ϕt,x0=x0.The solution for (27) is given as follows:(28)xt=x0+1−qFqN1t,xt+qFqΓq∫0tN1k,xkt−kq−1dk.Applying Lagrange’s interpolation polynomial ontk,tk+1 to equality N1s,xs,ys,zs,ϕs=β1y−x+4yz−0.02xtWu leads to(29)xk≈1hs−tk−1N1tk,xtk,ytk−s−tkN1tk−1,xtk−1,ytk−1,where h=tk−tk−1.Substituting (28) into (29), we obtain(30)xtn+1=x0+1−qFqN1tk,xtk,ytk+qFqΓq∑j=1nN1tj,xtj,ytjhϒj−1−N1tj−1,xtj−1,ytj−1hϒj,where(31)ϒj−1=∫tjtj+1s−tj−1tn+1−sq−1ds=−1qtj+1−tj−1tn+1−tj+1q−tj−tj−1tn+1−tjq−1qq+1tn+1−tj+1q+1tn+1−tj+1q−tn+1−tjq+1,(32)ϒj=∫tjtj+1y−tj−1tn+1−sq−1ds=−1qtj+1−tj−1tn+1−tj+1q−1qq+1tn+1−tj+1q+1−tn+1−tjq+1.Substitutingtj=jh in (31) and (32) results in(33)ϒj−1=hq+1qq+1n+1−jqn−j+2+q−n−jqn−j+2+2q,(34)ϒj=hq+1qq+1n+1−jq+1−n−jqn−j+1+q.The expression (30) can be expressed in terms of (31) and (34) as shown in the following equation:(35)xtn+1=x0+1−qFqN1tn,xtn,ytn+qFqΓq∑j=1nN1tj,xtj,ytjΓq+2×hqn+1−jqn−j+2+q−n−jqn−j+2+2q−N1tj−1,xtj−1,ytj−1Γq+2×hqn+1−jq+1−n−jqn−j+1+q.Similarly, we get the following equations for the rest of the state variables:(36)ytn+1=y0+1−qFqN2tn,xtn,ytn,ztn+qFqΓq∑j=1nN2tj,xtj,ytj,ztjΓq+2×hqn+1−jqn−j+2+q−n−jqn−j+2+2q−N2tj−1,xtj−1,ytj−1,ztj−1Γq+2×hqn+1−jq+1−n−jqn−j+1+q,(37)ztn+1=z0+1−qFqN3tn,ytn,ztn+qFqΓq∑j=1nN3tj,ytj,ztjΓq+2×hqn+1−jqn−j+2+q−n−jqn−j+2+2q−N3tj−1,ytj−1,ztj−1Γq+2×hqn+1−jq+1−n−jqn−j+1+q,(38)wtn+1=w0+1−qFqN4tn,ytn,ztn+qFqΓq∑j=1nN4tj,ytj,ztjΓq+2×hqn+1−jqn−j+2+q−n−jqn−j+2+2q−N4tj−1,ytj−1,ztj−1Γq+2×hqn+1−jq+1−n−jqn−j+1+q,(39)utn+1=u0+1−qFqN5tn,ytn,ztn+qFqΓq∑j=1nN5tj,ytj,ztjΓq+2×hqn+1−jqn−j+2+q−n−jqn−j+2+2q−N5tj−1,ytj−1,ztj−1Γq+2×hqn+1−jq+1−n−jqn−j+1+q.


## 6. Local Stability Analysis

In this section of the study, the local stability analysis of the fractional model represented in (25) is performed. It is known that the equilibrium points of chaotic systems are not generally stable. Some of the standard methods of stability analysis in fractional calculus are the Matignon criterion and the Laplace transform methods. In this work, the Matignon method is used for its simplicity and is most commonly used in the literature for the same purpose [1, 7, 13, 17].The Matignon criterion is given by(40)argλJ>qπ22,where J represents the Jacobian matrix, λJ is the set of the eigenvalues of J, and q is the fractional-order derivative. In the context of fractional derivative, an equilibrium point of (25) is said to be locally stable provided that the Matignon criterion (40) is satisfied for each of the eigenvalues of the Jacobian matrix.To determine if the equilibrium points of (25) are stable or not, we proceed as follows:(i)
The equilibrium points: the equilibrium points of (25) are given as 0,0,0,0,0, and line equilibrium is given by α0,0,0,0,1:α∈ℝ.(ii)
The Jacobian matrix: in this study, analysis is made regarding the trivial equilibrium point for simplicity. The parameter values used, unless otherwise mentioned, are set asa=0.1,b=0.01, β1=30,andβ2=8. Thus, the Jacobian matrix of (25) evaluated at the trivial equilibrium point is given by(41)Jo=−β1−0.002β1000−11601000−β2000−10000−10000.The eigenvalues of the Jacobian matrixJ are given by λ1=0, λ2=−29.3453, λ4=0.6981, λ4=0.6981, and λ5=−8.0000. To check the local stability of the trivial equilibrium point, using the Matignon criteria, we need to verify if all the eigenvalues satisfy condition (40). Indeed, argλ1,3,4=0>qπ/2 which is not possible since q∈0,1. argλ2,5=3.1416=π>qπ/2 is true for all q∈0,1. It can then be concluded that the equilibrium point Eeqpts=0,0,0 is locally unstable for the parameter values considered.Furthermore, since one of the eigenvalues of the Jacobian matrix is with positive real part, it can be inferred that system (25) satisfies the necessary condition for showing the double scroll attractor [19].


## 7. Lyapunov Exponents, Bifurcation, and Chaos with Different Fractional Orders ofq and Different Parameter Values

In this section, the level of chaos in system (25) is quantified using the Lyapunov exponent method. Bifurcation diagrams of the system (25) related to the fractional derivative order q and three of the parameters in the model named, β1,β2, and β3, are depicted.A Matlab code for Lyapunov exponents of fractional-order systems named the Danca algorithm [26] is used to quantify the chaos by calculating Lyapunov exponents for different fractional orders of model (25). The initial conditions used in this part of the work are given by 0.11,0.11,0.11,0.11. The parameter values used are similar to the ones used above for calculating the Jacobian matrix. The corresponding Lyapunov exponents (LEs) for different fractional orders q=0.94,0.96,0.98,1.00 are shown in Table 1.Table 1
LEs corresponding to different fractional orders of model (25) for a simulation time of 300 s.
q0.940.960.980.991LE119.832317.916316.185515.383714.6215LE20.93870.84900.76760.72980.6938LE30.00500.00500.00500.00500.0050LE4−10.7674−9.7560−8.8353−8.4066−7.9978LE5−38.8609−35.3857−32.1742−30.6645−29.2172Sum of all LEs−28.8523−26.3714−24.0496−22.9526−21.8497It is then possible to conclude that the 5D system considered in this study is dissipative since the sum of the LEs in each column of Table1 is negative and the system exhibits a hyperchaotic behavior since there are at least two positive LEs in each column of the table and the largest LE is positive. Moreover, Kaplan–Yorke dimension coresponding to the fractional derivatives considered in Table 1can be calculated. For instance, the dimension of two of the fractional orders is given as follows.Forq = 0.99,(42)dimLE=4+15.3837+0.7298+0.005−8.4128−30.6645=4.2511.Forq = 0.98,(43)dimLE=4+16.1855+0.7676+0.005−8.8488−32.1742=4.2520.


### 7.1. Bifurcations due to Variation of the Fractional Orderq

For obtaining bifurcation diagrams due to the variation of the fractional orderq, the values of all the parameters are kept fixed and the order of the fractional derivative q is varied in the interval 0.8,1 with an increment of 0.001. The other parameter values used in this simulation are a=0.1,b=0.01,β1=30,andβ2=8, and the initial condition is 0.11,0.11,0.11,0.11,0.11. The bifurcation diagram is shown in Figure 1.Figure 1
Bifurcation diagram for variation of the fractional derivative orderq in the interval (0.8, 1).As shown in Figure1, when q∈0.8,1, the system is hyperchaotic. As a verification of the observation, the phase portraits corresponding to q = 0.98 and q = 0.99 using the numerical approximations for fractional-order systems indicated by equations (35)–(39) are shown in Figures 2 and 3, respectively. In both Figures 2 and 3, the dependence of the hyperchaotic system on the fractional derivative is observable as the two figures are seen to be different from each other for a fractional order difference of 0.01. If we consider the trajectories in Figure 3, the orbit due to q = 0.98 converged in less than 13.5 seconds which is not the case for q = 0.99. That means different fractional derivative orders generate reasonably different dynamics of the system (25).Figure 2
Phase portraits of system (25) for derivative order of q = 0.98 projected on different planes. (a) x-y plane. (b) x-z plane. (c) y-z plane.
(a)(b)(c)Figure 3
Phase portraits of system (25) for derivative order of q = 0.99 projected on different planes. (a) x-y plane. (b) x-z plane. (c) y-z plane.
(a)(b)(c)


### 7.2. Bifurcations due to Variation of the Parameterβ1

The values of the parameters used area=0.1,b=0.01,andβ2=8, the fractional derivative order is q=0.99, and the initial condition is 0.11,0.11,0.11,0.11,0.11. The value of the parameter β1 is made to vary in the interval (29.5, 30.5) with an increment of 0.001. The bifurcation diagram is shown in Figure 4. It can be observed from the diagram that system (25) exhibits significant hyperchaos throughout the interval under consideration. In support of this conclusion, some of the phase portraits of hyperchaotic system (25) are depicted in Figure 5 projected on different planes.Figure 4
Bifurcation diagram for variation of the parameterβ1 in the interval (29.5, 30.5).Figure 5
Phase portraits of the system (25) for different values the parameter β1. (a) x-y plane. (b) x-u plane. (c) y-z plane. (d) z-u plane.
(a)(b)(c)(d)It can be inferred from Figures5(a)–5(d) that the orbits seem to overlap initially, but as time increases, their divergence from each other increases. That is, as time goes, the orbit due to β1=30.5 showed more contraction than the trajectory due to β1=30; this is because of the sum of all the Lyapunov exponents (LEs) of the respective parameter values considered in this case. The sum of LEs for β1=30 is −22.9526 (the LEs are shown in Table 1), the sum of LEs for β1=30.5 is −23.4693 (LEs are 15.3796, 0.7298, 0.0051, −8.4066, and −31.1772), and the sum of LEs for β1=30.25 is −23.2106 (LEs are 15.3817, 0.7298, 0.0051, −8.4066, and −30.9209). Since the system under investigation is dissipative, the contraction rate of volumes due to β1=30.5 is stronger than the contraction rate of volumes due to β1=30 which is inferred from the magnitude of the sum of all LEs for the corresponding parameter values.


### 7.3. Bifurcations due to Variation of the Parameterβ2

The values of the parameters used area=0.1,b=0.01,andβ2=8, the fractional derivative order is q=0.99, and the initial condition is 0.11,0.11,0.11,0.11,0.11. The value of the parameter β2 is made to vary in the interval [7, 8.5] with an increment of 0.001. The bifurcation diagram is shown in Figure 6. It can be observed from the diagram that system (25) exhibits significant hyperchaos throughout the interval under consideration. In support of this conclusion, some of the phase portraits of the hyperchaotic system (25) are depicted in Figure 7 projected on different planes and some of the time-series trajectories are depicted in Figure 8.Figure 6
Bifurcation diagram system (25) for the parameter β2 in the interval [7, 8.5].Figure 7
Phase portraits of the system (25) for different values parameter of β2. (a) x-y plane. (b) x-z plane. (c) y-z plane.
(a)(b)(c)Figure 8
Time series solutions of the state variablesy(t) and u(t) of system (25) for different parameter values of β2.
(a)(b)It can be inferred from Figures7(a)–7(c) that the orbits seem to overlap initially, but as time increases, their divergence from each other increases. That is, as time goes, the orbit due to β1=30.5 showed more contraction than the trajectory due to β1=30, as can be observed from the sum of all LEs of the respective parameter values considered in this case. The sum of LEs for β2=8 is −22.9526 (the LEs are shown in Table 1) and the sum of LEs for β2=7.5 is −22.4277 (LEs are 15.3837, 0.7298, −7.8817, 0.0050, and −30.6645). Since the system under investigation is dissipative, the contraction rate of volumes due to β2=8 is stronger than the contraction rate of volumes due to β1=7.5 which is inferred by the magnitude of the sum of all LEs corresponding to the parameter values.The trajectories in Figure8 seem to overlap for the first few seconds and then begin diverging from each other; the divergence increases with time.


## 7.1. Bifurcations due to Variation of the Fractional Orderq

For obtaining bifurcation diagrams due to the variation of the fractional orderq, the values of all the parameters are kept fixed and the order of the fractional derivative q is varied in the interval 0.8,1 with an increment of 0.001. The other parameter values used in this simulation are a=0.1,b=0.01,β1=30,andβ2=8, and the initial condition is 0.11,0.11,0.11,0.11,0.11. The bifurcation diagram is shown in Figure 1.Figure 1
Bifurcation diagram for variation of the fractional derivative orderq in the interval (0.8, 1).As shown in Figure1, when q∈0.8,1, the system is hyperchaotic. As a verification of the observation, the phase portraits corresponding to q = 0.98 and q = 0.99 using the numerical approximations for fractional-order systems indicated by equations (35)–(39) are shown in Figures 2 and 3, respectively. In both Figures 2 and 3, the dependence of the hyperchaotic system on the fractional derivative is observable as the two figures are seen to be different from each other for a fractional order difference of 0.01. If we consider the trajectories in Figure 3, the orbit due to q = 0.98 converged in less than 13.5 seconds which is not the case for q = 0.99. That means different fractional derivative orders generate reasonably different dynamics of the system (25).Figure 2
Phase portraits of system (25) for derivative order of q = 0.98 projected on different planes. (a) x-y plane. (b) x-z plane. (c) y-z plane.
(a)(b)(c)Figure 3
Phase portraits of system (25) for derivative order of q = 0.99 projected on different planes. (a) x-y plane. (b) x-z plane. (c) y-z plane.
(a)(b)(c)


## 7.2. Bifurcations due to Variation of the Parameterβ1

The values of the parameters used area=0.1,b=0.01,andβ2=8, the fractional derivative order is q=0.99, and the initial condition is 0.11,0.11,0.11,0.11,0.11. The value of the parameter β1 is made to vary in the interval (29.5, 30.5) with an increment of 0.001. The bifurcation diagram is shown in Figure 4. It can be observed from the diagram that system (25) exhibits significant hyperchaos throughout the interval under consideration. In support of this conclusion, some of the phase portraits of hyperchaotic system (25) are depicted in Figure 5 projected on different planes.Figure 4
Bifurcation diagram for variation of the parameterβ1 in the interval (29.5, 30.5).Figure 5
Phase portraits of the system (25) for different values the parameter β1. (a) x-y plane. (b) x-u plane. (c) y-z plane. (d) z-u plane.
(a)(b)(c)(d)It can be inferred from Figures5(a)–5(d) that the orbits seem to overlap initially, but as time increases, their divergence from each other increases. That is, as time goes, the orbit due to β1=30.5 showed more contraction than the trajectory due to β1=30; this is because of the sum of all the Lyapunov exponents (LEs) of the respective parameter values considered in this case. The sum of LEs for β1=30 is −22.9526 (the LEs are shown in Table 1), the sum of LEs for β1=30.5 is −23.4693 (LEs are 15.3796, 0.7298, 0.0051, −8.4066, and −31.1772), and the sum of LEs for β1=30.25 is −23.2106 (LEs are 15.3817, 0.7298, 0.0051, −8.4066, and −30.9209). Since the system under investigation is dissipative, the contraction rate of volumes due to β1=30.5 is stronger than the contraction rate of volumes due to β1=30 which is inferred from the magnitude of the sum of all LEs for the corresponding parameter values.


## 7.3. Bifurcations due to Variation of the Parameterβ2

The values of the parameters used area=0.1,b=0.01,andβ2=8, the fractional derivative order is q=0.99, and the initial condition is 0.11,0.11,0.11,0.11,0.11. The value of the parameter β2 is made to vary in the interval [7, 8.5] with an increment of 0.001. The bifurcation diagram is shown in Figure 6. It can be observed from the diagram that system (25) exhibits significant hyperchaos throughout the interval under consideration. In support of this conclusion, some of the phase portraits of the hyperchaotic system (25) are depicted in Figure 7 projected on different planes and some of the time-series trajectories are depicted in Figure 8.Figure 6
Bifurcation diagram system (25) for the parameter β2 in the interval [7, 8.5].Figure 7
Phase portraits of the system (25) for different values parameter of β2. (a) x-y plane. (b) x-z plane. (c) y-z plane.
(a)(b)(c)Figure 8
Time series solutions of the state variablesy(t) and u(t) of system (25) for different parameter values of β2.
(a)(b)It can be inferred from Figures7(a)–7(c) that the orbits seem to overlap initially, but as time increases, their divergence from each other increases. That is, as time goes, the orbit due to β1=30.5 showed more contraction than the trajectory due to β1=30, as can be observed from the sum of all LEs of the respective parameter values considered in this case. The sum of LEs for β2=8 is −22.9526 (the LEs are shown in Table 1) and the sum of LEs for β2=7.5 is −22.4277 (LEs are 15.3837, 0.7298, −7.8817, 0.0050, and −30.6645). Since the system under investigation is dissipative, the contraction rate of volumes due to β2=8 is stronger than the contraction rate of volumes due to β1=7.5 which is inferred by the magnitude of the sum of all LEs corresponding to the parameter values.The trajectories in Figure8 seem to overlap for the first few seconds and then begin diverging from each other; the divergence increases with time.


## 8. Impact of Initial Condition

In this section, the impact of different initial conditions on the dynamics of the system (25) is addressed. It is well known that one of the properties of chaotic systems is sensitivity to initial conditions, and thus it seems relevant to verify the impact of applying different initial conditions on the phase-space and time-series solutions of system (25). Accordingly, the parameter values and derivative order used are a=0.1,b=0.01,β1=30,β2=8, and q=0.99. In Figure 9, some of the phase portraits of hyperchaotic system (25) are depicted corresponding to initial conditions (0.11,0.11,0.11,0.11.0.11) by varying initial coordinates of u0=0.5, and1.0is made to run for 5 seconds with a time step of h = 0.001. The simulation result is shown in Figure 9.Figure 9
Some of the time series orbits of the system (25) due to small changes in the initial conditions.
(a)(b)As can be observed from the figures, hyperchaotic system (25) is sensitive to changes in initial conditions. As the variation in the initial condition increases, the dynamics of the system become different from each other at least for this example.As can be observed from Figures9 and 10, a small difference of the initial condition generates a significantly observable change in the dynamics of system (25).Figure 10
Some of the Phase portraits of the system (25) due to small changes in the initial condition. (a) z-y-u space. (b) y-u plane.
(a)(b)


## 9. Synchronization of the Hyperchaotic Model

Synchronization between two chaotic systems is one of the most interesting phenomena in the study of dynamic systems. Synchronization is an occurrence in which two or more chaotic/hyperchaotic systems express a strong correlation among themselves. It is very interesting to see that two hyperchaotic systems, being highly sensitive to initial conditions and exponential divergence of nearby orbits, get synchronized starting from two different initial conditions. Pecora and Caroll [33, 34] in 1990 showed for the first time the possibility of synchronizing two chaotic systems starting from two different initial conditions. Pecora and Caroll used a technique called replacement synchronization. Chaotic synchronization has a range of applications including securing communications based on phase synchronization and controlling insulin production by beta cells. In this study, we applied the master-slave synchronization technique. That is, we considered two identical copies of system (25) and related them with a coupling function. We considered different initial conditions for the coupled systems, showed that the error dynamics are asymptotically stable, and portrayed some of the simulation results. The coupling function is added to the slave system to make it respond to the master system.Now the master hyperchaotic fractional model is given by(44)D0ABCtqxt=30y−x+4yz−0.002x+0.0006xu2,D0ABCtqyt=−x+16y−xz+w,D0ABCtqzt=−8z+xy−xu−yw,D0ABCtqwt=−10y+0.15xz−0.3zu,D0ABCtqut=−x.Let the slave model be given by(45)D0ABCtqxst=30ys−xs+4yszs−0.002xs+0.0006xsus2+c1,D0ABCtqyst=−xs+16ys−xszs+ws+c2,D0ABCtqzst=−8zs+ysxs−xsus−ysws+c3,D0ABCtqwst=−10ys+0.15xszs−0.3uszs+c4,D0ABCtqust=−xs+c5.Define the error terms by(46)e1=xs−x,e2=ys−y,e3=zs−z,e4=ws−w,e5=us−u.The error dynamics are then given by(47)D0ABCtqe1t=β1e2−β1e1+4e3e2+y+e2z−0.002e1+0.0006−e1e52+2e1e5u−u2e1−xe52+2e5ux+c1,D0ABCtqe2t=−e1+16e2+−e1e3−e1z−xe3+e4+c2,D0ABCtqe3t=−β2e3+e1e2+e1y+xe2+−e1e5−e1u−xe5+−e2e4−e4y−e2w+c3,D0ABCtqe4t=−10e2+0.15e1e3+e1z+xe3+0.3−e5e4−e5z−ue4+c4,D0ABCtqe5t=−e1+c5.The vectorc1,c2,c3,c4,csT is chosen as(48)c1c2c3c4c5=−β1e2−β1e1+4e3e2+y+e2z−0.002e1+0.0006−e1e52+2e1e5u−u2e1−xe52+2e5ux,−e1+16e2+−e1e3−e1z−xe3,−β2e3+e1e2+e1y+xe2+−e1e5−e1u−xe5+−e2e4−e4y−e2w,−10e2+0.15e1e3+e1z+xe3+0.3−e5e4−e5z−ue4,−e1,+H.ζ,where(49)H.ζ=−30−14.5000.5−14.5−2104.5000−130004.50−500.5000−5e1e2e3e4e5.The choice of the matrixH is arbitrary except that all of its eigenvalues must satisfy the Matignon criteria of stability (40).Then, the error dynamics become(50)D0ABCtqe1D0ABCtqe2D0ABCtqe3D0ABCtqe4D0ABCtqe5=−30−14.5000.5−14.5−2104.5000−130004.50−500.5000−5e1e2e3e4e5.The eigenvalues of matrix H are −40.8876, −13.0000, −12.1143, −4.9900, and −3.0081, and all the eigenvalues satisfy Matignon criteria, and thus the error dynamics are asymptotically stable.Now the slave system is given by(51)D0ABCtqxst=44.5y−0.002x+4yz−0.5u−0.0006xu2−30xs−14.5ys+0.5us,D0ABCtqyst=13.5x+37y−xz−3.5w−14.5xs−21ys+4.5ws,D0ABCtqzst=7z+xy−xu−yw−13zs,D0ABCtqwst=5w+0.15xz−14.5y−0.3zu+4.5ys−5ws,D0ABCtqust=−1.5x+5u+0.5xs−5us.The numerical simulation that verified a strong correlation between the master and the slave system in the form of phase space, time-series orbits, and the error dynamics are depicted in Figures11–13. The numerical approximation for fractional order system (25) given in equations (35)–(39) is adapted to simulate each of the error dynamics (50), the master system (44), and slave system (51). The simulation results are shown in Figures 11–13. The parameter values and derivative order used are a=0.1,b=0.01,β1=30,β2=8, and q=0.99. The initial condition used for master systems is (0.11, 0.11, 0.11, 0.11, 0.11), and it is (−10, 10, −5, −0.5, −0.5) for the slave system. The error graphs indicate fast synchronization of the slave and master systems.Figure 11
Error dynamics of the synchronization process.Figure 12
Some of the time series graphs of synchronized orbits.Figure 13
Some of the synchronized phase spaces projected to different planes.


## 10. Conclusion

In this study, a 5-dimensional memristor-based hyperchaotic circuit in the context of the fractional operator was considered and different qualitative and quantitative analyses are made: numerical approximation of the solution, bifurcation diagrams, Lyapunov exponents for different fractional orders, and parameter values are used to investigate the nature of the solution and chaotic system. The system is found to be hyperchaotic. Sensitivity to initial conditions and sensitivity to parameter value changes were found to cause a significant effect on the dynamics of the system. Furthermore, in the course of making several simulation results on the dynamics of the system, the hyperchaotic system is sensitive to almost any change made in the system including time scale increment change used in the study. The sensitivity of the system to fractional derivative order is very strong in the sense that the simulation results forq = 0.98 and 0.99 are significantly different (see Figures 2, 3, and 14). This is perhaps related to the contraction rate of volumes due to the different fractional-order derivatives. The magnitude of the sum of all the LEs corresponding to different fractional-order derivatives decreases as the order of the derivative increases for the case of the hyperchaotic system considered in this study (see Table 1). This requires further investigation. Synchronization of the system is also established and exhibited a strong agreement. The dependence of synchronization on the fractional derivative orders and variation of parameters needs further investigation.Figure 14
Some of the time series orbits of the system (25) due to different fractional-order derivatives.
(a)(b)(c)(d)



---

*Source: 6027246-2021-10-25.xml*