\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{cite,amssymb}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 108, pp. 1--17.\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 - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/108\hfil  
Limit cycles for cubic polynomial differential systems]
{Limit cycles for piecewise smooth perturbations of
a cubic polynomial differential center}

\author[S. Li, T. Huang \hfil EJDE-2015/108\hfilneg]
{Shimin Li, Tiren Huang}

\address{Shimin Li (corresponding author)\newline
School of Mathematics and Statistics,
Guangdong University of Finance and Economics,
Guangzhou 510320,  China}
\email{lism1983@126.com}

\address{Tiren Huang \newline
Department of Mathematics, Zhejiang Sci-Tech University,
Hangzhou 310018, China}
\email{htiren@zstu.edu.cn}

\thanks{Submitted October 24, 2014. Published April 21, 2015.}
\makeatletter
\@namedef{subjclassname@2010}{\textup{2010} Mathematics Subject Classification}
\makeatother
\subjclass[2010]{34A36, 34C07, 37G15}
\keywords{Limit cycle; piecewise smooth system; averaging method}

\begin{abstract}
 In this article, we study the planar cubic polynomial differential system
 \begin{gather*}
 \dot{x}=-yR(x,y)\\
 \dot{y}=xR(x,y)
 \end{gather*}
 where $R(x,y)=0$ is a conic and $R(0,0)\neq 0$.
 We find a bound for the number of limit cycles which bifurcate
 from the period annulus of the center, under piecewise
 smooth cubic polynomial perturbations.
 Our results show that the piecewise smooth cubic system can have at
 least 1 more limit cycle than the smooth one.
\end{abstract}

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


\section{Introduction and statement of main results}

In the qualitative theory of real planar differential system, one of
the important problems is the determination and distribution of
limit cycles, such as the famous Hilbert's 16th problem and its weak
form \cite{Li}. Limit cycles can arise by bifurcation in several different
ways. One of the main methods is perturbing a system which has a
center via Poincar\'e bifurcation, in such a way that limit cycles
bifurcate in the perturbed system from the period annulus of the
center for the unperturbed system. This has been studied intensively
perturbing the limit cycles of the center for the Hamiltonian system
(e.g. \cite{GLV}) and non-Hamiltonian integrable system (e.g.\cite{LV}).

In this article, we consider the  non-Hamiltonian integrable system
\begin{equation} \label{e1}
\begin{gathered}
\dot{x}=-yR(x,y),\\
\dot{y}=xR(x,y),
\end{gathered}
\end{equation}
where $R(0,0)\neq 0$ and the dot denotes derivatives with respect to the
variable $t$.
Note that system \eqref{e1} has a first integral $H(x,y)=x^2+y^2$,
thus the origin is a center. The problem of studying the number of limit cycles which
bifurcate from the period annulus surrounding the origin of system
\eqref{e1} under smooth perturbation has been considered intensively, see for
instance \cite{LZL} and the references quoted therein.

Motivated by the non-smooth phenomena in the real world, piecewise
smooth systems were widely investigated in the past few decades
\cite{BD,F,K}. According to the extent of discontinuity, piecewise
smooth systems can be distinguished by the following three classes:
Piecewise smooth continuous systems, Filippov systems and Impacting systems,
see \cite{BD}.

There are several papers \cite{CGP,GT,LHR,LH,YHH} concerning the number of
limit cycles bifurcate from piecewise smooth continuous systems.
In  \cite{CGP,GT}, applying  Lyapunov constants, the
authors study the center-focus problem as well as the number of
limit cycles which bifurcate from a weak
focus for piecewise smooth continuous systems. Bifurcation of limit cycles by
perturbing piecewise smooth Hamiltonian systems has been studied in
\cite{LHR,LH,YHH}. Generally speaking, piecewise smooth differential
system can have more limit cycles than the smooth one.

For piecewise smooth continuous linear systems, the authors in \cite{HZ} showed
that piecewise linear differential system can have 2
limit cycles surrounding the origin. Later on, in \cite{HY}, the
authors proved that piecewise smooth linear systems can have 3
limit cycles surrounding a unique equilibrium. In
\cite{BPT}, the authors consider higher order piecewise smooth
perturbation of a linear center. For piecewise smooth quadratic
systems, the authors in \cite{CD} showed that there are
piecewise quadratic system with 9 small amplitude limit cycles. In
\cite{LM}, the authors showed that there are at least 5
limit cycles can bifurcate from quadratic isochronous centers under
piecewise smooth quadratic perturbations. Recently, the author
in \cite{X} showed that the piecewise smooth quadratic
isochronous systems can have at least 6 limit cycles. To the best of
our knowledge, there are only a few papers considering the number of limit
cycles for piecewise smooth cubic systems.

The objective of this article is to study the number of limit cycles
that bifurcate from the period annulus surrounding the origin of
system \eqref{e1} under piecewise smooth cubic polynomial perturbations.
More precisely, for $|\varepsilon|>0$ sufficiently small, we want to a
bound for the number of limit cycles for the  piecewise smooth system
\begin{equation} \label{e2}
\begin{pmatrix}\dot{x}\\ \dot{y}\end{pmatrix}
=\begin{cases}
\begin{pmatrix}-yR(x,y)+\varepsilon p_1(x,y)\\
xR(x,y)+\varepsilon q_1(x,y)\end{pmatrix}, &  y> 0,
\\[4pt]
\begin{pmatrix}-yR(x,y)+\varepsilon p_2(x,y)\\
xR(x,y)+\varepsilon q_2(x,y)\end{pmatrix}, &  y< 0,
\end{cases}
\end{equation}
where
\begin{equation} \label{e3}
\begin{gathered}
p_1(x,y)=a_1x+a_2y+a_3x^2+a_4xy+a_5y^2+a_6x^{3}
+a_{7}x^2y+a_{8}xy^2+a_{9}y^{3},\\
q_1(x,y)=b_1x+b_2y+b_3x^2+b_4xy+b_5y^2+b_6x^{3}
+b_{7}x^2y+b_{8}xy^2+b_{9}y^{3},\\
p_2(x,y)=c_1x+c_2y+c_3x^2+c_4xy+c_5y^2+c_6x^{3}+c_{7}x^2y
+c_{8}xy^2+c_{9}y^{3},\\
q_2(x,y)=d_1x+d_2y+d_3x^2+d_4xy+d_5y^2+d_6x^{3}
+d_{7}x^2y+d_{8}xy^2+d_{9}y^{3},
\end{gathered}
\end{equation}
and $R(x,y)=0$ is one of the following six conics:
\begin{itemize}
\item[(E)] Ellipse $R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1=0$
 with $a>0$, $b>0$.

\item[(CE)] Complex ellipse
$R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}+1=0$ with $a>0$, $b>0$.

\item[(H)] Hyperbola
$R(x,y)=\frac{x^2}{a^2}-\frac{y^2}{b^2}-1=0$ with $a>0$, %b>0$.

\item[(CL)] Two complex straight lines intersecting in a real point (the
origin) $R(x,y)=a^2x^2+y^2=0$ with $a\neq 0$.
In this case, we allow that $R(0,0)=0$.

\item[(RPL)] Two real parallel straight lines $R(x,y)=x^2-a^2=0$ with
$a>0$.

\item[(CPL)] Two complex parallel straight lines $R(x,y)=x^2+a^2=0$
with $a>0$.
\end{itemize}

Let $H(3)$ denote the maximum number of limit cycles bifurcating
from the period annulus surrounding the origin of piecewise smooth
system \eqref{e2}. From the first order averaging method,  we have the
following theorem.

\begin{theorem}\label{thm1}
Consider system \eqref{e2} with $|\varepsilon|>0$ sufficiently small.
\begin{itemize}
\item[(I)] If $R(x,y)=a^2x^2+y^2$, then $H(3)=2$.

\item[(II)] In the other cases above, $5\leqslant H(3)\leqslant 15$.
\end{itemize}
\end{theorem}

Note that if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then
$p_1(x,y)\equiv p_2(x,y),q_1(x,y)\equiv q_2(x,y)$, and
system \eqref{e2} is a smooth system. Let $\bar{H}(3)$ denote the maximum
number of limit cycles bifurcating from the period annulus
surrounding the origin of smooth system \eqref{e2}. From the first order
averaging method, then we have the following theorem.

\begin{theorem}\label{thm2}
Consider  \eqref{e2} with $|\varepsilon|>0$ sufficiently small,
suppose that $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$.
\begin{itemize}
\item[(i)] If $R(x,y)=a^2x^2+y^2$, then $\bar{H}(3)=1$.

\item[(ii)] In the other cases above, $\bar{H}(3)=3$.
\end{itemize}
\end{theorem}

\begin{remark} \label{rmk3} \rm
From the above two theorems, we know that the piecewise smooth system \eqref{e2}
have at least 1 (resp. 2) more limit cycles than the smooth one for the
case $R(x,y)=a^2x^2+y^2$ (resp. other cases above).

In \cite{LV}, the authors considered the case
for $a_i=c_i$, $b_i=d_i$, $i=1,2,6,7,8,9$ and $a_i=b_i=c_i=d_i=0$, $i=3,4,5$.
\end{remark}

The organizing of this paper is as follows.
In section 2, we introduce the first order averaging method for
discontinuous system derived from \cite{LNT}.
In section 3, we compute the averaged function according the different
kind of conics.
In section 4, we prove Theorems \ref{thm1} and  \ref{thm2}.
In section 5, we give a conclusion.

\section{Averaging method for discontinuous system}

In this section we summarize the theorems of first order averaging
method for discontinuous differential system as obtained in
\cite{LNT}. The original theorem is given for a system of
differential equations, but since we will use it only for one
differential equation, we state them in this case. For a general
introduction to averaging method, see the book \cite{SV}.

We consider the discontinuous differential equation
\begin{equation} \label{e4}
\frac{\mathrm{d}r}{\mathrm{d}\theta}=\varepsilon
F(\theta,r)+\varepsilon^2R(\theta,r,\varepsilon),
\end{equation}
with
\begin{equation} \label{e5}
\begin{gathered}
F(\theta,r)=F_1(\theta,r)+\operatorname{sign}(h(\theta,r))F_2(\theta,r),\\
R(\theta,r,\varepsilon)=R_1(\theta,r,\varepsilon)
+\operatorname{sign}(h(\theta,r))R_2(\theta,r,\varepsilon),
\end{gathered}
\end{equation}
where $F_1, F_2: \mathbb{R}\times D\to \mathbb{R},
R_1, R_2: \mathbb{R}\times D\times (-\varepsilon_{0},
+\varepsilon_{0})\to \mathbb{R}$ and $h:\mathbb{R}\times
D\to \mathbb{R}$ are continuous functions, T-periodic in the
variable in $\theta$ and $D$ is an open subset of $\mathbb{R}$.
The sign function is defined as
\[
\operatorname{sign}(u)=\begin{cases}
1 & \text{if } u>0,\\
0 & \text{if } u=0,\\
-1 & \text{if } u<0.
\end{cases} 
\]
We also assume that $h$ is a $C^{1}$ function having $0$ as a
regular value. Denote by $\mathcal{M}=h^{-1}(0)$, by
$\Sigma=\{0\}\times D \nsubseteqq \mathcal{M}$, by
$\Sigma_{0}=\Sigma\setminus\mathcal {M}\neq \emptyset $, and its elements
by $z\equiv (0,z)\notin \mathcal{M}$.
Define the averaged function $f:D\to \mathbb{R}$ as
\begin{equation} \label{e6}
f(r)=\int_{0}^{T}F(\theta,r)\mathrm{d}\theta.
\end{equation}

\begin{theorem}[\cite{LNT}]\label{thm4}
Assume the following three conditions ate satisfied.
\begin{itemize}
\item[(i)]  $F_1, F_2, R_1, R_2$ and $h$ are locally
$L$-Lipschitz with respect to $r$.

\item[(ii)] For $a\in \Sigma_{0}$ with $f(a)=0$, there exist a
neighborhood $V$ of $a$ such that $f(z)\neq 0$ for all 
$z\in \bar{V}\setminus \{a\}$ and the Brouwer degree function
$d_{B}(f,V,0)\neq 0$, see  \cite[Appendix A]{BL} for the
definition of the Brouwer degree.

\item[(iii)] If $\partial h/\partial \theta \neq 0$, then for
all $(\theta,r) \in \mathcal{M}$ we have that 
$(\partial h/\partial \theta)(\theta,r)\neq 0$; and if 
$\partial h/\partial \theta\equiv 0$, then 
$\langle \nabla_{r}h, F_1\rangle^2-\langle \nabla_{r}h,
F_2\rangle^2>0$ for all $(\theta,z)\in [0,T]\times\mathcal{M}$,
where $\nabla_{r}h$ denotes the gradient of the function $h$
restricted to variable $r$.
\end{itemize}
Then for $|\varepsilon|>0$ sufficiently small, there exists a
$T$-periodic solution $r(\theta,\varepsilon)$ of system $\eqref{e4}$ such
that $r(0,\varepsilon)\to a$ (in the sense of Hausdorff
distance) as $\varepsilon\to 0$.
\end{theorem}

To verify the hypothesis (ii) of Theorem \ref{thm4}, we have the
following remark, see for instance \cite{BL}.

\begin{remark} \label{rmk5} \rm
Let $f: D\to \mathbb{R}$ be a $C^{1}$ function with
$f(a)=0$, where $D$ is an open subset of $\mathbb{R}$ and $a\in D$.
Whenever the Jacobian matrix of $J_{f}(a)\neq 0$ , there exists a
neighborhood $V$ of $a$ such that $f(r)\neq 0$ for all $r\in
\bar{V}\setminus \{a\}$. Then $d_{B}(f,V,0)\neq 0$.
\end{remark}

By Theorem \ref{thm4} and Remark \ref{rmk5}, if system \eqref{e4} satisfies the
hypothesis (i) and (iii), then every simple zero of the averaged
function $f(r)$ defined by \eqref{e6} provides a limit cycle of system \eqref{e4}.
In the next section, we will compute the averaged function case by
case.

\section{Averaged function}

Taking polar coordinates $x=r\cos\theta$, $y=r\sin\theta,
\theta\in(0,2\pi)$ and viewing $\theta$ as the new independent
variable. For $|\varepsilon|>0$ sufficiently small, we can
transform the differential system \eqref{e2} into the equivalent
differential equation
\begin{equation} \label{e7}
\frac{\mathrm{d}r}{\mathrm{d}\theta}
=\begin{cases}
\varepsilon X^{+}(\theta,r)+\varepsilon^2 Y^{+}(\theta, r,
\varepsilon),  &\sin\theta>0,\\
\varepsilon X^{-}(\theta, r)+\varepsilon^2 Y^{-}(\theta, r,
\varepsilon), & \sin\theta<0,
\end{cases} 
\end{equation}
where
\begin{equation} \label{e8}
\begin{gathered}
X^{+}(\theta, r)=\frac{\cos\theta p_1(\theta,r)+\sin\theta
q_1(\theta,r)}{R(\theta,r)},\\
Y^{+}(\theta, r, \varepsilon)=\frac{\sin\theta
p_1(\theta,r)-\cos\theta q_1(\theta,r)}{r R(\theta,r)
-\varepsilon\left(\sin\theta p_1(\theta,r)-\cos\theta q_1(\theta,r)\right)},\\
X^{-}(\theta, r)=\frac{\cos\theta p_2(\theta,r)+\sin\theta
q_2(\theta,r)}{R(\theta,r)},\\
Y^{-}(\theta, r, \varepsilon)=\frac{\sin\theta
p_2(\theta,r)-\cos\theta q_2(\theta,r)}{r R(\theta,r)
-\varepsilon\left(\sin\theta p_2(\theta,r)-\cos\theta
q_2(\theta,r)\right)},
\end{gathered}
\end{equation}
with
\begin{gather*}
R(\theta,r)=R(r\cos\theta,r\sin\theta),\\
p_i(\theta,r)=p_i(r\cos\theta,r\sin\theta),\quad i=1,2,\\
q_i(\theta,r)=q_i(r\cos\theta,r\sin\theta),\quad i=1,2.
\end{gather*}
Recall that $p_i(x,y),q_i(x,y)$, $i=1,2$ are given by \eqref{e3}.

Denote
\begin{gather*}
F_1(\theta,r)=\frac{1}{2}\left(X^{+}(\theta,r)+X^{-}(\theta,r)\right),\\
F_2(\theta,r)=\frac{1}{2}\left(X^{+}(\theta,r)-X^{-}(\theta,r)\right),\\
R_1(\theta,r,\varepsilon)=\frac{1}{2}\left(Y^{+}(\theta,r,\varepsilon)
 +Y^{-}(\theta,r,\varepsilon)\right),\\
R_2(\theta,r,\varepsilon)=\frac{1}{2}\left(Y^{+}(\theta,r,\varepsilon)
 -Y^{-}(\theta,r,\varepsilon)\right),
\end{gather*}
then differential equation \eqref{e7} becomes
\begin{equation} \label{e9}
\frac{\mathrm{d}r}{\mathrm{d}\theta}=\varepsilon
F(\theta,r)+\varepsilon^2R(\theta,r,\varepsilon),
\end{equation}
where
\begin{equation} \label{e10}
\begin{gathered}
F(\theta,r)=F_1(\theta,r)+\operatorname{sign}(\sin\theta)F_2(\theta,r),\\
R(\theta,r,\varepsilon)=R_1(\theta,r,\varepsilon)
+\operatorname{sign}(\sin\theta)R_2(\theta,r,\varepsilon).
\end{gathered}
\end{equation}

It is obvious that $F(\theta,r)$, $R(\theta,r,\varepsilon)$,
$h(\theta,r)=\sin\theta$ are locally L-Lipschitz with respect $r$.
Since $\mathcal{M}=\{(r,\theta)|\theta=0,\pi\}$, the function
$\partial h(\theta,r)/\partial\theta=\cos\theta\neq 0$ when
$(r,\theta)\in\mathcal {M}$. By Theorem \ref{thm4}, we need to
estimate the number of simple zeros for the averaged function
\begin{equation} \label{e11}
\begin{aligned}
f(r)
&=\int_{0}^{2\pi}F(\theta,r)\mathrm{d}\theta\\
&=\int_{0}^{\pi}X^{+}(\theta,r)\mathrm{d}\theta
+\int_{\pi}^{2\pi}X^{-}(\theta,r)\mathrm{d}\theta.
\end{aligned}
\end{equation}

Note that for $r>0$, the zeros of the averaged function $f(r)$
coincide with the zeros of the function $F(r)=rf(r)$. Hence, in
order to simplify further computation, we will compute $F(r)$ as
follows:
\begin{equation} \label{e12}
\begin{aligned}
F(r)
&=\int_{0}^{\pi}\frac{r\cos\theta p_1(\theta,r)
 +r\sin\theta q_1(\theta,r)}{R(\theta,r)}\mathrm{d}\theta\\
&\quad+\int_{\pi}^{2\pi}\frac{r\cos\theta p_2(\theta,r)
 +r\sin\theta q_2(\theta,r)}{R(\theta,r)}\mathrm{d}\theta\\
&=\int_{0}^{\pi}\frac{r\cos\theta p_1(r\cos\theta,r\sin\theta)
 +r\sin\theta q_1(r\cos\theta,r\sin\theta)}{R(r\cos\theta,r\sin\theta)}
 \mathrm{d}\theta\\
&\quad-\int_{0}^{\pi}\frac{r\cos\theta p_2(-r\cos\theta,-r\sin\theta)
 +r\sin\theta q_2(-r\cos\theta,-r\sin\theta)}{R(-r\cos\theta,-r\sin\theta)}
 \mathrm{d}\theta.
\end{aligned}
\end{equation}

In the following subsections, we will deduce the averaged function \eqref{e12} case by
case.

\subsection{ Ellipse-case (E) }
In this subsection, we study system \eqref{e2} with
$R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1$. Without loss
of generality, we can assume that $0<a<b$, otherwise we can change the
variables $x_1=y$, $y_1=x$, $t_1=-t$.
Substituting
\begin{equation} \label{e13}
\begin{aligned}
R(\theta,r)
&=\frac{(r\cos\theta)^2}{a^2}+\frac{(r\sin\theta)^2}{b^2}-1\\
&=-\frac{a^2(b^2-r^2)-(b^2-a^2)r^2\cos^2\theta}{a^2b^2}
\end{aligned}
\end{equation}
in \eqref{e12}, we have
\begin{equation} \label{e14}
\begin{aligned}
F(r)
&=\int_{0}^{\pi}g_1(\theta,r)\mathrm{d}\theta
 +g_2(r)\int_{0}^{\pi}\frac{1}{a^2(b^2-r^2)-(b^2-a^2)r^2
 \cos^2\theta}\mathrm{d}\theta\\
&\quad +g_3(r)\int_{0}^{\pi}\frac{\cos\theta}{a^2(b^2-r^2)-(b^2
 -a^2)r^2\cos^2\theta}\mathrm{d}\theta\\
&\quad +g_4(r)\int_{0}^{\pi}\frac{\sin\theta}{a^2(b^2-r^2)-(b^2
 -a^2)r^2\cos^2\theta}\mathrm{d}\theta\\
&\quad +g_5(r)\int_{0}^{\pi}\frac{\cos\theta\sin\theta}{a^2(b^2-r^2)
 -(b^2-a^2)r^2\cos^2\theta}\mathrm{d}\theta,
\end{aligned}
\end{equation}
where
\begin{align*}
g_1(\theta,r)
&=\frac{a^2b^2}{(a^2-b^2)^2}\Big(b^2(a_1-b_2+c_1-d_2) \\
&\quad+a^2(-a_1+b_2-c_1+d_2+a^2b^2(a_6-a_{8}-b_{7}+b_{9}+c_6-c_{8}-d_{7}+d_{9})\\
&\quad(b^2(a_{8}+b_{7}-2b_{9}+c_{8}+d_{7}-2d_{9})+ a^2(-a_6+b_{9}-c_6+d_{9}))r^2\\
&\quad-(a^2-b^2)(a_3-a_5-b_4-c_3+c_5+d_4)r\cos\theta\\
&\quad-(a^2-b^2)(a_6-a_{8}-b_{7}+b_{9}+c_6-c_{8}-d_{7}+d_{9})r^2\cos^2\theta\\
&\quad-(a^2-b^2)(a_4+b_3-b_5-c_4-d_3+d_5)r\sin\theta\\
&\quad-(a^2-b^2)(a_{7}-a_{9}+b_6-b_{8}+c_{7}
-c_{9}+d_6-d_{8})r^2\cos\theta\sin\theta)\Big),
\end{align*}
\begin{align*}
g_2(r)&=\frac{-a^2 b^2}{(a^2-b^2)^2}\Big(a^2b^4(a_1-b_2+c_1-d_2)
 +a^4b^2(-a_1+b_2-c_1+d_2)\\
&\quad+a^4b^4(a_6-a_{8}-b_{7}+b_{9}+c_6-c_{8}-d_{7}+d_{9})\\
&\quad+\Big(a^4(a_1+c_1)+b^4(b_2+d_2)-a^2b^2(a_1+b_2+c_1+d_2)\\
&\quad-a^4 b^2(2a_6-a_{8}-b_{7}+2c_6-c_{8}-d_{7})\\
&\quad+a^2b^4(a_{8}+b_{7}-2b_{9}+c_{8}+d_{7}-2d_{9})r^2\\
&\quad+(a^4(a_6+c_6)-a^2b^2(a_{8}+b_{7}+c_{8}+d_{7})+b^4(b_{9}
+d_{9}))r^4\Big)\Big),
\end{align*}
\begin{align*}
g_3(r)
&=\frac{a^2b^2r}{a^2-b^2}\Big(a^2b^2(a_3-a_5-b_4
+c_3-c_5-d_4)\\
&\quad-(a^2(a_3+c_3)+b^2(a_5+b_4+c_5+d_4))r^2\Big),
\end{align*}
\begin{align*}
g_4(r)
&=\frac{a^2b^2r}{a^2-b^2}\Big(a^2b^2(a_4+b_3-b_5-c_4-d_3+d_5) \\
&\quad-(a^2(a_4+b_3-c_4-d_3)-b^2(b_5-d_5))r^2\Big),
\end{align*}
\begin{align*}
g_5(r)
&=\frac{a^2b^2r^2}{a^2-b^2}\Big(b^2(a_2+b_1+c_2+d_1)+a^2(-a_2-b_1-c_2-d_1)\\
&\quad+a^2b^2(a_{7}-a_{9}+b_6-b_{8}+c_{7}-c_{9}+d_6-d_{8})\\
&\quad+(-a^2a_{7}+a_{9}b^2-a^2b_6+b^2b_{8}-a^2c_{7}
+b^2c_{9}-a^2d_6+b^2d_{8})r^2\Big).
\end{align*}
Since the calculation of \eqref{e14} is tedious, the whole calculations 
can be done by using some algebraic manipulator such as Mathematica or Maple.

To compute the average function \eqref{e14}, we list some useful results 
on integrals in the following lemma.

\begin{lemma} \label{lem6}
The following equalities hold:
\begin{gather}
\int_{0}^{\pi}\frac{1}{a^2(b^2-r^2)-(b^2
-a^2)r^2\cos^2\theta}\mathrm{d}\theta
=\frac{\pi}{ab\sqrt{a^2-r^2}\sqrt{b^2-r^2}};
\\
\int_{0}^{\pi}\frac{\cos\theta}{a^2(b^2-r^2)-(b^2
-a^2)r^2\cos^2\theta}\mathrm{d}\theta=0;
\\
\int_{0}^{\pi}\frac{\sin\theta}{a^2(b^2-r^2)
-(b^2-a^2)r^2\cos^2\theta}\mathrm{d}\theta
=\frac{\ln \big(\frac{a\sqrt{b^2-r^2}
+\sqrt{b^2-a^2}r}{a\sqrt{b^2-r^2}-\sqrt{b^2
-a^2}r}\big)}{ar\sqrt{b^2-a^2}\sqrt{b^2-r^2}};
\\
\int_{0}^{\pi}\frac{\cos\theta\sin\theta}{a^2(b^2-r^2)
-(b^2-a^2)r^2\cos^2\theta}\mathrm{d}\theta=0.
\end{gather}
\end{lemma}


Note that for $0<r<a<b$, the proof of above lemma is straightforward, 
we omit it here.
Applying Lemma \ref{lem6} to \eqref{e14}, we  obtain the averaged function
\begin{equation} \label{e15}
F(r)=k_1f_1(r)+k_2f_2(r)+k_3f_3(r)+k_4f_4(r)+k_5f_5(r)+k_6f_6(r),
\end{equation}
where $r\in(0,a)$,
\begin{equation} \label{e16}
\begin{gathered}
f_1(r)=\frac{ab}{\sqrt{a^2-r^2}\sqrt{b^2-r^2}}-1,\quad
f_2(r)=\frac{r^2}{\sqrt{a^2-r^2}\sqrt{b^2-r^2}},\\
f_3(r)=\frac{r^{4}}{\sqrt{a^2-r^2}\sqrt{b^2-r^2}},\quad
f_4(r)=r^2,\\
f_5(r)=\frac{ab^2}{\sqrt{b^2-r^2}}\ln
\Big(\frac{a\sqrt{b^2-r^2}+\sqrt{b^2-a^2}r}{a\sqrt{b^2-r^2}
 -\sqrt{b^2-a^2}r}\Big)-2\sqrt{b^2-a^2}r,\\
f_6(r)=\frac{r^2}{\sqrt{b^2-r^2}}\ln
 \Big(\frac{a\sqrt{b^2-r^2}+\sqrt{b^2-a^2}r}{a\sqrt{b^2-r^2}
-\sqrt{b^2-a^2}r}\Big),
\end{gathered}
\end{equation}
and
\begin{equation} \label{e17}
\begin{aligned}
k_1&=\frac{-a^2b^2\pi}{(b^2-a^2)^2}\Big(b^2-a^2)(a_1+c_1-b_2-d_2)\\
&\quad+a^2b^2(a_6+c_6+b_{9}-d_{9}-a_{8}-c_{8}-b_{7}-d_{7})\Big),\\
k_2&= \frac{-ab\pi}{(b^2-a^2)^2}\Big(a^2b^2(a^2+b^2)(a_{8}
 +b_{7}+c_{8}+d_{7})-a^2(b^2-a^2)(a_1+c_1)\\
&\quad +b^2(b^2-a^2)(b_2+d_2)-2a^{4}b^2(a_6+c_6)
 -2a^2b^{4}(b_{9}+d_{9})\Big),\\
k_3&= \frac{-ab\pi}{(b^2-a^2)^2}
 \Big(a^{4}(a_6+c_6)+b^{4}(b_{9}+d_{9})-a^2b^2(a_{8}+b_{7}+c_{8}+d_{7})
 \Big),\\
k_4&= \frac{a^2b^2\pi}{2(b^2-a^2)^2}
 \Big((a^2+b^2)(a_{8}+c_{8}+b_{7}+d_{7})\\
&\quad -(3a^2-b^2)(a_6+c_6)-(3b^2-a^2)(b_{9}+d_{9})\Big),\\
k_5&= \frac{a^2b^2(a_4+b_3-b_5-c_4-d_3
 +d_5)}{(b^2-a^2)^{3/2}},\\
k_6&= \frac{ab^2}{(b^2-a^2)^{3/2}}
\big(a^2(a_4+b_3-c_4-d_3)-b^2(b_5-d_5)\big).
\end{aligned}
\end{equation}
From above analysis, we have the following proposition.

\begin{proposition} \label{prop7}
Suppose that $R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1$,
then the averaged function \eqref{e15} are linear combination of 6
linearly independent functions $f_i(r)$, $i=1,2,3,4,5,6$ defined by
\eqref{e16}, and the coefficients $k_i$, $i=1,2,3,4,5,6$ given by \eqref{e17}
can be chosen arbitrarily.

Moreover, if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then averaged function 
\eqref{e15} are linear combination of 4 linearly
independent functions $f_i(r)$, $i=1,2,3,4$ defined by \eqref{e16}, 
and the coefficients $k_i$, $i=1,2,3,4$ given by \eqref{e17} can be chosen arbitrarily.
\end{proposition}

\begin{proof}
It is obvious that the averaged functions $F(r)$ of \eqref{e15} is a linear
combination of the 6 functions defined by \eqref{e16}. In the following we
will prove these functions are linearly independent.
Suppose that
\begin{equation} \label{e18}
F(r)=\bar{k}_1f_1(r)+\bar{k}_2f_2(r)+\bar{k}_3f_3(r)+\bar{k}_4f_4(r)
+\bar{k}_5f_5(r)+\bar{k}_6f_6(r)=0.
\end{equation}
To prove that $f_i(r)$, $i=1,2,3,4,5,6$ defined by \eqref{e16} are linearly independent,
we need to prove that the coefficients $\bar{k}_i=0$ for $i=1,2,3,4,5,6$.

Derivative \eqref{e18} with respect to $r$ at $r=0$, we have
\begin{align*}
F^{(2)}(0)
&= \frac{(a^2+b^2)\bar{k}_1+2ab\bar{k}_2+2a^2b^2\bar{k}_4}{a^2b^2}=0,\\
F^{(3)}(0)
&= \frac{4\sqrt{b^2-a^2}\left((2a^2+b^2)\bar{k}_5+3a\bar{k}_6\right)}
 {a^2b^2}=0,\\
F^{(4)}(0)
 &= \frac{3\left((3a^{4}+2a^2b^2+3b^{4})\bar{k}_1+4(a^{3}b+ab^{3})\bar{k}_2
 +8a^{3}b^{3}\bar{k}_3\right)}{a^{4}b^{4}}=0,\\
F^{(5)}(0)&= \frac{16\sqrt{b^2-a^2}}{a^{4}b^{4}}
\left((8a^{4}+4a^2b^2+3b^{4})\bar{k}_5+(10a^{3}+5ab^2)\bar{k}_6\right)=0,\\
F^{(6)}(0)&= \frac{45(5a^{6}+3a^{4}b^2+3a^2b^{4}+5b^{6})}{a^{6}b^{6}}\bar{k}_1\\
&\quad +\frac{45\left((6a^{5}b+4a^{3}b^{3}+6ab^{5})\bar{k}_2
 +(8a^{5}b^{3}+a^{3}b^{5})\bar{k}_3\right)}{a^{6}b^{6}}=0,\\
F^{(8)}(0)&= \frac{315\left(35a^8+20a^6b^2+18a^4b^4+20a^2b^6+35b^8\right)}
{a^{8}b^{8}}\bar{k}_1\\
&\quad +\frac{315(40a^7b+24a^5b^3+24a^3b^5+40ab^7)}{a^{8}b^{8}}\bar{k}_2\\
&\quad +\frac{315(48a^7b^3 +32a^5b^5+48a^3b^7)}{a^{8}b^{8}}\bar{k}_3=0.
\end{align*}
Since $0<a<b$, the determinant of coefficient matrix,
\begin{equation*}\begin{array}{ll}
\det \frac{\partial(F^{(2)}(0),F^{(3)}(0),F^{(4)}(0),F^{(5)}(0),F^{(6)}(0),
F^{(8)}(0))}{\partial(\bar{k}_1,\bar{k}_2,\bar{k}_3,\bar{k}_4,
\bar{k}_5,\bar{k}_6)}=\frac{348364800(b^2-a^2)^9}{a^{20}b^{20}}\neq0,
\end{array}\end{equation*}
so that the only solution is the trivial one. From the above analysis we can obtain $\bar{k}_i=0,i=1,2,3,4,5,6$,
thus we have proved that these functions $f_i(r),i=1,2,3,4,5,6$ defined by \eqref{e16} are linearly independent.

In the following we will show that the coefficients $k_i,i=1,2,3,4,5,6$ given by \eqref{e17} can be chosen arbitrarily.

From \eqref{e17}, the determinant of Jacobian matrix
\[
\det \frac{\partial(k_1,k_2,k_3,k_4,k_5,k_6)}
{\partial(a_1,b_2,b_{7},b_{9},b_3,b_5)}
=\frac{a^9b^{12}\pi^{4}}{2(b^2-a^2)^{5}}\neq 0.
\]
The coefficients $a_1, b_2, b_3, b_5, b_{7}, b_{9}$ can be chosen
so that the $k_i$, $i=1,2,3,4,5,6$ can take any specified values.

Note that if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then
$k_5=k_6=0$. The averaged function \eqref{e15}
are linear combination of $f_i(r)$, $i=1,2,3,4$ defined by \eqref{e17}.
\end{proof}

\subsection{Complex-Ellipse case (CE)}
In this subsection, we study system \eqref{e2} with
$R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}+1$.
Substituting
\[
R(\theta,r)=\frac{(r\cos\theta)^2}{a^2}+\frac{(r\sin\theta)^2}{b^2}+1
\]
in \eqref{e12}, similar to the case (E), we have
\begin{equation} \label{e19}
F(r)
= k_1f_1(r)+k_2f_2(r)+k_3f_3(r)+k_4f_4(r)+k_5f_5(r)+k_6f_6(r),
\end{equation}
where $r\in(0,+\infty)$,
\begin{equation} \label{e20}
\begin{gathered}
f_1(r)=\frac{ab}{\sqrt{a^2+r^2}\sqrt{b^2+r^2}}-1,\quad
f_2(r)=r^2,\\
f_3(r)=\frac{r^2}{\sqrt{a^2+r^2}\sqrt{b^2+r^2}},\quad
f_4(r)=\frac{r^{4}}{\sqrt{a^2+r^2}\sqrt{b^2+r^2}},\\
f_5(r)=\frac{ab^2}{\sqrt{b^2+r^2}}\arctan
\Big(\frac{\sqrt{b^2-a^2}r}{a\sqrt{b^2+r^2}}\Big)-\sqrt{b^2-a^2}r,\\
f_6(r)=\frac{r^2}{\sqrt{b^2+r^2}}\arctan
\Big(\arctan\frac{\sqrt{b^2-a^2}r}{a\sqrt{b^2+r^2}}\Big),
\end{gathered}
\end{equation}
and
\begin{equation} \label{e21}
\begin{aligned}
k_1&= \frac{a^2b^2\pi}{(b^2-a^2)^2}\Big((b^2-a^2)(a_1+c_1-b_2-d_2) \\
&\quad +a^2b^2(a_6+c_6+b_{9}-d_{9}-a_{8}-c_{8}-b_{7}-d_{7})\Big),\\
k_2&= \frac{a^2b^2\pi}{2(b^2-a^2)^2}\Big((a^2+b^2)(a_{8}+c_{8}+b_{7}+d_{7}) \\
&\quad -(3a^2-b^2)(a_6+c_6)-(3b^2-a^2)(b_{9}+d_{9})\Big),\\
k_3&= \frac{-ab\pi}{(b^2-a^2)^2}\Big(a^2b^2(a^2+b^2)(a_{8}+b_{7}+c_{8}+d_{7})
 -a^2(b^2-a^2)(a_1+c_1) \\
&\quad +b^2(b^2-a^2)(b_2+d_2)-2a^{4}b^2(a_6+c_6)-2a^2b^{4}(b_{9}+d_{9})\Big),\\
k_4&= \frac{-ab\pi}{(b^2-a^2)^2}\Big(a^{4}(a_6+c_6)+b^{4}(b_{9}+d_{9})
 -a^2b^2(a_{8}+b_{7}+c_{8}+d_{7})\Big),\\
k_5&= \frac{-a^2b^2(a_4+b_3-b_5-c_4-d_3+d_5)}{(b^2-a^2)^{3/2}},\\
k_6&= \frac{-ab^2}{(b^2-a^2)^{3/2}}\Big(a^2(a_4+b_3-c_4-d_3)-b^2(b_5-d_5)\Big).
\end{aligned}
\end{equation}
From above analysis, we have the following proposition.

\begin{proposition} \label{prop8}
Suppose that $R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}+1$, then the averaged 
function \eqref{e19} are linear combination of the 6
linearly independent functions $f_i(r)$, $i=1,2,3,4,5,6$ defined by
\eqref{e20}, and the coefficients $k_i$, $i=1,2,3,4,5,6$ given by \eqref{e21} 
can be chosen arbitrarily.

Moreover, if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then the averaged 
function \eqref{e19} are linear combination of the 4 linearly independent 
functions $f_i(r)$, $i=1,2,3,4$ defined by \eqref{e20}, and the coefficients 
$k_i$, $i=1,2,3,4$ given by \eqref{e21} can be chosen arbitrarily.
\end{proposition}


The proof of the above proposition is similar to the proof of Proposition \ref{prop7},
 we omit it here.


\subsection{Hyperbolic case (H)}

In this subsection, we study system \eqref{e2} with
$R(x,y)=\frac{x^2}{a^2}-\frac{y^2}{b^2}-1$.
Substituting
\[
R(\theta,r)=\frac{(r\cos\theta)^2}{a^2}-\frac{(r\sin\theta)^2}{b^2}-1
\]
in \eqref{e12}, we have
\begin{equation} \label{e22}
F(r)= k_1f_1(r)+k_2f_2(r)+k_3f_3(r)+k_4f_4(r)+k_5f_5(r)+k_6f_6(r),
\end{equation}
where $r\in(0,a)$,
\begin{equation} \label{e23}
\begin{gathered}
f_1(r)=\frac{ab}{\sqrt{a^2-r^2}\sqrt{b^2+r^2}}-1,\quad f_2(r)=r^2,\\
f_3(r)=\frac{r^2}{\sqrt{a^2-r^2}\sqrt{b^2+r^2}},\quad
f_4(r)=\frac{r^{4}}{\sqrt{a^2-r^2}\sqrt{b^2+r^2}},\\
f_5(r)=\frac{ab^2}{\sqrt{b^2+r^2}}\ln 
\Big(\frac{a\sqrt{b^2+r^2}-\sqrt{a^2+b^2}r}{a\sqrt{b^2+r^2}
+\sqrt{a^2+b^2}r}\Big)+2\sqrt{b^2+a^2}r,\\
f_6(r)=\frac{r^2}{\sqrt{b^2+r^2}}\ln 
\Big(\frac{a\sqrt{b^2+r^2}-\sqrt{a^2+b^2}r}{a\sqrt{b^2+r^2}+\sqrt{a^2+b^2}r}\big),
\end{gathered}
\end{equation}
and
\begin{equation} \label{e24}
\begin{aligned}
k_1
&= \frac{-a^2b^2\pi}{(b^2+a^2)^2}\Big((b^2-a^2)(a_1+c_1-b_2-d_2) \\
&\quad +a^2b^2(a_6+c_6+b_{9}-d_{9}-a_{8}-c_{8}-b_{7}-d_{7})\Big),\\
k_2&= \frac{a^2b^2\pi}{2(b^2+a^2)^2}\Big((a^2+b^2)(a_{8}+c_{8}+b_{7}+d_{7}) \\
&\quad -(3a^2-b^2)(a_6+c_6)-(3b^2-a^2)(b_{9}+d_{9})\Big),\\
k_3&= \frac{-ab\pi}{(b^2+a^2)^2}\Big(a^2b^2(a^2+b^2)(a_{8}+b_{7}+c_{8}+d_{7})
 -a^2(b^2-a^2)(a_1+c_1) \\
&\quad +b^2(b^2-a^2)(b_2+d_2)-2a^{4}b^2(a_6+c_6)-2a^2b^{4}(b_{9}+d_{9})\Big),\\
k_4&= \frac{-ab\pi}{(b^2+a^2)^2}\left(a^{4}(a_6+c_6)+b^{4}(b_{9}+d_{9})
 -a^2b^2(a_{8}+b_{7}+c_{8}+d_{7})\right),\\
k_5&= \frac{a^2b^2(a_4+b_3-b_5-c_4-d_3+d_5)}{(b^2+a^2)^{3/2}},\\
k_6&= \frac{ab^2}{(b^2+a^2)^{3/2}}\left(a^2(a_4+b_3-c_4-d_3)-b^2(b_5-d_5)\right).
\end{aligned}
\end{equation}
From above analysis, we have the following proposition.

\begin{proposition} \label{prop9}
Suppose that $R(x,y)=\frac{x^2}{a^2}-\frac{y^2}{b^2}-1$, then the averaged 
function \eqref{e22} are linear combination of the 6
linearly independent functions $f_i(r)$, $i=1,2,3,4,5,6$ defined by \eqref{e23}, 
and the coefficients $k_i$, $i=1,2,3,4,5,6$ given by \eqref{e24} can be chosen 
arbitrarily.

Moreover, if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then the averaged function 
\eqref{e22} are linear combination of the 4 linearly independent functions 
$f_i(r)$, $i=1,2,3,4$ defined by \eqref{e23}, and the coefficients 
$k_i$, $i=1,2,3,4$ given by \eqref{e24} can be chosen arbitrarily.
\end{proposition}


The proof of the above proposition is similar to the proof of Proposition \ref{prop7}, 
we omit it here.


\subsection{Two complex straight lines intersecting in a real point case (CL)}

In this subsection, we study system \eqref{e2} with
$R(x,y)=a^2x^2+y^2$. Without loss of generality, we can assume
that $a>1$, Otherwise, we change the variables
$x=y_1,y=x_1,t_1=-a^2t$.

Substituting $R(\theta,r)=r^2(a^2\cos^2\theta+\sin^2\theta)$
in \eqref{e12}, we have
\begin{equation} \label{e25}
F(r)= k_1+k_2r+k_3r^2,
\end{equation}
where $r\in(0,+\infty)$ and
\begin{equation} \label{e26}
\begin{aligned}
k_1&= \frac{(a_1+c_1+ab_2+ad_2)\pi}{a(a+1)},\\
k_2&= \frac{-2(a_4+b_3-b_5-c_4-d_3+d_5)}{a^2-1}\\
&\quad +\frac{2(a_4+b_3-a^2b_5-c_4-d_3+a^2d_5)}{a^2-1}
\arctan (\sqrt{a^2-1}t),\\
k_3&= \frac{\pi((2+a)(a_6+c_6)+a(a_{8}+b_{7}+c_{8}+d_{7})+2a^2(b_{9}
+d_{9}))}{2a(a+1)^2}.
\end{aligned}
\end{equation}
From above analysis, we have the following proposition.

\begin{proposition} \label{prop10}
Suppose that $R(x,y)=a^2x^2+y^2$, then the averaged function \eqref{e25} 
are linear combination of the 3
linearly independent functions $1,r,r^2$, and the coefficients $k_i$,
$i=1,2,3$ given by \eqref{e26} can be chosen arbitrarily.

Moreover, if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then the averaged 
function \eqref{e25} are linear combination of the 2 linearly 
independent functions $1,r^2$, and the coefficients $k_i$, $i=1,3$ 
given by \eqref{e26} can be chosen arbitrarily.
\end{proposition}

\begin{proof}
It is obvious that functions $1,r,r^2$ are linearly independent.
The determinant of Jacobian matrix
\begin{equation*}\begin{array}{ll}
\det \frac{\partial(k_1,k_2,k_3)}{\partial(a_1,a_4,a_6)}
=\frac{(a+2)\pi^2(1-\arctan \sqrt{a^2-1})}{a^2(a-1)(a+1)^{4}}\neq
0.
\end{array}\end{equation*}
It is obvious that the coefficients $k_i$, $i=1,2,3$ given by \eqref{e26}
can be chosen arbitrarily since $a_1,a_4,a_6$ are arbitrary.
\end{proof}

\subsection{Two real parallel straight lines case (RPL)}

In this subsection, we study system \eqref{e2} with $R(x,y)=x^2-a^2$.
Substituting $R(\theta,r)=r^2\cos^2\theta-a^2$ in \eqref{e12}, we
have
\begin{equation} \label{e27}
F(r)= k_1f_1(r)+k_2f_2(r)+k_3f_3(r)+k_4f_4(r)+k_5f_5(r)+k_6f_6(r),
\end{equation}
where $r\in(0,a)$,
\begin{equation} \label{e28}
\begin{gathered}
f_1(r)= \frac{a}{\sqrt{a^2-r^2}}-1,\quad
f_2(r)= a-\sqrt{a^2-r^2},\\
f_3(r)= r^2,\quad
f_4(r)= r^2\sqrt{a^2-r^2},\\
f_5(r)= 2r-a\ln \Big(\frac{a+r}{a-r}\Big),\quad
f_6(r)= r^2\ln \Big(\frac{a+r}{a-r}\Big),
\end{gathered}
\end{equation}
and
\begin{equation} \label{e29}
\begin{aligned}
k_1&= -(a_1+c_1+a^2a_6+a^2c_6)\pi,\\
k_2&= \frac{-\pi}{a}(b_2+d_2+a^2(a_{8}+b_{7}-b_{9}+c_{8}+d_{7}-d_{9})),\\
k_3&= \frac{a_6+a_{8}+b_{7}-3b_{9}+c_6+c_{8}+d_{7}-3d_{9}}{2},\\
k_4&= \frac{b_{9}+d_{9}}{a},\\
k_5&= a_4+b_3-b_5-c_4-d_3+d_5,\\
k_6&= \frac{-b_5+d_5}{a}.
\end{aligned}
\end{equation}
From above analysis, we have the following proposition.

\begin{proposition} \label{prop11}
Suppose that $R(x,y)=x^2-a^2$, then the averaged function \eqref{e27}
 are linear combination of the 6 linearly independent functions 
$f_i(r)$, $i=1,2,3,4,5,6$ defined by \eqref{e28}, and the coefficients 
$k_i$, $i=1,2,3,4,5,6$ given by \eqref{e29} can be chosen arbitrarily.

Moreover, if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then the averaged 
function \eqref{e27} are linear combination of the 4 linearly independent
 functions $f_i(r)$, $i=1,2,3,4$ defined by \eqref{e28}, and the coefficients 
$k_i,i=1,2,3,4$ given by \eqref{e29} can be chosen arbitrarily.
\end{proposition}

The proof of the above proposition is similar to that of Proposition \ref{prop7}, 
we omit it here.

\subsection{Two complex parallel straight lines case (CPL)}
In this subsection, we study system \eqref{e2} with $R(x,y)=x^2+a^2$.
Substituting $R(\theta,r)=r^2\cos^2\theta+a^2$ in \eqref{e12}, we
have
\begin{equation} \label{e30}
F(r)= k_1f_1(r)+k_2f_2(r)+k_3f_3(r)+k_4f_4(r)+k_5f_5(r)+k_6f_6(r),
\end{equation}
where $r\in(0,+\infty)$,
\begin{equation} \label{e31}
\begin{gathered}
f_1(r)= \frac{a}{\sqrt{a^2+r^2}}-1,\quad
f_2(r)= \frac{r^2}{\sqrt{a^2+r^2}},\\
f_3(r)= r^2,\quad
f_4(r)= \frac{r^{4}}{\sqrt{a^2+r^2}},\\
f_5(r)= r-a\arctan \big(\frac{r}{a}\big),\quad 
f_6(r)= r^2\arctan \big(\frac{r}{a}\big),
\end{gathered}
\end{equation}
and
\begin{equation} \label{e32}
\begin{aligned}
k_1&= -(a_1-b_2+c_1-d_2-a^2(a_6-a_{8}-b_{7}+b_{9}+c_6-c_{8}-d_{7}+d_{9}))\pi,\\
k_2&= \frac{\pi}{a}(b_2+d_2-a^2(a_{8}+b_{7}-2b_{9}+c_{8}+d_{7}-2d_{9})),\\
k_3&= \frac{a_6+a_{8}+b_{7}-3b_{9}+c_6+c_{8}+d_{7}-3d_{9}}{2},\\
k_4&= \frac{b_{9}+d_{9}}{a},\\
k_5&= 2(a_4+b_3-b_5-c_4-d_3+d_5),\\
k_6&= \frac{2(b_5-d_5)}{a}.
\end{aligned}
\end{equation}
From above analysis, we have the following proposition.

\begin{proposition} \label{prop12}
Suppose that $R(x,y)=x^2+a^2$, then the averaged function \eqref{e30} are linear
combination of the 6 linearly independent functions $f_i(r)$, $i=1,2,3,4,5,6$ 
defined by \eqref{e31}, and the coefficients $k_i$, $i=1,2,3,4,5,6$ given by 
\eqref{e32} can be chosen arbitrarily.

Moreover, if $a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$, then the averaged function 
\eqref{e30} are linear combination of the 4 linearly independent 
functions $f_i(r)$, $i=1,2,3,4$ defined by \eqref{e31}, and the coefficients 
$k_i$, $i=1,2,3,4$ given by \eqref{e32} can be chosen arbitrarily.
\end{proposition}

The proof of the above proposition is to the proof of Proposition \ref{prop7}, 
we omit it here.

\section{Proof of the main results}

The tools for estimating the number of zeros for the averaged
function $F(r)$ are the following two lemmas.

\begin{lemma}[\cite{CGP1}] \label{lem13}
 Consider $p+1$ linearly independent analytical
functions $f_i: U \to \mathbb{R}$, $i=0,1,\dots,p$, where
$U \in \mathbb{R}$ is an interval. Suppose that there exists 
$j \in \{0,1,\dots,p\}$ such that $f_j$ has constant sign. then there
exists $p+1$ constants $C_i, i=0,1,\dots,p$  such that
$f(x)=\sum^{p}_{i=0}C_if_i(x)$ has at least $p$ simple
zeros in $U$.
\end{lemma}

\begin{lemma}[\cite{LZLZ}] \label{lem14}
Suppose that the real function $F(h)\in C^1(U)$ satisfies
\begin{equation} \label{e33}
\alpha(h)F'(h)-\beta(h)F(h)=\gamma(h),
\end{equation}
where $\alpha(h),\beta(h),\gamma(h)\in C^0(U)$ and
$U \in \mathbb{R}$ is an interval, then
\begin{equation} \label{e34}
\#\{F(h)\}\leqslant 1+\#\{\alpha(h)\}+ \#\{\gamma(h)\}.
\end{equation}
Here $\#\{F(h)\}$ is the number of zeros of $F(h)$ in $U$.
\end{lemma}

Although the proof of Lemma \ref{lem14} can be found in \cite{LZLZ},
 we include it here for the sake of convenience.

\begin{proof}
Let $h_1$ and $h_2$ are two consecutive zeros of $F(h)$. From \eqref{e33}, we have
\begin{equation*}\label{eq1}\begin{array}{ll}
\alpha(h_1)F'(h_1)=\gamma(h_1),\quad 
\alpha(h_2)F'(h_2)=\gamma(h_2),
\end{array}\end{equation*}
which shows $\gamma(h_1)\gamma(h_2)\leqslant 0$ if $\alpha(h)\neq0$ in 
$(h_1-\varepsilon,h_2+\varepsilon)$, $0<\varepsilon<<1$. Therefore, between
any two consecutive zeros of $F(h)$, there must exist at least one zero 
of $\alpha(h)$ or $\gamma(h)$. The result of this lemma follows.
\end{proof}


\begin{proof}[Proof of Theorem \ref{thm1}]
(I). From \eqref{e25}, $F(r)$ is a polynomial in variable $r$
with the degree of $2$, since the coefficients are arbitrary, then
we can conclude that $F(r)$ has at most $2$ zeros in the interval
$(0,+\infty)$ and this bound is sharp. From the first order
averaging method of Theorem \ref{thm4}, system \eqref{e2} has at most $2$ limit
cycles bifurcating from the period annulus surrounding the origin.
Moreover, there are system \eqref{e2} has $2$ limit cycles bifurcating from
the period annulus surrounding the origin. Thus we have $H(3)=2$,
this complete the proof of statement (I).

(II). For the case
$R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1$. First we
consider the lower bound of zeros for $F(r),r\in(0,r_{0})$, in this
case $r_{0}=a$. By the Proposition \ref{prop7}, $F(r)$ is an linear
combination of the $6$ independent generating functions
$f_i(r),i=1,2,3,4,5,6$ defined by \eqref{e16} with arbitrary coefficients.
All these functions are analytic in $(0, r_{0})$, and some of them
are strictly positive in this interval, hence $F(r)$ can have at least
$5$ simple zeros in $(0, r_{0})$ by Lemma \ref{lem13}.

Now we  find an upper bound for the zeros for $F(r)$ in
$(0,r_{0})$. From \eqref{e15}, we have
\begin{align*}
F(r)&=k_1f_1(r)+k_2f_2(r)+k_3f_3(r)+k_4f_4(r)+k_5f_5(r)+k_6f_6(r)\\
&=\frac{1}{\sqrt{a^2-r^2}\sqrt{b^2-r^2}}G(r),
\end{align*}
where
\begin{equation} \label{e35}
\begin{aligned}
G(r)&= abk_1+k_2r^2+k_3r^{4}\\
&\quad +\sqrt{a^2-r^2}\sqrt{b^2-r^2}\left(-k_1-2\sqrt{b^2-a^2}k_5r+k_4r^2\right)\\
&\quad +\sqrt{a^2-r^2}\ln 
\Big(\frac{a\sqrt{b^2-r^2}+\sqrt{b^2-a^2}r}{a\sqrt{b^2-r^2}-\sqrt{b^2-a^2}r}\Big)
\left(ab^2k_5+k_6r^2\right),
\end{aligned}
\end{equation}
Recall that $r\in(0,a)$ and $a<b$, then
$(a^2-r^2)(b^2-r^2)\neq 0$. The zeros of $F(r)=0$ is
coincide with $G(r)=0$ in the interval $(0,a)$. In the following we
will discuss the upper bound of zeros for $G(r)=0, r\in(0,a)$.

Computing the derivative on both sides of \eqref{e35} with respect to $r$, we
obtain
\begin{equation} \label{e36}
\begin{aligned}
G'(r)
&= \frac{r}{\sqrt{a^2-r^2}\sqrt{b^2-r^2}}P_4(r)+2k_2+4k_3r^2\\
&\quad -\frac{r}{\sqrt{b^2-r^2}}\ln 
\Big(\frac{a\sqrt{b^2-r^2}+\sqrt{b^2-a^2}r}{a\sqrt{b^2-r^2}-\sqrt{b^2-a^2}r}\Big)
(ab^2k_5-2a^2k_6+3k_6r^2),
\end{aligned}
\end{equation}
where
\begin{equation} \label{e37}
\begin{aligned}
P_4(r)&= a^2k_1+b^2k_1+2a^2b^2k_4+\sqrt{b^2-a^2}(4a^2k_5+4b^2k_5+2ak_6)r\\
&\quad -(2k_1+3a^2k_4+3b^2k_4)r^2-6\sqrt{b^2-a^2}k_5r^3+4k_4r^4
\end{aligned}
\end{equation}
is a polynomial in the variable $r$, of degree 4.

Applying \eqref{e35}, we can eliminate
$\ln \Big(\frac{a\sqrt{b^2-r^2}+\sqrt{b^2-a^2}r}{a\sqrt{b^2-r^2}-\sqrt{b^2-a^2}r}
\Big)$ in \eqref{e36}, then we have
\begin{equation} \label{e38}
\alpha(r)G'(r)=\beta(r)G(r)+\gamma(r),
\end{equation}
where
\begin{gather*}
\alpha(r)= \sqrt{a^2-r^2}\sqrt{b^2-r^2}(ab^2k_5+k_6r^2),\\
\beta(r)= -r(ab^2k_5-2a^2k_6+3k_6r^2),\\
\gamma(r)= P_{7}(r)+\sqrt{a^2-r^2}\sqrt{b^2-r^2}P_5(r)
\end{gather*}
with
\begin{gather*}
P_{7}(r)= r(ab^2k_5+k_6r^2)P_4(r)+r(ab^2k_5-2a^2k_6+3k_6r^2)(abk_1-k_2r^2-k_3r^{4}),\\
\begin{aligned}
P_5(r)&= (2k_2+4k_3r^2)(ab^2k_5+k_6r^2)\\
&\quad-r(ab^2k_5-2a^2k_6+3k_6r^2)(k_1+2\sqrt{b^2-a^2}k_5r-k_4r^2).
\end{aligned}
\end{gather*}

From Lemma \ref{lem14}, we obtain an upper bound for the zeros of
$F(r)$ in $(0,r_{0})$:
\begin{equation} \label{e39}
\#\{F(r)\}\leqslant 1+\#\{\alpha(r)\}+\#\{\gamma(r)\}
\leqslant 1+1+14 =16.
\end{equation}
Note that $F(0)=0$, then the number of zeros in $(0,r_{0})$ is at
most $15$.

From the first order averaging method of Theorem \ref{thm4},  system \eqref{e2}
with $R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1$ can
bifurcate at most $15$ limit cycles from the period annulus
surrounding the origin, and there are system \eqref{e2} with at least $5$
limit cycles bifurcating from such annulus. Thus we have 
$5\leqslant H(3)\leqslant 15$.
The proof of the other statements are similar.
\end{proof}

\begin{proof}[Proof of Theorem \ref{thm2}]

(I) From \eqref{e26}, we know that $k_2=0$ since
$a_i=c_i,b_i=d_i,i=1,2,\dots,9$, then
$F(r)=k_1+k_3r^2$. It is obvious that $F(r)$ has at most $1$
zero in the interval $(0,+\infty)$ and this bound is sharp. From the
first order averaging method of Theorem \ref{thm4}, system \eqref{e2} has at most
$1$ limit cycles bifurcating from the period annulus surrounding the
origin. Moreover, there are system \eqref{e2} has $1$ limit cycles
bifurcating from the period annulus surrounding the origin.

(II) For the case
$R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1$. From \eqref{e17}, we
know that $k_5=k_6=0$ since
$a_i=c_i$, $b_i=d_i$, $i=1,2,\dots,9$. By the Proposition \ref{prop7},
$F(r)$ is an linear combination of the $4$ independent generating
functions $f_i(r)$, $i=1,2,3,4$ given by \eqref{e16} with arbitrary
coefficients. All these functions are analytic in $(0, r_{0})$, and
some of them are strictly positive in this interval, hence $F(r)$
has at least $3$ simple zeros in $(0, r_{0})$ by Lemma \ref{lem13}.

Replacing $k_5=k_6=0$ in \eqref{e15}, we have
\begin{equation} \label{e40}
F(r)=\frac{1}{\sqrt{a^2-r^2}\sqrt{b^2-r^2}}(abk_1+k_2r^2+k_3r^{3})-(k_1-k_4r^2),
\end{equation}
where $k_i$, $i=1,2,3,4$ are given by \eqref{e14}. From \eqref{e40}, it is obvious
that $F(r)$ has at most 4 zeros in $[0,+\infty$. Note that $F(0)=0$,
then $F(r)$ has at 3 zeros in $(0,+\infty)$.

From the first order averaging method of Theorem \ref{thm4},  systems \eqref{e2}
with $R(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1$ can
bifurcate at most $3$ limit cycles from the period annulus
surrounding the origin, Moreover, there are system \eqref{e2} have $3$
limit cycles bifurcating from such annulus.
The proof of the other case is similar.
\end{proof}

\section{Conclusion}
Our work described bounds for the number of limit cycles for 
piecewise smooth perturbation of cubic polynomial differential center. 
Our results confirm that piecewise smooth system can have more limit cycles 
than the smooth one. One of the major difficulties encountered in our study 
is how to estimate the number of zeros for the averaged function. 
Some of the basic tools we have used are the methods of Lemmas \ref{lem13}
 and \ref{lem14}. 
 We can obtain a upper bound of the number of zeros for averaged function 
by Lemma \ref{lem14}, but this bound is not sharp. More advanced methods, such as 
Argument principle, are required to attack this kind of problems. 
We leave this as the future research problems.


\section{Acknowledgements}
We want to thank the anonymous referees for their useful suggestions and valuable
comments. The first author is supported by the NSFC (No. 11401111). 
The second author is supported by NSFC (No. 11401531).

\begin{thebibliography}{99}

\bibitem{BD}  M. di Bernardo, C. J. Budd,  A. R. Champneys, 
P. Kowalczyk, A. B. Nordmark, G. Olivar Tost, P. T. Piiroinen;
\emph{Bifurcations in nonsmooth dynamical systems}.
 SIAM Rev. 50 (2008), no. 4, 629–701. 

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

\bibitem{BPT}  C. Buzzi, C. Pessoa, J. Torregrosa;
\emph{Piecewise linear perturbations of a linear center}, 
Discrete Contin. Dyn. Syst., 33 (2013), 3915--3936.

\bibitem{CD}  X. Chen, Z. Du;
\emph{Limit cycles bifurcate from centers of
discontinuous quadratic systems}, Computers and Mathematics with
Applications, 59 (2010), 3836--3848.

\bibitem{CGP}  B. Coll, A. Gasull, R. Prohens;
\emph{Degenerate hopf bifurcations in discontinuous planar system}, 
J. Math. Anal. Appl., 253 (2001), 671--690.

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

\bibitem{F}  A. F. Filippov;
\emph{Differential Equation with Discontinuous
Right-Hand Sides}, Kluwer Academic, Amsterdam, 1988.

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

\bibitem{GT}  A. Gasull, J. Torregrosa;
\emph{Center-focus problems for discontinuous planar differential equations}, 
Internat. J. Bifurcation and Chaos, 13 (2003), 1755--1765.

\bibitem{HY}  M. Han, P. Yu;
\emph{Normal Forms, Melnikov Functions and
Bifurcations of Limit cycles}, Springer, NewYork, 2012.

\bibitem{HY2}  S. Huan, X. Yang;
\emph{On the number of limit cycles in general
planar piecewise linear systems}, Discrete Contin. Dyn. Syst. 32
(2012), 2147--2164.

\bibitem{HZ}  M. Han, W. Zhang;
\emph{On hopf bifurcations in nonsmooth planar systems}, 
J. Differential Equations. 248 (2010), 2399--2416.

\bibitem{K}  M. Kukucka;
\emph{Non-Smooth Dynamical Systems}, Springer-Verlag,
Berlin, Heidelberg, 2000.

\bibitem{Li}  J. Li;
\emph{Hilbert's 16th problem and bifurcations of planar
polynomial vector fields}, Int. J. Bifurcation and Chaos, 13 (2003),
47--106.

\bibitem{LZL}  S. Li, Y. Zhao, J. Li;
\emph{On the number of limit cycles of a
perturbed cubic polynomial defferential center}, J. Math. Anal.
Appl., 404 (2013), 212--220.

\bibitem{LZLZ}  W. Li, Y. Zhao, C. Li, Z. Zhang;
\emph{Abelian integrals for quadratic centres having almost
all their orbits formed by quartics}, Nonlinearity, 15 (2002), 863--885.

\bibitem{LHR}  F. Liang, M. Han, R. G. Romanovski;
\emph{Bifurcation of limit cycles by perturbing a piecewise 
linear Hamiltonian system with a homoclinic loop},  
Nonlinear Analysis, 75 (2012), 4355--4374.

\bibitem{LH}  X. Liu, M.Han;
\emph{Bifurcation of limit cycles by perturbing
piecewise Hamiltonian systems}, Internat. J. Bifurcation and Chaos, 20 (2010),
1379--1390.

\bibitem{LM}  J. Llibre, A. C. Mereu;
\emph{Limit cycles for discontinuous quadratic differential systems with 
two zones}, J. Math. Anal. Appl., 413 (2014), 763--775.

\bibitem{LNT}  J. Llibre, D. D. Novaes, M. A. Teixeira;
\emph{Averaging methods for studying the periodic orbits of discontinuous 
differential systems},
arXiv:1205.4211.

\bibitem{LV}  J. Llibre, C. Valls;
\emph{On the limit cycles of polynomial
differential systems with homogeneous nonlinearities of degree 3 via
the averaging method}, Dyn. Contin. Discrete Impuls Syst. Ser. A Math
Anal., 17 (2010), 453--473.

\bibitem{SV}  J. A. Sanders, F. Verhulst;
\emph{Averaging methods in nonlinear dynamic systems}. 
Applied Mathematical Science, vol. 9,
Springer-Verlag, New York, 1985.

\bibitem{X}  Y. Xiong;
\emph{Limit cycles bifurcations by perturbing piecewise
smooth Hamiltonian systems with multiple parameters}, J. Math. Anal.
Appl., 421 (2015), 260--275.

\bibitem{YHH}  J. Yang, M. Han, W. Huang;
\emph{On Hopf bifurcations of piecewise planar Hamiltonian systems}, 
J. Differential Equations. 250 (2011), 1026--1051.

\end{thebibliography}

\end{document}