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.

1. Introduction

Initial-dependent multistability [15] phenomenon is an intrinsic phenomenon of nonlinear dynamical systems, which provides great flexibility or crisis [69] 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 [1013] 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 [1315]. However, the special extreme multistability phenomenon [1625] 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 [3335]. 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 Section 2, 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 asFrom 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 parameters a and b are fixed as 0.7 and 10, respectively, in the later sections. When the initial conditions are specified as , , , and , this memristive jerk system generates a double-scroll chaotic attractor with Chua’s attractor structure, as shown in Figure 1.

When the memristor initial condition x4(0) is varied, the exhibited dynamical behaviors are greatly changed, which has been reported in various memristive circuits or systems [16, 17, 34, 38]. Taking x4(0) as a varying parameter, the bifurcation diagram of x1 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. When x4(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 for x4(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.

2.2. Equilibrium Points and Stabilities

The equilibrium points of system (1) are given bywhere is any real constant. Obviously, system (1) has three-line equilibrium points parallel to the x4 axis. For , the equilibrium points characterized by (2) do exist. The Jacobian matrix at the equilibrium point is given by For , the characteristic equation of (3) is formulated asWhile for S±, the corresponding characteristic equation is described asIt can be seen that the locations of three-line equilibrium points are determined by system parameter b and memristor initial condition μ, and their stabilities are determined by system parameter a and memristor initial condition μ.

When a = 0.7 and b = 10 are specified, the line equilibrium points described in (1) are expressed as S0 = (0, 0, 0, μ) and S± = (±0.32, 0, 0, μ). Their eigenvalues are calculated and depicted in Figure 3. The detailed eigenvalue distributions of S0 and S± 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 and S± are two index-2 saddle-foci; thus a double-scroll chaotic attractor could be evolved in the neighborhood of S±. While for 1 < μ < 29.6, S0 becomes an index-2 saddle-focus and S± 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, while S± 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.

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 the x2(0)-x4(0) and x3(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.

Ten typical color regions along x3(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, where x4(0) is varied while the other three initial conditions are fixed as x1(0) = 0, x2(0) = 0.1, and x3(0) = 0. Compared with the initial conditions used in Figures 1 and 2, i.e., x1(0) = 0, x2(0) = 10−6, and x3(0) = 0, this initial condition setting can rigorously achieved in the hardware circuit and is benefit for experimental verifications performed in the later section.

According to Figure 5, 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 points S+ and S, symmetric left-half attractors could be found for x4(0) = −22, −19, −16, and −13.5 through setting x2(0) = −0.1.

Moreover, when x4(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 the x4 coordinate and disconnected in the phase space. Taking x1(0) = 0, x2(0) = 0.1, x3(0) = 0, and x4(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). When x4(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.

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 and integrating the four equations of (1) from 0 to t [37], there areBased 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

According to (6) and the fourth equation of (7), (8) is formulated asDue 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 [3335, 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 .

To solve this problem, a new state variable is defined as follows:Substituting (9) and (10) into (7), a 4D system is constructed as follows:According 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 Figure 1, 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 variables X2 and X3 exercise in the bounded regions, but the amplitudes of X1 and W are slowly increased. This phenomenon is abnormal and not benefit for experimental measurements.

3.2. State Variable Linear Transformation

Considering that X1 and W 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 asWith (12), system (11) is linearly transformed asand simplified as a 3D model, i.e.,

The state variable correspondences of systems (14) and (1) are given as follows:Since the initials of system (11) are fixed as X1(0) = X2(0) = X3(0) = W(0) = 0, the values of Y1(0), Y2(0), and Y3(0) are all settled down to 0.

When the initial-related system parameters are set as those adopted in Figure 7, the typical phase portraits in the Y1 Y2 and Y1 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.

3.3. Transformed Equilibrium Points and Stabilities

Letting , the equilibrium points of system (14) are given bywhere

When and b > 0, the equilibrium points characterized in (16) do exist, and the Jacobian matrix at the equilibrium point is given by For E1 and E2, the characteristic equation is formulated as andrespectively. Further, for E3 and E5, the characteristic equation is expressed asFinally, E4 and E6 are considered, and the characteristic equation is written as

When δ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 as E1 = (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), and E6 = (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 of E1, E3, and E5 remain unchanged with the variation of δ4, while the Y1 and Y3 coordinates of E2, E4, and E6 are linearly proportional to δ4. In the region of δ4 < 0, E2, E4, and E6 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-foci E3 and E5. But for 0.6 < δ4 < 1, the repelling affections of E4 and E6 tend to be preponderant, and the orbit is pushed into infinite.

As for δ4 > 1, the six equilibriums are changed as E1 = (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), and E6 = (−0.22, −0.32, 0), respectively. With reference to Figure 10(b), it can be found that E2 is an index-2 saddle-focus located at the initiating point of system (14) and E4 and E6 are two symmetrical saddle points or index-1 saddle-foci situated at the two sides of E2, while E1, E3, and E5 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 Figure 11, based on the largest Lyapunov exponent, the two-parameter plots of system (1) are depicted in the x2(0)-x4(0) and x3(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.

With reference to Figures 11 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:Thus, system (14) is rewritten as

The physical implementation circuit of system (21) is synthesized as shown in Figure 13. The circuit equations are formulated asThree capacitor voltages , , and indicate the three state variables of system (21) and RC is the time constant. Denote the time constant as RC = 1 ms, i.e., R = 10 kΩ and C = 100 nF. The element parameters of Figure 13 are determined as

Based on this circuit, the desired operation mode of system (1) can be generated through tuning R3 and the three extra DC control signals V1 = −4δ1 V, V2 = 4δ2 V, and V3 = 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 capacitor C of the ideal integrators implemented by U1 and U3 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 tuning R3 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 as V1 = 0, V2 = 0.4 V, and V3 = 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 resistors R3 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 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 resistors R3 are slightly different from their theoretical values due to the existence of the parasitic circuit parameters and measurement errors.

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 parameter b and memristor initial condition x4(0). Consequently, when the system parameters a and b 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.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundations of China under Grant Nos. 61601062, 51777016, 51607013, and 61801054.