Fast-Slow Coupling Dynamics Behavior of the van der Pol-Rayleigh System

Abstract

In this paper, the dynamic behavior of the van der Pol-Rayleigh system is studied by using the fast–slow analysis method and the transformation phase portrait method. Firstly, the stability and bifurcation behavior of the equilibrium point of the system are analyzed. We find that the system has no fold bifurcation, but has Hopf bifurcation. By calculating the first Lyapunov coefficient, the bifurcation direction and stability of the Hopf bifurcation are obtained. Moreover, the bifurcation diagram of the system with respect to the external excitation is drawn. Then, the fast subsystem is simulated numerically and analyzed with or without external excitation. Finally, the vibration behavior and its generation mechanism of the system in different modes are analyzed. The vibration mode of the system is affected by both the fast and slow varying processes. The mechanisms of different modes of vibration of the system are revealed by the transformation phase portrait method, because the system trajectory will encounter different types of attractors in the fast subsystem.

1. Introduction

Nonlinear oscillators appear in many applied sciences, involving mechanical engineering, mechanics, chemistry, biology and physics, and have a wide range of engineering backgrounds. The Duffing [1,2], van der Pol [3,4] and Rayleigh [5,6] oscillators are well-known models of nonlinear systems in physics and mechanics.

In the wake of the rapid development of the economy and technology, the continuous emergence of new structures and construction technologies and the increasing aesthetic needs of the public, pedestrian bridges, stadiums, floors and other structures are becoming increasingly slender. However, these structures are subjected to walking load, so it is easy to produce a large-amplitude vibration. The modified hybrid van der Pol-Rayleigh (MHVR) system proposed by Erlicher et al. [7] is a kind of typical self-excited oscillation system, which can be suitably applied to the pedestrian walking lateral force model.

As a combination of the two oscillators, the dynamic behavior of the van der Pol-Rayleigh system has attracted extensive attention. Tang et al. [8] studied the vibration response and its generation mechanism in the van der Pol-Rayleigh system under slow-varying periodic excitation, and analyzed the excitation hysteresis behavior and its generation mechanism of the system. Saha et al. [9] classified the van der Pol oscillator and Rayleigh oscillator into the general form of the Liénard-Levinson-Smith (LLS) system, and thus design a generalized van der Pol and Rayleigh oscillator family system with multiple limit cycles. Hasegawa [10] studied the Jarzynski equality in van der Pol and Rayleigh oscillators to which a ramp force with a duration τ is applied. In the τ range, the Jarzynski equality (JE) is basically valid, but not completely satisfied. The work distribution function (WDF) uses a u-shaped structure with large damping parameters. Jerzy [11] studied the interactions between two parametrically coupled, self-excited oscillators. The different types of motion were classified by Lyapunov’s exponent criterion. Veskos and Demiris [12] studied the effects of different lower-level building blocks of a robotic swinging system using the van der Pol and Rayleigh oscillators. The similarities and differences between these oscillators were analyzed. Chen et al. [13] studied the dynamics of a hybrid van der Pol-Rayleigh oscillator and discovered new dynamics that are different from van der Pol and Rayleigh oscillators. At the same time, the van der Pol-Rayleigh system is also widely used in mechanics, electronic circuits, biomechanics and other engineering areas. For example, Bartkowiak and Woernle [14] adjusted the parameters of the van der Pol-Rayleigh oscillator to the given harmonic base excitation, and then added them to the linear multi-body system. After adjusting the amplitude and phase shift, the vibration absorption and energy collection tasks were completed. Chen et al. [15] used the multi-scale method to decouple the van der Pol-Rayleigh equation and applied the discussion results to the dynamic model of the side-coupling system of a footbridge—a flexible footbridge.

The fast–slow coupling systems can often cause bursting oscillation, which is an important dynamic behavior. The fast–slow analysis method [16] first proposed by Rinzelis can be used to reveal the generation mechanism of the bursting response. The bursting oscillation [17,18,19,20,21] is composed of large and small oscillations alternately in a given period. Generally speaking, the large-amplitude oscillation is called the spiking state motion, and its trajectory moves in the large-amplitude limit cycle vector field. The small-amplitude oscillation is called the quiescent state motion and its trajectory moves within the small-amplitude limit cycle vector field or the equilibrium attraction region.

In this paper, we will analyze the dynamic response behavior of the van der Pol-Rayleigh system, reveal the mechanism of vibration generation in different modes and discuss the transition process among the different types of attractors in the process of rapid change. In practical engineering, the results of the discussion can be used to reasonably utilize or avoid these large vibrations in advance, so as to improve the stability of the equipment, prolong the service cycle and increase the economic and social benefits.

The remainder of this paper is organized as follows. In Section 2, the stability and bifurcation behavior of the equilibrium point of the system are analyzed. In Section 3, the fast subsystem is simulated numerically and analyzed with or without external excitation. In Section 4, the vibration behavior and its generation mechanism of the system in different modes are analyzed. In Section 5, we draw the conclusions.

2. The System Model and Bifurcation Analysis

2.1. The System Model

The model equation in [7] is:

$$\ddot{x}-\epsilon \dot{x}\left(1-\beta x^2-\alpha x\dot{x}-\mu {\dot{x}}^2\right)+{\omega }_0{}^{2}x=A\mathit{cos}(\omega t),$$

(1)

where $x\left(t\right)$ is a real function, ${\omega }_0$ is the natural frequency of the system, $\omega$ is the periodic excitation frequency, $A$ is the periodic excitation amplitude, $\epsilon$, $\beta$, $\mu$, $\alpha$ are parameters, and ${\omega }_0>0$, $\omega >0$, $\alpha >0$, $\beta >0$, $\mu >0$, $\epsilon \in R$. Obviously, the system is the Rayleigh oscillator when $\beta =\alpha =0$; the system is the van der Pol oscillator when $\alpha =\mu =0$. Therefore, this is a coupling van der Pol-Rayleigh system.

Let $f=A\mathit{cos}(\omega t)$, and Equation (1) can become:

$$\left\{\begin{array}{c}\dot{x}=y, \\ \dot{y}=\epsilon y\left(1-\beta x^2-\alpha xy-\mu y^2\right)-{\omega }_0{}^{2}x+f.\end{array}\right.$$

(2)

When $\omega \ll {\omega }_0$ and $\omega \le 1$, we can view $f=A\mathit{cos}(\omega t)$ as the slow variable. In any period $t\in \left[t_0, t_0+2\pi /\Omega \right]$ corresponding to the natural frequency, the excitation term $f=A\mathit{cos}(\omega t)$ changes between ${\Omega }_A=Acos\left(\omega t_0\right)$ and ${\Omega }_B=Acos\left(\omega t_0+2\pi \omega /\Omega \right)$, obviously ${\Omega }_A\approx {\Omega }_B$, so the change in $f$ in any period of the natural frequency is very small. Regarding $f$ as a generalized parameter of Equation (2), Equation (2) can be called a generalized autonomous system with respect to $f$ accordingly. We regard the generalized autonomous system (2) as the fast subsystem corresponding to the fast variable $x$ and $y$, while treating $f=A\mathit{cos}(\omega t)$ as the corresponding slow subsystem.

2.2. Bifurcation Analysis

To further reveal the mixed-mode vibration behavior of Equation (2), its bifurcation behavior with respect to the slow variable $f$ will be discussed.

2.2.1. Equilibrium Analysis and Fold Bifurcation

The equilibrium point of the system is $E_0\frac{f}{{\omega }_0{}^2},0$, and the Jacobian matrix of the linear system derived at $E_0$ is

$$J\epsilon =\left(\begin{array}{cc}0& 1\\ -{\omega }_0{}^2& \left(1-\beta x^2\right)\epsilon \end{array}\right)$$

(3)

The corresponding characteristic equation is

$${\lambda }^2+p\lambda +q=0,$$

(4)

where $p=\left(\beta x^2-1\right)\epsilon$, $q={\omega }_0{}^2$. Its eigenvalues are ${\lambda }_{1,2}=\frac{-p\pm \sqrt{p^2-4q}}{2}$.

When $p>0$, namely $f>\frac{1}{\sqrt{\beta }}{\omega }_0{}^2$ or $f<-\frac{1}{\sqrt{\beta }}{\omega }_0{}^2$, the equilibrium point is stable; conversely, when $p<0$, namely $-\frac{1}{\sqrt{\beta }}{\omega }_0{}^2<f<\frac{1}{\sqrt{\beta }}{\omega }_0{}^2$, the equilibrium point is unstable. At the same time, when $p>2\sqrt{q}$, namely $f<-{\omega }_0{}^2\sqrt{\frac{2{\omega }_0+\epsilon }{\epsilon \beta }}$ or $f>{\omega }_0{}^2\sqrt{\frac{2{\omega }_0+\epsilon }{\epsilon \beta }}$, the equilibrium point is a stable node; when $p<-2\sqrt{q}$, namely $-{\omega }_0{}^2\sqrt{\frac{\epsilon -2{\omega }_0}{\epsilon \beta }}<f<{\omega }_0{}^2\sqrt{\frac{\epsilon -2{\omega }_0}{\epsilon \beta }}$, the equilibrium point is an unstable node. When $0<p<2\sqrt{q}$, namely $-{\omega }_0{}^2\sqrt{\frac{2{\omega }_0+\epsilon }{\epsilon \beta }}<f<-\frac{1}{\sqrt{\beta }}{\omega }_0{}^2$ or $\frac{1}{\sqrt{\beta }}{\omega }_0{}^2<f<{\omega }_0{}^2\sqrt{\frac{2{\omega }_0+\epsilon }{\epsilon \beta }}$, the equilibrium point is a stable focus; when $-2\sqrt{q}<p<0$, namely $-\frac{1}{\sqrt{\beta }}{\omega }_0{}^2<f<-{\omega }_0{}^2\sqrt{\frac{\epsilon -2{\omega }_0}{\epsilon \beta }}$ or ${\omega }_0{}^2\sqrt{\frac{\epsilon -2{\omega }_0}{\epsilon \beta }}<f<\frac{1}{\sqrt{\beta }}{\omega }_0{}^2$, the equilibrium point is an unstable focus.

Because $q={\omega }_0{}^2>0$, the characteristic Equation (4) has no zero eigenvalue, and the system (2) has no fold bifurcation.

2.2.2. Hopf Bifurcation

Theorem 1.

If$f$is the bifurcation parameter of the fast subsystem (2), the sufficient and necessary conditions for the equilibrium point $E_0\frac{f}{{\omega }_0{}^2},0$ of the fast subsystem to produce Hopf bifurcation at$f=f_c$are

(1)$f_c=\pm \frac{1}{\sqrt{\beta }}{\omega }_0{}^2$;

(2)$\epsilon \ne 0$.

Proof of Theorem 1.

The condition of the characteristic Equation (4) has a pair of pure imaginary roots: $p=0$, $q>0$.

Because $q={\omega }_0{}^2>0$, we can obtain that $q>0$.

When $p=\left(\beta x^2-1\right)\epsilon =0$, we can determine that $f_c=\pm \frac{1}{\sqrt{\beta }}{\omega }_0{}^2$. Thus, Equation (2) may cause the Hopf bifurcation.

The existence of the Hopf bifurcation must also satisfy the transversal condition $Re{\left(\frac{d\lambda }{df}\right)}_{f=f_c}\ne 0$. The transversal condition of the Hopf bifurcation is proven below.

Taking the partial derivative of both sides of the characteristic Equation (4) with respect to $f$, we obtain

$$\frac{d\lambda }{df}=\frac{-2\epsilon \beta f\lambda }{\epsilon \beta f^2+{\omega }_0^4\left(2\lambda -\epsilon \right)}.$$

(5)

Let $\lambda =i{\omega }_1$, and substitute it into Equation (5) to obtain

$$\frac{d\lambda }{df}=\frac{-4\epsilon \beta f{\omega }_0^4{\omega }_1^2-2{\epsilon }^2\beta f{\omega }_1\left(\beta f^2-{\omega }_0^4\right)i}{{\epsilon }^2{\left(\beta f^2-{\omega }_0^4\right)}^2+4{\omega }_0^8{\omega }_1^2}.$$

(6)

Separate the real part of $\frac{d\lambda }{df}$, and substitute $f=\pm \frac{1}{\sqrt{\beta }}{\omega }_0{}^2$ into Equation (6) to obtain

$$Re{\left(\frac{d\lambda }{df}\right)}_{f=f_c}=\mp \frac{\epsilon \sqrt{\beta }}{{\omega }_0^2}.$$

(7)

Thus, when $\epsilon \ne 0$, determine $Re{\left(\frac{d\lambda }{df}\right)}_{f=f_c}\ne 0$. Therefore, Equation (2) will cause the Hopf bifurcation when $\epsilon \ne 0$. □

2.2.3. The First Lyapunov Coefficient

In order to judge the bifurcation direction and the stability of the Hopf bifurcation, the first Lyapunov coefficient $l_1$ [22] of the system is considered. Firstly, we will shift $E_0$ to the original point $O\left(0,0\right)$, and the transformation of coordinates is

$$\left\{\begin{array}{c}X=x-\frac{f}{{\omega }_0{}^2},\\ Y=y. \end{array}\right.$$

(8)

Substitute Equation (8) into Equation (2), and then use $x$ and $y$ to replace $X$ and $Y$, so Equation (2) can be converted into

$$\left\{\begin{array}{c}\dot{x}=y, \\ \dot{y}=-{\omega }_0{}^{2}x+\left(\epsilon -\frac{\epsilon \beta f^2}{{\omega }_0^4}\right)y+\epsilon y\left(-\beta x^2-\frac{2f\beta }{{\omega }_0^2}x-\mu y^2-\frac{\alpha f}{{\omega }_0^2}y-\alpha xy\right).\end{array}\right.$$

(9)

Next, transform Equation (9) into

$$\dot{X}=AX+F\left(X\right),X\in R^2,$$

(10)

where

$$X=\left(\begin{array}{c}x\\ y\end{array}\right),A=\left(\begin{array}{cc}0& 1\\ -{\omega }_0^2& \epsilon -\frac{\epsilon \beta f^2}{{\omega }_0^4}\end{array}\right),F\left(X\right)=\left(\begin{array}{c}0\\ \epsilon y\left(-\beta x^2-\frac{2f\beta }{{\omega }_0^2}x-\mu y^2-\frac{\alpha f}{{\omega }_0^2}y-\alpha xy\right)\end{array}\right).$$

(11)

Substitute $\lambda =\pm i{\omega }_1$ into characteristic Equation (5), and obtain ${\omega }_0^4-\beta f^2=0$, ${\omega }_1^2={\omega }_0^2$. Thus,

$$A=\left(\begin{array}{cc}0& 1\\ -{\omega }_1^2& 0\end{array}\right).$$

(12)

The eigenvector of $A$ corresponding to $i{\omega }_1$ is called $q$, and the eigenvector of $A^T$ corresponding to $-i{\omega }_1$ is called $p$. In other words, $Aq=i{\omega }_{1}q$, $A^{T}p=-i{\omega }_{1}p$, $〈p,q〉=1$. Thus, we can calculate the following:

$$q=\left(\begin{array}{c}1\\ i{\omega }_1\end{array}\right),p=\left(\begin{array}{c}\frac{1}{2}\\ -\frac{1}{2i{\omega }_1}\end{array}\right).$$

(13)

$F\left(x\right)=O\left({\left|\left|x\right|\right|}^2\right)$ is a smooth function that can be expanded as

$$F\left(x\right)=\frac{1}{2}B\left(x,x\right)+\frac{1}{6}C\left(x,x,x\right)+O\left({\left|\left|x\right|\right|}^4\right),$$

(14)

where $B\left(x,x\right)$ and $C\left(x,x,x\right)$ are bilinear functions, and

$$B_i\left(x,y\right)={\left.{\displaystyle \sum }_{j,k=1}^2\frac{{\partial }^2{F}_i\left(\xi \right)}{\partial {\xi }_j\partial {\xi }_k}\right|}_{\xi =0}{x}_j{y}_k,i=1,2.$$

(15)

$$C_i\left(x,y,z\right)={\left.{\displaystyle \sum }_{j,k,l=1}^2\frac{{\partial }^3{F}_i\left(\xi \right)}{\partial {\xi }_j\partial {\xi }_k\partial {\xi }_l}\right|}_{\xi =0}{x}_j{y}_k{z}_l,i=1,2.$$

(16)

After calculation, we can obtain

$$B\left(x,y\right)=\left(\begin{array}{c}0\\ -\frac{2\epsilon \beta f}{{\omega }_0^2}{x}_1{y}_2-\frac{2\epsilon \beta f}{{\omega }_0^2}{x}_2{y}_1-\frac{2\epsilon \alpha f}{{\omega }_0^2}{x}_2{y}_2\end{array}\right),$$

(17)

$$C\left(x,y,z\right)=\left(\begin{array}{c}0\\ -2\epsilon \beta \left(x_1{y}_1{z}_2+x_1{y}_2{z}_1+x_2{y}_1{z}_1\right)-2\epsilon \alpha \left(x_1{y}_2{z}_2+x_2{y}_1{z}_2+x_2{y}_2{z}_1\right)-6\epsilon \mu x_2{y}_2{z}_2\end{array}\right).$$

(18)

Then, after a series of calculations, we obtain

$$\begin{gathered} l_1=\frac{1}{2{\omega }_1}Re\left[〈p,C\left(q,q,\overline{q}\right)〉-2〈p,B\left(q,A^{-1}B\left(q,\overline{q}\right)\right)〉+〈p,B\left(\overline{q},{\left(2i{\omega }_{1}E-A\right)}^{-1}B\left(q,q\right)\right)〉\right] \\ =\frac{1}{2{\omega }_0}\left(2{\epsilon }^2\alpha -\epsilon \beta -3\epsilon \mu {\omega }_0^2\right). \end{gathered}$$

(19)

In [22], we can see that when $l_1>0$, the system experiences the subcritical Hopf bifurcation and an unstable limit cycle is obtained; when $l_1<0$, the system experiences the supercritical Hopf bifurcation and a stable limit cycle is obtained; when $l_1=0$, the system experiences the codimensional 2 degenerate Hopf bifurcation.

Now, fixing the parameters $\alpha =0.1$, $\beta =0.36$, $\mu =1.25$, ${\omega }_0=0.6$, we can obtain

$$l_1=\frac{1}{6}\left({\epsilon }^2-8.55\epsilon \right).$$

(20)

Thus, when $0<\epsilon <8.55$, $l_1<0$, the system experiences the supercritical Hopf bifurcation and the resulting limit cycle is stable. When $\epsilon <0$ and $\epsilon >8.55$, $l_1>0$, the system experiences the subcritical Hopf bifurcation and the resulting limit cycle is unstable.

2.3. Bifurcation Curve

Fixing the parameters $\alpha =0.1$, $\beta =0.36$, $\mu =1.25$, ${\omega }_0=0.6$, $\epsilon =1.5$, $\omega =0.001$, the bifurcation diagram of the fast subsystem (2) with respect to $f$ can be drawn, as shown in Figure 1. The red solid line indicates that the equilibrium is stable, the black solid line indicates that the equilibrium point is unstable, and the green solid line indicates that the system has a stable limit cycle; these are consistent with the actual characteristics of the van der Pol-Rayleigh system.

In Figure 1, $H{B}_1$ and $H{B}_2$ are Hopf bifurcation points; $E{P}_1$ and $E{P}_2$ are limit points, and also are the saddle of the system, called the turning point. There are four different equilibrium lines distributed on the equilibrium line $N_1{N}_2$. The first kind of equilibrium lines are from the point $N_1$ to the point $E{P}_1$ and from the point $E{P}_2$ to the point $N_2$, and the points on this kind of line are stable nodes. The second kind of equilibrium lines are from the point $E{P}_1$ to the point $H{B}_1$ and from the point $H{B}_2$ to the point $E{P}_2$, and the points on this kind of line are stable focuses. The third kind of equilibrium lines are from the point $H{B}_1$ to the point $F_1$ and from the point $F_2$ to the point $H{B}_2$, and the points on this kind of line are unstable focuses. The fourth kind of equilibrium line is from the point $F_1$ to the point $F_2$, and the points on this kind of line are unstable nodes.

The images drawn by the software XPPAUT are consistent with the theoretical analysis, and the type and stability of the system equilibrium point discussed by the linear characteristic equation derived from the system are consistent with the results obtained by the actual software application.

3. Numerical Simulation of the System

Fixing the parameters $\alpha =0.1$, $\beta =0.36$, $\mu =1.25$, ${\omega }_0=0.6$, the vibration behavior of the system before and after external excitation is given by numerical simulation.

When $A=0$, the system degenerates to the original MHVDR system. The numerical simulation shows that the system causes a simple periodic vibration, as shown in Figure 2. The phase plane portrait is a closed elliptical trajectory. Thus, it is also verified that the system displays periodic motion from another aspect. From Figure 2, we can see that the system trajectory hardly changes when $\epsilon$ changes.

When the slow-varying external excitation $f$ is applied to the system, the system presents a relatively complex vibration mode, which is not a single vibration mode within a period, as shown in Figure 3 and Figure 4. The vibration mode is composed of large and small oscillations alternately in a period.

Comparing (a) with (b) in Figure 3, the system trajectory hardly changes when $\epsilon$ changes. Comparing (b) with (c) in Figure 3, the period becomes short when $\omega$ changes. Finally, comparing (c) with (d) in Figure 3, when $A$ changes from $2$ to $0.8$, its vibration mode is still in the form of bursting oscillation, but it is more violent than the spiking state in Figure 3c. Therefore, it can be seen that the parameters have a certain influence on the bursting oscillation of the system.

4. Mechanism of Vibration Generation in Different Modes

Regarding the external excitation $Acos\left(\omega t\right)$ as the slow-varying bifurcation parameter of the fast subsystem, the transformation phase portrait of the fast variables $x$ and the slow variables $Acos\left(\omega t\right)$ is drawn. To further reveal the bifurcation mechanism of the bursting oscillation, the transformation phase portrait is superimposed with the bifurcation diagram of the fast subsystem about $f$. Thus, we can further study the vibration response behavior of the system visually, analyze the transfer process between multiple attractors in the system, discuss the mechanism of the bursting phenomenon generated and further explore the influence of external excitation amplitude $A$ on the fast subsystem.

Fixing the parameters $\alpha =0.1$, $\beta =0.36$, $\mu =1.25$, ${\omega }_0=0.6$, $\epsilon =1.5$, $\omega =0.001$, we take the excitation amplitudes $A$ as $0.2$, $0.6$, $1.8$, $12$ and $20,$ respectively. Thus, the time history diagrams and the superposition diagrams of the system are drawn.

4.1. The Mode of $A=0.2$

Letting $A=0.2$, we draw the time history diagram and the superposition diagram of the system, as shown in Figure 5.

In Figure 5a, the system presents a large vibration mode of periodic vibration with two frequency coupling. One of the frequencies is the external excitation frequency $\omega =0.001$, and the other frequency is close to ${\omega }_0=0.6$. The system vibrates at constant amplitude and high frequency in accordance with the natural frequency of the fast subsystem in each slow-varying period, which is presented as the spiking state motion.

Figure 5b is the superposition graph drawn, and $f=0.2cos\left(0.001t\right)$ changes within the interval $\left[-0.2,0.2\right]$. At this time, the system trajectory is only within the attraction domain of the limit cycle attractor, the system presents the spiking state motion of high frequency vibration, and the vibration amplitude is the same as that of the limit cycle.

4.2. The Mode of $A=0.6$

When $A$ increases, the vibration state of the system becomes different. When $A=0.6$, the time history diagram and the superposition diagram of the system are drawn, as shown in Figure 6.

In Figure 6a, we show that the system still presents the spiking state. However, in each slow-varying period, the amplitudes of the system are not the same, but there are distinct differences. Thus, it can be judged that the system changes from constant-amplitude, high-frequency vibration to variable-amplitude, high-frequency vibration.

It can be seen from Figure 6b that the system trajectory spreads towards both sides of the limit cycle. When $f$ is the Hopf bifurcation point $H{B}_2$, the system trajectory covers the whole region. Moreover, the amplitude of the system also changes. As the amplitude of the limit cycle decreases, the amplitude of the system is no longer the original high-frequency vibration amplitude.

4.3. The Mode of $A=1.8$

When $A$ continues to increase, the system can display the situation of the coexistence of the limit cycle attractor and stable equilibrium attractor, resulting in the “bistable state” situation. Figure 7 plots the time history diagram and the superposition diagram of the system when the excitation amplitude $A$ increases to $1.8$. At this point, the system has a burst response. Moreover, the system presents the superposition of two vibration modes of high-frequency, large-amplitude vibration and low-frequency vibration in each slow-varying period, namely the mixed vibration mode of the spiking state and the quiescent state alternately in a period.

In Figure 7b, we can see that the system has two Hopf bifurcation points $f_{H{B}_1}=-0.6$ and $f_{H{B}_2}=0.6$. When the system trajectory is within the attraction domain of the limit cycle attractor $\left[-0.6,0.6\right]$, it presents the spiking state mode of high-frequency vibration. When the system trajectory is within the attractor domain $\left[-1.8,-0.6\right]$ and $\left[0.6,1.8\right]$ of stable equilibrium, the system converges to the stable equilibrium line and presents the quiescent state mode. In one period, the system trajectory passes through four Hopf bifurcation points and switches back and forth between the limit cycle attractor and stable equilibrium attractor, resulting in two spiking states and two quiescent states.

4.4. The Mode of $A=12$

When $A$ continues to increase to $12$, the time history diagram and the superposition diagram of the system are drawn, as shown in Figure 8. From Figure 8a, another new trajectory emerges. At this time, the system trajectory is still affected by the limit cycle attractor and the stable equilibrium attractor, and the system has the tendency to move away from and close to the attractor. The spiking state completely disappears, and the system only presents the quiescent state mode, as shown in Figure 8b.

4.5. The Mode of $A=20$

When $A$ continues to increase to $20$, the time history diagram and the superposition diagram of the system are drawn, as shown in Figure 9. Figure 9b plots the superposition diagram, from which it can be seen that the system trajectory slowly moves around the limit cycle and the equilibrium point. As ${\omega }_0$ weakens, the attractivity of the attractor also weakens, and the influence of periodic excitation becomes more and more significant.

4.6. Summary

To sum up, both fast- and slow-varying processes affect the vibration mode of the system. When $A$ is small, the system has two kinds of frequencies and is more affected by the fast-varying process. When $A$ gradually increases, $f$ increases, and the influence of the slow-varying process is gradually reflected. At this time, the bursting phenomenon will occur. When $A$ reaches a certain value, the system is only affected by the slow-varying process and presents a simple periodic vibration mode.

In Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9b, we reveal the vibration mode mechanism of (a), and $f$ vibrates in $\left[-A, A\right]$. The system trajectory encounters different types of attractors in the fast subsystem, resulting in different vibration modes. When $A$ is small, the system trajectory presents the spiking state under the influence of stable limit cycle attractor. As $A$ increases gradually, the system trajectory will be affected by both the limit cycle attractor and stable equilibrium attractor, which shows a "bistable state" bursting situation. Moreover, the spiking state and the quiescent state alternately appear in the same period. Finally, the influence of periodic excitation is increasing, the attraction of the fast subsystem is becoming weaker and weaker, the system appears a single periodic motion, and the trajectory moves around the limit cycle and the equilibrium point.

5. Conclusions

In this paper, the dynamic behavior of the van der Pol-Rayleigh system with the double Hopf bifurcations generated by external excitation is studied by the fast–slow analysis method.

Firstly, the stability and bifurcation behavior of the equilibrium point of the system (2) are analyzed. By analyzing the characteristic equation, we know that the equilibrium point of the system has different properties under different conditions. In addition, the system has no fold bifurcation point, but may have a Hopf bifurcation point. By calculating the Hopf bifurcation’s transversality condition, we find that the Hopf bifurcation theorem holds and a Hopf bifurcation point exists in the system. By calculating the first Lyapunov coefficient, the bifurcation direction and stability of Hopf bifurcation are obtained.

Then, the fast subsystem is simulated numerically and analyzed with or without external excitation. By comparing the time history diagrams with different values, it is found that the parameters have a certain influence on the bursting oscillation of the system.

Finally, the vibration mechanism of the system under different modes is analyzed. The vibration mode of the system is affected by both the fast- and slow-varying processes. The mechanism of different modes vibration of the system is revealed by the transformation phase portrait method, because the system trajectory will encounter different types of attractors in the fast subsystem. When the external excitation amplitude is small, the system trajectory presents the spiking state under the influence of a stable limit cycle attractor. As the external excitation amplitude increases gradually, the system trajectory will be affected by both the limit cycle attractor and stable equilibrium attractor, which shows a “bistable state” bursting situation. When the external excitation amplitude reaches a certain value, the system is only affected by the slow-varying process and presents a simple periodic vibration mode.