An Algorithmic Approach to Limit Cycles of Nonlinear Differential Systems: the Averaging Method Revisited

TL;DR

This paper presents three algorithms for analyzing bifurcations of limit cycles in nonlinear differential systems using an improved averaging method, implemented in Maple, and identifies errors in prior literature.

Contribution

It introduces a systematic, algorithmic framework for the averaging method to analyze limit cycles, including transformation, computation, and exact expression derivation.

Findings

Algorithms successfully implemented in Maple

Effective in analyzing bifurcations of limit cycles

Identifies errors in previous publications on averaging method

Abstract

This paper introduces an algorithmic approach to the analysis of bifurcation of limit cycles from the centers of nonlinear continuous differential systems via the averaging method. We develop three algorithms to implement the averaging method. The first algorithm allows to transform the considered differential systems to the normal formal of averaging. Here, we restricted the unperturbed term of the normal form of averaging to be identically zero. The second algorithm is used to derive the computational formulae of the averaged functions at any order. The third algorithm is based on the first two algorithms that determines the exact expressions of the averaged functions for the considered differential systems. The proposed approach is implemented in Maple and its effectiveness is shown by several examples. Moreover, we report some incorrect results in published papers on the averaging method.

Tables

Table 1. Table 1. Number of limit cycles for quadratic isochronous centers

Averaging order	${𝕊}_1$	${𝕊}_2$	${𝕊}_3$	${𝕊}_4$
1,2	0	0	0	0
3,4	1	1	1	1
5	1	2	2	2
6	2	2	2	2
7	2	2	2	2
8	-	-	-	-

Equations

$$\begin{array}{r}\dot{x}=f_n(x,y),\dot{y}=g_n(x,y).\end{array}$$

$$\begin{array}{r}\dot{x}=P(x,y),\dot{y}=Q(x,y).\end{array}$$

$$\begin{array}{rl}\dot{x}& ={\bar{P}}_{\bar{\alpha }}(x,y)=-y+\sum _{m=2}^{n_1}{P}_m(x,y),\\ \dot{y}& ={\bar{Q}}_{\bar{\beta }}(x,y)=x+\sum _{m=2}^{n_1}{Q}_m(x,y),\end{array}$$

$$\begin{array}{rl}\dot{x}& ={\bar{P}}_{\bar{\alpha }}(x,y)+p_{\alpha }(x,y,\epsilon ),\\ \dot{y}& ={\bar{Q}}_{\bar{\beta }}(x,y)+q_{\beta }(x,y,\epsilon )\end{array}$$

$$\begin{array}{r}{p}_{\alpha }(x,y,\epsilon )=\sum _{j=1}^k{\epsilon }^j{\tilde{p}}_j(x,y),q_{\beta }(x,y,\epsilon )=\sum _{j=1}^k{\epsilon }^j{\tilde{q}}_j(x,y),\end{array}$$

$$\begin{array}{r}\frac{dr}{d\theta }=\sum _{i=0}^k{\epsilon }^i{F}_i(\theta ,r)+{\epsilon }^{k+1}R(\theta ,r,\epsilon ),\end{array}$$

$$\begin{array}{r}{f}_i(z)=\frac{y_i(2\pi ,z)}{i!},\end{array}$$

$$\begin{split}y_{i}(\theta,z)&=i!\int_{0}^{\theta}\Bigg{[}F_{i}(s,\varphi(s,z))+\sum_{\ell=1}^{i}\sum_{S_{\ell}}\frac{1}{b_{1}!b_{2}!2!^{b_{2}}\cdots b_{\ell}!\ell!^{b_{\ell}}}\\ &\quad\cdot\partial^{L}F_{i-\ell}(s,\varphi(s,z))\prod_{j=1}^{\ell}y_{j}(s,z)^{b_{j}}\Bigg{]}ds,\end{split}$$

$$\begin{split}y_{1}(\theta,z)&=\int_{0}^{\theta}F_{1}(s,z)ds,\\ y_{i}(\theta,z)&=i!\int_{0}^{\theta}\Bigg{[}F_{i}(s,z)+\sum_{\ell=1}^{i-1}\sum_{S_{\ell}}\frac{1}{b_{1}!b_{2}!2!^{b_{2}}\cdots b_{\ell}!\ell!^{b_{\ell}}}\\ &\quad\cdot\partial^{L}F_{i-\ell}(s,z)\prod_{j=1}^{\ell}y_{j}(s,z)^{b_{j}}\Bigg{]}ds.\end{split}$$

$${\bar{f}}_i(r)=\sum _{j_2=0}^i\sum _{j_1=0}^{N_i}{c}_{j_1,j_2}{r}^{j_1}{\pi }^{j_2}.$$

$${\tilde{c}}_j\in \mathbb{Q}[\pi ,\bar{\alpha },\bar{\beta },\alpha ,\beta ]/{\Sigma }_k$$

$$F_{i,j}=f_{i,j}(r,S)C+g_{i,j}(r,S),$$

$$\begin{array}{r}{B}_{\ell ,m}(x_1,\dots ,x_{\ell -m+1})=\sum _{{\tilde{S}}_{\ell ,m}}\frac{\ell !}{b_1!b_2!\cdots b_{\ell -m+1}!}\prod _{j=1}^{\ell -m+1}{\left(\frac{x_j}{j!}\right)}^{b_j},\end{array}$$

$$\begin{split}y_{1}(\theta,z)&=\int_{0}^{\theta}F_{1}(s,z)ds,\\ y_{i}(\theta,z)&=i!\int_{0}^{\theta}\Bigg{[}F_{i}(s,z)+\sum_{\ell=1}^{i-1}\sum_{m=1}^{\ell}\frac{1}{\ell!}\partial^{m}F_{i-\ell}(s,z)\\ &\quad\cdot B_{\ell,m}(y_{1}(s,z),\ldots,y_{\ell-m+1}(s,z))\Bigg{]}ds.\end{split}$$

$$\begin{array}{r}\dot{x}=-y+x^{2}y,\dot{y}=x+x{y}^2.\end{array}$$

$$\begin{array}{rl}\dot{x}& =-y+x^{2}y+\sum _{s=1}^7{\epsilon }^s{p}_s(x,y),\\ \dot{y}& =x+x{y}^2+\sum _{s=1}^7{\epsilon }^s{q}_s(x,y),\end{array}$$

$$\begin{array}{rl}{p}_s(x,y)& ={\alpha }_{s,1}x+{\alpha }_{s,2}y+{\alpha }_{s,3}{x}^2+{\alpha }_{s,4}xy+{\alpha }_{s,5}{y}^2+{\alpha }_{s,6}{x}^3\\ & +{\alpha }_{s,7}{x}^{2}y+{\alpha }_{s,8}x{y}^2+{\alpha }_{s,9}{y}^3,\\ q_s(x,y)& ={\beta }_{s,1}x+{\beta }_{s,2}y+{\beta }_{s,3}{x}^2+{\beta }_{s,4}xy+{\beta }_{s,5}{y}^2+{\beta }_{s,6}{x}^3\\ & +{\beta }_{s,7}{x}^{2}y+{\beta }_{s,8}x{y}^2+{\beta }_{s,9}{y}^3,\end{array}$$

$$\begin{array}{r}\frac{dr}{d\theta }=\sum _{i=1}^7{\epsilon }^i{F}_i(\theta ,r)+\mathcal{O}({\epsilon }^8).\end{array}$$

$$F_1(\theta ,r)=r({\alpha }_{1,2}+{\beta }_{1,1})SC+r(-{\alpha }_{1,1}+{\beta }_{1,2})S^2+r{\alpha }_{1,1}$$

$$\begin{split}f_{2}(r)&=\frac{\pi r}{2}\big{(}\pi\alpha_{1,1}^{2}+2\pi\alpha_{1,1}\beta_{1,2}+\pi\beta_{1,2}^{2}+\alpha_{1,1}\alpha_{1,2}-\alpha_{1,1}\beta_{1,1}\\ &\quad+\alpha_{1,2}\beta_{1,2}-\beta_{1,1}\beta_{1,2}+2\alpha_{2,1}+2\beta_{2,2}\big{)}.\end{split}$$

$$f_3(r)=\frac{1}{4}\pi r\left(A_2{r}^2+A_0\right),$$

$$A_2=4{\alpha }_{1,1}+3{\alpha }_{1,6}+{\alpha }_{1,8}+{\beta }_{1,7}+3{\beta }_{1,9},A_0=4({\alpha }_{3,1}+{\beta }_{3,2}).$$

$$f_4(r)=\frac{1}{4}\pi r\left(B_2{r}^2+B_0\right),$$

$$\begin{array}{rl}{B}_2& =4{\alpha }_{1,1}{\alpha }_{1,2}+2{\alpha }_{1,1}{\alpha }_{1,7}+2{\alpha }_{1,1}{\beta }_{1,8}+{\alpha }_{1,2}{\alpha }_{1,8}\\ & +3{\alpha }_{1,2}{\beta }_{1,9}+{\alpha }_{1,3}{\alpha }_{1,4}-2{\alpha }_{1,3}{\beta }_{1,3}+{\alpha }_{1,4}{\alpha }_{1,5}\\ & +2{\alpha }_{1,5}{\beta }_{1,5}+{\alpha }_{1,8}{\beta }_{1,1}+3{\beta }_{1,1}{\beta }_{1,9}-{\beta }_{1,3}{\beta }_{1,4}\\ & -{\beta }_{1,4}{\beta }_{1,5}+4{\alpha }_{2,1}+3{\alpha }_{2,6}+{\alpha }_{2,8}+{\beta }_{2,7}+3{\beta }_{2,9},\\ B_0& =4({\alpha }_{4,1}+{\beta }_{4,2}).\end{array}$$

$$f_5(r)=\frac{1}{4}\pi r\left(C_4{r}^4+C_2{r}^2+C_0\right),$$

$$\begin{array}{r}{C}_4=2{\alpha }_{1,1}+2{\alpha }_{1,6}+{\alpha }_{1,8}+{\beta }_{1,9},C_0=4({\alpha }_{5,1}+{\beta }_{5,2}).\end{array}$$

$$f_6(r)=\frac{1}{24}\pi r\left(D_4{r}^4+D_2{r}^2+D_0\right),$$

$$\begin{array}{rl}{D}_4& =12{\alpha }_{1,1}{\alpha }_{1,7}-6{\alpha }_{1,1}{\alpha }_{1,9}-12{\alpha }_{1,1}{\beta }_{1,1}-18{\alpha }_{1,1}{\beta }_{1,6}\\ & -12{\alpha }_{1,2}{\alpha }_{1,6}+7{\alpha }_{1,3}{\alpha }_{1,4}-18{\alpha }_{1,3}{\beta }_{1,3}-20{\alpha }_{1,3}{\beta }_{1,5}\\ & +7{\alpha }_{1,4}{\alpha }_{1,5}-4{\alpha }_{1,5}{\beta }_{1,3}-6{\alpha }_{1,5}{\beta }_{1,5}-18{\alpha }_{1,6}{\alpha }_{1,9}\\ & -12{\alpha }_{1,6}{\beta }_{1,1}-18{\alpha }_{1,6}{\beta }_{1,6}-6{\alpha }_{1,8}{\alpha }_{1,9}-6{\alpha }_{1,8}{\beta }_{1,6}\\ & +{\beta }_{1,3}{\beta }_{1,4}+{\beta }_{1,4}{\beta }_{1,5}+12{\alpha }_{2,1}+12{\alpha }_{2,6}+6{\alpha }_{2,8}+6{\beta }_{2,9},\\ D_0& =24({\alpha }_{6,1}+{\beta }_{6,2}).\end{array}$$

$$f_7(r)=-\frac{1}{48}\pi r\left(E_6{r}^6+E_4{r}^4+E_2{r}^2+E_0\right),$$

$$\begin{array}{r}{E}_6=-3({\alpha }_{1,1}+{\alpha }_{1,6}+{\alpha }_{1,8}),E_0=-48({\alpha }_{7,1}+{\beta }_{7,2}).\end{array}$$

Full text

Algorithmic approach; averaging method; center; limit cycle; nonlinear differential systems

††journalyear: 2019††copyright: acmcopyright††conference: 2019 International Symposium on Symbolic and Algebraic Computation; July 15–18, 2019; Beijing, China††booktitle: 2019 Int’l Symposium on Symbolic $\&$ Algebraic Computation (ISSAC’19), July 15–18, 2019, Beijing, China††price: 15.00††doi: 10.1145/3326229.3326234††isbn: 978-1-4503-6084-5/19/07††ccs: Computing methodologies Symbolic and algebraic manipulation††ccs: Symbolic and algebraic algorithms Symbolic calculus algorithms

1. Introduction

Bounding the number of limit cycles for systems of polynomial differential equations is a long standing problem in the field of dynamical systems. As is well known, the second part of the 16th Hilbert’s problem (Hilbert, 1902; Ilyashenko, 2002) asks about “the maximal number $H(n)$ and relative configurations of limit cycles” for planar polynomial differential systems of degree $n$:

[TABLE]

Solving this problem, even in the case $n=2$, at the present state of knowledge seems to be hopeless. While it has not been possible to obtain uniform upper bounds for $H(n)$ in the near future, there has been success in finding lower bounds. Some known results are as follows: it is shown in (Chen and Wang, 1979; Shi, 1980) that $H(2)\geq 4$ and $H(3)\geq 13$ in (Li et al., 2009). In (Christopher and Lloyd, 1995), it is proved that $H(n)$ grows at least as rapidly as $n^{2}\log n$. For the latest development about $H(n)$, we refer the reader to (Christopher and Li, 2007; Li, 2003).

Recall that a limit cycle of system (1) is an isolated periodic orbit. It is the $\omega$-(forward) or $\alpha$-(backward) limit set of nearby orbits. One classical way of producing limit cycles is by perturbing a differential system which has a center. In this case the perturbed system displays limit cycles that bifurcate, either from the center (having the so-called Hopf bifurcation), or from some of the periodic orbits of the period annulus surrounding the center, see for instance Pontrjagin (Pontrjagin, 1934), the book of Christopher-Li (Christopher and Li, 2007), and the hundreds of references quoted there.

In this paper we study the maximal number of limit cycles that bifurcate from the centers of the unperturbed systems (the so-called small-amplitude limit cycles). The main technique is based on the averaging method. We point out that the method of averaging is a classic and mature tool for studying isolated periodic solutions of nonlinear differential systems in the presence of a small parameter. The method has a long history that started with the classical works of Lagrange and Laplace, who provided an intuitive justification of the method. The first formalization of this theory was done in 1928 by Fatou. Important practical and theoretical contributions to the averaging method were made in the 1930s by Bogoliubov-Krylov, and in 1945 by Bogoliubov. The ideas of averaging method have extended in several directions for finite and infinite dimensional differentiable systems. We refer to the books of Sanders-Verhulst-Murdock (Sanders et al., 2007) and Llibre-Moeckel-Simó (Llibre et al., 2015b) for a modern exposition of this subject.

We remark that most of these previous results developed the averaging method up to first order in a small parameter $\varepsilon$, and at most up to third order. In (Giné et al., 2013; Llibre et al., 2014) the averaging method at any order was developed to study isolated periodic solutions of nonsmooth but continuous differential systems. Recently, the averaging method has also been extended to study isolated periodic solutions of discontinuous differential systems; see (Llibre et al., 2015a; Itikawa et al., 2017). In practice, the evaluation of the averaged functions is a computational problem that require powerful computerized resources. Moreover, the computational complexity grows very fast with the averaging order. In view of this, our objective in this paper is to present an algorithmic approach to develop the averaging method at any order and to further study periodic solutions of nonlinear continuous differential systems.

It is known that the Liapunov constants are a good tool for studying the number of small-amplitude limit cycles which can bifurcate from a singular point, i.e., a Hopf bifurcation. Over the years, a number of algorithms for efficient computation of Liapunov constants have been developed (see (Gasull and Torregrosa, 2001; Wang, 1991, 2004) for instance). But a disadvantage of such an approach is that there is no clear geometry of the bifurcated limit cycles. In contrast, using the expressions of the averaged functions, we can estimate the size of the bifurcated limit cycles as a function of $\varepsilon$ for $|\varepsilon|>0$ sufficiently small, see (Giné and Llibre, 2007; Benterki and Llibre, 2017) for instance.

Overview of Paper. The structure of our paper is as follows. In Section 2, we introduce the basic results on the averaging method for planar differential systems before presenting our main results in Section 3. We give our algorithms and briefly describe their implementation in Maple in Section 4. Its application is illustrated in Section 5 using several examples including a cubic polynomial differential system known as Collins First Form and a class of generalized Kukles polynomial differential systems of degree 6. We end with some discussions in Section 6.

In view of space limitation, we moved the proof of Theorem 3.1 to Appendix A. Two of the examples are found in Appendices B and C. The version with appendices may be found at our website https://cs.nyu.edu/exact/papers/ as well as in the arXiv.

2. Basic Theory of the Averaging Method

In this section we introduce the basic results on the averaging method that we shall use for studying the limit cycles which bifurcate from the centers of polynomial differential systems of degree $n_{1}$ in the form of

[TABLE]

An accessible reference is (Chicone, 2006) (see also (Sanders et al., 2007)). The following definition is due to Poincaré (see (Chavarriga and Sabatini, 1999), Section 2).

Definition 2.1.

We say that an isolated singular point $O$ of (2) is a center if there exists a punctured neighbourhood $V$ of $O$, such that every orbit in $V$ is a cycle surrounding $O$.

Without loss of generality we can assume that the center $O$ of system (2) is the origin of coordinates. In this case, after a linear change of variables and a rescaling of time variable, we can write system (2) in the form

[TABLE]

where $P_{m}$, $Q_{m}$ are homogeneous polynomials of degree $m$ in $x$ and $y$ with $\bar{\alpha}$ and $\bar{\beta}$ are parameters appearing as coefficients of $\bar{P},\bar{Q}$ satisfying that system (3) has a center at the origin. It is well known since Poincaré (Poincaré, 1881) and Liapunov (Liapunov, 1949) that system (3) has a center at the origin if and only if there exists a local analytic first integral of the form $H(x,y)=x^{2}+y^{2}+F(x,y)$ defined in a neighborhood of the origin, where $F$ starts with terms of order higher than 2. For the well known center problem, see (Romanovski and Shafer, 2009; Mahdi et al., 2017).

We now consider the perturbations of (3) of the form

[TABLE]

with

[TABLE]

where the polynomials $\tilde{p}_{j},\tilde{q}_{j}$ are of degree at most $n_{2}$ (usually $n_{2}\geq n_{1}\geq 2$) in $x$ and $y$ with $\alpha$ and $\beta$ are free parameters appearing as coefficients of $\tilde{p}_{j},\tilde{q}_{j}$, and $\varepsilon$ is a small parameter. Note that by “free parameters” we mean that the coefficient of each monomial in $p_{\alpha}$ and $q_{\beta}$ is a distinct parameter in $\alpha$ or $\beta$. We are interested in the maximum number of small-amplitude limit cycles of (4) for $|\varepsilon|>0$ sufficiently small, which bifurcate at $\varepsilon=0$ from the center of (3).

Usually, the averaging method deals with planar differential systems in the following normal form

[TABLE]

where $F_{i}:\mathbb{R}\times D\rightarrow\mathbb{R}$ for $i=0,1,\ldots,k$, and $R:\mathbb{R}\times D\times(-\varepsilon_{0},\varepsilon_{0})\rightarrow\mathbb{R}$ are $\mathcal{C}^{k}$ functions, $2\pi$-periodic in the first variable, being $D$ an open and bounded interval of $(0,\infty)$, and $\varepsilon_{0}$ is a small parameter. As one of the main hypotheses, it is assumed that the solution $\varphi(\theta,z)$ of the unperturbed differential system, $dr/d\theta=F_{0}(\theta,r)$, is $2\pi$-periodic in the variable $\theta$ for every initial condition $\varphi(0,z)=z\in D$.

The averaging method consists in defining a collection of functions $f_{i}:D\rightarrow\mathbb{R}$, called the $i$-th order averaged function, for $i=1,2,\ldots,k$, which control (their simple zeros control), for $\varepsilon$ sufficiently small, the isolated periodic solutions of the differential system (5). In Llibre-Novaes-Teixeira (Llibre et al., 2014) it has been established that

[TABLE]

where $y_{i}:\mathbb{R}\times D\rightarrow\mathbb{R}$, for $i=1,2,\ldots,k$, is defined recursively by the following integral equation

[TABLE]

where $S_{\ell}$ is the set of all $\ell$-tuples of non-negative integers $[b_{1},b_{2},\ldots,b_{\ell}]$ satisfying $b_{1}+2b_{2}+\cdots+\ell b_{\ell}={\ell}$ and $L=b_{1}+b_{2}+\cdots+b_{\ell}$. Here, $\partial^{L}F(\theta,r)$ denotes the Fréchet’s derivative of order $L$ with respect to the variable $r$.

We remark that, in practical terms, the evaluation of the recurrence (7) is a computational problem. Recently in (Novaes, 2017) the Bell polynomials were used to provide a relatively simple alternative formula for the recurrence. In this paper, we will exploit this new formula in our algorithmic approach for solving this problem (see Section 4.2).

Related to the averaging functions (6) there exist two fundamentally different cases in (5), namely, when $F_{0}=0$ and when $F_{0}\neq 0$. We see that when $F_{0}\neq 0$, the formula for $y_{i}(\theta,z)$ in (7) requires the solution of a Cauchy problem because $y_{i}(\theta,z)$ appears on both sides of the equation (see Remark 3 in (Llibre et al., 2014)). The investigation in this paper is restricted to the case where $F_{0}=0$. In this case, we have $\varphi(\theta,z)=z$ for each $\theta\in\mathbb{R}$. Then the integral in equation (7) simplifies to

[TABLE]

The following $k$-th order averaging theorem gives a criterion for the existence of limit cycles. Its proof can be found in Section 2 of (Itikawa et al., 2017).

Theorem 2.2.

(Itikawa et al., 2017)* Assume that $f_{i}\equiv 0$ for $i=1,2,\ldots,j-1$ and $f_{j}\neq 0$ with $j\in\{1,2,\ldots,k\}$. If there exists $\bar{r}\in D$ such that $f_{j}(\bar{r})=0$ and $f^{\prime}_{j}(\bar{r})\neq 0$, then for $|\varepsilon|>0$ sufficiently small, there exists a $2\pi$-periodic solution $r(\theta,\varepsilon)$ of (5) such that $r(0,\varepsilon)\rightarrow\bar{r}$ when $\varepsilon\rightarrow 0$.*

We remark that in order to analyze the Hopf bifurcation for system (4), applying Theorem 2.2, we introduce a small parameter $\varepsilon$ doing the change of coordinates $x=\varepsilon X$, $y=\varepsilon Y$. After that we perform the polar change of coordinates $X=r\cos\theta$, $Y=r\sin\theta$, and by doing a Taylor expansion truncated at $k$-th order in $\varepsilon$ we obtain an expression for $dr/d\theta$ similar to (5) up to $k$-th order in $\varepsilon$. In doing so, the variable $\theta$ appears through sines and cosines, the differential equation in the form $dr/d\theta$ is $2\pi$-periodic. It suffices to take $D=\{r:0<r<r_{0}\}$ with $r_{0}>0$ is arbitrary, since we restrict $F_{0}=0$, the unperturbed system has periodic solutions passing through the points $(0,r)$ with $0<r<r_{0}$.

In general, it is not an easy thing to determine the exact number of simple zeros of the averaged functions (6), since the averaged functions may be too complicated, such as including square root functions, logarithmic functions, and the elliptic integrals. In the literature there is an abundance of papers dealing with zeros of the averaged functions (see for instance (Liang et al., 2016; Novaes and Torregrosa, 2017; Huang, 2019) and references therein). The techniques and arguments to tackle this kind of problem are usually very long and technical.

As a summary of this section, we remark that, using the expressions of the averaged functions, one can estimate the size of bifurcated limit cycles. In fact we know that if the averaged function $f_{j}=0$ for $j=1,\ldots,k-1$ and $f_{k}\neq 0$, and $\bar{r}$ is a simple zero of $f_{k}$, then by Theorem 2.2 there is a limit cycle $r(\theta,\varepsilon)$ of the differential system (5) such that $r(0,\varepsilon)=\bar{r}+\mathcal{O}(\varepsilon)$. Then, going back through the changes of variables we have for the differential system $(\dot{X},\dot{Y})$ the limit cycle $(X(t,\varepsilon),Y(t,\varepsilon))=(\bar{r}\cos\theta,\bar{r}\sin\theta)+\mathcal{O}(\varepsilon)$. Now due to the scaling $x=\varepsilon X,y=\varepsilon Y$ the limit cycles that we find for the differential system (5) coming from our system (4), are in fact limit cycles of the form $(x(t,\varepsilon),y(t,\varepsilon))=\varepsilon(\bar{r}\cos\theta,\bar{r}\sin\theta)+\mathcal{O}(\varepsilon^{2})$ for system (4), which tends to the origin from the origin, i.e., are limit cycles coming by a Hopf bifurcation, for more details on these kind of bifurcations see (Kuznetsov, 2004) for instance.

3. Main Results

Denote the exact upper bound for the number of positive simple zeros of the $i$-th order averaged function $f_{i}(r)$ associated to system (4) by $H_{i}(n_{1},n_{2})$ for $i=1,\ldots,k$. Applying Theorem 2.2, we know that the maximal number of small-amplitude limit cycles of (4) is $H_{i}(n_{1},n_{2})$ and this number can be reached. In this work, we attempt to prove upper bounds on the number of zeros of the $k$-th order averaged function. Our main theorem is the following:

Theorem 3.1.

Assume that $F_{0}=0$ in the normal form (5) associated to the system (4), then there exist a non-negative integer $\nu_{i}\leq i-1$ and a polynomial function $\bar{f}_{i}(r)=\sum_{j=0}^{N_{i}}c_{j}r^{j}$ with $N_{i}\leq in_{2}$, such that $r^{\nu_{i}}f_{i}(r)=\bar{f}_{i}(r)$ for $i=1,\ldots,k$, where the coefficients $c_{j}\in\mathbb{Q}[\pi]$ with degree no more than $i$ in $\pi$.

A detailed proof of it can be found in Appendix A. This result is the first work that deals with the bifurcation of limit cycles of system (4) in the general class of perturbations (see (Llibre and Itikawa, 2015; Llibre and da Silva, 2016) for a few results on some systems of special form). This theorem tells us that the maximum number of small-amplitude limit cycles of (4), which bifurcate from the center of (3) is always finite ( $H_{k}(n_{1},n_{2})\leq N_{k}$). But, for a given system (4), how can we determine the exact value of $N_{i}$ for $i=1,\ldots,k$? In this paper, we provide an algorithmic approach to the solution (see Algorithm 3 in Section 4.2).

Applying Theorems 2.2 and 3.1, we obtain the Theorem 3.2 on $f_{k}(r)$. We first introduce some notations based on Theorems 2.2 and 3.1 before we state this result. Let $R^{*}$ be the real polynomial ring $\mathbb{Q}[\bar{\alpha},\bar{\beta},\alpha,\beta]$. Then for each $\bar{f}_{i}(r)\in R^{*}[\pi][r]$, we define $\mbox{coeffs}(\bar{f}_{i};r,\pi)=\{c_{j_{1},{j_{2}}}:j_{1}=0,\ldots,N_{i};j_{2}=0,\ldots,i\}$, where

[TABLE]

Then $\Sigma_{k}=\cup_{i=1}^{k-1}\mbox{coeffs}(f_{i};r,\pi)=0\subseteq R^{*}$. Now taking the above notations into account and applying Theorems 2.2 and 3.1, we obtain the following theorem on $f_{k}(r)$.

Theorem 3.2.

Assume that $F_{0}=0$ in the normal form (5) associated to the system (4). Then there exist non-negative integers $\tilde{\nu}_{k}\leq k-1$, $\tilde{N}_{k}\leq kn_{2}$ and a polynomial function $\tilde{f}_{k}(r)=\sum_{j=0}^{N_{k}}\tilde{c}_{j}r^{j}$, such that $f_{k}(r)$ in (6) has the form $r^{\tilde{\nu}_{k}}f_{k}(r)=\tilde{f}_{k}(r)$ and the coefficients

[TABLE]

with degree no more than $k$ in $\pi$.

Proof.

The conclusion follows directly from the conditions $\bar{f}_{1}=\bar{f}_{2}=\cdots=\bar{f}_{k-1}=0$. ∎

We remark that, the study of the number of zeros of $f_{k}(r)$ is currently not-algorithmic. Below we give our analysis on this.

Let $\bar{N}=|\bar{\alpha}|+|\bar{\beta}|+|\alpha|+|\beta|$ be the number of parameters in system (4), and $V(\Sigma_{k})\subseteq\mathbb{R}^{\bar{N}}$ is the variety defined by $\Sigma_{k}$. For any point $p^{*}\in V(\Sigma_{k})$, let $\bar{f}_{k}(r;p^{*})\in\mathbb{R}[r]$ be the real polynomial when the parameter are instantiated by $p^{*}$. Finally let $\#(p^{*})$ denote the maximal number of zeros (counted with multiplicity) of $\bar{f}(r;p^{*})$. It follows that $H_{k}(n_{1},n_{2})\leq\mbox{max}\{\#(p^{*}:p\in V(\Sigma_{k}))\}$.

In order to study the number of zeros of function $f_{k}(r)$, according to our Theorem 3.1, it suffices to consider the number of zeros of a polynomial function. Here we provide the Descartes theorem (see (Berezin and Zhidkov, 1965)) to obtain the upper bound of the number of zeros for the polynomial functions.

Lemma 3.3.

(Descartes theorem). Consider the real polynomial $m(x)=a_{s_{1}}x^{s_{1}}+a_{s_{2}}x^{s_{2}}+\cdots+a_{s_{m}}x^{s_{m}}$ with $0=s_{1}<s_{2}<\cdots<s_{m}$ and $a_{s_{j}}\neq 0$ real constants for $j\in\{1,2,\ldots,m\}$. When $a_{s_{j}}a_{s_{j+1}}<0$, we say that $a_{s_{j}}$ and $a_{s_{j+1}}$ have a variation of sign. If the number of variations of signs is $m^{*}$, then $m(x)$ has at most $m^{*}$ positive real roots. Moreover, it is always possible to choose the coefficients of $m(x)$ in such a way that $m(x)$ has exactly $m-1$ positive real roots.

4. Algorithms for the $k$-th Order Averaging Theorem

In this section we will provide an algorithmic approach to revisit the averaging method. According to the averaging method described in Section 2, it is necessary to take the following steps to study the bifurcation of limit cycles for system (4).

STEP 1. Write the perturbed system (4) in the normal form of averaging (5) up to $k$-th order in $\varepsilon$.

STEP 2. (i) Compute the exact formula for the $k$-th order integral function $y_{k}(\theta,z)$ in (8). (ii) Derive the symbolic expression of the $k$-th order averaged function $f_{k}(z)$ by (6).

STEP 3. Determine the exact upper bound for number of positive simple zeros of $f_{k}(z)$.

In the following subsections we will present algorithms to implement the first two steps. We use “Maple-like” pseudo-code, based on our Maple implementation. Using these algorithms we reduce the problem of studying the number of limit cycles of system (4) to the problem of detecting STEP 3.

4.1. Algorithm for STEP 1

In this subsection we will devise an efficient algorithm which can be used to transform system (4) into the form (5). Our algorithm can derive (5) at any order in $\varepsilon$.

Now refer to (4), making the change of variables $x=\varepsilon\cdot r\cdot C$ and $y=\varepsilon\cdot r\cdot S$ with $C=\cos\theta$ and $S=\sin\theta$, we present the algorithm Normalize below based on the above analysis.

The if hypothesis in line 6 is to make sure that $F_{0}=0$. In line 9 the function $\mbox{prem}(a,b,x)$ is the pseudo-remainder of $a$ with respect to $b$ in the variable $x$. The following lemma is obtained directly by the property of the pseudo-remainder.

Lemma 4.1.

The expressions $F_{i,j}$ with $j\in\{1,2\}$ in the algorithm Normalize have the following properties:

[TABLE]

where $f_{i,j}$ and $g_{i,j}$ are polynomials in the variables $r$ and $S$.

4.2. Algorithms for STEP 2

This subsection is devoted to provide effective algorithms to compute the formula and exact expression of the $k$-th order averaged function.

According to (8), we should take the following substeps to compute the $k$-th order averaged function of system (5):

Substep 1. Compute the exact formula for the $k$-th order integral function $y_{k}(\theta,z)$.

Substep 2. Output the symbolic expression for the $k$-th order averaged function $f_{k}(r)$ (not simplified by using $f_{1}\equiv f_{2}\equiv\cdots\equiv f_{k-1}\equiv 0$) for a given differential system (4).

We first recall the partial Bell polynomials which can be used to implement the first substep. For $\ell$ and $m$ positive integers, the Bell polynomials:

[TABLE]

where $\tilde{S}_{\ell,m}$ is the set of all $(\ell-m+1)$-tuples of nonnegative integers $[b_{1},b_{2},\ldots,b_{\ell-m+1}]$ satisfying $b_{1}+2b_{2}+\cdots+(\ell-m+1)b_{\ell-m+1}=\ell$, and $b_{1}+b_{2}+\cdots+b_{\ell-m+1}=m$.

Then the integral in equation (8) reads ((Novaes, 2017), Theorem 2)

[TABLE]

The algorithm Averformula, presented below, is based on (9) that can be used to derive the formula of the $k$-th order integral function $y_{k}(\theta,z)$ (Substep 1).

For the generation of the Bell polynomials (lines 4 and 5) we use the routine IncompleteBellB built-in Maple. We give the outputs $Y_{k}$ of the algorithm for $k=1,2$ (see (32) in Appendix B). Note that the formula of $y_{k}(\theta,z)$ is the first element in the set $Y_{k}$. The second element in $Y_{k}$ (where $F_{k}$ without the dependence on $(s,z)$) can be used to derive an exact expression of $f_{k}$ if we give a concrete differential system (4) (then $F_{k}$ can be assigned to values by the algorithm Normalize), see next algorithm AverFun. We also remark that the formula for the $k$-th order averaged function $f_{k}$ can be obtained directly from $y_{k}$ using (6), so we omit the formula for $f_{k}$ in our algorithm Averformula. We deduce explicitly the formulae of $y_{k}$’s up to $k=5$ (see (33) in Appendix B); one can verify that the outputs of our algorithm are consistent with the results given in (Llibre et al., 2014; Itikawa et al., 2017). In fact our algorithm can compute arbitrarily high order formulae of $y_{k}$’s. In Section 4, we will study a cubic differential system (Collins First Form) and a class of generalized Kukles systems to show the feasibility of our algorithm.

In the last subsection, we provide an algorithm Normalize to transform system (4) into the form of $dr/d\theta$ (normal form of averaging). The algorithm AverFun, presented below, is based on the algorithms Normalize and Averformula, which provides a straightforward calculation method to derive the exact expression of the $k$-th order averaged function for a given differential system in the form (4) (Substep 2).

According to our Theorem 3.1, we know that the output of the algorithm AverFun has the property $f_{k}\in\mathbb{Q}[\pi][r,r^{-1}]$. And the numerator of the expression $f_{k}$ is a polynomial function with degree $N_{k}$. In practice, the calculation of $f_{k}$ typically requires powerful computer resources as the computational complexity grows exponentially with order $k$. It turns out that we can greatly improve the speed by updating the obtained $dr/d\theta$ by using the conditions $f_{1}\equiv f_{2}\equiv\cdots\equiv f_{k-1}\equiv 0$.

We implemented all the algorithms presented in this section in Maple. In the next section, we will apply our general algorithmic approach to analyze the bifurcation of limit cycles for several concrete differential systems in order to show its feasibility.

5. Experiments

In this section, we present the bifurcation of limit cycles for a cubic polynomial differential system as an illustration of our approach explained above. In addition, the bifurcation of limit cycles from the centers of a class of generalized Kukles polynomial differential systems of degree 6 is studied when it is perturbed inside the class of all polynomial differential systems of the same degree, and as an application of our method, we also report some results on quadratic differential systems with isochronous centers. The obtained results of our experiments show the feasibility of our approach.

5.1. Illustrative Example

In this subsection, we consider a cubic center of the following polynomial system

[TABLE]

This system is known as Collins First Form, see (Llibre and Itikawa, 2015) for more details.

More concretely, we consider the perturbations of (10) in the form of

[TABLE]

where

[TABLE]

being $\alpha_{s,j}$ and $\beta_{s,j}$, for $s=1,\ldots,7$ and $j=1,\ldots,9$, real constants.

Next, we use our algorithms to study the maximum number of limit cycles of (11) that bifurcate from the center of (10). Applying our algorithm Normalize by taking $k=7$ we obtain

[TABLE]

Here we give only the expression of $F_{1}(\theta,r)$, the explicit expressions of $F_{i}(\theta,r)$ for $i=2,\ldots,7$ are quite large so we omit them.

[TABLE]

with $C=\cos\theta$ and $S=\sin\theta$.

Using our algorithm AverFun in Section 4 and computing $f_{1}$ we obtain $f_{1}(r)=\pi r(\alpha_{1,1}+\beta_{1,2})$. Clearly equation $f_{1}(r)$ has no positive zeros. Thus the first averaged function does not provide any information about the limit cycles that bifurcate from the center of (10) when we perturb it.

Computing $f_{2}$ we obtain

[TABLE]

According to our Theorem 3.1, we take $\bar{f}_{2}(r)=f_{2}(r)$ with degree $N_{2}=1$, and $c_{1}$ is a polynomial in $\pi$ with degree 2. Note that $f_{1}(r)=0$ means that $\beta_{1,2}=-\alpha_{1,1}$. Using this condition we can simplify $f_{2}(r)$ into the form $f_{2}(r)=\pi r(\alpha_{2,1}+\beta_{2,2})$. As for the first averaged function, the second one also does not provide information on the bifurcating limit cycles. From now on, for each $k=3,\ldots,7$, we will perform the calculation of the averaged function $f_{k}$ under the hypothesis $f_{j}\equiv 0$ for $j=1,\ldots,k-1$.

Doing $\beta_{2,2}=-\alpha_{2,1}$ and computing $f_{3}$ we obtain

[TABLE]

where

[TABLE]

Therefore $f_{3}(r)$ can have at most one positive real root. From Theorem 2.2 it follows that the 3-th order averaging provides the existence of at most one small-amplitude limit cycle of system (11) and this number can be reached by Lemma 3.3, since $A_{i}$ for $i=0,2$ are independent constants ($\partial(A_{2},A_{0})/\partial(\beta_{1,7},\beta_{3,2})=4\neq 0$).

To consider the 4-th order averaging theorem we take $\beta_{1,7}=-A_{2}+\beta_{1,7}$ and $\beta_{3,2}=-A_{0}/4+\beta_{3,2}$. Computing $f_{4}$ we obtain

[TABLE]

where

[TABLE]

It is obvious that $f_{4}(r)$ can have at most one positive real root. From Theorem 2.2 it follows that the 4-th order averaging provides the existence of at most one small-amplitude limit cycle of system (11) and this number can be reached ($B_{2}$ and $B_{0}$ are independent constants).

Letting $\beta_{2,7}=-B_{2}+\beta_{2,7}$ and $\beta_{4,2}=-B_{0}/4+\beta_{4,2}$ we obtain $f_{4}(r)=0$. Computing $f_{5}$ we obtain

[TABLE]

where

[TABLE]

We do not explicitly provide the expression of $C_{2}$, because it is very long. It is not hard to check that $C_{4}$, $C_{2}$ and $C_{0}$ are independent constants. Therefore $f_{5}(r)$ can have at most two positive real roots. Then the 5-th order averaging provides the existence of at most two small-amplitude limit cycle of system (11) and this number can be reached.

To consider the 6-th order averaging theorem we let $\beta_{1,9}=-C_{4}+\beta_{1,9}$, $\beta_{3,7}=-C_{2}+\beta_{3,7}$ and $\beta_{5,2}=-C_{0}/4+\beta_{5,2}$. Computing $f_{6}$ we obtain

[TABLE]

where

[TABLE]

Here we do not provide the explicit expression of $D_{2}$ because it is quite long. Moreover $D_{4}$, $D_{2}$ and $D_{0}$ are independent constants. In fact only $D_{4}$ presents the parameter $\alpha_{2,6}$, only $D_{2}$ has the parameter $\alpha_{2,2}$, and $D_{0}$ is the only one with parameters $\alpha_{6,1}$ and $\beta_{6,2}$. Hence $f_{6}(r)$ has at most two positive simple roots. Then the 6-th order averaging provides the existence of at most two small-amplitude limit cycle of system (11) and this number can be reached.

To consider the 7-th order averaging theorem we take $\beta_{2,9}=-D_{4}/6+\beta_{2,9}$, $\beta_{4,7}=-D_{2}/6+\beta_{4,7}$ and $\beta_{6,2}=-D_{0}/24+\beta_{6,2}$. Computing $f_{7}$ we obtain

[TABLE]

where

[TABLE]

Here we do not provide the explicit expressions of $E_{2}$ and $E_{4}$ because they are quite long. Moreover $E_{j}$ for $j=0,2,4,6$ are independent constants. Hence $f_{7}(r)$ has at most three positive simple roots. Then the 7-th order averaging provides the existence of at most three small-amplitude limit cycle of system (11) and this number can be reached.

We remark that our averaged functions $f_{j}(r)$ for $j=1,\ldots,5$ are consistent with the forms in (Llibre and Itikawa, 2015). However, our averaged function $f_{6}(r)$ looks much simpler than the form given in (Llibre and Itikawa, 2015), this is because we rigorously simplify the function $f_{6}(r)$ under the conditions $f_{1}\equiv f_{2}\equiv\cdots\equiv f_{5}\equiv 0$. The averaged function $f_{6}(r)$ in (Llibre and Itikawa, 2015) is not correct, because the authors do not simplify this expression in a right way (in fact one should note that the isolated parameter $\beta_{3,7}$ contains the parameter $\beta_{1,9}$). As a consequence, the maximum number $\{3\}$ of limit cycles of system (10) up to the 6-th order averaging they obtained can not be reached. Thus, some calculations of the averaged functions in (Llibre and Itikawa, 2015) need to be reconsidered algorithmically, using the algorithm and exact formula of the averaged function in this paper.

Here we restate the result related to the Collins First Form as follows.

Theorem 5.1.

For $|\varepsilon|>0$ sufficiently small the maximum number of small-amplitude limit cycles of system (11) is 3 using the 7-th order averaging method, and this number can be reached.

5.2. A Class of Generalized Kukles Differential Systems

In this subsection we consider the perturbations

[TABLE]

of system $\eqref{eq4.2.1}_{\varepsilon=0}$, where $a_{s,j,i}$ and $b_{s,j,i}$ are real parameters, for $s=1,\ldots,6$, $0\leq i\leq j\leq 6$, and $a,b,c$ are real coefficients satisfying $a^{2}+b^{2}+c^{2}\neq 0$. We note that the bifurcation of limit cycles of $\eqref{eq4.2.1}$ has been studied in (Llibre and da Silva, 2016) up to 6-th order averaging theorem ((Llibre and da Silva, 2016), Section 7.3). Here restudy it by using our algorithmic approach to illustrate its feasibility.

We remark that taking $k=6$ our algorithm Normalize can not pass the if hypothesis, this is because the unperturbed term (i.e., the constant term of $H$ in the algorithm Normalize)

[TABLE]

does not vanish. So we have to exclude the perturbed terms $\varepsilon a_{1,0,0}$ and $\varepsilon b_{1,0,0}$ in $\eqref{eq4.2.1}$. However, the authors in (Llibre and da Silva, 2016) obtained a wrong expression of $F_{0}$ in the form

[TABLE]

In fact one can easily check this mistake by manual calculation. So the calculations of the averaged functions of system $\eqref{eq4.2.1}$ in (Llibre and da Silva, 2016) must be redone.

Now consider system $\eqref{eq4.2.1}$, letting $a_{1,0,0}=b_{1,0,0}=0$ and using our algorithm AverFun in Section 4 we obtain the averaged functions up to 6-th order as follows. Since the calculations and arguments are quite similar to those used in the previous subsection we do not explicitly present the process here.

[TABLE]

The expressions of $E_{i},G_{i}$ for $i=0,2$ and $H_{j},I_{j}$ for $j=0,2,4$ are quite long so we omit them for brevity.

In view of these expressions in (14), we verified that (Theorem 3 in (Llibre and da Silva, 2016)) the averaging theorem up to sixth order provides the existence of at most two small-amplitude limit cycles of system (13).

5.3. Quadratic Systems

In order to save space, we put the results in Appendix C.

6. Discussions

In this paper we present a systematical approach to study the maximum number of limit cycles of differential system (4) for $|\varepsilon|>0$ sufficiently small, which bifurcate from the centers of differential systems in the form of (3). In general, we give three algorithms to analyze the averaging method. Then with the aid of these algorithms, we reduce the study of the number of limit cycles of system (4) to the problem of estimating the number of simple zeros of the obtained averaged functions. Theoretically, we show that the maximum number of limit cycles of system (4) has no more than $kn_{2}$ (a rough bound) by using the $k$-th order averaging method. We believe that the first averaged function $f_{k}$ which is not identically zero is a polynomial in $r$ with odd terms. However, we cannot prove this, we leave this as a future research problem.

We remark that, though in the present paper, we focus our attention on the study of bifurcation of limit cycles of the continuous differential system (4), the developed algorithmic approach admits a generalization to the case of studying the bifurcation of limit cycles for discontinuous differential systems. It is of great interest to employ our approach to analyze the bifurcation of limit cycles for differential systems in many different fields (biology, chemistry, economics, engineering, mathematics, physics, etc.). It will be beneficial to generalize our approach to the case of higher dimension differential systems by using the general form of the averaging method. We leave this as the future research problems. Furthermore, how to simplify and optimize the steps of the computations of the averaged functions is also worthy of further study.

Acknowledgements.

Huang’s work is partially supported by China Scholarship Council under Grant No.: 201806020128. Yap’s work is partially supported by NSF Grants #CCF-1423228 and #CCF-1564132, and also a Chinese Academy of Science (Beijing) President’s International Fellowship Initiative (2018), and Beihang International Visiting Professor Program No. Z2018060. The first author is grateful to Professor Dongming Wang for his encouragement and helpful suggestions, and to Chee Yap for inviting him to visit NYU Courant. Both authors thank the anonymous referees for their valuable comments on improving the presentation.

Appendix A Proof of Theorem 3.1

We first give some lemmas before we prove the Theorem 3.1. The following lemma plays a key role in determining the numbers $\nu_{i}$ and $N_{i}$.

Lemma A.1.

If $k\geq 2$ and $p\in\mathbb{N}$, then for any polynomial $g_{n}(x)$ of degree $n$,

[TABLE]

where $\bar{g}_{n}(x)$ is a polynomial of degree no more than $n$. Here $g^{(p)}$ denotes the $p$-order derivative of a function $g$.

Proof.

The lemma follows directly from the following equality

[TABLE]

∎

The lemma described below can be used to determine the expression form of the averaged function $f_{i}(r)$ in Theorem 3.1.

Lemma A.2.

Define the integral function

[TABLE]

Then we have the following recursive formula for $M_{i,j,k}$

[TABLE]

Moreover, when $k=0$, we have

[TABLE]

Proof.

Doing integration by parts for (15), we have

[TABLE]

Then we find the recursive integral formula (16).

When $k=0$, doing integration by parts for $M_{i,j,0}$, we obtain

[TABLE]

On the other hand, by doing integration by parts for $M_{i-1,j-1,1}$, in a similar way we have

[TABLE]

Using equations (18) and (19), we obtain (17).

∎

Proof of Theorem 3.1. Now refer to system (4), we define the perturbed terms

[TABLE]

with $p_{t}^{j}$, $q_{t}^{j}$ homogeneous polynomials of degree $t$. The change of coordinates

[TABLE]

carries system (4) into

[TABLE]

In polar coordinates $X=rC$ and $Y=rS$ with $C=\cos\theta$, $S=\sin\theta$, system (20) has the form

[TABLE]

Then

[TABLE]

where

[TABLE]

Computing the first-order Taylor expansion of $dr/d\theta$ in $\varepsilon$ we obtain

[TABLE]

Since we assume that $F_{0}=0$, we need to let $p_{0}^{1}=q_{0}^{1}=0$. Then the resulting expression of $dr/d\theta$ is of the form

[TABLE]

where

[TABLE]

and the expressions of $A_{i}$ and $B_{i}$ for $i=2,\ldots,n_{2}+k-1$ are summation of a kind of polynomial functions in the form $r^{i_{1}}\bar{H}_{i_{1}}(C,S)$ with $i_{1}$ non-negative integer and $\bar{H}_{i_{1}}$ polynomial function in the variables $C$ and $S$. Moreover, by observing (21) we know that $B_{i}$ for $i=2,\ldots,k-1$ is a polynomial in $r$ of degree at most $n_{2}+1$ without constant term; and $A_{i}$ is a polynomial in $r$ of degree at most $n_{2}$ in the form:

[TABLE]

where $\bar{A}_{i,1}(r,C,S)$ is a polynomial in $r$ of degree at most $n_{2}$ without constant term, and

[TABLE]

We recall that, given any real value $|\eta|<1$, the following expansion holds:

[TABLE]

Thus, equation (22) can be written as

[TABLE]

where

[TABLE]

and the expressions of $F_{i}$ for $i=4,\ldots,k$ are linear combination of a kind of functions in the form $r^{\alpha_{1}}A_{j_{1}}^{\alpha_{2}}B_{j_{2}}$ with $-i\leq\alpha_{1}\leq-1$, $1\leq j_{1},\alpha_{2}\leq i-1$, and $1\leq j_{2}\leq i$ (here we have avoided the dependence on $(r,C,S)$ to simplify the notation). Recalling the property that $B_{i}$ is a polynomial in $r$ of degree at most $n_{2}+1$ without constant term and $A_{i}$ is a polynomial in $r$ of degree at most $n_{2}$ with constant term, we find that $F_{1}$ is a polynomial in $r$ of degree at most $n_{2}\geq 2$ and $F_{i}$ is a rational function in $r$ of the form

[TABLE]

where $\bar{F}_{i}(r,C,S)$ is a polynomial in $r$ of degree at most $in_{2}$.

In what follows, we first prove that there exist a non-negative integer $\nu_{i}$ and a polynomial function $\bar{f}_{i}(r)=\sum_{j=0}^{N_{i}}c_{j}r^{j}$, such that $r^{\nu_{i}}f_{i}(r)=\bar{f}_{i}(r)$ for $i=1,\ldots,k$, then we provide the bounds for the numbers $\nu_{i}$ and $N_{i}$.

Let $R^{SC}=\{r^{\lambda_{1}}\sin^{\lambda_{2}}\theta\cos^{\lambda_{3}}\theta:\lambda_{1}\in\mathbb{Z},\lambda_{2},\lambda_{3}\in\mathbb{N}\}$ be a set of functions. It is obvious that each $F_{i}$ in (25) (or (26)) is a function generated by linear combination of elements of $R^{SC}$. Note that the explicit expression of $F_{1}$ is of the form

[TABLE]

Now refer to (8), it is easy to check that $y_{1}(\theta,r)$ is a function generated by linear combination of elements of the set of functions in the form $\{\theta r,r^{j_{1}}\sin^{j_{2}}\theta\cos^{j_{3}}\theta\}$ with $0\leq j_{1}\leq 2$ and $0\leq j_{2},j_{3}\leq 3$.

Let $\bar{R}=R^{SC}\times\Theta=\{\theta^{\lambda_{0}}r^{\lambda_{1}}\sin^{\lambda_{2}}\theta\cos^{\lambda_{3}}\theta\}$ be a set of functions with $\lambda_{0}\in\mathbb{N}_{+}$. We denote by $\mbox{Span}(\bar{R})$ be the set of functions generated by linear combination of elements of $\bar{R}$. Next, we will show that the integral function $y_{i}(\theta,r)\in\mbox{Span}(\bar{R})$ for $i=2,\ldots,k$.

First, it is critical to observe that, the resulting form of $\partial^{L}F_{i}(\theta,r)$ is a function generated by linear combination of elements of $R^{SC}$. Since $y_{1}(\theta,r)$ contains $\theta$, the function in the square bracket of (8) is in $\mbox{Span}(\bar{R})$. In order to prove $y_{i}(\theta,r)\in\mbox{Span}(\bar{R})$, we need to consider the following integral equation:

[TABLE]

Second, we claim that $M_{i,j,k}\in\mbox{Span}(\bar{R})$. It follows from Lemma A.2 that $M_{i,j,k}\in\mbox{Span}(\bar{R})$ if and only if $M_{i-1,j+1,k-1}\in\mbox{Span}(\bar{R})$ and $M_{i,j+2,k-2}\in\mbox{Span}(\bar{R})$.

Reuse the recursive formula (16) until the subscript $i=0$ or $k=0$. In this way it suffices to consider $M_{0,j,k}\in\mbox{Span}(\bar{R})$ and $M_{i,j,0}\in\mbox{Span}(\bar{R})$. It is easy to judge that $M_{0,j,k}\in\mbox{Span}(\bar{R})$, so we focus on the proof of $M_{i,j,0}\in\mbox{Span}(\bar{R})$. By using Lemma A.2 and reusing the recursive formula (17), we conclude that $M_{i,j,0}\in\mbox{Span}(\bar{R})$. Then $M_{i,j,k}\in\mbox{Span}(\bar{R})$. Thus the desired result $y_{i}(\theta,r)\in\mbox{Span}(\bar{R})$ holds.

Finally, letting $\theta=2\pi$ in $y_{i}(\theta,r)$ (equation (6)) and taking into accounting the following formulae

[TABLE]

we prove that there exist a non-negative integer $\nu_{i}$ and a polynomial function $\bar{f}_{i}(r)=\sum_{j=0}^{N_{i}}c_{j}r^{j}$, such that $r^{\nu_{i}}f_{i}(r)=\bar{f}_{i}(r)$ for $i=1,\ldots,k$.

Next we provide the bounds for the numbers $\mu_{i}$ and $N_{i}$.

Case $i=1$, since $F_{1}$ is a polynomial in $r$ of degree 2 (at most $n_{2}$), we have by equation (8) that $y_{1}$ is a polynomial in $r$ of degree 2 (at most $n_{2}$).

We assume, by the induction hypothesis, that $y_{i}$ is a rational function in $r$ of the form

[TABLE]

where $\bar{y}_{i-1}(r,C,S)$ is a polynomial in $r$ of degree at most $(i-1)n_{2}$.

In the expression of $y_{i}$ given in (8), there only appear the previous functions $y_{j}$, for $1\leq j\leq i-1$. Now by using equation (26) and Lemma A.1, for a given integer $\ell$ with $1\leq\ell\leq i-1$, we have the following summation function

[TABLE]

We have used the equalities $L=b_{1}+b_{2}+\cdots+b_{\ell}$ and $b_{1}+2b_{2}+\cdots+\ell b_{\ell}={\ell}$ to simplify (28). Combining equations (26) and (27), we know that the numerator of the expression (28) is a polynomial in $r$ with degree at most

[TABLE]

Thus $y_{i}$ is a rational function in $r$ of the form

[TABLE]

where $\bar{y}_{i}(r,C,S)$ is a polynomial in $r$ of degree at most $in_{2}$.

Herewith, we prove that there exist a non-negative integer $\nu_{i}\leq i-1$ and a polynomial function $\bar{f}_{i}(r)=\sum_{j=0}^{N_{i}}c_{j}r^{j}$ with $N_{i}\leq in_{2}$, such that $r^{\nu_{i}}f_{i}(r)=\bar{f}_{i}(r)$ for $i=1,\ldots,k$.

Next we will show that the coefficients $c_{j}$ of $\bar{f}_{i}(r)$ is a polynomial in $\pi$ of degree at most $i$. To do this, we just consider the dependence on $\theta$ for brevity.

We define $\Delta Y_{i}=\{\theta^{\Delta_{i}}\sin^{k_{1}}\theta\cos^{k_{2}}\theta:0\leq\Delta_{i}\leq i,k_{1},k_{2}\in\mathbb{N}\}$ be a set of functions. We claim that the following property holds

[TABLE]

We begin to prove this by induction.

Case $i=1$, we recall that $y_{1}(\theta)$ is a function generated by linear combination of elements of the set of functions $\{\theta,\sin^{j_{2}}\theta\cos^{j_{3}}\theta\}$ with $0\leq j_{2},j_{3}\leq 3$. It obvious that $y_{1}(2\pi)$ is a polynomial in $\pi$ of degree at most $1$.

Suppose that for $\bar{k}\leq i-1$, property (30) holds, then for $\bar{k}=i$, using the integral equation (8), for a given integer $\ell$ with $1\leq\ell\leq i-1$, we have

[TABLE]

Note that by the induction hypothesis, we have $y_{j}(\theta)\in\mbox{Span}(\Delta Y_{j})$ for $1\leq j\leq\ell\leq i-1$. Then the degree of $\theta$ in (31) is at most $b_{1}+2b_{2}+\cdots+\ell b_{\ell}=\ell\leq i-1$.

By using Lemma A.2 (the integral formulae (16) and (17)) and noting also that $\int_{0}^{\theta}\theta^{i-1}d\theta=\theta^{i}/i$, we find that $y_{i}(\theta)\in\mbox{Span}(\Delta Y_{i})$. Finally, letting $\theta=2\pi$, we prove that $y_{i}(2\pi)$ is a polynomial in $\pi$ of degree at most $i$. That is to say, the coefficients $c_{j}$ of $\bar{f}_{i}(r)$ is a polynomial in $\pi$ of degree at most $i$.

Up to now, we finish the proof of Theorem 3.1.

Appendix B Fifth Order Averaging Formulae

We present some formulas computed by Averformula (Algorithm 2).

[TABLE]

[TABLE]

where

[TABLE]

Appendix C Quadratic Systems

This appendix is an overflow from Subsection 5.3. In this subsection we report some results on quadratic differential systems with centers of the form

[TABLE]

C.1. Isochronous Quadratic Centers

We recall that the classification of such quadratic system having an isochronous center at the origin is due to Loud (Loud, 1964). He proved that after an affine change of variables and a rescaling of time any quadratic isochronous center can be written as one of the following four systems.

[TABLE]

In the case of limit cycles bifurcating from the periodic orbits surrounding such quadratic isochronous centers, Chicone and Jacobs in (Chicone and Jacobs, 1991) proved that, under all quadratic polynomial perturbations, at most 1 limit cycle bifurcate from the periodic orbits of $\mathbb{S}_{1}$, and at most 2 limit cycles bifurcate from the periodic orbits of $\mathbb{S}_{2}$, $\mathbb{S}_{3}$ and $\mathbb{S}_{4}$. Iliev obtained in (Iliev, 1998) that the cyclicity of the period annulus surrounding the center $\mathbb{S}_{1}$ is also 2.

Here we focus on the study of limit cycles that bifurcate from such quadratic isochronous centers, and the perturbation terms in (4) are taken as follows:

[TABLE]

Since the calculations and arguments are quite similar to those used in the previous proofs, we just summarize our results in the following Table 1.

We remark that the computation of the 8-th order averaged functions would be too demanding (Maple was consuming too much of the CPU during a calculation). Since we are providing lower bounds for the maximum number of limit cycles that bifurcate from the origin of such quadratic systems, the results could be improved using higher orders of the averaging theorem. Thus, we have a conjecture that some of the numbers 2 obtained in Table 1 may could be increased to 3 as Bautin (Bautin, 1954) proved that in a sufficiently small neighborhood $\Omega$ of a quadratic center, all sufficiently small quadratic perturbations of the given system have at most three limit cycles in $\Omega$, and that three arbitrarily small-amplitude limit cycles can be produced.

C.2. Reversible System with Two Centers

Next, we study the following reversible quadratic system

[TABLE]

with two centers $(0,0)$ and $(1,0)$, where $a_{1}$ and $a_{4}$ are real coefficients satisfying $a_{1}<-1$ ((Yu and Han, 2012), Theorem 1). The authors in (Yu and Han, 2012) proved that 3 limit cycles can bifurcate from the center $(0,0)$ under the case $a_{4}=(a_{1}-5)/3$ based on the Melnikov function method by adding perturbed terms $p(x,y,\varepsilon)=\varepsilon a_{10}x$ and $q(x,y,\varepsilon)=\varepsilon(b_{01}y+b_{11}xy)$ (see Section 3 of (Yu and Han, 2012) for more details). Here using the averaging method we study system (35) by choosing a similar kind of perturbed terms, and then we give some remarks on the relations between these two methods.

First, introducing $x=-\bar{x}$, $y=\bar{y}$ into (35) results in

[TABLE]

which is similar to system (3), now has centers $(0,0)$ and $(-1,0)$. We then consider the perturbations

[TABLE]

of (36), where $c_{s,i,j}$ and $d_{s,i,j}$ are real parameters for $s=1,\ldots,10$.

Computing the averaged functions under the case $a_{4}+1\neq 0$, we obtain the expressions of $f_{k}$’s up to $k=5$ as follows:

[TABLE]

and for $k=6,\ldots,10$, we have

[TABLE]

where

[TABLE]

In view of these expressions of the obtained averaged functions, we find that the $k$-th ($k=5,\ldots,10$) order averaging provides the existence of at most two small-amplitude limit cycles of the perturbed system (37) and this number can be reached under the condition $(a_{1}-a_{4})(a_{1}+2a_{4})(a_{1}-3a_{4}-5)\neq 0$. We conjecture that the maximal number of small-amplitude limit cycles of the perturbed system (37) is 2 up to the $k$-th order averaging for any $k\geq 6$. This problem might be proved by using the recursive integral equation (8).

Our result on the quadratic system (35) describes the different mechanisms between the averaging method and the Melnikov function method when studying the number of limit cycles that can appear in a Hopf bifurcation from centers. The number of limit cycles obtained by the averaging method in some cases (under a similar kind of perturbations) seems to be less than the number obtained by the Melnikov function method. We want to say that in the study of the limit cycles which bifurcate from a period annulus surrounding the center, the equivalence between the averaging method and the Melnikov function method at any order has been proved in (Buicǎ, 2017; Han et al., 2016).