%ͬʱ ö ƪ ο

# Dynamic properties of a class of van der Pol-Duffing oscillators

Yelei Kuang, Xuemei Li 111 Corresponding author 000 E-mail:
yelei.k@hunnu.edu.cn(Y.Kuang), lixuemei_1@sina.com(X.Li).

Key Laboratory of High Performance Computing and Stochastic Information
Processing,

Department of Mathematics, Hunan Normal University, Changsha, Hunan 410081,
P.R. China.

Abstract. In this paper, we study the existence of bifurcation of a van der
Pol-Duffing oscillator with quintic terms and its quasi-periodic solutions by
means of qualitative and bifurcation theories. Firstly, we analyze the
autonomous system and find that it has two kinds of local bifurcations and a
global bifurcation: pitchfork bifurcation, Hopf bifurcation, homoclinic
bifurcation. It is worth noting that the disappearance of the homoclinic orbit
is synchronized with the emergence of a large limit cycle. Then, by discussing
the stability of equilibria at infinity and the orientation of the trajectory,
the existence and stability of limit circles of the autonomous system are
analyzed by combining the Poincaré-Bendixson theorem and the index theory. The
global phase portrait and the numerical simulation of the autonomous system in
different parameter values are given. Finally, the existence of periodic and
quasi-periodic solutions to periodic forced system is proved by a KAM theorem.

Key words: Homoclinic bifurcation; Hopf bifurcation; Limit circle; Stability;
Quasi-periodic solution.

##  1 Introduction and main result

In 1927, the physicist Balthasar van der Pol[1] proposed a oscillator model

$\displaystyle\ddot{x}-k(1-x^{2})\dot{x}+x=0,$ (1.1)

named after van der Pol oscillator, which is used to simulate the limit circle
oscillation phenomenon of vacuum tube amplifier, or with a periodic forcing
term

$\displaystyle\ddot{x}-k(1-x^{2})\dot{x}+x=b\lambda\sin(\lambda t),$ (1.2)

where $k$ is the damping coefficient. When $k=0$, the equation becomes an
ordinary simple harmonic vibration equation. When $k<0$, it is damped
vibration, and the vibration will gradually decay to 0; when $k>0$, self-
excited motion occurs. van der Pol’s work has become the basis for many modern
theories of nonlinear oscillations. van der Pol oscillator has no analytical
solution, but its approximate solutions are obtained by many methods, such as
perturbation method[3], multi-scale method[4], KBM method[5] and so on.

The Duffing oscillator[2] is a vibration oscillator describing forced
vibration, which is expressed by nonlinear differential equation in the
following

$\displaystyle\ddot{x}+2r\dot{x}+\alpha x+\beta x^{3}=\delta\cos(\omega t),$
(1.3)

where $\gamma$ is the damping coefficient, $\alpha$ is the toughness, $\beta$
is the nonlinearity of the controlling force, $\delta$ is the amplitude of the
driving force, and $\omega$ is the angular frequency of the driving force.

When studying some practical problems, van der Pol equation (1.1) and Duffing
equation(1.3) are usually combined together, called van der Pol-duffing
equation

$\displaystyle\ddot{x}+k(1-x^{2})\dot{x}+\alpha x+\beta
x^{3}=\delta\cos(\omega t).$ (1.4)

van der Pol-Duffing oscillator(1.4) is one of the most classical mathematical
physics models. Because of its rich dynamical properties, it is widely used in
physics, biology, engineering and even economics. There are many generalized
models of (1.4) also called van der Pol-Duffing oscillators, which have been
extensively studied. Yu, Murza et al.[10] investigated the dynamics of
autonomous ODE systems with D4-symmetry, and applied the obtained results to
two coupled van der Pol-Duffing oscillators with D4-symmetry. Xu and Luo[11]
studied the periodic driven van der Pol-Duffing oscillator and showed the
independent bifurcation tree from 2-periodic motion to the coexistence of
chaotic and 1-periodic motion by semi-analytical approach. Monsivais,
Velazquez et al. [12] investigated the dynamics of hierarchical weighted
networks of van der Pol oscillators. In addition, many scholars have
theoretically analyzed the global dynamical properties [13, 14, 15, 16, 17,
18, 19], and found attractors and chaotic phenomena [20, 21, 22] in van der
Pol-Duffing oscillator.

Based on the cell mapping method Han and Xu et al.[23] investigated the
dynamical behaviors of a model of van der Pol-Duffing oscillator with a
quintic term

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=y,\\\
\dot{y}=(\mu+x^{2}-x^{4})y+\omega_{0}^{2}x-\lambda x^{3}-\alpha x\cos(\omega
t).\end{array}\right.$ (1.7)

where $x,y$ are the state variable, $\mu,\omega_{0},\lambda$ and $\alpha$ are
real parameters. It is obvious that the system (1.7) is symmetric because the
two tracks generated by $(x,y)$ and $(-x,-y)$ are center symmetric about the
origin at all times. For the fixed parameters
$\mu=-1,\omega_{0}=1,\lambda=3,\omega=1$, though numerical simulations the
authors showed the change of dynamical behaviors of (1.7) in the region
$D=\\{(x,y):~{}|x|\leq 1.5,|y|\leq 2\\}$ as varying the amplitude parameter
$\alpha$.

The purpose of this paper is to theoretically analyze the global dynamic
properties of (1.7). For the convenience of writing, $\omega_{0}^{2}$ and
$\lambda$ are marked as $\beta$ and $\epsilon$, respectively. The system (1.7)
is written as

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=y,\\\
\dot{y}=(\mu+x^{2}-x^{4})y+\beta x-\epsilon x^{3}-\alpha x\cos(\omega
t),\end{array}\right.$ (1.10)

where $x,y$ are the state variable, $\mu,\beta,\epsilon,\alpha$ are real
parameter and $\beta\geq 0$. We will first analyze global dynamical behaviors
of the autonomous system without the external force $(\alpha=0)$, i.e

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=y,\\\
\dot{y}=(\mu+x^{2}-x^{4})y+\beta x-\epsilon x^{3}.\end{array}\right.$ (1.13)

Obviously, orbits of the system (1.13) are center symmetric about the origin.
By the transformation

$\displaystyle x\rightarrow x,~{}~{}y\rightarrow y-F(x),$

where $F(x)=-\mu x-\frac{1}{3}x^{3}+\frac{1}{5}x^{5}$. This equation (1.13)
can be changed into the $\mathbb{Z}_{2}-$equivariant Liénard equation

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=y-(-\mu
x-\frac{1}{3}x^{3}+\frac{1}{5}x^{5})\triangleq y-F(x),\\\ \dot{y}=-(-\beta
x+\epsilon x^{3})\triangleq-g(x).\end{array}\right.$ (1.16)

We will discuss the existence of limit cycles of (1.13) by the Liénard
equation (1.16).

In Section 2, we analyze the stability of equilibrium points $O(0,0)$,
$E_{1}(x_{1},0)$, $E_{2}(x_{2},0)$, where
$x_{1}=-\sqrt{\frac{\beta}{\epsilon}}$, $x_{2}=\sqrt{\frac{\beta}{\epsilon}}$,
of the system (1.13) depending on different parameter values, and find that
the pitchfork bifurcation and Hopf bifurcation occur at the equilibrium point
$O$ and the equilibrium point $E_{1}$ (or $E_{2}$), respectively. Then the
system (1.13) is transformed into a near Hamnilton system, and the branch
curve of the homoclinic orbit is obtained by the Melnikov function method. In
Section 3, we analyze the dynamical behavior of the system (1.13) at infinity
by projecting the system (1.13) onto the Poincaré disk. Then the existence and
non-existence of limit cycles of the system (1.13) in different parameter
ranges are proven based on the direction of orbits at infinity, and the global
phase diagram is also given. In Section 4, it is first proved that the system
(1.10) has a periodic solution with period $\frac{2\pi}{\omega}$ for some
parameter regions. Then by using a KAM theorem we obtain that the system
(1.10) possesses a quasi-periodic solution with two fundamental frequencies
generated by the Hopf bifurcation of the system (1.13) with periodic
perturbations (i.e., $|\alpha|$ is sufficiently small) near the points
$E_{1,2}$. Numerical simulations are also carried out in Section 5.

Our main results are as follows.

######  Theorem 1.

When $\epsilon>0$, System (1.13) goes through the following types of
bifurcations in the parameter plane $(\beta,\mu)$ with $\beta\geq 0$

(a) pitchfork bifurcation

supercritical pitchfork bifurcation curve:
$\mu_{-}=\\{(\beta,\mu):~{}\beta=0,~{}\mu<0\\}$,

subcritical pitchfork bifurcation curve:
$\mu_{+}=\\{(\beta,\mu):~{}\beta=0,~{}\mu>0\\}$;

(b) Hopf bifurcation

bifurcation curve:
$H=\\{(\beta,\mu):\mu=\frac{\beta^{2}-\epsilon\beta}{\epsilon^{2}},~{}\beta>0\\}$;

(c) Homoclinic orbit bifurcation

bifurcation curve:
$P=\\{(\beta,\mu):\mu=\frac{32\beta^{2}}{35\epsilon^{2}}-\frac{4\beta}{5\epsilon}\triangleq\mu_{3},~{}\beta>0\\}$.
The bifurcation diagram is presented in Figure 1.

See Propositions 1-2 and Lemmas 2-5 for more details.

Figure 1: when $\epsilon>0$, the bifurcation diagram of System (1.13)

######  Theorem 2.

When $\epsilon>0$ and $|\alpha|$ is sufficiently small, for most parameters in
the sufficiently small right neighborhood of the Hopf bifurcation point
$\mu=\mu_{c}(\triangleq\frac{\beta^{2}-\epsilon\beta}{\epsilon^{2}})$, the
system (1.10) has a quasi-periodic solution with two fundamental frequencies
near the point $E_{1}$ (or $E_{2}$).

See Theorem 3 for more details.

##  2 Local dynamics

In this section, we analyze the pitchfork bifurcation, Hopf bifurcation and
homoclinic orbit bifurcation of the system (1.13). Since when
$\beta=\epsilon=0$, all points on the line $y=0$ are the equilibrium points of
the system (1.13), we assume $\beta$ and $\epsilon$ are not 0 at the same time
in the sequel.

######  Lemma 1.

The system (1.13) has only one equilibrium, $O$ in the case with $\beta=0$ or
$\epsilon\leq 0$; has three equilibrium points, $O$ and $E_{1,2}$ in the case
with $\beta>0$ and $\epsilon>0$. Their properties are shown in the Table 1,
Where
$\mu_{1}=\frac{\beta^{2}-\beta\epsilon-2\epsilon^{2}\sqrt{2\beta}}{\epsilon^{2}}$,
$\mu_{2}=\frac{\beta^{2}-\beta\epsilon+2\epsilon^{2}\sqrt{2\beta}}{\epsilon^{2}}$,
$\mu_{c}=\frac{-\epsilon\beta+\beta^{2}}{\epsilon^{2}}$.

Table 1: The type and stability of equilibrium points in different parameter
value ranges

possibilities of $(\epsilon,\beta,\mu)$ | type and stability
---|---
$\epsilon\leq 0,\beta>0,\mu\in R$ |  $O(0,0)$, saddle
$\epsilon<0,\beta=0,\mu\in R$ |  $O(0,0)$, saddle
$\epsilon>0$ | $\beta=0$ | $\mu<0$ |  $O(0,0)$, stable node
$\mu=0$ |  $O(0,0)$, unstable focus
$\mu>0$ |  $O(0,0)$, unstable node
$\beta>0$ | $\mu\leq\mu_{1}$ |  $O(0,0)$, saddle; $E_{1}(-\sqrt{\frac{\beta}{\epsilon}},0)$, $E_{2}(\sqrt{\frac{\beta}{\epsilon}},0)$, stable nodes
$\mu_{1}<\mu\leq\mu_{c}$ |  $O(0,0)$, saddle; $E_{1}(-\sqrt{\frac{\beta}{\epsilon}},0)$, $E_{2}(\sqrt{\frac{\beta}{\epsilon}},0)$, stable foci
$\mu_{c}<\mu<\mu_{2}$ |  $O(0,0)$, saddle; $E_{1}(-\sqrt{\frac{\beta}{\epsilon}},0)$, $E_{2}(\sqrt{\frac{\beta}{\epsilon}},0)$, unstable foci
$\mu\geq\mu_{2}$ |  $O(0,0)$, saddle; $E_{1}(-\sqrt{\frac{\beta}{\epsilon}},0)$, $E_{2}(\sqrt{\frac{\beta}{\epsilon}},0)$, unstable nodes

###### Proof.

(I) We first prove the case where the system(1.13) has only one equilibrium.
When parameter $\beta=0$ or $\beta>0$, $\epsilon\leq 0$, solving

$\displaystyle y=0,$ $\displaystyle(\mu+x^{2}-x^{4})y+\beta x-\epsilon
x^{3}=0,$

it is easy to see system(1.13) has a unique equilibrium $O(0,0)$. Its
corresponding Jacobian matrix $A_{0}$ for system(1.13) at $O(0,0)$ is

$A_{0}=\left(\begin{array}[]{cc}0&1\\\ \beta&\mu\end{array}\right).$

For the characteristic equation is,

$\lambda^{2}-\mu\lambda-\beta=0,$

its characteristic root can be solved as

$\lambda_{1}=\frac{1}{2}(\mu+\sqrt{\mu^{2}+4\beta}),~{}~{}\lambda_{2}=\frac{1}{2}(\mu-\sqrt{\mu^{2}+4\beta}).$

Since $\beta\geq 0$, $\mu^{2}+4\beta\geq 0$ is constant. Thus, (1) when the
$\beta>0$, $\epsilon\leq 0$, $\mu\in\mathbb{R}$, $\lambda_{1},\lambda_{2}$ is
real number and $\lambda_{1}>0$, $\lambda_{2}<0$. So the equilibrium $O(0,0)$
is the saddle of the system.

(2) When the $\beta=0$, $\epsilon\neq 0$, $\mu\neq 0$,
$\lambda_{1},\lambda_{2}$ are real numbers and $\lambda_{1}=0$,
$\lambda_{2}=\mu$. So the equilibrium $O(0,0)$ is degenerate. We will analyze
its stability using [[25], Theorem 7.1 of Chapter 2], where the system is

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=y,\\\
\dot{y}=(\mu+x^{2}-x^{4})y-\epsilon x^{3}.\end{array}\right.$ (2.3)

make $x\rightarrow x+\frac{y}{\mu}$, $y\rightarrow y$,
$t\rightarrow\frac{t}{\mu}$ substitute into the system (2.3), we have

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=&\frac{\epsilon}{\mu^{2}}x^{3}+\frac{3\epsilon-\mu}{\mu^{3}}x^{2}y+\frac{3\epsilon-\mu}{\mu^{4}}xy^{2}+\frac{\epsilon-\mu}{\mu^{5}}y^{3}\\\
&+\frac{1}{\mu^{2}}x^{4}y+\frac{4}{\mu^{3}}x^{3}y^{2}+\frac{6}{\mu^{4}}x^{2}y^{3}+\frac{4}{\mu^{5}}xy^{4}+\frac{1}{\mu^{6}}y^{5}\triangleq
P(x,y),\\\
\dot{y}=&y-\frac{\epsilon}{\mu}x^{3}-\frac{3\epsilon-\mu}{\mu^{2}}x^{2}y-\frac{3\epsilon-\mu}{\mu^{3}}xy^{2}-\frac{\epsilon-\mu}{\mu^{4}}y^{3}\\\
&-\frac{1}{\mu}x^{4}y-\frac{4}{\mu^{2}}x^{3}y^{2}-\frac{6}{\mu
3}x^{2}y^{3}-\frac{4}{\mu^{4}}xy^{4}-\frac{1}{\mu^{5}}y^{5}\triangleq
y+Q(x,y).\end{array}\right.$ (2.8)

Make $y+Q(x,y)=0$ obtained $y=\varphi(x)$, since $Q(x,y)$is analytic function,
and $\varphi(0)={\varphi}\dot{(}0)=0$, so the method of undetermined
coefficients can be used for the $\varphi(x)$. Assuming
$\varphi(x)=a_{2}x^{2}+a_{3}x^{3}+a_{4}x^{4}+\cdots\cdots$, in
$\varphi(x)+Q(x,\varphi(x)=0$ and compare the coefficient get

$\displaystyle
a_{2}=0,~{}a_{3}-\frac{\epsilon}{\mu}=0,~{}a_{4}-\frac{-\mu+3\epsilon}{\mu^{2}}a_{2}=0,$
$\displaystyle
a_{5}-\frac{-\mu+3\epsilon}{\mu^{2}}a_{3}-\frac{3\epsilon-2\mu}{\mu^{3}}a^{2}_{2}=0,\cdots$

solve for

$\displaystyle
a_{2}=0,~{}a_{3}=\frac{\epsilon}{\mu},~{}a_{4}=0,a_{5}=\frac{-\mu\epsilon+3\epsilon^{2}}{\mu^{3}},\cdots,$

So we can get

$\displaystyle\varphi(x)=\frac{\epsilon}{\mu}x^{3}+\frac{-\mu\epsilon+3\epsilon^{2}}{\mu^{3}}x^{5}+\cdots.$

Substitute $y=\varphi(x)$ into $P(x,y)$ to get

$\displaystyle P(x,\varphi(x))=\frac{\epsilon}{\mu^{2}}x^{3}+\cdots,$

That is, when $m=3$, $a_{3}=\frac{\epsilon}{\mu^{2}}$. Considering the time
transformation $t=\frac{t}{\mu}$, when $\mu>0$, time transformation is
positive, when $\mu<0$, time transformation is the reverse. Hence, by [[25],
Theorem 7.1 of Chapter 2], we can get when $\mu>0$, $\epsilon>0$, the
equilibrium $O(0,0)$ is an unstable node of the system; when $\mu>0$,
$\epsilon<0$, the equilibrium $O(0,0)$ is a saddle of the system; when
$\mu<0$, $\epsilon>0$, the equilibrium $O(0,0)$ is a stable node of the
system; when $\mu<0$, $\epsilon>0$, the equilibrium $O(0,0)$ is a saddle of
the system.

(3) When $\beta=0$, $\epsilon\neq 0$, $\mu=0$, two eigenvalues are equal to
zero. In this case, the system is

$\displaystyle\left\\{\begin{array}[]{ll}\dot{x}=y,\\\
\dot{y}=(x^{2}-x^{4})y-\epsilon x^{3}.\end{array}\right.$ (2.11)

Obviously the corresponding [[25], Theorem 7.2 of Chapter 2] the
$P_{2}(x,y)\equiv 0$, and $Q_{2}(x,y(x))=-\epsilon x^{3}$,
$Q_{2y}(x,y(x))=x^{2}-x^{4}$, which in this case $m=1$, $n=2$,
$a_{2m+1}=-\epsilon$. So by [[25], Theorem 7.2 of Chapter 2], when
$\epsilon<0$, $a_{2m+1}=-\epsilon>0$, the equilibrium $O(0,0)$ is a saddle;
when $\epsilon>0$, $a_{2m+1}=-\epsilon<0$ and $n>m$, the equilibrium $O(0,0)$
is a center or focus. Let’s use the Lyapnouv function method to determine the
type of equilibrium $O(0,0)$ in the case of $\epsilon>0$.

For any $x\in(-1,1)$, construct the Lyapnouv function as

$\displaystyle V(x,y)=\frac{\epsilon}{4}x^{4}+\frac{1}{2}y^{2}.$ (2.12)

Since $\epsilon>0$, $V(x,y)\geq 0$ is a positive definite function, and

$\displaystyle\frac{dV}{dt}\mid_{(\ref{2.1.4})}=x^{2}y^{2}(1-x^{2})\geq 0.$
(2.13)

Therefore, the equilibrium $O(0,0)$ is unstable. Also because

$\displaystyle\mathrm{div}\\{y,(x^{2}-x^{4})y-\epsilon x^{3}\\}=x^{2}-x^{4},$

at $-1<x<1$ is invariant and not constant zero, so there is no closed orbit.
So in this case $O(0,0)$ is an unstable focus.

(II) When $\epsilon>0$, $\beta>0$, we know from

$\displaystyle y=0,$ $\displaystyle(\mu+x^{2}-x^{4})y+\beta x-\epsilon
x^{3}=0,$

that the system has three equilibria $O(0,0)$, $E_{1}(x_{1},0)$ and
$E_{2}(x_{2},0)$, where $x_{1}=-\sqrt{\frac{\beta}{\epsilon}}$,
$x_{2}=\sqrt{\frac{\beta}{\epsilon}}$. According to the previous analysis, the
origin $O(0,0)$ is a saddle of the system(1.13). Let’s analyze the stability
of the other two equilibria, $E_{1}(x_{1},0)$ and $E_{2}(x_{2},0)$.
$E_{1}(x_{1},0)$ and $E_{2}(x_{2},0)$ corresponding Jacobian matrix, the same
as

$B=\left(\begin{array}[]{cc}0&1\\\
-2\beta&\mu+\frac{\beta}{\epsilon}-\frac{\beta^{2}}{\epsilon^{2}}\\\
\end{array}\right).$

By Jacobian matrix to get the equilibrium $E_{1}(x_{1},0)$ and
$E_{2}(x_{2},0)$ corresponding characteristic equation(2.14)

$\lambda^{2}-(\mu+\frac{\beta}{\epsilon}-\frac{\beta^{2}}{\epsilon^{2}})\lambda+2\beta=0.$
(2.14)

The eigenvalue can be calculated as

$\lambda_{1}=\frac{\epsilon^{2}\mu+\epsilon\beta-\beta^{2}}{2\epsilon^{2}}+\frac{\sqrt{\Delta}}{2},$
$\lambda_{2}=\frac{\epsilon^{2}\mu+\epsilon\beta-\beta^{2}}{2\epsilon^{2}}-\frac{\sqrt{\Delta}}{2},$

where
$\Delta=(\mu+\frac{\beta}{\epsilon}-\frac{\beta^{2}}{\epsilon^{2}})^{2}-8\beta.$
For simplicity, let’s say
$\mu_{c}=\frac{-\epsilon\beta+\beta^{2}}{\epsilon^{2}}$,

$\displaystyle\mu_{1}$
$\displaystyle=\frac{\beta^{2}-\beta\epsilon-2\epsilon^{2}\sqrt{2\beta}}{\epsilon^{2}}=\mu_{c}-2\sqrt{2\beta},$
$\displaystyle\mu_{2}$
$\displaystyle=\frac{\beta^{2}-\beta\epsilon+2\epsilon^{2}\sqrt{2\beta}}{\epsilon^{2}}=\mu_{c}+2\sqrt{2\beta}.$

Through simple calculations we can see: when $\mu\leq\mu_{1}$, be able to get
$\Delta\geq 0$, and the characteristic values are negative real numbers, so
$E_{1}(x_{1},0)$ and $E_{2}(x_{2},0)$ are stable nodes; when $\mu\geq\mu_{2}$,
be able to get $\Delta\geq 0$, and the characteristic values are positive real
numbers, so $E_{1}(x_{1},0)$ and $E_{2}(x_{2},0)$ are unstable nodes; when
$\mu_{1}<\mu<\mu_{c}$, you can get $\Delta<0$, imaginary eigenvalues for a
negative real component, so $E_{1}(x_{1},0)$ and $E_{2}(x_{2},0)$ are stable
foci; when $\mu_{c}<\mu<\mu_{2}$, you can get $\Delta<0$, imaginary
eigenvalues for positive real part, so $E_{1}(x_{1},0)$ and $E_{2}(x_{2},0)$
are unstable foci. ∎

Based on Lemma 1, we investigate bifurcations from finite equilibria in the
following two propositions.

######  Proposition 1.

When $\beta$ passes through $0$, the autonomous system(1.13) will generate two
pitchfork bifurcations in the small neighborhood of the equilibrium $O(0,0)$:
(i) If $\epsilon>0,\mu>0$, a subcritical pitchfork bifurcation is generated at
$\beta=0$.
(ii) If $\epsilon>0,\mu<0$, a supercritical pitchfork bifurcation is generated
at $\beta=0$.

###### Proof.

By restricting the discussion of autonomous system(1.13) to the central
manifold, taking $\beta$ as the bifurcation parameter and $\beta$ as the new
independent variable, system(1.13) is rewritten as

$\left\\{\begin{array}[]{l}\dot{x}=y,\\\ \dot{y}=\mu y+\beta x-\epsilon
x^{3}+(x^{2}-x^{4})y-\epsilon x^{3},\\\ \dot{\beta}=0,\end{array}\right.$
(2.15)

where $\beta x$ is a nonlinear term. Obviously $(x,\beta,y)=(0,0,0)$ is a
fixed point with eigenvalues of $0$ and $\mu$. The corresponding eigenvectors
are

$\displaystyle v_{1}=(1,0)^{T},~{}v_{2}=(1,\mu)^{T},$

by invertible linear transformation

$X\rightarrow TX,$

where $X=(x,y)^{T}$, $T=(v_{1},v_{2}),$ system(2.15) into

$\begin{array}[]{l}\left(\begin{array}[]{c}\dot{x}\\\ \dot{y}\\\
\end{array}\right)=\left(\begin{array}[]{cc}0&0\\\ 0&\mu\\\
\end{array}\right)\left(\begin{array}[]{c}x\\\ y\\\
\end{array}\right)+\left(\begin{array}[]{c}-\frac{\beta}{\mu}x+\frac{\epsilon}{\mu}x^{3}+yF_{1}(x,y),\\\
\frac{\beta}{\mu}x-\frac{\epsilon}{\mu}x^{3}+yF_{2}(x,y)\\\
\end{array}\right),\\\ ~{}~{}~{}\dot{\beta}~{}~{}=~{}~{}0,\end{array}$ (2.16)

where $F_{1}(x,y),~{}F_{2}(x,y)$ are function of $x,~{}y$. For sufficiently
small $x$, $\beta$ has

$\displaystyle
W^{c}(0)=\\{(x,y,\beta)\in\mathbb{R}^{3}:y=h(x,\beta),h(0,0)=0,Dh(0,0)=0\\}.$

Can calculate the center manifold is $y=h(x,\beta)=0$.

So system(2.16) restrictions on the center manifold is

$\displaystyle\begin{split}\dot{x}&=-\frac{\beta}{\mu}x+\frac{\epsilon}{\mu}x^{3},\\\
\dot{\beta}&=0.\ \end{split}$ (2.17)

Make $\bar{\beta}=-\frac{\beta}{\mu}$, then the system(2.17) to write to

$\displaystyle\begin{split}\dot{x}&=\bar{\beta}x+\frac{\epsilon}{\mu}x^{3},\\\
\dot{\bar{\beta}}&=0.\ \end{split}$

Since $\beta\geq 0$, $\epsilon>0$, when $\mu>0$, there is
$\frac{\epsilon}{\mu}>0$, then the system(1.13) has a subcritical pitchfork
bifurcation; when $\mu<0$, $\frac{\epsilon}{\mu}<0$, this system (1.13) there
is supercritical pitchfork bifurcation.

∎

######  Proposition 2.

If $\epsilon>0$, $\beta>0$ and $\mu_{1}<\mu<\mu_{2}$, the equilibrium $E_{1}$
(or $E_{2}$) is a focus of the system(1.13). When $\mu>\mu_{c}$ the
system(1.13) yields a stable limit cycle in a small area of the equilibrium
$E_{1}$(or $E_{2}$).

###### Proof.

Since Norm-form is required later in the analysis of the existence of quasi-
periodic solutions for periodic forced systems, here we derive the standard
form of Hopf bifurcation. From Lemma(1) , we can see that in the parameter
plane

$H=\\{(\mu,\beta):\mu=\frac{-\epsilon\beta+\beta^{2}}{\epsilon^{2}},\epsilon>0,\beta>0\\}$

. The equilibrium of the system(1.13) $E_{1}$ (or $E_{2}$) has a pair of pure
virtual eigenvalues, $\lambda_{1,2}=\pm{\rm i}\sqrt{2\beta}$. Since the
equilibrium $E_{1}$ and $E_{2}$ are symmetric about the axis of $y$, only one
equilibrium $E_{2}$ can be discussed.

First through coordinate transformation

$\left\\{\begin{array}[]{l}x\rightarrow\sqrt{\frac{\beta}{\epsilon}}+x,\\\
y\rightarrow y,\end{array}\right.$ (2.18)

equilibrium $E_{2}$ move to the origin, get it

$\left\\{\begin{array}[]{l}\dot{x}=y;\\\ \dot{y}=-2\beta
x+(\mu+\frac{\beta}{\epsilon}-\frac{\beta^{2}}{\epsilon^{2}})y-3\epsilon\sqrt{\frac{\beta}{\epsilon}}x^{2}+2\sqrt{\frac{\beta}{\epsilon}}(1-2\frac{\beta}{\epsilon})xy\\\
\qquad-\epsilon
x^{3}+(1-6\frac{\beta}{\epsilon})x^{2}y-4\sqrt{\frac{\beta}{\epsilon}}x^{3}y-x^{4}y.\end{array}\right.$
(2.19)

Make

$\sigma(\mu)=\mu+\frac{\beta}{\epsilon}-\frac{\beta^{2}}{\epsilon^{2}},~{}\zeta(\mu)=2\beta,~{}\sigma(\mu_{c})=0,~{}\zeta(\mu_{c})=\omega^{2}_{0}=2\beta.$

Let $\lambda(\mu)=\delta(\mu)+{\rm i}\omega(\mu)$ be the eigenvalue of $E_{2}$
and $\bar{\lambda}(\mu)$ be the other eigenvalue. For small $|\mu|$, we have

$\delta(\mu)=\frac{1}{2}\sigma(\mu),~{}\omega(\mu)=\frac{1}{2}\sqrt{4\zeta(\mu)-\sigma^{2}(\mu)},$

and $\delta(\mu_{c})=0,~{}\omega(\mu_{c})=\omega_{0}=\sqrt{2\beta}>0.$
Eigenvalue $\lambda,~{}\bar{\lambda}$ corresponding eigenvector is
$v_{1},~{}v_{2}$, where

$v_{1}=(\frac{\sigma(\mu)-{\rm
i}\sqrt{4\zeta(\mu)-\sigma^{2}(\mu)}}{4\beta},1)^{T}=(\frac{\bar{\lambda}}{2\beta},1)^{T},$
$v_{2}=(\frac{\sigma(\mu)+{\rm
i}\sqrt{4\zeta(\mu)-\sigma^{2}(\mu)}}{4\beta},1)^{T}=(\frac{\lambda}{2\beta},1)^{T}.$

Take the non-degenerate matrix

$T=(v_{1},v_{2})^{T},$

and take a linear transformation

$X=TZ,$ (2.20)

where $X=(x,y)^{T}$, $Z=(z,\bar{z})^{T}$ and $z$ is conjugated with $\bar{z}$.
Then the system(2.19) can be turned into

$\left(\begin{array}[]{cc}\dot{z}\\\
\dot{\bar{z}}\end{array}\right)=\left(\begin{array}[]{cc}\lambda&0\\\
0&\bar{\lambda}\end{array}\right)\left(\begin{array}[]{cc}z\\\
\bar{z}\end{array}\right)+\left(\begin{array}[]{cc}N_{1}(Z)\\\
N_{2}(Z)\end{array}\right),$ (2.21)

where

$\displaystyle N_{1}(Z)=\sum\limits_{2\leq k+l\leq
3}g_{kl}z^{k}\bar{z}^{l}+O(|z|^{4}),$ $\displaystyle
N_{2}(Z)=\sum\limits_{2\leq k+l\leq
3}\bar{g}_{kl}z^{l}\bar{z}^{k}+O(|z|^{4}),$ $\displaystyle
g_{20}=-\frac{(\frac{\beta}{\epsilon})^{\frac{3}{2}}\lambda\bar{\lambda}(8\beta^{2}-4\beta\epsilon+3\epsilon^{2}\bar{\lambda})}{4\beta^{3}(\lambda-\bar{\lambda})},$
$\displaystyle
g_{11}=-\frac{(\frac{\beta}{\epsilon})^{\frac{3}{2}}\lambda(3\epsilon^{2}\lambda\bar{\lambda}+4\beta^{2}(\lambda+\bar{\lambda})-2\beta\epsilon(\lambda+\bar{\lambda}))}{2\beta^{3}(\lambda-\bar{\lambda})},$
$\displaystyle
g_{02}=-\frac{(\frac{\beta}{\epsilon})^{\frac{3}{2}}\lambda^{2}(8\beta^{2}-4\beta\epsilon+3\epsilon^{2}\lambda)}{4\beta^{3}(\lambda-\bar{\lambda})},$
$\displaystyle
g_{30}=-\frac{\lambda\bar{\lambda}^{2}(12\beta^{2}-2\beta\epsilon+\epsilon^{2}\bar{\lambda})}{8\beta^{3}\epsilon(\lambda-\bar{\lambda})},$
$\displaystyle
g_{21}=-\frac{\lambda\bar{\lambda}(3\epsilon^{2}\lambda\bar{\lambda}+12\beta
6(2\lambda+\bar{\lambda})-2\beta\epsilon(\bar{\lambda}+2\lambda))}{8\beta^{3}\epsilon(\lambda-\bar{\lambda})},$
$\displaystyle
g_{12}=-\frac{\lambda^{2}(3\epsilon^{2}\lambda\bar{\lambda}+12\beta
6(\lambda+2\bar{\lambda})-2\beta\epsilon(2\bar{\lambda}+\lambda))}{8\beta^{3}\epsilon(\lambda-\bar{\lambda})},$
$\displaystyle
g_{03}=-\frac{\lambda^{3}(12\beta^{2}-2\beta\epsilon+\epsilon^{2}\lambda)}{8\beta^{3}\epsilon(\lambda-\bar{\lambda})},$

$\bar{g}_{kl}$ represents the complex conjugation of $g_{kl}$.

Because system(2.21) in the first equation and the second equation is
conjugate relation, so we just need to discuss the first equation, namely

$\displaystyle\dot{z}=\lambda z+\sum\limits_{2\leq k+l\leq
3}g_{kl}z^{k}\bar{z}^{l}+O(|z|^{4}),$ (2.22)

Next, the Hopf bifurcation of the system is analyzed by the method of finding
the Nomal-Form, and the quadratic term in system(2.22) is eliminated by the
transformation

$\displaystyle z\rightarrow
z+\frac{h_{20}}{2}z^{2}+h_{11}z\bar{z}+\frac{h_{02}}{2}\bar{z}^{2},$ (2.23)

where $h_{20}=\frac{g_{20}}{\lambda}$, $h_{11}=\frac{g_{11}}{\bar{\lambda}}$,
$h_{02}=\frac{g_{02}}{2\bar{\lambda}-\lambda}$. The transformed equation is

$\dot{z}=\lambda
z+\frac{k_{30}}{6}z^{3}+\frac{k_{21}}{2}z^{2}\bar{z}+\frac{k_{12}}{2}z\bar{z}^{2}+\frac{k_{03}}{6}\bar{z}^{3}+\cdots\cdots$

where

$\displaystyle
k_{30}=g_{30}+\frac{3g^{2}_{20}}{\lambda}+\frac{3g_{11}\bar{g_{02}}}{2\lambda-\bar{\lambda}},$
$\displaystyle
k_{21}=g_{21}+\frac{g_{02}\bar{g_{02}}}{2\lambda-\bar{\lambda}}+\frac{2g_{11}\bar{g_{11}}}{\lambda}+\frac{g_{11}g_{20}(2\lambda+\bar{\lambda})}{\lambda\bar{\lambda}},$
$\displaystyle
k_{12}=g_{12}+\frac{2g_{02}\bar{g_{11}}}{\lambda}+\frac{2g_{11}^{2}}{\bar{\lambda}}+\frac{g_{11}\bar{g_{20}}}{\bar{\lambda}}-\frac{g_{02}g_{20}}{\lambda-2\bar{\lambda}},$
$\displaystyle
k_{03}=g_{03}+\frac{3g_{02}\bar{g_{20}}}{\bar{\lambda}}-\frac{3g_{11}g_{02}}{\lambda-2\bar{\lambda}},$

a further transformation is made to eliminate the cubic off-resonance term,
and the transformation is given by

$\displaystyle z\rightarrow
z+\frac{h_{30}}{6}z^{3}+\frac{h_{12}}{6}z\bar{z}^{2}+\frac{h_{03}}{2}\bar{z}^{3},$
(2.24)

where $h_{30}=\frac{g_{30}}{2\lambda}$,
$h_{12}=\frac{g_{12}}{2\bar{\lambda}}$,
$h_{03}=\frac{g_{03}}{3\bar{\lambda}-\lambda}$. Substitute it into the system,
we have

$\displaystyle\dot{z}=(\delta(\mu)+{\rm
i}\omega(\mu))z+c_{1}z^{2}\bar{z}+O(|z|^{4}),$ (2.25)

where $c_{1}$ is the function related to the parameter $\mu$,

$\displaystyle
c_{1}=\frac{g_{21}}{2}+\frac{g_{02}\bar{g}_{02}}{2(2\lambda-\bar{\lambda})}+\frac{g_{11}\bar{g}_{11}}{\lambda}+\frac{g_{11}g_{20}(2\lambda+\bar{\lambda})}{2\lambda\bar{\lambda}}.$

So it can be calculated at $\mu=\mu_{c}$

$\displaystyle c_{1}(\mu_{c})=-\frac{1}{2\beta}<0.$

Also because
$\delta(\mu)=\frac{1}{2}(\mu+\frac{\beta}{\epsilon}-\frac{\beta^{2}}{\epsilon^{2}})$,
so

$\displaystyle\frac{d\delta(\mu)}{d\mu}\mid_{\mu=\mu_{c}}=\frac{1}{2}>0.$

This shows that both the nondegeneracy and transversality conditions are
satisfied, so there is a coordinate transformation and a time transformation

$\displaystyle
z=\frac{z}{\sqrt{|l_{1}(\alpha)|}},~{}dt=(1+e_{1}(\alpha)|z|^{2})\omega(\mu)dt,$

where

$\displaystyle\alpha=\alpha(\mu)=\frac{\delta(\mu)}{\omega(\mu)},~{}e_{1}(\alpha)=\operatorname{Im}\frac{c_{1}(\mu(\alpha))}{\omega(\mu(\alpha))},~{}l_{1}(\alpha)=\operatorname{Re}\frac{c_{1}(\mu(\alpha))}{\omega(\mu(\alpha))}-\alpha~{}e_{1}(\alpha),$

and

$\displaystyle
l_{1}(0)=\frac{\operatorname{Re}c_{1}(\mu_{c})}{\omega(\mu_{c})}.$

Normalized the system(2.25)

$\displaystyle\dot{z}=(\alpha+{\rm i})z-z^{2}\bar{z}+O(|z|^{4}).$ (2.26)

Substitute $z=\rho e^{{\rm i}\varphi}$ into the equation(2.26) and separate
the real and imaginary parts to get the system in polar coordinates

$\left\\{\begin{array}[]{l}\dot{\rho}=\rho(\alpha-\rho^{2}),\\\
\dot{\varphi}=1.\end{array}\right.$ (2.27)

Can be seen when $\alpha\leq 0$, or
$\mu\leq\frac{\beta^{2}}{\epsilon^{2}}-\frac{\beta}{\epsilon}$, the
system(2.27) is only an equilibrium solution $\rho=0$. When $\alpha>0$, that
is, $\mu>\frac{\beta^{2}}{\epsilon^{2}}-\frac{\beta}{\epsilon}$, the system
has two equilibrium solutions: $\rho=0$ and $\rho=\sqrt{\alpha}$. However,
$\rho=\sqrt{\alpha}$ is a limit cycle, and since $c_{1}(\mu_{c})<0$, the first
Lyapunov coefficient $l_{1}<0$, it can be known that the limit cycle is stable
by Hopf bifurcation theory[33]. ∎

######  Lemma 2.

If $\epsilon>0$, $\beta>0$ and $\mu=\mu_{3}$, there are three equilibria in
the system(1.13), where the origin $O(0,0)$ is a saddle, and $E_{1}$ , $E_{2}$
are unstable foci. In this case, there are two homoclinic orbits that co-
reside at the saddle.

###### Proof.

In order to more conveniently study the homologous orbital bifurcation, we
introduce the following transformation:

$\displaystyle
x\rightarrow\epsilon_{1}x,~{}y\rightarrow\epsilon^{2}_{1}y,~{}\beta\rightarrow\epsilon^{2}_{1}\beta,~{}\mu\rightarrow\epsilon^{2}_{1}\mu,~{}t\rightarrow\frac{1}{\epsilon_{1}}t.$
(2.28)

The system(1.13) can be converted to

$\left\\{\begin{array}[]{l}\dot{x}=y,\\\ \dot{y}=\beta x-\epsilon
x^{3}+\epsilon_{1}(\mu y+x^{2}y-\epsilon^{2}_{1}x^{4}y).\end{array}\right.$
(2.29)

When $\epsilon_{1}=0$, the system(2.29) is

$\left\\{\begin{array}[]{l}\dot{x}=y,\\\ \dot{y}=\beta x-\epsilon
x^{3},\end{array}\right.$

this is a Hamilton system. Its energy function is

$\displaystyle
H(x,y)=\frac{1}{2}y^{2}-\frac{\beta}{2}x^{2}+\frac{\epsilon}{4}x^{4}=h.$

when $h=0$, there is a pair of homologous tracks, denoted as
$\Gamma^{0}=\\{\Gamma^{0}_{-}(t)|t\in
R\\}\bigcup\\{O\\}\bigcup\\{\Gamma^{0}_{+}(t)|t\in R\\}$, It consists of a
hyperbolic saddle $O(0,0)$ and two homoclinic orbital lines homoclinic at the
saddle $\\{\Gamma^{0}_{-}(t)|t\in R\\}$, $\\{\Gamma^{0}_{+}(t)|t\in R\\}$.
Make

$\displaystyle x(t)=$ $\displaystyle l_{1}\mathrm{sech}(t),$ $\displaystyle
y(t)=$ $\displaystyle l_{2}\mathrm{sech}(t)\tanh(t),$

by substituting the energy function equation $H(x,y)=0$, when $h=0$, the
equation for the homoclinic orbits can be written as

$\left\\{\begin{array}[]{l}x(t)=\pm\sqrt{\frac{2\beta}{\epsilon}}\mathrm{sech}(t),\\\
y(t)=\mp\frac{\sqrt{2}\beta}{\sqrt{\epsilon}}\mathrm{sech}(t)\tanh(t).\end{array}\right.$

Since the perturbation term of the system(2.29) is $\epsilon_{1}(\mu
y+x^{2}y-\epsilon^{2}_{1}x^{4}y)$, and depends on the time $t$. So the
Melnikov function can be written as:

$\displaystyle M(\beta,\epsilon,\mu)=$
$\displaystyle\int_{-\infty}^{\infty}y(t)(\mu
y(t)+x^{2}(t)y(t)-\epsilon^{2}_{1}x^{4}(t)y(t))dt$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}y(t)^{2}(\mu+x^{2}(t)-\epsilon^{2}_{1}x^{4}(t))dt$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{2\beta^{2}}{\epsilon}\mathrm{sech}^{2}(t)\tanh^{2}(t)(\mu+\frac{2\beta}{\epsilon}\mathrm{sech}^{2}(2)-\frac{4\beta^{2}\epsilon^{2}_{1}}{\epsilon^{2}}\mathrm{sech}^{4}(t))dt$
$\displaystyle=$
$\displaystyle\frac{4\beta^{2}\mu}{3\epsilon}+\frac{16\beta^{3}}{15\epsilon^{2}}-\frac{128\beta^{4}\epsilon^{2}_{1}}{105\epsilon^{3}}.$
(2.30)

In order for the homoclinic orbits to survive in the disturbance, this
requires $M(\beta,\epsilon,\mu)=0$. Therefore, from (2), the bifurcation curve
of the homoclinic orbits can be obtained:

$\displaystyle\mu=\frac{32\beta^{2}\epsilon^{2}_{1}}{35\epsilon^{2}}-\frac{4\beta}{5\epsilon}.$
(2.31)

Then substituting the inverse transformation of the parameter transformation
in (2.28) $\beta\rightarrow\frac{\beta}{\epsilon_{1}^{2}}$,
$\mu\rightarrow\frac{\mu}{\epsilon_{1}^{2}}$ into (2.31), the homoclinic orbit
bifurcation curve is obtained as follows:

$\displaystyle\mu=\frac{32\beta}{35\epsilon^{2}}-\frac{4\beta}{5\epsilon}\triangleq\mu_{3},$

That is, when $\mu=\mu_{3}$, the system(1.13) generated a homoclinic orbit. ∎

In this way, we have obtained all the bifurcations of the autonomous system
(1.13) and the corresponding bifurcation curves. In order to prove the Theorem
1, it is necessary to discuss the system(1.13) when the parameter lies in the
range. The problem of large limit rings exists when \scriptsize{1}⃝. Before
analyzing large limit cycles, it is necessary to obtain the stability of the
equilibrium at infinity of the system(1.13).

##  3 Limit cycles and homoclinic loops

In this section, we will analyze the behavior of the system (1.13) at infinity
and the existence and non-existence of large limit cycles, and give a
representative global phase diagram. First, for equilibria at infinity of the
unforced system(1.13), we conclude as follows.

######  Lemma 3.

System(1.13) in Poincaré the disc has four infinity equilibria $\bar{B}(1,0)$,
$\bar{\bar{B}}(1,0)$, $\bar{(}0,1)$C, $\bar{\bar{C}}(0,1)$, where $\bar{C}$
and $\bar{\bar{C}}$ are unstable nodes; when $\epsilon\leq 0$, $\bar{B}$ and
$\bar{\bar{B}}$ are stable nodes; when $\epsilon>0$, $\bar{B}$ and
$\bar{\bar{B}}$ are saddle.

###### Proof.

In order to analyze equilibria at infinity for the unforced system(1.13), it
is first necessary to map infinity onto the Poincaré disk. By the Poincaré
transformation

$\displaystyle x=\frac{1}{z},~{}~{}y=\frac{u}{z},$

mapping points on the equator onto the plane $(u,z)$, we get

$\displaystyle\dot{x}=-\frac{1}{z^{2}}\dot{z},~{}\dot{y}=\frac{\dot{u}z-\dot{z}u}{z^{2}},$

substituting into the system(1.13), we can get

$\left\\{\begin{array}[]{ll}\dot{u}=-u^{2}+(\mu
z^{2}+1-\frac{1}{z^{2}})\frac{u}{z^{2}}+\beta-\frac{\epsilon}{z^{2}},\\\
\dot{z}=-zu.\end{array}\right.$ (3.1)

The space of the variable $z$ at the position of the denominator, when $z=0$
is meaningless, so we need to travel through the time transformation
$dt=z^{4}d\tau$, get an equivalent system

$\left\\{\begin{array}[]{ll}\frac{du}{d\tau}=-u-\epsilon z^{2}+z^{2}u+\beta
z^{4}+\mu z^{4}u-u^{2}z^{4},\\\ \frac{dz}{d\tau}=-z^{5}u.\end{array}\right.$
(3.2)

Obviously the system(3.2) only has one equilibrium $B(0,0)$, which corresponds
to the Jacobian matrix is

$\left(\begin{array}[]{cc}-1&0\\\ 0&0\\\ \end{array}\right),$

this is a degenerate equilibrium with eigenvalues $\lambda_{1}=-1$ and
$\lambda_{2}=0.$

In order to facilitate the analysis, the variable in the equation(3.2) is
still represented by $x,y$ through the time transformation
$d\tau\rightarrow-d\tau$ to get

$\left\\{\begin{array}[]{l}\frac{dx}{d\tau}=x^{5}y\triangleq\varphi(x,y),\\\
\frac{dy}{d\tau}=y+\epsilon x^{2}-x^{2}y-\beta x^{4}-\mu
x^{4}y+x^{4}y^{2}\triangleq y+\psi(x,y).\end{array}\right.$

Implicit function $y=\phi(x)$ can be solved form $y+\psi(x,y)=0$ by the
implicit function theorem. And because $\psi(x,y)$ is an analytic function,
and $\phi(0)=\dot{\phi}(0)=0$, so the method of undetermined coefficients can
be used for the $\phi(x)$, we can assume

$\displaystyle\phi(x)=a_{2}x^{2}+a_{3}x^{3}+a_{4}x^{4}+\cdots.$

Substituting $\phi+\psi(x,\phi)=0$, compare the coefficients to determine the
implicit function

$\phi(x)=-\epsilon x^{2}+(\beta-\epsilon)x^{4}+\cdots.$

Put in $\varphi(x,y)$

$\varphi(x,\phi(x))=-\epsilon x^{7}+(\beta-\epsilon)x^{9}+\cdots,$

when$\epsilon\neq 0$, $m=7,~{}a_{7}=-\epsilon,$ when$\epsilon=0$,
$m=9,~{}a_{9}=\beta>0.$

Therefore, according to [[25], Theorem 7.1 of Chapter 2], when $\epsilon\leq
0$, $B(0,0)$ is an unstable node; when $\epsilon>0$, $B(0,0)$ is a saddle.
Notice that we did the time transform $d\tau\rightarrow-d\tau$, so we need to
return to the original time $\tau$, that is, when $\epsilon\leq 0$, $B(0,0)$
is a stable node; when $\epsilon>0$, $B(0,0)$ is a saddle.

Since the time transformation $dt=z^{4}d\tau$ has an even number of degrees,
therefore, the corresponding point of $B(0,0)$ on the Poincaré disk
$\bar{B}(1,0)$, and $\bar{\bar{B}}(-1,0)$ is consistent with the stability and
trajectory trend of $B(0,0)$.

To make $x=\frac{v}{z},~{}~{}y=\frac{1}{z}$, Plug into the system(1.13) to get

$\left\\{\begin{array}[]{ll}\dot{v}=1-\beta
v^{2}-(\mu+\frac{v^{2}}{z^{2}}-\frac{v^{4}}{z^{4}})v+\epsilon\frac{v^{4}}{z^{2}},\\\
\dot{z}=-\mu z-\beta
vz-\frac{v^{2}}{z}+\epsilon\frac{v^{3}}{z}+\frac{v^{4}}{z^{3}}.\end{array}\right.$

Similarly, you need to do the time transformation $dt=z^{4}d\tau$, which is

$\left\\{\begin{array}[]{ll}\frac{dv}{d\tau}=z^{4}+v^{5}-\mu z^{4}v-\beta
v^{2}z^{4}-v^{3}z^{2}+\epsilon v^{4}z^{2},\\\ \frac{dv}{d\tau}=-\mu
z^{5}-v^{2}z^{3}+v^{4}z-\beta vz^{5}+\epsilon v^{3}z^{3}.\end{array}\right.$
(3.3)

Obviously, $C(0,0)$ is a equilibrium of the system (3.3), and its
corresponding Jacabian matrix is the zero matrix, and both eigenvalues are 0.
Since the equilibrium is degenerate, the type of equilibrium can not be
directly judged, so it is necessary to discuss the characteristic direction
and trajectory direction. By the polar coordinate transformation
$v=r\cos\theta$, $z=r\sin\theta$, (3.3) is reduced to a system in polar
coordinates

$\left\\{\begin{array}[]{ll}\dot{r}=r^{4}\cos\theta\sin^{4}\theta+O(r^{5}),\\\
r\dot{\theta}=-r^{4}\sin^{5}\theta.\end{array}\right.$

That is

$\displaystyle\frac{1}{r}\frac{dr}{d\theta}=\frac{\cos\theta\sin^{4}\theta+O(r)}{-\sin^{5}\theta}=\frac{H(\theta)+O(1)}{G(\theta)+O(1)}.$

So the characteristic equation $G(\theta)=-\sin^{5}\theta$, the value that
satisfies $G(\theta)=0$ has $0,\pi$. While
$H(\theta)=\cos\theta\sin^{4}\theta$, and $H(0)=H(\pi)=0$. It cannot be
directly determined whether $\theta=0$ and $\theta=\pi$ are one characteristic
direction.

Below Briot-Bouquet transformation in rail line, the $v=v$, $z=z_{1}v$,
$dt_{1}=vd\tau$, system (3.3) as

$\left\\{\begin{array}[]{l}\frac{dv}{dt_{1}}=v^{4}-v^{4}z^{2}_{1}+v^{3}z^{4}_{1}-\mu
v^{4}z^{4}_{1}+\epsilon v^{5}z^{2}_{1}-\beta v^{5}z^{4}_{1},\\\
\frac{dz_{1}}{dt_{1}}=v^{2}z^{5}_{1}.\end{array}\right.$ (3.4)

At this point, the origin is still a degenerate equilibrium. By the polar
coordinate transformation $v=r\cos\theta$, $z_{1}=r\sin\theta$, get

$\displaystyle\frac{dr}{rd\theta}=\frac{\cos^{5}\theta+O(r^{1})}{-\cos^{4}\theta\sin\theta+O(r^{1})}.$

So $G(\theta)=-\cos^{4}\theta\sin\theta$, $H(\theta)=\cos^{5}\theta$, such
that the solution of $G(\theta)=0$ is
$\theta_{0}=0,\pi,\frac{\pi}{2},\frac{3\pi}{2}$.

Because of $G^{{}^{\prime}}(0)H(0)=G^{{}^{\prime}}(\pi)H(\pi)=-1<0$,
therefore, $t_{1}\rightarrow\infty$ has a unique path along $\theta=0$,
$\theta=\pi$ to enter the origin. Because
$H(\frac{\pi}{2})=H(\frac{3\pi}{2})=0$, so not sure whether the two direction
rail line into the origin. Repeat the above steps to make $v=v$,
$z_{1}=z_{2}v$, $dt_{2}=vdt_{1}$ into system(3.4)

$\left\\{\begin{array}[]{l}\frac{dv}{dt_{2}}=v^{3}-v^{5}z^{2}_{2}+\epsilon
v^{6}z^{2}_{2}+v^{6}z^{4}_{2}-\mu v^{7}z^{4}_{2}-\beta v^{8}z^{4}_{2},\\\
\frac{dz_{2}}{dt_{2}}=-v^{2}z_{2}+v^{4}z^{3}_{2}-\epsilon v^{5}z^{3}_{2}+\mu
v^{6}z^{5}_{2}+\beta v^{7}z^{5}_{2}.\end{array}\right.$ (3.5)

By polar coordinate transformation $v=r\cos\theta$, $z_{2}=r\sin\theta$, get
$G(\theta)=-2\cos^{3}\theta\sin\theta$,
$H(\theta)=\cos^{2}\theta(\cos^{2}\theta-\sin^{2}\theta)$, such that the
solution of $G(\theta)=0$ is $\theta_{0}=0,\pi,\frac{\pi}{2},\frac{3\pi}{2}$.
Due to $G^{{}^{\prime}}(0)H(0)=G^{{}^{\prime}}(\pi)H(\pi)=-2<0$, therefore,
when $t_{2}\rightarrow\infty$, there is a unique path along $\theta=0$,
$\theta=\pi$ to enter the origin. Because
$H(\frac{\pi}{2})=H(\frac{3\pi}{2})=0$, so not sure whether the two direction
rail line into the origin. Repeat the above steps to make $v=v$,
$z_{2}=z_{3}v$, $dt_{3}=vdt_{2}$ into system(3.5)

$\left\\{\begin{array}[]{l}\frac{dv}{dt_{3}}=v^{2}-v^{6}z^{2}_{3}+\epsilon
v^{7}z^{2}_{3}+v^{9}z^{4}_{3}-\mu v^{10}z^{4}_{3}-\beta v^{11}z^{4}_{3},\\\
\frac{dz_{3}}{dt_{3}}=-2vz_{3}+2v^{5}z^{3}_{3}-2\epsilon
v^{6}z^{3}_{3}-v^{8}z^{5}_{3}+2\mu v^{9}z^{5}_{3}+2\beta
v^{11}z^{5}_{3}.\end{array}\right.$ (3.6)

By polar coordinate transformation $v=r\cos\theta$, $z_{3}=r\sin\theta$, get
$G(\theta)=-3\cos^{2}\theta\sin\theta$,
$H(\theta)=\cos^{3}\theta-2\cos\theta\sin^{2}\theta$, such that the solution
of the $G(\theta)=0$ is
$\theta_{0}=0,~{}\pi,~{}\frac{\pi}{2},~{}\frac{3\pi}{2}$. Because of
$G^{{}^{\prime}}(0)H(0)=G^{{}^{\prime}}(\pi)H(\pi)=-3<0$, so when
$t_{3}\rightarrow\infty$ along $\theta=0$, $\theta=\pi$ has a unique track
into the origin. Because $H(\frac{\pi}{2})=H(\frac{3\pi}{2})=0$, so not sure
whether the two direction rail line into the origin. Although the result of
this step is the same as the previous step, you can see that the minimum
number of times has been reduced once again. To make $v=v$, $z_{3}=z_{4}v$,
$dt_{4}=vdt_{3}$ into system(3.6)

$\left\\{\begin{array}[]{l}\frac{dv}{dt_{4}}=v-v^{7}z^{2}_{4}+\epsilon
v^{8}z^{2}_{4}+v^{9}z^{4}_{4}-\mu v^{13}z^{4}_{4}-\beta v^{14}z^{4}_{4},\\\
\frac{dz_{4}}{dt_{4}}=-3z_{4}+v^{4}z^{2}_{4}+2v^{6}z^{3}_{4}-v^{8}z^{5}_{4}-v^{11}z^{5}_{4}+3\beta
v^{15}z^{5}_{4}+3\mu v^{13}z^{5}_{4}-3\epsilon
v^{7}z^{3}_{4}.\end{array}\right.$ (3.7)

In this case $O(0,0)$ is a saddle of the system(3.7). According to the
characteristics of the Briot-Bouquet transformation, it can be concluded that
$C(0,0)$ is an unstable node. And since the time transformation
$dt=z^{4}d\tau$ has an even degree, So corresponding to Poincaré disc two
equilibria $\bar{C}(0,1)$ and $\bar{\bar{C}}(0,1)$ stability and rail line to
the same as $C(0,0)$. ∎

Next, we will discuss when the autonomous system(1.13) generates limit cycle.
Below we discuss the existence of limit cycles for autonomous systems(1.13) in
different parameter ranges. For simplicity, the whole parameter space is
divided into the following five subsets:

$(c1):\left\\{\begin{array}[]{lll}\epsilon\leq 0,\\\ \beta>0,\\\
\mu\in\mathbb{R},\end{array}\right.(c2):\left\\{\begin{array}[]{lll}\epsilon<0,\\\
\beta=0,\\\
\mu\in\mathbb{R},\end{array}\right.(c3):\left\\{\begin{array}[]{lll}\epsilon>0,\\\
\beta=0,\\\ \mu\leq-\frac{5}{36},\end{array}\right.$
$(c4):\left\\{\begin{array}[]{lll}\epsilon\in\mathbb{R},\\\
\beta\in\mathbb{R},\\\
\mu\leq-\frac{1}{4},\end{array}\right.(c5):\left\\{\begin{array}[]{lll}\epsilon>0,\\\
\beta=0,\\\ \mu\geq
0,\end{array}\right.(c6):\left\\{\begin{array}[]{lll}\epsilon>0,\\\
\beta>0,\\\
\mu=\mu_{3},\end{array}\right.(c7):\left\\{\begin{array}[]{lll}\epsilon>0,\\\
\beta>0,\\\ \mu>\mu_{3}.\end{array}\right.$

######  Lemma 4.

When one of conditions (c1),(c2),(c3) and (c4) holds, System (1.13) exhibits
neither limit cycles nor homoclinic loops.

###### Proof.

When $\mu\leq-\frac{1}{4}$, we can calculate the system (1.13) divergence

$\displaystyle\mathrm{div}\\{y,(\mu+x^{2}-x^{4})y+\beta x-\epsilon x^{3}\\}=$
$\displaystyle\mu+x^{2}-x^{4}$ $\displaystyle=$
$\displaystyle-(x^{2}-\frac{1}{2})^{2}+(\mu+\frac{1}{4})$ $\displaystyle\leq$
$\displaystyle~{}0,$

Therefore, according to the Bendixon-Dulac criterion, when
$\mu\leq-\frac{1}{4}$, the system (1.13) has no closed orbit.

When $\epsilon\leq 0$, $\beta>0$ and $\epsilon<0$, $\beta=0$, the finite
equilibrium has only the origin $O(0,0)$, which is a saddle, if there is a
limit cycle must contain this saddle, and according to the index theory, The
system(1.13) has no limit cycle.

When $\epsilon>0$, $\beta=0$, $\mu\leq-\frac{5}{36}$, equilibrium $O(0,0)$ is
a stable nodes

$\displaystyle E(x,y)=\int_{0}^{x}g(s)ds+\frac{y^{2}}{2}.$

Therefore

$\displaystyle\frac{dE(x,y)}{dt}\mid_{(1.7)}=-g(x)F(x)=\epsilon
x^{4}(\mu+\frac{1}{3}x^{2}-\frac{1}{5}x^{4})\leq 0.$

If there is a limit cycle $\gamma$, then there is

$\displaystyle\oint_{\gamma}\frac{dE(x,y)}{dt}\mid_{(1.7)}=\oint_{\gamma}\epsilon
x^{4}(\mu+\frac{1}{3}x^{2}-\frac{1}{5}x^{4})dt<0.$

With the $\oint_{\gamma}\frac{dE(x,y)}{dt}=0$ contradiction, so the system
(1.13) no limit cycle.

∎

######  Lemma 5.

When the condition (c5) holds, there is a unique stable limit cycle for
system(1.13).

###### Proof.

When $\epsilon>0,\beta=0,\mu\geq 0$, by Lemma 3 we can know that
$\bar{B}(1,0)$ and $\bar{\bar{B}}(1,0)$ are saddle, $\bar{C}(0,1)$ and
$\bar{\bar{C}}(0,-1)$ are unstable nodes, so we choose the equator and these
two pairs of diametral points together to form the outer boundary of the ring
domain, and the orbitals point to the inner boundary. At this time, only one
equilibrium at the origin on the finite plane is an unstable equilibrium
(unstable node or unstable focus), and this equilibrium can be used as an
inner boundary line, thus from the Poincaré-Bendixson ring domain theorem we
can know that the system(1.13) has a stable limit cycle. ∎

Thus we prove that the system (1.13) has a large limit cycle with three
equilibria when the parameter is in the range \scriptsize{1}⃝ in the Theorem
1. Through the above analysis, we can give a representative global phase
diagram of the unforced system(1.13).

When $\epsilon\leq 0$, the system(1.13) has only one equilibrium as the
saddle, then there is no limit cycle, as shown in Figure 2(a).

When $\epsilon>0,\beta=0,\mu\geq 0$, system(1.13) has only one unstable
equilibrium (unstable node or unstable focus), in which case there is a stable
limit cycle, as shown in Figure 2(b).

When $\epsilon>0,\beta>0,\mu\leq\mu_{c}$, system(1.13) has three equilibrium,
where the origin $O(0,0)$ is the saddle and $E_{1,2}$ are stable foci, and
there is no limit cycle, as shown in Figure 2(c).

When $\epsilon>0,\beta>0,\mu_{c}<\mu<\mu_{3}$, system(1.13) has three
equilibrium, in which $O(0,0)$ is a saddle, $E_{1,2}$ are unstable foci, and
there is a stable limit cycle near $E_{1,2}$, as shown in Figure 2(d).

When $\epsilon>0,\beta>0,\mu=\mu_{3}$, system(1.13) has three equilibrium,
where the origin $O(0,0)$ is the saddle, $E_{1,2}$ are unstable foci, and
there is a homohoming orbit at the saddle point, as shown in Figure 2(e).

When $\epsilon>0,\beta>0,\mu>\mu_{3}$, system(1.13) has three equilibrium,
where the origin $O(0,0)$ is a saddle, $E_{1,2}$ are unstable foci, and there
is a stable large limit cycle containing three equilibria, as shown in Figure
2(f).

(a) $\epsilon\leq 0$

(b) $\epsilon>0,\beta=0,\mu\geq 0$

(c) $\epsilon>0,\beta>0,\mu\leq\mu_{c}$

(d) $\epsilon>0,\beta>0,\mu_{c}<\mu<\mu_{3}$

(e) $\epsilon>0,\beta>0,\mu=\mu_{3}$

(f) $\epsilon>0,\beta>0,\mu>\mu_{3}$

Figure 2: Global phase diagram for an unforced system (1.13)

##  4 The existence of periodic and quasi-periodic solutions for forced
systems

In the previous section, we studied the dynamics of autonomous systems. In
this section, we will apply a KAM theory to discuss whether the original
system has two-dimensional invariant torus. In the previous chapter, we know
that the autonomous system(1.13) will generate Hopf bifurcation when the
parameter $\mu=\mu_{c}$ is near the equilibria $E_{1,2}$, so we need to first
shift the bifurcation parameter to the critical value, and then move the
equilibrium to the origin. Then the system (1.10) is reduced to a standard
form suitable for KAM theoretical analysis by a series of transformations
similar to the Hopf bifurcation canonical form derived from autonomous
systems. Finally, the existence of a two-dimensional quasi-periodic invariant
torus of the system(1.10) is discussed by this KAM theorem. For convenience,
we have relisted(1.10) here

$\left\\{\begin{array}[]{ll}\dot{x}=y,\\\ \dot{y}=(\mu+x^{2}-x^{4})y+\beta
x-\epsilon x^{3}-\alpha x\cos(\omega t).\end{array}\right.$ (4.1)

In the previous section, we know that the system(1.13) when the parameter
$\mu=\mu_{c}$, the Hopf bifurcation will be generated near two equilibria:
$E_{1},E_{2}$. Since $E_{1}$ is symmetric with the orbit near $E_{2}$, So just
analyze the equilibrium $E_{2}(x_{2},0)$, where
$x_{2}=\sqrt{\frac{\beta}{\epsilon}}$. We get the standard form of autonomous
system by a series of transformations. Due to the addition of periodic
external forces, we guess that the system(4.1) will produce quasi-periodic
solutions of two fundamental frequencies near $E_{2}$. The standard form will
be obtained by a series of transformations similar to the previous section.
First, translate the bifurcation parameters so that $\mu=\mu_{c}+\xi$, $\xi$
are small perturbed parameters, then the system(1.13) becomes

$\left\\{\begin{array}[]{ll}\dot{x}=y,\\\
\dot{y}=(\mu_{c}+\xi+x^{2}-x^{4})y+\beta x-\epsilon x^{3}.\end{array}\right.$
(4.2)

Then we apply the translation transformation (2.18) to move $E_{2}(x_{2},0)$
to the origin,

$\left\\{\begin{array}[]{ll}\dot{x}=&y,\\\ \dot{y}=&-2\beta x+\xi y-3\epsilon
x_{1}x^{2}+\frac{2x_{1}(-2\beta+\epsilon)}{\epsilon}xy\\\ &-\epsilon
x^{3}+(1-6x^{2}_{1})x^{2}y-4x_{1}x^{3}y-x^{4}y.\end{array}\right.$ (4.3)

The corresponding Jacabian matrix of the origin is

$\left(\begin{array}[]{cc}0&1\\\ -2\beta&\xi\\\ \end{array}\right),$

Its corresponding eigenvalue is

$\displaystyle\lambda=$ $\displaystyle\frac{1}{2}(\xi+{\rm
i}\sqrt{4\beta-\xi^{2}}),$ $\displaystyle\bar{\lambda}=$
$\displaystyle\frac{1}{2}(\xi-{\rm i}\sqrt{4\beta-\xi^{2}}).$

After a series of transformations (2.20), (2.23), (2.24) similar to the
transformation in the previous chapter (4.3) can be written

$\dot{z}=\lambda(\xi)z+c_{1}(\xi)z^{2}\bar{z}+O(|z|^{4}).$ (4.4)

Next order

$\displaystyle z=\sqrt{\rho}e^{{\rm i}\varphi},$

then the system(4.4) can be reduced to

$\left\\{\begin{array}[]{ll}\dot{\rho}&=2\operatorname{Re}\lambda(\xi)\rho+2\operatorname{Re}c_{1}(\xi)\rho^{2}+P_{1}(\rho,\varphi,\xi),\\\
\dot{\varphi}&=\operatorname{Im}\lambda(\xi)+\operatorname{Im}c_{1}(\xi)\rho+P_{2}(\rho,\varphi,\xi),\end{array}\right.$
(4.5)

where

$\displaystyle\operatorname{Re}\lambda(\xi)=$ $\displaystyle\frac{1}{2}\xi,$
$\displaystyle\operatorname{Im}\lambda(\xi)=$
$\displaystyle\frac{1}{2}(\sqrt{4\beta-\xi^{2}})=\sqrt{\beta}+\frac{1}{4\sqrt{\beta}}\xi+O(\xi^{2}),$

and $P_{1}(\rho,\varphi,\xi)$ is a smooth function about $\rho$ starting from
the cubic term, and $P_{2}(\rho,\varphi,\xi)$ is a smooth function about
$\rho$ starting from the quadratic term.

######  Remark 1.

We provide the following concrete example to show that the following scaling
transformation of the parameters and the system can be achieved. Fixed
$\beta=1$, $\epsilon=1$ then $\mu_{c}=0$, can be calculated

$\displaystyle x_{1}=-1,$ $\displaystyle\operatorname{Re}\lambda(\xi)=0.5\xi,$
$\displaystyle\operatorname{Im}\lambda(\xi)=1+0.25\xi+O(\xi^{2}),$
$\displaystyle
c_{1}(\xi)\approx(2-1.65885\mathrm{i})+(2.94835-0.625\mathrm{i})\xi+O(\xi^{2}).$

For the truncation equation
$2\operatorname{Re}\lambda(\xi)\rho+2\operatorname{Re}c_{1}(\xi)\rho^{2}$ of
the polar radius of the system(4.5), it has a non-zero equilibrium solution

$\displaystyle\rho_{0}=-\frac{\operatorname{Re}\lambda(\xi)}{\operatorname{Re}c_{1}(\xi)}.$

Since $\xi$ is a small perturbation parameter, we can set
$\xi\rightarrow\varepsilon\xi$. Obviously
$\rho_{0}=\frac{1}{4}\varepsilon\xi+O(\varepsilon^{2}\xi^{2})=O(\varepsilon)$,
make a translation transformation

$\displaystyle\rho=\varepsilon^{\frac{3}{2}}I+\rho_{0},$

then the system(4.5) can be written as

$\left\\{\begin{array}[]{ll}\dot{I}=&-2\operatorname{Re}\lambda(\xi)I\varepsilon+2\operatorname{Re}c_{1}(\xi)I^{2}\varepsilon^{\frac{3}{2}}+\varepsilon^{\frac{5}{2}}\widetilde{P_{1}}(I,\varphi,\xi),\\\
\dot{\varphi}=&\operatorname{Im}\lambda(\xi)+\operatorname{Im}c_{1}(\xi)\rho_{0}+\varepsilon^{\frac{3}{2}}\operatorname{Im}c_{1}(\xi)I+\varepsilon^{\frac{3}{2}}\widetilde{P_{2}}(I,\varphi,\xi),\end{array}\right.$
(4.6)

denote as

$\left\\{\begin{array}[]{ll}\dot{I}&=\varepsilon(H_{1}(\xi)I+\varepsilon^{\frac{1}{2}}\widetilde{G_{1}}(I,\varphi,\xi)),\\\
\dot{\varphi}&=w_{1}(\xi)+\varepsilon^{\frac{3}{2}}\widetilde{G_{2}}(I,\varphi,\xi),\end{array}\right.$
(4.7)

where

$\displaystyle H_{1}(\xi)$ $\displaystyle=-\xi,$ $\displaystyle w_{1}(\xi)$
$\displaystyle=\operatorname{Im}\lambda(\varepsilon\xi)-\frac{\operatorname{Re}\lambda(\varepsilon\xi)\operatorname{Im}(\varepsilon\xi)}{\operatorname{Re}c_{1}(\varepsilon\xi)}$
$\displaystyle=w_{10}+w_{11}(\xi)\varepsilon+O(\varepsilon^{2}),$

and $\widetilde{G_{1}}(I,\varphi,\xi)$, $\widetilde{G_{2}}(I,\varphi,\xi)$ all
about $I,\varphi$ is real analytic, and about $\xi$ is sufficiently smooth on
some bounded closed set.

######  Remark 2.

In order to facilitate the proof of the following theorem, specific examples
are given: fixed $\beta=1$, $\epsilon=1$, then $\mu_{c}=0$, we can calculate
it

$\displaystyle\rho_{0}\approx-0.5\varepsilon\xi,$ $\displaystyle
H_{1}(\xi)=-\xi,$ $\displaystyle w_{1}(\xi,\varepsilon)\approx
1+0.82943\varepsilon\xi+O((\varepsilon\xi)^{2}).$

Now consider the original system(4.1), because $\alpha$is a small quantity,
can make $\alpha\rightarrow\varepsilon^{3}\alpha$, $\psi=\omega t$, after the
above similar series of transformations, The system(4.1) can be written as

$\left\\{\begin{array}[]{ll}\dot{I}&=\varepsilon[H_{1}(\xi)I+\varepsilon^{\frac{1}{2}}G_{1}(I,\varphi,\psi,\xi,\alpha,\varepsilon)],\\\
\dot{\varphi}&=w_{1}(\xi,\varepsilon)+\varepsilon^{\frac{3}{2}}G_{2}(I,\varphi,\psi,\xi,\alpha,\varepsilon),\\\
\dot{\psi}&=\omega.\end{array}\right.$ (4.8)

and $G_{1}(I,\varphi,\psi,\xi)$, $G_{2}(I,\varphi,\psi,\xi)$ are all about
$I,~{}\varphi,~{}\psi$ real analytic, and about $\xi$ is sufficiently smooth
in the region of $\xi>0$. For reduced equation(4.8) we have the following
theorem.

######  Theorem 3.

Hypothesis $\xi\in\Pi=[\frac{1}{16},1]$, So for the given $0<\gamma_{0}\ll 1$,
there are sufficiently small positive numbers $\varepsilon_{0}^{*}$, such that
for $0<\varepsilon<\varepsilon_{0}^{*}$,
$\varepsilon_{0}^{*}=o(\gamma^{4}_{0})$, there is a Cantor subset
$\Pi_{\gamma_{0}}\subset[\frac{1}{16},1]$, for any of
$\xi\in\Pi_{\gamma_{0}}$, the system(4.8) exist two quasi-periodic solutions
with the fundamental frequency $(\omega_{1}^{*},\omega)$,
$|\omega_{1}^{*}-\omega_{1}|=O(\varepsilon^{\frac{3}{2}})$ and when
$\gamma_{0}\rightarrow 0$, we have ${\rm
meas}(\Pi\setminus\Pi_{\gamma_{0}})\rightarrow 0$.

###### Proof.

In section 2 we know that the autonomous system(1.13) has a small limit cycle
in the sufficiently small right neighborhood of the bifurcation point
$\mu=\mu_{c}$ (i.e. $\xi=0$), and because we have scaled
$\xi\rightarrow\varepsilon\xi$ for the parameter $\xi$, so long as
$\varepsilon$ is small enough to value $\xi$ in a closed set, $\varepsilon\xi$
is in the sufficiently small right neighborhood of $\xi=0$. Since the KAM
theorem is used to analyze the existence of quasi-periodic solutions, it is
usually necessary to restrict $\xi$ to a closed interval leaving $\xi=0$, we
may as well take $\xi\in[\frac{1}{16},1]$ and let $\Pi=[\frac{1}{16},1]$.
System(4.8) is a special case that takes $n_{11}=n_{21}=0$, $n_{12}=1$,
$n_{22}=2$, $q_{1}=q_{7}=\frac{3}{2}$, $q_{2}=q_{5}=\frac{1}{2}$,
$q_{3}=q_{4}=q_{6}=1$ , make

$\displaystyle I_{2}=I,~{}~{}\varphi_{2}=(\varphi,\psi)^{T},$

generation into the reduced equation(4.8) can be written

$\left\\{\begin{array}[]{ll}\dot{I}_{2}&=\varepsilon[A_{2}(\xi)I_{2}+\varepsilon^{\frac{1}{2}}g_{2}(I_{2},\varphi_{2},\xi,\varepsilon)],\\\
\dot{\varphi}_{2}&=w(\xi,\varepsilon)+\varepsilon^{\frac{3}{2}}g_{4}(I_{2},\varphi_{2},\xi,\varepsilon),\\\
\end{array}\right.$ (4.9)

where

$\displaystyle A_{2}(\xi)=-\xi,$
$w^{0}(\xi,\varepsilon)=\left(\begin{array}[]{cc}w_{1}(\xi,\varepsilon)&\\\
\omega\end{array}\right).$

The assumptions conditions (H1)-(H3) and the nondegenerativity condition with
respect to frequency are verified below. Since the system(4.9) is a special
case when [[31],Theorem 2.2] takes $n_{11}=n_{21}=0$, $n_{12}=1$, $n_{22}=2$,
$q_{1}=q_{7}=\frac{3}{2}$, $q_{2}=q_{5}=\frac{1}{2}$, $q_{3}=q_{4}=q_{6}=1$,
it clearly satisfies (H1).

Since $g_{2}$, $g_{4}$ is continuously differentiable with respect to the
arguments and smooth with respect to the coordinate variables of any order,
and $n_{2}=n_{21}+n_{22}=2$, so it is desirable to $l=30$, $\alpha=1$,
$\iota=3$. Obviously $g_{i}\in
C^{l,\alpha}(\Omega\times\mathbb{T}^{n_{2}},~{}\Pi_{0})(i=2,4)$,
$l>2(\alpha+1)(\iota+2)+\alpha\iota,~{}\iota>\alpha n_{2}-1$, that is, the
hypothesis (H3) is satisfied.

For $\xi\in\Pi=[\frac{1}{16},1]$,
$\inf\limits_{\xi\in\Pi}|A_{2}(\xi)|=\frac{1}{16}$, therefore, according to
the existence theorem of implicit functions, there exist positive constants
$c_{0},~{}c_{1}$ and $\varepsilon^{*}$ such that for any
$\varepsilon\in(0,\varepsilon^{*}]$ satisfying

$\displaystyle\inf\mid\lambda(\xi)\mid=$ $\displaystyle\inf\mid
A_{2}(\xi)\mid\geq 0,$ $\displaystyle\parallel
B_{2}\parallel_{1;\Pi},~{}\parallel B^{-1}_{2}\parallel_{1;\Pi},$
$\displaystyle\parallel w^{0}\parallel_{0;\Pi}:=\sup\limits_{\xi\in\Pi}\mid
w^{0}\mid\leq c_{1},$

It is known from Remark 2

$\displaystyle\parallel\partial_{\xi}w^{0}\parallel_{\Pi}:=\sup\limits_{\xi\in\Pi}\mid\partial_{\xi}w^{0}(\xi,\varepsilon)\mid\leq
c_{1}\varepsilon.$

So hypothesis (H2) is true.

That is, assumption conditions (H1)-(H3) in [[31],Theorem 2.2], standard
type(4.9) are satisfied, so for a given $0<\gamma_{0}\ll 1$, there exists a
sufficiently small positive number $\varepsilon^{*}$ such that for
$0<\varepsilon\leq\varepsilon^{*}$, $\varepsilon^{*}=o(\gamma^{4}_{0})$, there
is a subset of Cantor $\Pi_{\gamma_{0}}\in\Pi$, For any
$\xi\in\Pi_{\gamma_{0}}$, the system(4.9) has a quasi-periodic solution and an
estimate is given. Since the frequency mapping does not satisfy the conditions
in the measure estimation of [[31],Theorem 2.3], so [[31],Theorem 2.3] cannot
be directly applied. It is shown below that the measure estimate of
$\Pi_{\gamma 0}$, $\Pi_{\gamma 0}$ is the set of $\Pi$ by removing some
parameters that make the denominator too small, because

$\displaystyle\omega^{v}$
$\displaystyle=w^{0}+O(\varepsilon^{\frac{3}{2}}),\Pi_{v}=\Pi_{v-1}\backslash\bigcup_{k}\mathfrak{R}^{v}_{k}(\gamma_{0}),$
$\displaystyle k\in$
$\displaystyle\mathbb{Z}^{2}\setminus\\{0\\},K_{v-1}<|k|_{2}\leq
K_{v},~{}v=1,2,\cdots$

such that $\Pi_{\gamma 0}=\bigcap\limits_{v=0}^{\infty}\Pi_{v}$, where
$\Pi_{0}=\Pi$,

$\displaystyle\mathfrak{R}_{k}^{v}(\gamma_{0})=\\{\xi:|(k,\omega^{v})|<\frac{\varepsilon\gamma_{0}}{|k|^{3}_{1}}\\}.$
(4.10)

$\gamma_{0}=\varepsilon^{\kappa},~{}0<\kappa\leq\frac{1}{4}$. Dig out the
parameter set $\bigcup_{k}\mathfrak{R}^{v}_{k}(r_{0})$, measure estimation of
$\mathfrak{R}^{v}_{k}(r_{0})$ is performed below, $k=(k_{1},~{}k_{2})$,
$|k|_{1}=|k_{1}|+|k_{2}|$, Substituting the expression $\omega^{v}$ into
(4.10) yields

$\displaystyle|(k,\omega^{v})|=|k_{1}w_{1}(\xi,\varepsilon)+k_{2}\omega|=|k_{1}(w_{10}+w_{11}(\xi)\varepsilon+O(\varepsilon^{2}))+k_{2}\omega|<\frac{\varepsilon\gamma_{0}}{|k|^{3}_{1}},$

where $w_{10}=w_{1}(\xi,\varepsilon)$,
$w_{11}=\partial_{\xi}w_{1}|_{\varepsilon=0}$, and $w_{10}$ is a non-zero
constant, $\inf\limits_{\xi\in\Pi}|w_{11}|>0$. make

$\phi(\xi)=|k_{1}(w_{10}+w_{11}(\xi)\varepsilon+O(\varepsilon^{2}))+k_{2}\omega|.$

The following analysis needs to be discussed in several different cases:

(i)If $k_{1}w_{10}+k_{2}\omega=0$, Since $k\neq 0$, $k_{1}\neq 0$, and
$w_{10},~{}\omega$ are constants, so for sufficiently small $\gamma_{0}$,
there is

$\displaystyle|\phi(\xi)|=|k_{1}(w_{11}(\xi)\varepsilon+O(\varepsilon^{2}))|\geq\frac{\varepsilon\gamma_{0}}{|k|^{3}_{1}},$

so $\mathfrak{R}^{v}_{k}(\gamma_{0})={\emptyset}$.

(ii)If $k_{1}w_{10}+k_{2}\omega\neq 0$;

(a)If $k_{1}=0$, then $k_{2}\neq 0$, then there is

$\displaystyle|\phi(\xi)|=|k_{2}\omega|\geq\frac{\varepsilon\gamma_{0}}{|k|^{3}_{1}},$

so $\mathfrak{R}^{v}_{k}(\gamma_{0})={\emptyset}$.

(b)If $k_{1}\neq 0$, since $w_{10}$ and $\omega$ have nothing to do with
$\xi$, the derivative of $\phi(\xi)$ with respect to $\xi$ can be obtained as

$\displaystyle\frac{d\phi(\xi)}{d\xi}=k_{1}\varepsilon\frac{dw_{11}(\xi)}{d\xi}+O(\varepsilon^{2}),$

so, by

$\displaystyle\mathrm{meas}\\{\xi\in\Pi:~{}|\phi(\xi)|<a\\}\leq\frac{2a}{\inf\limits_{\xi\in\Pi}|\phi^{{}^{\prime}}(\xi)|},$

and $k_{1}\geq 1$, Remark 3 It follows that $\frac{dw_{11}(\xi)}{d\xi}$ has a
positive lower bound, so there exists $c_{2}$ such that

$\displaystyle\mathrm{meas}~{}\mathfrak{R}^{v}_{k}(\gamma_{0})\leq
c_{2}\frac{r_{0}}{|k|^{3}_{1}}.$

Combined with the above analysis, we can get

$\displaystyle\mathrm{meas}~{}(\bigcup_{v,k}\mathfrak{R}^{v}_{k}(\gamma_{0}))\leq
c_{2}r_{0}\sum_{0\neq k\in\mathbb{Z}^{2}}\frac{1}{|k|^{3}_{1}}.$

And because $\sum\limits_{0\neq k\in\mathbb{Z}^{2}}\frac{1}{|k|^{3}_{1}}$
converges, there exists $c_{3}$, such that

$\displaystyle\mathrm{meas}~{}(\bigcup_{v,k}\mathfrak{R}^{v}_{k}(\gamma_{0}))\leq
c_{3}\gamma_{0}.$

Then there is

$\displaystyle\Pi_{\gamma
0}=\Pi\setminus(\bigcup_{v,k}\mathfrak{R}^{v}_{k}(\gamma_{0})),$

so

$\displaystyle\mathrm{meas}~{}\Pi_{\gamma_{0}}=\mathrm{meas}~{}\Pi-O(\gamma_{0}).$

Then for sufficiently small $\gamma_{0}$, the Cantor set $\Pi_{\gamma_{0}}$
defined in theorem3 has a positive Lebesgue measure, and when
$\gamma_{0}\rightarrow 0$, there is $\mathrm{meas}~{}(\Pi\setminus\Pi_{\gamma
0})\rightarrow 0$.
∎

Since the previous transformations are invertible, it can be seen from the
Theorem 3 that the system(1.10) has quasi-periodic solutions. Obviously, a
periodic solution is generated near the saddle, $O(0,0)$.

##  5 Numerical simulations

In this section we will give the phase diagram of system(1.13) by numerical
simulation and show the bifurcation. For simplicity, we use $UM$ and $SM$ to
represent unstable and stable manifolds in the simulation phase diagram. The
qualitative properties of the system(1.13) at infinity cannot be reflected in
numerical simulations. Since when $\epsilon\leq 0$ the system has only one
equilibrium that origin is a saddle and the phase diagram structure is
relatively simple, we consider the case where $\epsilon>0$ in the following.

Example 1: When $\epsilon=2$, $\beta=0$, $\mu=0$, the system has only one
equilibrium, which the origin is a unstable focus, and there is a stable limit
cycle. See the Figure 2(a).

Example 2: When $\epsilon=2$, $\beta=0$, $\mu=1$, the system has only one
unstable equilibrium that origin, and there is a stable limit cycle. See the
Figure 2(b).

Example 3: When $\epsilon=2$, $\beta=1$, $\mu=-0.25$, the system has three
equilibria, the origin is saddle, $E_{1,2}$ are stable foci. See the Figure
2(c).

Example 4: When $\epsilon=2$, $\beta=1$, $\mu=0.2$, the system has three
equilibria, the origin is saddle, $E_{1,2}$ are unstable foci, and in the
$E_{1,2}$ near have a stable limit cycle. See the Figure 2(d).

Example 5: When $\epsilon=2$, $\beta=1$, $\mu=-0.171$, the system has three
equilibria, where the origin is saddle, $E_{1,2}$ are the unstable foci, and
there is a homocyclic ring that resides at the saddle. See the Figure 2(e).

Example 6: When $\epsilon=2$, $\beta=1$, $\mu=-0.1$, the system has three
equilibria, where the origin is a saddle, $E_{1,2}$ are unstable foci, and
there is a large limit cycle containing three equilibria. See the Figure 2(f).

(a) $\epsilon=2$, $\beta=0$, $\mu=0$

(b) $\epsilon=2$, $\beta=0$, $\mu=1$

(c) $\epsilon=2$, $\beta=1$, $\mu=-0.25$

(d) $\epsilon=2$, $\beta=1$, $\mu=-0.2$

(e) $\epsilon=2$, $\beta=1$, $\mu=-0.171$

(f) $\epsilon=2$, $\beta=1$, $\mu=-0.1$

Figure 3: Numerical simulation of the unforced system (1.13)

Next, we will show the phase diagram of the system(4.1) by numerical
simulation. Take $\beta=1$, $\epsilon=3$, $\omega=1$, $\alpha=-0.3$. When
$\mu\leq-\frac{2}{9}$, there is no quasi-periodic solution for the system
(4.1), as shown in Figure 4(a) is $\mu=-0.3$. When $\mu>-\frac{2}{9}$, the
autonomous system(1.13) has a limit cycle caused by Hopf bifurcation, and when
the amplitude $\alpha$ changes in a small range, the system(4.1) has a two-
dimensional torus. The phase diagram of the system(4.1) when Figure 4(b) is
$\mu=-0.1$.

(a) Phase diagram of the system (4.1) in the $(x,y)$ plane when $\mu=?0.3$,
$\beta=1$, $\epsilon=3$, $\omega=1$, and $\alpha=?0.3$

(b) Phase diagram of the system (4.1) in the $(x,y)$ plane when $\mu=-0.1$,
$\beta=1$, $\epsilon=3$, $\omega=1$, $\alpha=-0.3$

Figure 4: Numerical simulation of forced system (4.1)

Figure 5 is the oscillatory graph of $x$ and $y$ changing with time $t$
respectively when $\mu=-0.3$. We can see that the changes of $x$ and $y$ are
irregular at this time, that is, there is no quasi-periodic solution at this
time.

(a) $x$ changes track with time $t$

(b) $y$ changes track with time $t$

Figure 5: When $\mu=-0.3$, the system (4.1) takes the initial value
$(x_{0},y_{0})=(0,1.2)$as the oscillation diagram of the change of $x$and
$y$with time $t$, respectively

Figure 6 shows the oscillation diagram of $x$ and $y$ changing with time $t$
respectively when $\mu=-0.1$. We can see that the changes of $x$ and $y$ show
a certain rule, that is, there is a quasi-periodic solution at this time.

(a) $x$ changes track with time $t$

(b) $y$ changes track with time $t$

Figure 6: When $\mu=-0.1$, the system (4.1) takes the initial value
$(x_{0},y_{0})=(0,1.2)$as the oscillation diagram of the change of $x$and
$y$with time $t$, respectively

##  6 Conclusions and Prospects

In this paper, we mainly study the dynamic properties of van der Pol-Duffing
unforced systems with quintic terms and the existence of periodic and quasi-
periodic solutions for systems with periodic external forces. Through
analysis, it is found that the unforced system(1.13) has rich dynamic
properties. The unforced system(1.13) has a total of three parameters, and
there may be 1 to 3 equilibria in the finite plane: $O(0,0)$, $E(x_{1},0)$,
$E(x_{2},0)$, where the origin exists for any value of the argument. Choosing
different parameters as bifurcation will produce different bifurcation at
different equilibrium. By numerical simulation, we find that the autonomous
system(1.13) can produce pitchfork bifurcation, Hopf bifurcation and
homoclinic orbit bifurcation, and the numerical simulation results are
strictly proved theoretically.

Pitchfork bifurcation and Hopf bifurcation are local partial bifurcation, and
the higher order term does not affect the bifurcation phase diagram. Through
the stability analysis of the equilibria, it is found that if $\beta$ is
regarded as the branching parameter, pitchfork bifurcation will be generated
near the origin $O$. If $\mu$ is regarded as a bifurcation parameter, the Hopf
bifurcation will be generated near two equilibria, $E_{1,2}$. In the analysis
of pitchfork bifurcation, the center direction and non-center direction are
obtained by reversible linear transformation of coordinate variables, and then
the central manifold of the system(1.13) is calculated to limit the system to
the central manifold and the standard form of pitchfork bifurcation is
obtained. It is found that the system(1.13) generates a subcritical pitchfork
bifurcation when $\mu>0$, and a supercritical pitchfork bifurcation when
$\mu<0$. When analyzing Hopf bifurcation, since the system(1.13) is symmetric,
we only need to analyze the Hopf bifurcation generated by an equilibrium.
Firstly, the equilibrium point $E_{2}$ is moved to the origin by coordinate
translation transformation. Then the standard form of the Hopf bifurcation is
obtained by calculating the first Lyapunov coefficient and Normal-Form method.
It is found that when $\mu>\mu_{c}$, the system (1.13) branches off a stable
small limit cycle near the equilibrium points $E_{1}$ and $(E_{2})$,
respectively.

Homoclinic orbit bifurcation is a global bifurcation. The system(1.13) is
transformed into an approximate Hamilton system by transformation, and then
the Melnikov function is calculated to obtain the homoclinic orbit when
$\mu=\mu_{3}$, and the homoclinic orbit bifurcation line is obtained. Next, we
analyze the homoclinic orbit generated by the system(1.13) by mapping method,
and prove the existence of the homoclinic orbit again from another angle, and
prove the existence of the large limit cycle (including three finite
equilibria) by using the ring domain theorem. Combined with numerical
simulation, it can be seen that there is only one large limit cycle and it is
stable. In the global analysis, we focus on when the limit cycle of the system
occurs. By calculating the stability of equilibria at infinity, we know that
when $\epsilon>0$, $\beta=0$, $\mu\geq 0$, the unforced system(1.13) has a
stable limit cycle. Then, through the Bendion-Dulac criterion and index
theory, we know that when $\epsilon\leq 0$, $\beta>0$ and $\epsilon<0$,
$\beta=0$, system(1.13) haven’t limit cycle, and there is no limit cycle for
$\mu\leq-\frac{1}{4}$. The global phase diagram of the system(1.13) is given
based on the above analysis.

Finally, the existence of quasi-periodic solutions of system(1.10) with
periodic disturbances is analyzed by KAM theorem. We reduce the system (1.10)
to normal form in polar coordinates by a shift transformation of the
parameters and a series of coordinate transformations similar to the canonical
form of the Hopf bifurcation. It is proved by KAM theorem that there are
quasi-periodic solutions with two fundamental frequencies for this canonical
form, and thus the existence of quasi-periodic solutions caused by Hopf
bifurcation near the equilibrium points $E_{1}$ and $E_{2}$ of the original
system(1.10) is proved. There is no theoretical analysis and numerical
simulation in the reference[23] for the above analysis.

In this paper, the system(1.10) has been systematically analyzed, but there is
still some work to be further studied on the unforced system (1.13), mainly in
the following two aspects:

1\. The existence of the large limit cycle of the system(1.13) when $\mu<0$
has not been discussed, we obtain from the Bendion-Dulac criterion that the
system(1.13) has no closed orbit when $\mu\leq-\frac{1}{4}$ has no closed
orbit. It is speculated that the existence of the limit cycle of the system
(1.13) when $\mu<0$ may be the same as the case of $\mu\leq-\frac{1}{4}$, but
no theoretical proof is given.

2\. When analyzing the existence and uniqueness of the large limit cycle of
the $\epsilon>0,~{}\beta>0,~{}\mu>\zeta(\mu)$system (1.13), the uniqueness of
the large limit cycle is obtained by combining numerical simulation. This is
partly because of the difficulty of proof, no rigorous theoretical proof has
been given.

Acknowledgments

This work is supported by the NNSF(11371132) of China, by Key Laboratory of
High Performance Computing and Stochastic Information Processing.

References

## References

  * [1] B. van der Pol. The nonlinear theory of electric oscillations. Proceedings of the Institute of Radio Engineers, 1934, 22(9):1051-1086.
  * [2] M. L. Cartwbight. Balthazar van der Pol. Journal of the London Mathematical Society, 1960, 35(3):367-376.
  * [3] M. Van Dyke. Perturbation Methods in Fluid Mechanics. Academic Press, 1964.
  * [4] A. H. Nayfeh and D. T. Mool. Nonlinear Oscillations. New York: John Wiley and Sons, 1979.
  * [5] Y.Yang. KBM method of analyzing limit cycle flutter of a wing with an external store and comparison with a wind-tunnel test. Journal of Sound and Vibration, 1995, 187(2):271-280.
  * [6] F. Battelli, K. J. Palmer. Chaos in the Duffing equation. Journal of Differential Equations, 1993, 101:276-301.
  * [7] I. Lopes, D. Passos, M. Nagy, et al. Oscillator models of the solar cycle: Towards the development of inversion methods. Physics, 2015, 186(1-4):535-559.
  * [8] Y. Qian, W. Zhang, B.Lin, et al. Analytical approximate periodic solutions for two-degree-of-freedom coupled van der Pol-Duffing oscillators by extended homotopy analysis method. Acta Mechanica, 2011, 219:1-14.
  * [9] P. Kumar, A. Kumar, S. Erlicher. A modified hybrid van der Pol-Duffing-Rayleigh oscillator for modelling the lateral walking force on a rigid floor. Physica D: Nonlinear Phenomena, 2017, 358:1-14.
  * [10] P. Yu, A. C. Murza, et al. Coupled oscillatory systems with D-4 symmetry and application to van der Pol oscillators. International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, 2016, 26(8):1650141.
  * [11] Y. Xu, A. C. J. Luo. Independent period-2 motions to chaos in a van der Pol-Duffing oscillator. International Journal of Bifurcation and Chaos, 2020, 2030045.
  * [12] D. Monsivais-Velazquez, K. Bhattacharya, R. A. Barrio, et al. Dynamics of hierarchical weighted networks of van der Pol oscillators[J]. Chaos, 2020, 30(12):123146.
  * [13] F. Dumortier, C. Rousseau. Cubic Lienard equations with linear damping[J]. Nonlinearity, 1990, 3:1015-1039.
  * [14] F. Dumortier, C. Herssens. Polynomial Liénard equations near infinity[J]. Journal of Differential Equations, 1999, 153:1-29.
  * [15] H.Chen, X.Chen. Dynamical analysis of a cubic Liénard system with global parameters. Nonlinearity, 2015, 28:3535-3562.
  * [16] M. R. Cndido, J. Llibre, C. Valls. Non-existence, existence, and uniqueness of limit cycles for a generalization of the Van der Pol-Duffing and the Rayleigh-Duffing oscillators. Physica D: Nonlinear Phenomena, 2020, 407:132458.
  * [17] H.Chen, M. Han, Y. Xia. Limit cycles of a Liénard system with symmetry allowing for discontinuity. Journal of Mathematical Analysis and Applications, 2018, 468:799-816.
  * [18] H.Chen, H. Yang, R. Zhang. Global dynamics of two classes of Liénard systems with $Z_{2}$-equivariance. Mathematical theory and application, 2020, 40(4):1-14.
  * [19] Z. Wang, H. Chen. A nonsmooth van der Pol-Duffing oscillator (I): The sum of indices of equilibria is 1. Discrete and Continuous Dynamical Systems-B, 2021, 27(3):1549-1589.
  * [20] J. Kengne, J. C. Chedjou, M. Kom, et al. Regular oscillations, chaos, and multistability in a system of two coupled van der Pol oscillators: numerical and experimental studies. Nonlinear Dynamics, 2014, 76(2):1119-1132.
  * [21] H. Zhao, Y.Lin, Y.Dai. Hidden attractors and dynamics of a general autonomous van der Pol-Duffing oscillator*. International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, 2014, 24(6):1450080.
  * [22] K. Rajagopal, A. J. M. Khalaf, Z. Wei, et al. Hyperchaos and coexisting attractors in a modified van der Pol-Duffing oscillator. International Journal of Bifurcation and Chaos, 2019, 1950067.
  * [23] Q. Han, W. Xu, X. Yue. Global bifurcation analysis of a Duffing-van der Pol oscillator with parametric excitation. International Journal of Bifurcation and Chaos, 2014, 24(4):333-126.
  * [24] S.Cai, X.Qiang. Introduction to Qualitative Theory of Ordinary Differential Equations. Higher Education Press. 1994, (in Chinese).
  * [25] Z.Zhang, T.Ding, W.Huang. Qualitative Theory of Differential Equations. Science Press. 1985, (in Chinese).
  * [26] F. Dumortier, J. Llibre, C. Joan. Qualitative Theory of Planar Differential Systems, Springer-Verlag Berlin Heidelberg, 2006.
  * [27] B.Liu, L.Huang. Existence of periodic solutions for nonlinear $n$th order ordinary differential equations. Acta Mathematica Sinica, 2004, 47(6):1133-1140.
  * [28] A. N. Kolmogorov. On conservation of conditionally periodic motions for a small change in Hamilton’s function. Doklady Akademii nauk SSSR, 1954, 98: 527-530.
  * [29] V. I. Arnold. Proof of a theorem of A. N. Kolmogorov on the invariance of quasi-periodic motions under small perturbations of the Hamiltonian. Russian Mathematical Surveys, 2007, 18(5):9-36.
  * [30] J. Moser. On invariant curves of area-preserving mappings of an annulus. Matematika, 1962.
  * [31] X. Li, Z. Shang. On the existence of invariant tori in non-conservative dynamical systems with degeneracy and finite differentiability. Discrete and Continuous Dynamical Systems, 2019, 39(7):4225-4257.
  * [32] H. Rüssmann. Invariant tori in non-degenerate nearly integrable Hamiltonian systems. Regular and Chaotic Dynamics, 2001, 6(2):119-204.
  * [33] Y. A. Kuznetsov. Elements of applied bifurcation theory. Applied Mathematical Sciences, 2004, 288(2): 715-730.
  * [34] C. Li, C. Rousseau. Codimension 2 symmetric homoclinic bifurcations and application to 1:2 resonance. Canadian Journal of Mathematics, 1990, 42(2): 279-287.
  * [35] W. Zhang, P. Yu. Degenerate bifurcation analysis on a parametrically and externally excited mechanical system. International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, 2001, 689-709.