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

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

\begin{document}
\title[\hfilneg EJDE-2015/271\hfil Cubic homogeneous polynomial centers]
{Limit cycles bifurcating from the periodic annulus of cubic 
homogeneous polynomial centers}

\author[J. Llibre, B.D. Lopes, J.R. de Moraes \hfil EJDE-2015/271\hfilneg]
{Jaume Llibre, Bruno D. Lopes, Jaime R. de Moraes}  

\address{Jaume Llibre \newline
Departament de Matem\`atiques,
Universitat Aut\`onoma de Barcelona, 08193 Bellaterra, Barcelona,
Catalonia, Spain}
\email{jllibre@mat.uab.cat}

\address{Bruno D. Lopes \newline
 Departamento de Matem\'atica - IBILCE-UNESP,
Rua C. Colombo, 2265, CEP 15054-000 S. J. Rio Preto, S\~ao Paulo,
Brazil}
\email{brunodomicianolopes@gmail.com}

\address{Jaime R. de Moraes \newline
Curso de Matem\'atica - UEMS,
Rodovia Dourados-Itaum Km 12, CEP 79804-970 Dourados, 
Mato Grosso do Sul, Brazil}
\email{jaime@uems.br}


\thanks{Submitted October 8, 2014. Published October 21, 2015.}
\subjclass[2010]{37C10, 34H99}
\keywords{Polinomial vector field; limit cycle; averaging method;
periodic orbit; 
\hfill\break\indent isochronous center; homogeneous cubic centers} 

\begin{abstract}
 We obtain an explicit polynomial whose simple positive real roots
 provide the limit cycles which bifurcate from the periodic orbits
 of any cubic homogeneous polynomial center when it is perturbed
 inside the class of all polynomial differential systems of degree $n$.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section] 
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks


\section{Introduction and statement of main results}

One of the main goals in the qualitative theory of real planar
differential systems is the determination of their limit cycles.
It is well known that perturbing the periodic orbits of a center
often produces limit cycles, see for instance
\cite{ALGM-biurcation, BloPer94, Yi4}. One of the first in studying
these perturbations was Pontrjagin \cite{Po}. These last years
this problem has been studied by many authors see the second part
of the book \cite{Chr} and the hundreds of references quoted there.


Hilbert in 1900 was interested in the maximum number of the limit
cycles that a polynomial differential system of a given degree can
have. This problem is the well-known 16-th Hilbert problem,
which together with the Riemann conjecture are the two problems
of the famous list of 23 problems of Hilbert which  remain open.
See for more details \cite{HI,SS}.


There exist several methods to study the number of limit cycles
that bifurcate from the periodic annulus of a center, such as the
\emph{Poincar\'e return map}, the \emph{Poincar\'e-Melnikov
integrals}, the \emph{Abelian integrals}, the \emph{inverse
integrating factor}, and the \emph{averaging theory}. In the plane
all of them are essentially equivalent.


There are few works trying to study this problem for homogeneous
cubic polynomial differential systems. Our main objetive will be
to solve this problem for the cubic  homogeneous polynomial differential systems.


In \cite{Cima} the authors classified all the cubic homogeneous
polynomial differential systems. In \cite{Li2} the authors proved
that any real planar cubic homogeneous polynomial differential
system having a center can be written as
\begin{equation}
\begin{gathered}
\dot{x}=ax^3+(b-3\alpha\mu)x^2y-axy^2-\alpha y^3 = P(x,y),\\
\dot{y}=\alpha x^3+ax^2y+(b+3\alpha\mu)xy^2-ay^3 = Q(x,y),
\end{gathered}\label{sis1}
\end{equation}
with $\alpha\in\{-1,1\}$, $a,b,\mu\in\mathbb{R}$ and $\mu>-1/3$, after doing
an affine change of variables and a rescaling of the time.

It is known that the maximum number of limit cycles which bifurcate
from the periodic orbits of a cubic homogeneous center \eqref{sis1}
using perturbations of first order inside the class of all polynomial
differential systems of degree $n$ is $[(n-1)/2]$, see for details statement (c)
of Theorem A of \cite{Li1}. Here $[x]$ denotes the integer part function of $x$.

The objetive of this work is \emph{to provide an explicit
polynomial whose real positive simple zeros gives the exact number
of limit cycles which bifurcate, at first order in the perturbation
parameter, from the periodic orbits of any cubic homogeneous
center \eqref{sis1}}.
More precisely consider the system.
\begin{equation}
\begin{gathered}
\dot{x}=ax^3+(b-3\alpha\mu)x^2y-axy^2-\alpha y^3 + \varepsilon p(x,y),\\
\dot{y}=\alpha x^3+ax^2y+(b+3\alpha\mu)xy^2-ay^3 + \varepsilon q(x,y),
\end{gathered} \label{sis2}
\end{equation}
where
\begin{equation}\label{eq1}
p(x,y)=\sum_{i=0}^{n}p_i(x,y), \quad q(x,y)=\sum_{i=0}^{n}q_i(x,y),
\end{equation}
$p_i,q_i$ are homogeneous polynomials of degree $i$, and $\varepsilon$ is
a small parameter.

 Define the following functions
\begin{gather*}
\begin{aligned}
f_1({\theta})&= -a \sin ^4\theta +a \cos ^4\theta + (b-3 \alpha
\mu +\alpha )\sin \theta \cos ^3\theta \\
&\quad +(b+3\alpha  \mu -\alpha ) \sin ^3 \theta  \cos \theta,
\end{aligned}\\
g_1({\theta})=\alpha  \left(6 \mu  \sin ^2\theta  \cos ^2\theta +\sin
   ^4\theta +\cos ^4\theta \right),\\
k(\theta)=\exp\Big(\int_0^\theta\frac{f_1(s)}{g_1(s)}ds \Big),\\
B_i(\theta)=Q(\cos\theta,\sin\theta) p_i(\cos\theta,\sin\theta)
-P (\cos\theta,\sin\theta)q_i(\cos\theta,\sin\theta).
\end{gather*}

In sequel we state our main result where  the function
$M(\theta)$ is defined in \eqref{eqM}, we do not provide it here
due to its length.

\begin{theorem}\label{teo1}
For $|\varepsilon|>0$ sufficiently small and for every positive simple zero
$r_{0}^{*}$ of the polynomial
\begin{equation*}\label{eqS01}
r_0\mathcal{F}(r_0)= \frac{1}{2\pi}\displaystyle\sum_{k=0}^{[
\frac{n-1}{2}]}r_0^{2k}\int_0^{2\pi}A_{2k+1}(\theta)d\theta,
\end{equation*}
where
$$
A_i(\theta)=\frac{B_i(\theta) k(\theta)^{i-2}}{g_1(\theta)^2
M(\theta)},
$$
for $i=1,2,\dots,[\frac{n-1}{2}]$, the perturbed systems
\eqref{sis2} has a limit cycle bifurcating from the periodic orbit
$r(\theta, r_0^{*})=k(\theta)r_0^{*}$  of the period annulus of the center
\eqref{sis1} using the averaging theory of first order.
In particular the perturbed systems \eqref{sis2} has
at most $[\frac{n-1}{2}]$ limit cycles.
\end{theorem}

Theorem \ref{teo1} is proved in Section \ref{sec3}. In Section \ref{sec4} we
provide an example that illustrates Theorem \ref{teo1} with $n=5$. We obtain
two limit cycles.


\section{Preliminaries}

In this section we give some known results that we shall need for
proving  Theorem \ref{teo1}.
Consider a system in the form
\begin{equation}\label{eq2}
\dot{\mathrm{x}}=F_0(t,\mathrm{x})+\varepsilon F_1(t,\mathrm{x})+\varepsilon^2
F_2(t,\mathrm{x},\varepsilon),
\end{equation}
where $\varepsilon\neq0$ is sufficiently small and the functions
$F_0,F_1:\mathbb{R}\times\Omega\to \mathbb{R}^n$ and $F_2:\mathbb{R}\times
\Omega\times (-\varepsilon_0,\varepsilon_0)\to \mathbb{R}^n$
are $\mathcal{C}^2$ functions, $T-$periodic in the first variable
and $\Omega$ is an open subset of $\mathbb{R}^n$. We assume that the
unperturbed system
\begin{equation}\label{eq3}
\dot{\mathrm{x}}=F_0(t,\mathrm{x})
\end{equation}
has a submanifold of periodic solutions of dimension $n$.

Consider $x(t,z,\varepsilon)$ the solution of system \eqref{eq3} such
that $\mathrm{x}(0,\mathrm{z},\varepsilon)=z$. The linearization of the unperturbed
system along a periodic solution $\mathrm{x}(t,\mathrm{z},0)$ is given by
\begin{equation}\label{eq4}
\dot{\mathrm{y}}=D_{\mathrm{x}}F_0(t,\mathrm{x}(t,\mathrm{z},0))\mathrm{y}.
\end{equation}

In sequel we denote by $M_{\mathrm{z}}(t)$ the fundamental matrix of
the linearized system \eqref{eq4} such that $M_{\mathrm{z}}(0)$ is the identity.

We suppose that there exists an open set $V$ with
$\mathrm{Cl}(V)\subset\Omega$ such that for each
$z\in\mathrm{Cl}(V)$, $\mathrm{x}(t,\mathrm{z},0)$ is $T-$periodic, where
$\mathrm{x}(t,z,0)$ denotes the solution of the unperturbed system
\eqref{eq3}. Here $\mathrm{Cl}(V)$ denotes the closure of $V$.  We
have that the set $\mathrm{Cl}(V)$ is \emph{isochronous} for system
\eqref{eq3}, i.e. it is formed only by periodic orbits with period
$T$.

The next result is the averaging theorem for studying the
bifurcation of $T$-periodic solutions of system \eqref{eq2} from the
periodic solutions $\mathrm{x}(t,\mathrm{z},0)$ contained in
$\mathrm{Cl}(V)$ of system \eqref{eq3} when $|\varepsilon|>0$ is
sufficiently small. See \cite{Bui1} for a proof. For more details on
the averaging theory see \cite{BuiLli} and the book \cite{SanVer}.

\begin{theorem}[Perturbations of an isochronous set]\label{teo2}
We assume that there exists an open and bounded set $V$ with
$\mathrm{Cl}(V)\subset\Omega$ such that for each
$\mathrm{z}\in \mathrm{Cl}(V)$, the solution $\mathrm{x}(r,\mathrm{z},0)$
is $T-$periodic. Consider the function
$\mathcal{F}:\mathrm{Cl}(V)\to \mathbb{R}^n$
\begin{equation}\label{eq5}
\mathcal{F}(\mathrm{z})=\frac{1}{T}\int_0^T M_{\mathrm{z}}^{-1}(t)
F_1(t,\mathrm{x}(t,\mathrm{z},0))dt.
\end{equation}
Then the following statements hold.
\begin{itemize}
\item[(i)] If there exists $\mathbf{a}\in V$ with
$\mathcal{F}(\mathbf{a})=0$ and
$\det((\partial \mathcal{F}/\partial \mathrm{z})(\mathbf{a}))\neq0$
then there exists a $T-$periodic solution $\mathrm{x}(t,\varepsilon)$ of
system \eqref{eq2} such that $\mathrm{x}(0,\varepsilon)\to \mathbf{a}$
when $\varepsilon\to 0$.



\item[(ii)] The kind of the stability of the periodic solution
$\mathrm{x}(t,\varepsilon)$ is given by the eigenvalues of the Jacobian
matrix $((\partial \mathcal{F}/\partial \mathrm{z})(\mathbf{a}))$.
\end{itemize}
\end{theorem}

\section{Proof of Theorem \ref{teo1}}\label{sec3}
The next result follows easily.

\begin{lemma}\label{lema1}
Let $P_k(x,y)$ and $P_3(x,y)$ be homogeneous polynomials of degree
$k$ and\, $3$ respectively, where $(x,y)\in\mathbb{R}^2$. Thus in polar
coordinates $x=r\cos\theta$ and $y=r\sin\theta$ we have
\begin{align*}
&P_k(r\cos\theta,r\sin\theta)P_3 (r\cos\theta,r\sin\theta)\\
&=(-1)^{k+1}P_k(r\cos(\theta+\pi),r\sin(\theta+\pi)
 P_3 (r\cos(\theta+\pi),r\sin(\theta+\pi)).
\end{align*}
\end{lemma}

Now we pass system \eqref{sis2} to polar coordinates taking
$x=r \cos\theta$, $y=r \sin\theta$ and we obtain
\begin{gather*}
\dot{r} = r^3 f_1(\theta)+ \varepsilon(\cos\theta \ p(r\cos\theta,r
\sin\theta)+\sin\theta \ q(r\cos\theta,r \sin\theta)),\\
\dot{\theta} = r^2g_1(\theta)+  \varepsilon \frac{1}{r}(\cos\theta \
q(r\cos\theta,r \sin\theta)-\sin\theta \ p(r\cos\theta,r \sin\theta)).
\end{gather*}
Note that $g_1({\theta}) \neq 0$ for all $ \theta \in [0,2\pi]$.
Thus we take the quotient $\dot{r}/\dot{\theta}$  and  we get the
differential equation
\begin{equation}\label{eq100}
\frac{dr}{d\theta}= F_0(r,\theta)+\varepsilon F_1(r,\theta)+
\mathcal{O}(\varepsilon^2),
\end{equation}
in the standard form for applying the
averaging theory of first order, where
\begin{gather*}
F_0(r,\theta) = \frac{f_1(\theta)}{g_1(\theta)}r, \\
\begin{aligned}
F_1(r,\theta) 
&= \frac{1}{r^5 g_1(\theta)^2}(Q(r\cos\theta,r \sin\theta)
p(r\cos\theta,r \sin\theta) \\
&\quad - P(r\cos\theta,r \sin\theta)q(r\cos\theta,r \sin\theta)).
\end{aligned}
\end{gather*}
Note that the differential equation \eqref{eq100} satisfies the
assumptions of Theorem \ref{teo2}. Consider $r(\theta,r_0)$ the
periodic solution of the differential equation
$ \dot{r}=  r f_1(\theta)/g_1(\theta)$ such that $r(0,r_0)= r_0$. 
Solving this differential equation we obtain
$$
r(\theta,r_0)= k(\theta )r_0=r_0 e^{k_1(\theta )} k_2(\theta ),
$$
where
\begin{gather*}
\begin{aligned}
k_1(\theta)&= -\frac{a \alpha}{2R}  \bigg(\frac{(-3 \mu +R-1)
\tan ^{-1}\big(\frac{\tan \theta }{\sqrt{3 \mu -R}}\big)}
{\sqrt{3 \mu -R}} 
 +\frac{(3 \mu +R+1) \tan ^{-1}\big(\frac{\tan \theta }
{\sqrt{3 \mu +R}}\big)}{\sqrt{3 \mu +R}}\bigg),
\end{aligned}\\
\begin{aligned}
k_2(\theta)&= \sqrt{\sec ^2\theta } (3 \mu -R)^{\frac{R-\alpha  b}{4 R}}
(3 \mu +R)^{\frac{\alpha  b+R}{4 R}} \left(\tan ^2\theta+
3 \mu -\alpha  R\right)^{\frac{1}{4}  \left(\frac{b}{R}-1\right)}\\
&\quad\times  \left(\tan ^2\theta+3 \mu +\alpha  R\right)^{-\frac{b+R}{4 R}},
\end{aligned}\\
R=\sqrt{9\mu^2-1}.
\end{gather*}

Solving the variational equation \eqref{eq4} for our differential
equation \eqref{eq100} we get that the fundamental matrix of \eqref{eq4} is
\begin{equation} \label{eqM}
\begin{aligned}
M(\theta)
&=\bigg(-2 a \alpha  (-3 \mu +R-1) \sqrt{3 \mu +R}
\tan ^{-1}\Big(\frac{\tan \theta}{\sqrt{3 \mu -R}}\Big) \\
&\quad  -2 a \alpha  \sqrt{3 \mu -R} (3 \mu +R+1) \tan
   ^{-1}\Big(\frac{\tan \theta}{\sqrt{3 \mu +R}}\Big)+\alpha
\sqrt{3 \mu -R} \\
&\quad\times \sqrt{3 \mu +R} (b-\alpha  R) \log \left(\tan ^2\theta
+3 \mu -R\right)-\alpha  \sqrt{3 \mu -R} \\
&\quad\times \sqrt{3 \mu +R} (b+\alpha  R) \log \left(\tan ^2\theta
+3 \mu +R\right)+\sqrt{3 \mu -R} \sqrt{3 \mu +R}\\
&\quad\times  (R-\alpha  b) \log (3 \mu -R)+\sqrt{3
   \mu -R} \sqrt{3 \mu +R} (\alpha  b+R) \log (3 \mu +R) \\
&\quad +2 R \sqrt{3 \mu -R} \sqrt{3 \mu +R} \log \left(\sec ^2\theta\right) \\
&\quad +4 R \sqrt{3 \mu -R} \sqrt{3 \mu +R}\bigg)
\big/\Big(4 R \sqrt{3 \mu -R} \sqrt{3 \mu +R}\Big) .
\end{aligned}
\end{equation}
Note that $M(\theta)$ does not depend on $r_0$. Using
the polynomials $p$ and $q$ given in \eqref{eq1} and system
\eqref{sis1} we have that the integrant of the integral
\eqref{eq5} for our differential equation is
\begin{equation}\label{eq101}
M^{-1}(\theta)F_1(\theta ,r(\theta,r_0))=
\frac{F_1(\theta ,r(\theta,r_0))}{M(\theta)}=
\frac{h(r,\theta)}{r^5 g_1(\theta)^2 M(\theta)},
\end{equation}
where
$$
h(r,\theta)= Q(r\cos\theta,r\sin\theta) p(r\cos\theta,r\sin\theta)
-P(r\cos\theta,r\sin\theta) q(r\cos\theta,r\sin\theta).
$$
Since $p$ and $q$ are sum of homogeneous polynomials
(see eq. \eqref{eq1}), we can rewrite the equality \eqref{eq101}
as follows
\begin{align*}
M^{-1}(\theta)F_1(\theta ,r(\theta,r_0))
&= \sum_{i=0}^{n}\frac{B_i(\theta)}{g_1(\theta)^2M(\theta)}
r(\theta,r_0)^{i-2} \\
&= \sum_{i=0}^{n} r_0^{i-2} \frac{B_i(\theta) (e^{k_1(\theta)}
k_2(\theta))^{i-2}}{g_1(\theta)^2M(\theta)}\\
&= \sum_{i=0}^{n} r_0^{i-2}A_i(\theta).
\end{align*}
Computing the integral \eqref{eq5} we have
$$
\mathcal{F}(r_0)= \frac{1}{2\pi}\int_0^{2\pi} M^{-1}
(\theta)F_1(\theta ,r(\theta,r_0) d\theta =   \frac{1}{2\pi}
\displaystyle\sum_{i=0}^{n}r_0^{i-2}\int_0^{2\pi}A_i(\theta)d\theta.
$$

If $i$ is even then Lemma \ref{lema1} implies
$B_i(\theta )= -B_i(\theta+\pi)$. Thus since
$k_1(\theta)=k_1(\theta+\pi)$, $k_2(\theta)=k_2(\theta+\pi)$ and
$M(\theta)=M(\theta+\pi)$ we easily obtain
\begin{align*}
\int_\pi^{3\pi/2} A_i(\theta)d\theta 
&= \int_\pi^{3\pi/2} \frac{B_i(\theta) (e^{k_{1}(\theta) }
k_2(\theta))^{i-2}}{g_1(\theta)^2 M(\theta)}d\theta\\
&= \int_0^{\pi/2}\frac{B_i(\theta+\pi) (e^{k_1(\theta+\pi)}
k_2(\theta+\pi))^{i-2}}{g_1(\theta+\pi))^2 M(\theta+\pi))}d\theta\\
&= \int_0^{\pi/2}-\frac{B_i(\theta) (e^{k_{1}(\theta)}
k_2(\theta))^{i-2}}{g_1(\theta)^2 M(\theta)}d\theta\\
&=  - \int_0^{\pi/2}A_i(\theta)d\theta,
\end{align*}
\begin{align*}
\int_{3\pi/2}^{2\pi} A_i(\theta)d\theta 
&= 
\int_{3\pi/2}^{2\pi} \frac{B_i(\theta) (e^{k_{1}(\theta)}
k_2(\theta))^{i-2}}{g_1(\theta)^2 M(\theta)}d\theta\\
&= \int_{\pi/2}^{\pi} \frac{B_i(\theta+\pi) (e^{k_1(\theta+\pi)}
k_2(\theta+\pi))^{i-2}}{g_1(\theta+\pi))^2 M(\theta+\pi))}d\theta\\
&= \int_{\pi/2}^{\pi} - \frac{B_i(\theta) (e^{k_{1}(\theta) }
k_2(\theta))^{i-2}}{g_1(\theta)^2 M(\theta)}d\theta\\
&=  - \int_{\pi/2}^{\pi} A_i(\theta)d\theta.
\end{align*}
Therefore for $i$ even we have
$$ 
\int_{0}^{2\pi} A_i(\theta)d\theta = 0.
$$
Analogously if $i$ is odd then Lemma \ref{lema1} implies
$B_i(\theta )= B_i(\theta+\pi)$. Thus we easily can check that
$$
\int_0^{\pi/2}A_i(\theta)d\theta=\int_\pi^{3\pi/2}
A_i(\theta)d\theta \quad \text{and} \quad
\int_{\pi/2}^{\pi} A_i(\theta)d\theta= \int_{3\pi/2}^{2\pi} A_i(\theta)d\theta.
$$
So we have
$$ 
\int_{0}^{2\pi} A_i(\theta)d\theta =  2 \int_{0}^{\pi} A_i(\theta)
d\theta \neq 0,
$$
because $A_i$ is a $\pi-$periodic even function.
So the function $\mathcal{F}$ can be written in the following way
\begin{equation}\label{eqS0}
\mathcal{F}(r_0)= \frac{1}{2\pi}\sum_{k=0}^{[
\frac{n-1}{2}]}r_0^{2k-1}\int_0^{2\pi}A_{2k+1}(\theta)d\theta.
\end{equation}

Note that the coefficients $A_{2k+1}(\theta )$ in \eqref{eqS0}
are linearly independent because the polynomials $p_i$ and
$q_i$ are linearly independents. Thus the averaged function
$\mathcal{F}$ has at most $[(n-1)/2]$ simple zeros which correspond
to the limit cycles of system \eqref{sis2} and Theorem \ref{teo1}
is proved.


\section{Example}\label{sec4}

In this section we present an example that illustrates Theorem \ref{teo1}.
Consider the cubic polynomial homogeneous center
\[\dot{x}=-y^3,\quad \dot{y}=x^3,\]
and its perturbation
\begin{equation}\label{eq6}
\dot{x}=-y^3+\varepsilon(a_1 x+a_2 x^3+a_3 x^5),\quad \dot{y}=x^3.
\end{equation}
Passing system \eqref{eq6} to the polar coordinates we get
\begin{gather*}
\dot{r} =  r^3 (\sin \theta \cos ^3\theta
   -\sin ^3\theta \cos \theta)+ \varepsilon  r\cos ^2\theta 
(a_1+a_2 r^2 \cos ^2\theta+a_3
r^4
   \cos ^4\theta),\\
 \dot{\theta} = r^2 (\sin ^4\theta+\cos ^4\theta)-\varepsilon  \sin
   \theta \cos \theta (a_1+a_2 r^2 \cos ^2\theta+a_3 r^4
   \cos ^4\theta) .
\end{gather*}
Taking the quotient $\dot{r}/\dot{\theta}$ we obtain the following
system in the standard form of Theorem \ref{teo2} for applying the
averaging theory
\begin{equation}\label{abc}
\frac{dr}{d\theta}=F_0(r,\theta)+\varepsilon
F_1(r,\theta)+\mathcal{O}(\varepsilon^2),
\end{equation}
where
\begin{gather*}
F_0(r,\theta) = \frac{r \sin \theta \cos ^3\theta-r \sin ^3\theta \cos
\theta}{\sin ^4\theta+\cos ^4\theta},\\
F_1(r,\theta) = \frac{\cos ^4\theta
(a_1+a_2 r^2 \cos ^2\theta+a_3 r^4 \cos
 ^4\theta)}{r (\sin ^4\theta+\cos ^4\theta )^2}.
\end{gather*}
Thus for system  \eqref{abc} we have
\begin{gather*}
k(\theta) = \frac{\sqrt{2}}{\sqrt[4]{\cos (4 \theta )+3}},\\
M(\theta) = \frac{1}{4} (-\log (\cos (4 \theta )+3)+4+\log 4),
\end{gather*}
and the integrant of the integral \eqref{eq5} of system  \eqref{abc}
is
\[
\frac{A(\theta)+B(\theta)r_0^2+C(\theta)r_0^4}{r_0},
\]
where
\begin{gather*}
A(\theta) = a_1\frac{32 \sqrt{2} \cos ^4\theta}{(\cos (4 \theta )+3)^{7/4}
   (-\log (\cos (4 \theta )+3)+4+\log 4)},\\
B(\theta) = a_2\frac{64 \sqrt{2} \cos ^6\theta}{(\cos (4 \theta )+3)^{9/4}
   (-\log (\cos (4 \theta )+3)+4+\log 4)},\\
C(\theta) = a_3\frac{128 \sqrt{2} \cos ^8\theta}{(\cos (4 \theta )+3)^{11/4}
   (-\log (\cos (4 \theta )+3)+4+\log 4)}.
\end{gather*}
Computing numerically the integral \eqref{eq5} for system
\eqref{abc}  we obtain
\[
\mathcal{F}(r_0)=\frac{ 3.72731\dots a_1 + 3.34745\dots a_2 r_0^2+ 3.10284\dots
 a_3 r_0^4 }{r_0}.
\]
Taking
\[
a_1=\frac{-1}{3.72731\dots}, \quad 
a_2=\frac{2}{3.34745\dots}\quad\text{and}
a_3=\frac{-0.1}{3.10284\dots},
\]
it is easy to check that the function $\mathcal{F}$ has two positive
simple zeros given by
\[
r_0^{*}=0.716357\dots\quad \text{and}\quad r_0^{**}=4.41439\dots
\]
which correspond to two limit cycles of the perturbed system
\eqref{eq6} with $\varepsilon\ne 0$ sufficiently small.


\subsection*{Acknowledgements}
The first author is partially supported by a MINECO/FEDER grant
number MTM2009-03437, by an AGAUR grant number 2014SGR-568, by an
ICREA Academia, two FP7+PEOPLE+2012+IRSES numbers 316338 and 318999,
and FEDER-UNAB10-4E-378. The second author is supported
by CAPES/GDU - 7500/13-0.


\begin{thebibliography}{10}

\bibitem{ALGM-biurcation}
A.~A. Andronov, E.~A. Leontovich, I.~I. Gordon, A.~G. Ma{\u\i}er;
  \emph{Theory of bifurcations of dynamic systems on a plane}, Halsted Press [A
  division of John Wiley \& Sons], New York-Toronto, Ont., 1973, Translated
  from the Russian.

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

\bibitem{Bui1} A. Buic{\u{a}}, J. P. Fran\c{c}oise, J. Llibre;
 \emph{Periodic solutions of
  nonlinear periodic differential systems with a small parameter}, Commun. Pure
  Appl. Anal., \textbf{6} (2007), 103--111.

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

\bibitem{Chr} C.~Christopher, C.~Li;
\emph{Limit cycles in differential equations},   Birkhauser, Boston, 1999.

\bibitem{Cima} A.~Cima, J.~Llibre;
\emph{Algebraic and topological classification of the
  homogeneous cubic vector fields in the plane}, J. Math. Anal. Appl.,
  \textbf{147} (1990), 420--448.

\bibitem{HI} C.~Hilbert;
\emph{Mathematische probleme}, Lecture, Second Internat. Congr.
  Math. (Paris, 1900), \emph{Nachr. Ges. Wiss. G''ttingen Math. Phys. KL.}
  (1900), 253--297; English transl., \emph{Bull. Amer. Math. Soc.} \textbf{8}
  (1902), 437--479; \emph{Bull. (New Series) Amer. Math. Soc.} \textbf{37}
  (2000), 407--436.

\bibitem{Li2} C.~Li, J.~Llibre;
\emph{Cubic homogeneous polynomial centers}, Publ. Mat.,
  \textbf{58} (2014), 297--308.

\bibitem{Li1} W.~Li, J.~Llibre, J.~Yang, Z.~Zhang;
\emph{Limit cycles bifurcating from
  the period annulus of quasi--homogeneous centers}, J. Dynam. Differential
  Equations, \textbf{21} (2009), 133--152.

\bibitem{Po} L.~S. Pontrjagin;
\emph{\"Uber autoschwingungssysteme, die den
  hamiltonshen nahe liegen}, Physikalische Zeitschrift der Sowjetunion,
  \textbf{6} (1934), 25--28.

\bibitem{SanVer} J. A. Sanders, F.~Verhulst;
 \emph{Averaging methods in nonlinear dynamical   systems}, 
Applied Mathematical Sciences, vol.~59, Springer, New York, 2007.

\bibitem{Yi4} Y.~Shao, K.~Wu;
\emph{Bifurcation of limit cycles for cubic reversible
  systems}, Electron. J. Differential Equations, \textbf{2014} (2014), No. 96,
  10 pp.

\bibitem{SS} S.~Smale;
\emph{Mathematical problems for next century}, Math. Intelligencer,
\textbf{20} (1998), 7--15.

\end{thebibliography}


\end{document}