# Multistability and Formation of Spiral Waves in a Fractional-Order Memristor-Based Hyperchaotic Lü System with No Equilibrium Points

Authors: Bo Yan; Shaobo He; Shaojie Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Multistablity analysis and formation of spiral wave in the fractional-order nonlinear systems is a recent hot topic. In this paper, dynamics, coexisting attractors, complexity, and synchronization of the fractional-order memristor-based hyperchaotic Lü system are investigated numerically by means of bifurcation diagram, Lyapunov exponents (LEs), chaos diagram, and sample entropy (SampEn) algorithm. The results show that the system has rich dynamics and high complexity. Meanwhile, coexisting attractors in the system are observed and hidden dynamics are illustrated by changing the initial conditions. Finally, the network based on the system is built, and the emergence of spiral waves is investigated and chimera states are observed.


---

## Body

## 1. Introduction

A new two-terminal circuit element called memristor characterized by a relationship between the charge and the flux was postulated as the fourth basic circuit element by Chua [1] in 1971. Memristor is short for memory resistor, and its memristance depends on how much electric charge flows through it in a particular direction [2]. Until 2008, researchers at HP Labs firstly published a paper in Nature and reported the realization of the memristor [3]. Since then, applications of the memristor have attracted lots of attention of scholars and a big progress in the fields such as secure communications [4, 5], neural network [6, 7], chaos [8, 9], and simulation of neural network behaviour [10, 11] and the memristor oscillators’ behaviour [12, 13] has been made.Memristor is a nonlinear electronic element. When the memristor is introduced, the chaotic oscillation of the circuit could be more complex due to its nonlinearity. So the research of memristor-based chaotic systems and its dynamic characteristics has become a new hot spot. In recent years, relevant research results have been obtained [9, 14–16]. Fractional calculus is a topic which is more than 300 years old [17]. Lots of physical phenomena show that the fractional-order model can describe the characteristics of the actual system more accurately than the classical integer-order model [18]. Based on this consideration, the application of fractional-order system has attracted more and more researchers’ attention [19, 20]. When the fractional calculus is introduced into the memristor-based chaotic circuit, a more accurate model for the real situation is obtained. Nowadays, many research studies on chaotic dynamics of the fractional-order memristor-based chaotic system have been conducted [21–25]. Among these studies, Rajagopal et al. [24] analyzed the dynamic properties of the fractional-order memristor system with no equilibrium points through Lyapunov exponents and phase diagrams for the first time. Prakash et al. [25] found rich chaotic behaviours of a novel fractional-order memristor-based chaotic jerk system with no equilibrium points by bifurcation diagram, attractors, and instantaneous phase plots. As a result, the analysis of chaotic dynamical behaviours of the fractional-order memristor-based chaotic systems without equilibrium points is an interesting topic.Recently, Qiao et al. [26] designed a novel improved memristor-based hyperchaotic Lü circuit without any equilibrium points by adding a linear feedback term, a constant term, and a memristor to the classical Lü chaotic system [27]. The coexisting hidden dynamical behaviours of this system are analyzed. However, the fractional-order chaotic dynamics of the system have not been investigated. Moreover, collective behaviours in the coupled nonlinear network have aroused much research interests, and these networks usually consist of nonlinear chaotic systems [28, 29], neural models [30, 31], and neural models with memristors [32, 33]. In fact, these networks can present different kinds of dynamical behaviours such as synchronization and chimera states. However, there are few reports on the formation of spiral waves in the network of fractional-order memristor-based hyperchaotic systems with no equilibrium points. Thus, the dynamics in the network of the new fractional-order memristor-based Lü system need further investigations.The rest of this paper is organized as follows. Section2 presents the details of the solution algorithm and mathematical model of the fractional-order memristor-based hyperchaotic Lü system. Section 3 investigates dynamics properties and complexity of the system by using bifurcation, Lyapunov exponents, and SampEn algorithm. The multiple coexisting hidden attractors are also illustrated with different initial conditions. Section 4 investigates the synchronization of the system by the ring network structure of the neuron network. Section 5 gives the conclusion.


## 2. The Fractional-Order Memristor-Based Hyperchaotic Lü System

### 2.1. The Mathematical Model

The most commonly used memristor-based dynamical system [16] is described as(1)i=α+βφ2u,φ˙=u,where u and i are the input voltage variable and the input current variable of the memristor, respectively, φ is the internal magnetic flux, α and β are two positive constants, and α + βφ2 is a memductance function describing the flux-dependent rate of change of charge. Qiao et al. [26] designed the memristor-based hyperchaotic Lü system, and its equations are given by(2)x˙=ay−x+gα+βw2y,y˙=by−xz−dx+e,z˙=xy−cz,w˙=y.Here, by introducing the fractional calculus to this system, a novel fractional-order memristor-based hyperchaotic Lü system is constructed, and the mathematical model is given by(3)Dt0qx=ay−x+gα+βw2y,Dt0qy=by−xz−dx+e,Dt0qz=xy−cz,Dt0qw=y,where Dt0q is the q-order Caputo differential operator [34], x, y, z, and w are state variables, a, b, c, d, and e are the control parameters of classical Lü system [27], and g is the gain of the memristor. The Caputo fractional-order definition is given in Definition 1.Definition 1.
(see [34]). The Caputo fractional-order derivative definition is given by(4)Dt0qxt=1Γ1−q∫t0tx˙τt−τdτ,0<q<1,x˙t,q=1,where q ∈ R+ and Γ⋅ is the Gamma function.
In order to analyze the dynamic and complexity characteristics of this fractional-order memristor-based hyperchaotic Lü system, we fix the system parameters asa = 36, b = 20, c = 3, d = 5, e = 0.1, α = 4, and β = 0.18 and select the gain parameter g as the only adjustable parameter of the system.
It can be obtained by (3):(5)∇qV=∂Dt0qx∂x+∂Dt0qy∂y+∂Dt0qy∂z+∂Dt0qw∂w=−a+b−c=−19<0,where ∇qV<0 means that this fractional-order memristor-based hyperchaotic Lü system is dissipative. It also indicates that the improvement of Lü system and the introduction of memristor and fractional calculus do not change the dissipation proprieties of the original classical Lü system.
LetDt0qx=Dt0qy=Dt0qz=Dt0qw=0; obviously, (3) has no real solutions and the fractional-order memristor-based hyperchaotic Lü system has no equilibrium. As shown in Ref [26], the integer-order hyperchaotic Lü system has hidden attractors. According to the definition of the hidden attractors [35], the generated attractors and limit cycles of system (3) are all hidden because this fractional-order system investigated in this paper also has no equilibrium points.


### 2.2. The Numerical Analysis Methods

The predictor-corrector method is an effective method for solving fractional-order equations. According to Kai et al. [36], it is also called the Adams–Bashforth–Moulton algorithm. Suppose that the fractional-order nonlinear system is defined as(6)Dt0qxt=fx,xt,0≤t≤T,xk0=x0k,k=0,1,2,…,q−1,where xt=xt,yt,zt,wt, xk=x0k is the initial condition of the system, 0 < q ≤ 1, ⌈⋅⌉ is the ceil function, and Dt0qxt is the Caputo derivative. Formula (6) is equivalent to the Volterra integral equation [33, 37]:(7)xt=∑k=0q−1x0ktkk!+1Γα∫0Tt−τq−1fτ,xτdτ.Leth = T/N, tj = jh (j = 0, 1, …, N ∈ Z+), and T be the upper limit for solving this integral equation. Then, the correction formula of formula (7) is(8)xhtn+1=∑k=0q−1x0ktn+1kk!+hqΓq+2ftn+1,xhptn+1+hqΓq+2∑j=0naj,n+1ftj,xhtj,in which(9)xhptn+1=∑k=0q−1x0ktn+1kk!+1Γq∑j=0nbj,n+1ftj,xhtj,bj,n+1=hqqn−j+1q−n−jq,0≤j≤n,aj,n+1=nq+1−n−qn+1q,j=0,n−j+2q+1+n−jq+1−2n−j+1q+1,1≤j≤n.The utilized source Matlab code of the predictor-corrector method can be downloaded online (http://www.sourcecodeonline.com/details/predictorcorrector_pece_method_for_fractional_differential_equations.html); the function name is “FDE12.m.” In [38], the method about how to use “FDE12.m” to solve fraction-order system is described in detail. The LEs (Lyapunov exponents) [39] of system (6) can be calculated as follows. Firstly, suppose Λk (t) is the eigenvalue of fx,xt,k=1,2,…,d,whered is the dimension of the system. The orbit xt can be given in terms of the tangent linear extension:(10)Dt0qv=∇fxtT⋅v,where v is a vector evolving in the tangent bundle of the manifold and the superscript “T” denotes the transpose. Then, the LEs can be defined as the limit eigenvalues of the symmetric operator:(11)Λkt=1tAt⋅ATt,where A (t) is the exponential of fx,xt computed along the trajectory xt. Thus, we can obtain Λ1t≤Λ2t≤⋯≤Λdt. Finally, The LEs of the trajectory xt obtained by system (6) can be determined by(12)λk=limt⟶∞sup1tlnΛkt.


## 2.1. The Mathematical Model

The most commonly used memristor-based dynamical system [16] is described as(1)i=α+βφ2u,φ˙=u,where u and i are the input voltage variable and the input current variable of the memristor, respectively, φ is the internal magnetic flux, α and β are two positive constants, and α + βφ2 is a memductance function describing the flux-dependent rate of change of charge. Qiao et al. [26] designed the memristor-based hyperchaotic Lü system, and its equations are given by(2)x˙=ay−x+gα+βw2y,y˙=by−xz−dx+e,z˙=xy−cz,w˙=y.Here, by introducing the fractional calculus to this system, a novel fractional-order memristor-based hyperchaotic Lü system is constructed, and the mathematical model is given by(3)Dt0qx=ay−x+gα+βw2y,Dt0qy=by−xz−dx+e,Dt0qz=xy−cz,Dt0qw=y,where Dt0q is the q-order Caputo differential operator [34], x, y, z, and w are state variables, a, b, c, d, and e are the control parameters of classical Lü system [27], and g is the gain of the memristor. The Caputo fractional-order definition is given in Definition 1.Definition 1.
(see [34]). The Caputo fractional-order derivative definition is given by(4)Dt0qxt=1Γ1−q∫t0tx˙τt−τdτ,0<q<1,x˙t,q=1,where q ∈ R+ and Γ⋅ is the Gamma function.
In order to analyze the dynamic and complexity characteristics of this fractional-order memristor-based hyperchaotic Lü system, we fix the system parameters asa = 36, b = 20, c = 3, d = 5, e = 0.1, α = 4, and β = 0.18 and select the gain parameter g as the only adjustable parameter of the system.
It can be obtained by (3):(5)∇qV=∂Dt0qx∂x+∂Dt0qy∂y+∂Dt0qy∂z+∂Dt0qw∂w=−a+b−c=−19<0,where ∇qV<0 means that this fractional-order memristor-based hyperchaotic Lü system is dissipative. It also indicates that the improvement of Lü system and the introduction of memristor and fractional calculus do not change the dissipation proprieties of the original classical Lü system.
LetDt0qx=Dt0qy=Dt0qz=Dt0qw=0; obviously, (3) has no real solutions and the fractional-order memristor-based hyperchaotic Lü system has no equilibrium. As shown in Ref [26], the integer-order hyperchaotic Lü system has hidden attractors. According to the definition of the hidden attractors [35], the generated attractors and limit cycles of system (3) are all hidden because this fractional-order system investigated in this paper also has no equilibrium points.


## 2.2. The Numerical Analysis Methods

The predictor-corrector method is an effective method for solving fractional-order equations. According to Kai et al. [36], it is also called the Adams–Bashforth–Moulton algorithm. Suppose that the fractional-order nonlinear system is defined as(6)Dt0qxt=fx,xt,0≤t≤T,xk0=x0k,k=0,1,2,…,q−1,where xt=xt,yt,zt,wt, xk=x0k is the initial condition of the system, 0 < q ≤ 1, ⌈⋅⌉ is the ceil function, and Dt0qxt is the Caputo derivative. Formula (6) is equivalent to the Volterra integral equation [33, 37]:(7)xt=∑k=0q−1x0ktkk!+1Γα∫0Tt−τq−1fτ,xτdτ.Leth = T/N, tj = jh (j = 0, 1, …, N ∈ Z+), and T be the upper limit for solving this integral equation. Then, the correction formula of formula (7) is(8)xhtn+1=∑k=0q−1x0ktn+1kk!+hqΓq+2ftn+1,xhptn+1+hqΓq+2∑j=0naj,n+1ftj,xhtj,in which(9)xhptn+1=∑k=0q−1x0ktn+1kk!+1Γq∑j=0nbj,n+1ftj,xhtj,bj,n+1=hqqn−j+1q−n−jq,0≤j≤n,aj,n+1=nq+1−n−qn+1q,j=0,n−j+2q+1+n−jq+1−2n−j+1q+1,1≤j≤n.The utilized source Matlab code of the predictor-corrector method can be downloaded online (http://www.sourcecodeonline.com/details/predictorcorrector_pece_method_for_fractional_differential_equations.html); the function name is “FDE12.m.” In [38], the method about how to use “FDE12.m” to solve fraction-order system is described in detail. The LEs (Lyapunov exponents) [39] of system (6) can be calculated as follows. Firstly, suppose Λk (t) is the eigenvalue of fx,xt,k=1,2,…,d,whered is the dimension of the system. The orbit xt can be given in terms of the tangent linear extension:(10)Dt0qv=∇fxtT⋅v,where v is a vector evolving in the tangent bundle of the manifold and the superscript “T” denotes the transpose. Then, the LEs can be defined as the limit eigenvalues of the symmetric operator:(11)Λkt=1tAt⋅ATt,where A (t) is the exponential of fx,xt computed along the trajectory xt. Thus, we can obtain Λ1t≤Λ2t≤⋯≤Λdt. Finally, The LEs of the trajectory xt obtained by system (6) can be determined by(12)λk=limt⟶∞sup1tlnΛkt.


## 3. Numerical Analysis

### 3.1. Bifurcations and Lyapunov Exponents

Dynamics of the system with the variation ofg and q are analyzed by means of bifurcation diagram, LEs, and phase diagrams. Three cases are presented.Case 1. Fixq = 0.998: let the gain parameter g vary from 0 to 20 with step size of 0.04 and the initial condition be x0;y0;z0;w0=1;0;1;0. Bifurcation diagram and LEs are illustrated in Figure 1. It is shown in Figure 1 that the system is periodic for g∈0.3,2.1 and [3.5, 7.8], is chaotic for g∈0,0.3, [2.1, 3.5], and [14.5, 15.2], is hyperchaotic for g∈7.8,14.5, and is convergent for g∈15.2,20. In this case, chaotic and periodic states are alternate in the interval g∈0,7.8. It shows the rich dynamics in the system.Case 2. Fixg=10: vary the derivative order q from 0.94 to 1 with step size of 1.2 × 10−4. Set the initial condition as x0;y0;z0;w0=1;0;1;0. Dynamical analysis results are shown in Figure 2. It is shown in Figure 2 that the system is hyperchaotic when q > 0.972. Meanwhile, when q ∈ [0.965, 0.972], the system vibrates between chaotic and hyperchaotic states. For the rest values of q, the system is convergent.Case 3: Vary parameterg from 1 to 16 with step size of 0.04 and increase derivative order q from 0.94 to 1 with step size 1.2 × 10−4. Thus, the parameter space is divided as a 400 × 500 grid. As with above cases, the initial condition is given by x0;y0;z0;w0=1;0;1;0. Maximum Lyapunov exponent-based contour plot in the parameter plane g-q is shown in Figure 3. Chaotic region is observed in a triangular area located in the top left of the q-g parameter plane. When the derivative order q takes different values, the system has relative larger maximum Lyapunov exponents with smaller g. Meanwhile, the maximum Lyapunov exponents of the system decrease with the increase of derivative order q. Finally, it can be found out that the minimum order for chaos is q = 0.927.Figure 1
Dynamical behaviours with derivative orderq = 0.998 and varying memristor gain g. (a) Bifurcation diagram. (b) LEs.
(a)(b)Figure 2
Dynamical behaviours with memristor gaing=10 and varying derivative order q. (a) Bifurcation diagram. (b) LEs.
(a)(b)Figure 3
Maximum Lyapunov exponent contour plot in theg-q parameter plane.Phase portraits of the system with the initial conditionx0;y0;z0;w0=1;0;1;0 and different g and q are plotted in Figure 4. Four sets of parameters are given, and they are q = 0.995 and g=0.3 for Figure 4(a), q = 0.995 and g=1 for Figure 4(b), q = 0.995 and g=3 for Figure 4(c), and q = 0.995 and g=14 for Figure 4(d). Different states including periodic circle, chaotic attractor, and hyperchaotic attractor are observed in the system. According to the above analyses, the conclusion can be drawn that the fractional-order memristor-based hyperchaotic Lü system has rich dynamical behaviours.Figure 4
Phase diagrams with derivative orderq = 0.995 and different memristor gain g: (a) g=0.3; (b) g=1; (c) g=3; (d) g=14.
(a)(b)(c)(d)


### 3.2. SampEn Complexity Analysis

Sample entropy (SampEn) [40, 41] is a measure of the probability size of a new model generated by a time series when its dimension changes. The greater the complexity of the time series is, the greater the probability for a new model that the time series will produce and the greater the entropy holds. The SampEn algorithm can be briefly described as follows:Step 1. Given a time series {x (i), i = 0, 1, 2, …, N − 1}, if m represents embedding dimension, then the m-dimensional vectors can be expressed as(13)Xi=xi,xi+1,…,xi+m−1,i=1,2,…,N−m.Step 2. Define the Euclidean distance betweenX (i) and X (j) as the maximum absolute difference between their corresponding scalar elements, i.e.,(14)dXi,Xj=Maxk=0⟶m−1Xi+k−Xj+k.Step 3. Define the criterion of similarityr as a positive real number. Count the number of dXi,Xj<r and its ratio to total distance N − m denoted as Brm (i), i = 0, 1, 2, …, N − m, j = 0, 1, 2, …, N − m, and i ≠ j. Then, the average value of Bm (i) can be calculated as(15)Bmr=1N−m+1∑i=0N−mBimr.Step 4. Similarly, changem to m + 1 and repeat Steps 1 to 3; B(m+1) (r) also can be obtained as(16)Bm+1r=1N−m+1∑i=0N−mBim+1r.Finally, theoretically, the SampEn is defined as(17)SampEnm,r=−lnBm+1rBmr.The value of SampEn is related to the values ofm and r. Usually, m ranges from 2 to 5, and r takes 0.1 to 0.25 times of the standard deviation (SD) of the original time series. This article takes m = 2 and r = 0.2SD.SampEn complexity of the fractional-order memristor-based hyperchaotic Lü system with parameterg and q is analyzed, and the results are shown in Figure 5. Time series z is used for complexity measure and length of the time series is 10000 with the first 3000 points of data removed. The parameters used for Figure 5(a) are the same as that used for Figure 1. And the parameters used for Figure 5(b) are the same as that used for Figure 2. From Figure 5, we can see that higher complexity is found when the system is chaotic or hyperchaotic. This result agrees well with the analysis results by the corresponding maximum Lyapunov exponents. Comparing Figure 5(a) with Figure 1(b), it can be seen that the difference of the SampEn complexity is more obvious in Figure 5(a) than the difference of corresponding maximum Lyapunov exponents in Figure 1(b), which indicates that SampEn complexity is more convenient in distinguishing chaotic with nonchaotic states of the system.Figure 5
SampEn complexity analysis results with derivative orderq parameter g with initial condition x0;y0;z0;w0=1;0;1;0. (a) q = 0.998 and varying gain g. (b) g=10 and varying derivative order q.
(a)(b)Let parameterg vary from 0 to 16 with step size of 0.004, derivative order q vary from 0.94 to 1 with step size 1.2 × 10−4, and the initial condition be x0;y0;z0;w0=1;0;1;0. SampEn complexity in the parameter plane g-q is calculated and shown in Figure 6. Obviously, the higher complexity region in Figure 6 matches well with the larger maximum Lyapunov exponents area in Figure 3. However, unlike the maximum Lyapunov exponent plane, the SampEn complexity plane does not change gradually but is divided into three regions: chaos, periodic cycle, and convergence. The upper right corner in Figure 6 is periodic region. So, the SampEn complexity analysis results are more intuitive and obvious.Figure 6
SampEn complexity contour plots in theg-q parameter plane.


### 3.3. Coexisting Hidden Dynamics

Coexisting attractors can be observed in this system. Fixq = 0.998 and g=5; four phase diagrams are shown in different forms in Figure 7 because of the different initial conditions. When x0=y0=z0=w0=0 and x0=y0=z0=w0=1, the system is periodic but has different number of cycles. When x0=y0=z0=w0=3 and x0=y0=z0=w0=5, the system is chaotic, but Figure 7(d) shows the weak chaotic state. When g varies from 4.8 to 5, the bifurcation diagrams for the corresponding 4 different initial conditions are shown in Figure 8. The blue points are for x0=y0=z0=w0=0, red points are for x0=y0=z0=w0=1, blue points are for x0=y0=z0=w0=3, and magenta points are for x0=y0=z0=w0=5. It can be seem from Figure 8 that the system trajectory shows the multiple stable phenomena of coexistence for different periodic limit cycles, chaos, and weak chaos due to the different initial conditions.Figure 7
Coexisting attractors withq = 0.998, g=5, and different initial conditions. (a) x0=y0=z0=w0=0. (b) x0=y0=z0=w0=1. (c) x0=y0=z0=w0=3. (d) x0=y0=z0=w0=5.
(a)(b)(c)(d)Figure 8
Bifurcation diagram in the rangeg∈4.85 with q = 0.998 and different initial conditions.In order to further study the effect of initial conditions on the dynamics of the system, the SampEn complexity chaos diagram for different parameters is illustrated in Figure9. Comparing Figure 9(a) with Figure 9(b), it can be seen that the effect of initial conditions on system dynamics is multidimensional. Each initial component will play a role in the dynamics of the system. Comparing Figure 9(c) with Figure 9(d), we can find that when the derivative order q is different, the effect of initial value on system dynamics is also different. By comparing the four figures in Figure 9, it can be found that the SampEn complexity diagrams are different with different initial conditions and SampEn algorithm is an effective method for the analysis of the hidden dynamics in the system.Figure 9
SampEn complexity-based chaos diagram. (a)q = 0.998, g=5, w0=z0=1, x0 − y0 plane. (b) q = 0.998, g=5, x0=w0=1, y0 − z0 plane. (c) q = 0.94, g=0.1, w0=z0=1, x0 − y0 plane. (d) q = 0.97, g=0.1, w0=z0=1, x0 − y0 plane.
(a)(b)(c)(d)Compared with the basin attraction plots, the SampEn-based chaos diagrams illustrate the dynamics and complexity of the system with different initial conditions. If the complexity is different, the state of the system is different. Thus, it can show the multistability of the system in the initial plane through a complexity of view. Since complexity measure just needs a segment of time series, it is very convenient to detect the coexisting attractors in the multistability systems. However, it can only detect the coexisting states with different complexity. Usually, it means that there exists coexisting chaotic attractors and periodic circles.


## 3.1. Bifurcations and Lyapunov Exponents

Dynamics of the system with the variation ofg and q are analyzed by means of bifurcation diagram, LEs, and phase diagrams. Three cases are presented.Case 1. Fixq = 0.998: let the gain parameter g vary from 0 to 20 with step size of 0.04 and the initial condition be x0;y0;z0;w0=1;0;1;0. Bifurcation diagram and LEs are illustrated in Figure 1. It is shown in Figure 1 that the system is periodic for g∈0.3,2.1 and [3.5, 7.8], is chaotic for g∈0,0.3, [2.1, 3.5], and [14.5, 15.2], is hyperchaotic for g∈7.8,14.5, and is convergent for g∈15.2,20. In this case, chaotic and periodic states are alternate in the interval g∈0,7.8. It shows the rich dynamics in the system.Case 2. Fixg=10: vary the derivative order q from 0.94 to 1 with step size of 1.2 × 10−4. Set the initial condition as x0;y0;z0;w0=1;0;1;0. Dynamical analysis results are shown in Figure 2. It is shown in Figure 2 that the system is hyperchaotic when q > 0.972. Meanwhile, when q ∈ [0.965, 0.972], the system vibrates between chaotic and hyperchaotic states. For the rest values of q, the system is convergent.Case 3: Vary parameterg from 1 to 16 with step size of 0.04 and increase derivative order q from 0.94 to 1 with step size 1.2 × 10−4. Thus, the parameter space is divided as a 400 × 500 grid. As with above cases, the initial condition is given by x0;y0;z0;w0=1;0;1;0. Maximum Lyapunov exponent-based contour plot in the parameter plane g-q is shown in Figure 3. Chaotic region is observed in a triangular area located in the top left of the q-g parameter plane. When the derivative order q takes different values, the system has relative larger maximum Lyapunov exponents with smaller g. Meanwhile, the maximum Lyapunov exponents of the system decrease with the increase of derivative order q. Finally, it can be found out that the minimum order for chaos is q = 0.927.Figure 1
Dynamical behaviours with derivative orderq = 0.998 and varying memristor gain g. (a) Bifurcation diagram. (b) LEs.
(a)(b)Figure 2
Dynamical behaviours with memristor gaing=10 and varying derivative order q. (a) Bifurcation diagram. (b) LEs.
(a)(b)Figure 3
Maximum Lyapunov exponent contour plot in theg-q parameter plane.Phase portraits of the system with the initial conditionx0;y0;z0;w0=1;0;1;0 and different g and q are plotted in Figure 4. Four sets of parameters are given, and they are q = 0.995 and g=0.3 for Figure 4(a), q = 0.995 and g=1 for Figure 4(b), q = 0.995 and g=3 for Figure 4(c), and q = 0.995 and g=14 for Figure 4(d). Different states including periodic circle, chaotic attractor, and hyperchaotic attractor are observed in the system. According to the above analyses, the conclusion can be drawn that the fractional-order memristor-based hyperchaotic Lü system has rich dynamical behaviours.Figure 4
Phase diagrams with derivative orderq = 0.995 and different memristor gain g: (a) g=0.3; (b) g=1; (c) g=3; (d) g=14.
(a)(b)(c)(d)


## 3.2. SampEn Complexity Analysis

Sample entropy (SampEn) [40, 41] is a measure of the probability size of a new model generated by a time series when its dimension changes. The greater the complexity of the time series is, the greater the probability for a new model that the time series will produce and the greater the entropy holds. The SampEn algorithm can be briefly described as follows:Step 1. Given a time series {x (i), i = 0, 1, 2, …, N − 1}, if m represents embedding dimension, then the m-dimensional vectors can be expressed as(13)Xi=xi,xi+1,…,xi+m−1,i=1,2,…,N−m.Step 2. Define the Euclidean distance betweenX (i) and X (j) as the maximum absolute difference between their corresponding scalar elements, i.e.,(14)dXi,Xj=Maxk=0⟶m−1Xi+k−Xj+k.Step 3. Define the criterion of similarityr as a positive real number. Count the number of dXi,Xj<r and its ratio to total distance N − m denoted as Brm (i), i = 0, 1, 2, …, N − m, j = 0, 1, 2, …, N − m, and i ≠ j. Then, the average value of Bm (i) can be calculated as(15)Bmr=1N−m+1∑i=0N−mBimr.Step 4. Similarly, changem to m + 1 and repeat Steps 1 to 3; B(m+1) (r) also can be obtained as(16)Bm+1r=1N−m+1∑i=0N−mBim+1r.Finally, theoretically, the SampEn is defined as(17)SampEnm,r=−lnBm+1rBmr.The value of SampEn is related to the values ofm and r. Usually, m ranges from 2 to 5, and r takes 0.1 to 0.25 times of the standard deviation (SD) of the original time series. This article takes m = 2 and r = 0.2SD.SampEn complexity of the fractional-order memristor-based hyperchaotic Lü system with parameterg and q is analyzed, and the results are shown in Figure 5. Time series z is used for complexity measure and length of the time series is 10000 with the first 3000 points of data removed. The parameters used for Figure 5(a) are the same as that used for Figure 1. And the parameters used for Figure 5(b) are the same as that used for Figure 2. From Figure 5, we can see that higher complexity is found when the system is chaotic or hyperchaotic. This result agrees well with the analysis results by the corresponding maximum Lyapunov exponents. Comparing Figure 5(a) with Figure 1(b), it can be seen that the difference of the SampEn complexity is more obvious in Figure 5(a) than the difference of corresponding maximum Lyapunov exponents in Figure 1(b), which indicates that SampEn complexity is more convenient in distinguishing chaotic with nonchaotic states of the system.Figure 5
SampEn complexity analysis results with derivative orderq parameter g with initial condition x0;y0;z0;w0=1;0;1;0. (a) q = 0.998 and varying gain g. (b) g=10 and varying derivative order q.
(a)(b)Let parameterg vary from 0 to 16 with step size of 0.004, derivative order q vary from 0.94 to 1 with step size 1.2 × 10−4, and the initial condition be x0;y0;z0;w0=1;0;1;0. SampEn complexity in the parameter plane g-q is calculated and shown in Figure 6. Obviously, the higher complexity region in Figure 6 matches well with the larger maximum Lyapunov exponents area in Figure 3. However, unlike the maximum Lyapunov exponent plane, the SampEn complexity plane does not change gradually but is divided into three regions: chaos, periodic cycle, and convergence. The upper right corner in Figure 6 is periodic region. So, the SampEn complexity analysis results are more intuitive and obvious.Figure 6
SampEn complexity contour plots in theg-q parameter plane.


## 3.3. Coexisting Hidden Dynamics

Coexisting attractors can be observed in this system. Fixq = 0.998 and g=5; four phase diagrams are shown in different forms in Figure 7 because of the different initial conditions. When x0=y0=z0=w0=0 and x0=y0=z0=w0=1, the system is periodic but has different number of cycles. When x0=y0=z0=w0=3 and x0=y0=z0=w0=5, the system is chaotic, but Figure 7(d) shows the weak chaotic state. When g varies from 4.8 to 5, the bifurcation diagrams for the corresponding 4 different initial conditions are shown in Figure 8. The blue points are for x0=y0=z0=w0=0, red points are for x0=y0=z0=w0=1, blue points are for x0=y0=z0=w0=3, and magenta points are for x0=y0=z0=w0=5. It can be seem from Figure 8 that the system trajectory shows the multiple stable phenomena of coexistence for different periodic limit cycles, chaos, and weak chaos due to the different initial conditions.Figure 7
Coexisting attractors withq = 0.998, g=5, and different initial conditions. (a) x0=y0=z0=w0=0. (b) x0=y0=z0=w0=1. (c) x0=y0=z0=w0=3. (d) x0=y0=z0=w0=5.
(a)(b)(c)(d)Figure 8
Bifurcation diagram in the rangeg∈4.85 with q = 0.998 and different initial conditions.In order to further study the effect of initial conditions on the dynamics of the system, the SampEn complexity chaos diagram for different parameters is illustrated in Figure9. Comparing Figure 9(a) with Figure 9(b), it can be seen that the effect of initial conditions on system dynamics is multidimensional. Each initial component will play a role in the dynamics of the system. Comparing Figure 9(c) with Figure 9(d), we can find that when the derivative order q is different, the effect of initial value on system dynamics is also different. By comparing the four figures in Figure 9, it can be found that the SampEn complexity diagrams are different with different initial conditions and SampEn algorithm is an effective method for the analysis of the hidden dynamics in the system.Figure 9
SampEn complexity-based chaos diagram. (a)q = 0.998, g=5, w0=z0=1, x0 − y0 plane. (b) q = 0.998, g=5, x0=w0=1, y0 − z0 plane. (c) q = 0.94, g=0.1, w0=z0=1, x0 − y0 plane. (d) q = 0.97, g=0.1, w0=z0=1, x0 − y0 plane.
(a)(b)(c)(d)Compared with the basin attraction plots, the SampEn-based chaos diagrams illustrate the dynamics and complexity of the system with different initial conditions. If the complexity is different, the state of the system is different. Thus, it can show the multistability of the system in the initial plane through a complexity of view. Since complexity measure just needs a segment of time series, it is very convenient to detect the coexisting attractors in the multistability systems. However, it can only detect the coexisting states with different complexity. Usually, it means that there exists coexisting chaotic attractors and periodic circles.


## 4. Wave Formations of the Fractional-Order Ring Network

### 4.1. The Structure of the Ring Network

To further capture the complex dynamics of the system, we construct a ring network of nonlinear systems which are coupled locally to the nearby 3 nodes with a coupling strength ofD. The ring network structure of the neuron network is depicted in Figure 10 and can be numerically defined as(18)Dt0qxi=ayi−xi+gα+βwi2yi+D3∑x3∈Nodes−xi,Dt0qyi=byi−xizi−dzi+e+D3∑y3∈Nodes−yi,Dt0qzi=xiyi−czi+D3∑z3∈Nodes−zi,Dt0qwi=yi+D3∑w3∈Nodes−wi,where i = 1 to N. N is the number of nodes. The set nodes contain the three nodes near the target node. As shown in Figure 10, the nearby nodes of node 1 are node 2, node 3, and node 49, while the nearby nodes of node 36 are node 34, node 35, and node 38. To find the nodes, Algorithm 1 is used.Algorithm 1: The algorithm to find the nodes near node i.
Input:i, NOutput:Nodesifi = 1Nodes = [2, 3, N − 1];elseifi = 2Nodes = [1, N, 4];elseifi = N − 1Nodes = [1, N, N − 3];elseifi = NNodes = [2, N − 1, N − 2];elseifi%2 = 1Nodes = [i + 1, i − 2, i + 2];elseNodes = [i − 1, i − 2, i + 2];end ifFigure 10
Model of the ring network structure of the network (by Pajek 64 5.08).In order to derive the dynamics of the neuronal ring network, we also fix the system parameters asa = 36, b = 20, c = 3, d = 5, e = 0.1, α = 4, and β = 0.18 and vary the gain parameter g. In Figure 10, the built network has 50 nodes, and thus N = 50.


### 4.2. Numerical Simulation of the Fractional-Order Network

Spatiotemporal patterns of the model with the variation ofg and q are analyzed for different coupling strengths. Three cases are presented.Case 1. Fixq = 0.998, g=3, and various coupling strengths. The spatiotemporal patterns are shown in Figure 11. In these figures, the level of excitability of the network is illustrated by colours. We can see that when the coupling strength D = 2 and 4, only unsynchronized nodes are shown in Figures 11(a) and 11(b). When the coupling strength D increases to 6, some synchronized neurons are found in Figure 11(c), but there are still many unsynchronized nodes. By further increasing the coupling strength to D = 8, the systems in the network are synchronized locally.Case 2. Fixq = 0.998, g=8, and various coupling strengths. It can be seen from Figure 12 that the spatiotemporal behaviours for D = 2 show clear indication of coherent and incoherent nodes. But when the coupling strength is D = 4 and 6, some of the systems have synchronization with their nearby nodes. However, by further increasing the coupling strength to D = 8, all of the systems are synchronized.Case 3. Fixq = 0.97, g=3, and various coupling strengths. Figure 13 shows that the systems are not synchronized for the coupling strength D = 1. With the coupling strength increasing gradually, the number of synchronized systems increases. When the coupling strength D = 6, all nodes are in a local synchronization state between different nodes.Figure 11
Spatiotemporal patterns of the ring network for different coupling strengths of Case 1. (a)D = 2. (b) D = 4. (c) D = 6. (d) D = 8.
(a)(b)(c)(d)Figure 12
Spatiotemporal patterns of the ring network for different coupling strengths of Case 2. (a)D = 2. (b) D = 4. (c) D = 6. (d) D = 8.
(a)(b)(c)(d)Figure 13
Spatiotemporal patterns of the neuronal ring network for different coupling strengths of Case 3. (a)D = 1. (b) D = 2. (c) D = 4. (d) D = 6.
(a)(b)(c)(d)As shown in Figures11–13, the behaviours of the network are different with different coupling strengths. When the coupling strength is small, the network does not show synchronization. For instance, it is shown in Figures 11(a) and 11(b) that the network is not synchronized. When the coupling strength is large, the network can be synchronized. For example, Figure 12(d) shows that there exists synchronization. The chimera state means coexistence of synchronization and asynchronization in the network. According to Figures 11–13, there are different kinds of chimera state in the network. The network with different fractional-order derivative order has different dynamics. As for the coupling strength for chimera states in the network, it should be proper. According to the analysis results, when the coupling strength D∈2,4, the chimera states are observed. Moreover, the derivative order also affects the dynamics of the network. In Case 1, we fix q = 0.998 and in Case 3 we let q = 0.97, and it is shown in Figures 11 and 13 that the dynamics of the network are different. It is easier for the network with q = 0.97 to synchronize or to show chimera states.


## 4.1. The Structure of the Ring Network

To further capture the complex dynamics of the system, we construct a ring network of nonlinear systems which are coupled locally to the nearby 3 nodes with a coupling strength ofD. The ring network structure of the neuron network is depicted in Figure 10 and can be numerically defined as(18)Dt0qxi=ayi−xi+gα+βwi2yi+D3∑x3∈Nodes−xi,Dt0qyi=byi−xizi−dzi+e+D3∑y3∈Nodes−yi,Dt0qzi=xiyi−czi+D3∑z3∈Nodes−zi,Dt0qwi=yi+D3∑w3∈Nodes−wi,where i = 1 to N. N is the number of nodes. The set nodes contain the three nodes near the target node. As shown in Figure 10, the nearby nodes of node 1 are node 2, node 3, and node 49, while the nearby nodes of node 36 are node 34, node 35, and node 38. To find the nodes, Algorithm 1 is used.Algorithm 1: The algorithm to find the nodes near node i.
Input:i, NOutput:Nodesifi = 1Nodes = [2, 3, N − 1];elseifi = 2Nodes = [1, N, 4];elseifi = N − 1Nodes = [1, N, N − 3];elseifi = NNodes = [2, N − 1, N − 2];elseifi%2 = 1Nodes = [i + 1, i − 2, i + 2];elseNodes = [i − 1, i − 2, i + 2];end ifFigure 10
Model of the ring network structure of the network (by Pajek 64 5.08).In order to derive the dynamics of the neuronal ring network, we also fix the system parameters asa = 36, b = 20, c = 3, d = 5, e = 0.1, α = 4, and β = 0.18 and vary the gain parameter g. In Figure 10, the built network has 50 nodes, and thus N = 50.


## 4.2. Numerical Simulation of the Fractional-Order Network

Spatiotemporal patterns of the model with the variation ofg and q are analyzed for different coupling strengths. Three cases are presented.Case 1. Fixq = 0.998, g=3, and various coupling strengths. The spatiotemporal patterns are shown in Figure 11. In these figures, the level of excitability of the network is illustrated by colours. We can see that when the coupling strength D = 2 and 4, only unsynchronized nodes are shown in Figures 11(a) and 11(b). When the coupling strength D increases to 6, some synchronized neurons are found in Figure 11(c), but there are still many unsynchronized nodes. By further increasing the coupling strength to D = 8, the systems in the network are synchronized locally.Case 2. Fixq = 0.998, g=8, and various coupling strengths. It can be seen from Figure 12 that the spatiotemporal behaviours for D = 2 show clear indication of coherent and incoherent nodes. But when the coupling strength is D = 4 and 6, some of the systems have synchronization with their nearby nodes. However, by further increasing the coupling strength to D = 8, all of the systems are synchronized.Case 3. Fixq = 0.97, g=3, and various coupling strengths. Figure 13 shows that the systems are not synchronized for the coupling strength D = 1. With the coupling strength increasing gradually, the number of synchronized systems increases. When the coupling strength D = 6, all nodes are in a local synchronization state between different nodes.Figure 11
Spatiotemporal patterns of the ring network for different coupling strengths of Case 1. (a)D = 2. (b) D = 4. (c) D = 6. (d) D = 8.
(a)(b)(c)(d)Figure 12
Spatiotemporal patterns of the ring network for different coupling strengths of Case 2. (a)D = 2. (b) D = 4. (c) D = 6. (d) D = 8.
(a)(b)(c)(d)Figure 13
Spatiotemporal patterns of the neuronal ring network for different coupling strengths of Case 3. (a)D = 1. (b) D = 2. (c) D = 4. (d) D = 6.
(a)(b)(c)(d)As shown in Figures11–13, the behaviours of the network are different with different coupling strengths. When the coupling strength is small, the network does not show synchronization. For instance, it is shown in Figures 11(a) and 11(b) that the network is not synchronized. When the coupling strength is large, the network can be synchronized. For example, Figure 12(d) shows that there exists synchronization. The chimera state means coexistence of synchronization and asynchronization in the network. According to Figures 11–13, there are different kinds of chimera state in the network. The network with different fractional-order derivative order has different dynamics. As for the coupling strength for chimera states in the network, it should be proper. According to the analysis results, when the coupling strength D∈2,4, the chimera states are observed. Moreover, the derivative order also affects the dynamics of the network. In Case 1, we fix q = 0.998 and in Case 3 we let q = 0.97, and it is shown in Figures 11 and 13 that the dynamics of the network are different. It is easier for the network with q = 0.97 to synchronize or to show chimera states.


## 5. Conclusions

In this paper, a fractional-order memristor-based hyperchaotic system with no equilibrium points has been numerically analyzed by employing the predictor-corrector method. It shows that the system has rich dynamical behaviours, and different states including periodic circle, chaotic attractor, and hyperchaotic attractor are plotted in the system. The minimum derivative order of generating chaos is found, and it isq = 0.927. Meanwhile, through numerical simulation, the coexisting hidden attractors are observed with different initial conditions. Moreover, the SampEn complexity verifies the hidden dynamic characteristics of the system and it shows that this complexity analysis algorithm can be employed to detect the multiple coexisting attractors. Finally, by constructing the ring network structure of the system, it can be found that with increase of the coupling strength, the state of the network is different. The network can have the global synchronization state with large coupling strength, while in most cases, the network shows chimera states with proper coupling strength.



---

*Source: 2468134-2020-06-29.xml*