Limit Cycle Bifurcations Near a Cuspidal Loop

Abstract

In this paper, we study limit cycle bifurcation near a cuspidal loop for a general near-Hamiltonian system by using expansions of the first order Melnikov functions. We give a method to compute more coefficients of the expansions to find more limit cycles near the cuspidal loop. As an application example, we considered a polynomial near-Hamiltonian system and found 12 limit cycles near the cuspidal loop and the center.

1. Introduction

As we know, in the qualitative theory of planar differential systems, one of the most important problems is to study the number of limit cycles for a near-Hamiltonian system, which is closely related to the Hilbert’s 16th problem. Many studies focused on the number of limit cycles for the following near-Hamiltonian system:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{x}& =H_y+\epsilon f(x,y,\delta ),\hfill \\ \hfill \dot{y}& =-H_x+\epsilon g(x,y,\delta ),\hfill \end{array}\end{array}$$

(1)

where H, f and g are analytic functions, $\epsilon >0$ is a small parameter, and $\delta \in D\subset {\mathbb{R}}^m$ with $m\in Z^{+}$ and D compact; see [1,2,3,4,5]. When $\epsilon =0$, system (1) becomes the following Hamiltonian system:

$$\begin{array}{c}\hfill \dot{x}=H_y,\dot{y}=-H_x.\end{array}$$

(2)

We suppose that the equation $H(x,y)=h$ defines a family periodic orbits $L_h$ of system (2), where $h\in J$ with J an open interval. The boundary of the family of periodic orbits $L_h$ may be a center, a homoclinic loop or a heteroclinic loop, among other possibilities. To study limit cycle bifurcations of system (1), the following Melnikov function

$$\begin{array}{c}\hfill M(h,\delta )={\oint }_{L_h}gdx-fdy,h\in J,\end{array}$$

(3)

plays an important role. In order to study the maximal number of limit cycles of system (1), we convert to study the maximal number of isolated zeros of M which is called the weak Hilbert’s 16th problem posed by Arnold. Some interesting advances can be found in Li et al. [6] around the weak Hilbert’s 16th problem. Recently, there are many new results have been obtained about the problem; see [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21] for example.

Suppose that the Hamiltonian system (2) has a homoclinic loop $L_0$ defined by the equation $H(x,y)=0$, passing through the origin O which is a cusp, and two families of periodic orbits

$$\begin{array}{c}\hfill L_h^{\pm }:H(x,y)=h,h\in J^{\pm },\end{array}$$

where $J^{+}=(0,h^{+})$, $J^{-}=(h^{-},0)$ for some constants $h^{+}>0$ and $h^{-}<0$. Correspondingly, there are two Melnikov functions as following

$$\begin{array}{c}\hfill M^{\pm }(h,\delta )={\oint }_{L_h^{\pm }}gdx-fdy,h\in J^{\pm }.\end{array}$$

(4)

For the case of cuspidal loop, Han et al. [15] studied the property of $M^{\pm }$ for general near-Hamiltonian system and obtained the expansions of $M^{\pm }$ near $h=0$. Later, based on Han et al. [15], Atabigi et al. [16] and Xiong [17] gave the formulas of the first eight and eleven coefficients in the expansions of $M^{\pm }$ as the cusp has order 2 and order 3, respectively.

As we know, for the purpose of estimating the maximal number of limit cycles, we need to find the maximal number of zeros of the first order Melnikov function. We can do this by using the expansion of it. In [15], the authors gave formulas of the first five coefficients in the expansions of the functions $M^{\pm }$ as the cusp point has order 1. Then using these formulas one can find up to six limit cycles. If one wants to find more limit cycles near the cuspidal loop, one needs to establish formulas of more coefficients in the expansions of the functions $M^{\pm }$ than [15]. For that purpose, in this paper we develop the idea for computing coefficients in the the expansion of M near a homoclinic loop passing through a hyperbolic saddle used in [18] to the case of a cuspidal loop. We give some conditions under which formulas of the first 11 coefficients in the expansions of $M^{\pm }$ can be obtained, as the cusp has order 1; see Theorem 1. Thus, in our case we can obtain up to 19 limit cycles near the cuspidal loop; see Theorems 2 and 3. As an application example, we prove that the following Liénard system

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \dot{x}& =y,\dot{y}\hfill & =x^2(1-x)-\epsilon \sum _{i\ge 0}^9{a}_i{x}^{i}y\end{array}\end{array}$$

can have 12 limit cycles. We can see that this paper is a continuation of [15].

We organize the paper as follows. In Section 2, we present some preliminary lemmas. In Section 3, we give the formulas for some coefficients in the expansions of $M^{\pm }$ and the conditions to obtain limit cycles near the cuspidal loop of order 1. In Section 4, as an application example, we consider a Liénard system and find 12 limit cycles near the cuspidal loop and the center.

2. Preliminary

Suppose that system (2) has a nilpotent singular point and it is at the origin. In other words, the function $H(x,y)$ satisfies $H_x(0,0)=0$, $H_y(0,0)=0$, and

$$\begin{array}{c}\hfill \frac{\partial (H_y,-H_x)}{\partial (x,y)}\ne 0,det\frac{\partial (H_y,-H_x)}{\partial (x,y)}=0.\end{array}$$

Further, we can suppose that

$$\begin{array}{c}\hfill H_{yy}(0,0)=1,H_{xy}(0,0)=H_{xx}(0,0)=0.\end{array}$$

Then, the function H at the origin has the following form

$$\begin{array}{c}\hfill H(x,y)=\frac{1}{2}{y}^2+\sum _{i+j\ge 3}{h}_{ij}{x}^i{y}^j.\end{array}$$

(5)

Applying the implicit function theorem, there is a unique $C^{\infty }$ function $\phi \left(x\right)={\displaystyle \sum _{j\ge 2}}{e}_j{x}^j$ exists such that $H_y(x,\phi \left(x\right))=0$ for $0<\left|x\right|\ll 1$. By (5), we can write

$$\begin{array}{c}\hfill H(x,\phi \left(x\right))=\sum _{j\ge 3}{h}_j{x}^j,0<\left|x\right|\ll 1.\end{array}$$

(6)

Let $k\ge 3$ be an integer such that

$$\begin{array}{c}\hfill h_k\ne 0,h_j=0,\mathrm{for}j<k.\end{array}$$

(7)

Regarding the order of the nilpotent singular point, Han et al. [15] gave the following definition.

Definition 1

([15]). Let (5)–(7) be satisfied. Then the origin is called a cusp of order m if $k=2m+1$.

We assume that the equation $H(x,y)=0$ defines a cuspidal loop $L_0$ of the unperturbed system (2), and $L_0$ passes through a cusp of order 1 at the origin, and surrounds an elementary center $C(x_c,y_c)$. We further suppose that the cuspidal loop $L_0$ is clockwise oriented. Then, two families of periodic orbits $L_h^{+}$ and $L_h^{-}$ are defined by the equation $H(x,y)=h$ for $0<h<h^{+}$ and $h_c<h<0$, respectively, where $h_c=H(x_c,y_c)$ and $h^{+}$ is a positive constant as given before. The phase portrait is shown in Figure 1.

The aim of our paper is to estimate the number of limit cycles of system (1) in the region V, where

$$\begin{array}{c}\hfill V=\mathrm{Cl}.(V_1\bigcup V_2),V_1=\bigcup _{h_c<h<0}{L}_h^{-},V_2=\bigcup _{0<h<h^{+}}{L}_h^{+},\end{array}$$

(8)

and $\mathrm{Cl}.(V_1\bigcup V_2)$ denotes the closure of $V_1\bigcup V_2$. Under the above assumptions, two Melnikov functions of system (1) have the following form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M^{+}(h,\delta )& ={\oint }_{L_h^{+}}gdx-fdy,0<h<h^{+},\hfill \\ \hfill M^{-}(h,\delta )& ={\oint }_{L_h^{-}}gdx-fdy,h_c<h<0.\hfill \end{array}\end{array}$$

(9)

By Han et al. [15] we have the following lemma.

Lemma 1

([15]). Consider the analytic system (1). Let (7) hold with $h_3<0$ and

$$\begin{array}{c}\hfill f(x,y,\delta )=\sum _{i+j\ge 1}{a}_{ij}{x}^i{y}^j,g(x,y,\delta )=\sum _{i+j\ge 1}{b}_{ij}{x}^i{y}^j.\end{array}$$

Then,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M^{-}(h,\delta )& =c_0+B_{00}{c}_1{\left|h\right|}^{\frac{5}{6}}+(c_2+O\left(\right|c_1\left|\right))h+B_{10}{c}_3{\left|h\right|}^{\frac{7}{6}}-\frac{B_{00}}{11}{c}_4{\left|h\right|}^{\frac{11}{6}}+c_5{h}^2+c_6{\left|h\right|}^{\frac{13}{6}}\hfill \\ & +c_7{\left|h\right|}^{\frac{17}{6}}+c_8{h}^3+c_9{\left|h\right|}^{\frac{19}{6}}+c_{10}{\left|h\right|}^{\frac{23}{6}}+O\left(h^4\right),0<-h\ll 1,\hfill \\ \hfill M^{+}(h,\delta )& =c_0+B_{00}^{*}{c}_1{h}^{\frac{5}{6}}+(c_2+O\left(\right|c_1\left|\right))h+B_{10}^{*}{c}_3{h}^{\frac{7}{6}}+\frac{B_{00}^{*}}{11}{c}_4{h}^{\frac{11}{6}}+c_5^{*}{h}^2+c_6^{*}{h}^{\frac{13}{6}}\hfill \\ & +c_7^{*}{h}^{\frac{17}{6}}+c_8^{*}{h}^3+c_9^{*}{h}^{\frac{19}{6}}+c_{10}^{*}{h}^{\frac{23}{6}}+O\left(h^4\right),0<h\ll 1,\hfill \end{array}\end{array}$$

(10)

where $B_{00}>0$, $B_{10}>0$, $B_{00}^{*}<0$, $B_{10}^{*}<0$ are constants and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill c_0=& {\oint }_{L_0}gdx-fdy,c_1=2\sqrt{2}{h}_{30}^{-\frac{1}{3}}(a_{10}+b_{01}),\hfill \\ \hfill c_2=& {\oint }_{L_0}(f_x+g_y-a_{10}-b_{01})dt,\hfill \\ \hfill c_3=& 2\sqrt{2}{h}_{30}^{-\frac{5}{3}}\{h_{30}[2a_{20}+b_{11}-h_{12}(a_{10}+b_{01})]+\frac{1}{3}[h_{21}^2-2h_{40}](a_{10}+b_{01})\},\hfill \\ \hfill c_4=& 9{\mu }_1^{-1}{\alpha }_{01}-2{\mu }_1^{-7}[(20{\mu }_2^3-20{\mu }_1{\mu }_2{\mu }_3+4{\mu }_1^2{\mu }_4){\alpha }_{00}+(4{\mu }_1^2{\mu }_3-10{\mu }_1{\mu }_2^2){\alpha }_{10}\hfill \\ & +4{\mu }_1^2{\mu }_2{\alpha }_{20}-{\mu }_1^3{\alpha }_{30}]\hfill \end{array}\end{array}$$

(11)

with

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mu }_1=& h_{30}^{\frac{1}{3}},{\mu }_2=-\frac{1}{6}{h}_{30}^{-\frac{2}{3}}[-2h_{40}+h_{21}^2],\hfill \\ \hfill {\mu }_3=& \frac{1}{36}{h}_{30}^{-\frac{5}{3}}\{12h_{30}[h_{50}-h_{31}{h}_{21}+h_{12}{h}_{21}^2]-4h_{40}^2+4h_{40}{h}_{21}^2-h_{21}^4\},\hfill \\ \hfill {\mu }_4=& -\frac{1}{648}{h}_{30}^{-\frac{8}{3}}\{5h_{21}^6-40h_{40}^3+30h_{40}{h}_{21}^2[2h_{40}-h_{21}^2]+144h_{30}{h}_{40}[h_{50}-h_{31}{h}_{21}+h_{12}{h}_{21}^2]\hfill \\ & +72h_{30}{h}_{21}^2[-h_{50}+h_{21}{h}_{31}-h_{21}^2{h}_{12}]+216h_{30}^2[-h_{60}-h_{22}{h}_{21}^2+h_{41}{h}_{21}+h_{03}{h}_{21}^3]\hfill \\ & +108h_{30}^2{h}_{31}^2+432h_{30}^2{h}_{12}{h}_{21}(h_{12}{h}_{21}-h_{31})\},\hfill \\ \hfill {\alpha }_{00}=& 2\sqrt{2}(a_{10}+b_{01}),{\alpha }_{10}=2\sqrt{2}[-h_{12}(a_{10}+b_{01})+2a_{20}+b_{11}],\hfill \\ \hfill {\alpha }_{20}=& 2\sqrt{2}\{(a_{10}+b_{01})[3h_{03}{h}_{21}-h_{22}+\frac{3}{2}{h}_{12}^2]-2h_{12}{a}_{20}-h_{12}{b}_{11}+3a_{30}+b_{21}-a_{11}{h}_{21}-2b_{02}{h}_{21}\},\hfill \\ \hfill {\alpha }_{30}=& 2\sqrt{2}\{(a_{10}+b_{01})[3h_{13}{h}_{21}+3h_{03}{h}_{31}+3h_{12}{h}_{22}-15h_{12}{h}_{03}{h}_{21}-\frac{5}{2}{h}_{12}^3-h_{32}]\hfill \\ & +(3a_{11}+6b_{02})h_{12}{h}_{21}-2(b_{12}+a_{21})h_{21}-(a_{11}+2b_{02})h_{31}+4a_{40}+b_{31}\hfill \\ & +(3b_{11}+6a_{20})h_{03}{h}_{21}-(2a_{20}+b_{11})h_{22}+(3a_{20}+\frac{3}{2}{b}_{11})h_{12}^2-(3a_{30}+b_{21})h_{12}\},\hfill \\ \hfill {\alpha }_{01}=& 2\sqrt{2}\{\frac{2}{3}{a}_{12}+2b_{03}-2h_{03}{a}_{11}-4h_{03}{b}_{02}+(a_{10}+b_{01})[5h_{03}^2-2h_{04}]\}.\hfill \end{array}\end{array}$$

(12)

We suppose that

$$\begin{array}{c}\hfill H(x,y)=h_c+\frac{1}{2}[{(x-x_c)}^2+{(y-y_c)}^2]+O\left(\right|x-x_c,y-y_c{|}^3),\end{array}$$

for $(x,y)$ near $(x_c,y_c)$. About the expansion of $M^{-}$ near the center $C(x_c,y_c)$, we have (see [19])

$$\begin{array}{c}\hfill M^{-}(h,\delta )=\sum _{j\ge 0}{b}_j\left(\delta \right){(h-h_c)}^{j+1},0<h-h_c\ll 1.\end{array}$$

(13)

The formulas of $b_j$, $j=0,1,2,3$, can be found in [20] and more coefficients can be obtained by using the programs in [21]. The formula of $b_0\left(\delta \right)$ is

$$\begin{array}{c}\hfill b_0\left(\delta \right)=T_0{\overline{b}}_0\left(\delta \right),\end{array}$$

(14)

where ${\overline{b}}_0\left(\delta \right)=(f_x+g_y)(C,\delta )$ and $T_0>0$ is a constant.

3. Main Results

In this section, we use Lemma 1 and the method given by Tian and Han in [18] to obtain more coefficients in the expansions of $M^{\pm }$. Then, we will get more limit cycles in the neighborhood of the cuspidal loop $L_0$ and the center. The following lemma was obtained in Han [22].

Lemma 2

([22]). For the Melnikov functions $M^{\pm }$ defined by (4) we have

$$\begin{array}{c}\hfill \frac{\partial M^{\pm }}{\partial h}={\oint }_{L_h^{\pm }}(f_x+g_y)dt.\end{array}$$

(15)

We assume that there are analytic functions $P_1(x,y,\delta )$ and $Q_1(x,y,\delta )$ defined on V such that

$$\begin{array}{c}\hfill f_x+g_y=H_x(x,y)P_1(x,y,\delta )+H_y(x,y)Q_1(x,y,\delta ),(x,y)\in V\end{array}$$

(16)

for ${\overline{b}}_0=c_1=c_3=0$, where V is defined in (8). Let

$$\begin{array}{c}\hfill P_1(x,y,\delta )=\sum _{i+j\ge 1}{\tilde{a}}_{ij}{x}^i{y}^j,Q_1(x,y,\delta )=\sum _{i+j\ge 1}{\tilde{b}}_{ij}{x}^i{y}^j\end{array}$$

(17)

for $(x,y)$ near the origin. Differentiating (10) and (13) with respect to h, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial M^{-}}{\partial h}& =2b_1(h-h_c)+3b_2{(h-h_c)}^2+4b_3{(h-h_c)}^3+\cdots ,0<h-h_c\ll 1,\hfill \\ \hfill \frac{\partial M^{-}}{\partial h}& =c_2+\frac{B_{00}}{6}{c}_4{\left|h\right|}^{\frac{5}{6}}+2c_{5}h-\frac{13}{6}{c}_6{\left|h\right|}^{\frac{7}{6}}-\frac{17}{6}{c}_7{\left|h\right|}^{\frac{11}{6}}+\cdots ,0<-h\ll 1,\hfill \\ \hfill \frac{\partial M^{+}}{\partial h}& =c_2+\frac{B_{00}^{*}}{6}{c}_4{h}^{\frac{5}{6}}+2c_5^{*}h+\frac{13}{6}{c}_6^{*}{h}^{\frac{7}{6}}+\frac{17}{6}{c}_7^{*}{h}^{\frac{11}{6}}+\cdots ,0<h\ll 1.\hfill \end{array}\end{array}$$

(18)

Next, applying Lemma 2 and (16) we have

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \frac{\partial M^{\pm }}{\partial h}& ={\oint }_{L_h^{\pm }}(H_x{P}_1+H_y{Q}_1)dt\hfill & ={\oint }_{L_h^{\pm }}{Q}_{1}dx-P_{1}dy\equiv M_1^{\pm }(h,\delta ).\end{array}\end{array}$$

(19)

By (10), (11), (13) and (14), it can be seen that $M_1^{-}(h,\delta )$ and $M_1^{+}(h,\delta )$ have the form

$$\begin{array}{c}\hfill M_1^{-}(h,\delta )=b_{1,0}\left(\delta \right)(h-h_c)+b_{1,1}\left(\delta \right){(h-h_c)}^2+\cdots ,\end{array}$$

(20)

for $0<h-h_c\ll 1$,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M_1^{-}(h,\delta )& =c_{1,0}+B_{00}{c}_{1,1}{\left|h\right|}^{\frac{5}{6}}+(c_{1,2}+O\left(\right|c_{1,1}\left|\right))h+B_{10}{c}_{1,3}{\left|h\right|}^{\frac{7}{6}}-\frac{B_{00}}{11}{c}_{1,4}{\left|h\right|}^{\frac{11}{6}}+\cdots ,\hfill \end{array}\end{array}$$

(21)

for $0<-h\ll 1$, and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M_1^{+}(h,\delta )& =c_{1,0}+B_{00}^{*}{c}_{1,1}{h}^{\frac{5}{6}}+(c_{1,2}+O\left(\right|c_{1,1}\left|\right))h+B_{10}^{*}{c}_{1,3}{h}^{\frac{7}{6}}+\frac{B_{00}^{*}}{11}{c}_{1,4}{h}^{\frac{11}{6}}+\cdots ,\hfill \end{array}\end{array}$$

(22)

for $0<h\ll 1$, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b_{1,0}=& T_0{\overline{b}}_1,{\overline{b}}_1=(P_{1x}+Q_{1y})(C,\delta ),\hfill \\ \hfill c_{1,1}=& 2\sqrt{2}{h}_{30}^{-\frac{1}{3}}({\tilde{a}}_{10}+{\tilde{b}}_{01}),c_{1,2}={\oint }_{L_0}(P_{1x}+Q_{1y}-{\tilde{a}}_{10}-{\tilde{b}}_{01})dt,\hfill \\ \hfill c_{1,3}=& 2\sqrt{2}{h}_{30}^{-\frac{5}{3}}\{h_{30}[2{\tilde{a}}_{20}+{\tilde{b}}_{11}-h_{12}({\tilde{a}}_{10}+{\tilde{b}}_{01})]+\frac{1}{3}[h_{21}^2-2h_{40}]({\tilde{a}}_{10}+{\tilde{b}}_{01})\},\hfill \\ \hfill c_{1,4}=& 9{\mu }_1^{-1}{\alpha }_{01}-2{\mu }_1^{-7}[(20{\mu }_2^3-20{\mu }_1{\mu }_2{\mu }_3+4{\mu }_1^2{\mu }_4){\alpha }_{00}+(4{\mu }_1^2{\mu }_3-10{\mu }_1{\mu }_2^2){\alpha }_{10}\hfill \\ & +4{\mu }_1^2{\mu }_2{\alpha }_{20}-{\mu }_1^3{\alpha }_{30}],\hfill \end{array}\end{array}$$

(23)

and $T_0>0$ is the same as before; ${\mu }_i$ and ${\alpha }_{01},{\alpha }_{i0},i=1,2,3,4$ satisfy (12), with $a_{ij}$ and $b_{ij}$ replaced by ${\tilde{a}}_{ij}$ and ${\tilde{b}}_{ij}$, respectively.

Let

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\overline{c}}_4=c_{1,1},{\overline{c}}_5=c_{1,2},{\overline{c}}_6=c_{1,3},{\overline{c}}_7=c_{1,4}.\end{array}\end{array}$$

(24)

Now by comparing the three expansions of $M_1^{\pm }$ above with (18), we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b_1& =\frac{T_0}{2}{\overline{b}}_1,c_4=6{\overline{c}}_4,\hfill \\ \hfill c_5& =c_5^{*}=\frac{1}{2}({\overline{c}}_5+O\left(\right|{\overline{c}}_4\left|\right)),\hfill \\ \hfill c_6& =-\frac{6}{13}{B}_{10}{\overline{c}}_6,c_7=\frac{6}{187}{B}_{00}{\overline{c}}_7,\hfill \\ \hfill c_6^{*}& =\frac{6}{13}{B}_{10}^{*}{\overline{c}}_6,c_7^{*}=\frac{6}{187}{B}_{00}^{*}{\overline{c}}_7.\hfill \end{array}\end{array}$$

(25)

Applying Lemma 2 to (19) again, we have

$$\begin{array}{c}\hfill \frac{{\partial }^2{M}^{\pm }}{\partial h^2}=\frac{\partial M_1^{\pm }}{\partial h}={\oint }_{L_h^{\pm }}(P_{1x}+Q_{1y})dt.\end{array}$$

(26)

Suppose there are analytic functions $P_2(x,y,\delta )$ and $Q_2(x,y,\delta )$ defined on V such that for ${\overline{b}}_0=c_1=c_3={\overline{b}}_1={\overline{c}}_4={\overline{c}}_6=0,$

$$\begin{array}{c}\hfill P_{1x}+Q_{1y}=H_x(x,y)P_2(x,y,\delta )+H_y(x,y)Q_2(x,y,\delta ),(x,y)\in V,\end{array}$$

(27)

where

$$\begin{array}{c}\hfill \begin{array}{c}\hfill P_2(x,y,\delta )=\sum _{i+j\ge 1}{\overline{a}}_{ij}{x}^i{y}^j,Q_2(x,y,\delta )=\sum _{i+j\ge 1}{\overline{b}}_{ij}{x}^i{y}^j\end{array}\end{array}$$

(28)

for $(x,y)$ near the origin. By Lemma 2, we further have

$$\begin{array}{c}\hfill \frac{{\partial }^2{M}^{\pm }}{\partial h^2}={\oint }_{L_h^{\pm }}{Q}_{2}dx-P_{2}dy\equiv M_2^{\pm }(h,\delta ).\end{array}$$

Then, by applying formulas (10), (11), (13) and (14) to the functions $M_2^{\pm }(h,\delta )$ again, we obtain

$$\begin{array}{c}\hfill M_2^{-}(h,\delta )=b_{2,0}\left(\delta \right)(h-h_c)+b_{2,1}\left(\delta \right){(h-h_c)}^2+\cdots ,\end{array}$$

(29)

for $0<h-h_c\ll 1$,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M_2^{-}(h,\delta )& =c_{2,0}+B_{00}{c}_{2,1}{\left|h\right|}^{\frac{5}{6}}+(c_{2,2}+O\left(\right|c_{2,1}\left|\right))h+B_{10}{c}_{2,3}{\left|h\right|}^{\frac{7}{6}}-\frac{B_{00}}{11}{c}_{2,4}{\left|h\right|}^{\frac{11}{6}}+\cdots ,\hfill \end{array}\end{array}$$

(30)

for $0<-h\ll 1$, and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M_2^{+}(h,\delta )& =c_{2,0}+B_{00}^{*}{c}_{2,1}{h}^{\frac{5}{6}}+(c_{2,2}+O\left(\right|c_{2,1}\left|\right))h+B_{10}^{*}{c}_{2,3}{h}^{\frac{7}{6}}+\frac{B_{00}^{*}}{11}{c}_{2,4}{h}^{\frac{11}{6}}+\cdots ,\hfill \end{array}\end{array}$$

(31)

for $0<h\ll 1$, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b_{2,0}=& T_0{\overline{b}}_2,{\overline{b}}_2=(P_{2x}+Q_{2y})(C,\delta ),\hfill \\ \hfill c_{2,1}=& 2\sqrt{2}{h}_{30}^{-\frac{1}{3}}({\overline{a}}_{10}+{\overline{b}}_{01}),c_{2,2}={\oint }_{L_0}(P_{1x}+Q_{1y}-{\overline{a}}_{10}-{\overline{b}}_{01})dt,\hfill \\ \hfill c_{2,3}=& 2\sqrt{2}{h}_{30}^{-\frac{5}{3}}\{h_{30}[2{\overline{a}}_{20}+{\overline{b}}_{11}-h_{12}({\overline{a}}_{10}+{\overline{b}}_{01})]+\frac{1}{3}[h_{21}^2-2h_{40}]({\overline{a}}_{10}+{\overline{b}}_{01})\},\hfill \\ \hfill c_{2,4}=& 9{\mu }_1^{-1}{\alpha }_{01}-2{\mu }_1^{-7}[(20{\mu }_2^3-20{\mu }_1{\mu }_2{\mu }_3+4{\mu }_1^2{\mu }_4){\alpha }_{00}+(4{\mu }_1^2{\mu }_3-10{\mu }_1{\mu }_2^2){\alpha }_{10}\hfill \\ & +4{\mu }_1^2{\mu }_2{\alpha }_{20}-{\mu }_1^3{\alpha }_{30}],\hfill \end{array}\end{array}$$

(32)

and $T_0>0$ is the constant as given before; ${\mu }_i$ and ${\alpha }_{01},{\alpha }_{i0},i=1,2,3,4$ satisfy (12), with $a_{ij}$ and $b_{ij}$ replaced by ${\overline{a}}_{ij}$ and ${\overline{b}}_{ij}$, respectively. Let

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\overline{c}}_8=c_{2,2},{\overline{c}}_9=c_{2,3},{\overline{c}}_{10}=c_{2,4}.\end{array}\end{array}$$

(33)

Under the conditions of ${\overline{b}}_0=c_1=c_3={\overline{b}}_1={\overline{c}}_4={\overline{c}}_6=0$, we take the derivative of (18) with respect to h again, and then have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{{\partial }^2{M}^{-}}{\partial h^2}& =6b_2(h-h_c)+12b_3{(h-h_c)}^2+\cdots ,0<h-h_c\ll 1,\hfill \\ \hfill \frac{{\partial }^2{M}^{-}}{\partial h^2}& =2c_5+\frac{187}{36}{c}_7{\left|h\right|}^{\frac{5}{6}}+6c_{8}h+\frac{247}{36}{c}_9{\left|h\right|}^{\frac{7}{6}}+\frac{391}{36}{c}_{10}{\left|h\right|}^{\frac{11}{6}}+\cdots ,0<-h\ll 1,\hfill \\ \hfill \frac{{\partial }^2{M}^{+}}{\partial h^2}& =2c_5^{*}+\frac{187}{36}{c}_7^{*}{h}^{\frac{5}{6}}+6c_8^{*}h+\frac{247}{36}{c}_9^{*}{h}^{\frac{7}{6}}+\frac{391}{36}{c}_{10}^{*}{h}^{\frac{11}{6}}+\cdots ,0<h\ll 1.\hfill \end{array}\end{array}$$

(34)

By (30) and (34) we can obtain $c_7=\frac{36}{187}{B}_{00}{c}_{2,1}$. Thus, by the formula of $c_7$ in (25) it can be seen that ${\overline{c}}_7$ can be rewritten as ${\overline{c}}_7=6c_{21}$. By comparing (34) with (29)–(31) respectively, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b_2& =\frac{T_0}{6}{\overline{b}}_2,c_8=c_8^{*}=\frac{1}{6}({\overline{c}}_8+O\left(\right|{\overline{c}}_7\left|\right)),\hfill \\ \hfill c_9& =\frac{36}{247}{B}_{10}{\overline{c}}_9,c_{10}=-\frac{36}{4301}{B}_{00}{\overline{c}}_{10},\hfill \\ \hfill c_9^{*}& =\frac{36}{247}{B}_{10}^{*}{\overline{c}}_9,c_{10}^{*}=\frac{36}{4301}{B}_{00}^{*}{\overline{c}}_{10}.\hfill \end{array}\end{array}$$

(35)

Summarizing the results above, we obtain the following theorem.

Theorem 1.

(i) Suppose there exist analytic functions $P_1$, $Q_1$ such that (16) holds for ${\overline{b}}_0=c_1=c_3$. Then (25) holds. (ii) Under the conditions of (i), further suppose there are analytic functions $P_2$ and $Q_2$ such that (27) holds for ${\overline{b}}_0=c_1=c_3={\overline{b}}_1={\overline{c}}_4={\overline{c}}_6=0$. Then (35) holds.

Now we apply Theorem 1 to study the number of limit cycles bifurcating both from the cuspidal loop and the center of system (1).

Theorem 2.

Consider the system (1). Suppose that the conditions of Theorem 1 (i) are satisfied, and let ${\overline{c}}_i\left(\delta \right)=c_i\left(\delta \right)$, $0\le i\le 3$. If there exist ${\delta }_0\in {\mathbb{R}}^m$ and $5\le l\le 7$ such that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{b}}_0\left({\delta }_0\right)& =0,{\overline{c}}_j\left({\delta }_0\right)=0,j=0,\cdots ,l-1,\hfill \\ \hfill \mu & ={(-1)}^{(l+1)}{\overline{b}}_1\left({\delta }_0\right){\overline{c}}_l\left({\delta }_0\right)<0,(resp.,>0),\hfill \\ \hfill \mathrm{rank}& \frac{\partial ({\overline{b}}_0,{\overline{c}}_0,\cdots ,{\overline{c}}_{l-1})}{\partial ({\delta }_1,{\delta }_2,\cdots ,{\delta }_m)}\left({\delta }_0\right)=l+1,\hfill \end{array}\end{array}$$

(36)

then system (1) can have $\left[\frac{3l+6}{2}\right](resp.,\left[\frac{3l+4}{2}\right])$ limit cycles for some $(\epsilon ,\delta )$ near $(0,{\delta }_0)$ as ε sufficiently small.

Proof. 

We proof the case of $l=7$. Let $\mu <0$. For definiteness, assume that ${\overline{b}}_1\left({\delta }_0\right)<0$ and ${\overline{c}}_7\left({\delta }_0\right)>0$. Then by (10), (13) and (25) we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M^{-}(h,{\delta }_0)& =c_7\left({\delta }_0\right){\left|h\right|}^{\frac{17}{6}}+{O\left(\right|h|}^3)>0,0<-h\ll 1,\hfill \\ \hfill M^{-}(h,{\delta }_0)& =b_1\left({\delta }_0\right){(h-h_c)}^2+O\left({(h-h_c)}^3\right)<0,0<h-h_c\ll 1.\hfill \end{array}\end{array}$$

Thus, there exists $h_0\in (h_c,0)$ such that

$$\begin{array}{c}\hfill M^{-}(h_0,{\delta }_0)=0,M^{-}(h_0-{\epsilon }_0,{\delta }_0)M^{-}(h_0+{\epsilon }_0,{\delta }_0)<0\end{array}$$

for ${\epsilon }_0>0$ sufficiently small. From (10), (13) and (25) we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M^{-}(h,\delta )=& b_0(h-h_c)+b_1{(h-h_c)}^2+O\left(h^3\right),0<h-h_c\ll 1,\hfill \\ \hfill M^{-}(h,\delta )=& {\overline{c}}_0+B_{00}{\overline{c}}_1{\left|h\right|}^{\frac{5}{6}}+({\overline{c}}_2+O\left(\right|{\overline{c}}_1\left|\right))h+B_{10}{\overline{c}}_3{\left|h\right|}^{\frac{7}{6}}-\frac{6}{11}{B}_{00}{\overline{c}}_4{\left|h\right|}^{\frac{11}{6}}+\frac{1}{2}({\overline{c}}_5+O\left(\right|{\overline{c}}_4\left|\right))h^2\hfill \\ & -\frac{6}{13}{B}_{10}{\overline{c}}_6{\left|h\right|}^{\frac{13}{6}}+\frac{6}{187}{B}_{00}{\overline{c}}_7{\left|h\right|}^{\frac{17}{6}}+{O\left(\right|h|}^3),0<-h\ll 1,\hfill \\ \hfill M^{+}(h,\delta )=& {\overline{c}}_0+B_{00}^{*}{\overline{c}}_1{h}^{\frac{5}{6}}+({\overline{c}}_2+O\left(\right|{\overline{c}}_1\left|\right))h+B_{10}^{*}{\overline{c}}_3{h}^{\frac{7}{6}}+\frac{6}{11}{B}_{00}^{*}{\overline{c}}_4{h}^{\frac{11}{6}}+\frac{1}{2}({\overline{c}}_5+O\left(\right|{\overline{c}}_4\left|\right))h^2\hfill \\ & +\frac{6}{13}{B}_{10}^{*}{\overline{c}}_6{h}^{\frac{13}{6}}+\frac{6}{187}{B}_{00}^{*}{\overline{c}}_7{h}^{\frac{17}{6}}+{O\left(\right|h|}^3),0<h\ll 1.\hfill \end{array}\end{array}$$

By (36), the coefficients ${\overline{b}}_0$, ${\overline{c}}_j,j=0,\cdots ,6$, can be taken as free parameters. First, we take ${\overline{c}}_i=0(0\le i\le 5)$, $0<{\overline{c}}_6\ll 1$. In this step, there exists a simple zero $h_1$ of the function $M^{-}$ for $0<-h\ll 1$. Next, if we take ${\overline{c}}_i=0(0\le i\le 4)$ and $0<{\overline{c}}_5\ll {\overline{c}}_6\ll 1$, then there exists a simple zero $h_2$ of $M^{-}$ for $0<-h\ll 1$ and a simple zero ${\overline{h}}_1$ of $M^{+}$ for $0<h\ll 1$. Further taking ${\overline{c}}_i=0(0\le i\le 3)$ and $0<{\overline{c}}_4\ll {\overline{c}}_5$, in this step, there exists a simple zero $h_3$ of $M^{-}$ and a simple zero ${\overline{h}}_2$ of $M^{+}$ for $0<h\ll 1$. Similarly, we can change the sign of ${\overline{c}}_0$, ${\overline{c}}_1$, ${\overline{c}}_2$, ${\overline{c}}_3$ in turn with

$$\begin{array}{c}\hfill 0<-{\overline{c}}_0\ll {\overline{c}}_1\ll {\overline{c}}_2\ll {\overline{c}}_3\ll {\overline{c}}_4.\end{array}$$

This ensures that there exist four more simple zeros of $M^{-}$, $h_7$, $h_6$, $h_5$ and $h_4$; and two simple zeros of $M^{+}$, ${\overline{h}}_4$ and ${\overline{h}}_3$. Thus, if

$$\begin{array}{c}\hfill 0<-{\overline{c}}_0\ll {\overline{c}}_1\ll {\overline{c}}_2\ll {\overline{c}}_3\ll {\overline{c}}_4\ll {\overline{c}}_5\ll {\overline{c}}_6\ll 1,\end{array}$$

then there exist seven simple zeros of $M^{-}$ for $0<-h\ll 1$ and four zeros of $M^{+}$ for $0<h\ll 1$.

Finally, let ${\overline{b}}_0{\overline{b}}_1\left({\delta }_0\right)<0$, $|{\overline{b}}_0|\ll |{\overline{b}}_1\left({\delta }_0\right)|$; there exists a simple zero of $M^{-}$ for $0<h-h_c\ll 1$.

Thus, we obtain 13 zeros of $M^{-}$ and $M^{+}$ altogether which lead to 13 limit cycles for $\epsilon$ sufficiently small.

By using the same method, if ${\overline{b}}_1\left({\delta }_0\right){\overline{c}}_7\left({\delta }_0\right)>0$, we can prove that system (1) can have 12 limit cycles.

Other cases are similar to proof. This ends the proof. □

Theorem 3.

Consider the system (1). Suppose that the conditions of Theorem 1 (ii) are satisfied, and let ${\overline{c}}_i\left(\delta \right)=c_i\left(\delta \right)$, $0\le i\le 3$. If there exist ${\delta }_0\in {\mathbb{R}}^m$ and $8\le l\le 10$ such that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{b}}_0\left({\delta }_0\right)& ={\overline{b}}_1\left({\delta }_0\right)=0,{\overline{c}}_j\left({\delta }_0\right)=0,j=0,\cdots ,l-1,\hfill \\ \hfill \mu & ={(-1)}^{(l+1)}{\overline{b}}_2\left({\delta }_0\right){\overline{c}}_l\left({\delta }_0\right)<0,(resp.,>0),\hfill \\ \hfill \mathrm{rank}& \frac{\partial ({\overline{b}}_0,{\overline{c}}_0,\cdots ,{\overline{c}}_{l-1})}{\partial ({\delta }_1,{\delta }_2,\cdots ,{\delta }_m)}\left({\delta }_0\right)=l+1,\hfill \end{array}\end{array}$$

(37)

then system (1) can have $\left[\frac{3l+9}{2}\right](resp.,\left[\frac{3l+7}{2}\right])$ limit cycles for some $(\epsilon ,\delta )$ near $(0,{\delta }_0)$.

The proof is similar to Theorem 2.

4. Application

Consider the following Liénard cubic oscillator:

$$\dot{x}=y,\dot{y}=(x^2+b)(1-x),$$

where b is a small parameter. Note that the system above is Hamiltonian for all b; there does not exist a limit cycle. However, its portrait varies with the change of b and a cuspidal loop occurs in the case of $b=0$. See Figure 2 below.

As an application of our main results, we consider the following Liénard system:

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \dot{x}& =y,\dot{y}\hfill & =x^2(1-x)-\epsilon \sum _{i\ge 0}^9{a}_i{x}^{i}y,\end{array}\end{array}$$

(38)

and obtain the following theorem.

Theorem 4.

If $a_9\ne 0$, the system (38) can have 12 limit cycles for some $(\epsilon ,a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8)$.

Proof. 

System ${\left(38\right)|}_{\epsilon =0}$ has a nilpotent cusp $O(0,0)$ and an elementary center $C(1,0)$, and has the following Hamiltionian function

$$\begin{array}{c}\hfill H(x,y)=\frac{1}{2}{y}^2-\frac{1}{3}{x}^3+\frac{1}{4}{x}^4.\end{array}$$

It is obvious that

$$\begin{array}{c}\hfill H(0,0)=0,H(1,0)=-\frac{1}{12}.\end{array}$$

From (11) and (14) we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{b}}_0& =-\sum _{i=0}^9{a}_i,c_0=-\sum _{i=0}^9{a}_i{I}_i,c_1=2^{\frac{3}{2}}{3}^{\frac{1}{3}}{a}_0,\hfill \\ \hfill c_2& =-\sum _{i=1}^9{a}_i{J}_i,c_3=-3^{\frac{2}{3}}{2}^{\frac{1}{2}}(a_0+2a_1),\hfill \end{array}\end{array}$$

(39)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill I_i& ={\oint }_{L_0}{x}^{i}ydx=2{\int }_0^{\frac{4}{3}}{x}^{i+\frac{3}{2}}\sqrt{\frac{2}{3}-\frac{1}{2}x}dx,0\le i\le 9,\hfill \\ \hfill J_i& ={\oint }_{L_0}\frac{x^i}{y}dx=2{\int }_0^{\frac{4}{3}}\frac{x^{i-\frac{3}{2}}}{\sqrt{\frac{2}{3}-\frac{1}{2}x}}dx,1\le i\le 9,\hfill \end{array}\end{array}$$

and $L_0$ is given by $y=\pm x^{\frac{3}{2}}\sqrt{\frac{2}{3}-\frac{1}{2}x}$, $0\le x\le \frac{4}{3}$. Simple computions give

$$\begin{array}{c}\hfill \begin{array}{cccccc}\hfill I_0& =\frac{4\sqrt{2}}{27}\pi ,\hfill & I_1& =\frac{10\sqrt{2}}{81}\pi ,& I_2& =\frac{28\sqrt{2}}{243}\pi ,\\ \hfill I_3& =\frac{28\sqrt{2}}{243}\pi ,\hfill & I_4& =\frac{88\sqrt{2}}{729}\pi ,& I_5& =\frac{286\sqrt{2}}{2187}\pi ,\\ \hfill I_6& =\frac{2860\sqrt{2}}{19683}\pi ,\hfill & I_7& =\frac{9724\sqrt{2}}{59049}\pi ,& I_8& =\frac{33592\sqrt{2}}{177147}\pi ,\\ \hfill I_9& =\frac{117572\sqrt{2}}{531441}\pi ,\hfill & J_1& =2\sqrt{2}\pi ,& J_2& =\frac{4\sqrt{2}}{3}\pi ,\\ \hfill J_3& =\frac{4\sqrt{2}}{3}\pi ,\hfill & J_4& =\frac{40\sqrt{2}}{27}\pi ,& J_5& =\frac{140\sqrt{2}}{81}\pi ,\\ \hfill J_6& =\frac{56\sqrt{2}}{27}\pi ,\hfill & J_7& =\frac{616\sqrt{2}}{243}\pi ,& J_8& =\frac{2288\sqrt{2}}{729}\pi ,\\ \hfill J_9& =\frac{2860\sqrt{2}}{729}\pi .\hfill \end{array}\end{array}$$

From ${\overline{b}}_0=c_1=c_3=0$, we obtain that $a_0=a_1=0$, $a_2=-a_3-a_4-a_5-a_6-a_7-a_8-a_9$. Hence, if we let $F(x,y)=-{\displaystyle \sum _{i\ge 0}^9}{a}_i{x}^i$, there exist $P_1(x,y)$ and $Q_1(x,y)$ as follows

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_1(x,y)& =\frac{F(x,0)}{x^2(x-1)}=-\sum _{j=3}^9\sum _{i=j}^9{a}_i{x}^{j-3},\hfill \\ \hfill Q_1(x,y)& =\frac{F(x,y)-F(x,0)}{y}=0,\hfill \end{array}\end{array}$$

such that (16) holds. Then by (23) and (24) we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{b}}_1& =-a_4-3a_5-6a_6-10a_7-15a_8-21a_9,\hfill \\ \hfill {\overline{c}}_4& =2^{\frac{3}{2}}{3}^{\frac{1}{3}}(a_4+a_5+a_6+a_7+a_8+a_9),\hfill \\ \hfill {\overline{c}}_5& =-4\sqrt{2}\pi (a_5+2a_6+\frac{10}{3}{a}_7+\frac{140}{27}{a}_8+\frac{630}{81}{a}_9),\hfill \\ \hfill {\overline{c}}_6& =-2^{\frac{1}{2}}·3^{\frac{2}{3}}(a_4+5a_5+5a_6+5a_7+5a_8+5a_9).\hfill \end{array}\end{array}$$

(40)

From the equations ${\overline{b}}_0\left(\delta \right)=c_0\left(\delta \right)=c_1\left(\delta \right)=c_2\left(\delta \right)=c_3\left(\delta \right)={\overline{c}}_4\left(\delta \right)={\overline{c}}_5\left(\delta \right)=0$, where $\delta =(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9)$, we can solve

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill a_0& =a_1=0,a_2=-a_3,a_4=\frac{303763}{202392}{a}_9,a_8=\frac{35}{14992}{a}_9,\hfill \\ \hfill a_5& =\frac{4}{3}{a}_7+\frac{20865}{7496}{a}_9,a_6=-\frac{7}{3}{a}_7-\frac{2139965}{404784}{a}_9.\hfill \end{array}\end{array}$$

(41)

Let ${\delta }_0=(0,0,-a_3,a_3,\frac{303763}{202392}{a}_9,\frac{4}{3}{a}_7+\frac{20865}{7496}{a}_9,-\frac{7}{3}{a}_7-\frac{2139965}{404784}{a}_9,a_7,\frac{35}{14992}{a}_9,a_9)$. Then ${\overline{b}}_1\left({\delta }_0\right)=\frac{337495a_9}{404784}$, ${\overline{c}}_6\left({\delta }_0\right)=\frac{303763\sqrt{2}·3^{\frac{2}{3}}{a}_9}{50598}$, and it can be seen that

$$\begin{array}{c}\hfill \det \frac{\partial ({\overline{b}}_0,c_0,c_1,c_2,c_3,{\overline{c}}_4,{\overline{c}}_5)}{\partial (a_0,a_1,a_2,a_4,a_5,a_6,a_8)}=-\frac{614072320{\pi }^3}{4782969·3^{\frac{2}{3}}}\ne 0,\end{array}$$

and $\mu =-{\overline{b}}_1\left({\delta }_0\right){\overline{c}}_6\left({\delta }_0\right)<0$ if $a_9\ne 0$. Thus, by Theorem 2 we know that system (38) can have 12 limit cycles for some $(\epsilon ,\delta )$ near $(0,{\delta }_0)$. The proof is completed. □

5. Conclusions and Perspectives

This paper considered the number of limit cycles bifurcating from a cuspidal loop with order 1. By establishing a method based on the idea in [18], we obtained more coefficients of the expansions of the two first order Melnikov functions near the loop. This enabled us to find more limit cycles than [15]. By the same method, we can also obtain more coefficients in the expansions of the first order Melnikov functions as the cusp has order 2 and order 3. The method can also be applied to cases of loops with nilpotent saddles. However, we were not able to study an upper bound of the number of limit cycles in Theorems 2 and 3. It should be possible to give an upper bound for the maximal number of zeros of the Melnikov functions $M^{\pm }$ under the conditions of Theorems 2 or 3. In Theorem 1 we give formulas for computing the coefficients $c_5$, $c_6$ and $c_7$ in (10) if the functions $P_1$ and $Q_1$ satisfying condition (16) exist. In fact, the formulas of these coefficients do not depend on the existence of $P_1$ and on $Q_1$ satisfying condition (16). Then an interesting problem is to give the formulas without assuming (16). For high-dimensional systems, the Melnikov function is a vector function, and it is difficult to study the number of its zeros in this case. However, it is possible to use the expansion of the Melnikov function to study the number of periodic orbits; see [23]. As we knew, the Melnikov method can be also used to predict the occurrence of subharmonic solutions, invariant tori and chaotic orbits in non-autonomous smooth nonlinear systems under periodic perturbation; see [22,24] for example.