Theory and Modern Applications

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

## Abstract

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 fractional-order 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 . 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 . 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 , modeling financial systems and representing circuit diagrams , and many others . 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 .

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 .

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, 4346]. Different chaotic systems including Chua’s electric circuit and several other chaotic systems were analyzed using fractional-based order mathematical models by Petras . Bifurcation and chaotic behaviors in fractional-order simplified Lorenz system using Adams–Bashforth–Moulton predictor-corrector scheme is considered in . A dynamical fractional order of HIV-1 in the Caputo fractional order sense that led to chaotic behavior is considered in .

Atangana–Baleanu fractional derivative operators were used for modeling and analysis of different chaotic and hyperchaotic systems, and solutions were approximated by implementing a two-step Adams–Bashforth numerical algorithm in . 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 .

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 $$(\phi =\phi (q))$$, or it models a relationship between charges and flux $$(q=q(\phi ))$$ 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  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 . 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. . 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 . Then Petras  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 . 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  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 originality 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

### Memristor-based circuit

A memristor is regarded as a 2-terminal tool in which the magnetic flux $$(\phi )$$ between the terminals is considered as a function of the electric charge q that passes through the device . Here, the memristor is of flux-controlled type specified by its incremental memductance function $$W(\phi )$$ given as follows:

$$W(\phi ) \triangleq \frac{\mathrm{d}q(\phi )}{\mathrm{d}\phi }.$$
(1)

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

$$I_{M}(t) \triangleq \frac{\mathrm{d}q}{\mathrm{d}t} = \frac{\mathrm{d}q}{\mathrm{d}\phi } \frac{\mathrm{d}\phi }{\mathrm{d}t} = W\bigl(\phi (t)\bigr)V_{1}(t).$$
(2)

By (2) and $$W(\phi (t)) = W(\int 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

$$W\bigl(\phi (t)\bigr) = W\biggl( \int V_{1}(t)\biggr) = \mathrm{const},$$

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

\begin{aligned} &\frac{\mathrm{d}V_{1}(t)}{\mathrm{d}t}= \frac{1}{C_{1}} \biggl[ \frac{(V_{2}(t)-V_{1}(t))}{R} + GV_{1}(t) - I_{M}(t) \biggr], \\ &\frac{\mathrm{d}V_{2}(t)}{\mathrm{d}t}= \frac{1}{C_{2}} \biggl[ \frac{(V_{2}(t)-V_{1}(t))}{R} + I_{L}(t) \biggr], \\ &\frac{\mathrm{d}I_{L}(t)}{\mathrm{d}t}= \frac{1}{L} \bigl[ -V_{2}(t)-R_{L}I_{L}(t) \bigr], \\ &\frac{\mathrm{d}\phi }{\mathrm{d}t}= V_{1}(t), \end{aligned}
(3)

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

$$q(\phi ) = a_{0} + a_{1} \phi + a_{2} \phi ^{3},\quad a_{1} >0, a_{2} >0.$$
(4)

Consequently, the memductance function is formulated as

$$W(\phi ) = a_{1} + 3a_{2}\phi ^{2}.$$
(5)

The system of differential equations (3) can be converted to a structure without the dimension given as follows:

$$\textstyle\begin{cases} \frac{\mathrm{d}x}{\mathrm{d}t} = \beta _{1} (- W(\phi )x(t) + y(t)-x(t)+ \varsigma x(t) ), \\\frac{\mathrm{d}y}{\mathrm{d}t} = x(t)-y(t) + z(t), \\\frac{\mathrm{d}z}{\mathrm{d}t} = - \beta _{2} y(t)-\beta _{3}z(t), \\\frac{\mathrm{d}\phi }{\mathrm{d}t} = x(t), \end{cases}$$
(6)

where

$$\begin{gathered}x:=V_{1},\qquad y:=V_{2},\qquad z:=I_{L},\qquad C_{2}=1,\qquad R=1,\\ \beta _{1}:= \frac{1}{C_{1}},\qquad \varsigma :=G,\qquad \beta _{2} :=\frac{1}{L},\qquad \beta _{3} := \frac{R_{L}}{L}. \end{gathered}$$
(7)

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

### Definition 2.1

([35, 58, 59])

Let $$\phi \in C^{1}(a,b)$$ and $$\alpha \in (0,1)$$. The ABC $$\alpha ^{th}$$-derivative of ϕ is defined as

$${}^{ABC}_{a} \mathcal{D}_{t}^{\alpha } \phi (t) = \frac{F(\alpha )}{1-\alpha } \int _{a}^{t} \frac{\mathrm{d}\phi }{\mathrm{d}k} E_{\alpha } \biggl[ - \frac{\alpha }{1-\alpha } (t-k)^{\alpha } \biggr] \,\mathrm{d}k,$$
(8)

where $$F(\alpha ) = 1-\alpha - \frac{\alpha }{\Gamma (\alpha )}$$, and $$E_{\alpha }[\cdot ]$$ is the Mittag-Leffler function.

### Definition 2.2

([35, 58, 59])

The Atangana–Baleanu (AB) $$\alpha ^{th}$$-integral of $$\phi \in C^{1}(a,b)$$ is given by

$${}^{AB}_{a} \mathcal{I}_{t}^{\alpha } \phi (t) = \frac{1-\alpha }{F(\alpha )} \phi (t) + \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{a}^{t} \phi (k) (t-k)^{ \alpha - 1} \,\mathrm{d}k.$$
(9)

### Lemma 2.3

()

The ABC $$\alpha ^{th}$$-derivative and AB $$\alpha ^{th}$$-integral of $$\phi \in C^{1}(a,b)$$ fulfill the Newton–Leibniz relation

$${}^{AB}_{a} \mathcal{I}_{t}^{\alpha } \bigl( {}^{ABC}_{a} \mathcal{D}_{t}^{\alpha } \phi (t) \bigr) = \phi (t)-\phi (a).$$
(10)

### Lemma 2.4

()

For $$\phi _{1},\phi _{2}\in C^{1}(a,b)$$, the AB $$\alpha ^{th}$$-derivative in the Caputo sense satisfies the Lipschitz condition

$$\bigl\vert {}^{ABC}_{a} \mathcal{D}_{t}^{\alpha }\phi _{1}(t) - {}^{ABC}_{a} \mathcal{D}_{t}^{\alpha }\phi _{2}(t) \bigr\vert \leq \Lambda \bigl\Vert \phi _{1}(t)- \phi _{2}(t) \bigr\Vert .$$
(11)

Now we follow to formulate model (6) in terms of ABC fractional derivatives.

### Lemma 2.5

()

Let $$G\in C^{1}(a,b)$$. Then the following fractional differential equation

$${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }\phi (t) = G(t),\quad f(0)=G_{0}$$

admits a solution uniquely given by

$$\phi (t) = G_{0} + \frac{1-\alpha }{F(\alpha )} G(t) + \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} G(k)\,\mathrm{d}k.$$

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

$$\mathcal{X}=\bigl\{ x\in C(\mathbb{I},\mathbb{R}), \mathbb{I} = [0,1]\bigr\} ,$$

with a norm defined as $$\Vert x\Vert _{\mathcal{X}} =\sup_{t\in \mathbb{I}}\vert x(t)\vert$$.

Now define operators $$O_{1}, O_{2},O_{3}, O_{4}: \mathcal{X} \to \mathcal{X}$$ by

\begin{aligned} O_{1}x(t) &= x_{0} + \frac{1-\alpha }{F(\alpha )} \bigl( \beta _{1} \bigl(-W( \phi )x(t) + y(t)+ (\varsigma -1)x(t)\bigr) \bigr) \\ &\quad{}+ \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl( \beta _{1} \bigl(y(k)+ (\varsigma -1)x(k)-W(\phi )x(k)\bigr) \bigr) \,\mathrm{d}k, \end{aligned}
(12)

and

\begin{aligned} O_{2}y(t) &= y_{0} + \frac{1-\alpha }{F(\alpha )} \bigl( x(t)-y(t)+z(t) \bigr) \\ &\quad{}+ \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl( x(k)-y(k)+z(k) \bigr) \,\mathrm{d}k, \end{aligned}
(13)

and

\begin{aligned} O_{3}z(t) &= z_{0} + \frac{1-\alpha }{F(\alpha )} \bigl( -\beta _{2} y(t) - \beta _{3} z(t) \bigr) \\ &\quad{}+ \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl( -\beta _{2} y(k) - \beta _{3} z(k) \bigr) \,\mathrm{d}k, \end{aligned}
(14)

and

\begin{aligned} O_{4}\phi (t) = \phi _{0} + \frac{1-\alpha }{F(\alpha )} x(t) + \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} x(k) \,\mathrm{d}k. \end{aligned}
(15)

We are now ready to describe the dynamic equation for the aforesaid memristor-based circuit given in (6) using the ABC fractional derivative

$$\textstyle\begin{cases} {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }x(t) = \beta _{1} (- W(\phi )x(t) + y(t)-x(t)+ \varsigma x(t) ), \\{}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }y(t) = x(t)-y(t) + z(t), \\{}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }z(t) = - \beta _{2} y(t)-\beta _{3}z(t), \\{}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }\phi (t) = x(t), \end{cases}$$
(16)

in which $$x(0)=x_{0}$$, $$y(0)=y_{0}$$, $$z(0)=z_{0}$$, and $$\phi (0)=\phi _{0}$$ are initial conditions.

### 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, )

Let $$\Pi \neq \emptyset$$ be any closed subset of a Banach space $$\mathcal{X}$$. Then any contraction mapping $$O: \Pi \to \Pi$$ admits a fixed point uniquely.

### Theorem 2.7

Assume that x, y, z, ϕ are continuous mappings satisfying the following assumptions:

1. (D1)

$$\exists d_{1},d_{2},d_{3},d_{4}>0$$ s.t. $$\vert x(t)\vert < d_{1}$$, $$\vert y(t)\vert < d_{2}$$, $$\vert z(t)\vert < d_{3}$$, $$\vert \phi (t)\vert < d_{4}$$;

2. (D2)

$$\beta _{1}((\varsigma -1)+ \vert a_{1}\vert + 3a_{2}d_{4}^{2})((1- \alpha )\Gamma (\alpha )+ 1) < \Gamma (\alpha )F(\alpha )$$ and

$$(1-\alpha )\Gamma (\alpha ) + 1 < \Gamma (\alpha )F(\alpha ) \quad \textit{and} \quad \beta _{3} \bigl((1-\alpha )\Gamma (\alpha ) + 1 \bigr) < \Gamma ( \alpha )F(\alpha ).$$

Then the ABC fractional derivative system given by (16) has a unique solution in the region $$\mathcal{X}$$.

### Proof

Let us show that the operator $$O_{1}$$ defined in (12) is well defined in the sense that $$O_{1}x(t) \in \Im _{r}$$ and $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{1}x)(t) \in C(\mathbb{I}=[0,1], \mathbb{R})$$, where

$$\Im _{r} := \bigl\{ w\in \mathcal{X}, \Vert w \Vert _{\mathcal{X}} \leq r\bigr\} ,$$
(17)

and

$$r> \frac{ \vert x_{0} \vert F(\alpha )\Gamma (\alpha ) + ((1-\alpha )\Gamma (\alpha ) +1) (\beta _{1}(d_{2}+d_{1}(\varsigma -1)+ \vert a_{1}+3a_{2}d_{4}^{2} \vert d_{1}) )}{F(\alpha )\Gamma (\alpha )} .$$

Now, for any $$x\in \Im _{r}$$ and from (12), we have

\begin{aligned} & \bigl\Vert O_{1}x(t) \bigr\Vert _{\mathcal{X}} \\ &\quad \leq \vert x_{0} \vert + \frac{1-\alpha }{F(\alpha )} \bigl( \beta _{1} \bigl(d_{2}+ (\varsigma -1)d_{1} + \bigl\vert a_{1}+3a_{2}d_{4}^{2} \bigr\vert d_{1}\bigr) \bigr) \\ &\qquad{}+ \sup_{t\in \mathbb{I}} \biggl( \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl\vert \beta _{1} \bigl(y(k)+ (\varsigma -1)x(k)-W(\phi )x(k)\bigr) \bigr\vert \,\mathrm{d}k \biggr) \\ &\quad \leq \vert x_{0} \vert + \frac{1-\alpha }{F(\alpha )} \bigl( \beta _{1} \bigl(d_{2}+ (\varsigma -1)d_{1} + \bigl\vert a_{1}+3a_{2}d_{4}^{2} \bigr\vert d_{1}\bigr) \bigr) \\ &\qquad{}+ \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \bigl( \beta _{1} \bigl(d_{2}+ ( \varsigma -1)d_{1} + \bigl\vert a_{1}+3a_{2}d_{4}^{2} \bigr\vert d_{1}\bigr) \bigr) \sup_{t\in \mathbb{I}=[0,1]} \biggl( \int _{0}^{t} (t-k)^{\alpha -1} \,\mathrm{d}k \biggr) \\ &\quad \leq \frac{ \vert x_{0} \vert F(\alpha )\Gamma (\alpha ) + ((1-\alpha )\Gamma (\alpha ) + 1) (\beta _{1}(d_{2}+d_{1}(\varsigma -1)+ \vert a_{1}+3a_{2}d_{4}^{2} \vert d_{1}) )}{F(\alpha )\Gamma (\alpha )} < r. \end{aligned}

Consequently, we have $$O_{1}x(t) \in \Im _{r}$$. To show continuous differentiability on $$\mathbb{I}$$, we proceed from

\begin{aligned} {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{1}x) (t) &= {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }{}^{AB}_{0} \mathcal{I}_{t}^{\alpha } \bigl( \beta _{1} \bigl(- W( \phi )x(t) + y(t)+(\varsigma -1) x(t) \bigr) \bigr) \\ &= \beta _{1} \bigl(- W(\phi )x(t) + y(t) + (\varsigma -1) x(t) \bigr), \end{aligned}

which is continuous on $$\mathbb{I}$$, and so we conclude that $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{1}x)(t)$$ is continuous on $$\mathbb{I}$$. As a result $$O_{1}\Im _{r} \subset \Im _{r}$$.

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} \in \mathcal{X}$$, $$t\in \mathbb{I}$$. Hence,

\begin{aligned} & \Vert O_{1}x_{1} -O_{1}x_{2} \Vert _{\mathcal{X}} \\ &\quad \leq \frac{1-\alpha }{F(\alpha )} \beta _{1} \bigl( (\varsigma -1) \Vert x_{1} -x_{2} \Vert _{\mathcal{X}} + \bigl\vert W(\phi ) \bigr\vert \Vert x_{1}-x_{2} \Vert _{\mathcal{X}} \bigr) \\ &\qquad{}+ \frac{\alpha }{F(\alpha ) \Gamma (\alpha )} \beta _{1} \bigl( ( \varsigma -1) \Vert x_{1}-x_{2} \Vert _{\mathcal{X}} + \bigl\vert W(\phi ) \bigr\vert \Vert x_{1}-x_{2} \Vert _{\mathcal{X}} \bigr) \sup_{t\in \mathbb{I}} \int _{0}^{t} (t-k)^{\alpha -1} \,\mathrm{d}k \\ &\quad \leq \biggl( \frac{\beta _{1}((\varsigma -1)+ \vert a_{1} \vert + 3a_{2}d_{4}^{2})(\Gamma (\alpha )(1-\alpha )+ 1)}{F(\alpha )\Gamma (\alpha )} \biggr) \Vert x_{1}-x_{2} \Vert _{\mathcal{X}} \\ &\quad = H \Vert x_{1}-x_{2} \Vert _{\mathcal{X}}, \end{aligned}

where $$H= \frac{\beta _{1}((\varsigma -1)+ \vert a_{1}\vert + 3a_{2}d_{4}^{2})(\Gamma (\alpha )(1-\alpha )+ 1)}{F(\alpha )\Gamma (\alpha )}$$. Since $$H<1$$, by the hypothesis of the theorem, we conclude that $$O_{1}$$ is a contraction.

To show that the remaining operators $$O_{2}$$, $$O_{3}$$, $$O_{4}$$ are contractions, we proceed in the same manner.

Let us verify that the operator $$O_{2}$$ defined in (13) is well defined in the sense that $$O_{2}y(t) \in \Im '_{r'}$$ and $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{2}y)(t) \in C(\mathbb{I}=[0,1], \mathbb{R})$$, where

$$\Im '_{r'} := \bigl\{ w\in \mathcal{X}, \Vert w \Vert _{\mathcal{X}} \leq r' \bigr\} ,$$
(18)

for

$$r'> \frac{ \vert y_{0} \vert F(\alpha )\Gamma (\alpha ) + ((1-\alpha )\Gamma (\alpha ) + 1) ( d_{1} + d_{2} + d_{3} )}{F(\alpha )\Gamma (\alpha )} .$$

Now, for any $$y\in \Im '_{r'}$$ and from (13), we have

\begin{aligned} \bigl\Vert O_{2}y(t) \bigr\Vert _{\mathcal{X}} &\leq \vert y_{0} \vert + \frac{1-\alpha }{F(\alpha )} ( d_{1} + d_{2} + d_{3} ) \\ &\quad{}+ \sup_{t\in \mathbb{I}} \biggl( \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl\vert x(k)-y(k)+z(k) \bigr\vert \,\mathrm{d}k \biggr) \\ &\leq \vert y_{0} \vert + \frac{1-\alpha }{F(\alpha )} ( d_{1} + d_{2} + d_{3} ) \\ &\quad{}+ \frac{\alpha }{F(\alpha )\Gamma (\alpha )} ( d_{1} + d_{2} + d_{3} ) \sup_{t\in \mathbb{I}=[0,1]} \biggl( \int _{0}^{t} (t-k)^{ \alpha -1} \,\mathrm{d}k \biggr) \\ &\leq \frac{ \vert y_{0} \vert F(\alpha )\Gamma (\alpha ) + ((1-\alpha )\Gamma (\alpha ) + 1) ( d_{1} + d_{2} + d_{3} )}{F(\alpha )\Gamma (\alpha )} < r'. \end{aligned}

Consequently, we have $$O_{2}y(t) \in \Im '_{r'}$$. To show continuous differentiability on $$\mathbb{I}$$, we proceed from

$${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{2}y) (t) = {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }{}^{AB}_{0} \mathcal{I}_{t}^{\alpha } \bigl( x(t)-y(t)+z(t) \bigr) = x(t)-y(t)+z(t),$$

which is continuous on $$\mathbb{I}$$, and so we conclude that $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{2}y)(t)$$ is continuous on $$\mathbb{I}$$. As a result, $$O_{2}\Im '_{r'} \subset \Im '_{r'}$$.

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} \in \mathcal{X}$$, $$t\in \mathbb{I}$$. Thus,

\begin{aligned} & \Vert O_{2}y_{1} -O_{2}y_{2} \Vert _{\mathcal{X}} \\ &\quad \leq \frac{1-\alpha }{F(\alpha )} \bigl( \Vert y_{1} -y_{2} \Vert _{\mathcal{X}} \bigr) + \frac{\alpha }{F(\alpha ) \Gamma (\alpha )} \bigl( \Vert y_{1} -y_{2} \Vert _{\mathcal{X}} \bigr) \sup _{t\in \mathbb{I}} \int _{0}^{t} (t-k)^{\alpha -1} \,\mathrm{d}k \\ &\quad \leq \biggl( \frac{(1-\alpha )\Gamma (\alpha ) + 1}{F(\alpha )\Gamma (\alpha )} \biggr) \Vert y_{1}-y_{2} \Vert _{\mathcal{X}} \\ &\quad = H' \Vert y_{1}-y_{2} \Vert _{\mathcal{X}}, \end{aligned}

where $$H'= \frac{(1-\alpha )\Gamma (\alpha ) + 1}{F(\alpha )\Gamma (\alpha )}$$. 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) \in \Im ''_{r''}$$ and $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{3}z)(t) \in C(\mathbb{I}=[0,1], \mathbb{R})$$, where

$$\Im ''_{r''} := \bigl\{ w\in \mathcal{X}, \Vert w \Vert _{\mathcal{X}} \leq r'' \bigr\}$$
(19)

for

$$r''> \frac{ \vert z_{0} \vert \Gamma (\alpha )F(\alpha ) + (\Gamma (\alpha )(1-\alpha ) + 1) ( \beta _{2} d_{2} + \beta _{3} d_{3} )}{\Gamma (\alpha )F(\alpha )} .$$

Now, for any $$z\in \Im ''_{r''}$$ and from (14), we have

\begin{aligned} \bigl\Vert O_{3}z(t) \bigr\Vert _{\mathcal{X}} &\leq \vert z_{0} \vert + \frac{1-\alpha }{F(\alpha )} ( \beta _{2} d_{2} + \beta _{3} d_{3} ) \\ &\quad{}+ \sup_{t\in \mathbb{I}} \biggl( \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl\vert - \beta _{2} y(k)- \beta _{3} z(k) \bigr\vert \,\mathrm{d}k \biggr) \\ &\leq \vert z_{0} \vert + \frac{1-\alpha }{F(\alpha )} ( \beta _{2} d_{2} + \beta _{3} d_{3} ) \\ &\quad{}+ \frac{\alpha }{\Gamma (\alpha )F(\alpha )} ( \beta _{2} d_{2} + \beta _{3} d_{3} ) \sup_{t\in \mathbb{I}=[0,1]} \biggl( \int _{0}^{t} (t-k)^{\alpha -1} \,\mathrm{d}k \biggr) \\ &\leq \frac{ \vert z_{0} \vert F(\alpha )\Gamma (\alpha ) + ((1-\alpha )\Gamma (\alpha ) + 1) ( \beta _{2} d_{2} + \beta _{3} d_{3} )}{F(\alpha )\Gamma (\alpha )} < r''. \end{aligned}

Consequently, we have $$O_{3}z(t) \in \Im ''_{r''}$$. To show continuous differentiability on $$\mathbb{I}$$, we proceed from

$${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{3}z) (t) = {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }{}^{AB}_{0} \mathcal{I}_{t}^{\alpha } \bigl( - \beta _{2} y(t)- \beta _{3} z(t) \bigr) = -\beta _{2} y(t)- \beta _{3} z(t),$$

which is continuous on $$\mathbb{I}$$, and so we conclude that $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{3}z)(t)$$ is continuous on $$\mathbb{I}$$. As a result, $$O_{3}\Im ''_{r''} \subset \Im ''_{r''}$$.

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} \in \mathcal{X}$$, $$t\in \mathbb{I}$$. Hence,

\begin{aligned} & \Vert O_{3}z_{1} -O_{3}z_{2} \Vert _{\mathcal{X}} \\ &\quad \leq \frac{1-\alpha }{F(\alpha )} \bigl( \beta _{3} \Vert z_{1} -z_{2} \Vert _{\mathcal{X}} \bigr) + \frac{\alpha }{F(\alpha ) \Gamma (\alpha )} \bigl( \beta _{3} \Vert z_{1} -z_{2} \Vert _{\mathcal{X}} \bigr) \sup_{t \in \mathbb{I}} \int _{0}^{t} (t-k)^{\alpha -1} \,\mathrm{d}k \\ &\quad \leq \biggl( \frac{\beta _{3} ((1-\alpha )\Gamma (\alpha ) + 1 ) }{F(\alpha )\Gamma (\alpha )} \biggr) \Vert z_{1}-z_{2} \Vert _{\mathcal{X}} \\ &\quad = H'' \Vert z_{1}-z_{2} \Vert _{\mathcal{X}}, \end{aligned}

where $$H''= \frac{\beta _{3} ((1-\alpha )\Gamma (\alpha ) + 1 ) }{F(\alpha )\Gamma (\alpha )}$$. Since $$H''<1$$, by the hypothesis of the theorem, we conclude that $$O_{3}$$ is a contraction.

Lastly, we verify that the operator $$O_{4}$$ defined in (15) is well defined in the sense that $$O_{4}\phi (t) \in \Im '''_{r'''}$$ and $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{4}\phi )(t) \in C( \mathbb{I}=[0,1], \mathbb{R})$$, where

$$\Im '''_{r'''} := \bigl\{ w\in \mathcal{X}, \Vert w \Vert _{\mathcal{X}} \leq r'''\bigr\}$$
(20)

for

$$r'''> \frac{ \vert \phi _{0} \vert \Gamma (\alpha )F(\alpha ) + (\Gamma (\alpha )(1-\alpha ) + 1) d_{1}}{\Gamma (\alpha )F(\alpha )} .$$

Now, for any $$\phi \in \Im '''_{r'''}$$ and from (15), we have

\begin{aligned} \bigl\Vert O_{4}\phi (t) \bigr\Vert _{\mathcal{X}} &\leq \vert \phi _{0} \vert + \frac{1-\alpha }{F(\alpha )} d_{1} + \sup _{t\in \mathbb{I}} \biggl( \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \int _{0}^{t} (t-k)^{ \alpha -1} \bigl\vert x(k) \bigr\vert \,\mathrm{d}k \biggr) \\ &\leq \vert \phi _{0} \vert + \frac{1-\alpha }{F(\alpha )} d_{1} + \frac{\alpha }{F(\alpha )\Gamma (\alpha )} d_{1} \sup_{t\in \mathbb{I}=[0,1]} \biggl( \int _{0}^{t} (t-k)^{\alpha -1} \,\mathrm{d}k \biggr) \\ &\leq \frac{ \vert \phi _{0} \vert F(\alpha )\Gamma (\alpha ) + ((1-\alpha )\Gamma (\alpha ) + 1) d_{1}}{F(\alpha )\Gamma (\alpha )} < r'''. \end{aligned}

Consequently, we have $$O_{4}\phi (t) \in \Im '''_{r'''}$$. To show continuous differentiability on $$\mathbb{I}$$, we proceed from

$${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{4} \phi ) (t) = {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }{}^{AB}_{0} \mathcal{I}_{t}^{\alpha }x(t) = x(t),$$

which is continuous on $$\mathbb{I}$$, and so we conclude that $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }(O_{4}\phi )(t)$$ is continuous on $$\mathbb{I}$$. As a result, $$O_{4}\Im '''_{r'''} \subset \Im '''_{r'''}$$. Immediately, it follows that $$O_{4}$$ is a contraction. Hence, by BFPT 2.6, system (16) admits a solution uniquely in $$\mathcal{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 . 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:

$$\textstyle\begin{cases} {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }x(t) = M_{1} ( t,x(t),y(t),z(t), \phi (t) ), \\{}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }y(t) = M_{2} ( t,x(t),y(t),z(t), \phi (t) ), \\{}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }z(t) = M_{3} ( t,x(t),y(t),z(t), \phi (t) ), \\{}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }\phi (t) = M_{4} ( t,x(t),y(t),z(t), \phi (t) ), \end{cases}$$
(21)

where

$$\textstyle\begin{cases} M_{1} ( t,x(t),y(t),z(t),\phi (t) ) = \beta _{1} ( - W(\phi )x(t) + y(t) + (\varsigma -1) x(t)), \\M_{2} ( t,x(t),y(t),z(t),\phi (t) ) = x(t)-y(t) + z(t), \\M_{3} ( t,x(t),y(t),z(t),\phi (t) ) = - \beta _{2} y(t)-\beta _{3}z(t), \\M_{4} ( t,x(t),y(t),z(t),\phi (t) ) = x(t). \end{cases}$$

Now from Lemma 2.5 and the first equation of (21) we have

$$\textstyle\begin{cases} {}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }x(t) = M_{1} ( t,x(t),y(t),z(t), \phi (t) ), \\x(0)=x_{0}. \end{cases}$$
(22)

The solution to FIVP (22) is given as

\begin{aligned} x(t) &= x_{0} + \frac{1-\alpha }{F(\alpha )} M_{1}{ \bigl( t,x(t),y(t),z(t), \phi (t) \bigr)} \\ &\quad{}+ \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \int _{0}^{t} M_{1}{ \bigl(k,x(k),y(k),z(k), \phi (k) \bigr)} (t-k)^{\alpha -1} \,\mathrm{d}k. \end{aligned}
(23)

With the help of Lagrange’s interpolation polynomial on $$[t_{k},t_{k+1}]$$, the equality

$$M_{1} \bigl( s,x(s),y(s),z(s),\phi (s) \bigr) = \beta _{1} \bigl(y(s) + ( \varsigma -1) x(s) - W(\phi )x(s)\bigr),$$

\begin{aligned} x_{k} &\approx \frac{1}{h} \bigl[ (s-t_{k-1})M_{1}{ \bigl( t_{k},x(t_{k}),y(t_{k}),z(t_{k}), \phi (t_{k}) \bigr)} \\ &\quad {}- (s-t_{k})M_{1} { \bigl( t_{k-1},x(t_{k-1}),y(t_{k-1}),z(t_{k-1}), \phi (t_{k-1}) \bigr)} \bigr], \end{aligned}
(24)

where $$h=t_{k}-t_{k-1}$$. Substituting (24) into (23), we obtain

\begin{aligned} x(t_{n+1}) &= x_{0} + \frac{1-\alpha }{F(\alpha )} M_{1}{ \bigl( t_{k},x(t_{k}),y(t_{k}),z(t_{k}), \phi (t_{k}) \bigr)} \\ &\quad{}+ \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \sum_{m=1}^{n} \biggl( \frac{M_{1}{ ( t_{m},x(t_{m}),y(t_{m}),z(t_{m}),\phi (t_{m}) )}}{h} \Upsilon _{m-1} \\ &\quad{}- \frac{M_{1}{ ( t_{m-1},x(t_{m-1}),y(t_{m-1}),z(t_{m-1}),\phi (t_{m-1}) )} }{h} \Upsilon _{m} \biggr), \end{aligned}
(25)

where

\begin{aligned} \Upsilon _{m-1}& = \int _{t_{m}}^{t_{m+1}} (k-t_{m-1}) (t_{n+1}-k)^{ \alpha -1} \,\mathrm{d}k \\ &= - \frac{1}{\alpha } \bigl[ (t_{m+1}-t_{m-1}) (t_{n+1} -t_{m+1})^{\alpha }- (t_{m}-t_{m-1}) (t_{n+1}-t_{m})^{\alpha } \bigr] \\ &\quad{}- \frac{1}{\alpha (\alpha +1)} \bigl[ (t_{n+1} -t_{m+1})^{\alpha +1} (t_{n+1} -t_{m+1})^{\alpha }- (t_{n+1} -t_{m})^{\alpha + 1} \bigr], \end{aligned}
(26)

and

\begin{aligned} \Upsilon _{m} &= \int _{t_{m}}^{t_{m+1}} (k-t_{m-1}) (t_{n+1}-k)^{ \alpha -1} \,\mathrm{d}k = - \frac{1}{\alpha } \bigl[ (t_{m+1}-t_{m-1}) (t_{n+1} -t_{m+1})^{\alpha } \bigr] \\ &\quad{}- \frac{1}{\alpha (\alpha +1)} \bigl[ (t_{n+1} -t_{m+1})^{\alpha +1} - (t_{n+1} -t_{m})^{\alpha + 1} \bigr]. \end{aligned}
(27)

Inserting $$t_{m} = m h$$ into (26) and (27), we get

\begin{aligned} &\Upsilon _{m-1}= \frac{h^{\alpha +1}}{\alpha (\alpha +1)} \bigl[ (n+1-m)^{\alpha }(n-m+2+ \alpha ) - (n-m)^{\alpha }(n-m+2+2\alpha ) \bigr], \end{aligned}
(28)
\begin{aligned} &\Upsilon _{m}= \frac{h^{\alpha +1}}{\alpha (\alpha +1)} \bigl[ (n+1-m)^{ \alpha +1} - (n-m)^{\alpha }(n-m+1+\alpha ) \bigr]. \end{aligned}
(29)

Expression (25) can be expressed in terms of (28) and (29) as shown in

\begin{aligned} x(t_{n+1}) &= x_{0} + \frac{1-\alpha }{F(\alpha )} M_{1}{ \bigl( t_{n},x(t_{n}),y(t_{n}),z(t_{n}), \phi (t_{n}) \bigr)} + \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \\ &\quad{}\times \sum_{m=1}^{n} \biggl( \frac{M_{1}{ ( t_{m},x(t_{m}),y(t_{m}),z(t_{m}),\phi (t_{m}) )} }{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha }(n-m+2+ \alpha ) - (n-m)^{\alpha }(n-m+2+2 \alpha ) \bigr] \\ &\quad{}- \frac{M_{1}{ ( t_{m-1},x(t_{m-1}),y(t_{m-1}),z(t_{m-1}),\phi (t_{m-1}) )} }{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha +1} - (n-m)^{\alpha }(n-m+1+ \alpha ) \bigr] \biggr). \end{aligned}
(30)

Accordingly, we obtain other equations for the rest of the variables

\begin{aligned} y(t_{n+1}) &= y_{0} + \frac{1-\alpha }{F(\alpha )} M_{2}{ \bigl( t_{n},x(t_{n}),y(t_{n}),z(t_{n}), \phi (t_{n}) \bigr)} + \frac{\alpha }{\Gamma (\alpha )F(\alpha )} \\ &\quad{}\times \sum_{m=1}^{n} \biggl( \frac{M_{2}{ ( t_{m},x(t_{m}),y(t_{m}),z(t_{m}),\phi (t_{m}) )}}{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha }(n-m+2+ \alpha ) - (n-m)^{\alpha }(n-m+2+2 \alpha ) \bigr] \\ &\quad{}- \frac{M_{2}{ ( t_{m-1},x(t_{m-1}),y(t_{m-1}),z(t_{m-1}),\phi (t_{m-1}) )}}{\Gamma (\alpha +2)} \end{aligned}
(31)
\begin{aligned} &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha +1} - (n-m)^{\alpha }(n-m+1+ \alpha ) \bigr] \biggr), \end{aligned}
(32)

and

\begin{aligned} z(t_{n+1}) &= z_{0} + \frac{1-\alpha }{F(\alpha )} M_{3}{ \bigl( t_{n},x(t_{n}),y(t_{n}),z(t_{n}), \phi (t_{n}) \bigr)} + \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \\ &\quad{}\times \sum_{m=1}^{n} \biggl( \frac{M_{3}{ ( t_{m},x(t_{m}),y(t_{m}),z(t_{m}),\phi (t_{m}) )}}{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha }(n-m+2+ \alpha ) - (n-m)^{\alpha }(n-m+2+2 \alpha ) \bigr] \\ &\quad{}- \frac{M_{3}{ ( t_{m-1},x(t_{m-1}),y(t_{m-1}),z(t_{m-1}),\phi (t_{m-1}) )}}{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha +1} - (n-m)^{\alpha }(n-m+1+ \alpha ) \bigr] \biggr), \end{aligned}
(33)

and

\begin{aligned} \phi (t_{n+1}) &= \phi _{0} + \frac{1-\alpha }{F(\alpha )} M_{4}{ \bigl( t_{n},x(t_{n}),y(t_{n}),z(t_{n}), \phi (t_{n}) \bigr)} + \frac{\alpha }{F(\alpha )\Gamma (\alpha )} \\ &\quad{}\times \sum_{m=1}^{n} \biggl( \frac{M_{4}{ ( t_{m},x(t_{m}),y(t_{m}),z(t_{m}),\phi (t_{m}) )}}{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha }(n-m+2+ \alpha ) - (n-m)^{\alpha }(n-m+2+2 \alpha ) \bigr] \\ &\quad{}- \frac{M_{4}{ ( t_{m-1},x(t_{m-1}),y(t_{m-1}),z(t_{m-1}),\phi (t_{m-1}) )}}{\Gamma (\alpha +2)} \\ &\quad{}\times h^{\alpha } \bigl[ (n+1-m)^{\alpha +1} - (n-m)^{\alpha }(n-m+1+ \alpha ) \bigr] \biggr). \end{aligned}
(34)

## 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

$$\bigl\vert \arg \lambda (J) \bigr\vert > \frac{1}{2}q \pi ,$$
(35)

in which J represents the Jacobian matrix, $$\lambda (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}_{0} \mathcal{D}_{t}^{\alpha }x(t) =0$$, $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }y(t) = 0$$, $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }z(t) = 0$$, and $${}^{ABC}_{0} \mathcal{D}_{t}^{\alpha }\phi (t) = 0$$, which is found to be $$E_{eqpts} = (0,0,0,\delta )$$ for any arbitrary constant δ. Hence, it is a trivial equilibrium point for $$\delta =0$$, and the nontrivial equilibrium point is a line of equilibrium points given by $$E_{eqpts} = \{ (0,0,0,\delta ): \delta \neq 0 \}$$.

(II) The Jacobian matrix:

In this part of the paper, analysis on the local stability is made for the trivial equilibrium point for simplicity. The parameter values used unless otherwise mentioned in this work are taken to be $$a_{1} =0.3$$, $$a_{2}=0.8$$, $$C_{2}=1$$, $$\varsigma = 1.4$$, $$R=1$$, $$\beta _{1} = \frac{1}{C_{1}} = 10$$, $$\beta _{2} = \frac{1}{L} = 13$$, $$\beta _{3} = \frac{R_{L}}{L}=0.1$$. Thus, the Jacobian matrix of (21) computed at the trivial EP equals

$$J_{0} = \begin{pmatrix} 1 & 10 & 0 & 0 \\ 1 & -1 & 1 & 0 \\ 0 & -13 & -0.1 & 0 \\ 1 & 0 & 0 & 0 \end{pmatrix} .$$

(III) Eigenvalues:

The eigenvalues of $$J_{0}$$ are computed as $$\lambda _{1} = 0$$, $$\lambda _{2,3} = - 1.1063 \pm 2.3346i$$ and $$\lambda _{4} = 2.1126 + 0.0000i$$. To check the local stability of the trivial EP, using the Matignon criteria, we need to verify if all the eigenvalues satisfy condition (35). Indeed, $$\vert \arg (\lambda _{1}) \vert = \vert \arg (\lambda _{4})\vert = 0 > q \frac{\pi }{2}$$, which is not possible since $$q\in (0,1)$$.

On the other side, $$\vert \arg (\lambda _{2,3}) \vert = 2.0133 > q \frac{\pi }{2}$$ is true for any $$q\in (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 .

## 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 $$\beta _{1}$$, $$\beta _{2}$$, and $$\beta _{3}$$ are performed. In addition, a Matlab pseudo-code for Lyapunov exponents of fractional systems called the Danca algorithm  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$$, $$\varsigma = 1.4$$, $$\beta _{1}=10$$, $$\beta _{2} = 13$$, and $$\beta _{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 \in (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

$$\dim (LE) = 3 + \frac{LE1 + LE2 + LE3}{ \vert LE4 \vert } = 3.9098.$$

The Kaplan–Yorke dimension corresponding to $$q=0.98$$ is given as

$$\dim (LE) = 3 + \frac{LE1 + LE2 + LE3}{ \vert LE4 \vert } = 3.9092.$$

Similarly, the other Kaplan–Yorke dimension corresponding to $$q=0.99$$ is given as

$$\dim (LE) = 3 + \frac{LE1 + LE2 + LE3}{ \vert LE4 \vert } = 3.9090.$$

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 $$\beta _{1}$$, $$\beta _{2}$$, and $$\beta _{3}$$.

### Bifurcation diagram due to variation of $$\beta _{1}$$

The values of parameters are $$a_{1} = 0.3$$, $$a_{2} = 0.8$$, $$\varsigma = 1.4$$, $$\beta _{1} = 10$$, $$\beta _{2} = 13$$, and $$\beta _{3} = 0.1$$, and the initial condition is given by $$X_{0} = (0.11, 0.11, 0.11, 0.11)$$, and $$q= 0.98$$. Let the value of the parameter $$\beta _{1}$$ vary in $$(9.5,10.5)$$ with an increment of 0.001. The bifurcation diagram for this case is illustrated in Fig. 4.

One can observe from Fig. 4 that for the parameter $$\beta _{1} \in [9.5,10.5]$$ system (21) gets significantly chaotic. To verify the observation, the phase portrait for $$\beta _{1} = 10.09$$ is depicted in Fig. 5.

### Bifurcation diagram due to variation of $$\beta _{2}$$

To obtain the bifurcation diagram as a result of variation of parameter $$\beta _{2}$$, the other parameter values are kept unchanged and the order is taken to be $$q=0.98$$. The parameter $$\beta _{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 $$\beta _{2} \in [12.5,13.5]$$, system (21) remains chaotic. To verify this observation, the phase portraits of system (21) for $$\beta _{2} = 12.5$$ and $$\beta _{2} = 13.02$$ are illustrated in Figs. 7 and 8, respectively.

Based on the bifurcation diagram (Fig. 6) and the figures for $$\beta _{2}=12.5$$ and $$\beta _{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 $$\beta _{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 $$\beta _{3}$$

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

In the figure illustrated by varying the parameter $$\beta _{3}$$ in $$(0, 1)$$, we find a chaotic behavior of system (21) for $$\beta _{3} \in (0,0.182)$$. When $$\beta _{3}\in (0.182,0.51)$$, all the trajectories $$\phi (t)$$ tend to a nontrivial EP. To verify these findings, the phase portraits of model (21) corresponding to $$\beta _{3} = 0.11$$ and $$\beta _{3} = 0.4$$ are depicted in Figs. 10 and 11.

Note that each of the phase portraits projected to different planes indicated in Fig. 11 converges accordingly to the nontrivial equilibrium point (EP); that is, the phase portrait on the y-ϕ, z-ϕ, and x-ϕ planes converges to a nontrivial EP as illustrated in Fig. 11. The eigenvalues of the trivial EP of model (21) for $$\beta _{3}=0.4$$, ($$a_{1}=0.3$$, $$a_{2}=0.8$$, $$\varsigma = 1.4$$, $$\beta _{1} = 10$$, $$\beta _{2} = 13$$) are given as $$\lambda _{1} =0$$, $$\lambda _{2} = - 1.3113 \pm 2.7417i$$, and $$\lambda _{3} = 2.2226$$, and hence the trivial EP of system (21) is unstable for $$\beta _{3} = 0.4$$. Therefore the chance for system (21) in this particular case is to tend to a nontrivial EP as illustrated in Fig. 11.

## 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$$, $$\beta _{1} = 10$$, $$\beta _{2} = 13$$, $$\beta _{3} = 0.1$$, $$a_{1} = 0.3$$, $$a_{2} = 0.8$$, $$\varsigma = 1.4$$. As well as, the used initial condition is $$X_{0} = [x_{0}, y_{0}, z_{0}, \phi _{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 is shown in Fig. 12 for the solution $$\phi (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 $$\phi (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 fractional-order 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.

## Availability of data and materials

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

## References

1. Tuan, H.T., Ganji, R.M., Jafari, H.: A numerical study of fractional rheological models and fractional Newell–Whitehead–Segel equation with non-local and non-singular kernel. Chin. J. Phys. 68, 224–234 (2020). https://doi.org/10.1016/j.cjph.2020.08.019

2. Baleanu, D., Etemad, S., Rezapour, S.: A hybrid Caputo fractional modeling for thermostat with hybrid boundary value conditions. Bound. Value Probl. 2020, 64 (2020). https://doi.org/10.1186/s13661-020-01361-0

3. Baleanu, D., Mohammadi, H., Rezapour, S.: Analysis of the model of HIV-1 infection of CD4+ T-cell with a new approach of fractional derivative. Adv. Differ. Equ. 2020, 71 (2020). https://doi.org/10.1186/s13662-020-02544-w

4. Mohammadi, H., Kumar, S., Rezapour, S., Etemad, S.: A theoretical study of the Caputo-Fabrizio fractional modeling for hearing loss due to Mumps virus with optimal control. Chaos Solitons Fractals 144, 110668 (2021). https://doi.org/10.1016/j.chaos.2021.110668

5. Nikan, O., Jafari, H., Golbabai, A.: Numerical analysis of the fractional evolution model for heat flow in materials with memory. Alex. Eng. J. 59(4), 2627–2637 (2020). https://doi.org/10.1016/j.aej.2020.04.026

6. Ganji, R.M., Jafari, H., Baleanu, D.: A new approach for solving multi variable orders differential equations with Mittag-Leffler kernel. Chaos Solitons Fractals 130, 109405 (2020). https://doi.org/10.1016/j.chaos.2019.109405

7. Abdeljawad, T., Agarwal, R.P., Karapinar, E., Kumari, P.S.: Solutions of the nonlinear integral equation and fractional differential equation using the technique of a fixed point with a numerical experiment in extended b-metric space. Symmetry 11(5), 686 (2019). https://doi.org/10.3390/sym11050686

8. Alqahtani, B., Fulga, A., Karapinar, E.: Common fixed point results on an extended b-metric space. J. Inequal. Appl. 2018, 158 (2018). https://doi.org/10.1186/s13660-018-1745-4

9. Nguyen, P.D.: Note on a Allen-Cahn equation with Caputo-Fabrizio derivative. Results Nonlinear Anal. 4(3), 179–185 (2021). https://doi.org/10.53006/rna.962068

10. Sabetghadam, F., Masiha, H.P.: Common fixed points for generalized φ-pair mappings on cone metric spaces. Fixed Point Theory Appl. 2010, Article ID 718340 (2010). https://doi.org/10.1155/2010/718340

11. Sabetghadam, F., Masiha, H.P.: Fixed-point results for multi-valued operators in quasi-ordered metric spaces. Appl. Math. Lett. 25(11), 1856–1861 (2012). https://doi.org/10.1016/j.aml.2012.02.046

12. Afshari, H., Shojaat, H., Moradi, M.S.: Existence of the positive solutions for a tripled system of fractional differential equations via integral boundary conditions. Results Nonlinear Anal. 4(3), 186–193 (2021). https://doi.org/10.53006/rna.938851

13. Shojaat, H., Afshari, H., Asgari, M.S.: A new class of mixed monotone operators with concavity and applications to fractional differential equations. TWMS J. Appl. Eng. Math. 11(1), 122–133 (2021)

14. Lazreg, J.E., Abbas, S., Benchohra, M., Karapinar, E.: Impulsive Caputo-Fabrizio fractional differential equations in b-metric spaces. Open Math. 19(1), 363–372 (2021). https://doi.org/10.1515/math-2021-0040

15. Karapinar, E., Fulga, A., Rashid, M., Shahid, L., Aydi, H.: Large contractions on quasi-metric spaces with an application to nonlinear fractional differential equations. Mathematics 7(5), 444 (2019). https://doi.org/10.3390/math7050444

16. Adiguzel, R.S., Aksoy, U., Karapinar, E., Erhan, I.M.: On the solution of a boundary value problem associated with a fractional differential equation. Math. Methods Appl. Sci. (2020). https://doi.org/10.1002/mma.6652

17. Baleanu, D., Jajarmi, A., Mohammadi, H., Rezapour, S.: A new study on the mathematical modelling of human liver with Caputo-Fabrizio fractional derivative. Chaos Solitons Fractals 134, 109705 (2020). https://doi.org/10.1016/j.chaos.2020.109705

18. Adiguzel, R.S., Aksoy, U., Karapinar, E., Erhan, I.M.: On the solutions of fractional differential equations via Geraghty type hybrid contractions. Int. J. Appl. Comput. Math. 20(2), 313–333 (2021)

19. Adiguzel, R.S., Aksoy, U., Karapinar, E., Erhan, I.M.: Uniqueness of solution for higher-order nonlinear fractional differential equations with multi-point and integral boundary conditions. Rev. R. Acad. Cienc. Exactas Fís. Nat., Ser. A Mat. 2021, 155 (2021). https://doi.org/10.1007/s13398-021-01095-3

20. Alqahtani, B., Aydi, H., Karapinar, E., Rakocevic, V.: A solution for Volterra fractional integral equations by hybrid contractions. Mathematics 7(8), 694 (2019). https://doi.org/10.3390/math7080694

21. Masiha, H.P., Sabetghadam, F., Shahzad, N.: Fixed point theorems in partial metric spaces with an application. Filomat 27(4), 617–624 (2013)

22. Sabetghadam, F., Masiha, H.P., Altun, I.: Fixed-point theorems for integral-type contractions on partial metric spaces. Ukr. Math. J. 68, 940–949 (2016). https://doi.org/10.1007/s11253-016-1267-5

23. Rezapour, S., Mohammadi, H., Jajarmi, A.: A new mathematical model for Zika virus transmission. Adv. Differ. Equ. 2020, 589 (2020). https://doi.org/10.1186/s13662-020-03044-7

24. Aydogan, S.M., Baleanu, D., Mousalou, A., Rezapour, S.: On high order fractional integro-differential equations including the Caputo-Fabrizio derivative. Bound. Value Probl. 2018, 90 (2018). https://doi.org/10.1186/s13661-018-1008-9

25. Baleanu, D., Mousalou, A., Rezapour, S.: On the existence of solutions for some infinite coefficient-symmetric Caputo-Fabrizio fractional integro-differential equations. Bound. Value Probl. 2017, 145 (2017). https://doi.org/10.1186/s13661-017-0867-9

26. Baleanu, D., Rezapour, S., Saberpour, Z.: On fractional integro-differential inclusions via the extended fractional Caputo-Fabrizio derivation. Bound. Value Probl. 2019, 79 (2019). https://doi.org/10.1186/s13661-019-1194-0

27. Samei, M.E., Rezapour, S.: On a system of fractional q-differential inclusions via sum of two multi-term functions on a time scale. Bound. Value Probl. 2020, 135 (2020). https://doi.org/10.1186/s13661-020-01433-1

28. Tuan, N.H., Mohammadi, H., Rezapour, S.: A mathematical model for COVID-19 transmission by using the Caputo fractional derivative. Chaos Solitons Fractals 140, 110107 (2020). https://doi.org/10.1016/j.chaos.2020.110107

29. Rezapour, S., Samei, M.E.: On the existence of solutions for a multi-singular pointwise defined fractional q-integro-differential equation. Bound. Value Probl. 2020, 38 (2020). https://doi.org/10.1186/s13661-020-01342-3

30. Alrabaiah, H., Ahmad, I., Shah, K., Ur Rahman, G.: Qualitative analysis of nonlinear coupled pantograph differential equations of fractional order with integral boundary conditions. Bound. Value Probl. 2020, 138 (2020). https://doi.org/10.1186/s13661-020-01432-2

31. Deressa, C.T., Duressa, G.F.: Modeling and optimal control analysis of transmission dynamics of COVID-19: the case of Ethiopia. Alex. Eng. J. 60(1), 719–732 (2021)

32. Thabet, S.T.M., Etemad, S., Rezapour, S.: On a coupled Caputo conformable system of pantograph problems. Turk. J. Math. 45(1), 496–519 (2021). https://doi.org/10.3906/mat-2010-70

33. Li, T., Wang, Y.: Stability of a class of fractional-order nonlinear systems. Discrete Dyn. Nat. Soc. 2014, Article ID 724270 (2014). https://doi.org/10.1155/2014/724270

34. Khan, H., Gomez-Aguilar, J.F., Alkhazzan, A., Khan, A.: A fractional order HIV-TB coinfection model with nonsingular Mittag-Leffler law. Math. Methods Appl. Sci. 43(6), 3786–3806 (2020). https://doi.org/10.1002/mma.6155

35. Deressa, C.T., Duressa, G.F.: Analysis of Atangana-Baleanu fractional-order SEAIR epidemic model with optimal control. Adv. Differ. Equ. 2021, 174 (2021). https://doi.org/10.1186/s13662-021-03334-8

36. Atangana, A., Gomez-Aguilar, J.F.: Hyperchaotic behaviour obtained via a nonlocal operator with exponential decay and Mittag-Leffler laws. Chaos Solitons Fractals 102, 285–294 (2017). https://doi.org/10.1016/j.chaos.2017.03.022

37. Atangana, A., Owolabi, K.M.: New numerical approach for fractional differential equations. Math. Model. Nat. Phenom. 13(1), 3 (2018). https://doi.org/10.1051/mmnp/2018010

38. Trzaska, Z.W.: Matlab solutions of chaotic fractional order circuits. Engineering education and research using MATLAB. IntechOpen (2011). https://doi.org/10.5772/23144

39. Li, Y., Huang, X., Guo, M.: The generation, analysis and circuit implementation of a new memristor based chaotic system. Math. Probl. Eng. 2013, Article ID 398306 (2013). https://doi.org/10.1155/2013/398306

40. Sene, N.: Study of a fractional-order chaotic system represented by the Caputo operator. Complexity 2021, Article ID 5534872, 1–20 (2021). https://doi.org/10.1155/2021/5534872

41. Barrow-Green, J.: Poincare and the Three-Body Problem. Am. Math. Soc., New York (1996)

42. Lorenz, E.: Deterministic non-periodic flow. J. Atmos. Sci. 20(2), 130–141 (1963). https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2

43. Sene, N.: Analysis of a fractional-order chaotic system in the context of the Caputo fractional derivative via bifurcation and Lyapunov exponents. J. King Saud Univ., Sci. 33(1), 101275 (2021). https://doi.org/10.1016/j.jksus.2020.101275

44. Sene, N.: Introduction to the fractional-order chaotic system under fractional operator in Caputo sense. Alex. Eng. J. 60(4), 3997–4014 (2021). https://doi.org/10.1016/j.aej.2021.02.056

45. Sene, N.: Qualitative analysis of class of fractional-order chaotic system via bifurcation and Lyapunov exponents notions. J. Math. 2021, Article ID 5548569 (2021). https://doi.org/10.1155/2021/5548569

46. Sene, N.: Mathematical views of the fractional Chua’s electrical circuit described by the Caputo-Liouville derivative. Rev. Mex. Fis. 67(1), 91–99 (2021). https://doi.org/10.31349/REVMEXFIS.67.91

47. Petras, I.: Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation. Springer, Berlin (2011)

48. Sun, K., Wang, X., Sprott, J.C.: Bifurcations and chaos in fractional-order simplified Lorenz system. Int. J. Bifurc. Chaos 20(4), 1209–1219 (2010). https://doi.org/10.1142/S0218127410026411

49. Naik, P.A., Owolabi, K.M., Yavuz, M., Zu, J.: Chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells. Chaos Solitons Fractals 140, 110272 (2020). https://doi.org/10.1016/j.chaos.2020.110272

50. Owolabi, K.M., Atangana, A.: On the formulation of Adams-Bashforth scheme with Atangana-Baleanu-Caputo fractional derivative to model chaotic problems. Chaos 29(2), 023111 (2019). https://doi.org/10.1063/1.5085490

51. Owolabi, K.M., Atangana, A.: Modelling and simulation of a dynamical system with the Atangana-Baleanu fractional derivative. Eur. Phys. J. Plus 133, 15 (2018). https://doi.org/10.1140/epjp/i2018-11863-9

52. Buscarino, A., Fortuna, L., Frasca, M., Gambuzza, L.V.: A chaotic circuit based on Hewlett-Packard memristor. Chaos, Interdiscip. J. Nonlinear Sci. 22, 023136 (2012). https://doi.org/10.1063/1.4729135

53. Ruan, J., Sun, K., Mou, J., He, S., Zhang, L.: Fractional-order simplest memristor-based chaotic circuit with new derivative. Eur. Phys. J. Plus 133, 3 (2018). https://doi.org/10.1140/epjp/i2018-11828-0

54. Wang, R., Li, M., Gao, Z., Sun, H.: A new memristor-based 5d chaotic system and circuit implementation. Complexity 2018, Article ID 6069401 1–12 (2018). https://doi.org/10.1155/2018/6069401

55. Toufik, M., Atangana, A.: New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models. Eur. Phys. J. Plus 132, 444 (2017). https://doi.org/10.1140/epjp/i2017-11717-0

56. Danca, M.F., Kuznetsov, N.: Matlab code for Lyapunov exponents of fractional-order systems. Int. J. Bifurc. Chaos 28(5), 1850067 (2018). https://doi.org/10.1142/S0218127418500670

57. Muthuswamy, B.: Implementing memristor based chaotic circuits. Int. J. Bifurc. Chaos 20(5), 1335–1350 (2010). https://doi.org/10.1142/S0218127410026514

58. Atangana, A., Araz, S.I.: New numerical method for ordinary differential equations: Newton polynomial. J. Comput. Appl. Math. 372, 112622 (2020). https://doi.org/10.1016/j.cam.2019.112622

59. Atangana, A., Baleanu, D.: New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Therm. Sci. 20(2), 763–769 (2016)

60. Abdeljawad, T., Baleanu, D.: Integration by parts and its applications of a new nonlocal fractional derivative with Mittag-Leffler nonsingular kernel. J. Nonlinear Sci. Appl. 10, 1098–1107 (2017). https://doi.org/10.22436/jnsa.010.03.20

61. Granas, A., Dugundji, J.: Fixed Point Theory. Springer, New York (2003)

## Acknowledgements

The first author was supported by Jimma University. The second and third authors were supported by Azarbaijan Shahid Madani University. The authors express their gratitude to dear unknown referees for their helpful suggestions which improved the final version of this paper.

Not applicable.

## Author information

Authors

### Contributions

The authors declare that the study was realized in collaboration with equal responsibility. All authors read and approved the final manuscript.

### Corresponding author

Correspondence to Shahram Rezapour.

## Ethics declarations

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

Not applicable.

## Rights and permissions
{"ft_lang_label":"__label__en","ft_lang_prob":0.82397735,"math_prob":0.9996437,"size":49116,"snap":"2022-27-2022-33","text_gpt3_token_len":13707,"char_repetition_ratio":0.16124369,"word_repetition_ratio":0.08184545,"special_character_ratio":0.30197898,"punctuation_ratio":0.17888299,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9999782,"pos_list":[0,1,2],"im_url_duplicate_count":[null,2,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-06-27T17:35:07Z\",\"WARC-Record-ID\":\"<urn:uuid:165e5971-169f-4f69-aa7d-3ab3d3cb836c>\",\"Content-Length\":\"426971\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:f834f16c-84be-4b9b-963c-8b5de7090963>\",\"WARC-Concurrent-To\":\"<urn:uuid:90902727-7bb3-49f2-9018-fb58312e84a4>\",\"WARC-IP-Address\":\"146.75.32.95\",\"WARC-Target-URI\":\"https://advancesindifferenceequations.springeropen.com/articles/10.1186/s13662-021-03600-9\",\"WARC-Payload-Digest\":\"sha1:BXDLBBXE32NG53PQFB233Z6LNXUNS55T\",\"WARC-Block-Digest\":\"sha1:VTYHYIRGPXNNAXCZ56X32ILA2TX6YOA3\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-27/CC-MAIN-2022-27_segments_1656103337962.22_warc_CC-MAIN-20220627164834-20220627194834-00119.warc.gz\"}"}