A physical memristor based Muthuswamy–Chua–Ginoux system

Abstract

In 1976, Leon Chua showed that a thermistor can be modeled as a memristive device. Starting from this statement we designed a circuit that has four circuit elements: a linear passive inductor, a linear passive capacitor, a nonlinear resistor and a thermistor, that is, a nonlinear “locally active” memristor. Thus, the purpose of this work was to use a physical memristor, the thermistor, in a Muthuswamy–Chua chaotic system (circuit) instead of memristor emulators. Such circuit has been modeled by a new three-dimensional autonomous dynamical system exhibiting very particular properties such as the transition from torus breakdown to chaos. Then, mathematical analysis and detailed numerical investigations have enabled to establish that such a transition corresponds to the so-called route to Shilnikov spiral chaos but gives rise to a “double spiral attractor”.

Introduction

Michael Faraday (1791-1867) is generally well known for his contributions to the study of electromagnetism and electrochemistry. However, according to Orton¹, while investigating the effect of temperature on the conductivity of “sulphuret of silver” (silver sulfide) in 1833 he discovered what is now considered as a thermistor². Thermistor, i.e., thermal resistor, is thus an electrical-resistance element made of a semiconducting material the resistance of which varies with temperature. Thermistor production was then difficult, and applications were limited. One century had to pass before thermistors became commonly used by commercial manufacturers. In the 1930s, Samuel Ruben (1900-1988) invented the first commercial thermistor that he called “electrical pyrometer resistance” (Patents N^o 2,021,491). He explained that:

This invention relates to an electrical pyrometer and specifically to one utilizing the resistance change of a metallic compound with heat to indicate temperature changes. The characteristic property of the material employed for the temperature indicating resistance element is one having a high negative resistance coefficient.

In the early 1940s, development of a technique improving the overall consistency and repeatability of many manufacturing processes enabled the large-scale commercialization of thermistor. At that time, thermistor was most frequently used for protection, regulation, and compensating for temperature in electronic circuits. By the 1960s, thermistors were being used in the aerospace industry. Then, John Steinhart and Stanley Hart found a function modeling thermistor characteristics, i.e., the resistance according to the temperature which proved to be suitable for a wide variety of thermistors for ranges of a few degrees to a few hundred degrees³. In 1976, five years after Chua⁴ had postulated a missing circuit element that he called memristor, Chua published a paper with Sun Kang⁵ in which they recalled that:

Thermistors have been widely used as a linear resistor whose resistance varies with the ambient temperature. In particular, a negative-temperature coefficient thermistor is characterized by

$$\begin{aligned} v_T = R_0 \left( T_0 \right) \exp \left[ \beta \left( \frac{1}{T} - \frac{1}{T_0} \right) \right] i {\mathop {=}\limits ^{\triangle }} R\left( T \right) i \end{aligned}$$

(1)

where $\beta$ is the material constant, T is the thermistor temperature and $T_0$ the room temperature both in kelvin. The constant $R_0(T_0)$ denotes the cold temperature resistance at $T = T_0$. The instantaneous temperature T, however, is known to be a function of the power dissipated in the thermistor and is governed by the heat transfer equation

$$\begin{aligned} p\left( T \right) = v_T \left( t \right) i \left( t \right) = \delta \left( T - T_0 \right) + c \frac{dT}{dt} \end{aligned}$$

(2)

where c is the heat capacitance and $\delta$ is the dissipation constant of the thermistor, defined as the ratio of a change in the power dissipation to the resultant change in the body temperature. Substituting (1) into (2) and by rearranging terms, we obtain

$$\begin{aligned} \frac{dT}{dt} = - \frac{\delta }{c} \left( T - T_0 \right) + \frac{R_0 \left( T_0 \right) }{c} \exp \left[ \beta \left( \frac{1}{T} - \frac{1}{T_0} \right) \right] i^2 \end{aligned}$$

(3)

We observe from (1) and (3) that a thermistor is in fact not a memoryless temperature-dependent linear resistor—as usually assumed—but rather a first-order time-invariant current-controlled memristive.

On May 1st, 2008, an ReRAM based memristor was discovered by Strukhov et al.⁶. Their work triggered a renewed interest in memristors and their applications across widely different fields which has consistently grown since. During these last five years, Rajamani et al.⁷ analyzed an electronic oscillator circuit designed by connecting an inductor in series with a “locally-active” Positive Temperature Coefficient (PTC) memristor and a battery. Then, Sah et al.⁸ have considered a “second order memristor which represents the model of a physical device called Positive Temperature Coefficient (PTC) and Negative Temperature Coefficient (NTC) thermistor connected in series.” The aim of this work is to investigate a circuit consisting of a linear passive capacitor, a linear passive inductor, a nonlinear resistor and a Negative Temperature Coefficient thermistor, that is, a nonlinear “locally active” volatile memristor. Notice that a few implementations of chaotic circuits have recently used ReRAM memristive models^9,10,11. However, memristive models for ReRAM devices such as the HP memristor is still the subject of research¹². In contrast, the memristive model for thermistors is well established and the device can be bought as an off-the-shelf component with a variety of parameter options¹³.

Circuit setup

Starting from the Muthuswamy–Chua circuit¹³, we designed a circuit (system) consisting of a linear passive capacitor of capacitance C, a linear passive inductor of inductance L, a nonlinear resistor $N_R$ and a thermistor, that is, a nonlinear “locally active” memristor M. A schematic is shown in Fig 1.

Figure 1

(a) The Muthuswamy–Chua–Ginoux (MCG) circuit is a generalization of the (b) Muthuswamy–Chua circuit, because a “nonlinear resistor in series with a memristor is a memristor” (see Muthuswamy and Banerjee for a proof of this theorem¹³).

By Kirchhoff’s laws, we have: $v_C + v_L + v_R + v_M = 0$ with $v_C = q/C, v_L = L di/dt, v_R = f(i)$ and $v_M=R(T)i$. $v_R$ is the voltage across the the nonlinear resistor modeled by a function f(i) defined below. Taking into account the Eq. (3) and since the current intensity across the capacitor is given by $i = dq/dt$, we obtain the following set of equations:

$$\begin{aligned} \begin{aligned} \dfrac{dv_C}{dt}&= \dfrac{i}{C}, \\ \dfrac{di}{dt}&= - \dfrac{1}{L} \left[ v_C + f\left( i \right) + R\left( T \right) i \right] , \\ \dfrac{dT}{dt}&= \dfrac{ R\left( T \right) }{c}i^2 - \frac{\delta }{c} \left( T - T_0 \right) , \\ \end{aligned} \end{aligned}$$

(4)

Following Balthasar Van de Pol¹⁴, we model the nonlinear resistor as a cubic function of the current, that is, $f(i) = ai + bi^3$ where a and b are constants. The thermistor characteristics (the curve representing the resistance R as a function of the temperature T) is modeled using the classical Steinhart–Hart equation ($T^{-1} = A + B \log R + C (\log R)^3$) which is a third-order approximation³. By considering that the third Steinhart-Hart coefficient $C= 0$ we obtain the $\beta$-parameter equation (1) which thus corresponds to a second-order approximation. Nevertheless, from both analytical and numerical point of view it is quite difficult to investigate the model (4) with exponential functions. Therefore, physically, we can restrict the thermistor to a small temperature change and hence approximate R(T) by its Taylor series expansion up to the second order:

$$\begin{aligned} R\left( T \right) = R_0 \left( T_0 \right) \exp \left[ \beta \left( \frac{1}{T} - \frac{1}{T_0} \right) \right] = R_0 \left[ 1 - \frac{\beta }{T_0^2} \left( T - T_0 \right) + \frac{\beta \left( \beta + 2T_0 \right) }{2 T_0^4} \left( T - T_0 \right) ^2 + O \left( (T-T_0)^3 \right) \right] \end{aligned}$$

(5)

Let us notice that the left hand side of Eq. (5) is a second-order approximation of the Steinhart-Hart equation. As a consequence, its Taylor series expansion must be limited to the second order. Moreover, numerical investigations have shown that if R(T) is approximated by its first order Taylor series expansion, the dynamical system (4) can only exhibit relaxation oscillations and not chaotic behavior. In order to simplify the study of system (4), let’s pose for the variables:

$$\begin{aligned} x = V_C \quad \text{; } \quad y = i \quad \text{; } \quad z = T - T_0. \end{aligned}$$

and for the parameters:

$$\begin{aligned} \alpha = C \quad \text{; } \quad \eta = L \quad \text{; } \quad \theta = \dfrac{R_0}{c} \quad \text{; } \quad \gamma = - \dfrac{R_0}{c}\frac{\beta }{T_0^2} \quad \text{; } \quad \mu = \dfrac{R_0}{c} \frac{\beta \left( \beta + 2T_0 \right) }{2 T_0^4} \quad \text{; } \quad \epsilon = \frac{\delta }{c}. \end{aligned}$$

In Fig. 2, we have plotted both thermistor characteristics and its Taylor series expansion (5) while using the NTC thermistor B57236S0250M000 (This device can be purchased and its specifications can be found from the following website: https://product.tdk.com/info/en/documents/data_sheet/50/db/icl_16/S236.pdf) the parameters of which are the following: $R_0 = 60 \Omega$ and $\beta = 3000$ K.

Figure 2

Thermistor characteristics in blue and its Taylor series expansion in red.

Then, we have computed the coefficient of determination $\mathcal {R}^2$ for $T \in [240, 300]$ to measure how well the Taylor series expansion fits with the thermistor characteristics and found $\mathcal {R}^2 = 0.9726$ (p value < 0.025) which indicates a quite good fit¹⁵.

Thus, we obtain what we call the Muthuswamy–Chua–Ginoux (MCG) system:

$$\begin{aligned} \begin{aligned} \dfrac{dx}{dt}&= \dfrac{y}{\alpha }, \\ \dfrac{dy}{dt}&= - \dfrac{1}{\eta } \left[ x + f\left( y \right) + R\left( z \right) y \right] , \\ \dfrac{dz}{dt}&= R\left( z \right) y^2 - \epsilon z, \\ \end{aligned} \end{aligned}$$

(6)

where $f(y) = ay + by^3$ and $R(z) = \mu z^2 + \gamma z + \theta$.

In order to justify the novelty of our system, we apply the criteria for publication standard for chaotic systems set forth by Sprott¹⁶.

• 1. The system should credibly model some important unsolved problem in nature and shed insight on that problem.

  Ever since the announcement of a memristor by Hewlett-Packard labs⁶, memristor based chaotic circuits abound¹⁷. Nevertheless, such circuits do not use readily available physical memristors such as thermistors, discharge tubes or junction diodes¹³. Our work hence resolves the unsolved question: “can a physical memristor based electronic circuit exhibit chaos”? Moreover, the insight gained is that a nonlinear resistor in series with a memristor (still a memristor, refer to Fig. 1) can be used to appropriately shape the nonlinear memristance for chaotic behaviour.

• 2. The system should exhibit some behavior previously unobserved.

  We observe, for the very first time in a memristor based chaotic circuit, the Shilnikov spiral chaos¹⁸ but with exotic phenomenon such as “multiple” reverse period-doubling cascade.

• 3. The system should be simpler than all other known examples exhibiting the observed behavior.

  Referring to Fig. 1, one can see that we only have four elements in series. Analog memristor based chaotic circuits that use emulators require quite a lot more components¹⁷. Although this paper uses a cubic characteristic for the nonlinear resistor, one could simply use a piecewise-linear characteristic and hence implement the nonlinear resistor in Fig. 1 using a single op-amp¹⁷. Please refer to the appendix for electronic circuit implementation.

Let’s notice that the flow of this system (6) is invariant under the symmetry $(-x,-y,z) \rightarrow (x,y,z)$. Hence if (x(t), y(t), z(t)) is a solution of system (6), so is $(-x(t), -y(t), z(t))$.

Parameters set

From Eq. (5), it follows that its right hand side is strictly positive since it is an exponential function modeling the resistance of the thermistor which is supposed to be positive. As a consequence, the same should be true for its left hand side, i.e., its Taylor series expansion. Thus, the second degree polynomial $R(z) = \mu z^2 + \gamma z + \theta$ must be positive. This implies that its discriminant $\gamma ^2 - 4 \mu \theta$ must be negative and that $\mu$ must be positive. Otherwise, the resistance could take negative values which is impossible from a thermodynamics point of view. Moreover, since $\theta = R_0 / c$ it must also be positive. Notice that the negativity of $\theta$ would imply the positivity of the discriminant which is here precluded. Such considerations lead to the following conditions:

$$\begin{aligned} \mu> 0 \quad \text{; } \quad -2 \sqrt{\mu \theta }< \gamma < 2 \sqrt{\mu \theta } \quad \text{; } \quad \theta > 0. \end{aligned}$$

So, first for the stability analysis of the dynamical system (6) presented in the next section, we will use the following parameters set which meet the previous conditions and which facilitate the identification of the transition from torus breakdown to “double spiral chaos”:

$$\begin{aligned} a = -6, b = 3, \eta = 12.2, \mu = 3, \gamma = -2, \theta = 3, \epsilon = 0.6. \end{aligned}$$

Then, in the Appendix, dynamical system (6) will be recast in dimensionless form which is more suitable for analog simulations with a parameters set including the NTC thermistor B57236S0250M000 features ($R_0 = 60\,\Omega ,\beta = 3000$ K). These analog simulations could be thus used to highlight such a transition.

Stability analysis

Fixed points

Fixed points are determined while using the classical nullclines method. MCG system (6) has the origin O(0, 0, 0) as the unique fixed point.

Jacobian matrix and eigenvalues

The Jacobian matrix of MCG system (6) reads:

$$\begin{aligned} J = \begin{pmatrix} 0 \ &{} \ \dfrac{1}{\alpha } \ &{} \ 0 \\ -\dfrac{1}{\eta } \ &{} \ -\dfrac{1}{\eta } \left( a + 3by^2 + \mu z^2 + \gamma z + \theta \right) \ &{} \ -\dfrac{y}{\eta } \left( 2 \mu z + \gamma \right) \\ 0 \ &{} \ 2 \left( \mu z^2 +\gamma z + \theta \right) y \ &{} \ \left( 2 \mu z + \gamma \right) y^2 - \epsilon \end{pmatrix} \end{aligned}$$

(7)

By replacing the coordinate of the fixed point, i.e., the origin in the Jacobian matrix (7) one obtains the Cayley-Hamilton third degree eigenpolynomial which can be easily factorized as follows:

$$\begin{aligned} \left( \lambda + \epsilon \right) \left[ \alpha \eta \lambda ^2 + \alpha \left( a + \theta \right) \lambda + 1 \right] = 0 \end{aligned}$$

(8)

Thus, there is one real eigenvalue:

$$\begin{aligned} \lambda _1 = - \epsilon \end{aligned}$$

(9)

which is negative since the parameter $\epsilon$ is positive and two eigenvalues:

$$\begin{aligned} \lambda _{2,3} = \frac{1}{2 \eta } \left[ - \left( a + \theta \right) \pm \sqrt{ \left( a + \theta \right) ^2 - 4 \frac{\eta }{\alpha } } \right] . \end{aligned}$$

(10)

Eigenvalues’ properties

The value $\alpha \left( a + \theta \right) ^2 - 4 \eta / \alpha$ depends on $\alpha$ which we will be chosen as bifurcation parameter. Since we have posed $\alpha = C >0$, this value is negative provided that: $0 < \alpha \leqslant 4 \eta / ( a + \theta )^2$. In this stability analysis, we use the following parameters set: $a = -6, \eta = 12.2, \mu = 3, \gamma = -2, \theta = 3, \epsilon = 0.6 \text{ and } b = 3$. So, we have $a + \theta < 0$. Thus, according to Poincaré¹⁹, the real parts of the eigenvalues (10) are positive and the fixed point, i.e., the origin is a saddle node or a saddle focus. More specifically, if $0 < \alpha \leqslant 4 \eta / ( a + \theta )^2$, it is a saddle focus. These results are summarized in the Table 1.

Table 1 Fixed point stability.

To investigate the effects of the changes of parameter $\alpha$ on the dynamics of the MCG system, a bifurcation diagram has been plotted (see Fig. 3). Then the information can be used to have a better understanding of the phase space orbits plotted in Fig. 4. According to I.M. Ovsyannikov and L.P. Shilnikov²⁰:

The transition to spiral chaos along the lines of the above scenario is usually preceded by either a cascade of period doubling bifurcations (if three-dimensional volumes are contracting) or the collapse of a two-dimensional torus arising from L (a stable periodic motion).

Figure 3

Bifurcation diagram $z_{max}$ as function of $\alpha$.

We observe, on the bifurcation diagram (Fig. 3), a “double” period-doubling cascade for $\alpha \in [0.001, 0.3]$ and $\alpha \in [0.9, 1.2]$. This “doubling” is due to the fact that, for some particular values of the parameters set, the trajectory curve take the shape of a “double spiral attractor” (see Fig. 4d). Hence, each spiral can perform a period-doubling cascade. That’s the reason why the bifurcation diagram exhibits two branches for the period-doubling cascade.

Figure 4

Phase portraits of Muthuswamy–Chua–Ginoux system (6) in the phase space for $\alpha = 0.05, 0.1, 0.5$ and 1.2.

The transition from torus breakdown to “double spiral chaos” has been highlighted by plotting the phase portrait of MCG system (6) for various values of $\alpha$ (see Fig. 4) and by comparing it with the bifurcation diagram (see Fig. 3). We observe that for $\alpha = 0.05$ (all other parameters are the same as above), the attractor is a torus (see Fig. 4a). As $\alpha$ increases between 0.08 and 0.1, a torus breakdown is observed and a limit cycle appears (see Fig. 4b). When $\alpha$ increases again between 0.11 and 0.13, a “spiral attractor” appears. Then, for $\alpha \in [0.14, 0.19]$ a period-doubling cascade occurs. For $\alpha = 0.2$ a “spiral chaos” attractor is observed while for $\alpha \in [0.21, 0.23]$, a second period-doubling cascade occurs. For $\alpha \in [0.24, 0.28]$ a limit cycle appears and for $\alpha \in [0.29, 0.30]$, a third period-doubling cascade occurs (see Fig. 3). For $\alpha \in [0.31, 0.98]$, a “spiral chaos” attractor is again observed (see Fig. 4c). For $\alpha = 1$ a limit cycle appears. Finally, after a short fourth period-doubling cascade for $\alpha = 1.1$, a “double spiral attractor” is observed for $\alpha = 1.2$ (see Fig. 4d). The corresponding time series x(t), y(t) and z(t) of Muthuswamy–Chua–Ginoux system (6) have been plotted in Fig. 5 for $\alpha = 0.5$ and 1.2.

Figure 5

Time series of Muthuswamy–Chua–Ginoux system (6) for $\alpha = 0.5$ and 1.2.

In order to confirm such scenario, Lyapunov Characteristic Exponents (LCE) have been computed in each case. The algorithm developed by Sandri²¹ for Mathematica has been used to perform the numerical calculation of the Lyapunov characteristics exponents (LCE) of the Muthuswamy–Chua–Ginoux system (6) in each case. LCEs values have been computed within each considered interval ($\alpha \in [0.001, 0.3]$ and $\alpha \in [0.9, 1.2]$). As an example, for $\alpha = 0.05, 0.1, 0.5$ and 1.2, Sandri’s algorithm has provided respectively the following LCEs $(0, 0, -0.13), (0, -0.08, -0.08), (0.08, 0, -0.4)$ and $(0.073, 0, -0.38)$. Then, following the works of Klein and Baier²², a classification of (autonomous) continuous-time attractors of dynamical system (6) on the basis of their Lyapunov spectrum, together with their Hausdorff dimension is presented in Table 2. LCEs values have been also computed with the Lyapunov Exponents Toolbox (LET) developed by Siu for MatLab and involving the two algorithms proposed by Wolf et al. ²³ and Eckmann and Ruelle²⁴ (see https://fr.mathworks.com/matlabcentral/fileexchange/233-let). Results obtained by both algorithms are consistent.

Table 2 Lyapunov characteristics exponents of dynamical system (6) for various values of $\alpha$.

On Fig. 6, the bifurcation diagram has been plotted for $\alpha \in [0.001, 12]$. Starting from $\alpha = 3$, we observe that a “multiple” reverse period-doubling cascade with 1, 2, 3, 4, 5, ...branches occurs. Such a number corresponds to the period of each spiral of the “double spiral attractor”. As an example, for $\alpha = 7.5$, there are three bifurcation points and so, each spiral has a period 3 while for $\alpha = 11$, there are five bifurcation points and so, each spiral has a period 5 (see Fig. 7).The time series y(t) and z(t) of of Muthuswamy–Chua–Ginoux system (6) plotted in Fig. 8 for $\alpha = 7.5$ and 11 highlights period three and five.

Figure 6

Bifurcation diagram $z_{max}$ as function of $\alpha$.

Figure 7

Phase portraits of Muthuswamy–Chua–Ginoux system (6) in the phase space for various values $\alpha$.

Figure 8

Time series of Muthuswamy–Chua–Ginoux system (6) for $\alpha = 7.5$ and 11.

For $\alpha \in [2.9, 12]$, the bifurcation diagram (see Fig. 6) consists in intervals within which the trajectory curve takes either the shape of a “double spiral attractor” or a “n-periodic limit cycle”.

Discussion

By using a four element circuit comprising a linear passive inductor, a linear passive capacitor, a nonlinear resistor and a thermistor, that is, a nonlinear “locally active” memristor, we designed a new three-dimensional autonomous dynamical system exhibiting very particular properties such as a transition from a limit cycle to double spiral chaos. Then, mathematical analysis and detailed numerical investigations have enabled to establish that such a transition corresponds to the route to Shilnikov spiral chaos originally described by L.P. Shilnikov²⁵ and then by I.M. Ovsyannikov and L.P. Shilnikov²⁶ (see also I.M. Ovsyannikov and L.P. Shilnikov²⁰ and Shilnikov et al. ²⁷) but gives rise to a “double spiral attractor”. Another scenario leading to spiral chaos, as well as double spiral chaos, and also triple spiral chaos has been proposed by B. Deng¹⁸ twenty-five years ago. Nevertheless, if the attractors he obtained seem to be topologically equivalent to those exhibited by the Muthuswamy–Chua–Ginoux system (6) for various values of the bifurcation parameter $\alpha$, the route leading to these spiral and double spiral chaotic attractors is different from the “Z-switches” described by Deng¹⁸ but corresponds to the route to Shilnikov spiral chaos²⁵ as highlighted by the bifurcation diagrams (see Figs. 3, 6) and the computation of Lyapunov Characteristic Exponents (see Table 2). Double scroll attractor (not double spiral attractor) has been used during the 1990s in cryptography²⁸. In 2017, this idea has been patented by Rainer Plaga, Ralph Breithaupt, Sven Freud & Stephan Gieseler, “Chaotic circuit having variable dynamic states as secure information memory,” (https://patents.google.com/patent/WO2017097909A1/en). So, it seems possible to imagine that the double spiral attractor could be used for the same application.