A Simple Parallel Chaotic Circuit Based on Memristor

This paper reports a simple parallel chaotic circuit with only four circuit elements: a capacitor, an inductor, a thermistor, and a linear negative resistor. The proposed system was analyzed with MATLAB R2018 through some numerical methods, such as largest Lyapunov exponent spectrum (LLE), phase diagram, Poincaré map, dynamic map, and time-domain waveform. The results revealed 11 kinds of chaotic attractors, 4 kinds of periodic attractors, and some attractive characteristics (such as coexistence attractors and transient transition behaviors). In addition, spectral entropy and sample entropy are adopted to analyze the phenomenon of coexisting attractors. The theoretical analysis and numerical simulation demonstrate that the system has rich dynamic characteristics.


Introduction
In 1971, Chua predicted the existence of the memristor based on the incompleteness of the relationship among the basic circuit elements composed of resistance, inductance, and capacitance [1]. However, for the next 30 years, memristor had not been taken seriously until D.B. Strukov and coauthors discovered a new TiO 2 nanoscale device that could be ascribed to memristive properties [2]. Then, the memristor was applied in many fields such as artificial neural networks [3], RFID security system, new memory [4], and chaotic circuits.
The memristor is one type of nonlinear element that can easily produce chaotic signals; thus, many mathematical methods have been proposed [5,6], and the design of chaotic circuits based on memristors has attracted extensive close attention [7][8][9][10][11][12][13][14][15][16][17][18][19]. Chua's circuit or its variants also have been examined by replacing the diode with memristive devices. For example, in 2013, Xu proposed a parallel chaotic system formed by an inductance, a capacitance, and a memristor [20]; in 2019, Yuan designed a parallel chaotic circuit made up of a memristor, a memcapacitor, and a meminductor [21]; in 2020, Ma constructed a parallel chaotic circuit containing a memristor, a memcapacitor, and an inductance [22]. Several kinds of filter circuits have been explored by replacing the resistance or capacitance with memristive devices [23,24]. Some chaotic circuits are based on Shinriki's circuit or jerk chaotic system by replacing the nonlinear element or a resistor with a memristive diode bridge [24,25].
In 1976, Chua elaborated that the thermistor together with ionic systems and discharge tubes can be modeled as memristive devices [26]. In fact, in 1833, Michael Faraday discovered that the resistance of a thermistor varies nonlinearly with temperature; in the early 1940s, the thermistor was widely used in business such as circuit protection as well as temperature sensor and regulation. In 2020, Ginoux introduced a series chaotic circuit based on a thermistor [27]. In this work, the thermistor is applied in a parallel chaotic circuit to construct a new three-dimensional chaotic system for the first time. Although the high-precision thermistor and inductor are difficult to obtain for producing chaos, we provide a possibility of producing chaos in theory. From the simulation results, 11 chaotic attractors and 4 periodic attractors have been found by altering the parameters of the proposed system. In addition, we investigate the coexisting attractors and use the spectral entropy and sample entropy algorithms to obtain the accurate initial values that can generate obviously different attractors. Furthermore, transient transition behavior has also been found in this system. Thus, the proposed system has complex dynamic behaviors.
The rest of this article is arranged as follows: Section 2 gives the mathematical model of the thermistor that is used to construct the parallel chaotic circuit. Section 3 is concerned with numerical analysis and describes the use of MATLAB R2018 to perform the numerical simulation. Bifurcation diagrams combined with their corresponding largest Lyapunov exponents (LLE) are presented to reveal the impacts of parameters. Spectral entropy (SE), sample entropy (SampE), and attractor phase diagrams are mainly used to analyze coexistence attractors. The transient transition behaviors are also discussed using Lyapunov exponents, time-domain waveforms, and attractor phase diagrams. Section 4 summarizes this work.

Model of Thermistor
A negative temperature coefficient (NTC) thermistor is described by Equation (1): In Equation (1), β is the material constant, T and T 0 separately show the thermistor temperature and room temperature with the unit of Kelvin, and R 0 (T 0 ) and R(T) are the resistance respectively at the temperature of T 0 and T. Thus, the conductance G(T) can be obtained as below.
One can obtain the approximate G(T) through Taylor series expansion to the second order shown in Equation (3): From [26], the instantaneous temperature T is known to be a function of the power dissipated in the thermistor and governed by the heat transfer equation in Equation (4): In Equation (4), c is the heat capacitance and δ is the dissipation constant [26] of the thermistor that is the ratio of the power dissipation to the temperature. Thus, based on Equations (1), (2) and (4), we obtained the first-order time-invariant equation shown in Equation (5): Chua and Kang have defined the mathematical model of memristor system, expressed as follows [26]: where u, y, and z respectively represent the memristor's input, output, and the state of the memristor. Based on Equations (3) and (5), the model of the thermistor can be obtained: where v M and i M respectively represent the voltage and current through the thermistor, while T denotes the state of the thermistor. To simplify the system, let T − T 0 = z, δ/c = k, and 1/c = ϕ; in Equation (7), let G 0 = θ, βG 0 /T 0 2 = η, G 0 β(β − 2T 0 )/T 0 4 = λ, and G(T) = θ + ηz + λz 2 ; then, the system can be expressed as In 1976, Chua and his team extended the definition of the memristor [26]. By definition, the relationship between the current and voltage of the memristor system has the characteristics of compact hysteresis loop, zero crossing (passivity criterion), and closure theorem.
Assuming the voltage of the thermistor to be a sinusoidal signal with 2 V and setting the initial value as z 0 = 0, the thermistor's i-v characteristic curve is displayed in Figure 1a. We can find that the trajectory of the i-v curve is similar to the number "8", and the area of the hysteresis loop decreases with the increase in the excitation frequency f. The phenomenon verifies that the thermistor is consistent with the characteristics of a compact hysteresis loop. Moreover, for a passive memristor system, its terminal voltage and current should cross the zero point at the same time, as shown in Figure 1b. Furthermore, the dimension of a thermistor is the same as that of a common resistance; that is, a thermistor can be used in series and parallel as a common resistance, which indicates that the thermistor conforms to the closure theorem. Thus, a thermistor can be considered as a memristor system. Chua and Kang have defined the mathematical model of memristor system, expressed as follows [26]: where u, y, and z respectively represent the memristor's input, output, and the state of the memristor. Based on Equations (3) and (5), the model of the thermistor can be obtained: where vM and iM respectively represent the voltage and current through the thermistor, while T denotes the state of the thermistor. To simplify the system, let T − T0 = z, /c = k, and 1/c = ; in Equation (7), let G0 = , G0/T0 2 = , G0( − 2T0)/T0 4 = , and G(T) =  + z + z 2 ; then, the system can be expressed as In 1976, Chua and his team extended the definition of the memristor [26]. By definition, the relationship between the current and voltage of the memristor system has the characteristics of compact hysteresis loop, zero crossing (passivity criterion), and closure theorem.
Assuming the voltage of the thermistor to be a sinusoidal signal with 2 V and setting the initial value as z0 = 0, the thermistor's i-v characteristic curve is displayed in Figure 1a. We can find that the trajectory of the i-v curve is similar to the number "8", and the area of the hysteresis loop decreases with the increase in the excitation frequency f. The phenomenon verifies that the thermistor is consistent with the characteristics of a compact hysteresis loop. Moreover, for a passive memristor system, its terminal voltage and current should cross the zero point at the same time, as shown in Figure 1b. Furthermore, the dimension of a thermistor is the same as that of a common resistance; that is, a thermistor can be used in series and parallel as a common resistance, which indicates that the thermistor conforms to the closure theorem. Thus, a thermistor can be considered as a memristor system.

A Simple Chaotic Oscillator
A simple chaotic circuit (see Figure 2) is only composed of four elements: a capacitor, an inductor, a thermistor, and a linear negative resistor. Based on Kirchhoff's laws and  (8)), we get the system described in Equation (9). To simplify the system, let i L (t) = x, u C (t) = y, 1/L = a, 1/C = b, and g = r; then, the system can be expressed as The parameters we have set have some practical significance, so we need to discuss the spans. Equation (3) indicates the conductance of the thermistor that is supposed to be positive; thus, G(T) = θ + ηz + λz 2 is positive, and hence, η 2 − 4λθ must be negative and λ must be positive. Otherwise, the conductance may have some negative values. Other parameters a, b, θ, k, ϕ, and γ are positive according to their physical meanings, and η = βG 0 /T 0 2 (β is the material constant) is also positive. Thus, we get the value ranges as follows:

of 14
Obviously, if the characteristic values meet the conditions mentioned above, the equilibrium point is a saddle focus or a saddle node, which is the necessary condition for generating chaotic attractors.
Based on Equations (11) and (14), we give the preliminary value range of the parameters, which may create chaos: According to the conditions discussed above, we choose a set of data a = 1, b = 1, r = 1.5, k = 2, ϕ = 1, η = 1, λ = 2, θ = 0.5 and the initial conditions (0, 0.1, 0); then, Equation (10) can be expressed by Equation (16): Figure 3 shows the chaotic characteristics of the system. Figure 3a-d separately displays the 2D and 3D phase portraits, exhibiting the chaotic behavior: instability and finiteness. Figure 3e shows the Lyapunov exponents (LEs) of the system obtained by the Wolf algorithm [28,29]. The results are LE 1 = 0.18728, LE 2 = −0.0044319, and LE 3 = −0.87559. LE 1 is positive, LE 3 is negative, LE 2 is approximately equal to zero, and the sum of these LEs is negative, so it is a stable chaotic system. Figure 3f is the Poincaré map [30][31][32][33][34] of the system that contains random locations of dots, which also indicates the properties of chaos. In brief, it is a stable chaotic system.

The Impacts of Parameters
We observe the chaotic characteristics of Equation (10) through alerting the parameters a, b, r and k with the initial values as (0, 0.1, 0). The calculated results of the largest Lyapunov exponent spectrum (LLE), bifurcation diagram, attractor phase diagrams, and dynamic maps are used to analyze the impact of parameters a, b, r, and k. First, we determine the state of the system by comparing the value of LLE and bifurcation diagram; then, attractor phase diagrams and dynamic maps are used to future verify the results. Through multiple analyses, we can obtain reliable results. Figure 4 shows the attractor phase diagrams that fully illustrate the state of the proposed system. Through analysis and generalization, 11 kinds of chaotic attractors and 4 kinds of periodic attractors have been found.
The largest Lyapunov exponent spectrum (LLE), bifurcation diagram, attractor phase diagrams, and dynamic maps are obtained by the fourth-order Runge-Kutta method with double precision using the time-step ∆t = 10 −2 s. It is well known that all the numerical algorithms will produce the numerical noises such as truncation error and round-off error. For a system of ODEs, set dϕ/dt = F (ϕ), where F is differentiable, if |∂F i / ∂ϕ j | < L (for all i and j), Then, Equation (17) can be obtained [35] within the time of [0, T]: where N is the order accuracy, D is a constant due to F, ∆t is the time step, E 0 is the initial error, and ϕ t n,∆t and ϕ t a respectively represent numerical solution and analytic solution.
For the proposed chaotic system, we have fixed the values of N, D, ∆t, and E 0 ; thus, as the simulation time T increases, the error will increase. We calculated the LLE with T = 500 s and T = 1000 s, as shown in Figures 5a-8a. From the results, we can find some differences between the two sets of values; moreover, the results of LLE with T = 500 s are better correlated with bifurcation diagrams in Figures 5b-8b. Then, through considering the numerical noises and referring to other articles [22], we set the simulation time to 500 s.      Figure 5a shows the largest Lyapunov spectrum versus the parameter a, and Figure 5b exhibits the bifurcation diagram of the variable z versus the parameter a. The parameter a is increasing from 0 to 2, during which there are three kinds of chaotic attractors and one periodic attractor specifically summarized in Table 1.  (10) when varying the parameter a for the selected set of b = 1, r = 1.5, k = 2, ϕ = 1, η = 1, λ = 2, θ = 0.5 and the initial conditions (0, 0.1, 0).

Range (a)
LEs The largest Lyapunov exponent spectrum and corresponding bifurcation diagram versus the parameter b with the range from 0 to 6 are in Figure 6. Five kinds of chaotic attractors and one periodic attractor have been found. During b ∈ [2.75, 6], the chaotic attractors are transformed from Type II to Type VI at b = 3.45, and during b ∈ [0.06, 0.19], the chaotic attractor is of Type V, which is the only time this type is found in our study. All the dynamic characteristics with varying b are presented in Table 2. Table 2. Corresponding states and LEs of Equation (10) when varying the parameter b for the selected set of a = 1, r = 1.5, k = 2, ϕ = 1, η = 1, λ = 2, θ = 0.5 and the initial conditions (0, 0.1, 0).

Range (b)
LEs  Figure 7, the largest Lyapunov exponent spectrum and corresponding bifurcation diagram versus the parameter r with the range from 1 to 6 are displayed. Especially when r ∈ [1.00, 1.03] and r ∈ [1.08, 1.15], the periodic attractors are separately of Type II and Type III that almost symmetric with respect to the z-axis. All the dynamic characteristics with varying r are presented in Table 3. Table 3. Corresponding states and LEs of Equation (10) when varying the parameter r for the selected set of a = 1, b = 1, k = 2, ϕ = 1, η = 1, λ = 2, θ = 0.5 and the initial conditions (0, 0.1, 0).

Range (r)
LEs State Attractor Type The largest Lyapunov exponent spectrum and corresponding bifurcation diagram versus the parameter k with the range from 0 to 0.5 are presented in Figure 8. Three kinds of chaotic attractors and one periodic attractor have been found. In particular, the chaotic attractor of Type X only exists at k = 0.13 in our discussion area. All the dynamic characteristics with varying k are presented in Table 4. Table 4. Corresponding states and LEs of Equation (10) when varying the parameter k for the selected set of a = 1, b = 1, ϕ = 1, η = 1, λ = 2, r = 2.5, θ = 0.5 and the initial conditions (0, 0.1, 0).

Range (k)
LEs To better illustrate the complex behavior of the circuit, we have given the dynamical maps of the proposed system between the parameters a, b, r, and k depicted in Figure 9. The dynamic maps have been calculated by the method of largest Lyapunov exponent spectrum (LLE). Recently, the relative integration error (RIE) algorithm based on semiimplicit integration schemes has been proposed for the calculation of dynamic maps [30]. In the dynamic maps below, the red regions indicate the chaotic state, the yellow regions indicate the weak chaotic state, the blue regions indicate the periodic or quasiperiodic state, and the violet regions indicate the unbounded state. From Figure 9, we can intuitively understand the changes of the system.

Coexistence of Attractors
Coexistence of attractors is a particular phenomenon in which one system's orbits may gradually change to other stable states with the same parameters but different initial values. In the section, some numerical methods have been used to analyze the coexistence of attractors, including spectral entropy (SE), algorithm sample entropy (SampE) algorithm, largest Lyapunov spectrum (LLE), and attractor phase diagrams. We think different values of SE and SampE imply different states which can be verified by attractor phase diagrams and LLE.
As we know, spectral entropy takes the maximum value when all the probability values of the relative power spectrum for one sequence are equal; that is, the more balanced the distribution of the power spectrum, the bigger the value of SE [32]. Thus, a bigger SE stands for more complex oscillation for the chaotic signal; otherwise, a smaller SE stands for more obvious oscillation for the chaotic signal. Sample entropy is measured by the probability of generating new patterns in signals. The greater the probability, the greater the complexity. Thus, different values of spectral entropy or sample entropy correspond to different chaotic signal trajectories.
From Figure 10a, it can be seen that the value of sample entropy (SampE) is about 0.1, higher than that of spectral entropy (SE), and both sample entropy and spectral entropy have the same pattern of changing. We chose two points (u = 0.1408, SE = 0.4968, SampE = 0.6083) and (u = 0.149, SE = 0.188, SampE = 0.3085), marked in Figure 10a. As expected, different attractors are obtained in Figure 10c,e. From Figure 10b, we discovered that the attractors are in the chaotic state. At the same time, we discovered that the attractors are almost symmetric with respect to the z−axis with setting the negative initial value (0, −0.1408, 0) and (0, −0.149, 0) shown in Figure 10d,f. In brief, two values of variable u can lead to obviously different attractors. In fact, the entropy diagrams shown in Figure 10a reveal that the values of SE and SampE are varying with the variable u, and the complexity of the system is altered, so we can infer that the system has a large number of coexistence attractors.

Transient Transition Behaviors
Transient transition behavior is a special phenomenon in which the state of the circuit is unstable and the system will depart from one dynamic region to another at some point. Sometimes, transient behaviors may include transitions, i.e., transition chaos, transition period, and transition hyperchaos. The proposed system also presents the unstable phenomenon.
Let the parameters a = 1, b = 1, r = 1.5, k = 2, ϕ = 1, η = 1, λ = 2, and θ = 0.5, selecting the initial values (0, 3.75, 0), setting the simulation time as 5000 s. The time-domain diagram displayed in Figure 11c shows that there are two time-period boundaries of t ∈ (0 s, 891.5 s) and t ∈ (891.5 s, 5000 s), which implies that the system should have two different transient transition states. Through numerical simulation, we observed that the chaotic attractor has a transition from Figure 11a  Setting the parameters a = 0.5, b = 1, r = 2.5, k = 0.23, ϕ = 1, η = 1, λ = 2, and θ = 0.5, the initial values (−1, −1, 0), and the simulation time as 5000 s, Figure 11f shows the time domain diagram in which there are two obviously different domains respectively in red and blue. When calculating the LEs within t∈(0 s, 1195 s) and t∈(2075 s, 4062 s), the corresponding values are (0.012051, −0.0059675, −0.058887) and (0.0073981, −0.0042146, −0.037354); thus, the system is in a chaotic state with different attractors. The phase diagrams on the y-z plane shown in Figure 11d,f fully illustrate the transient behavior of the system.
In summary, this section mainly introduces the transition state behavior of the system. The numerical methods of time-domain waveforms, attractor phase diagrams, and Lyapunov exponents are given in detail, and all the simulation results confirm the proposed system has transition chaos behaviors.

Conclusions
In this work, we have designed a chaotic system and found that the thermistor can produce complex nonlinear behaviors. The numerical simulation shows that the system has 11 chaotic attractors and 4 periodic attractors when the circuit parameters change; hence, the system is complex and sensitive. Furthermore, the phenomena of coexistence attractors and transient transition behaviors also exhibit that the system has very rich dynamic characteristics.