On a new four-dimensional model of memristor-based chaotic circuit in the context of nonsingular Atangana–Baleanu–Caputo operators

A memristor is naturally a nonlinear and at the same time memory element that may substitute resistors for next-generation nonlinear computational circuits that can show complex behaviors including chaos. A four-dimensional memristor system with the Atangana–Baleanu fractional nonsingular operator in the sense of Caputo is investigated. The Banach fixed point theorem for contraction principle is used to verify the existence–uniqueness of the fractional representation of the given system. A newly developed numerical scheme for fractional-order systems introduced by Toufik and Atangana is utilized to obtain the phase portraits of the suggested system for different fractional derivative orders and different parameter values of the system. Analysis on the local stability of the fractional model via the Matignon criteria showed that the trivial equilibrium point is unstable. The dynamics of the system are investigated using Lyapunov exponents for the characterization of the nature of the chaos and to verify the dissipativity of the system. It is shown that the supposed system is chaotic and it is significantly sensitive to parameter variation and small initial condition changes.


Introduction
In the last decade, fractional differential equations started gaining much attention in modeling several real-world problems in different areas including mathematical epidemiology, physics, engineering, and many others, in which fractional-order operators are either with the singular kernel (Caputo derivative and Riemann-Liouville fractional derivatives) or nonsingular kernels (Atangana-Baleanu and Caputo-Fabrizio derivatives) .
The difference between integer-and fractional-order derivatives is that the integer-order derivative indicates some properties at a particular time of a system, while the fractionalorder derivation operator describes a certain feature of a dynamical system for the whole time. Moreover, the integer-order derivative describes the local properties of a certain dynamic system, whereas the fractional-order derivative representation of a dynamic system involves the whole space of the process [33]. In other words, applying derivation operators via noninteger orders in modeling real-world problems is essential for describing the hereditary specifications and effectiveness of the memory as an important feature of different mechanisms in the problem [34,35].
The emergence of different definitions of a fractional derivative is interesting and is an opportunity to address the complexity of nature in the sense that some problems in nature follow the power law for the case of Riemann-Liouville fractional operator, others follow the Mittag-Leffler law for the case of Atangana-Baleanu fractional derivative operator, and others the exponential decay law l for the case of Caputo-Fabrizio fractional operator or a combination of the above [36,37].
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 in relation to complex systems are found to be nonlocal having long memory in time, and intrinsically fractional derivation operators can describe such systems more accurately than the integer derivatives [38]. In other words, important features of many physical systems are best described or exposed by using fractional-order operators.
Chaos theory has attracted many researchers and has applications in the fields of encryption and secure communication [39], modeling financial systems and representing circuit diagrams [40], and many others [37]. It is believed that the first chaotic phenomenon occurred when the 3-body problem was studied by Poincare in the 1980s. In his study Poincare realized that the 3-body problem is not integrable and the numerical solution used for the system is highly sensitive to initial conditions. Later on, in the middle of the twentieth century, Edward Lorenz modeled atmospheric problems using three coupled nonlinear differential equations, and he also noted sensitivity to initial conditions of his atmospheric model. The sensitivity to initial conditions made predicting the weather for a long enough amount of time very difficult or impossible [41].
Besides, Lorenz recognized that trajectories of a chaotic system are not scattered all over the places but approach an attractor of the state space. He also detected that sensitivity to initial conditions reveals itself by generating instability in the attractors and as a result, an attractor is named 'strange' attractor [42].
Nowadays there are many kinds 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 systems. Some of the literature is reviewed in what 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 different chaotic models investigated in [40,[43][44][45][46]. Different chaotic systems including Chua's electric circuit and several other chaotic systems were analyzed using fractionalbased order mathematical models by Petras [47]. Bifurcation and chaotic behaviors in fractional-order simplified Lorenz system using Adams-Bashforth-Moulton predictorcorrector scheme is considered in [48]. A dynamical fractional order of HIV-1 in the Caputo fractional order sense that led to chaotic behavior is considered in [49].
Atangana-Baleanu fractional derivative operators were used for modeling and analysis of different chaotic and hyperchaotic systems, and solutions were approximated by im-plementing a two-step Adams-Bashforth numerical algorithm in [50]. An intra-specific relation of predators and a prey model based on the Atangana-Baleanu fractional derivation operator, in which the numerical simulation led to a more chaotic dynamic system for different fractional derivative orders, were considered in [51].
It was in 1971 that Chua, circuit theorist, proposed memristor as a missing two-terminal nonlinear electrical component. The three basic components of a circuit are resistor, capacitor, and inductor. Memristor is famous for its memory effect and nonlinear specifications which is considered as the fourth component of a circuit. 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 [39,52].
At present, there have been several studies conducted on memristor-based chaotic circuits using both integer and fractional-order derivatives. In [53] a conformal model of the simplest fractional memristor-based chaotic circuit is designed and studied in direction of the conformable ADM (Adomian decomposition method), bifurcation diagram, Lyapunov exponent, and Poincare sections. Buscarino et al. introduced a chaotic circuit with the help of a realistic model of an HP memristor, and numerical results showed the generation of chaotic attractors [52]. A novel five-dimensional chaotic system via flux-controlled memristor and integer-order derivative, extracted from Wang's 4D hyperchaotic system, is proposed by Wang et al. [54]. A memristor-based chaotic circuit modified by putting a nonlinear resistor in the circuit due to Chua via a flux-controlled memristor and negative conductance is analyzed for its chaotic dynamics using Lyapunov exponent, bifurcation mapping, and Poincare mapping and power spectrum, along with laboratory experiments using integer-order derivatives in [39]. Then Petras [47] investigated a memristor-based Chua's oscillator in the framework of the fractional orders.
The objective of this work is to investigate the memory effect properties and detection of chaos in a four-dimensional memristor-based system. Accordingly, an integer model of memristor-based circuit is represented by Atangana-Baleanu fractional derivatives of the Caputo type (ABC) and the existence-uniqueness of the mentioned ABC fractional model are performed based on the Banach fixed point theorem (BFPT) for contraction principle. Numerical approximation of the ABC fractional model is made using the newly developed numerical approximation for fractional derivative by Toufik and Atanganain [55]. The local stability of given fractional model is accompanied using the Matignon stability criterion. The existence and nature of chaos in the fractional model are checked using Lyapunov exponents. A bifurcation diagram for different fractional derivatives and parameter variations is performed. Several phase portraits are depicted as a verification for the impact of different parameter values and different fractional 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. All the phase portraits and solution trajectories in this work are obtained from the numerical scheme of Toufik and Atangana adapted for the memristor chaotic model considered in this work. A computing software called Matlab 2019a is used for the simulation of different results.
A Matlab code in relation to Lyapunov exponents and bifurcation diagram of suggested fractional systems called Danca algorithm [56] 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 pieces of evidence for the orig-inality of this work include 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, and obtaining the phase portraits of the system from the numerical scheme. The Lyapunov exponents are calculated and bifurcation diagrams are depicted for different fractional orders and different range of parameter values of the model respectively. The impact of small changes in given initial conditions on the dynamics of the chaotic system is also investigated using simulation results for different initial values.
The remaining part of the research manuscript is arranged as follows: In the second part of the paper, the fractional representation of the memristor-based systems is made following some recapping of preliminary concepts and definitions of memristor circuit and basic properties of Atangana-Baleanu derivatives. The existence-uniqueness of the solution for the fractional representation of such a model are accomplished in this same section of the paper. The numerical algorithm applied to get the phase portrait of the memristor-based system is developed in the third part of the paper. The fourth part of the paper is concerned with the local stability analysis of the fractional model, and then, Lyapunov exponents, bifurcation diagrams with different fractional order and parameter variation are considered in the fifth section. Investigation of the impact of small change in the initial conditions on the dynamics of the chaotic system is considered in the sixth section, and the conclusions and references are presented in the subsequent sections.

Mathematical model describing memristor-based circuit 2.1 Memristor-based circuit
A memristor is regarded as a 2-terminal tool in which the magnetic flux (φ) between the terminals is considered as a function of the electric charge q that passes through the device [57]. Here, the memristor is of flux-controlled type specified by its incremental memductance function W (φ) given as follows: From (1) we reach a relation between the current I M (t) through the memristor and the voltage V 1 (t) which is presented as By (2) and W (φ(t)) = W ( V 1 (t)), the integral operator on the memductance function indicates that the mentioned function remembers the history of the voltage values. On the other side, if then a memristor will be a resistor. The memristor was derived from a circuit attributed to Chua by putting the Chua diode with the flux-controlled memristor [47,57]. Chua's electric circuit involves the resistor R, the inductor L, capacitors C j , j = 1, 2, and a nonlinear resistor (NR). The dynamic equation for the memristor-based chaotic circuit is formulated by in which V 1 , V 2 , and R L stand for the voltages over the capacitors C j , j = 1, 2, I L stands for the current through the existing inductance L, and I M (t) is defined in (2). Based on the motivation that smooth nonlinearity does give rise to chaos, let us choose a smooth continuous cubic monotone increasing nonlinearity presented by Consequently, the memductance function is formulated as The system of differential equations (3) can be converted to a structure without the dimension given as follows: where

Fractional representation of the memristor-based circuit
In this subsection, we recall the definitions and basic properties of Atangana-Baleanu derivative of Caputo type (ABC fractional derivative). where is the Mittag-Leffler function.
Now we follow to formulate model (6) in terms of ABC fractional derivatives.
. Then the following fractional differential equation admits a solution uniquely given by We can describe the result in Lemma 2.5 in the form of the Banach fixed point theorem as follows. We begin by defining a Banach space and and and We are now ready to describe the dynamic equation for the aforesaid memristor-based circuit given in (6) using the ABC fractional derivative in which

Results on the existence-uniqueness
Here, the existence-uniqueness of the solution for the ABC fractional model given in (16) are proved using Banach fixed point theorem (BFPT) for contraction mapping. The following two theorems are worth recalling before proceeding further.
Theorem 2.6 (BFPT, [61]) Let = ∅ be any closed subset of a Banach space X . Then any contraction mapping O : → admits a fixed point uniquely.

Theorem 2.7
Assume that x, y, z, φ are continuous mappings satisfying the following assumptions: Then the ABC fractional derivative system given by (16) has a unique solution in the region X .
Proof Let us show that the operator O 1 defined in (12) is well defined in the sense that and Now, for any x ∈ r and from (12), we have To show continuous differentiability on I, we proceed from which is continuous on I, and so we conclude that ABC To show that the operator O 1 has a fixed point, we apply Theorem 2.6. Based on this theorem, it is enough to show that O 1 is a contraction mapping. Indeed, let x 1 , x 2 ∈ X , t ∈ I. Hence, . for Now, for any y ∈ r and from (13), we have Consequently, we have O 2 y(t) ∈ r . To show continuous differentiability on I, we proceed from which is continuous on I, and so we conclude that ABC To show that the operator O 2 has a fixed point, we apply Theorem 2.6. Based on this theorem, it is enough to show that O 2 is a contraction mapping. Indeed, let y 1 , y 2 ∈ X , t ∈ I. Thus, Since H < 1, by the hypothesis of the theorem, we conclude that O 2 is a contraction.
In the sequel, we verify that the operator O 3 defined in (14) is well defined in the sense that O 3 z(t) ∈ r and ABC for Now, for any z ∈ r and from (14), we have To show continuous differentiability on I, we proceed from which is continuous on I, and so we conclude that ABC To show that the operator O 3 has a fixed point, we apply Theorem 2.6. Based on this theorem, it is enough to show that O 3 is a contraction mapping. Indeed, let z 1 , z 2 ∈ X , t ∈ I. Hence, . Since H < 1, by the hypothesis of the theorem, we conclude that O 3 is a contraction.
Now, for any φ ∈ r and from (15), we have Consequently, we have O 4 φ(t) ∈ r . To show continuous differentiability on I, we proceed from ABC 0 which is continuous on I, and so we conclude that ABC As a result, O 4 r ⊂ r . Immediately, it follows that O 4 is a contraction. Hence, by BFPT 2.6, system (16) admits a solution uniquely in X .

Numerical approximation of solutions of a given model
In this section, the numerical algorithm utilized to get the phase portrait of the dynamic ABC-system (16) is introduced. In the context of chaotic or hyperchaotic fractional differential equations, the use of analytical methods such as the Laplace transform method, the Sumudu transform, the HATM, and the homotopy perturbation technique cannot easily be applied because of the nonlinearities of the system [36,37]. 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 Toufic-Atangana numerical schemes [38]. Both of these methods are based on Lagrange interpolation polynomials.
In this study, the newly developed numerical approximation for fractional derivatives developed by Toufic and Atangana is employed. The numerical scheme is particularly developed for approximation of Atangana-Baleanu fractional-derivatives, and it is proved to be convergent, stable, and consistent.
For convenience, let us write (16) in the form shown in the following: where Now from Lemma 2.5 and the first equation of (21) we have The solution to FIVP (22) is given as

x(t), y(t), z(t), φ(t)
With the help of Lagrange's interpolation polynomial on [t k , t k+1 ], the equality leads to where h = t kt k-1 . Substituting (24) into (23), we obtain where and Inserting t m = mh into (26) and (27), we get Expression (25) can be expressed in terms of (28) and (29) as shown in Accordingly, we obtain other equations for the rest of the variables and and

Analysis on local stability
In the current part, the local stability analysis of the fractional model represented in (21) is performed. It is known in general that the equilibrium points (EPs) of chaotic systems are not 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 most commonly used in the literature for the same purpose [1,35,37,45]. The Matignon condition is presented as in which J represents the Jacobian matrix, λ(J) is the set of all existing eigenvalues of J, and q is the fractional derivative order. In relation to fractional derivatives, the EP of (21) is named locally stable provided that the Matignon criterion (35) is satisfied for each of the eigenvalues of J.
To determine if the equilibrium points (EPs) of (21) are stable or not, we proceed as follows. (

I) The equilibrium points (EPs):
The equilibrium points of (21) are obtained by solving ABC which is found to be E eqpts = (0, 0, 0, δ) for any arbitrary constant δ. Hence, it is a trivial equilibrium point for δ = 0, and the nontrivial equilibrium point is a line of equilibrium points given by E eqpts = {(0, 0, 0, δ) : δ = 0}.
On the other side, | arg(λ 2,3 )| = 2.0133 > q π 2 is true for any q ∈ (0, 1). It can then be concluded that the equilibrium point E eqpts = (0, 0, 0, 0) is locally unstable for the considered parameter values. Furthermore, since the real part of one of the eigenvalues of the Jacobian matrix is positive, so we infer that system (21) fulfills the necessary criterion for showing the double scroll attractor [1].

Lyapunov exponents, bifurcation, and chaos via different values of q and parameters
In this section, the level of chaos in system (21) is quantified using the Lyapunov exponent method. Bifurcation analysis of system (21) related to fractional derivative order q and three of the parameters in the model β 1 , β 2 , and β 3 are performed. In addition, a Matlab pseudo-code for Lyapunov exponents of fractional systems called the Danca algorithm [56] is used to quantify the chaos by calculating Lyapunov exponents for different fractional orders of model (21). The values for parameters are a 1 = 0.3, a 2 = 0.8, ς = 1.4, β 1 = 10, β 2 = 13, and β 3 = 0.1. These values are similar to the above ones used for calculating the Jacobian matrix. The initial conditions used in this part of the work are given by (0.11, 0.11, 0.11, 0.11). The final time of the simulation is t end = 300 seconds. The corresponding Lyapunov exponents (LE) for different values q = 0.95, 0.96, 0.98, 0.99 can be found in Table 1. As we can observe from Table 1, in each of the rows of Table 1, one of the Lyapunov exponents is positive, and as a result, the dynamic system in (21) is chaotic, and, since the sum of the LEs in each row is negative, we conclude that the system is dissipative.

Bifurcations due to variation of q
To obtain bifurcation diagrams due to the variation of the fractional order q, the values of all the parameters are kept fixed and the order of the fractional derivative q is made to vary in the interval (0.95, 1) with an increment of 0.01. The bifurcation diagram is illustrated in Fig. 1.
As shown in Fig. 1, the system started to bifurcate at q = 0.95, and then for q ∈ (0.95, 1] the system is chaotic. This observation is verified by the phase portrait of given model (21) depicted in Figs. 2 and 3 for q = 0.95 and q = 0.99, respectively. It is observed from  the figures that as the order of the derivative increases to 1, the chaotic nature of dynamic system (21) gets more and more significant. Based on the LE shown in Table 1 for the fractional order q = 0.95, the Kaplan-Yorke dimension of system (21) can be calculated as follows: Since the sum of the first three LEs are positive, i.e., LE1 + LE2 + LE3 > 0, and the sum of all the LEs is less than zero, thus this system is dissipative, and accordingly, there is the Kaplan-Yorke dimension which equals The Kaplan-Yorke dimension corresponding to q = 0.98 is given as Similarly, the other Kaplan-Yorke dimension corresponding to q = 0.99 is given as The overall conclusion from the LEs is that the system is a chaotic system with one positive LE. The Lyapunov dimension is found to be noninteger for all the fractional orders considered above. Another observation is that though the fractional order q increased from 0.95 to 0.99, the LEs are approximately equal to each other exactly to three decimal places.
In the next part of this section, the chaos and bifurcation with different values of parameters are investigated. The three parameters considered are β 1 , β 2 , and β 3 .
One can observe from Fig. 4 that for the parameter β 1 ∈ [9.5, 10.5] system (21) gets significantly chaotic. To verify the observation, the phase portrait for β 1 = 10.09 is depicted in Fig. 5.

Bifurcation diagram due to variation of β 2
To obtain the bifurcation diagram as a result of variation of parameter β 2 , the other parameter values are kept unchanged and the order is taken to be q = 0.98. The parameter   β 2 is made to vary in the interval (12.5, 13.5) with an increment of 0.001. The relevant diagram is illustrated in Fig. 6.
One can infer from Fig. 6 that when β 2 ∈ [12.5, 13.5], system (21) remains chaotic. To verify this observation, the phase portraits of system (21) for β 2 = 12.5 and β 2 = 13.02 are illustrated in Figs. 7 and 8, respectively.  Based on the bifurcation diagram (Fig. 6) and the figures for β 2 = 12.5 and β 2 = 13.02, it can be observed that the significance of the chaos gets weaker and weaker as the values of the parameter increase to 13.5 in the interval (12.5, 13.5).
This observation is also supported by the LEs and corresponding Kaplan-Yorke dimension of the attractors shown in Table 2. As is observed from Table 2, the positive LEs and the corresponding dimension decrease by increasing the value of β 2 , and because of this, the significance of the chaos decreases by increasing the value of the parameter in the interval (12.5, 13.5).

Bifurcation diagram due to variation of β 3
To obtain the bifurcation diagram as a result of variation of parameter β 3 , the values of other parameters are kept unchanged and the order is taken to be q = 0.98. The parameter β 3 is made to vary in the interval (0, 1) with an increment of 0.001. The bifurcation diagram is illustrated in Fig. 9.

Impact of initial conditions
The impact of applying different initial conditions on dynamic system (21) is addressed using different initial conditions and simulation results in this section. Since chaotic systems are famous to be sensitive to given initial conditions, it is better to investigate dynamic system (21) in terms of different initial conditions. The values of used parameters are q = 0.98, β 1 = 10, β 2 = 13, β 3 = 0.1, a 1 = 0.3, a 2 = 0.8, ς = 1.4. As well as, the used initial condition is X 0 = [x 0 , y 0 , z 0 , φ 0 ] = [0.11, 0.11, 0.11, 0.11].
The first case of investigation of the influence of initial condition on the dynamics of system (21), x 0 is made to vary from 0.11 to 0.001, 0.1, 0.2, and 0.12. The simulation result Figure 12 The solution trajectory φ(t) of system (21) for different initial conditions due to variation of x 0 Figure 13 The solution trajectory φ(t) of system (21) for different initial conditions due to variation of y 0 is shown in Fig. 12 for the solution φ(t) of system (21) corresponding to the different initial conditions obtained by varying the first coordinate of X 0 . The other solution graphs of system (21) have a similar pattern as Fig. 12 for this case of initial conditions. Secondly, keeping all the values of parameters fixed as in the first case, the initial value y 0 is made to vary from 0.11 to 0.001, 0.1, and 0.12. The simulation result is shown in Fig. 13 for the solution φ(t) of system (21) corresponding to the different initial conditions obtained by varying the second coordinate of X 0 . The other solution graphs of system (21) have a similar pattern as Fig. 13 for this case of the initial condition.
In general, we can conclude that variation of the initial conditions in the chaotic systems generates different properties of the solution trajectories of the system. This happens because the values of Lyapunov exponents are very sensitive to initial conditions.

Conclusions
In this study, a four-dimensional memristor-based system of a circuit along with a flux-controlled memristor is investigated. The integer-order system is represented by Atangana-Baleanu fractional derivative and the Banach contraction theorem is utilized to show the existence-uniqueness of solutions of the fractional representation of the mathematical model. The memristor fractional model exhibited chaotic behavior for different fractional-order derivatives and different values of the parameters. One can say that the fractional representation showed additional chaotic behaviors that may not be clear using integer derivatives. This is verified by simulation of the system for different fractionalorder derivatives. The Lyapunov exponent was applied to show the nature of the chaotic system, and it is found that there is one positive Lyapunov exponent and then the mentioned system is chaotic. The memristor-based chaotic system is found to be sensitive to small variations of the parameters as depicted by bifurcation diagrams. The dynamics of the system are also found to be sensitive to the initial conditions as it is depicted via different simulation results. For the future, other fractional operators such as Caputo-Fabrizio could be applied to the same system and compared to the result obtained in this work.