\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2014 (2014), No. 95, pp. 1--14.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2014 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2014/95\hfil Bifurcation of limit cycles]
{Bifurcation of limit cycles from quartic isochronous systems}

\author[L. Peng, Z. Feng \hfil EJDE-2014/95\hfilneg]
{Linping Peng, Zhaosheng Feng}  

\address{Linping Peng \newline
School of Mathematics and System Sciences, Beihang University,
LIMB of the Ministry of Education, Beijing, 100191, China}
\email{penglp@buaa.edu.cn, fax (86-10) 8231-7933}

\address{Zhaosheng Feng \newline
Department of Mathematics,
University of Texas-Pan American,
Edinburg, Texas 78539, USA}
\email{zsfeng@utpa.edu}

\thanks{Submitted December 2, 2013. Published April 10, 2014.}
\subjclass[2000]{34C07, 37G15, 34C05}
\keywords{Bifurcation; limit cycles; homogeneous perturbation; 
\hfill\break\indent averaging method; isochronous center; period annulus}

\begin{abstract}
 This article concerns the bifurcation of limit cycles for a quartic
 system with an isochronous center. By using the averaging theory,
 it shows that under any small quartic homogeneous perturbations,
 at most two limit cycles bifurcate from the period annulus of the
 considered system, and this upper bound can be reached.
 In addition, we study a family of perturbed isochronous systems
 and prove that there are at most three limit cycles
 bifurcating from the period annulus of the unperturbed one, and the
 upper bound is sharp.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks

\section{Introduction} \label{sect1}

There has been a longstanding problem, called the Hilbert 16th problem,
whose second part asks for the maximum $H(n)$
of the number of limit cycles and the relative positions for all
planar polynomial differential systems of degree $n$.
One of the most remarkable achievements, Ecalle-Ilyashenko Theorem,
claims that the number of limit cycles is finite for any
individual vector field \cite{BN, I1, E, I2}. However, the existence of
a uniform upper bound for the number even for quadratic vector
fields is still an open problem.


To attack the Hilbert 16th problem, many researchers investigate the
number of limit cycles of various planar polynomial differential
systems. Among them, the problem of the number of limit cycles by perturbing the
periodic orbits of a center has been extensively studied in the
literatures \cite{CLLZ, GI, GGI, I0, LL, LPR, L1} and the
references therein. In general, some useful methods have been proposed based on
the Poincar\'e map \cite{BP, CJ, LLZ}, the
Poincar\'e-Pontryagin-Melnikov integrals or the Abelian
integrals \cite{AI, AN1, AN2, BL2, CGP, GH, XH1, XH2},
the inverse integrating factor \cite{GLV1, GLV2, GLV3, VLG},
and the averaging method which is equivalent to the Abelian integrals
in the plane \cite{ BL1, CLP, GL, LL, LPR, L1}.

Although in the plane the methods based on the Abelian integrals and
the averaging theory are equivalent, each has its own advantages.
For example, when the associated Abelian integrals are complicated
or we need to study
the periodic orbits of the non-autonomous differential systems, the
averaging method displays more flexibility. Roughly speaking, the
averaging method gives a quantitative relation between the solutions
of a non-autonomous periodic differential system and the solutions
of its averaged differential system, which is autonomous. In
particular, for the averaging method of the first order, the number
of hyperbolic equilibrium points of the averaged differential system
can give a lower bound of the maximal number of limit cycles of the
non-autonomous periodic differential system \cite{SV, V}.

As mentioned above, by using the averaging method, the problem on
the number of limit cycles of the
non-autonomous periodic differential systems is equivalent to
the exploration of
the number of hyperbolic equilibrium points of the averaged
differential systems. Hence, the averaging theory has played a
crucial role in the study of limit cycles of differential
systems. Now there are quite many important results on
the number of limit cycles of the polynomial differential systems
by the averaging method, such as Llibre \cite{L1}, Buic\u a
and Llibre \cite{BL1, BL2}, Gine and Llibre \cite{GL} and so on. It
seems that among these results,  more are focused on
differential systems of lower degree. As far as we know, for the integrable
systems of higher degree, in some cases the first integrals may have
complicated expressions so that it is out of the reach to study the
bifurcation of limit cycles of these systems under small
perturbations.

In this article, we consider the  quartic system
\begin{equation}\label{eq1}
\begin{gathered}
\dot{x}=-y+x^3y+xy^3,\\
\dot{y}=x+x^2y^2+y^4,
\end{gathered}
\end{equation}
which has
$$
H(x,y)=\frac{1}{3(x^2+y^2)^{3/2}}-\frac{x}{(x^2+y^2)^{1/2}}=c
$$ 
as its first integral with the integrating factor $1/(x^2+y^2)^{5/2}$ 
and the unique finite singularity $(0, 0)$
as its isochronous center. The period annulus, denoted by
$$
\{(x, y)|H(x,y)=c, \; c\in(1, +\infty)\}
$$
starts at the center $(0, 0)$ and
terminates with the separatrix
passing the infinite degenerate singularity on the equator. The phase portrait of system \eqref{eq1} is shown in Fig.1.

 \begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1}
\end{center}
\caption{The  phase portrait of system \eqref{eq1} in the Poincar\'{e} disk.}
 \label{fig1}
 \end{figure}

By using the averaging method, we study the bifurcation of limit
cycles from system \eqref{eq1} under any small perturbations, and
prove the following main results.

\begin{theorem}\label{thm1}
For any sufficiently small parameter $|\varepsilon|$, and any real
constants $a_{ij}$ and $b_{ij}$ $(i, j=0, 1, 2, 3, 4)$, the following
quartic perturbation of system \eqref{eq1},
\begin{equation}\label{eq2}
\begin{gathered}
\dot{x}=-y+x^3y+xy^3+\varepsilon\sum_{i+j=4}a_{ij}x^iy^j,\\
\dot{y}=x+x^2y^2+y^4+\varepsilon\sum_{i+j=4}b_{ij}x^iy^j,
\end{gathered}
\end{equation}
has at most two limit cycles bifurcating from the period annulus
around the center $(0, 0)$ of the unperturbed one, and this upper
bound is sharp.
\end{theorem}

\begin{theorem}\label{thm2}
For the family of quartic perturbations
\begin{equation}\label{eq3}
\begin{split}
\dot{x}&=-y+x^3y+xy^3+\varepsilon(a_{10}x+a_{01}y+a_{11}xy+a_{21}x^2y+a_{03}y^3\\
&\quad +a_{40}x^4+a_{31}x^3y+a_{22}x^2y^2+a_{13}xy^3+a_{04}y^4),\\
\dot{y}&=x+x^2y^2+y^4+\varepsilon(b_{10}x+b_{01}y+b_{20}x^2+b_{02}y^2+b_{30}x^3+b_{12}xy^2\\
&\quad +b_{40}x^4+b_{31}x^3y+b_{22}x^2y^2+b_{13}xy^3+b_{04}y^4),
\end{split}
\end{equation}
where $|\varepsilon|$ is sufficiently small, $a_{i, j}$ and $b_{i,
j}$ $(i, j=0, 1, 2, 3, 4)$ are any real constants.  Then there are at most three
limit cycles bifurcating from the period annulus surrounding the
center $(0, 0)$ of the unperturbed system, and this upper bound is
sharp.
\end{theorem}

The rest of this paper is organized as follows. In Section
$2$, we give an introduction on the averaging theory,
including some technical lemmas and methods employed in the
averaging theory. Section $3$ is dedicated to the proof of Theorem
\ref{thm1} by computing the averaged equations corresponding to the
equivalent system of system \eqref{eq2} and exploring the number of
its hyperbolic equilibriums. In Section $4$, after making a
transformation to system \eqref{eq3}, theorem \ref{thm2} is proven
through analyzing an equivalent system
and a corresponding averaged system. In addition,
some examples are illustrated to verify the obtained results.

\section{Preliminary results} \label{sect2}

In this section, we introduce some preliminary results on the averaging theory
that will be used in our quartic polynomial systems.

The following lemma provides a first order approximation for the periodic
solution of a periodic differential equation. For the proof, we refer the reader to
\cite[Theorem 2.6.1]{SV} and \cite[Theorems 11.5 and 11.6]{V}.

\begin{lemma}\label{lem1}
Consider the two initial value problems
\begin{equation}\label{eq4}
\dot{x}=\varepsilon f(t, x)+\varepsilon^2h(t, x, \varepsilon),
\quad x(0)=x_0,
\end{equation} 
and
\begin{equation}\label{eq5}
\dot{y}=\varepsilon f^0(y), \quad y(0)=x_0,
\end{equation}
 where $x, y, x_0\in D$, here $D$ is an open subset of 
$R, t\in[0, +\infty), \varepsilon\in(0, \varepsilon_0], f$ and $h$ 
are periodic with period $T$ in $t$, and
\begin{equation}\label{eq6}
f^0(y)=\frac{1}{T}\int_0^Tf(t, y)dt.
\end{equation}
We suppose that
\begin{enumerate}
\item $f$, $\partial f/\partial x$, $\partial^2f/\partial x^2$ and
$\partial h/\partial x$ are continuous and bounded by a
constant independent on $\varepsilon$ in $[0, +\infty)\times D$ and
$\varepsilon\in(0, \varepsilon_0]$;

\item
$T$ is independent on $\varepsilon$; and

\item
$y(t)$ belongs to $D$ on the time-scale $1/\varepsilon$.
\end{enumerate}
Then the following statements hold.
\begin{itemize}
\item[(a)] On the time-scale $1/\varepsilon$, we have that 
$$
x(t)-y(t)=O(\varepsilon), \quad\text{as }\varepsilon\to 0.
$$
\item[(b)] If $p$ is an equilibrium point of the averaged system \eqref{eq5}
such that
\begin{equation}\label{eq7}
(df^0/dy)(p)\neq 0, 
\end{equation} 
then there exists a $T$-periodic solution $\phi(t, \varepsilon)$ of 
equation \eqref{eq4}
which is close to $p$ such that $\phi(t, \varepsilon)\to p$
as $\varepsilon\to 0$.

\item[(c)] If \eqref{eq7} is negative, then the corresponding periodic
solution $\phi(t, \varepsilon)$ in the plane $(t, x)$ is
asymptotically stable for any sufficiently small $|\varepsilon|$. If
\eqref{eq7} is positive, then it is unstable.
\end{itemize}
\end{lemma}

Let us consider another integrable system of the form
\begin{equation}\label{eq8}
\begin{gathered}
\dot{x}=P(x, y),\\
\dot{y}=Q(x, y),
\end{gathered}
\end{equation}
with a first integral $H$ and  a continuous family of ovals
$$
\{\gamma_{h}\}\subset\{(x, y)|H(x, y)=h, h_1<h<h_2\}.
$$
We consider a perturbed system:
\begin{equation}\label{eq9}
\begin{gathered}
\dot{x}=P(x, y)+\varepsilon p(x, y),\\
\dot{y}=Q(x, y)+\varepsilon q(x, y).
\end{gathered}
\end{equation}

To study the number of limit cycles for any sufficiently
small $|\varepsilon|$ by using the above averaging theory, we
need to transform system \eqref{eq9} to the canonical form in Lemma
\ref{lem1}. The following lemma \cite{BL1} provides us a useful
transformation.

\begin{lemma}\label{lem2}
For system \eqref{eq8}, assume $xQ(x, y)-yP(x, y)\neq 0$ for all
$(x, y)$ in the period annulus formed by the ovals $\gamma_{h}$. Let
$$
\rho: (\sqrt{h_1}, \sqrt{h_2})\times [0, 2\pi)\to [0,+\infty)
$$
be a continuous function such that
$$
H(\rho(R, \varphi)\cos\varphi, \rho(R, \varphi)\sin\varphi)=R^2,
$$
for all $R\in(\sqrt{h_1}, \sqrt{h_2})$ and
$\varphi \in[0,2\pi)$. Then the differential equation which describes the
dependence between the square root of energy, $R=\sqrt{h}$, and the
angle $\varphi$ for system \eqref{eq9} is
\begin{equation}\label{eq90}
\frac{dR}{d\varphi}=\varepsilon\frac{\mu(x^2+y^2)(Qp-Pq)}{2R(Qx-Py)}
\Big(1-\varepsilon\frac{qx-py}{Qx-Py}\Big)+O(\varepsilon^3),
\end{equation}
where $x=\rho(R, \varphi)\cos\varphi$ and $y=\rho(R,
\varphi)\sin\varphi$.
\end{lemma}


The following lemma presents the version of the formula of the first
order Melnikov integral associated with system \eqref{eq9} in the
polar coordinates \cite{BL1}.

\begin{lemma}\label{lem3}
Under the conditions of Lemma \ref{lem2}, we define
\begin{equation}\label{eq91}
\begin{gathered}
d(R,\varepsilon)=\int_0^{2\pi}\Big[\varepsilon\frac{\mu(x^2+y^2)(Qp-Pq)}{2R(Qx-Py)}
\Big(1-\varepsilon\frac{qx-py}{Qx-Py}\Big)+O(\varepsilon^3)\Big]d\varphi,\\
M_1(R)=\int_0^{2\pi}\frac{\mu(x^2+y^2)(Qp-Pq)}{2R(Qx-Py)}d\varphi,
\end{gathered}
\end{equation}
for system \eqref{eq9}, where $\mu=\mu(x, y)$ is the integrating
factor of system \eqref{eq8} corresponding to the first integral
$H$, and $x=\rho\cos\varphi$ and $y=\rho\sin\varphi$.
Then $d(R, \varepsilon)$ and $M_1(R)$ expressed by \eqref{eq91} are the
displacement function and the first order Melnikov integral of
system \eqref{eq9}, respectively.
\end{lemma}

Based on Lemmas \ref{lem1}, \ref{lem2} and \ref{lem3}, we can obtain

\begin{corollary}\label{coro1}
If $d^0(R)$ represents the averaged function of the first
approximation in $\varepsilon$ of the right side of system
\eqref{eq90}, then the following relation holds,
\[
2\pi d^0(R)=M_1(R),
\]
where $M_1(R)$ is defined by \eqref{eq91}.
\end{corollary}


Corollary \ref{coro1} provides a relation between the averaged
function and the first order Melnikov integral associated with the
same differential system, which enables us to
explore the maximal number of limit cycles of system \eqref{eq9}
bifurcating from the period annulus of system \eqref{eq8} via the
averaging method.


\section{Proof of Theorem \ref{thm1}} \label{sect3}

For 
$$
 H(x, y)=\frac{1}{3(x^2+y^2)^{3/2}}-\frac{x}{(x^2+y^2)^{1/2}},
$$
we choose the function 
\begin{equation}\label{eq10}
\begin{split}
\rho(R, \varphi)=\frac{1}{(R^2+3\cos\varphi)^{1/3}}
\end{split}
\end{equation}
 such
that $ H(\rho\cos\varphi, \rho\sin\varphi)=R^2/3$.
Let
\begin{equation}\label{eq11}
\begin{gathered}
x=\rho(R, \varphi)\cos\varphi,\\
y=\rho(R, \varphi)\sin\varphi,
\end{gathered}
\end{equation}
for $\varphi\in[0, 2\pi)$ and $R>\sqrt{3}$. By using Lemma \ref{lem2}, we
can transform system \eqref{eq2} as
\begin{equation}\label{eq12}
\begin{aligned}
\frac{dR}{d\varphi}
&=\Big(\varepsilon\frac{3(Qp-Pq)}{2R(x^2+y^2)^{5/2}}
-\varepsilon^2\frac{3(Qp-Pq)(qx-py)}{2R(x^2+y^2)^{7/2}}\Big)\Big|_{x=\rho(R,
\varphi)\cos\varphi,  y=\rho(R, \varphi)\sin\varphi}\\
&\quad +O(\varepsilon^3),
\end{aligned}
\end{equation}
where
\begin{align*}
Qp-Pq&=-b_{40}x^{7}y+(a_{40}-b_{31})x^{6}y^2+(a_{31}-b_{40}-b_{22})x^{5}y^3\\
&\quad +(a_{40}+a_{22}-b_{31}-b_{13})x^4y^4+(a_{31}+a_{13}-b_{22}-b_{04})x^3y^{5}\\
&\quad +(a_{22}+a_{04}-b_{13})x^2y^{6}+(a_{13}-b_{04})xy^{7}+a_{04}y^{8}\\
&\quad +a_{40}x^{5}+(a_{31}+b_{40})x^4y+(a_{22}+b_{31})x^3y^2\\
&\quad +(a_{13}+b_{22})x^2y^3+(a_{04}+b_{13})xy^4+b_{04}y^{5},
\end{align*}
\begin{align*}
qx-py&=b_{40}x^4+(b_{31}-a_{40})x^4y+(b_{22}-a_{31})x^3y^2+(b_{13}-a_{22})x^2y^3\\
&\quad +(b_{04}-a_{13})xy^4-a_{04}y^{5}.
\end{align*}
The averaged equation corresponding to system \eqref{eq12} is
\begin{equation}\label{eq13}
\dot{R}=\varepsilon f^0(R),
\end{equation}
where
\begin{equation}
\begin{split}
f^0(R)&=\frac{1}{2\pi}\int_0^{2\pi}\Big(\frac{3(Qp-Pq)}{2R(x^2+y^2)^{5/2}}\Big)
 \Big|_{x=\rho(R,\varphi)\cos\varphi, y=\rho(R, \varphi)\sin\varphi}d\varphi\\
&=\frac{1}{4\pi
R}\Big[M_1\int_0^{2\pi}\frac{\cos^{6}\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi
+M_2\int_0^{2\pi}\frac{\cos^4\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi\\
&\quad +M_3\int_0^{2\pi}
\frac{\cos^2\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi
-(M_1+M_2+M_3)\int_0^{2\pi}
\frac{1}{\cos\varphi+\frac{R^2}{3}}d\varphi\Big],
\end{split}
\end{equation}
and
\begin{equation}\label{eq13*}
 \begin{gathered}
 M_1=-a_{40}+a_{22}-a_{04}+b_{31}-b_{13},\\
 M_2=a_{40}-2a_{22}+3a_{04}-b_{31}+2b_{13},\\
 M_3=a_{22}-3a_{04}-b_{13}.
 \end{gathered}
 \end{equation}
Straightforward computations give
\begin{gather*}
\int_0^{2\pi}\frac{\cos^{6}\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi=-\frac{\pi
R^2}{4}-\frac{\pi R^{6}}{27}-\frac{2\pi R^{10}}{243}+\frac{2\pi
R^{12}}{243\sqrt{R^4-9}},\\
\int_0^{2\pi}\frac{\cos^4\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi=-\frac{\pi
R^2}{3}-\frac{2\pi R^{6}}{27}+\frac{2\pi R^{8}}{27\sqrt{R^4-9}},\\
\int_0^{2\pi}\frac{\cos^2\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi=-\frac{2\pi
R^2}{3}+\frac{2\pi R^4}{3\sqrt{R^4-9}}.
\end{gather*}
From these expressions, we obtain
\begin{equation}\label{eq14}
\begin{split}
f^0(R)
&=\frac{1}{4R}\Big\{\Big[-\frac{2M_1}{243}R^{10}-\frac{M_1+2M_2}{27}
R^{6}-\frac{3M_1+4M_2+8M_3}{12}R^2\Big]\\
&\quad +\Big[\frac{2M_1}{243}R^{12}+\frac{2M_2}{27}R^{8}
 +\frac{2M_3}{3}R^4-6(M_1+M_2+M_3)\Big]\frac{1}{\sqrt{R^4-9}}\Big\}\\
&=\frac{1}{4R}\Big\{\Big[-\frac{2M_1}{243}S^{5}-\frac{M_1+2M_2}{27}
S^3-\frac{3M_1+4M_2+8M_3}{12}S\Big]\\
&\quad +\Big[\frac{2M_1}{243}S^{6}+\frac{2M_2}{27}S^4
 +\frac{2M_3}{3}S^2-6(M_1+M_2+M_3)\Big]\frac{1}{\sqrt{S^2-9}}\Big\},
\end{split}
\end{equation}
where $S=R^2$.
Let 
$$
S=\frac{3(1+w^2)}{1-w^2}.
$$ 
For $0<w<1$, formula \eqref{eq14} becomes
\begin{align*}
f^0(R)&=\frac{(w-1)}{16R(w+1)^{5}} g(w)\\
&=-\frac{\sqrt{3}(1-w)^{3/2}}{48(1+w^2)^{1/2}(1+w)^{9/2}}
[N_1w^4+N_2w^3+N_3w^2+N_2w+N_1],
\end{align*}
where
\begin{gather*}
g(w)=N_1w^4+N_2w^3+N_3w^2+N_2w+N_1, \\
N_1=15M_1+12M_2+8M_3,\quad N_2=42M_1+40M_2+32M_3,\\
N_3=62M_1+56M_2+48M_3.
\end{gather*}

As a result of the symmetry of coefficients of $g(w)$, we know
that if $w_0\neq 0$ is one root of $g(w)=0$, so is $1/w_0$.
Hence, the fact that $g(w)$ has at most two zeros in $w\in(0, 1)$
implies that there exist at most two zeros for $f^0(R)$ in
$R\in(\sqrt{3}, +\infty)$. By  Lemma \ref{lem1} and Corollary
\ref{coro1}, we get that system \eqref{eq12} has at most two periodic
solutions which tend to the corresponding hyperbolic equilibriums,
respectively. That is, for system \eqref{eq2} with any sufficiently
small $|\varepsilon|$, at most two limit cycles bifurcate from the
period annulus around the center $(0, 0)$ of system \eqref{eq1}.

In fact, there exist  many systems expressed like \eqref{eq2} which
have exactly two limit cycles emerging from the period annulus
 of the unperturbed system. In the following, we not only provide
some examples satisfying this property, but also introduce a method
to construct such systems.

Suppose that
\begin{align*}
\tilde{g}(w)&=\big(w-\frac{1}{10}\big)\big(w-\frac{1}{5}\big)(w-10)(w-5)\\
&=w^4-\frac{153}{10}w^3+\frac{1363}{25}w^2-\frac{153}{10}w+1.
\end{align*}
Take the constants
\[
C_1=1,\quad C_2=-\frac{153}{10},\quad C_3=\frac{1363}{25},
\]
then we can choose 
\begin{equation}\label{eq15}
M_1=\frac{1089}{100}, \quad M_2=-\frac{2209}{100},\quad M_3=\frac{10273}{800},
\end{equation}
such that 
\begin{gather*} 
15M_1+12M_2+8M_3=1,\\
42M_1+40M_2+32M_3=-\frac{153}{10},\\
62M_1+56M_2+48M_3=\frac{1363}{25}.
\end{gather*}
From \eqref{eq13*} and \eqref{eq15}, we have
$$
a_{40}=b_{31}-\frac{213}{160},\quad
a_{22}=b_{13}+\frac{3167}{400},\quad
a_{04}=-\frac{1313}{800}.
$$
Hence, for the sufficiently small $|\varepsilon|$, we obtain a family of systems
\begin{equation}\label{eq16}
\begin{gathered}
\begin{aligned}
\dot{x}
&=-y+x^3y+xy^3+\varepsilon\Big[\big(b_{31}-\frac{213}{160}\big)x^4+a_{31}x^3y\\
&\quad +\big(b_{13}+\frac{3167}{400}\big)x^2y^2+a_{13}xy^3-\frac{1313}{800}y^4\Big],
\end{aligned}\\
\dot{y}=x+x^2y^2+y^4+\varepsilon
\big[b_{40}x^4+b_{31}x^3y+b_{22}x^2y^2+b_{13}xy^3+b_{04}y^4\big],
\end{gathered}
\end{equation}
where $a_{13},\ a_{31}$ and $b_{ij}\ (i, j=0, 1, 2, 3, 4)$ are any
real constants.

By using polar coordinates $x=\rho\cos\varphi$ and $y=\rho\sin\varphi$,
system \eqref{eq16} can be rewritten as
\begin{equation}\label{eq17}
\frac{dR}{d\varphi}=\varepsilon F(\varphi, R)+O(\varepsilon^2),
\end{equation}
where
\begin{align*}
F(\varphi,R)
&=\rho^3\Big[-b_{40}\cos^{7}\varphi\sin\varphi-\frac{213}{160}\cos^{6}\varphi\sin^2\varphi+\left(a_{31}-b_{40}-b_{22}\right)
\cos^{5}\varphi\sin^3\varphi\\
&\quad +\frac{5269}{800}\cos^4\varphi\sin^4\varphi+\left(a_{31}+a_{13}-b_{22}-b_{04}\right)\cos^3\varphi\sin^{5}\varphi\\
&\quad +\frac{5021}{800}\cos^2\varphi\sin^{6}\varphi+\left(a_{13}-b_{04}\right)\cos\varphi\sin^{7}\varphi-\frac{1313}{800}sin^{8}\varphi\Big]\\
&\quad +\Big[\big(b_{31}-\frac{213}{160}\big)\cos^{5}\varphi+(a_{31}+b_{40})\cos^4\varphi\sin\varphi\\
&\quad +\big(b_{31}+b_{13}+\frac{3167}{400}\big)\cos^3\varphi\sin^2\varphi
+\big(a_{13}+b_{22}\big)\cos^2\varphi\sin^3\varphi\\
&\quad +\big(b_{13}-\frac{1313}{800}\big)\cos\varphi\sin^4\varphi
 +b_{04}\sin^{5}\varphi\Big].
\end{align*}
The averaged equation of system \eqref{eq17} is given by
\begin{equation}\label{eq18}
\frac{dR}{d\varphi}=\varepsilon f_{*}^0(R)+O(\varepsilon^2),
\end{equation}
where 
\begin{equation}\label{eq19}
\begin{split}
f_{*}^0(R)
&=\frac{1}{2\pi}\int_0^{2\pi}F(\varphi, R)d\varphi\\
&=\frac{1}{4\pi
R}\Big[\frac{1089}{100}\int_0^{2\pi}
 \frac{\cos^{6}\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi
-\frac{2209}{100}\int_0^{2\pi}\frac{\cos^4\varphi}{\cos\varphi
 +\frac{R^2}{3}}d\varphi\\
&\quad +\frac{10273}{800}\int_0^{2\pi}\frac{\cos^2\varphi}{\cos\varphi+\frac{R^2}{3}}d\varphi
-\frac{1313}{800}\int_0^{2\pi}\frac{1}{\cos\varphi+\frac{R^2}{3}}d\varphi\Big]\\
&=-\frac{\sqrt{3}(1-w)^{3/2}}{48(1+w^2)^{1/2}(1+w)^{9/2}}
 (w-\frac{1}{10})(w-\frac{1}{5})(w-10)(w-5).
\end{split}
\end{equation}
Apparently, $f_{*}^0(R)$ has exactly two positive zeros, denoted
by
\[
R_1=\frac{\sqrt{29997}}{99}\approx1.749458791,\quad
R_2=\frac{\sqrt{1872}}{24}\approx1.802775638,
\]
corresponding to
$w_1=1/10$ and $w_2=1/5$ in $R\in(\sqrt{3}, +\infty)$.
Moreover, we have
\begin{gather*}
\frac{df_{*}^0(R_1)}{dR}=\frac{107163}{387200}\approx0.5260835926>0,\\
\frac{df_{*}^0(R_2)}{dR}=-\frac{49}{675}\approx-0.07259259259<0.
\end{gather*}

It follows from Lemma \ref{lem1} and Corollary \ref{coro1} that for the
sufficiently small $|\varepsilon|$, system \eqref{eq16} has just two
limit cycles emerging from the period annulus of the corresponding
unperturbed system: one is unstable and the other is stable.
This completes the proof of Theorem \ref{thm1}.

As a byproduct, we obtain

\begin{theorem} \label{thm3}
For the sufficiently small $|\varepsilon|$, system \eqref{eq17} has exactly
two periodic solutions, denoted  by $l_1$ and $l_2$ respectively,
such that $l_1$ shrinks to $R_1$ and $l_2$ shrinks to $R_2$ as
$\varepsilon$ goes to $0$.
Moreover,  $l_1$ is unstable while $l_2$ is stable.
\end{theorem}

\section{Proof of Theorem \ref{thm2}} \label{sect4}

After using the transformation \eqref{eq11}, system \eqref{eq3} can
be re-expressed as
\begin{equation}\label{eq20}
\frac{dR}{d\varphi}=\varepsilon
\Big(\frac{3}{2R}\frac{Q\tilde{p}-P\tilde{q}}{\rho^{5}}\Big)
\Big|_{x=\rho\cos\varphi, y=\rho\sin\varphi}+O(\varepsilon^2),
\end{equation}
where $\rho$ is defined as \eqref{eq10}, and
\begin{align*}
&Q\tilde{p}-P\tilde{q}\\
&=[a_{10}x^2+(a_{01}+b_{10})xy+b_{01}y^2]+[(a_{11}+b_{20})x^2y+b_{02}y^3]\\
&\quad +[(a_{21}+b_{30})x^3y+(a_{03}+b_{12})xy^3]+[a_{40}x^{5}+(a_{31}+b_{40}-b_{10})x^4y\\
&\quad +(a_{22}+a_{10}+b_{31}-b_{01})x^3y^2+(a_{13}+a_{01}+b_{22}-b_{10})x^2y^3\\
&\quad +(a_{04}+a_{10}+b_{13}-b_{01})xy^4+(a_{01}+b_{04})y^{5}]\\
&\quad +[-b_{20}x^{5}y+(a_{11}-b_{20}-b_{02})x^3y^3+(a_{11}-b_{02})xy^{5}]\\
&\quad +[-b_{30}x^{6}y+(a_{21}-b_{30}-b_{12})x^4y^3+(a_{21}+a_{03}-b_{12})x^2y^{5}+a_{03}y^{7}]\\
&\quad +[-b_{40}x^{7}y+(a_{40}-b_{31})x^{6}y^2+(a_{31}-b_{40}-b_{22})x^{5}y^3\\
&\quad +(a_{40}+a_{22}-b_{31}-b_{13})x^4y^4+(a_{31}+a_{13}-b_{22}-b_{04})x^3y^{5}\\
&\quad +(a_{22}+a_{04}-b_{13})x^2y^{6}
   +(a_{13}-b_{04})xy^{7}+a_{04}y^{8}].
\end{align*}
The averaged equation associated with system \eqref{eq20} is
\begin{equation}\label{eq21}
\frac{dR}{d\varphi}=\varepsilon g^0(R)+O(\varepsilon^2),
\end{equation}
where
\begin{equation}\label{eq21*}
\begin{split}
g^0(R)
&=\frac{1}{2\pi}\int_0^{2\pi}\frac{3}{2R}
\big(\frac{Q\tilde{p}-P\tilde{q}}{\rho^{5}}\big)\Big|_{x=\rho\cos\varphi,
y=\rho\sin\varphi}d\varphi\\
&=\frac{3}{4\pi R}\int_0^{2\pi}\Big\{\frac{a_{10}\cos^2\varphi
 +b_{01}\sin^2\varphi}{\rho^3}\\
&\quad +\Big[a_{40}\cos^{5}\varphi+(a_{22}+a_{10}
 +b_{31}-b_{01})\cos^3\varphi\sin^2\varphi\\
&\quad +(a_{04}+a_{10}+b_{13}-b_{01})\cos\varphi\sin^4\varphi\Big]\\
&\quad +\rho^3\Big[(a_{40}-b_{31})\cos^{6}\varphi\sin^2\varphi+(a_{40}
 +a_{22}-b_{31}-b_{13})\cos^4\varphi\sin^4\varphi\\
&\quad +(a_{22}+a_{04}-b_{13})\cos^2\varphi\sin^{6}\varphi
 +a_{04}\sin^{8}\varphi\Big]\Big\}d\varphi.
\end{split}
\end{equation}
Using a similar transformation as in the preceding section to \eqref{eq21*}, the
function $g^0(R)$ can be simplified as
\begin{equation}\label{eq22}
\begin{split}
g^0(R)
&=\frac{3}{4R}\Big[-\frac{2M_1}{729}R^{10}-\frac{M_1+2M_2}{81}R^{6}
+\Big(-\frac{3M_1+4M_2+8M_3}{36}+M_4\Big) R^2\\
&\quad +\Big(\frac{2M_1}{729}R^{12}+\frac{2M_2}{81}R^{8}
 +\frac{2M_3}{9}R^4-2(M_1+M_2+M_3)\Big)\frac{1}{\sqrt{R^4-9}}\Big]\\
&=\frac{3}{4R}\Big[-\frac{2M_1}{729}S^{5}-\frac{M_1+2M_2}{81}S^3
 +\Big(-\frac{3M_1+4M_2+8M_3}{36}+M_4\Big) S\\
&\quad +\Big(\frac{2M_1}{729}S^{6}+\frac{2M_2}{81}S^4
+\frac{2M_3}{9}S^2-2(M_1+M_2+M_3)
 \Big)\frac{1}{\sqrt{S^2-9}}\Big]\\
&=-\frac{\sqrt{3}}{48(1-w^2)^{1/2}(1+w^2)^{1/2}(1+w)^4}\\
&\quad\times [\tilde{N}_1w^{6}+\tilde{N}_2w^{5}+\tilde{N}_3w^4 
+\tilde{N}_4w^3+\tilde{N}_3w^2+\tilde{N}_2w+\tilde{N}_1],
\end{split}
\end{equation}
where $M_{i}$ $(i=1, 2, 3)$ are defined as \eqref{eq13*}, and
\begin{gather*}
M_4=a_{10}+b_{01},\\
\tilde{N}_1=15M_1+12M_2+8M_3-36M_4,\\
\tilde{N}_2=12M_1+16M_2+16M_3-144M_4,\\
\tilde{N}_3=-7M_1-12M_2-8M_3-252M_4,\\
\tilde{N}_4=-40M_1-32M_2-32M_3-288M_4.
\end{gather*}
Similarly, from \eqref{eq22}, we get that $g^0(R)$ has at most
three zeros in $R\in(\sqrt{3}, +\infty)$. Using this fact together with
Lemma \ref{lem1} and Corollary \ref{coro1}, it follows that system
\eqref{eq20} has at most three periodic solutions tending to the
corresponding hyperbolic equilibriums, respectively. This means that
the maximal number of limit cycles of system \eqref{eq3} emerging
from the period annulus of the unperturbed one is three. Moreover,
the upper bound can be reached.

As an example, we consider the system
\begin{equation}\label{eq23}
\begin{split}
\dot{x}
&=-y+x^3y+xy^3+\varepsilon\Big[\big(-b_{01}+\frac{9}{800}\big)x
 +a_{01}y+a_{11}xy+a_{21}x^2y +a_{03}y^3\\
&\quad +\big(b_{31}-\frac{109}{80}\big)x^4+a_{31}x^3y
 +\big(b_{13}+\frac{28279}{3200}\big)x^2y^2
 +a_{13}xy^3-\frac{1313}{640}y^4\Big],\\
\dot{y}
&=x+x^2y^2+y^4+\varepsilon\Big[b_{10}x+b_{01}y+b_{20}x^2+b_{02}y^2+b_{30}x^3+b_{12}xy^2+b_{40}x^4\\
&\quad +b_{31}x^3y+b_{22}x^2y^2+b_{13}xy^3+b_{04}y^4\Big],
\end{split}
\end{equation}
where $|\varepsilon|$ is sufficiently small, $a_{ij}$ 
$(i=0, 1, 2, 3, j=1, 3)$ and $b_{ij}\ (i, j=0, 1, 2, 3, 4)$ are any 
real constants.

By polar coordinates in \eqref{eq11}, system \eqref{eq23} is
equivalent to
\begin{equation}\label{eq24}
\frac{dR}{d\varphi}=\varepsilon G(\varphi, R)+O(\varepsilon^2),
\end{equation}
where
\begin{align*}
&G(\varphi, R)\\
&=\frac{\big(-b_{01}+\frac{9}{800}\big)\cos^2\varphi+(a_{01}+b_{10})\cos\varphi\sin\varphi+b_{01}\sin^2\varphi}{\rho^3}\\
&\quad+\frac{(a_{11}+b_{20})\cos^2\varphi\sin\varphi+b_{02}\sin^3\varphi}{\rho^2}\\
&\quad +\frac{(a_{21}+b_{30})\cos^3\varphi\sin\varphi+(a_{03}+b_{12})\cos\varphi\sin^3\varphi}{\rho}\\
&\quad +\Big[\big(b_{31}-\frac{109}{80}\big)\cos^{5}\varphi+(a_{31}+b_{40}-b_{10})\cos^4\varphi\sin\varphi\\
&\quad +\big(b_{31}+b_{13}-2b_{01}+\frac{5663}{640}\big)\cos^3\varphi\sin^2\varphi\\
&\quad +(a_{13}+a_{01}+b_{22}-b_{10})\cos^2\varphi\sin^3\varphi\\
&\quad +\big(b_{13}-2b_{01}-\frac{6529}{3200}\big)\cos\varphi\sin^4\varphi
 +(a_{01}+b_{04})\sin^{5}\varphi \Big]\\
&\quad +\rho\Big[-b_{20}\cos^{5}\varphi\sin\varphi
 +(a_{11}-b_{02})\cos\varphi\sin^{5}\varphi\\
&\quad +(a_{11}-b_{20}-b_{02})\cos^3\varphi\sin^3\varphi\Big]\\
&\quad +\rho^2\Big[-b_{30}\cos^{6}\varphi\sin\varphi
 +(a_{21}-b_{30}-b_{12})\cos^4\varphi\sin^3\varphi\\
&\quad +(a_{21}+a_{03}-b_{12})\cos^2\varphi\sin^{5}\varphi
 +a_{03}\sin^{7}\varphi\Big]\\
&\quad +\rho^3\Big[-b_{40}\cos^{7}\varphi\sin\varphi
 -\frac{109}{80}\cos^{6}\varphi\sin^2\varphi+(a_{31}
 -b_{40}-b_{22})\cos^{5}\varphi\sin^3\varphi\\
&\quad +\frac{23919}{3200}\cos^4\varphi\sin^4\varphi
+(a_{31}+a_{13}-b_{22}-b_{04})\cos^3\varphi\sin^{5}\varphi\\
&\quad +\frac{10857}{1600}\cos^2\varphi\sin^{6}\varphi
+(a_{13}-b_{04})\cos\varphi\sin^{7}\varphi-\frac{1313}{640}\sin^{8}\varphi\Big].
\end{align*}
The averaged equation of system \eqref{eq24} is
\begin{equation}\label{eq25}
\frac{dR}{d\varphi}=\varepsilon g_{*}^0(R)+O(\varepsilon^2),
\end{equation}
where
\begin{equation}
\begin{split}
g_{*}^0(R)
&=\frac{1}{2\pi}\int_0^{2\pi}G(\varphi, R)d\varphi\\
&=\frac{3}{4\pi R}\int_0^{2\pi}\Big\{\rho^3
 \Big[-\frac{109}{80}\cos^{6}\varphi\sin^2\varphi+\frac{23919}{3200}\cos^4\varphi\sin^4\varphi\\
&\quad +\frac{10857}{1600}\cos^2\varphi\sin^{6}\varphi-\frac{1313}{640}\sin^{8}\varphi\Big]\\
&\quad +\Big[\big(b_{31}-\frac{109}{80}\big)\cos^{5}\varphi
 +\big(b_{31}+b_{13}-2b_{01}+\frac{5663}{640}\big)\cos^3\varphi\sin^2\varphi\\
&\quad +\big(b_{13}-2b_{01}-\frac{6529}{3200}\big)\cos\varphi\sin^4\varphi\Big]\\
&\quad +\frac{1}{\rho^3}\Big[\big(-b_{01}+\frac{9}{800}\big)
 \cos^2\varphi+b_{01}\sin^2\varphi\Big]\Big\}d\varphi\\
&=-\frac{\sqrt{3}}{48(1-w^2)^{1/2}(1+w^2)^{1/2}(1+w)^4}\\
&\quad\times (w-\frac{1}{10})(w-\frac{1}{5})(w-\frac{1}{2})(w-10)(w-5)(w-2).
\end{split}
\end{equation}
Hence, $g_{*}^0(R)$ has exactly three positive zeros, denoted by
\begin{gather*}
\tilde{R}_1=\frac{\sqrt{29997}}{99}\approx1.749458791,\quad
\tilde{R}_2=\frac{\sqrt{1872}}{24}\approx1.802775638,\\
\tilde{R}_3=\sqrt{5}\approx2.236067977,
\end{gather*}
which correspond to
$$
\tilde{w}_1=\frac{1}{10},\quad \tilde{w}_2=\frac{1}{5},\quad
\tilde{w}_3=\frac{1}{2}
$$
in $R\in(\sqrt{3}, +\infty)$, respectively. Moreover, we have
\begin{gather*}
\frac{dg_{*}^0(\tilde{R}_1)}{dR}=\frac{25137}{96800}\approx0.2596797521>0,\\
\frac{dg_{*}^0(\tilde{R}_2)}{dR}=-\frac{49}{800}\approx-0.06125<0,\\
\frac{dg_{*}^0(\tilde{R}_3)}{dR}=\frac{19}{800}\approx0.02375>0.
\end{gather*}
According to Lemma \ref{lem1} and Corollary \ref{coro1}, we obtain that
for the sufficiently small $|\varepsilon|$, system \eqref{eq23} has
exactly three limit cycles emerging from the period annulus of
the unperturbed system. Hence, we completes the proof of Theorem \ref{thm2}.

\begin{theorem}\label{thm4}
For the sufficiently small $|\varepsilon|$, system \eqref{eq24} has
just three periodic solutions, denoted  by
$\tilde{l}_1, \tilde{l}_2$ and $\tilde{l}_3$ respectively, such that
$\tilde{l}_1$ shrinks to $\tilde{R}_1,$ $\tilde{l}_2$ shrinks to
$\tilde{R}_2$ and $\tilde{l}_3$ shrinks to $\tilde{R}_3$ as
$\varepsilon$ goes to $0$.
Moreover, $\tilde{l}_1$ and $\tilde{l}_3$ are unstable while
$\tilde{l}_2$ is stable.
\end{theorem}

\subsection*{Acknowledgments} 

This work is supported by the National
Science Foundation of China under contracts No. 11371046 and No.11290141.
The first author would like to thank the Department of
Mathematics at the University of Texas-Pan American for its
hospitality and generous support during her visiting from January
2013 to January 2014.

\begin{thebibliography}{00}

\bibitem{AI} V. I. Arnold, Y. S. Ilyashenko;
\emph{Dynamical systems I:
Ordinary differential equations, Encyclopaedia Math. Sci., Vol. 1},
Springer, Berlin, 1986.

\bibitem{AN1} A. Atabaigi, N. Nyamoradi, H. R. Z. Zangeneh;
\emph{The number of limit cycles of a quintic polynomial system}, 
Comput. Math. Appl., 57 (2009), 677-684.

\bibitem{AN2} A. Atabaigi, N. Nyamoradi, H. R. Z. Zangeneh;
\emph{The number of limit cycles of a quintic polynomial system with a center,
Nonlinear Anal.}, 71 (2009), 3008-3017.


\bibitem{BL1} A. Buic\u{a}, J. Llibre;
\emph{Averaging methods for finding periodic orbits via Brouwer degree}, 
Bull. Sci. Math., 128 (2004), 7-22.

\bibitem{BL2} A. Buic\u{a}, J. Llibre;
Limit cycles of a perturbed cubic polynominal
differential center, Chaos Solitons Fractals, 32 (2007), 1059-1069.

\bibitem{BP} T. R. Blows, L. M. Perko;
Bifurcation of limit cycles from centers and separatrix cycles of planar 
analytic systems, SIAM Rev., 36 (1994), 341-376.

\bibitem{BN} S. Benditkis, D. Novikov;
On the number of zeros of Melnikov fuctions, arXiv: 1007.0672vl, 
[math. DS] 5, July, 2010.

\bibitem{CLLZ} F. D. Chen, C. Li, J. Llibre, Z. H. Zhang;
\emph{A unified proof on the weak Hilbert 16th problem for $n=2$}, 
J. Differential Equations, 221 (2006), 309-342.

\bibitem{CLP} B. Coll, J. Llibre, R. Prohens;
\emph{Limit cycles bifurcating from a perturbed quartic center}, 
Chaos Solitons Fractals, 44 (2011), 317-334.

\bibitem{CGP} B. Coll, A. Gasull, R. Prohens;
\emph{Bifurcation of limit cycles from two families of centers}, 
Dyn. Contin. Discrete Impuls. Syst., Ser. A (Math. Anal.), 12 (2005), 275-287.

\bibitem{CJ} C. Chicone, M. Jacobs;
\emph{Bifurcation of limit cycles from quadratic isochrones}, 
J. Differential Equations, 91 (1991), 268-326.

 \bibitem{E} J. \'Ecalle;
\emph{Introduction aux fonctions analysables et  preuve constructive de la 
conjecture de Dulac}, Actualiti\'ees Math., Hermann, Paris,  1992.

\bibitem{GH}J. Guckenheimer, P. Holmes;
\emph{Nonlinear oscillations, dynamical systems, and bifurcations of vector fields}, 
Appl. Math. Sci., Vol. 42, New York, Springer-Verlag, 1986.

\bibitem{GI} L. Gavrilov, I. D. Iliev;
\emph{Quadratic perturbations of quadratic codimension-four centers}, 
J. Math. Anal. Appl., 357 (2009), 69-76.

 \bibitem{GGI} S. Gautier, L. Gavrilov, I. D. Iliev;
\emph{Perturbations of quadratic center of genus one},
Discrete Contin. Dyn. Syst., 25 (2009), 511-535.

\bibitem{GLV1} H. Giacomini, J. Llibre, M. Viano;
\emph{On the nonexistence, existence and uniqueness of limit cycles},
Nonlinearity, 9(1996), 501-516.

\bibitem{GLV2} H. Giacomini, J. Llibre, M. Viano;
\emph{On the shape of limit cycles that bifurcate from Hamiltonian centers}, 
Nonlinear Anal., 41 (2000), 523-537.

\bibitem{GLV3} H. Giacomini, J. Llibre, M. Viano;
\emph{On the shape of limit cycles that bifurcate from non-Hamiltonian centers}, 
Nonlinear Anal., 43 (2001), 837-859.

\bibitem{GL} J. Gin\'e, J. Llibre;
\emph{Limit cycles of cubic polynomial vector fields via the averaging theory}, 
Nonlinear Anal., 66 (2007), 1707-1721.

 \bibitem{I0} I. D. Iliev;
\emph{Perturbations of quadratic centers},
 Bull. Sci. Math., 122 (1998), 107-161.

 \bibitem{I1} Y. S. Ilyashenko;
\emph{Finiteness theorems for limit cycles},
 Uspekhi Mat. Nauk, 45 (1990), 143-200 (Russian); English transl. Russian Math. 
Survey, 45 (1990), 129-203.

 \bibitem{I2} Y. S. Ilyashenko;
\emph{Centennial history Hilbert's 16th problem},
    Bull. Amer. Math. Soc. (N.S.), 39 (2002), 301-354.

\bibitem{LLZ}C. Li, J. Llibre, Z. Zhang;
\emph{Weak focus, limit cycles and  bifurcations for bounded quadrtic systems},
 J. Differential Equations, 115 (1995), 193-223.

 \bibitem{LL} C. Li, J. Llibre;
\emph{Quadratic perturbations of a  quadratic reversible Lotka-Volterra system}, 
Qual. Theory Dyn.  Syst., 9 (2010), 235-249.

 \bibitem{LPR} J. Llibre, J. S. P\'erez del R\'io, J.  A. Rodr\'iguez;
\emph{Averaging analysis of a perturbed quadratic  center}, 
Nonlinear Anal., 46 (2001), 45-51.

 \bibitem{L1} J. Llibre;
\emph{Averaging theory and limit cycles for
 quadratic systems}, Radovi Matemati\u{c}ki, 11 (2002), 1-14.

\bibitem{SV} J. A. Sanders, F. Verhulst;
\emph{Averaging methods in nonlinear dynamical systems}, 
Appl. Math. Sci., Vol. 59, New York, Springer-Verlag, 1985.

\bibitem{V} F. Verhulst;
\emph{Nonlinear differential equations and dynamical systems}, 
Universitext, Berlin, Springer-Verlag, 1996.

\bibitem{VLG} M. Viano, J. Llibre, H. Giacomini;
\emph{Arbitrary order bifurcations for perturbed Hamiltonian planar 
systems via the reciprocal of an integrating factor}, Nonlinear Anal., 48(2002),
117-136.

\bibitem{XH1} G. Xiang, M. Han;
\emph{Global bifurcation of limit cycles in a family of polynomial systems}, 
J. Math. Anal. Appl., 295(2004), 633-644.

\bibitem{XH2} G. Xiang, M. Han;
\emph{Global bifurcation of limit cycles in a family of multiparameter systems}, 
Int. J. Bifur. Chaos Appl. Sci. Engrg., 14 (2004), 3325-3335.

\end{thebibliography}

\end{document}