Initial-Condition Effects on a Two-Memristor-Based Jerk System

: Memristor-based systems can exhibit the phenomenon of extreme multi-stability, which results in the coexistence of inﬁnitely many attractors. However, most of the recently published literature focuses on the extreme multi-stability related to memristor initial conditions rather than non-memristor initial conditions. In this paper, we present a new ﬁve-dimensional (5-D) two-memristor-based jerk (TMJ) system and study complex dynamical effects induced by memristor and non-memristor initial conditions therein. Using multiple numerical methods, coupling-coefﬁcient-reliant dynamical behaviors under different memristor initial conditions are disclosed, and the dynamical effects of the memristor initial conditions under different non-memristor initial conditions are revealed. The numerical results show that the dynamical behaviors of the 5-D TMJ system are not only dependent on the coupling coefﬁcients, but also dependent on the memristor and non-memristor initial conditions. In addition, with the analog and digital implementations of the 5-D TMJ system, PSIM circuit simulations and microcontroller-based hardware experiments validate the numerical results.

In recent years, memristor-based chaotic circuits and systems have been broadly investigated, since memristor is a special nonlinear circuit component with memory effect and synaptic plasticity [19,28,29]. This particular type of chaotic circuit and system is conductive to deriving coexisting infinitely many attractors. Such coexistence of infinitely many attractors means initial-condition-related extreme multi-stability [30][31][32]. To implement the initial-condition-related extreme multi-stability, an effective and simple method is to introduce memristor into an existing dynamical system to construct a new memristor-based dynamical system, which is different from the method of using periodic trigonometric function to realize the initial condition-boosted infinitely many attractors in some special boostable systems [33,34]. In fact, the memristor-based dynamical system has a total bifurcation route to chaos with the evolution of the initial conditions [19] and can display the coexistence of various types of attractors [27]. However, most of the recently published literature only focuses on the extreme multi-stability related to the initial conditions of memristors [35][36][37][38], and little on the extreme multi-stability related to the initial conditions of non-memristors. In this paper, we present a new 5-D TMJ system and emphatically study complex dynamical effects induced by the initial conditions of memristors and non-memristors therein. Thus, the dynamical effects of the initial conditions on the 5-D TMJ system are disclosed comprehensively, which have not been wholly reported in the previously published literature.
The rest of this paper is arranged as follows. Section 2 first presents a 5-D TMJ system and then studies complex dynamics related to the coupling coefficients. Section 3 focuses on the dynamical effects of memristor and non-memristor initial conditions on the TMJ system. With the analog and digital implementations, PSIM circuit simulations and experimental measurements are carried out to validate the numerical simulations in Section 4. Finally, the paper is summarized in Section 5.

The TMJ System with Complex Dynamics
This section first presents a new 5-D TMJ system and then studies complex dynamics related to the coupling coefficients using multiple numerical methods.

Mathematical Modeling
A jerk system is a third-order ordinary differential equation. It is of the following form: J(·) is a nonlinear function. Since the mathematical model in (1) represents the cubic time derivative of variable x, it is named the 'jerk'. A simple jerk system used in this paper is improved from [39] and its general jerk form can be described as follows: where a and b are two positive control parameters. Denote x = x 1 ,ẋ = x 2 , andẍ = ax 3 , respectively. The simple jerk system in (2) can be rewritten as: where x 1 , x 2 , and x 3 are three state variables. Therefore, the simple jerk system is a threedimensional nonlinear dynamical system. The simple jerk system in (3) only has a cubic polynomial nonlinearity of state variable x 1 . On this basis, and referring to [35], a new 5-D TMJ system is presented in this paper, which is achieved by replacing the cubic polynomial in brackets with one memristor and introducing another memristor into the first equation.
For an input x and an output y, a flux-controlled memristor with state variable ϕ can be described as: The memductance W(ϕ) in (4) is selected as a threshold nonlinear function bounded above and below, which can be expressed by: Thus, the 5-D TMJ system is established using the above two memristors as: where x 4 and x 5 represent the inner state variables of two memristors, and k 1 and k 2 represent two changeable coupling coefficients. Therefore, the system modeled by (6) is a 5-D nonlinear dynamical system. Besides, the initial conditions are defined as ICs = (x 1 (0), x 2 (0), x 3 (0), x 4 (0), x 5 (0)), in which the first three initial conditions IC 1 = (x 1 (0), x 2 (0), x 3 (0)) are called the non-memristor initial conditions and the last two initial conditions IC 2 = (x 4 (0), x 5 (0)) are called the memristor initial conditions. To aim at the revelation of extreme multi-stability, two positive control parameters in (6) are kept unchanged as a = 2.5 and b = 0.8. When the initial conditions are fixed as ICs = (10 −9 , 0, 0, 1, 0) and the two changeable coupling coefficients are determined as k 1 = 1 and k 2 = 2.3, the 5-D TMJ system can generate a representative chaotic attractor, whose phase portraits in two different planes are depicted in Figure 1. Here, a MATLAB ODE45 algorithm with time-step 0.01 and time-span (400, 1000) is utilized. The results show that chaotic dynamics can be established in the 5-D TMJ system.

Coupling Coefficient-Reliant Complex Dynamics
Bifurcation diagram and Lyapunov exponent (LE) spectra of a dynamical system are taken as two dynamical indicators to characterize the type of bifurcation scenario giving rise to chaos. To show complex dynamical behaviors in system (6) intuitively, the twoparameter bifurcation diagram and dynamical map are numerically simulated in the k 1 -k 2 plane, as shown in Figure 2, where the initial conditions are fixed as ICs = (10 −9 , 0, 0, 1, 0). The two-parameter bifurcation diagram is plotted by measuring the periodicities of the state variable x 1 based on a MATLAB ODE45 algorithm with time-step 0.01 and time-span (400, 1000) and the two-parameter dynamical map is depicted by calculating the maximal Lyapunov exponent (LE) based on Wolf's Jacobi method.
In Figure 2a, the red block labeled CH represents chaos, the black block labeled UB represents unbounded behavior, the white block labeled P0 represents stable point, the yellow block labeled MP represents multi-period with the periodicity greater than 8, and the other color blocks labeled P1 to P8 represent period-1 to period-8. Meanwhile, in Figure 2b, the pink-red-yellow-cyan blocks with positive maximal LE denote chaotic behaviors and the blue blocks with zero maximal LE denote stable-point behaviors with different positions or periodic behaviors with different periodicities. Besides, the dynamical distributions described by Figure 2a,b are consistent with each other, which clarifies how the dynamical behaviors evolve in the k 1 -k 2 plane. As can be seen from Figure 2, chaotic behaviors with several periodic windows appear in the center left part of the pictures, implying the existence of tangent bifurcations and chaos crises; the dynamical behaviors have the transitions from P1 to P2, to P4, and to P8, demonstrating the occurrence of period-doubling bifurcations. To demonstrate the dynamical effects of the initial conditions, the initial conditions of the 5-D TMJ system are considered as ICs = (1, 0, 0, 1, 0) and (0, 1, 0, 1, 0). In other words, the memristor initial conditions are kept as IC 2 = (1, 0), while the non-memristor initial conditions are set to IC 1 = (1, 0, 0) and (0, 1, 0), respectively. For these two different cases of the initial conditions, the dynamical distributions of the 5-D TMJ system vary greatly in the parameter intervals concerned. Under these two sets of initial conditions, the two-parameter bifurcation diagrams are numerically simulated in the k 1 -k 2 plane, as shown in Figure 3. Looking at the dynamical distributions in Figures 2a and 3, it is easy to see that they are quite different from each other. This shows that the initial conditions of different state variables have great influence on the dynamical behaviors of the 5-D TMJ system. To describe the evolution of its dynamical behaviors with respect to the singleparameter interval, the single-parameter bifurcation diagram and LE spectra are employed. The initial conditions are set as ICs = (10 −9 , 0, 0, 1, 0) for simplicity.
First, the coupling coefficient k 2 is set as a fixed value of 2.3, and the coupling coefficient k 1 varies within the interval (0.2, 2). Second, k 1 = 1, and k 2 varies within the interval (1.8, 2.6). The single-parameter bifurcation diagrams and the corresponding LE spectra are simulated by MATLAB software and their results are depicted in the lower half and the upper half of Figure 4, respectively. Note that the single-parameter bifurcation diagrams in the lower half of Figure 4 are acquired by plotting the local maxima of state variable x 1 according to the coupling coefficients that are increased in slight steps. When increasing k 1 or k 2 , it is found that the 5-D TMJ system under consideration can undergo abundant bifurcation scenarios and exhibit rich dynamical behaviors. Take the case in Figure 4a as an example to illustrate the bifurcation route to chaos of the 5-D TMJ system. When increasing k 1 , the running orbit begins with unbounded behavior at k 1 = 0.2, transfers into bounded chaotic behavior at k 1 = 0.24, and rapidly degenerates into periodic behaviors via a serial of period-doubling bifurcation. Afterwards, the running orbit breaks into chaos by tangent bifurcation at k 1 = 0.29. In the relatively wide interval k 1 ∈ (0.29, 1.25) the 5-D TMJ system mainly works in chaos state, but accompanied by some periodic window behaviors. The largest periodic window with period-4 appears in the interval k 1 ∈ (0.78, 0.88), during which the chaos crisis and tangent bifurcation happen naturally. When further increasing k 1 , the running orbit ultimately degenerates into periodic behaviors via serial period-doubling bifurcation. Meanwhile, for positive maximal LE, the 5-D TMJ system evolves irregularly in the folded space of a chaotic attractor, while for zero maximal LE, the 5-D TMJ system oscillates regularly on a limit cycle.
Consequently, the 5-D TMJ system can exhibit complex dynamical behaviors that are reliant on to the coupling coefficients, and its dynamical distributions in the parameter plane are greatly affected by the initial conditions of the system.

Memristor and Non-Memristor Initial-Condition Effects
This section focuses on the dynamical effects of the memristor and non-memristor initial conditions on the 5-D TMJ system through theoretically exploring the stability distribution of the plane equilibrium point and numerically investigating the initial-conditionrelated extreme multi-stability.

Stability Distribution for Plane Equilibrium Point
The equilibrium point plays an essential role in the dynamical behaviors of the system and it can be solved by setting the right-hand side of the system equations to zero. Thus, a plane equilibrium point can be easily obtained from system (6) and it can be described by: where µ, η are two arbitrary constants, representing the initial positions of the equilibrium point S 0 .
Evaluated at the equilibrium point S 0 , the Jacobian matrix J of system (6) can be derived as follows: It is obvious that the equilibrium point S 0 has two zero roots and three nonzero roots, similar to that reported in [36]. For the nonzero roots, the corresponding cubic characteristic equation can be derived from the Jacobian matrix of system (6) at S 0 as: where For the three nonzero roots, Routh-Hurwitz criteria for the above cubic characteristic equation are given by: If the three criteria of (10) are satisfied, the equilibrium point S 0 is stable and a point attractor appears; otherwise, S 0 is unstable and periodic or chaotic behaviors may occur in system (6).
The two changeable parameters are determined as k 1 = 1 and k 2 = 2.3, the initial parameter µ is varied in the region (−0.5, 1), and the initial parameter η is evolved in the region (−3, 3). According to the conditions given in (10), the stability distribution identified by the three nonzero roots in the µ-η plane can be depicted, as shown in Figure 5. The different color blocks are used to divide three categories of stability regions consisting of the unstable region I, stable region II, and unstable region III. The equilibrium point S 0 located in the region I and region III is the unstable saddle-focus (USF), whereas S 0 in the region II is the stable node-focus (SNF). Besides, the three nonzero roots marked by 0P3N represent the three negative real parts with no positive real part, and the three nonzero roots marked by 1P2N and 2P1N represent the one and two positive real parts, respectively. Particularly, the boundary marked by HB line between the region II and region III is the Hopf bifurcation (HB) line. Therefore, different types of attractors' behaviors may emerge in the different stability regions, such as point, periodic, chaotic, and unbounded behaviors, demonstrating that the dynamical behaviors are extremely dependent on the initial positions of the equilibrium point S 0 , i.e., the initial conditions of two memristors. Due to the existence of two zero roots, the aforementioned SNF for S 0 is critically stable, indicating that the stability of system (6) cannot be effectively judged only by the three nonzero roots [19]. Just for intuition's sake, the initial conditions are determined as ICs = (10 −9 , 0, 0, µ, η). For different values of µ and η located in different regions of Figure 5, the three nonzero roots, stability regions, and attractor types are numerically simulated, as listed in Table 1. Consequently, for different values of µ and η, various kinds of attractors appear in system (6), leading to the emergence of extreme multi-stability. The results from Figure 5 and Table 1 demonstrate that, with the variations of the memristor initials µ and η, system (6) derives multiple stability distributions and presents complex dynamical behaviors. Particularly, the unbounded behavior appears in system (6) under a very large region of memristor initial conditions, meaning that system (6) is less robust to the initial conditions. Besides, the stable-point attractors also appear in the unstable region I, and are triggered by the two zero roots of S 0 .

Initial-Condition-Related Extreme Multi-Stability
The local basin of attraction is the attracting region of a steady-state attractor in the plane of two initial conditions, within which any initial conditions will settle down to the attractor. The above analysis results indicate that the 5-D TMJ system can produce extreme multi-stability dependent on the memristor initial conditions. Thus, the memristor initial conditions are taken as two invariant measures for classifying the coexisting infinitely many attractors' behaviors [27,40]. Note that the local basin of attraction is colorfully drawn by calculating the periodicities of the state variable x 1 based on a MATLAB ODE45 algorithm with time-step 0.01 and time-span (400, 1000), regardless of the structure and position of the generated attractor.
The representative parameters are determined as a = 2.5, b = 0.8, k 1 = 1, and k 2 = 2.3, and the memristor initial conditions are denoted as IC 2 = (x 4 (0), x 5 (0)) = (µ, η). To depict the local basin of attraction, we examine the periodicities of the state variable x 1 at each point in the µ-η plane. For two sets of tiny non-memristor initial conditions of IC 1 = (10 −9 , 0, 0) and IC 1 = (−10 −9 , 0, 0), the local basins of attraction in the µ-η plane are measured and the results are shown in Figure 6. The color blocks represent the memristor initial conditions that trigger the trajectories of different dynamical behaviors, which are exactly the same as those used in Figure 2a. As a result, the 5-D TMJ system can exhibit complex dynamical behaviors closely dependent on the memristor initial conditions, indicating the emergence of extreme multi-stability therein.
Particularly, as revealed in Figure 6, the local basins of attraction with IC 1 = (10 −9 , 0, 0) and (−10 −9 , 0, 0) are very different. The attracting regions in Figure 6a are completely symmetric with respect to µ, while the right half of the attracting regions in Figure 6b is exactly the same as that in Figure 6a, but the left half is unbounded. This is shocking that completely different dynamical behaviors appear in the µ-η plane when the initial condition x 1 (0) is slightly varied. Unfortunately, the dynamical mechanism is quite difficult to be identified based on the mathematical model of system (6). Moreover, because there are extremely small differences in the initial condition x 1 (0), the dynamical distributions depicted by Figure 6a are different from those depicted by Figure 2, but the dynamical distributions depicted by Figure 6b are similar to those depicted by Figure 5. This is an unexpected and very interesting result.
For manifesting dynamical effects of the non-memristor initial conditions, two sets of large non-memristor initial conditions of IC 1 = (1, 0, 0) and IC 1 = (0, 1, 0) are employed. Accordingly, the local basins of attraction in the µ-η plane are measured and their results are shown in Figure 7. As can be seen, with the changes of the non-memristor initial conditions, the left and right attracting regions appeared in the cases of tiny non-memristor initial conditions are connected as a whole, and the bounded behaviors are mainly confined to the lower plane regions of the memristor initial conditions. In brief, the simulated results further demonstrate that the non-memristor initial conditions have a great influence on the dynamical behaviors of the 5-D TMJ system. Referring to the local basin of attraction revealed in Figure 6a, the memristor initial condition-induced attractors are measured to verify the coexistence of infinitely many attractors in the 5-D TMJ system. For some representative memristor initial conditions IC 2 , the phase portraits of the induced attractors in the x 4 -x 5 plane are depicted in Figure 8, where the stable point attractors are marked by the symbol "*". For better visual effect, only partial phase portraits of the induced attractors are simulated in Figure 8, but they sufficiently demonstrate the memristor initial-condition-related extreme multi-stability. In fact, more phase portraits of the induced attractors with different structures, periodicities, and positions can be numerically measured from the 5-D TMJ system, indicating the appearance of the coexisting infinitely many attractors' behaviors.

Analog and Digital Implementations
This section designs an analog circuit and a digital microcontroller-based hardware platform for implementing the 5-D TMJ system. The PSIM circuit simulations and experimental measurements acquire the memristor initial condition-induced attractors to validate the numerical simulations.

Analog Circuit Design and PSIM Circuit Simulations
PSIM (Power Simulation) circuit simulations are employed to verify the dynamical effects of the initial conditions on the 5-D TMJ system using a physical circuit. Based on the mathematical model given in (6), an analog implementation circuit in a jerk-circuit form is designed for the 5-D TMJ system and its circuit schematic is shown in Figure 9. Here, the memristor used in system (6) is implemented equivalently by the circuit module shown at the top of Figure 9, where the −tanh(·) module can be defined by a function module in PSIM software. When the two capacitor voltages v 3 and v 1 are taken as the input voltages of two memristors, respectively, the generated currents i 1 and i 2 can be expressed as: where τ 0 = RC is the integration time constant, v 4 and v 5 stand for the capacitor voltages of the memristor circuit modules, g 1 and g 2 represent the multiplier gains of the memristor circuit modules, and R k1 and R k2 are used for adjusting the two coupling coefficients. Usually, the multiplier gains g 1 and g 2 are fixed as 1, and the integration time constant τ 0 = RC set to 10 kΩ × 100 nF = 1 ms. Thus, the resistances R a and R b in the main circuit are determined by R a = R/a and R b = R/b, and the resistances R k1 and R k2 in the memristor circuit modules are obtained by R k1 = R/k 1 and R k2 = R/k 2 . For the representative control parameters a = 2.5 and b = 0.8, as well as the coupling coefficients k 1 = 1 and k 2 = 2.3, the four resistances can be calculated as R a = 4 kΩ, R b = 12.5 kΩ, R k1 = 10 kΩ, and R k2 = 4.3478 kΩ. Besides, two inverters should be connected to the capacitor voltage terminals of the memristor circuit modules to achieve the same results of circuit simulations and numerical simulations. The three initial capacitor voltages in the main circuit are preset as v 1 (0) = 1 nV, v 2 (0) = 0 V, and v 3 (0) = 0 V, and only the two initial capacitor voltages v 4 (0) and v 5 (0) of the memristor circuit modules are adjustable to acquire the desired phase portraits. Corresponding to the numerical simulations shown in Figure 8, the phase portraits of the induced attractors in the plane of v 4 and v 5 are simulated by PSIM software and the measured outputs under different values of v 4 (0) and v 5 (0) are displayed in Figure 10, where IC 2 = (v 4 (0), v 5 (0)). Therefore, the circuit simulations in Figure 10 validate the memristor initial condition-induced attractors revealed in Figure 8, indicating the feasibility of the designed circuit.

Digital Hardware Platform and Experiments
Based on a pony-size and low-cost microcontroller, a digitally circuit-implemented hardware platform can be developed for the 5-D TMJ system, in which STM32F407VET6 microcontroller with 32-bit RISC core and DAC8563 digital-to-analog (D/A) converter are utilized. According to the mathematical model given in (6), a discrete-time mathematical model is built by means of the fourth-order Runge-Kutta algorithm, which can be coded in C language and uploaded to the microcontroller. All the system parameters and nonmemristor initial conditions are determined in advance, and only the memristor initial conditions are adjusted to acquire the desired phase portraits using a digital oscilloscope. Figure 11 displays the digitally circuit-implemented hardware platform for the 5-D TMJ system and the memristor initial condition-induced attractors acquired by the digital oscilloscope. Similarly, following the numerical simulations shown in Figure 8, the phase portraits of the induced attractors in the plane of v 4 and v 5 are acquired, as shown in Figure 12, where the memristor initial conditions IC 2 = (v 4 (0), v 5 (0)) are changed as different values according to those given in Figure 8. Therefore, the experimentally acquired results in Figure 12 validate the numerical results well, reflecting the feasibility of the digitally circuit-implemented hardware platform.

Conclusions
In this paper, a new 5-D nonlinear dynamical system based on two memristors was presented by introducing two memristors to a simple jerk system, and its complex dynamics related to the coupling coefficients under different memristor initial conditions were studied using multiple numerical methods. Emphatically, the dynamical effects of the memristor and non-memristor initial conditions on the 5-D TMJ system were explored by theoretical analyses and numerical simulations, which were validated by PSIM circuit simulations and microcontroller-based hardware experiments. However, because of the two zero roots caused by the plane equilibrium point, the dynamical mechanism of the initialcondition effects on the 5-D TMJ system cannot be effectively explained only by the stability distributions of the three nonzero roots. To address this issue, an incremental integral reconstitution model can be used for elaborating the dynamical mechanism of the initialcondition effects [19,35], which deserves further study.

Data Availability Statement:
Data generated during the current study will be made available at reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.