Symplectic Dynamics and Simultaneous Resonance Analysis of Memristor Circuit Based on Its van der Pol Oscillator

Abstract

A non-autonomous memristor circuit based on van der Pol oscillator with double periodically forcing term is presented and discussed. Firstly, the differences of the van der Pol oscillation of memristor model between Euler method and symplectic Euler method, four-order Runge–Kutta method (RK4) and four-order symplectic Runge–Kutta–Nyström method (SRKN4), symplectic Euler method and RK4 method, and symplectic Euler method and SRKN4 method in preserving structure are compared from theoretical and numerical simulations, the symmetry and structure preserving and numerical stability of symplectic scheme are demonstrated. Moreover, the analytic solution of the primary and subharmonic simultaneous resonance of this system is obtained by using the multi-scale method. Finally, based on the resonance relation of the system, the chaotic dynamics behaviors with different parameters are studied.

1. Introduction

In 1927, the illustrious van der Pol equation was proposed to describe the oscillation effect of triodes in electronic circuits [1]. Then, in physics [2,3,4], biology [5,6] and economics [7], the van der Pol equation has become a fundamental model to describe the oscillation process. Meanwhile, the memristor was defined in 1971 as the fourth basic circuit device by Academician Chua, following the three basic Double-terminal circuit elements of resistors, capacitors and inductors [8], it can be regarded as a time-varying resistor, and its resistance depends on the amount of charge passing through it. As a fundamental passive circuit element, memristor will be widely used in circuit design as a general RLC element [9,10,11]. It is well-known that Chua’s circuit can be realized by only three elements: inductors, capacitors, and voltage-controlled memristors [12], and van der Pol equation can be realized by two elements: inductors and memristor [13]. Additionally, the time-varying characteristics of memristor give it rich dynamical behaviors [14], and a series of chaotic oscillators are established by using memristors as nonlinear components [15,16]. However, most memristor chaotic circuits developed in recent years are autonomous systems [17,18,19,20], there are few non-autonomous systems with periodically forcing terms, which introduce chaos through forcing terms and other parameters, increase regulations and are easier to be applied in practice. Therefore, it is fascinating to study the non-autonomous memristor circuit based on a van der Pol oscillator with a periodically forcing term.

All along, the numerical method has been using difference approximation to approximate calculus, but its disadvantage is that it focuses only on the extent to which the numerical solution approximates the exact solution, without giving more consideration to the geometric properties of the mechanical system itself as described by the mathematical model. Therefore, a new direction in numerical analysis is that geometric numerical integration is paid more and more attention [21,22]. The well known unpublished report by De Vogelaere is the first demonstration of the existence of numerical integrators (now called symplectic integrators) [23], and it was not until later that Feng proposed a new method to numerically integrate the Hamiltonian system based on symplectic geometry [24], and this implicit method was later applied by Channell [25]. Meanwhile, Sanz-Sema [26] and Lasagni [27] found a condition for implicit Runge–Kutta methods to be symplectic, and Ruth developed an idea of explicit symplectic schemes for the Hamiltonian of the form $H=T\left(p\right)+V\left(q\right)$ [28]. Along the line of explicit symplectic schemes, higher order integrators were presented by Yoshida later on [29]. Then, Qin proposed a multi-step symplectic method and applied it to the numerical analysis of wave problems [30], Forest and Ruth constructed the symplectic maps for nonlinear motion of particles in accelerators [31], Yoshida modified non-existence of first integral by symplectic integration methods [32], Cieśliński compares several discretizations of the simple pendulum equation in a series of numerical experiments and puts forward a new numerical scheme of improved discrete gradient method [33]. Subsequently, the development, analysis and use of various numerical solution algorithms for differential equations concerned in geometric numerical integration [34,35,36,37,38,39]. These algorithms preserve the geometric or qualitative properties of the exact solutions, such as an integral or symmetry, or preservation of a differential invariant such as symplecticity or phase space volume. In engineering, especially for practical systems such as van der Pol and Duffing oscillator, classic numerical method are used commonly, but it has dissipation mechanisms, and errors will accumulate. Therefore, it is very urgent to study the structure preserving methods for these systems.

The path to chaos is also one of the fundamental problems in the study of dynamic behaviors of specific systems [40,41,42,43,44]. However, most of its studies are inclined towards Duffing oscillators, and the studies on periodic oscillation solutions of these systems can be divided into two categories. One is to establish a dynamic model more in line with engineering practice by considering the complexity of structure from the perspective of system structure [45]. The other is to study the dynamical behaviors of the system under complex excitation and the path leading to chaos [46,47]. Nonlinear systems are prone to resonance. In addition to the primary resonance [48], there also have some special resonance phenomena, such as super-harmonic resonance [49], sub-harmonic resonance [50], etc. However, most studies are still limited to the response law of nonlinear system under single harmonic excitation. In practical engineering, combination resonance [51] or simultaneous resonance [52,53] often occurs under multi-frequency excitation, but they are mentioned rarely, especially for van der Pol equation. The paper is organized as follows: a non-autonomous memristor circuit based on van der Pol oscillator with a double periodically forcing term is introduced, and the zero intersection and frequency-dependence characteristics of the memristor are verified by Multisim in Section 2. In Section 3, the stability and dynamic characteristics are studied, and the differences between Euler method and symplectic Euler method, four-order Runge–Kutta method and four-order Symplectic Runge–Kutta–Nyström method, symplectic Euler method and four-order Runge–Kutta method in preserving symmetry and structure are compared from theoretical and numerical simulation. In Section 4, the analytic solution of the primary and subharmonic (1/3 order) simultaneous resonance of this system is obtained by using the multi-scale method, then the chaotic dynamical behaviors are studied from the resonance relationship. Finally, concluding remarks are given in Section 5.

2. Memristor Circuit Based on van der Pol Oscillation Model

The memristor is a device used to describe the relationship between the magnetic flux $\phi$ and the charge q in electronic circuits, and it can be expressed by $f(\phi ,q)=0$. Memristors include charge-controlled memristor and flux-controlled memristor [54]. If the $f(\phi ,q)$ is expressed as the single-valued function of charge, the charge-controlled memristor model is $M\left(q\right)=\frac{d\phi \left(q\right)}{dq}$, and it is expressed as the single-valued function of voltage, the voltage-controlled memristor model is $W\left(\phi \right)=\frac{dq\left(\phi \right)}{d\phi }$, similarly. Where, $R\left(q\right)$ and $W\left(q\right)$ represent the memristance and memductance of the memristor, respectively, and they are nonlinear functions that satisfy $M\left(q\right)=\frac{1}{W\left(\phi \right)}$, similar to the resistance and conductance of the resistor.

Van der Pol first derived the van der Pol equation in order to describe the oscillating effect of triodes in electronic circuits. The van der Pol system is a kind of classical self-excited oscillation system. As an important mathematical model, the van der Pol system can be widely used in complex dynamic system modeling. It provides solutions for many practical engineering problems, and its standard mathematical form is as follows:

$$\left\{\begin{array}{c}\dot{x}=y\hfill \\ \dot{y}=-x+\mu (1-x^2)y+Acos\omega t\hfill \end{array}\right.$$

(1)

where, $\mu$ is nonlinear parameter, A is amplitude of applied excitation and $\omega$ is frequency of applied excitation. On the basis of the form of the system (1), Lu designed the second-order nonautonomous van der Pol oscillator in the Chua’s diode [55]. So, based on that, we add excitation source, and built a van der Pol oscillator based on this cubic memristor, as is shown in Figure 1.

The capacitor voltage $u_c$ and the magnetic flux $\phi$ of the memristor are selected as state variables, and the equation of state of the circuit can be written as

$$\left\{\begin{array}{c}\dot{\phi }=u_c\hfill \\ {\dot{u}}_c=\frac{-u_1\cos {\omega }_{1}t-u_2\cos {\omega }_{2}t-u_c}{RC}-\frac{W\left(\phi \right)u_c}{C}-\frac{\phi }{LC}\hfill \end{array}\right.$$

(2)

where, $W\left(\phi \right)$ is the memductance, and its form is $W\left(\phi \right)=\alpha {\phi }^2+\beta \phi +\gamma$, which implies $q\left(\phi \right)=\frac{\alpha }{3}{\phi }^3+\frac{\beta }{2}{\phi }^2+\gamma \phi$. Additionally, to explore the circuit properties of the memristor, a corresponding emulator is designed, as shown inside the frame in Figure 2.

The voltage between the capacitor electrodes is defined as the internal state of the memristor, and the applied voltage and output current are set as $v_i$ and i, respectively. $U1$ is inverting amplifier; $A1,A2,A3$ are multipliers; $U2,U3,U4$ are inverse integrators and $R1,R2,R3,R4,R5,R6,R7,R8,R9$ are resistors of the memristor, the simulation results of the hysteresis curves and time series are observed, as shown in Figure 3.

As can be seen from Figure 3a, the zero intersection and frequency-dependence characteristics of the memristor are verified by Multisim, and there is an asymmetric pinched hysteresis loop with frequency f = 30 Hz and amplitude voltage $v_i=5$ V, which verifies the existence of non-ideal memrisitve behavior of the memristor.

Then, let $x=\phi$, $y=u_c$, we have

$$\left\{\begin{array}{c}\dot{x}=y\hfill \\ \dot{y}=\frac{1-\gamma R}{RC}\left(1-\frac{\alpha R}{1-\gamma R}{x}^2-\frac{\beta R}{1-\gamma R}x\right)y-\frac{1}{LC}x+\frac{u_1}{RC}\cos {\omega }_{1}t+\frac{u_2}{RC}\cos {\omega }_{2}t\hfill \end{array}\right.$$

(3)

The Equation (3) can be rewritten as a second-order nonautonomous differential equation:

$$\ddot{x}-\xi (1-{\Delta }_1{x}^2-{\Delta }_{2}x)\dot{x}+{\omega }_0^{2}x=A_{1}cos{\omega }_{1}t+A_{2}cos{\omega }_{2}t$$

(4)

where, $\xi =\frac{1-\gamma R}{RC}$, ${\Delta }_1=\frac{\alpha R}{1-\gamma R}$, ${\Delta }_2=\frac{\beta R}{1-\gamma R}$, $A_1=\frac{u_1}{RC}$, $A_2=\frac{u_2}{RC}$, ${\omega }_0^2=\frac{1}{LC}$, here $\xi$ denotes stiffness parameter and $\xi >0$, ${\omega }_0$ is the natural frequency of this model. Additionally, the circuit model has the form of a van der Pol equation from Figure 1.

3. Symplectic Dynamic Analysis of van der Pol Self-Excited Oscillator

3.1. Equilibria and Stability

Considering remove the external excitation, the Equation (4) can be expressed as

$$\ddot{x}-\xi (1-{\Delta }_1{x}^2-{\Delta }_{2}x)\dot{x}+{\omega }_0^{2}x=0$$

(5)

and Equation (5) can be rewritten as

$$\left\{\begin{array}{c}\dot{x}=y\hfill \\ \dot{y}=\xi (1-{\Delta }_1{x}^2-{\Delta }_{2}x)y-{\omega }_0^{2}x\hfill \end{array}\right.$$

(6)

Hence, the equilibrium of the system is $O(0,0)$, and the Jacobian matrix at this equilibrium is

$$J=\left[\begin{array}{cc}0& 1\\ -{\omega }_0^2& \xi \end{array}\right]$$

and the corresponding eigenvalue equation can be written as ${\lambda }^2-\epsilon \lambda +{\omega }_0^2=0$, Then, it is easy to obtain the eigenvalues as

$${\lambda }_{1,2}=\frac{1}{2}\xi \pm \frac{1}{2}\sqrt{{\xi }^2-{\left(2{\omega }_0\right)}^2}$$

(7)

Lemma 1

(Lyapunov’s theorem of stability [56]). If $X^{\prime }=f(t,x)$, $f(t,x)=f(t+T,x)$, $f(t,x)$ and $\frac{\partial f}{\partial x}$ are continuous. $x\subset R^n$ has a periodic solution $X^{*}\left(t\right)$. If the eigenvalue of the linear variational equation $Y^{\prime }=\frac{\partial f}{\partial x}(t,X^{*}\left(t\right))Y$ has strictly negative real parts, then $X^{*}\left(t\right)$ is uniformly asymptotically stable. If there is an eigenvalue with strictly positive real parts, then $X^{*}\left(t\right)$ is unstable.

Lemma 2

(Poincaré-Andronov-Hopf bifurcation existence theorem [57]). Assumption $X^{\prime }=A\left(\lambda \right)+F(\lambda ,x)$ is a $C^k$, where $k\ge 3$, and the planar vector field depends on the scalar parameter λ, makes $F(\lambda ,0)=0$, $D_{x}F(\lambda ,0)=0$, where $\left|\lambda \right|$ is small sufficiently, if the equation of the linear part of $A\left(\lambda \right)$ has a pair of pure imaginary root $\alpha \left(\lambda \right)+i\beta \left(\lambda \right)$ at the $O(0,0)$, and eigenvalue with nonzero velocity across the imaginary axis, namely $\frac{d\alpha \left(0\right)}{d\lambda }\ne 0$. Then, there is $\left|\overline{\lambda }\right|<{\lambda }_0$, for any $U\in R^2$, $O(0,0)\in U$, for any ${\lambda }_0>0$, ${\lambda }_0$ is relatively fixed, subject to the equation $X^{\prime }$ has a nontrivial periodic orbits in U.

According to Lemma 1, the relationship between equilibrium state and $\xi$ can be obtained as shown in the

Table 1. The equilibrium state of system (6).

Conditions of $\mathit{\xi }$	Equilibria	Equilibria Properties
$\xi >2{\omega }_0$	$O\left(0,0\right)$	Unstable Node
$\xi =2{\omega }_0$	$O\left(0,0\right)$	Unstable Degenerate Node
$0<\xi <2{\omega }_0$	$O\left(0,0\right)$	Unstable Focus
$\xi =0$	$O\left(0,0\right)$	Center
$-2{\omega }_0<\xi <0$	$O\left(0,0\right)$	Stable Focus
$\xi \le -2{\omega }_0$	$O\left(0,0\right)$	Stable Node

.

From Table 1, when $\xi >0$, equilibria are unstable, and the phase trajectory is asymptotic to the equilibrium, the movement of phase points along the phase trace deviates from the equilibrium to generate an isolated limit cycle. Meanwhile, when $\xi$ changes from positive to negative, the equilibrium loses stability, because $\frac{d\mathrm{Re}{\lambda }_{1,2}}{d\xi }=\frac{1}{2}\ne 0$. According to Lemma 2, Hopf bifurcation is generated at $\xi =0$, and bifurcation produces a family of periodic solutions near the equilibrium $O(0,0)$.

3.2. Hamiltonian and Exact Solution of Oscillator

3.2.1. Exact Solution Method

Consider $\xi$ is a disturbance term, In order to study the exactness of the solution, here let $\xi =0$, and the system (6) becomes

$$\left\{\begin{array}{c}\dot{x}=y\hfill \\ \dot{y}=-{\omega }_0^{2}x\hfill \end{array}\right.$$

(8)

Obviously, the system (8) is the Hamilton canonical equation, and the corresponding Hamiltonian can be defined as

$$H(x,y)=\frac{{\omega }_0^2{x}^2+y^2}{2}$$

(9)

According to Table 1, when $\xi =0$, its phase trajectory is a closed-loop orbit centered at $O(0,0)$. Here, we take the starting point $(0,Y_0)$ on y-axis, then, the undisturbed parameter equation with time [58] is

$$\left\{\begin{array}{c}x\left(t\right)=Y_{0}sin{\omega }_{0}t\hfill \\ y\left(t\right)={\omega }_0{Y}_{0}cos{\omega }_{0}t\hfill \end{array}\right.$$

(10)

For Equation (10), when $t=0$, we have $x\left(t\right)=0$, $y\left(0\right)={\omega }_0{Y}_0$. Additionally, the starting point of the exact solution (symplectic scheme) is $(0,{\omega }_0{Y}_0)$, the Hamiltonian becomes $H(x,y)=\frac{{\omega }_0^2{Y}_0^2}{2}$. When $\xi \ne 0$, the system (6) can be written in the form of a Hamiltonian with periodic perturbation:

$$\left\{\begin{array}{c}\dot{x}=\xi f_1(x,y)+\frac{\partial H(x,y)}{\partial y}\hfill \\ \dot{y}=\xi f_2(x,y)-\frac{\partial H(x,y)}{\partial x}\hfill \end{array}\right.$$

(11)

where, $f_1(x,y)=0$ and $f_2(x,y)=y(1-{\Delta }_1{x}^2-{\Delta }_{2}x)$.

Then, by the Melnikov method, we have Melnikov function

$$M\left(Y_0\right)={\int }_0^{T\left(Y_0\right)}(\frac{\partial H(x,y)}{\partial x}{f}_1(x,y)+\frac{\partial H(x,y)}{\partial y}{f}_2(x,y))={\int }_0^{\frac{2\pi }{{\omega }_0}}{y}^2(1-{\Delta }_1{x}^2-{\Delta }_{2}x)dt$$

(12)

substitute Equation (10) into the Equation (12) to obtain

$$\begin{array}{c}M\left(Y_0\right)={\int }_0^{\frac{2\pi }{{\omega }_0}}{({\omega }_0{Y}_{0}cos{\omega }_{0}t)}^2(1-{\Delta }_1{(Y_{0}sin{\omega }_{0}t)}^2-{\Delta }_2{Y}_{0}sin{\omega }_{0}t)dt\hfill \\ ={\int }_0^{\frac{2\pi }{{\omega }_0}}{({\omega }_0{Y}_{0}cos{\omega }_{0}t)}^2(1-{\Delta }_1{(Y_{0}sin{\omega }_{0}t)}^2)dt-{\int }_0^{\frac{2\pi }{{\omega }_0}}{({\omega }_0{Y}_{0}cos{\omega }_{0}t)}^2({\Delta }_2{Y}_{0}sin{\omega }_{0}t)dt\hfill \\ =\pi {\omega }_0{Y}_0^2(1-\frac{1}{4}{\Delta }_1{Y}_0^2)-0=\pi {\omega }_0{Y}_0^2(1-\frac{1}{4}{\Delta }_1{Y}_0^2)\hfill \end{array}$$

(13)

When $\xi \ne 0$ and smaller, the Poincaré mapping will have a hyperbolic fixed point, and the homoclinic orbital of the Hamilton system (11) will break, so as to generate stable invariant manifold $W_1$ and unstable invariant manifold $W_2$ of this fixed point. Melnikovs’ theory explains that whether $W_1$ and $W_2$ intersect determines the existence of chaos in the system. Hence, when $M\left(Y_0\right)=0$, that is $Y_0=0$ or $Y_0=\frac{2}{\sqrt{{\Delta }_1}}$, there is a sufficiently small $\xi$, so that stable manifolds $W_1$ and unstable manifolds $W_2$ intersect. Further, $\frac{dM\left(Y_0\right)}{dt}=\frac{-4\pi {\omega }_0}{\sqrt{{\delta }_1}}<0\ne 0$, there must be a transversal homoclinic point on the Poincaré section.

3.2.2. Numerical Simulation

For the above, assume that $\xi$ is extremely close to 0, that is $\xi =0$, given ${\Delta }_1=16$, ${\omega }_0=10$, and obtain $Y_0=\frac{2}{\sqrt{{\Delta }_1}}=0.5$, for the Equation (10), the starting point of the exact solution is $(0,{\omega }_0{Y}_0)=(0,5)$, then the Hamiltonian is $H(x,y)=12.5$, and the phase diagram and time series under exact solution are shown in Figure 4.

Figure 4 shows the self-excited oscillation output under the initial condition $(0,5)$ of the exact solution. The response is sinusoidal waveform with a frequency approaching ${\omega }_0$. By solving the exact solution, the comparison between symplectic and non-symplectic methods are prepared.

3.3. Solutions of Numerical Scheme

3.3.1. Euler Scheme

There are many numerical solutions for differential equations, one of these is discretization. Firstly, the approximate value $y_0,y_1,\cdots y_{n-1},y_n,y_{n+1},\cdots$ of $y\left(t\right)$ on discrete node is obtained. Then, the recurrence formula of the approximate value $y_n$ of $y\left(t_n\right)$ is established to obtain the approximate value of $y\left(t\right)$ at each node. The spacing between two adjacent nodes is called step $h=t_{n+1}-t_n$, assuming that h is constant, then node $t_n=t_0+nh$, $n=0,1,2,\cdots ,n$. Then, defined by the derivative, Euler discretization of the equation, have

$$\frac{y\left(t_{n+1}\right)-y\left(t_n\right)}{h}=\frac{dy\left(t_n\right)}{dt}=f(t_n,y\left(t_n\right))$$

(14)

Through it, we can re describe as

$$y_{n+1}=y_n+hf(t_n,y_n)$$

(15)

According to a memristor circuit based on its van der Pol oscillator without external excitation, the Equation (6) is discretized and obtained

$$\left\{\begin{array}{c}{x}_{n+1}=x_n+h{f}_1(t_n,x_n,y_n)\\ y_{n+1}=y_n+h{f}_2(t_n,x_n,y_n)\end{array}\right.$$

(16)

where, $f_1(t,x,y)=\dot{x}=\frac{dx}{dt}$, $f_2(t,x,y)=\dot{y}=\frac{dy}{dt}$. Hence, the iteration format is

$$\left\{\begin{array}{c}{x}_{n+1}=x_n+h{y}_n\hfill \\ y_{n+1}=y_n+h\left(\xi (1-{\Delta }_1{x}_n^2-{\Delta }_2{x}_n)y_n-{\omega }_0^2{x}_n\right)\hfill \end{array}\right.$$

(17)

3.3.2. Symplectic Euler Scheme

The algebraic relation obtained by Euler method does not depend on initial and boundary conditions. However, no matter how small the step size is, the vibration will be unbounded, and with the evolution of time, the vibration amplitude will tend to infinity, destroying the periodicity of the system. Therefore, the symplectic Euler scheme is constructed to improve the structure-preserving performance of the traditional Euler scheme, and the symplectic Euler discrete scheme is

$$y_{n+1}=y_n+hf(t_{n+1},y_n)$$

(18)

According to memristor circuit based on its van der Pol oscillator without external excitation, the Equation (6) is discretized and obtained

$$\left\{\begin{array}{c}{x}_{n+1}=x_n+h{f}_1(t_n,x_{n+1},y_n)\\ y_{n+1}=y_n+h{f}_2(t_n,x_{n+1},y_n)\end{array}\right.$$

(19)

Hence, the iteration format of memristor circuit based on its van der Pol oscillator without external excitation is

$$\left\{\begin{array}{c}{x}_{n+1}=x_n+h{y}_n\hfill \\ y_{n+1}=y_n+h\left(\xi (1-{\Delta }_1{x}_{n+1}^2-{\Delta }_2{x}_{n+1})y_n-{\omega }_0^2{x}_{n+1}\right)\hfill \end{array}\right.$$

(20)

It is easy to obtain that the Jacobian matrix of Equation (20) is symplectic.

3.3.3. Four-Order Runge–Kutta Scheme

The four-order Runge–Kutta method is more widely used than Euler method in practical applications, Runge–Kutta method is a favorable tool for numerical solution of ordinary differential nonlinear equations, with high calculation accuracy. By shortening the step distance and increasing the order, the error range can be further controlled, and better results can be obtained in the case of calculation convergence. The classical method is as follows:

$$\left\{\begin{array}{c}\begin{array}{c}{\mathit{y}}_{\mathit{n}+1}={\mathit{y}}_{\mathit{n}}+\frac{h}{6}({\mathit{K}}_1+2{\mathit{K}}_2+2{\mathit{K}}_3+{\mathit{K}}_4)\hfill \\ {\mathit{K}}_1=f(x_n,{\mathit{y}}_{\mathit{n}})\hfill \\ {\mathit{K}}_2=f(x_n+\frac{h}{2},{\mathit{y}}_{\mathit{n}}+\frac{h}{2}{\mathit{K}}_1)\hfill \\ {\mathit{K}}_3=f(x_n+\frac{h}{2},{\mathit{y}}_{\mathit{n}}+\frac{h}{2}{\mathit{K}}_2)\hfill \\ {\mathit{K}}_4=f(x_n+h,{\mathit{y}}_{\mathit{n}}+h{\mathit{K}}_3)\hfill \end{array}\end{array}\right.$$

(21)

where, h is step, and let $y_n=y\left(t_n\right)$, then the expression of ${\mathit{y}}_{n+1}$ is exactly identical to the first four terms of Taylor expansion of $\mathit{y}\left(x_{n+1}\right)$ at $x_{n+1}$, and the local truncation error is $O\left(h^4\right)$.

3.3.4. Four-Order Symplectic Runge–Kutta–Nyström Scheme

For the solution of system (8), there is the following four-order symplectic Runge–Kutta–Nyström scheme as

$$\left\{\begin{array}{c}{g}_i=x_n+c_{i}h{y}_n+h^2{\displaystyle \sum _{j=1}^s}{a}_{ij}f\left(g_j\right),i\in (1,s)\hfill \\ x_{n+1}=x_n+h{y}_n+h^2{\displaystyle \sum _{j=1}^s}{d}_{j}f\left(g_j\right)\hfill \\ y_{n+1}=y_n+h{\displaystyle \sum _{j=1}^s}{b}_{j}f\left(g_j\right)\hfill \end{array}\right.$$

(22)

where, $b_j=d_j\left(1-c_j\right),1⩽j⩽s$ and $d_i{d}_j\left(c_j-c_i\right)=d_j{a}_{ji},1⩽i<j⩽s$, and this scheme is symplectic.

The coefficient of Equation (22) are

$$\left\{\begin{array}{c}A=\left(\begin{array}{cc}\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}& \begin{array}{cc}{a}_{13}& a_{14}\\ a_{23}& a_{24}\end{array}\\ \begin{array}{cc}{a}_{31}& a_{32}\\ a_{41}& a_{42}\end{array}& \begin{array}{cc}{a}_{33}& a_{34}\\ a_{43}& a_{44}\end{array}\end{array}\right)=\left(\begin{array}{cccc}0& 0& 0& 0\\ \frac{5}{12}-\frac{\alpha }{2}& 0& 0& 0\\ 0& -\frac{1}{2}+\alpha & 0& 0\\ 0& -\frac{1}{2}+\alpha & 0& 0\end{array}\right),\hfill \\ B=\left(\begin{array}{cc}{b}_1& \begin{array}{cc}{b}_2& \begin{array}{cc}{b}_3& b_4\end{array}\end{array}\end{array}\right)=\left(\begin{array}{cc}\frac{3}{4}-\alpha & \begin{array}{ccc}\frac{1}{2}& -\frac{1}{4}+\alpha -b_4& b_4\end{array}\end{array}\right),\hfill \\ C=\left(\begin{array}{cc}{c}_1& \begin{array}{cc}{c}_2& \begin{array}{cc}{c}_3& c_4\end{array}\end{array}\end{array}\right)=\left(\begin{array}{cc}\alpha & \begin{array}{cc}1-\alpha & \begin{array}{cc}\alpha & \alpha \end{array}\end{array}\end{array}\right),\hfill \\ D=\left(\begin{array}{cc}{d}_1& \begin{array}{cc}{d}_2& \begin{array}{cc}{d}_3& d_4\end{array}\end{array}\end{array}\right)\hfill \\ =\left(\begin{array}{cc}\frac{7-9\alpha }{12}& \begin{array}{cc}\frac{\alpha }{2}& \begin{array}{cc}(1-\alpha )(\alpha -\frac{1}{4}-b_4)& (1-\alpha )b_4\end{array}\end{array}\end{array}\right).\hfill \end{array}\right.$$

where, $b_4$ is arbitrary parameter and $6{\alpha }^2-6\alpha +1=0$.

3.4. Numerical Simulation

3.4.1. Euler Method and Symplectic Euler Method

a.$\xi =0.1$

Here, according to Section 3.2.2, also set the starting point to $(0,5)$, given $\xi =0.1$, ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, and set the step $h=0.001$, the number of iterations is n = 10,000. Then, the comparison of Euler and symplectic Euler phase trajectories and displacement component changes of time series were obtained in Figure 5.

From Figure 5, it can be found that the non-symplectic method leads to divergence of phase trajectory and time displacement components, it does not preserve the structure, while the symplectic method has a good effect of symmetry and structure preserving. Additionally, the comparison diagram of the exact solution and symplectic Euler method as shown in Figure 6.

Through the comparison, the symplectic method is shown to have a higher coincidence with the exact solution and the structure is preserved.

According to Section 3.2.2, we know when $\xi =0$, the Hamiltonian of the system (6) is $H(x,y)=12.5$. Then, we calculate the Hamiltonian under the symplectic method and the non-symplectic method when $\xi$ is small enough, as shown in

Table 2. Exact solution, symplectic and non-symplectic Hamiltonian of $\xi =0.1$.

Iterations	Euler	Symplectic Euler
1	12.5000	12.5000
2	12.5038	12.5011
3	12.5074	12.5022
⋮	⋮	⋮
9999	22.1613	12.5596
10,000	22.1693	12.5627

.

Therefore, a comparison diagram of error trends of different methods is presented, as shown in Figure 7.

When $\xi$ is sufficiently small, the Hamiltonian of the system (6) remains around $H(x,y)=12.5$ by symplectic Euler method, however, the classic Euler method diverges and does not preserve the symmetry and structure.

b.$\xi =20$

In order to further illustrate the superiority of the symplectic method, we set $\xi =20$ and adjust the step h based on this, which is described in two cases.

Case 1:$h=0.001$.

When the value of $\xi$ is very large, given $\xi =20$, likewise ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, the number of iterations n = 10,000, set the step $h=0.001$, Figure 8 shows the phase plane and time traces between the Euler method and the symplectic Euler method.

Obviously, from Figure 8, symplectic method is more reliable in structure preserving performance and convergent faster than non-symplectic method.

Case 2:$h=0.01$.

As in case 1, we only increase the step h, and let $h=0.01$, other conditions remain unchanged, resulting in Figure 9.

According to Figure 9, when the step increased, the trajectory trend obtained by Euler method will be quickly decoupled from the periodic one and instantaneous collapsed, but the stability properties of the periodic orbit are preserved. by the symplectic Euler method.

3.4.2. Runge–Kutta Method and Symplectic Runge–Kutta–Nyström Method

a.$\xi =0.1$

As in Section 3.3.2, when $\xi =0.1$, we analyzed the classic four-order Runge–Kutta method and four-order symplectic Runge–Kutta–Nyström method.

Case 1:$h=0.01$.

Here, set the starting point to $(0,5)$, given $\xi =0.1$, ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, and set the step $h=0.01$, the number of iterations is n = 10,000. Then, the comparison of classic four-order Runge–Kutta method and four-order symplectic Runge–Kutta–Nyström phase trajectories and displacement component changes of time series were obtained in Figure 10.

From Figure 10, it can be found that both methods preserved the structure and symmetry of system well, but the computation time of symplectic method is faster than that of non-symplectic method. For further analysis, we increased the step.

Case 2:$h=0.09$.

As in case 1, we only increase the step h, let $h=0.09$, the number of iterations is n = 10,000, and other conditions remain unchanged, resulting in Figure 11.

According to Figure 11, when the step increased, the trajectory trend obtained by Runge–Kutta method will be quickly decoupled from the periodic one and instantaneously collapse, and the structure cannot be preserved. However, the stability properties of the periodic orbit and symmetry are preserved by the symplectic Runge–Kutta–Nyström method. Then, we calculate the Hamiltonian under classic four-order Runge–Kutta method and four-order symplectic Runge–Kutta–Nyström method when $\xi$ is small enough, as shown in

Table 3. Exact solution, symplectic and non-symplectic Hamiltonian of $\xi =0.1$ and $h=0.09$.

Iterations	RK4	SRKN4
1	12.5000	12.5000
2	12.2044	12.3024
3	12.0098	12.1816
⋮	⋮	⋮
9999	3.3516	12.5290
10,000	3.3554	12.7716

.

Therefore, a comparison diagram of error trends of different methods is presented, as shown in Figure 12.

When $\xi$ is sufficiently small, the Hamiltonian of symplectic Runge–Kutta method still remain within a certain range, however, the classic Runge–Kutta method diverges and does not preserve the structure. Through the comparison, the symplectic method is shown to have a higher coincidence with the exact solution and the structure is preserved.

b.$\xi =20$

In order to further illustrate the superiority of the symplectic method, we set $\xi =20$ and adjust the step h based on this, also consider the description of two cases.

Case 1:$h=0.001$.

When the value of $\xi$ is very large, given $\xi =20$, likewise ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, the number of iterations n = 10,000, set the step $h=0.001$, Figure 13 shows the phase plane and time traces between the classic four-order Runge–Kutta method and four-order symplectic Runge–Kutta–Nyström method.

Obviously, from Figure 13, symplectic method is more reliable in structure preserving performance and convergent faster than non-symplectic method.

Case 2:$h=0.02$.

As in case 1, we only increase the step h, and let $h=0.02$, other conditions remain unchanged, resulting in Figure 14.

According to Figure 14, when the step increased, the trajectory trend obtained by classic Runge–Kutta method will be quickly decoupled from the periodic one and instantaneous collapsed.

Through case 1 and case 2 of Section 3.4.1 and Section 3.4.2, under $\xi =0.1$ and $\xi =20$, in the case of long-time and large step numerical integration, symplectic method can still better maintain the characteristics of the system than non-symplectic method. That is, the symplectic method can preserve the structure well. In particular, when the perturbation term is small, the system tends to be a Hamiltonian system, and the symplectic method can maintain the symmetry of the system.

3.4.3. Symplectic Euler Method and Four-Order Runge–Kutta Method

In the previous sections, we compare the symplectic and classical methods of Euler method and Runge–Kutta method, respectively. In this section, we consider comparing the symplectic Euler method with the classical four-order Runge–Kutta method.

a.$\xi =0.1$

As in Section 3.3.2, when $\xi =0.1$, we analyzed the classic four-order Runge–Kutta method and symplectic Euler method.

Case 1:$h=0.001$.

Set the starting point to $(0,5)$, given $\xi =0.1$, ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, and set the step $h=0.001$, the number of iterations is n = 10,000. Then, the comparison of classic four-order Runge–Kutta method and symplectic Euler phase trajectories and displacement component changes of time series were obtained in Figure 15.

From Figure 15, it can be found that both methods preserved the structure and symmetry of the system well, but the computation time of symplectic Euler method is faster than that of classic Runge–Kutta method. For further analysis, we increased the step.

Case 2:$h=0.1$.

As in case 1, we only increase the step h, let $h=0.1$, the number of iterations is n = 10,000 and other conditions remain unchanged, resulting in Figure 16.

According to Figure 16, when the step increases, the trajectory trend obtained by the classical Runge–Kutta method will be decoupled from the periodic trend, resulting in instantaneous collapse and the structure cannot be preserved. Meanwhile, the symplectic Euler method did not structure preserved well in the end, but it did not collapse, nor did it collapse in time. Then, we calculate the Hamiltonian under classic four-order Runge–Kutta method and symplectic Euler method when $\xi$ is small enough, as shown in

Table 4. Exact solution, symplectic and non-symplectic Hamiltonian of $\xi =0.1$ and $h=0.1$.

Iterations	RK4	Symplectic Euler
1	12.5000	12.5000
2	12.0852	12.5613
3	11.8489	23.2033
⋮	⋮	⋮
9999	$5.7256\times {10}^{-11}$	8.1131
10,000	$5.6771\times {10}^{-11}$	18.6570

.

Then, comparison diagram of error trends of different methods is presented, resulting in Figure 17.

The $\xi$ is sufficiently small, the Hamiltonian of the classic Runge–Kutta method diverges and does not preserve the structure, and the errors produced by symplectic Euler method are also larger, but compared with Runge–Kutta method, it is still remain within a certain range.

b.$\xi =20$, $h=0.001$

In order to further illustrate the contrast between symplectic Euler and classical Runge–Kutta method. When the value of $\xi$ is very large, given $\xi =20$, likewise ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, the number of iterations n = 10,000, set the step $h=0.001$, Figure 18 shows the phase plane and time traces between the classic four-order Runge–Kutta method and symplectic Euler method.

From Figure 18, symplectic Euler method is more reliable in structure preserving performance and convergent faster than classic four-order Runge–Kutta method. Furthermore, In this case, we did not compare the methods when the steps increased, so it was unnecessary to preserve the structure of the system. Moreover, the above method has obvious advantages in speed and has also achieved victory in structure preserving.

3.4.4. Symplectic Euler Method and Symplectic Runge–Kutta–Nyström Method

a.$\xi =0.1$

When $\xi =0.1$, we analyzed the symplectic Euler method and four-order symplectic Runge–Kutta–Nyström method.

Case 1:$h=0.01$.

Similarly, set the starting point to $(0,5)$, given $\xi =0.1$, ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, and set the step $h=0.01$, the number of iterations is n = 10,000. Then, the comparison of the symplectic Euler and four-order symplectic Runge–Kutta–Nyström phase trajectories and displacement component changes of time series were obtained in Figure 19.

From Figure 19, it can be found that both methods preserved the structure and symmetry of system well, but the computation time of symplectic Euler method is faster than that of symplectic Runge–Kutta–Nyström method. For further analysis, we increased the step.

Case 2:$h=0.1$.

As in case 1, we only increase the step h, let $h=0.1$, the number of iterations is n = 10,000, and other conditions remain unchanged, resulting in Figure 20.

According to Figure 20, when the step increases, the trajectory trend of the symplectic Euler method is decoupled from the periodic one, and the structure cannot be preserved. However, symplectic Runge–Kutta–Nyström preserves symmetry and structure, while the symplectic Euler method has an advantage in calculating speed. Then, we calculate the Hamiltonian under symplectic Euler method and four-order symplectic Runge–Kutta–Nyström method when $\xi$ is small enough, as shown in

Table 5. Symplectic Euler and four-order symplectic Runge–Kutta–Nyström methods’ Hamiltonian of $\xi =0.1$ and $h=0.1$.

Iterations	Symplectic Euler	SRKN4
1	12.5000	12.5000
2	12.5613	12.2914
3	23.2033	12.2471
⋮	⋮	⋮
9999	8.1131	12.7378
10,000	18.6570	12.5994

.

Therefore, a comparison diagram of error trends of different methods is presented, as shown in Figure 21.

From above, when $\xi$ is sufficiently small and step increases, the symplectic Runge–Kutta–Nyström method’s Hamiltonian remains within a certain range, while the symplectic Euler method diverges, but remains within a certain range (albeit a large one). By comparison, the symmetry and structure of the symplectic Runge–Kutta–Nyström method are better under the condition of small perturbation and large step.

b.$\xi =20$

To further illustrate the conditions of symplectic method and the characteristics of each methods, we set $\xi =20$ and adjust the step h based on this, also consider the description of two cases.

Case 1:$h=0.001$.

When the value of $\xi$ is large, given $\xi =20$, likewise ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, the number of iterations n = 10,000, set the step $h=0.001$, Figure 22 shows the phase plane and time traces between the symplectic Euler method and four-order symplectic Runge–Kutta–Nyström method.

From Figure 22, the symplectic Euler method and symplectic Runge–Kutta–Nyström method both preserve the structure of the system well, but the symplectic Euler method is superior in computing speed by comparison.

Case 2: Large step.

As in case 1, we only increase the step h, and make the step change constantly, other conditions remain unchanged, resulting in Figure 23.

According to Figure 23, it can be seen from the figure that when the step increases, $h=0.01$, the symplectic Euler method and the four-order symplectic Runge–Kutta–Nyström method do not well preserve the system structure. When the step continues to increase and $h=0.0159$, the trajectory trend obtained by the symplectic Euler method is decoupled from the periodic trend first, and the system structure is no longer preserved. Meanwhile, if the step size continues to increase, both methods can no longer preserve the system structure.

Through the discussion of various cases in Section 3.4, we conclude that the symplectic method maintains the characteristics of the system better than the non-symplectic method in the case of large and long numerical integrals over long time. In particular, when the perturbation is very small, the system tends to be Hamiltonian, and symplectic methods can preserve the symmetry and structure of the system better than classic methods. Meanwhile, when the method is not confined to the comparison of the two symplectic methods, the perturbation is very small, the symplectic methods do not differ much. When the perturbation becomes larger, the step becomes the key to influence the numerical integration, that is, when the step is larger, the symplectic Runge–Kutta–Nyström method wins, but the computational speed of the symplectic Euler method is better than that of the symplectic Runge–Kutta–Nyström method.

4. Primary and Subharmonic Simultaneous Resonance of Forced van der Pol Oscillator

Due to the interaction between self-excited oscillation and forced oscillation, complex dynamic phenomena appear in the van der Pol non-autonomous system with a sinusoidal forcing term. Many resonance phenomena occur in non-autonomous dynamical systems under forced excitation, especially under multi-frequency excitation. In this section, the primary and subharmonic simultaneous resonances of van der Pol oscillators are studied.

4.1. First-Order Approximate Solution of Primary and Subharmonic Simultaneous Resonance

A memristor circuit based on forced van der Pol oscillator system excited by multiple frequencies can be described as system (4). The parameters are specified as follows: ${\omega }_1={\omega }_0+\xi {\psi }_1$, ${\omega }_2=3{\omega }_0+\xi {\psi }_2$, $A_1=\xi f$, $f=O\left(1\right)$, ${\psi }_1=O\left(1\right)$, ${\psi }_2=O\left(1\right)$, the system (4) changes into

$$\ddot{x}+{\omega }_0^{2}x=\xi \left((1-{\Delta }_1{x}^2-{\Delta }_{2}x)\dot{x}+fcos{\omega }_{1}t\right)+A_{2}cos{\omega }_{2}t$$

(23)

where, ${\omega }_0$ is the natural frequency of the linearized system (2). A multi-scale method [59] is applied to study the first-order approximate solution of the system. Suppose that the solution form of the system (23) is as follows:

$$x\left(t\right)=x_0(T_0,T_1)+\xi x_1(T_0,T_1)$$

(24)

where, time scales are $T_0=t$, $T_1=\xi t$, and substituting Equation (24) into Equation (23), by comparing the same power of, a system of partial differential equations is obtained:

$$D_0^2{x}_0+{\omega }_0^2{x}_0=A_{2}cos(3{\omega }_0{T}_0+{\psi }_2{T}_1)$$

(25)

$$D_0^2{x}_1+{\omega }_0^2{x}_1=-2D_0{D}_1{x}_0+fcos({\omega }_0{T}_0+{\psi }_1{T}_1)+2D_0(1-{\Delta }_1{x}_0^2-{\Delta }_2{x}_0)x_0$$

(26)

Then, the general solution of Equation (25) is $x_0(T_0,T_1)=a\left(T_1\right)cos\left({\omega }_0{T}_0+b\left(T_1\right)\right)+\frac{A_2}{{\omega }_0^2-{\omega }_1^2}cos(3{\omega }_0{T}_0+{\psi }_2{T}_1)$ and it can also be written as plural form:

$$x_0(T_0,T_1)=K\left(T_1\right)e^{j{\omega }_0{T}_0}+P{e}^{j(3{\omega }_0{T}_0+{\psi }_2{T}_1)}+c$$

(27)

where, $K\left(T_1\right)=\frac{a\left(T_1\right)e^{jb\left(T_1\right)}}{2}$, $P=\frac{A_2}{2({\omega }_0^2-{\omega }_1^2)}$, c is the conjugate of all the preceding terms. $a\left(T_1\right)$, $b\left(T_1\right)$ are slowly varying amplitude and phase, respectively. Substitute Equation (27) into Equation (26), eliminate the secular term and separating real and imaginary parts. The differential equations satisfied by slowly varying amplitude $a\left(T_1\right)$ and phase $b\left(T_1\right)$ are obtained

$$\left\{\begin{array}{c}{D}_{1}a=\frac{\xi {\psi }_1{a}^{2}B}{24{\omega }_0}sin(3b-{\psi }_2{T}_1)-\frac{f}{2{\omega }_0}sin(b-{\psi }_1{T}_1)-\frac{\xi {\psi }_{1}a}{6}\hfill \\ a{D}_{1}b=\frac{\xi a}{6{\omega }_0}+\frac{\xi {\Delta }_{1}a{P}^2}{12{\omega }_0}+\frac{\xi {\Delta }_1{a}^{2}P}{6{\omega }_0}cos(3b-{\psi }_2{T}_1)-\frac{f}{2{\omega }_0}cos(b-{\psi }_1{T}_1)+\frac{\xi {\Delta }_1{a}^3}{6{\omega }_0}\hfill \end{array}\right.$$

(28)

Consequently, the first-order approximate solution of system (23) can be expressed as

$$x_0\left(t\right)=2Pcos\left({\omega }_{2}t\right)+acos({\omega }_{0}t+b)$$

(29)

where, a and b determined by the Equation (28).

4.2. Steady Solution and Its Stability Conditions

From above, the necessary conditions for the existence of steady solutions are $b-{\psi }_1{T}_1$ and $3b-{\psi }_2{T}_1$ are constants from Equation (27), where ${\psi }_2=3{\psi }_1$, ${\omega }_2=3{\omega }_1$. In other words, the steady solution of primary and subharmonic simultaneous resonance can be obtained only when the two excitation frequencies are sufficient to a specific multiple relationship. Let $3\psi =3{\psi }_1={\psi }_2$, $b-{\psi }_1{T}_1=\gamma$, then Equation (28) can be written as autonomous differential equations:

$$\left\{\begin{array}{c}{D}_{1}a=\frac{\xi {\psi }_1{a}^{2}P}{24{\omega }_0}sin\left(3\gamma \right)-\frac{f}{2{\omega }_0}sin\left(\gamma \right)-\frac{\xi {\psi }_{1}a}{6}\hfill \\ a{D}_{1}b=\frac{\xi a}{6{\omega }_0}+\frac{\xi {\Delta }_{1}a{P}^2}{12{\omega }_0}+\frac{\xi {\Delta }_1{a}^{2}P}{6{\omega }_0}cos\left(3\gamma \right)-\frac{f}{2{\omega }_0}cos\left(\gamma \right)+\frac{\xi {\Delta }_1{a}^3}{6{\omega }_0}\hfill \end{array}\right.$$

(30)

The corresponding first-order approximate solution becomes

$$x\left(t\right)=acos({\omega }_{1}t+\gamma )+2Pcos\left({\omega }_{2}t\right)$$

(31)

where, ${\omega }_2=3{\omega }_1$. To verify the accuracy of the first-order approximate solution of Equation (31). The Runge–Kutta method with variable step length is used to calculate the system (4), and it is compared with the first-order approximate solution of Equation (31). According to the above section, let $\xi =0.1$, ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=2.5$, ${\omega }_2=3{\omega }_1=7.5$, $A_1=0.2$, $A_2=2$ and the simulation duration $t=50$. The transient response is called when the initial time $t\in (0,20)$. Figure 24 shows the transient responses times.

Then, let the last time $t\in (20,50)$, and the response amplitude as the steady-state amplitude, take the initial value $(a_0,{\gamma }_0)=(0.72,-1.51)$, the steady-state response is shown in Figure 25.

As can be seen from Figure 24 and Figure 25, when ${\omega }_2=3{\omega }_1$, the approximation between numerical solution and analytical solution is better, and the multiple solutions and unstability phenomena of the system are mainly dominated by $acos({\omega }_{1}t+\gamma )$ of system (31). Therefore, in the subsequent study of the dynamical behaviors of system (4), only ${\omega }_1$ in the external excitation needs to be studied.

4.3. Analysis of Chaotic Dynamics

In the resonance analysis in Section 4.1 and Section 4.2, we get ${\omega }_2=3{\omega }_1$, then the system (4) becomes

$$\left\{\begin{array}{c}\dot{x}=y\hfill \\ \dot{y}=\xi (1-{\Delta }_1{x}^2-{\Delta }_{2}x)y-{\omega }_0^{2}x+A_{1}cos\left({\omega }_{1}t\right)+A_{2}cos\left(3{\omega }_{1}t\right)\hfill \end{array}\right.$$

(32)

Then, the system (32) will have different and complex dynamical behaviors when changing $\xi$, and we will analyze these dynamical behaviors under different $\xi$.

4.3.1. Dynamical Behaviors of $\xi =0.1$

The 2D nonlinear non-autonomous systems can produce chaotic behaviors, and chaos will occur when the nonlinearity of a memristor circuit based on its van der Pol oscillators is strong. When $\xi =0.1$, with ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, $A_1=12$, $A_2=8$ and ${\omega }_1=7.1$, the initial value $\left(x,y\right)=\left(0,0.5\right)$ is selected. The chaotic attractor and time series are shown in Figure 26.

Under the same parameters, when control parameters ${\omega }_1\in \left(6,14\right)$, the bifurcation diagram of the system (32) is shown in Figure 27.

Figure 27 shows a bifurcation diagram of x as a function of ${\omega }_1$, and the system undergoes a variety of dynamical behaviors. In addition, to further verify the existence of chaos, we give the Lyapunov exponential spectrum and the largest Lyapunov exponential trend diagram in Figure 28.

The Lyapunov exponential spectrum in Figure 28 shows rich dynamical behaviors and shows that the largest Lyapunov exponential is positive. The largest Lyapunov exponential diagram has the same results as the bifurcation diagram. These results also agree with the analysis of van der Pol oscillators.

4.3.2. Dynamical Behaviors of $\xi =5$

Similarly, when $\xi =5$, with ${\Delta }_1=16$, ${\Delta }_2=8$, ${\omega }_0=10$, $A_1=10$, $A_2=5$ and ${\omega }_1=3.7$, the initial value $\left(x,y\right)=\left(0,0.5\right)$ is selected. The chaotic attractor and zoom in on Poincaré section are shown in Figure 29.

From Figure 29, the chaotic attractor generates a type relaxation oscillator configuration, and the continuously dense point shows chaotic behavior. Then, we change the value of ${\omega }_1$ to get the phase diagram and Poincaré section with different parameters and the states of periodic 1, periodic 2 and periodic3 can be obtained, as shown in Figure 30.

When control parameters ${\omega }_1\in \left(0,9\right)$, the bifurcation diagram of the system (32) is shown in Figure 31.

From above, it can be seen that many narrow periodic windows are embedded in chaotic regions, and it contains rich dynamical behaviors. Thus, the forced parameter $\xi$ and ${\omega }_1$ can be used to tune the memristor circuit to different chaotic or periodic oscillations.

5. Conclusions

A memristor circuit based on its van der Pol oscillator is introduced, and the stability and dynamic characteristics of memristor circuit are studied. Meanwhile, the differences of the van der Pol oscillation model of memristor between the Euler method and symplectic Euler method, four-order Runge–Kutta method and four-order symplectic Runge–Kutta–Nyström method, symplectic Euler method and four-order Runge–Kutta method, symplectic Euler method and four-order symplectic Runge–Kutta–Nyström method in preserving symmetry and structure are compared from theoretical and numerical simulation: The symplectic scheme can better simulate the oscillate process of the oscillator and keep the conservation of the harmonic oscillator accurately in the time domain for a long time. Secondly, the analytic solution of the primary and subharmonic (1/3 order) simultaneous resonance of the van der Pol oscillation system is obtained by using the multi-scale method, and the stability condition of the steady solution is obtained by using Lyapunov stability theory. Finally, the dynamical behaviors of the system is studied from the resonance relationship, which is helpful to understand the oscillation characteristics of other nonlinear oscillation problems in the form of van der Pol equation later. In the future, we will consider comparing more structure preserving symplectic methods with classical integral methods to demonstrate the advantages of symplectic methods, and which will be used in aerospace, nanomaterials and other applications. Meanwhile, we will further study the memristor system on the electronic circuit, and the global dynamics and infinite behaviors of the system are considered.