# Hybrid State Variable Incremental Integral for Reconstructing Extreme Multistability in Memristive Jerk System with Cubic Nonlinearity

Authors: Mo Chen; Yang Feng; Han Bao; Bocheng Bao; Huagan Wu; Quan Xu
Journal: Complexity
 (2019)

---

## Abstract

Memristive system with infinitely many equilibrium points has attracted much attention for the generation of extreme multistability, whose initial-dependent dynamics can be interpreted in a reduced-order model through incremental integral transformation of state variables. But, the memristive system with any extra nonlinear terms besides the memristor ones cannot be handled directly using this method. In addition, the transformed state variables could be divergent due to the asymmetry of the original system. To solve these problems, a hybrid state variable incremental integral (HSVII) method is proposed in this paper. With this method, the extreme multistability in a four-dimensional (4D) memristive jerk system with cubic nonlinearity is successfully reconstituted in a three-dimensional (3D) model and the divergent state variables are eliminated through ingenious linear state variable mapping. Thus, mechanism analysis and physical control of the special extreme multistability can readily be performed. A hardware circuit is finally designed and fabricated, and the theoretical and numerical results are verified by the experimental measurements. It is demonstrated that this HSVII method is effective for the analysis of multistable system with high-order nonlinearities.


---

## Body

## 1. Introduction

Initial-dependent multistability [1–5] phenomenon is an intrinsic phenomenon of nonlinear dynamical systems, which provides great flexibility or crisis [6–9] for seeking potential uses in chaos-based engineering applications. Dynamics of these multistable systems are significantly affected by the initial conditions, and efficient detection and control strategies [10–13] should be applied to keep them operating on the desired modes. Lots of researches have been performed to develop control techniques for multistable systems through changing the coexisting states and driving these systems from multistability to monostability [13–15]. However, the special extreme multistability phenomenon [16–25] still requires some novel control strategies, especially when we want to keep the flexibility of coexisting infinite many attractors and achieve the controllability in the meantime.Due to the special constitutive relation of ideal memristor, the memristor-based chaotic systems commonly exhibit the extreme multistability phenomenon with no [23] or line [17, 20, 21] or plane [16, 18] equilibrium points related to the initial conditions; thus they can hardly have long-term steady operation. Considering that the memristor is described as nondynamic element in the flux-charge domain [26], memristive circuit can be modeled in a reduced order form through taking flux and charge as main state variables [18, 27], which simplifies the analysis difficulties. Moreover, when the initial values are considered during the integration of voltage and current, the initial states of all dynamic elements can be expressed as standalone system parameters [28, 29] and the initial-dependent dynamics of the original memristive circuit can be reflected in the flux-charge domain. With this improvement, the incremental flux-charge analysis method can readily be applied for the investigations of relative complex memristive circuits [30], e.g., the memristive cellular neural networks [31] and arrays of memristor oscillators [32]. Additionally, it provides a new solution for the control of extreme multistability, which has been verified by the investigations of fourth-order one-memristor-based chaotic circuits [33–35]. It is revealed that, in the flux-charge domain, a reduced-order model is constructed, the original sensitive initial-dependent dynamics is reconstructed as the easily controlled system-parameters-related dynamics, and the dynamical effects of the initial values are readily interpreted [33, 35].For a memristive circuit consisting of linear discrete electronic components, such as resistor, inductor, capacitor, and operational-amplifier-based linear calculating module, the flux-charge model is easily formulated based on the duality of circuit variables [36]. Considering that the flux and charge are the time integrals of voltage and current, respectively, it can be found that the intrinsic quality of the flux-charge analysis method is the incremental integral transformation of the state variables. Inspired by this, the memristive system with normalized state variables can also be modeled in the reduced order form using incremental integral transformation [37]. However, since the integrations of high-order nonlinear terms are hard to be explicitly formulated, only memristive systems with no extra complex nonlinear terms besides the memristor term can be handled using this method. Moreover, when the original system is not symmetric to the origin point, the newly mapped state variables could be divergent and the newly formulated system may operate abnormally.To overcome the aforementioned two drawbacks, a hybrid state variable incremental integral (HSVII) method consisting of incremental integral transformation and linear transformation of the state variables is proposed in this paper. A memristive jerk system with a cubic nonlinearity is proposed as a typical example to demonstrate the effectiveness of the improved method. This 4D memristive jerk system possesses three-line equilibrium points and demonstrates an enormous number of coexisting attractors. Thus, it cannot be stably measured in the hardware circuit. Utilizing the proposed HSVII method, incremental integral transformation is first performed. The integral of the cubic nonlinear term is represented by a newly introduced state variable; thus a 4D system is obtained. Unfortunately, the newly mapped first and fourth state variables are slightly divergent due to the integration process. Thus, in the second step, the difference of the first and fourth state variables is chosen as a new state variable. Through this linear transformation, the divergence is eliminated and a 3D reduced-order model is finally obtained, which facilitates mechanism analyses and physical measurements of the original initial-dependent dynamics. Numerical simulations and experimental measurements prove that the random changing coexisting attractors of the original 4D memristive jerk system can be stably visualized in the 3D transformed system.The rest parts are organized as follows. In Section2, the memristive jerk system with cubic nonlinearity is introduced, whose extreme multistability phenomenon is numerically revealed by bifurcation diagrams, finite time Lyapunov exponent spectra, attraction basins, and phase portraits. In Section 3, based on the HSVII method, dimensionality reduction model of this memristive system is established. Formation mechanism of the initial-dependent dynamics is explored through evaluating the stabilities of the equilibrium points. Besides, dynamical correspondences between these two systems are demonstrated using two-parameter plots. In Section 4, hardware experiments for the verifications of the theoretical and numerical results are executed. The conclusions are drawn in the last section.


## 2. Memristive Jerk System with Three-Line Equilibrium Points

### 2.1. Mathematical Model

According to the construction method reported in [38], a 4D memristive jerk system is proposed by introducing an ideal memristor defined in [39] into a 3D chaotic jerk system described in [40]. The mathematical model is expressed as(1)x˙1=x2x˙2=1-x4x3x˙3=x1-ax2-x3-bx13x˙4=-x3From the physical circuit realization of view, the ideal memristor is introduced to replace the input resistor of the realization circuit for the second differential equation of the original 3D chaotic jerk system [38, 41]. The two control parametersa andb are fixed as 0.7 and 10, respectively, in the later sections. When the initial conditions are specified as x1(0)=0, x2(0)=10-6, x3(0)=0, and x4(0)=-0.5, this memristive jerk system generates a double-scroll chaotic attractor with Chua’s attractor structure, as shown in Figure 1.Figure 1
Double-scroll chaotic attractor generated by the memristive jerk system witha=0.7, b=10, x1(0)=0, x2(0)=10-6, x3(0)=0, and x4(0)=-0.5. (a) Phase portrait in thex1(t)−x2(t) plane. (b) Phase portrait in thex3(t)−x4(t) plane.
(a) (b)When the memristor initial conditionx4(0) is varied, the exhibited dynamical behaviors are greatly changed, which has been reported in various memristive circuits or systems [16, 17, 34, 38]. Takingx4(0) as a varying parameter, the bifurcation diagram ofx1 and the corresponding finite time Lyapunov exponent spectra are presented in Figure 2. For the investigations of memristor-initial-dependent dynamics, the other three initial conditions are set almost equal to 0. Whenx4(0) is gradually increased, system (1) starts from the periodic mode and goes into the chaotic mode via period doubling bifurcations and then breaks into the period-3 mode via tangent bifurcation. Thereafter, the dynamics evolves into chaotic mode again through period doubling bifurcations, and breaks into the period-3 mode via tangent bifurcation and then turns into chaos mode through chaos crisis, and finally jumps into divergent mode forx4(0) > 0.59. The minor nondivergent region of 0 <x4(0) < 0.59 is not presented in Figure 2. Within the three chaotic regions, there also exist some narrow periodic windows.Figure 2
The memristor-initial-dependent dynamics in the range of −30 <x4(0) < 0 witha = 0.7,b = 10,x1(0) = 0,x2(0) = 10−6, andx3(0) = 0. (a) Bifurcation diagram ofx1(t). (b) Finite time Lyapunov exponent spectra.
(a) (b)


### 2.2. Equilibrium Points and Stabilities

The equilibrium points of system (1) are given by(2)S+=1b,0,0,μS0=0,0,0,μS-=-1b,0,0,μwhere μ is any real constant. Obviously, system (1) has three-line equilibrium points parallel to thex4 axis. For b≥0, the equilibrium points characterized by (2) do exist. The Jacobian matrix at the equilibrium point S=(x-1,0,0,μ) is given by (3)JS=0100001-μ01-3bx-12-a-1000-10For S0, the characteristic equation of (3) is formulated as(4)λλ3+λ2-aμ-1λ+μ-1=0While forS±, the corresponding characteristic equation is described as(5)λλ3+λ2-aμ-1λ-2μ+2=0It can be seen that the locations of three-line equilibrium points are determined by system parameterb and memristor initial conditionμ, and their stabilities are determined by system parametera and memristor initial conditionμ.Whena = 0.7 andb = 10 are specified, the line equilibrium points described in (1) are expressed asS0 = (0, 0, 0,μ) andS± = (±0.32, 0, 0,μ). Their eigenvalues are calculated and depicted in Figure 3. The detailed eigenvalue distributions ofS0 andS± are summarized in Tables 1 and 2, respectively. With reference to these results, it can be found that, forμ < 1,S0 is an index-1 saddle-focus andS± are two index-2 saddle-foci; thus a double-scroll chaotic attractor could be evolved in the neighborhood ofS±. While for 1 <μ< 29.6,S0 becomes an index-2 saddle-focus andS± turn into two repulsive index-1 saddle-foci, and the trajectory quickly runs into infinite. Asμfurther increased,S0 change to be a saddle point with two positive real roots, whileS± still are two index-1 saddle-focus within 29.6 <μ< 60.5 and evolve into two saddle points with one positive real root for 60.5 <μ< 90. Consequently, when the memristor initial condition is set in the range of 1 <μ< 90, system (1) will remain operating in divergent mode.Table 1
Eigenvalues of line equilibrium pointS0 in different initial condition regions of μ.
μ λ1 λ2 λ3 λ4 μ < 1 Positive real root Zero Complex root with negative real part Complex root with negative real part 1 < μ < 29.6 Complex root with positive real part Complex root with positive real part Zero Negative real root μ > 29.6 Positive real root Positive real root Zero Negative real rootTable 2
Eigenvalues of line equilibrium pointS± in different initial condition regions of μ.
μ λ1 λ2 λ3 λ4 μ < 1 Complex root with positive real part Complex root with positive real part Zero Negative real root 1 < μ < 60.5 Positive real root Zero Complex root with negative real part Complex root with negative real part μ > 60.5 Positive real root Zero Negative real root Negative real rootFigure 3
Eigenvalues distributions of three-line equilibrium points with the variation ofμ, where the system parameters are set asa = 0.7 andb = 10. (a) Line equilibrium pointS0. (b) Line equilibrium pointsS±.
(a) (b)


### 2.3. Extreme Multistability with Complex Attraction Basin

Through numerical simulations, it is revealed that all the four initial conditions play vital roles in the formation of the dynamical behaviors. When different initial conditions are specified, various kinds of attractors are generated in system (1), leading to the emergence of extreme multistability. Two typical attraction basins in thex2(0)-x4(0) andx3(0)-x4(0) planes are presented in Figures 4(a) and 4(b), respectively. The regions marked with different colors represent different coexisting attractor types. To get a better visual effect, only partial attractors, which exhibit distinct topologies and possess relatively large attraction basins, are distinguished and exhibited in Figure 4.Figure 4
Extreme multistability demonstrated by attraction basins in two different planes. (a) Attraction basin in thex2(0)-x4(0) plane withx1(0) = 0 andx3(0) = 0.2. (b) Attraction basin in thex3(0)-x4(0) plane withx1(0) = 0 andx2(0) = 0.1.
(a) (b)Ten typical color regions alongx3(0) = 0 of the local attraction basins depicted in Figure 4(b) are summarized in Table 3. Phase portraits of the corresponding attractors are depicted in Figure 5, wherex4(0) is varied while the other three initial conditions are fixed asx1(0) = 0,x2(0) = 0.1, andx3(0) = 0. Compared with the initial conditions used in Figures 1 and 2, i.e.,x1(0) = 0,x2(0) = 10−6, andx3(0) = 0, this initial condition setting can rigorously achieved in the hardware circuit and is benefit for experimental verifications performed in the later section.Table 3
Attractor types in different attraction regions of Figure4.
Color Attractor type Phase portrait Color Attractor type Phase portrait Right-half period-4 limit cycle Figure5(a)  Integrated period-5 limit cycle Figure5(f) Right-spiral chaotic attractor Figure5(b)  Integrated period-9 limit cycle Figure5(g) Right-half period-3 limit cycle Figure5(c)  Double-scroll chaotic attractor Figure5(h) Right-spiral chaotic attractor Figure5(d)  Integrated period-3 limit cycle Figure5(i) Double-scroll chaotic attractor Figure5(e)  Right-half period-1 limit cycle Figure5(j)Figure 5
Phase portraits of typical coexisting attractors with different memristor initial conditions, wherex1(0) = 0,x2(0) = 0.1, andx3(0) = 0. (a)x4(0) = −22. (b)x4(0) = −19. (c)x4(0) = −16. (d)x4(0) = −13.5. (e)x4(0) = −10. (f)x4(0) = −8. (g)x4(0) = −5.8. (h)x4(0) = −1.6. (i)x4(0) = −1.4. (j)x4(0) = −1.05.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)According to Figure5, it is demonstrated that there are right-half and integrated limit cycles with different periodicities, as shown in Figures 5(a), 5(c), 5(f), 5(g), 5(i), and 5(j), as well as chaotic right-spiral and double-scroll attractors with different topologies, as depicted in Figures 5(b), 5(d), 5(e), and 5(h). Due to the symmetric line equilibrium pointsS+ andS−, symmetric left-half attractors could be found forx4(0) = −22, −19, −16, and −13.5 through settingx2(0) = −0.1.Moreover, whenx4(0) is settled down, the equilibrium points of system (1) are determined as three fixed points (0.32, 0, 0,x4(0)), (0, 0, 0,x4(0)), and (−0.32, 0, 0,x4(0)). Consequently, the generated attractors are centralized around different points along thex4 coordinate and disconnected in the phase space. Takingx1(0) = 0,x2(0) = 0.1,x3(0) = 0, andx4(0) = −21.1/−21/−20.9 as typical examples, three period-4 limit cycles with similar topologies but different locations are generated, as predicted in Figure 6(a). Whenx4(0) = −8.2/−8/−7.8 are selected, chaotic and periodic attractors with different topologies and positions are revealed, as illustrated in Figure 6(b). Thus, it can be concluded that the extreme multistability phenomenon, i.e., the coexistence of infinitely many disconnected attractors, is emerged in the 4D memristive jerk system.Figure 6
Disconnected phase portraits in thex1(t)-x4(t) plane with differentx4(0), wherex1(0) = 0,x2(0) = 0.1, andx3(0) = 0. (a) Period-4 limit cycles with similar topologies but different locations underx4(0) = −21.1/−21/−20.9. (b) Chaotic and periodic attractors with different topologies and locations underx4(0) = −8.2/−8/−7.8.
(a) (b)


## 2.1. Mathematical Model

According to the construction method reported in [38], a 4D memristive jerk system is proposed by introducing an ideal memristor defined in [39] into a 3D chaotic jerk system described in [40]. The mathematical model is expressed as(1)x˙1=x2x˙2=1-x4x3x˙3=x1-ax2-x3-bx13x˙4=-x3From the physical circuit realization of view, the ideal memristor is introduced to replace the input resistor of the realization circuit for the second differential equation of the original 3D chaotic jerk system [38, 41]. The two control parametersa andb are fixed as 0.7 and 10, respectively, in the later sections. When the initial conditions are specified as x1(0)=0, x2(0)=10-6, x3(0)=0, and x4(0)=-0.5, this memristive jerk system generates a double-scroll chaotic attractor with Chua’s attractor structure, as shown in Figure 1.Figure 1
Double-scroll chaotic attractor generated by the memristive jerk system witha=0.7, b=10, x1(0)=0, x2(0)=10-6, x3(0)=0, and x4(0)=-0.5. (a) Phase portrait in thex1(t)−x2(t) plane. (b) Phase portrait in thex3(t)−x4(t) plane.
(a) (b)When the memristor initial conditionx4(0) is varied, the exhibited dynamical behaviors are greatly changed, which has been reported in various memristive circuits or systems [16, 17, 34, 38]. Takingx4(0) as a varying parameter, the bifurcation diagram ofx1 and the corresponding finite time Lyapunov exponent spectra are presented in Figure 2. For the investigations of memristor-initial-dependent dynamics, the other three initial conditions are set almost equal to 0. Whenx4(0) is gradually increased, system (1) starts from the periodic mode and goes into the chaotic mode via period doubling bifurcations and then breaks into the period-3 mode via tangent bifurcation. Thereafter, the dynamics evolves into chaotic mode again through period doubling bifurcations, and breaks into the period-3 mode via tangent bifurcation and then turns into chaos mode through chaos crisis, and finally jumps into divergent mode forx4(0) > 0.59. The minor nondivergent region of 0 <x4(0) < 0.59 is not presented in Figure 2. Within the three chaotic regions, there also exist some narrow periodic windows.Figure 2
The memristor-initial-dependent dynamics in the range of −30 <x4(0) < 0 witha = 0.7,b = 10,x1(0) = 0,x2(0) = 10−6, andx3(0) = 0. (a) Bifurcation diagram ofx1(t). (b) Finite time Lyapunov exponent spectra.
(a) (b)


## 2.2. Equilibrium Points and Stabilities

The equilibrium points of system (1) are given by(2)S+=1b,0,0,μS0=0,0,0,μS-=-1b,0,0,μwhere μ is any real constant. Obviously, system (1) has three-line equilibrium points parallel to thex4 axis. For b≥0, the equilibrium points characterized by (2) do exist. The Jacobian matrix at the equilibrium point S=(x-1,0,0,μ) is given by (3)JS=0100001-μ01-3bx-12-a-1000-10For S0, the characteristic equation of (3) is formulated as(4)λλ3+λ2-aμ-1λ+μ-1=0While forS±, the corresponding characteristic equation is described as(5)λλ3+λ2-aμ-1λ-2μ+2=0It can be seen that the locations of three-line equilibrium points are determined by system parameterb and memristor initial conditionμ, and their stabilities are determined by system parametera and memristor initial conditionμ.Whena = 0.7 andb = 10 are specified, the line equilibrium points described in (1) are expressed asS0 = (0, 0, 0,μ) andS± = (±0.32, 0, 0,μ). Their eigenvalues are calculated and depicted in Figure 3. The detailed eigenvalue distributions ofS0 andS± are summarized in Tables 1 and 2, respectively. With reference to these results, it can be found that, forμ < 1,S0 is an index-1 saddle-focus andS± are two index-2 saddle-foci; thus a double-scroll chaotic attractor could be evolved in the neighborhood ofS±. While for 1 <μ< 29.6,S0 becomes an index-2 saddle-focus andS± turn into two repulsive index-1 saddle-foci, and the trajectory quickly runs into infinite. Asμfurther increased,S0 change to be a saddle point with two positive real roots, whileS± still are two index-1 saddle-focus within 29.6 <μ< 60.5 and evolve into two saddle points with one positive real root for 60.5 <μ< 90. Consequently, when the memristor initial condition is set in the range of 1 <μ< 90, system (1) will remain operating in divergent mode.Table 1
Eigenvalues of line equilibrium pointS0 in different initial condition regions of μ.
μ λ1 λ2 λ3 λ4 μ < 1 Positive real root Zero Complex root with negative real part Complex root with negative real part 1 < μ < 29.6 Complex root with positive real part Complex root with positive real part Zero Negative real root μ > 29.6 Positive real root Positive real root Zero Negative real rootTable 2
Eigenvalues of line equilibrium pointS± in different initial condition regions of μ.
μ λ1 λ2 λ3 λ4 μ < 1 Complex root with positive real part Complex root with positive real part Zero Negative real root 1 < μ < 60.5 Positive real root Zero Complex root with negative real part Complex root with negative real part μ > 60.5 Positive real root Zero Negative real root Negative real rootFigure 3
Eigenvalues distributions of three-line equilibrium points with the variation ofμ, where the system parameters are set asa = 0.7 andb = 10. (a) Line equilibrium pointS0. (b) Line equilibrium pointsS±.
(a) (b)


## 2.3. Extreme Multistability with Complex Attraction Basin

Through numerical simulations, it is revealed that all the four initial conditions play vital roles in the formation of the dynamical behaviors. When different initial conditions are specified, various kinds of attractors are generated in system (1), leading to the emergence of extreme multistability. Two typical attraction basins in thex2(0)-x4(0) andx3(0)-x4(0) planes are presented in Figures 4(a) and 4(b), respectively. The regions marked with different colors represent different coexisting attractor types. To get a better visual effect, only partial attractors, which exhibit distinct topologies and possess relatively large attraction basins, are distinguished and exhibited in Figure 4.Figure 4
Extreme multistability demonstrated by attraction basins in two different planes. (a) Attraction basin in thex2(0)-x4(0) plane withx1(0) = 0 andx3(0) = 0.2. (b) Attraction basin in thex3(0)-x4(0) plane withx1(0) = 0 andx2(0) = 0.1.
(a) (b)Ten typical color regions alongx3(0) = 0 of the local attraction basins depicted in Figure 4(b) are summarized in Table 3. Phase portraits of the corresponding attractors are depicted in Figure 5, wherex4(0) is varied while the other three initial conditions are fixed asx1(0) = 0,x2(0) = 0.1, andx3(0) = 0. Compared with the initial conditions used in Figures 1 and 2, i.e.,x1(0) = 0,x2(0) = 10−6, andx3(0) = 0, this initial condition setting can rigorously achieved in the hardware circuit and is benefit for experimental verifications performed in the later section.Table 3
Attractor types in different attraction regions of Figure4.
Color Attractor type Phase portrait Color Attractor type Phase portrait Right-half period-4 limit cycle Figure5(a)  Integrated period-5 limit cycle Figure5(f) Right-spiral chaotic attractor Figure5(b)  Integrated period-9 limit cycle Figure5(g) Right-half period-3 limit cycle Figure5(c)  Double-scroll chaotic attractor Figure5(h) Right-spiral chaotic attractor Figure5(d)  Integrated period-3 limit cycle Figure5(i) Double-scroll chaotic attractor Figure5(e)  Right-half period-1 limit cycle Figure5(j)Figure 5
Phase portraits of typical coexisting attractors with different memristor initial conditions, wherex1(0) = 0,x2(0) = 0.1, andx3(0) = 0. (a)x4(0) = −22. (b)x4(0) = −19. (c)x4(0) = −16. (d)x4(0) = −13.5. (e)x4(0) = −10. (f)x4(0) = −8. (g)x4(0) = −5.8. (h)x4(0) = −1.6. (i)x4(0) = −1.4. (j)x4(0) = −1.05.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)According to Figure5, it is demonstrated that there are right-half and integrated limit cycles with different periodicities, as shown in Figures 5(a), 5(c), 5(f), 5(g), 5(i), and 5(j), as well as chaotic right-spiral and double-scroll attractors with different topologies, as depicted in Figures 5(b), 5(d), 5(e), and 5(h). Due to the symmetric line equilibrium pointsS+ andS−, symmetric left-half attractors could be found forx4(0) = −22, −19, −16, and −13.5 through settingx2(0) = −0.1.Moreover, whenx4(0) is settled down, the equilibrium points of system (1) are determined as three fixed points (0.32, 0, 0,x4(0)), (0, 0, 0,x4(0)), and (−0.32, 0, 0,x4(0)). Consequently, the generated attractors are centralized around different points along thex4 coordinate and disconnected in the phase space. Takingx1(0) = 0,x2(0) = 0.1,x3(0) = 0, andx4(0) = −21.1/−21/−20.9 as typical examples, three period-4 limit cycles with similar topologies but different locations are generated, as predicted in Figure 6(a). Whenx4(0) = −8.2/−8/−7.8 are selected, chaotic and periodic attractors with different topologies and positions are revealed, as illustrated in Figure 6(b). Thus, it can be concluded that the extreme multistability phenomenon, i.e., the coexistence of infinitely many disconnected attractors, is emerged in the 4D memristive jerk system.Figure 6
Disconnected phase portraits in thex1(t)-x4(t) plane with differentx4(0), wherex1(0) = 0,x2(0) = 0.1, andx3(0) = 0. (a) Period-4 limit cycles with similar topologies but different locations underx4(0) = −21.1/−21/−20.9. (b) Chaotic and periodic attractors with different topologies and locations underx4(0) = −8.2/−8/−7.8.
(a) (b)


## 3. Dimensionality Reduction Modeling through HSVII Method

System (1) can easily be implemented using multipliers and op-amps linked with resistors and/or capacitors [4, 9]. However, the memristor initial condition and other initial conditions are hard to be rigorously assigned in the hardware circuit. Moreover, the large memristor initial condition easily pushes the op-amps into saturation and makes the implementation circuit to operate in unpredicted mode. As a result, the initial-dependent extreme multistability phenomena cannot be captured in the implementation circuit of system (1). To stably achieve a desired operation mode in the hardware circuit, the extreme multistability is reconstituted using the HSVII method consisting by incremental integral transformation and liner transformation of state variables.


### 3.1. State Variable Incremental Integral Transformation

Denoting(6)Xit=∫0txiτdτ,δi=xi0i=1,2,3,4 and integrating the four equations of (1) from 0 tot [37], there are(7)x1t-δ1=X2x2t-δ2=X3-∫0tx3τx4τdτx3t-δ3=X1-aX2-X3-b∫0tx13τdτx4t-δ4=-X3Based on the fourth equation of (1), we have dx4(t) = −x3(t)dt, and then the integration of the quadric term in the second equation of (7) can be expressed as(8)∫0tx3τx4τdτ=-∫0tx4τdx4τ=-0.5x42t+0.5x420According to (6) and the fourth equation of (7), (8) is formulated as(9)∫0tx3τx4τdτ=-0.5-X3+δ42+0.5δ42=-0.5X32t+δ4X3tDue to the special constitutive relation of memristor, the integration of memristor related high-order term can easily be expressed by the newly mapped state variables [33–35, 37]. The main difficulty lies in the integration of the cubic term in the third equation of (7), which cannot be explicitly formulized using the newly mapped state variables Xi(t).To solve this problem, a new state variable is defined as follows:(10)W=b∫0tx13dτSubstituting (9) and (10) into (7), a 4D system is constructed as follows:(11)X˙1=X2+δ1X˙2=0.5X32+1-δ4X3+δ2X˙3=X1-aX2-X3-W+δ3W˙=bX2+δ13According to the first equation of (6), we can determined that the initial conditions of (11) are fixed as zero, i.e.,X1(0) =X2(0) =X3(0) =W(0) = 0. In another word, system (11) should be initiated from the origin point (0, 0, 0, 0).However, when the initial-related system parameters are assigned according to the initial conditions used in Figure1, i.e.δ1 =δ3=0,δ2 = 10−6, andδ4 = −0.5, the trajectory of system (11) gradually spreads into infinite, as shown in Figure 7. It can be seen that the state variablesX2 andX3 exercise in the bounded regions, but the amplitudes ofX1 andW are slowly increased. This phenomenon is abnormal and not benefit for experimental measurements.Figure 7
Divergent attractors of system (11) witha = 0.7,b = 10,δ1 =δ3=0,δ2 = 10−6, andδ4 = −0.5. (a) Phase portrait in theX1-X2 plane. (b) Phase portrait in theX3-W plane.
(a) (b)


### 3.2. State Variable Linear Transformation

Considering thatX1 andW are all gradually increased, if we take their difference as a new state variable, the divergence can be eliminated. Inspired by this idea, three new state variables are introduced, which are defined as(12)Y1=X1-W,Y2=X2,Y3=X3With (12), system (11) is linearly transformed as(13)Y˙1=X˙1-W˙=Y2+δ1-bY2+δ13Y˙2=X˙2=0.5Y32+1-δ4Y3+δ2Y˙3=X˙3=X1-W-aY2-Y3+δ3and simplified as a 3D model, i.e.,(14)Y˙1=Y2+δ1-bY2+δ13Y˙2=0.5Y32+1-δ4Y3+δ2Y˙3=Y1-aY2-Y3+δ3The state variable correspondences of systems (14) and (1) are given as follows:(15)x1t=Y2t+δ1,x2t=Y˙2t,x3t=Y˙3t,x4t=-Y3t+δ4Since the initials of system (11) are fixed asX1(0) =X2(0) =X3(0) =W(0) = 0, the values ofY1(0),Y2(0), andY3(0) are all settled down to 0.When the initial-related system parameters are set as those adopted in Figure7, the typical phase portraits in theY1 –Y2 andY1 –Y3 planes are presented in Figures 8(a) and 8(b), respectively. Based on the relations described in (15), this chaotic attractor can be transformed back into its original state variable domain via the state variable derivation. The transformed phase portraits are plotted in Figures 9(a) and 9(b), respectively, which demonstrate the same topologies as those plotted in Figure 1.Figure 8
Typical chaotic phase portraits of system (14) witha = 0.7,b = 10,δ1 =δ3= 0,δ2 = 10−6, andδ4 = −0.5. (a) Phase portrait in theY1(t) –Y2(t) plane. (b) Phase portrait in theY1(t) –Y3(t) plane.
(a) (b)Figure 9
Transformed phase portraits of system (14) witha = 0.7,b = 10,δ1 =δ3= 0,δ2 = 10−6, andδ4 = −0.5. (a) Phase portrait in the Y2t+δ1-Y˙2(t) plane. (b) Phase portrait in the Y˙3t-[-Y3(t)+δ4] plane.
(a) (b)


### 3.3. Transformed Equilibrium Points and Stabilities

LettingY˙1=Y˙2=Y˙3=0, the equilibrium points of system (14) are given by(16)E1=-aδ1+ξ1-δ3,-δ1,ξ1,E2=-aδ1+ξ2-δ3,-δ1,ξ2,E3=aζ1+ξ1-δ3,ζ1,ξ1,E4=aζ1+ξ2-δ3,ζ1,ξ2,E5=aζ2+ξ1-δ3,ζ2,ξ1,E6=aζ2+ξ2-δ3,ζ2,ξ2where(17)ξ1=-1-δ4+1-δ42-2δ2,ξ2=-1-δ4-1-δ42-2δ2,ζ1=-δ1+1b,ζ2=-δ1-1bWhen(1-δ4)2≥2δ2 andb > 0, the equilibrium points characterized in (16) do exist, and the Jacobian matrix at the equilibrium point E=(Y-1,Y-2,Y-3) is given by (18)JE=01-3bY-2+δ12000Y-3+1-δ41-a-1ForE1 andE2, the characteristic equation is formulated as (19a)λ3+λ2+a1-δ42-2δ2λ-1-δ42-2δ2=0and(19b)λ3+λ2-a1-δ42-2δ2λ+1-δ42-2δ2=0respectively. Further, forE3 andE5, the characteristic equation is expressed as(19c)λ3+λ2+a1-δ42-2δ2λ+21-δ42-2δ2=0Finally,E4 andE6 are considered, and the characteristic equation is written as(19d)λ3+λ2-a1-δ42-2δ2λ-21-δ42-2δ2=0Whenδ4 is chosen as a changing variable and the other initial-related system parameters are settled asδ1 =δ3=0 andδ2 = 10−6 ≈ 0, locations and stabilities of the six equilibrium points are evaluated as follows. Forδ4 < 1, the equilibrium points defined in (16) are expressed asE1 = (0, 0, 0),E2 = (2δ4 − 2, 0, 2δ4 − 2),E3 = (0.22, 0.32, 0),E4 = (2δ4 − 1.78, 0.32, 2δ4 − 2),E5 = (−0.22, −0.32, 0), andE6 = (2δ4 − 2.22, −0.32, 2δ4 − 2), respectively. Their locations and stabilities are illustrated in Figure 10, in which the red, blue, and cyan lines represent the index-2 saddle-foci, index-1 saddle-foci, and saddle points, respectively. Obviously, the locations ofE1,E3, andE5 remain unchanged with the variation ofδ4, while theY1 andY3 coordinates ofE2,E4, andE6 are linearly proportional toδ4. In the region ofδ4 < 0,E2,E4, andE6 locates relatively far away from the initiating point (0, 0, 0); thus they do not affect the evolution of the moving trajectory. Chaotic or periodic attractors are generated in the neighborhood of the two index-2 saddle-fociE3 andE5. But for 0.6 <δ4 < 1, the repelling affections ofE4 andE6 tend to be preponderant, and the orbit is pushed into infinite.Figure 10
Locations and stabilities of the six determined equilibrium points, where the red lines indicate index-2 saddle-foci, the blue lines represent index-1 saddle-foci, and the cyan lines indicate saddle points. (a) Locations and stabilities of six equilibrium points for −28 <δ4< −27. (b) Locations and stabilities of equilibrium points for 0.5 <δ4< 1.5.
(a) (b)As forδ4 > 1, the six equilibriums are changed asE1 = (2δ4 − 2, 0, 2δ4 − 2),E2 = (0, 0, 0),E3 = (2δ4 − 1.78, 0.32, 2δ4 − 2),E4 = (0.22, 0.32, 0),E5 = (2δ4 − 2.22, −0.32, 2δ4 − 2), andE6 = (−0.22, −0.32, 0), respectively. With reference to Figure 10(b), it can be found thatE2is an index-2 saddle-focus located at the initiating point of system (14) andE4 andE6 are two symmetrical saddle points or index-1 saddle-foci situated at the two sides ofE2, whileE1,E3, andE5 locate relatively far away from the initiating point (0, 0, 0). Consequently, in this region, orbits initiating from the origin point quickly run into infinite.


### 3.4. Reconstruction of Extreme Multistability

Systems (1) and (14) describe the same nonlinear system in different state variable domains, and the extreme multistability of system (1) is reconstructed as the system-parameters-related dynamics of system (14). This dynamical correspondence is numerically revealed through the simulations of two-parameter plots.In Figure11, based on the largest Lyapunov exponent, the two-parameter plots of system (1) are depicted in thex2(0)-x4(0) andx3(0)-x4(0) initial spaces, where ODE45 Runge-Kutta algorithm with time-step 0.5 and time-end 1×104 is adopted. Correspondingly, the two-parameter plots of system (14) in theδ2-δ4 andδ3-δ4 parameter spaces are presented in Figure 12, in which the simulation parameters and labeled colors are kept the same as those used in Figure 11.Figure 11
Extreme multistability of the 4D system (1) demonstrated by two-parameter plots in the initial spaces. (a) Two-parameter plot in thex2(0)-x4(0) plane withx1(0) = 0 andx3(0) = 0.2. (b) Two-parameter plot in thex3(0)-x4(0) plane withx1(0) = 0 andx2(0) = 0.1.
(a) (b)Figure 12
Two-parameter plots of 3D system (14) in initial-related system parameter spaces. (a) Two-parameter plot in theδ2-δ4 plane withδ1 = 0 andδ3 = 0.2. (b) Two-parameter plot in theδ3-δ4 plane withδ1 = 0 andδ2 = 0.1.
(a) (b)With reference to Figures11 and 12, we can see that through the incremental integral transformation and linear transformation of state variables, dynamics in the initial spaces of the original 4D system (1) consists with that exhibited in the initial-related system parameter spaces of the reconstructed 3D system (14). It provides major advantages for handling extreme multistability in the transformed state variable domain [33, 35, 37].


## 3.1. State Variable Incremental Integral Transformation

Denoting(6)Xit=∫0txiτdτ,δi=xi0i=1,2,3,4 and integrating the four equations of (1) from 0 tot [37], there are(7)x1t-δ1=X2x2t-δ2=X3-∫0tx3τx4τdτx3t-δ3=X1-aX2-X3-b∫0tx13τdτx4t-δ4=-X3Based on the fourth equation of (1), we have dx4(t) = −x3(t)dt, and then the integration of the quadric term in the second equation of (7) can be expressed as(8)∫0tx3τx4τdτ=-∫0tx4τdx4τ=-0.5x42t+0.5x420According to (6) and the fourth equation of (7), (8) is formulated as(9)∫0tx3τx4τdτ=-0.5-X3+δ42+0.5δ42=-0.5X32t+δ4X3tDue to the special constitutive relation of memristor, the integration of memristor related high-order term can easily be expressed by the newly mapped state variables [33–35, 37]. The main difficulty lies in the integration of the cubic term in the third equation of (7), which cannot be explicitly formulized using the newly mapped state variables Xi(t).To solve this problem, a new state variable is defined as follows:(10)W=b∫0tx13dτSubstituting (9) and (10) into (7), a 4D system is constructed as follows:(11)X˙1=X2+δ1X˙2=0.5X32+1-δ4X3+δ2X˙3=X1-aX2-X3-W+δ3W˙=bX2+δ13According to the first equation of (6), we can determined that the initial conditions of (11) are fixed as zero, i.e.,X1(0) =X2(0) =X3(0) =W(0) = 0. In another word, system (11) should be initiated from the origin point (0, 0, 0, 0).However, when the initial-related system parameters are assigned according to the initial conditions used in Figure1, i.e.δ1 =δ3=0,δ2 = 10−6, andδ4 = −0.5, the trajectory of system (11) gradually spreads into infinite, as shown in Figure 7. It can be seen that the state variablesX2 andX3 exercise in the bounded regions, but the amplitudes ofX1 andW are slowly increased. This phenomenon is abnormal and not benefit for experimental measurements.Figure 7
Divergent attractors of system (11) witha = 0.7,b = 10,δ1 =δ3=0,δ2 = 10−6, andδ4 = −0.5. (a) Phase portrait in theX1-X2 plane. (b) Phase portrait in theX3-W plane.
(a) (b)


## 3.2. State Variable Linear Transformation

Considering thatX1 andW are all gradually increased, if we take their difference as a new state variable, the divergence can be eliminated. Inspired by this idea, three new state variables are introduced, which are defined as(12)Y1=X1-W,Y2=X2,Y3=X3With (12), system (11) is linearly transformed as(13)Y˙1=X˙1-W˙=Y2+δ1-bY2+δ13Y˙2=X˙2=0.5Y32+1-δ4Y3+δ2Y˙3=X˙3=X1-W-aY2-Y3+δ3and simplified as a 3D model, i.e.,(14)Y˙1=Y2+δ1-bY2+δ13Y˙2=0.5Y32+1-δ4Y3+δ2Y˙3=Y1-aY2-Y3+δ3The state variable correspondences of systems (14) and (1) are given as follows:(15)x1t=Y2t+δ1,x2t=Y˙2t,x3t=Y˙3t,x4t=-Y3t+δ4Since the initials of system (11) are fixed asX1(0) =X2(0) =X3(0) =W(0) = 0, the values ofY1(0),Y2(0), andY3(0) are all settled down to 0.When the initial-related system parameters are set as those adopted in Figure7, the typical phase portraits in theY1 –Y2 andY1 –Y3 planes are presented in Figures 8(a) and 8(b), respectively. Based on the relations described in (15), this chaotic attractor can be transformed back into its original state variable domain via the state variable derivation. The transformed phase portraits are plotted in Figures 9(a) and 9(b), respectively, which demonstrate the same topologies as those plotted in Figure 1.Figure 8
Typical chaotic phase portraits of system (14) witha = 0.7,b = 10,δ1 =δ3= 0,δ2 = 10−6, andδ4 = −0.5. (a) Phase portrait in theY1(t) –Y2(t) plane. (b) Phase portrait in theY1(t) –Y3(t) plane.
(a) (b)Figure 9
Transformed phase portraits of system (14) witha = 0.7,b = 10,δ1 =δ3= 0,δ2 = 10−6, andδ4 = −0.5. (a) Phase portrait in the Y2t+δ1-Y˙2(t) plane. (b) Phase portrait in the Y˙3t-[-Y3(t)+δ4] plane.
(a) (b)


## 3.3. Transformed Equilibrium Points and Stabilities

LettingY˙1=Y˙2=Y˙3=0, the equilibrium points of system (14) are given by(16)E1=-aδ1+ξ1-δ3,-δ1,ξ1,E2=-aδ1+ξ2-δ3,-δ1,ξ2,E3=aζ1+ξ1-δ3,ζ1,ξ1,E4=aζ1+ξ2-δ3,ζ1,ξ2,E5=aζ2+ξ1-δ3,ζ2,ξ1,E6=aζ2+ξ2-δ3,ζ2,ξ2where(17)ξ1=-1-δ4+1-δ42-2δ2,ξ2=-1-δ4-1-δ42-2δ2,ζ1=-δ1+1b,ζ2=-δ1-1bWhen(1-δ4)2≥2δ2 andb > 0, the equilibrium points characterized in (16) do exist, and the Jacobian matrix at the equilibrium point E=(Y-1,Y-2,Y-3) is given by (18)JE=01-3bY-2+δ12000Y-3+1-δ41-a-1ForE1 andE2, the characteristic equation is formulated as (19a)λ3+λ2+a1-δ42-2δ2λ-1-δ42-2δ2=0and(19b)λ3+λ2-a1-δ42-2δ2λ+1-δ42-2δ2=0respectively. Further, forE3 andE5, the characteristic equation is expressed as(19c)λ3+λ2+a1-δ42-2δ2λ+21-δ42-2δ2=0Finally,E4 andE6 are considered, and the characteristic equation is written as(19d)λ3+λ2-a1-δ42-2δ2λ-21-δ42-2δ2=0Whenδ4 is chosen as a changing variable and the other initial-related system parameters are settled asδ1 =δ3=0 andδ2 = 10−6 ≈ 0, locations and stabilities of the six equilibrium points are evaluated as follows. Forδ4 < 1, the equilibrium points defined in (16) are expressed asE1 = (0, 0, 0),E2 = (2δ4 − 2, 0, 2δ4 − 2),E3 = (0.22, 0.32, 0),E4 = (2δ4 − 1.78, 0.32, 2δ4 − 2),E5 = (−0.22, −0.32, 0), andE6 = (2δ4 − 2.22, −0.32, 2δ4 − 2), respectively. Their locations and stabilities are illustrated in Figure 10, in which the red, blue, and cyan lines represent the index-2 saddle-foci, index-1 saddle-foci, and saddle points, respectively. Obviously, the locations ofE1,E3, andE5 remain unchanged with the variation ofδ4, while theY1 andY3 coordinates ofE2,E4, andE6 are linearly proportional toδ4. In the region ofδ4 < 0,E2,E4, andE6 locates relatively far away from the initiating point (0, 0, 0); thus they do not affect the evolution of the moving trajectory. Chaotic or periodic attractors are generated in the neighborhood of the two index-2 saddle-fociE3 andE5. But for 0.6 <δ4 < 1, the repelling affections ofE4 andE6 tend to be preponderant, and the orbit is pushed into infinite.Figure 10
Locations and stabilities of the six determined equilibrium points, where the red lines indicate index-2 saddle-foci, the blue lines represent index-1 saddle-foci, and the cyan lines indicate saddle points. (a) Locations and stabilities of six equilibrium points for −28 <δ4< −27. (b) Locations and stabilities of equilibrium points for 0.5 <δ4< 1.5.
(a) (b)As forδ4 > 1, the six equilibriums are changed asE1 = (2δ4 − 2, 0, 2δ4 − 2),E2 = (0, 0, 0),E3 = (2δ4 − 1.78, 0.32, 2δ4 − 2),E4 = (0.22, 0.32, 0),E5 = (2δ4 − 2.22, −0.32, 2δ4 − 2), andE6 = (−0.22, −0.32, 0), respectively. With reference to Figure 10(b), it can be found thatE2is an index-2 saddle-focus located at the initiating point of system (14) andE4 andE6 are two symmetrical saddle points or index-1 saddle-foci situated at the two sides ofE2, whileE1,E3, andE5 locate relatively far away from the initiating point (0, 0, 0). Consequently, in this region, orbits initiating from the origin point quickly run into infinite.


## 3.4. Reconstruction of Extreme Multistability

Systems (1) and (14) describe the same nonlinear system in different state variable domains, and the extreme multistability of system (1) is reconstructed as the system-parameters-related dynamics of system (14). This dynamical correspondence is numerically revealed through the simulations of two-parameter plots.In Figure11, based on the largest Lyapunov exponent, the two-parameter plots of system (1) are depicted in thex2(0)-x4(0) andx3(0)-x4(0) initial spaces, where ODE45 Runge-Kutta algorithm with time-step 0.5 and time-end 1×104 is adopted. Correspondingly, the two-parameter plots of system (14) in theδ2-δ4 andδ3-δ4 parameter spaces are presented in Figure 12, in which the simulation parameters and labeled colors are kept the same as those used in Figure 11.Figure 11
Extreme multistability of the 4D system (1) demonstrated by two-parameter plots in the initial spaces. (a) Two-parameter plot in thex2(0)-x4(0) plane withx1(0) = 0 andx3(0) = 0.2. (b) Two-parameter plot in thex3(0)-x4(0) plane withx1(0) = 0 andx2(0) = 0.1.
(a) (b)Figure 12
Two-parameter plots of 3D system (14) in initial-related system parameter spaces. (a) Two-parameter plot in theδ2-δ4 plane withδ1 = 0 andδ3 = 0.2. (b) Two-parameter plot in theδ3-δ4 plane withδ1 = 0 andδ2 = 0.1.
(a) (b)With reference to Figures11 and 12, we can see that through the incremental integral transformation and linear transformation of state variables, dynamics in the initial spaces of the original 4D system (1) consists with that exhibited in the initial-related system parameter spaces of the reconstructed 3D system (14). It provides major advantages for handling extreme multistability in the transformed state variable domain [33, 35, 37].


## 4. Circuit Synthesis and Experimental Measurements

### 4.1. Circuit Synthesis

To get high-quality experimental results, the state variables of system (14) are enlarged 4 times through taking the following transformations:(20)Z1=4Y1,Z2=4Y2,Z3=4Y3,Thus, system (14) is rewritten as (21)Z˙1=Z2+4δ1-b16Z2+4δ13Z˙2=18Z32+1-δ4Z3+4δ2Z˙3=Z1-aZ2-Z3+4δ3The physical implementation circuit of system (21) is synthesized as shown in Figure 13. The circuit equations are formulated as(22)RCdv1dt=v2+4δ1-RR1v2+4δ13RCdv2dt=RR2v32+RR3v3+4δ2RCdv3dt=v1-RR4v2-v3+4δ3Three capacitor voltages v1, v2, and v3 indicate the three state variables of system (21) andRC is the time constant. Denote the time constant asRC = 1 ms, i.e.,R = 10 kΩ andC = 100 nF. The element parameters of Figure 13 are determined as(23)R1=16Rb=16kΩ,R2=8R=80kΩ,R3=R1-δ4,R4=Ra=14.3kΩ.Figure 13
Schematic diagram of the implementation circuit for the third-order reconstructed memristive system.Based on this circuit, the desired operation mode of system (1) can be generated through tuningR3 and the three extra DC control signalsV1 = −4δ1 V,V2 = 4δ2 V, andV3 = 4δ3 V. The zero initial states of three capacitors are achieved through shortening them out for a while before powering on the circuit. The topologies of the generated attractors are reformed using differential circuits with proper time constants.


### 4.2. Breadboard Experiments

The breadboard prototype circuit is made using discrete components including operational amplifiers AD711KN and multiplier AD633JNZ with ±15 V DC bias voltages, ceramic capacitors, metal film resistors, and precision potentiometers. During measurements, a 5 MΩ parallel resistor is added to the capacitorC of the ideal integrators implemented byU1 andU3 to eliminate the DC voltage integral drift [42, 43]. Due to the large resistance value, the added parallel resistors do not affect the operation state of the implementation circuit. Thus, the desired phase portraits can be stably measured by a digital oscilloscope in X-Y mode.Through tuningR3 and the three extra DC control signals, any coexisting operation state of system (1) can be reproduced in the constructed hardware circuit. Five typical operation modes of Figures 5(a)–5(c), 5(e), and 5(f) are selected for the experimental verifications. Correspondingly, in system (14), the initial-related system parameters are specified asδ1 = 0,δ2 = 0.1, andδ3 = 0, while in the hardware circuit, the three DC control signals are fixed asV1 = 0,V2 = 0.4 V, andV3 = 0.Whenδ4 is set to −22, −19, −16, −10, and −8, respectively, the phase portraits of system (14) and their transformed phase portraits are numerically simulated by MATLAB, as shown in Figures 14(a)–14(e), respectively. Accordingly, when synchronously turning the values of the precision adjustable resistorsR3 as 475 Ω, 531 Ω, 614 Ω, 1.06 kΩ, and 1.27 kΩ, respectively, the phase portraits are experimentally measured and presented in Figures 15(a)–15(e), respectively. The waveform of v˙2(t) is obtained through using a differential circuit. The time scale of the differential circuit is specified as 10 kΩ × 100 nF = 1 ms. Thus, we can observe the phase portraits with exactly the same topological structures as those generated in the original 4D memristive system. These results are both consistent with each other. It is demonstrated that the initial-dependent coexisting attractors can be stably visualized in the reconstructed system, reflecting the efficiency of the HSVII method. In addition, the experimental values of resistorsR3 are slightly different from their theoretical values due to the existence of the parasitic circuit parameters and measurement errors.Figure 14
Phase portraits numerically simulated by MATLAB for different values ofδ4, in whichδ1 = 0,δ2 = 0.1, andδ3 = 0. (a)δ4 = −22. (b)δ4 = −19. (c)δ4 = −16. (d)δ4 = −10. (e)δ4 = −8.Figure 15
Phase portraits experimentally measured by breadboard circuit for different values of precision adjustable resistorR3. (a)R3 = 475 Ω. (b)R3 = 531 Ω. (c)R3 = 614 Ω. (d)R3 = 1.06 kΩ. (e)R3 = 1.27 kΩ.


## 4.1. Circuit Synthesis

To get high-quality experimental results, the state variables of system (14) are enlarged 4 times through taking the following transformations:(20)Z1=4Y1,Z2=4Y2,Z3=4Y3,Thus, system (14) is rewritten as (21)Z˙1=Z2+4δ1-b16Z2+4δ13Z˙2=18Z32+1-δ4Z3+4δ2Z˙3=Z1-aZ2-Z3+4δ3The physical implementation circuit of system (21) is synthesized as shown in Figure 13. The circuit equations are formulated as(22)RCdv1dt=v2+4δ1-RR1v2+4δ13RCdv2dt=RR2v32+RR3v3+4δ2RCdv3dt=v1-RR4v2-v3+4δ3Three capacitor voltages v1, v2, and v3 indicate the three state variables of system (21) andRC is the time constant. Denote the time constant asRC = 1 ms, i.e.,R = 10 kΩ andC = 100 nF. The element parameters of Figure 13 are determined as(23)R1=16Rb=16kΩ,R2=8R=80kΩ,R3=R1-δ4,R4=Ra=14.3kΩ.Figure 13
Schematic diagram of the implementation circuit for the third-order reconstructed memristive system.Based on this circuit, the desired operation mode of system (1) can be generated through tuningR3 and the three extra DC control signalsV1 = −4δ1 V,V2 = 4δ2 V, andV3 = 4δ3 V. The zero initial states of three capacitors are achieved through shortening them out for a while before powering on the circuit. The topologies of the generated attractors are reformed using differential circuits with proper time constants.


## 4.2. Breadboard Experiments

The breadboard prototype circuit is made using discrete components including operational amplifiers AD711KN and multiplier AD633JNZ with ±15 V DC bias voltages, ceramic capacitors, metal film resistors, and precision potentiometers. During measurements, a 5 MΩ parallel resistor is added to the capacitorC of the ideal integrators implemented byU1 andU3 to eliminate the DC voltage integral drift [42, 43]. Due to the large resistance value, the added parallel resistors do not affect the operation state of the implementation circuit. Thus, the desired phase portraits can be stably measured by a digital oscilloscope in X-Y mode.Through tuningR3 and the three extra DC control signals, any coexisting operation state of system (1) can be reproduced in the constructed hardware circuit. Five typical operation modes of Figures 5(a)–5(c), 5(e), and 5(f) are selected for the experimental verifications. Correspondingly, in system (14), the initial-related system parameters are specified asδ1 = 0,δ2 = 0.1, andδ3 = 0, while in the hardware circuit, the three DC control signals are fixed asV1 = 0,V2 = 0.4 V, andV3 = 0.Whenδ4 is set to −22, −19, −16, −10, and −8, respectively, the phase portraits of system (14) and their transformed phase portraits are numerically simulated by MATLAB, as shown in Figures 14(a)–14(e), respectively. Accordingly, when synchronously turning the values of the precision adjustable resistorsR3 as 475 Ω, 531 Ω, 614 Ω, 1.06 kΩ, and 1.27 kΩ, respectively, the phase portraits are experimentally measured and presented in Figures 15(a)–15(e), respectively. The waveform of v˙2(t) is obtained through using a differential circuit. The time scale of the differential circuit is specified as 10 kΩ × 100 nF = 1 ms. Thus, we can observe the phase portraits with exactly the same topological structures as those generated in the original 4D memristive system. These results are both consistent with each other. It is demonstrated that the initial-dependent coexisting attractors can be stably visualized in the reconstructed system, reflecting the efficiency of the HSVII method. In addition, the experimental values of resistorsR3 are slightly different from their theoretical values due to the existence of the parasitic circuit parameters and measurement errors.Figure 14
Phase portraits numerically simulated by MATLAB for different values ofδ4, in whichδ1 = 0,δ2 = 0.1, andδ3 = 0. (a)δ4 = −22. (b)δ4 = −19. (c)δ4 = −16. (d)δ4 = −10. (e)δ4 = −8.Figure 15
Phase portraits experimentally measured by breadboard circuit for different values of precision adjustable resistorR3. (a)R3 = 475 Ω. (b)R3 = 531 Ω. (c)R3 = 614 Ω. (d)R3 = 1.06 kΩ. (e)R3 = 1.27 kΩ.


## 5. Conclusions

In this paper, a 4D memristive jerk system is present, which is obtained by introducing an ideal memristor into a 3D chaotic jerk system. This system possesses three-line equilibrium points restricted by the system parameterb and memristor initial conditionx4(0). Consequently, when the system parametersa andb are kept unchanged, coexisting attractors with different topologies or locations are generated with the variations of any initial conditions, leading to the generation of extreme multistability.In order to perform inner interpretation and stable control of the sensitive extreme multistability phenomenon, this memristive jerk system is remodeled using a novel HSVII method consisting of incremental integral transformation and linear transformation of state variables. With this method, the formation of incremental integral of cubic nonlinear term and the divergent of the newly mapped state variables can be resolved. A dimensionality reduction system, which maintaining the dynamics of the original 4D system, is successfully constructed. Based on this reduced-order model, extreme multistability with an enormous number of coexisting attractors is interpreted through evaluating the mapped determined equilibrium points, depicting the two-parameter plots, and simulating the phase portraits. A hardware implementation circuit is successfully designed, in which the sensitive initial conditions of the original memristive system are reflected by a resistor and three extra DC control signals. Consequently, the initial-dependent dynamics is stably verified. These results demonstrate that this HSVII method is effective for the handling of memristive system with high-order nonlinear terms. It lays a good foundation for future investigations of ordinary nonlinear systems [5, 21, 22] or even complex memristive neural systems or circuits [44, 45] with multistability or extreme multistability.



---

*Source: 8549472-2019-06-23.xml*