Algorithmic Averaging for Studying Periodic Orbits of Planar Differential Systems
aa r X i v : . [ c s . S C ] M a y Algorithmic Averaging for Studying Periodic Orbits of PlanarDifferential Systems
Bo Huang
LMIB-School of Mathematical Sciences, Beihang UniversityBeijing, ChinaCourant Institute of Mathematical Sciences, New York UniversityNew York, [email protected]
ABSTRACT
One of the main open problems in the qualitative theory of real pla-nar differential systems is the study of limit cycles. In this article,we present an algorithmic approach for detecting how many limitcycles can bifurcate from the periodic orbits of a given polynomialdifferential center when it is perturbed inside a class of polynomialdifferential systems via the averaging method. We propose foursymbolic algorithms to implement the averaging method. The firstalgorithm is based on the change of polar coordinates that allowsone to transform a considered differential system to the normalform of averaging. The second algorithm is used to derive the solu-tions of certain differential systems associated to the unperturbedterm of the normal of averaging. The third algorithm exploits thepartial Bell polynomials and allows one to compute the integral for-mula of the averaged functions at any order. The last algorithm isbased on the aforementioned algorithms and determines the exactexpressions of the averaged functions for the considered differen-tial systems. The implementation of our algorithms is discussedand evaluated using several examples. The experimental resultshave extended the existing relevant results for certain classes ofdifferential systems.
CCS CONCEPTS • Computing methodologies → Symbolic and algebraic manip-ulation; •
Symbolic and algebraic algorithms → Symbolic cal-culus algorithms.
KEYWORDS
Algorithmic approach; averaging method; limit cycles; planar dif-ferential systems; periodic orbits
ACM Reference Format:
Bo Huang. 2020. Algorithmic Averaging for Studying Periodic Orbits of Pla-nar Differential Systems. In & AlgebraicComputation (ISSAC’20), July 20–23, 2020, Kalamata, Greece.
ACM, NewYork, NY, USA, 13 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn
Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full cita-tion on the first page. Copyrights for components of this work owned by others thanACM must be honored. Abstracting with credit is permitted. To copy otherwise, or re-publish, to post on servers or to redistribute to lists, requires prior specific permissionand/or a fee. Request permissions from [email protected].
ISSAC ’20, July 20–23, 2020, Kalamata, Greece © 2020 Association for Computing Machinery.ACM ISBN 978-x-xxxx-xxxx-x/YY/MM...$15.00https://doi.org/10.1145/nnnnnnn.nnnnnnn
We deal with polynomial differential systems in R of the form dxdt = Û x = f n ( x , y ) , dydt = Û y = д n ( x , y ) , (1.1)where n is the maximum degree of the polynomials f and д . Aswe knew, the second part of the 16th Hilbert’s problem [19, 25]asks for “the maximal number H ( n ) and relative configurations oflimit cycles” for the differential system (1.1). Here H ( n ) is called theHilbert number. The problem is still open even for n =
2. However,there have been many interesting results on the lower bound of H ( n ) for n ≥
2: it is shown in [5, 53] that H ( ) ≥ H ( ) ≥ H ( n ) grows at least as rapidly as n log n . For the latest development about H ( n ) , we refer the readerto [6, 16, 32].We recall that a limit cycle of the differential system (1.1) is anisolated periodic orbit of the system. One of the best ways of pro-ducing limit cycles is by perturbing a differential system which hasa center. In this case the perturbed system displays limit cycles thatbifurcate, either from the center (having the so-called Hopf bifur-cation), or from some of the periodic orbits surrounding the center,see the book of Christopher-Li [6] and the references cited therein.Usually, a limit cycle which bifurcates from a center equilibriumpoint is called a small amplitude limit cycle , and a medium ampli-tude limit cycle is one which bifurcates from a periodic orbit sur-rounding a center (see [35, 38]). Note that the notation of “large”limit cycle may occur in several situations in the literature, see [18,58]. In the past seven decades, many researchers have consideredthe small amplitude limit cycles and obtained many results (e.g.,[1, 30, 43, 54, 57]). Over the years, a number of algebraic methodsand algorithms have been developed (e.g., [13, 17, 44, 50, 55, 56])based on the tools of Liapunov constants or Melnikov function.In our recent work [24], we provide an algorithmic approach tosmall amplitude limit cycles of nonlinear differential systems bythe averaging method, and give an upper bound of the number ofzeros of the averaged functions for the general class of perturbeddifferential systems ([24], Thm. 3.1). The goal of this paper is toextend our algorithmic approach to study the maximal numberof medium amplitude limit cycles that bifurcate from some peri-odic orbits surrounding the centers of the unperturbed systems.The main technique is based on the general form of the averagingmethod for planar differential systems.The method of averaging is an important tool to study the ex-istence of isolated periodic solutions of nonlinear differential sys-tems in the presence of a small parameter. It can be used to find aower bound of the Hilbert number H ( n ) for certain differential sys-tems. The method has a long history that started with the classicalworks of Lagrange and Laplace, who provided an intuitive justifica-tion of the method. The first formalization of this theory was donein 1928 by Fatou. Important practical and theoretical contributionsto the averaging method were made in the 1930’s by Bogoliubov-Krylov, and in 1945 by Bogoliubov. The ideas of averaging methodhave extended in several directions for finite and infinite dimen-sional differentiable systems. For a modern exposition of this sub-ject, see the books of Sanders-Verhulst-Murdock [52] and Llibre-Moeckel-Simó [39].The averaging method provides a straightforward calculationapproach to determine the number of limit cycles that bifurcatefrom some periodic orbits of the regarded particular class of dif-ferential systems. However, in practice, the evaluation of the aver-aged functions is a computational problem that requires powerfulcomputerized resources. Moreover, the computational complexitygrows very fast with the averaging order. Our objective in this pa-per is to present an algorithmic approach to develop the averagingmethod at any order and to further study the number of mediumamplitude limit cycles for nonlinear differential systems.In general, to obtain analytically periodic solutions of a differ-ential system is a very difficult problem, many times a problemimpossible to solve. As we shall see when we can apply the aver-aging method, this difficult problem is reduced to finding the zerosof a nonlinear function in an open interval of R , i.e., now the prob-lem has the same difficulty as the problem of finding the singularor equilibrium points of a differential system.The structure of our paper is as follows. In Section 2, we intro-duce the basic results on the averaging method for planar differ-ential systems. We give our algorithms and briefly describe theirimplementation in Maple in Section 3. Its application is illustratedin Section 4 using several examples including a class of generalizedKukles polynomial differential systems and certain differential sys-tems with uniform isochronous centers of degrees 3 and 4. Finally,a conclusion is provided in Section 5. The Maple code of the algo-rithms can be download from https://github.com/Bo-Math/limit-cycle.In view of space limitation, we present the explicit formulae ofthe k -th order averaged function up to k = In this section we introduce the basic theory of the averaging method.We consider the following polynomial differential system of degree n Û x = P ( x , y ) , Û y = Q ( x , y ) (2.1)having a center at the point ¯ x ∈ R . Without loss of generalitywe can assume that the center ¯ x of system (2.1) is the origin ofcoordinates. The following definition is due to Poincaré ([4], Sect.2). Definition 2.1.
We say that an isolated singular point ¯ x of (2.1) isa center if there exists a neighbourhood of ¯ x , such that every orbitin this neighbourhood is a cycle surrounding ¯ x . Remark . Determining the conditions on the parameters underwhich the origin for system (2.1) is a center is the well-known cen-ter problem , see [45, 51]. There are many partial results for the cen-ters of system (2.1) of degree n ≥
2. Unfortunately, at present, weare very far from obtaining the classification of all the centers ofcubic polynomial differential systems. In general, the huge numberof computations necessary for obtaining complete classification be-comes the central problem which is computationally intractable,see for instances, [15] and the references cited therein.Now consider the perturbations of (2.1) of the form Û x = P ( x , y ) + p ( x , y , ε ) , Û y = Q ( x , y ) + q ( x , y , ε ) , (2.2)where the polynomials p , q are of degree at most n (usually n ≥ n ) in x and y , and ε is a small parameter. We are interested inthe maximum number of medium amplitude limit cycles of (2.2)for | ε | > drdθ = k Õ i = ε i F i ( θ , r ) + ε k + R ( θ , r , ε ) , (2.3)where F i : R × D → R for i = , , . . . , k , and R : R × D ×(− ε , ε ) → R are C k functions, 2 π -periodic in the first variable, being D anopen and bounded interval of ( , ∞) , and ε is a small parameter.As one of the main hypotheses, it is assumed that r ( θ , z ) is a 2 π -periodic solution of the unperturbed differential system dr / dθ = F ( θ , r ) , for every initial condition r ( , z ) = z ∈ D .The averaging method consists in defining a collection of func-tions f i : D → R , called the i -th order averaged functions, for i = , , . . . , k , which control (their simple zeros control), for ε sufficiently small, the isolated periodic solutions of the differentialsystem (2.3). In Llibre-Novaes-Teixeira [41] it has been establishedthat f i ( z ) = y i ( π , z ) i ! , (2.4)where y i : R × D → R , for i = , , . . . , k , is defined recursively bythe following integral equations y ( θ , z ) = ∫ θ (cid:16) F ( s , r ( s , z )) + ∂ F ( s , r ( s , z )) y ( s , z ) (cid:17) ds , y i ( θ , z ) = i ! ∫ θ F i ( s , r ( s , z )) + i Õ ℓ = Õ S ℓ b ! b !2! b · · · b ℓ ! ℓ ! b ℓ · ∂ L F i − ℓ ( s , r ( s , z )) ℓ Ö j = y j ( s , z ) b j ! ds , (2.5)where S ℓ is the set of all ℓ -tuples of nonnegative integers [ b , b , . . . , b ℓ ] satisfying b + b + · · · + ℓ b ℓ = ℓ and L = b + b + · · · + b ℓ . Here, ∂ L F ( θ , r ) denotes the Fréchet’s derivative of order L with respectto the variable r . n [14, 41] the averaging method at any order was developedto study isolated periodic solutions of nonsmooth but continuousdifferential systems. Recently, the averaging method has also beenextended to study isolated periodic solutions of discontinuous dif-ferential systems; see [28, 37, 40]. The following k -th order averag-ing theorem is proved in Llibre-Novaes-Teixeira [41]. Theorem 2.3.
Assume that the following conditions hold:(a) for each i = , , . . . , k and θ ∈ R , the function F i ( θ , ·) is ofclass C k − i , ∂ k − i F i is locally Lipschitz in the second variable, and R ( θ , · , ε ) is a continuous function locally Lipschitz in the second vari-able;(b) f i ≡ for i = , , . . . , j − and f j , with j ∈ { , , . . . , k } ;(c) for some z ∗ ∈ D with f j ( z ∗ ) = , there exists a neighbor-hood V ⊂ D of z ∗ such that f j ( z ) , for all z ∈ ¯ V \{ z ∗ } , and that d B ( f j ( z ) , V , ) , .Then, for | ε | > sufficiently small, there exists a π -periodic so-lution r ε ( θ ) of (2.3) such that r ε ( ) → z ∗ when ε → .Remark . The above symbol d B denotes the Browder degree; seeBrowder [2] for a general definition. When f j is a C function andthe derivative of f j at z ∈ V is distinct from zero (i.e., f ′ j ( z ) , f ′ j ( z ∗ ) , d B ( f j ( z ) , V , ) , ℓ and m positive integers, werecall the Bell polynomials: B ℓ, m ( x , . . . , x ℓ − m + ) = Õ ˜ S ℓ, m ℓ ! b ! b ! · · · b ℓ − m + ! ℓ − m + Ö j = (cid:18) x j j ! (cid:19) b j , where ˜ S ℓ, m is the set of all ( ℓ − m + ) -tuples of nonnegative integers [ b , b , . . . , b ℓ − m + ] satisfying b + b + · · · + ( ℓ − m + ) b ℓ − m + = ℓ ,and b + b + · · · + b ℓ − m + = m .The following result is an equivalent formulation of the integralequation (2.5) via above Bell polynomials, its proof can be foundin [48]. Theorem 2.5.
For i = , , . . . , k the recursive equation (2.5) reads y ( θ , z ) = Y ( θ , z ) ∫ θ Y ( s , z ) − F ( s , r ( s , z )) ds , y i ( θ , z ) = Y ( θ , z ) ∫ θ Y ( s , z ) − i ! F i ( s , r ( s , z )) + i Õ m = ∂ m F ( s , r ( s , z )) B i , m (cid:0) y ( s , z ) , . . . , y i − m + ( s , z ) (cid:1) + i − Õ ℓ = ℓ Õ m = i ! ℓ ! ∂ m F i − ℓ ( s , r ( s , z )) B ℓ, m (cid:0) y ( s , z ) , . . . , y ℓ − m + ( s , z ) (cid:1)! ds , (2.6) where Y ( θ , z ) is the fundamental solution of the variational equation Y ′ = ∂ F ( θ , r ( θ , z )) Y satisfying the initial condition Y ( , z ) = . The general study of the exact number of simple zeros of the av-eraged functions (2.4) up to every order is also very difficult to bedone, since the averaged functions may be too complicated, suchas including square root functions, logarithmic functions, and theelliptic integrals. In the literature there is an abundance of papersdealing with zeros of the averaged functions (see, for instance, [20–22, 34, 42] and references therein). Note that one can estimate thesize of bifurcated limit cycles by using the expressions of the aver-aged functions. In fact we know that if the averaged functions f j = j = , . . . , k − f k ,
0, and ¯ z ∈ D is a simple zero of f k ,then by Theorem 2.3 there is a limit cycle r ε ( θ ) of differential sys-tem (2.3) such that r ε ( θ ) = r ( θ , ¯ z ) + O( ε ) . Then, going back throughthe changes of variables we have for the differential system (2.2)the limit cycle ( x ( t , ε ) , y ( t , ε )) = ( r ( θ , ¯ z ) cos θ , r ( θ , ¯ z ) sin θ ) + O( ε ) . The process of using the averaging method for studying limit cy-cles of differential systems can be divided into three steps ([24],Sect. 4).
STEP 1 . Write the perturbed system (2.2) in the normal form ofaveraging (2.3) up to k -th order in ε . STEP 2 . (i) Compute the exact formula for the k -th order inte-gral function y k ( θ , z ) in (2.6). (ii) Derive the symbolic expressionof the k -th order averaged function f k ( z ) by (2.4). STEP 3 . Determine the exact upper bound of the number ofsimple zeros of f k ( z ) for z ∈ D .In the following subsections we will present algorithms to im-plement the first two steps. We use “Maple-like” pseudo-code, basedon our Maple implementation. Using these algorithms we reducethe problem of studying the number of limit cycles of system (2.2)to the problem of detecting STEP 3 . In this subsection we will devise an efficient algorithm which canbe used to transform system (2.2) into the form (2.3).We first describe the underlying equations before presenting thealgorithm. Doing the change of polar coordinates x = rC , y = rS with C = cos θ and S = sin θ , then we can transform system (2.2)into the following form drdθ = dr / dtdθ / dt = r ( x Û x + y Û y ) x Û y − y Û x (cid:12)(cid:12)(cid:12) x = rC , y = r S = r C ( P ( x , y ) + p ( x , y , ε )) + S ( Q ( x , y ) + q ( x , y , ε )) C ( Q ( x , y ) + q ( x , y , ε )) − S ( P ( x , y ) + p ( x , y , ε )) (cid:12)(cid:12)(cid:12) x = rC , y = r S = r CP ( x , y ) + SQ ( x , y ) CQ ( x , y )− SP ( x , y ) + Cp ( x , y , ε ) + Sq ( x , y , ε ) CQ ( x , y )− SP ( x , y ) + Cq ( x , y , ε )− Sp ( x , y , ε ) CQ ( x , y )− SP ( x , y ) (cid:12)(cid:12)(cid:12) x = rC , y = r S = F ( θ , r ) + εF ( θ , r ) + . . . + ε k F k ( θ , r ) + O( ε k + ) . (3.1)The last equality is obtained by carrying the order k + ariable ε , around the point ε =
0. The first algorithm
Normal-Form , presented below, is a direct implementation of the formuladerivation in (3.1).
Algorithm 1 NormalForm ( P , Q , p , q , k ) Input: a perturbed system (2.2) with an order k ≥ Output: an expression of dr / dθ up to k -th order in ε d = normal ( subs ( x = r · C , y = r · S , x · ( P + p ) + y · ( Q + q ))/ r ) ; d = normal ( subs ( x = r · C , y = r · S , x · ( Q + q ) − y · ( P + p ))/ r ) ; T : = taylor ( d / d , ε = , k + ) ; H : = convert ( T , polynom ) ; F : = coeff ( ε · H , ε ) ; for i from to k do c i : = coeff ( H , ε i ) ; F i , : = prem (cid:0) numer ( c i ) , C + S − , S (cid:1) ; F i , : = prem (cid:0) denom ( c i ) , C + S − , S (cid:1) ; F i : = F i , / F i , ; dr / dθ : = subs ( C = cos θ , S = sin θ , F + Í kj = F j ε j ); return dr / dθ ; In line 8 the function prem ( a , b , x ) is the pseudo-remainder of a with respect to b in the variable x . By the property of the pseudo-remainder we know that the degree in S is at most 1 of the polyno-mials F i , and F i , . This subsection is devoted to provide effective algorithms to com-pute the formula and exact expression of the k -th order averagedfunction. According to Theorem 2.5, we should take the followingsubsteps to compute the k -th order averaged function of system(2.3). Substep 1 . Determine the open and bounded interval D , the2 π -periodic solution r ( θ , z ) of the unperturbed system dr / dθ = F ( θ , r ) with initial condition r ( , z ) = z ∈ D , and the fundamentalsolution Y ( θ , z ) of the variational equation Y ′ = ∂ F ( θ , r ( θ , z )) Y with initial condition Y ( , z ) = Substep 2 . Compute the exact formula for the k -th order inte-gral function y k ( θ , z ) . Substep 3 . Output the symbolic expression for the k -th order av-eraged function f k ( z ) (simplified by using the conditions for f ≡ f ≡ · · · ≡ f k − ≡
0) for a given differential system (2.2).We provide each of the substep an algorithm. For the Substep 1we first derive the 2 π -periodic solution r ( θ , z ) , and then use it tofurther obtain the interval D and the fundamental solution Y ( θ , z ) . Algorithm 2 DSolutions ( F ) Input: the unperturbed term F in (2.3) Output: r ( θ , z ) , a set of inequalities ( SIs ) with respect to z , and Y ( θ , z ) de1 : = diff ( r ( θ ) , θ ) = subs ( r = r ( θ ) , F ) ; ans1 : = dsolve ({ de1 , r ( ) = z } , r ( θ )) ; r ( θ , z ) : = op ( , ans1 ) ; minvalue : = minimize ( r ( θ , z ) , θ = .. π ) ; m : = nops ([ op ( minvalue )]) ; SIs : = {} ; for i from to m do SIs : = SIs union { op ( i , minvalue ) > } ; de2 : = diff ( Y ( θ ) , θ ) = subs ( r = r ( θ , z ) , diff ( F , r )) · Y ( θ ) ; ans2 : = dsolve ({ de2 , Y ( ) = } , Y ( θ )) ; Y ( θ , z ) : = op ( , ans ) ; return [ r ( θ , z ) , SIs , Y ( θ , z ) ]; Remark . The output results of r ( θ , z ) and Y ( θ , z ) can be reducedby using the identity sin θ + cos θ = θ is at most 1. We use the routine dsolve built-inMaple for solving an ordinary differential equation. We remarkthat the unperturbed term F ( θ , r ) is usually a rational trigonomet-ric function in r , sin θ and cos θ . As far as we know, we do not havea systematic approach to the solution of the differential equation dr / dθ = F ( θ , r ) in the general case. In Section 4 we will considercertain classes of differential systems with uniform isochronouscenters to illustrate the effectiveness of our algorithm. It is impor-tant to emphasize that the interval D can be determined by usingthe output set SIs . Since the original system may contain someparameters, the resulting set
SIs could be parametric. In order toderive the interval D in this case, we will construct an equivalentsolution set SIs of SIs that contains only the rational polynomialinequalities, and then use the
SemiAlgebraic command in Mapleto compute the solutions. Below we provide a concrete example toshow the feasibility this algorithm, one may check the results in[34]. More experiments can be found in Section 4.
Example . Consider the following quintic polynomial differentialsystem Û x = − y + x y ( x + y ) , Û y = x + xy ( x + y ) . (3.2)Applying our algorithm NormalForm for p = q = k =
0, we have dr / dθ = r cos θ sin θ . Then applying the algorithm DSolutions we obtain a list [ r ( θ , z ) , SIs , Y ( θ , z ) ], where r ( θ , z ) = z ( z ( cos θ − ) + ) / , SIs = (cid:26) < z , < z (− z + ) / (cid:27) , Y ( θ , z ) = ( z ( cos θ − ) + ) / . To obtain the interval D in this case, we construct an equivalentsolution set SIs of SIs that contains only the rational polynomials:
SIs : = { < z , < z − z + } . Then using the Maple command Sol-veTools[SemiAlgebraic] , we compute the solution of the set
SIs ,and obtain that D = { < z < − / } . e want to say that the expressions of the returned results on r ( θ , z ) and Y ( θ , z ) may be complicated, such as including squareroot functions, and exponential functions. Below we give a simpleexample to show this, one may find the related results in [33]. Example . Consider the following polynomial differential system Û x = − y ( x + y ) , Û y = x ( x − y ) . (3.3)The normal form dr / dθ = − r cos θ sin θ can be obtained by thealgorithm NormalForm , and the algorithm
DSolutions returns alist (cid:2) ze cos θ − , (cid:8) < z , < ze − (cid:9) , e cos θ − (cid:3) .For the Substep 2, we present our algorithm AveragingFor-mula . This algorithm can be used to compute the exact formulaof the k -th order integral function y k ( θ , z ) . Correctness of it fol-lows from Theorem 2.5. Algorithm 3 AveragingFormula ( k ) Input: an order k ≥ Output: the integral function y k ( θ , z ) T : = T : = for m from to k do T : = T + Diff ( F ( s , r ( s , z )) , r $ m ) · IncompleteBellB ( k , m , y ( s , z ) , . . . , y k − m + ( s , z )) ; for ℓ from to k − do for m from to ℓ do T : = T + k ! ℓ ! · Diff ( F k − ℓ ( s , r ( s , z )) , r $ m ) · IncompleteBellB ( ℓ, m , y ( s , z ) , . . . , y ℓ − m + ( s , z )) ; y k ( θ , z ) : = Y ( θ , z )· Int (cid:0) Y − ( s , z ) · ( k ! · F k ( s , r ( s , z )) + T + T ) , s = .. θ (cid:1) ; return y k ( θ , z ) ;We deduce explicitly the formulae of y k ’s up to k = y k ’s. In Section 4, we will study several differentialsystems to show the feasibility of our algorithm.In the last subsection, we provide an algorithm NormalForm to transform system (2.2) into the form dr / dθ . The algorithm DSo-lutions admits one to obtain the fundamental solutions r ( θ , z ) , Y ( θ , z ) and the interval D (Substep 1). The algorithm Averaged-Function , presented below, is based on the algorithms
Normal-Form , DSolutions and Theorem 2.5, which provides a straight-forward calculation method to derive the exact expression of the k -th order averaged function for a given differential system in theform (2.2) (Substep 3). Algorithm 4 AveragedFunction ( dr / dθ , r ( θ , z ) , Y ( θ , z ) , k ) Input: a normal formal of averaging (3.1) with an order k ≥ r ( θ , z ) , Y ( θ , z ) Output: a list of expressions of the averaged functions F : = coeff ( ε · ( dr / dθ ) , ε ) ; for j from to k do F j : = coeff ( dr / dθ , ε j ) ; A j : = AFormula ( j ) ; H j : = normal (cid:16) Y ( θ , z ) · expand ( subs ( r = r ( θ , z ) , value ( A j ))) (cid:17) ; H j , : = collect ( expand ( numer ( H j )) , { cos θ , sin θ } , distributed ) ; H j , : = denom ( H j ) ; for h from to nops ( H j , ) do д j , h : = int (cid:16) op ( h , H j , ) H j , , θ = .. θ , AllSolutions (cid:17) ; s j , h : = int (cid:16) op ( h , H j , ) H j , , θ = .. π (cid:17) ; y j : = Y ( θ , z ) · sum ( д j , t , t = .. nops ( H j , )) ; f j : = j ! · sum ( s j , t , t = .. nops ( H j , )) ; return [ y k , f k ]; In line 4, the routine
AFormula is a subalgorithm we use for thegeneration of the expression in the parenthesis of equation (2.6)without dependence on ( s , z ) . The detailed information of this sub-algorithm is as follows. Subalgorithm: AFormula
INPUT: An averaging order k ≥ ( s , z ) .STEP 0. U = V = For m from to k do U : = U + Diff ( F , r $ m ) · IncompleteBellB ( k , m , seq ( y i , i = .. k − m + )) ; end do ;STEP 2. For ℓ from to k − for m from to ℓ do V : = V + k ! ℓ ! · Diff ( F k − ℓ , r $ m ) · IncompleteBellB ( ℓ, m , seq ( y i , i = ..ℓ − m + )) ; end do ; end do ;STEP 3. Output k ! F k + U + V . Remark . In order to obtain an exact and simplified expressionof the averaged function, one should make some assumptions (e.g.,the interval D on z , and possible conditions on the parameters thatmay appear in the original differential systems) before preformingthe algorithm AveragedFunction . For more details see our exper-iments in Section 4. We also remark that, throughout the compu-tation, an assumption on θ (i.e., θ ∈ ( π − ϵ , π + ϵ ) with ϵ a smallnumber) was made to identify a valid branch of the possible re-turned piecewise functions (in line 9), since the integral functions y i ( θ , z ) for i = , . . . , k evaluate at the point θ = π in (2.4).We implemented all the algorithms presented in this section inMaple. In the next section, we will apply our general algorithmicapproach to analyze the bifurcation of limit cycles for several con-crete differential systems. In this section we demonstrate our algorithmic tests using severalexamples. We present the bifurcation of limit cycles for a class ofgeneralized Kukles polynomial differential systems as an illustra-tion of our approach explained in previous sections. In addition, we tudy the number of limit cycles that bifurcate from some periodicsolutions surrounding the isochronous centers for certain differen-tial systems by the first and second order averaging method. Theobtained results of our experiments extend the existing relevantresults and show the feasibility of our approach. In this subsection we consider a very particular case of the 16thHilbert problem; we study the number of limit cycles of the gener-alized Kukles polynomial differential system Û x = − y , Û y = x + Q ( x , y ) , (4.1)where Q ( x , y ) is a polynomial with real coefficients of degree n .This system was introduced by Kukles in [30], examining the con-ditions under which the origin of the system Û x = − y , Û y = x + a x + a xy + a y + a x + a x y + a xy + a y is a center. For long time, it had been thought that the conditionsgiven by Kukles were necessary and sufficient conditions, but somenew cases have been found, see [7, 29].Here we are interested in studying the maximum number oflimit cycles that bifurcate from the periodic orbits of the linearcenter Û x = − y , Û y = x , perturbed inside the following class of gen-eralized Kukles polynomial differential systems Û x = − y + Õ k ≥ ε k l km ( x ) , Û y = x − Õ k ≥ ε k (cid:16) f kn ( x ) + д kn ( x ) y + h kn ( x ) y + d k y (cid:17) , (4.2)where for every k the polynomials l km ( x ) , f kn ( x ) , д kn ( x ) , and h kn ( x ) have degree m , n , n , and n respectively, d k , ε is a small parameter. This question has been studied in [46]for k = ,
2, and the authors obtained the following result.
Theorem 4.1.
Assume that for k = , the polynomials l km ( x ) , f kn ( x ) , д kn ( x ) , and h kn ( x ) have degree m , n , n , and n respectively,with m , n , n , n ≥ , and d k , is a real number. Then for ε sufficiently small the maximum number of limit cycles of the Kuklespolynomial system (4.2) bifurcating from the periodic orbits of thelinear center Û x = − y , Û y = x , (1) is max (cid:8) (cid:2) m − (cid:3) , (cid:2) n (cid:3) , (cid:9) by using the first order averagingmethod; (2) is max (cid:8) (cid:2) n (cid:3) + h n − i , (cid:2) n (cid:3) + (cid:2) m (cid:3) − , h n + i , h n + i , (cid:2) n (cid:3) + (cid:2) m (cid:3) , h n + i + (cid:2) n (cid:3) , (cid:2) n (cid:3) , (cid:2) m − (cid:3) , h n − i + µ , h n + i + µ , (cid:9) by using the second order averaging method,where µ = min (cid:8) (cid:2) m − (cid:3) , (cid:2) n (cid:3) (cid:9) . Here, [ · ] denotes the integer part function. Remark that, manyresearchers have discussed the bifurcation of limit cycles for gener-alized Kukles polynomial differential system in the form (4.1). Werefer the readers to [36, 47] for some interesting results on thissubject. The next result extends Theorem 4.1 to arbitrary order ofaveraging. Lemma 4.2.
Let max { m , n , n + , n + } = N ≥ , then theKukles polynomial system (4.2) for ε sufficiently small has no morethan [ k ( N − )/ ] limit cycles bifurcating from the periodic orbits ofthe linear center Û x = − y , Û y = x , using the averaging method up toorder k . Proof.
This result follows directly from Theorem 6 in [14]. (cid:3)
In what follows, using our algorithms we will do some experi-mental results by fixing some values of the degrees in system (4.2).Note that the maximum numbers of limit cycles in Theorem 4.1and Lemma 4.2 may not be reached. The following corollary showsthat these maximum numbers can be reached for some orders ofaveraging.
Corollary 4.3. (i) When m = , n = , n = , and n = , themaximum number of limit cycles of the Kukles polynomial system (4.2) bifurcating from the periodic orbits of the linear center Û x = − y , Û y = x , using the fifth order averaging method is five and it is reached.(ii) When m = , n = , n = , and n = , the maximum num-ber of limit cycles of the Kukles polynomial system (4.2) bifurcatingfrom the periodic orbits of the linear center Û x = − y , Û y = x , using thefourth order averaging method is five and it is reached. The detailed proof of the first statement of Corollary 4.3 can befound in Appendix B. Since the calculations and arguments of thesecond part are quite similar to those used in the first one, we omitthe proof of statement (ii) in Corollary 4.3. More concretely, weprovide in Table 1 the maximum number of limit cycles for system(4.2) in each case of Corollary 4.3 up to the k -th order averagingmethod for k = , . . . , Table 1: Number of limit cycles of system (4.2) in Corollary4.3
Averaging order Statement (i) Statement (ii)1 1 22 2 23 3 44 4 55 5 -The number of limit cycles in statement (i) can be reached foreach order of averaging. That is to say, the bound given in Lemma4.2 is sharp for the case in statement (i). However, for the statement(ii), the bound given in Lemma 4.2 is only sharp for the first orderof averaging. We note also that for each statement in Corollary4.3, the bound provided in Theorem 4.1 can be reached up to thesecond order.
Remark . The calculation of the high order averaged function f k involves heavy computations with complicated expressions. It maynot work effectively when one of the degrees ( m and n i , i = , , dr / dθ using the conditions on the parameters of f ≡ f ≡ · · · ≡ f k − = .2 Limit cycles for certain differential systemswith uniform isochronous centers Recall that a center ¯ x of system (2.1) is an isochronous center if ithas a neighborhood such that in this neighborhood all the periodicorbits have the same period. An isochronous center is uniform ifin polar coordinates x = r cos θ , y = r sin θ , it can be written as Û r = G ( θ , r ) , Û θ = η , η ∈ R \{ } , see Conti [10] for more details. Thenext result on the uniform isochronous center (UIC) is well-known,a proof of it can be found in [27]. Proposition 4.5.
Assume that system (2.1) has a center at theorigin ¯ x . Then ¯ x is a UIC if and only if by doing a linear change ofvariables and a rescaling of time the system can be written as Û x = − y + x f ( x , y ) , Û y = x + y f ( x , y ) , (4.3) where f is a polynomial in x and y of degree n − , and f ( , ) = . In what follows, we recall some important results on the UICsof planar cubic and quartic differential systems. The following re-sult due to Collins [9] in 1997, also obtained by Devlin, Lloyd andPearson [11] in 1998, and by Gasull, Prohens and Torregrosa [12]in 2005 characterizes the UICs of cubic polynomial systems.
Theorem 4.6.
A planar cubic differential system has a UIC at theorigin if and only if it can be written as system (4.3) with f ( x , y ) = a x + a y + a x + a xy − a y satisfying that a a − a a + a a a = . Moreover, this planar cubic differential system can be reduced toeither one of the following two forms: Û x = − y + x y , Û y = x + xy , (4.4) Û x = − y + x + Ax y , Û y = x + xy + Axy , (4.5) where A ∈ R . Systems (4.4) and (4.5) are known as
Collins First Form and
CollinsSecond Form , respectively. See ([35], Thm. 9) for more details of theglobal phase portraits of the Collins forms.In order to save space, we put the remaining results in AppendixC.
We have presented a systematical approach to analyze how manylimit cycles of differential system (2.2) can bifurcate from the pe-riodic orbits of an unperturbed one via the averaging method. Wedesigned four algorithms to analyze the averaging method andshown that the general study of the number of limit cycles of sys-tem (2.2) can be reduced to the problem of estimating the numberof simple zeros of the obtained averaged functions with the aid ofthese algorithms.Our algorithms admit a generalization to the case of studyingthe bifurcation of limit cycles for discontinuous differential sys-tems. It would be interesting to employ our approach to analyzethe bifurcation of limit cycles for differential systems in many dif-ferent fields, which are of high interest in nature sciences and en-gineering. It will be beneficial to generalize our current approachto the case of higher dimension differential systems by using thegeneral form of the averaging method. We leave this as the futureresearch problems. In addition, we noticed the phenomenon of tremendous growthof expressions in intermediate calculations while we done exper-iments for the linear center Û x = − y , Û y = x by using the high or-der of averaging. For the nonlinear polynomial differential centers,the evaluation of the high order averaged functions is highly non-trivial; the main difficulty exists in the technical and cumbersomecomputations of some complicated integral equations. How to sim-plify and optimize the steps of the computations of the averagedfunctions is also a question that remains for further investigation. ACKNOWLEDGMENTS
Huang’s work is partially supported by China Scholarship Councilunder Grant No.: 201806020128. The author is grateful to ProfessorChee Yap and Professor Dongming Wang for their profound con-cern and encouragement.
REFERENCES [1] Nikolai N. Bautin. 1954. On the number of limit cycles which appear with thevariation of the coefficients from an equilibrium position of focus or center type.
Amer. Math. Soc. Transl.
100 (1954), 397–413.[2] Felix E. Browder. 1983. Fixed point theory and nonlinear problems.
Bull. Amer.Math. Soc.
Int. J. Bifur. Chaos
11 (2001), 711–722.[4] Javier Chavarriga and Marco Sabatini. 1999. A survey of isochronous centers.
Qual. Theory Dyn. Syst.
Acta Math. Sinica
22 (1979),751–758.[6] Colin J. Christopher and Cheng Z. Li. 2007.
Limit Cycles of Differential Equations .Birkhäuser, Boston.[7] Colin J. Christopher and Noel G. Lloyd. 1990. On the paper of Jin and Wangconcerning the conditions for a centre in certain cubic systems.
Bull. LondonMath. Soc.
22 (1990), 5–12.[8] Colin J. Christopher and Noel G. Lloyd. 1995. Polynomial systems: A lowerbound for the Hilbert numbers.
Proc. R. Soc. Lond. Ser. A
450 (1995), 218–224.[9] Christopher B. Collins. 1997. Conditions for a centre in a simple class of cubicsystems.
Differential Integral Equations
10 (1997), 333–356.[10] Roberto Conti. 1994. Uniformly isochronous centers of polynomial systems in R . In Differential equations, dynamical systems, and control science. LectureNotes in Pure and Appl. Math. 152, New York, 21–31.[11] James Devlin, Noel G. Lloyd, and Jane M. Pearson. 1998. Cubic systems and Abelequations. J. Differ. Equ.
147 (1998), 435–454.[12] Armengol Gasull, Rafel Prohens, and Joan Torregrosa. 2005. Limit cycles forrigid cubic systems.
J. Math. Anal. Appl.
303 (2005), 391–404.[13] Armengol Gasull and Joan Torregrosa. 2001. A new algorithm for the compu-tation of the Lyapunov constants for some degenerated critical points.
Nonlin.Anal.
47 (2001), 4479–4490.[14] Jaume Giné, Maite Grau, and Jaume Llibre. 2013. Averaging theory at any orderfor computing periodic orbits.
Phys. D
250 (2013), 58–65.[15] Jaume Giné and Xavier Santallusia. 2004. Implementation of a new algorithmof computation of the Poincaré-Liapunov constants.
J. Comput. Appl. Math.
J. Differ. Equ.
252 (2012), 3278–3304.[17] Mao A. Han, Jun M. Yang, and Pei Yu. 2009. Hopf bifurcation for near-Hamiltonian.
Int. J. Bifur. Chaos
19 (2009), 4117–4130.[18] Mao A. Han, Tong H. Zhang, and Hong Zang. 2004. On the number and distri-bution of limit cycles in a cubic system.
Int. J. Bifur. Chaos
14 (2004), 4285–4292.[19] David Hilbert. 1902. Mathematical problems.
Bull. Am. Math. Soc.
Int. J. Bifur. Chaos
27 (2017), 1750072–1–16.[21] Bo Huang. 2019. Limit cycles for a discontinuous quintic polynomial differentialsystem.
Qual. Theory Dyn. Syst.
18 (2019), 769–792.[22] Bo Huang. 2020. On the limit cycles for a class of discontinuous piecewise cubicpolynomial differential systems.
Electron. J. Qual. Theory Differ. Equ.
25 (2020),1–24.[23] Bo Huang and Wei Niu. 2019. Limit cycles for two classes of planar polynomialdifferential systems with uniform isochronous centers.
J. Appl. Anal. Comput. Bull. Am.Math. Soc.
39 (2002), 301–354.[26] Jackson Itikawa and Jaume Llibre. 2015. Limit cycles bifurcating from the periodannulus of a uniform isochronous center in a quartic polynomial differentialsystem.
Electron. J. Differential Equations
246 (2015), 1–11.[27] Jackson Itikawa and Jaume Llibre. 2015. Phase portraits of uniform isochronousquartic centers.
J. Comput. Appl. Math.
287 (2015), 98–114.[28] Jackson Itikawa, Jaume Llibre, and Douglas D. Novaes. 2017. A new result onaveraging theory for a class of discontinuous planar differential systems withapplications.
Rev. Mat. Iberoam.
33 (2017), 1247–1265.[29] Xiao F. Jin and Dong M. Wang. 1990. On the conditions of Kukles for the exis-tence of a centre.
Bull. London Math. Soc.
22 (1990), 1–4.[30] Isaak S. Kukles. 1944. Sur quelques cas de distinction entre un foyer et un centre.
Dokl. Akad. Nauk. SSSR.
42 (1944), 208–211.[31] Cheng Z. Li, Chang J. Liu, and Jia Z. Yang. 2009. A cubic system with thirteenlimit cycles.
J. Differ. Equ.
246 (2009), 3609–3619.[32] Ji B. Li. 2003. Hilbert’s 16th problem and bifurcations of planar polynomialvector fields.
Int. J. Bifur. Chaos
13 (2003), 47–106.[33] Shi M. Li, Yu L. Zhao, and Zhao H. Sun. 2015. On the limit cycles of planarpolynomial system with non-rational first integral via averaging method at anyorder.
Appl. Math. Comput.
256 (2015), 876–880.[34] Hai H. Liang, Jaume Llibre, and Joan Torregrosa. 2016. Limit cycles coming fromsome uniform isochronous centers.
Adv. Nonlinear Stud.
16 (2016), 197–220.[35] Jaume Llibre and Jackson Itikawa. 2015. Limit cycles for continuous and discon-tinuous perturbations of uniform isochronous cubic centers.
J. Comput. Appl.Math.
277 (2015), 171–191.[36] Jaume Llibre and Ana C. Mereu. 2011. Limit cycles for generalized Kukles poly-nomial differential systems.
Nonlin. Anal.
74 (2011), 1261–1271.[37] Jaume Llibre, Ana C. Mereu, and Douglas D. Novaes. 2015. Averaging theoryfor discontinuous piecewise differential systems.
J. Differ. Equ.
258 (2015), 4007–4032.[38] Jaume Llibre, Ana C. Mereu, and Marco A. Teixeira. 2010. Limit cycles of thegeneralized polynomial Liénard differential equations.
Math. Proc. Camb. Phil.Soc.
148 (2010), 363–383.[39] Jaume Llibre, Richard Moeckel, and Carles Simó. 2015.
Central Configurations,Periodic Orbits, and Hamiltonian Systems . Birkhäuser, Basel.[40] Jaume Llibre, Douglas D. Novaes, and Camila A.B. Rodrigues. 2017. Averagingtheory at any order for computing limit cycles of discontinuous piecewise dif-ferential systems with many zones.
Phys. D
Nonlinearity
27 (2014), 563–583.[42] Jaume Llibre and Grzegorz Świrszcz. 2011. On the Limit cycles of polynomialvector fields.
Dyn. Contin. Discrete Impuls. Syst. Ser A: Math Anal.
18 (2011),203–214.[43] Noel G. Lloyd. 1988. Limit cycles of polynomial systems-some recent develop-ments.
London Math. Soc. Lecture Note Ser.
127 (1988), 192–234.[44] Noel G. Lloyd and Stephen Lynch. 1988. Small-amplitude limit cycles of certainLiénard systems.
Proc. R. Soc. Lond. A
418 (1988), 199–208.[45] Adam Mahdi, Claudio Pessoa, and Jonathan D. Hauenstein. 2017. A hybridsymbolic-numerical approach to the center-focus problem.
J. Symb. Comput.
82 (2017), 57–73.[46] Nawal Mellahi, Amel Boulfoul, and Amar Makhlouf. 2019. Maximum number oflimit cycles for generalized Kukles polynomial differential systems.
Differ. Equ.Dyn. Syst.
27 (2019), 493–514.[47] Ana C. Mereu, Regilene Oliveira, and Camila A.B. Rodrigues. 2018. Limit cyclesfor a class of discontinuous piecewise generalized Kukles differential systems.
Nonlin. Dyn.
93 (2018), 2201–2212.[48] Douglas D. Novaes. 2017.
An Equivalent Formulation of the Averaged Functionsvia Bell Polynomials . Springer, New York, 141–145.[49] Lin P. Peng and Zhao S. Feng. 2014. Bifurcation of limit cycles from quarticisochronous systems.
Electron. J. Differential Equations
95 (2014), 1–14.[50] Valery G. Romanovski. 1993. Calculation of Lyapunov numbers in the case oftwo pure imaginary roots.
Differ. Equ.
29 (1993), 782–784.[51] Valery G. Romanovski and Douglas S. Shafer. 2009.
The Center and CyclicityProblems: A Computational Algebra Approach . Birkhäuser, Boston.[52] Jan A. Sanders, Ferdinand Verhulst, and James Murdock. 2007.
Averaging Meth-ods in Nonlinear Dynamical Systems . Springer, New York.[53] Song L. Shi. 1980. A concrete example of the existence of four limit cycles forquadratic system.
Sci. Sinica
23 (1980), 153–158.[54] Dong M. Wang. 1990. A class of cubic differential systems with 6-tuple focus.
J.Differ. Equ.
87 (1990), 305–315.[55] Dong M. Wang. 1991. Mechanical manipulation for a class of differential systems.
J. Symb. Comput.
12 (1991), 233–254.[56] Pei Yu and Guan R. Chen. 2008. Computation of focus values with applications.
Nonlin. Dyn.
51 (2008), 409–427.[57] Pei Yu and Mao A. Han. 2005. Twelve limit cycles in a cubic case of the 16thHilbert problem.
Int. J. Bifur. Chaos
15 (2005), 2191–2205.[58] Pei Yu and Mao A. Han. 2012. Four limit cycles from perturbing quadratic inte-grable systems by quadratic polynomials.
Int. J. Bifur. Chaos
22 (2012), 1250254–1–28. FIFTH ORDER AVERAGING FORMULAE
The explicit formulae of the functions y k ( θ , z ) for k = , , . . . , y ( θ , z ) = Y ( θ , z ) ∫ θ F ( s , r ( s , z )) Y ( s , z ) ds , y ( θ , z ) = Y ( θ , z ) ∫ θ Y ( s , z ) F ( s , r ( s , z )) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) ! ds , y ( θ , z ) = Y ( θ , z ) ∫ θ Y ( s , z ) F ( s , r ( s , z )) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) ! ds , y ( θ , z ) = Y ( θ , z ) ∫ θ Y ( s , z ) F ( s , r ( s , z )) + ∂ ∂ r F ( s , r ( s , z )) (cid:16) y ( s , z ) y ( s , z ) + y ( s , z ) (cid:17) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) ! ds , y ( θ , z ) = Y ( θ , z ) ∫ θ Y ( s , z ) F ( s , r ( s , z )) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) (cid:16) y ( s , z ) y ( s , z ) + y ( s , z ) (cid:17) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) y ( s , z ) + ∂∂ r F ( s , r ( s , z )) y ( s , z ) + ∂ ∂ r F ( s , r ( s , z )) y ( s , z ) ! ds . (A.1) B PROOF OF COROLLARY 4.3
Let l k ( x ) = Õ i = e k , i x i , f k ( x ) = Õ i = a k , i x i , д k ( x ) = Õ i = b k , i x i , h k ( x ) = Õ i = c k , i x i . (B.1)We consider the following perturbed system Û x = − y + Õ k = ε k l k ( x ) , Û y = x − Õ k = ε k (cid:16) f k ( x ) + д k ( x ) y + h k ( x ) y + d k y (cid:17) . (B.2)Now we study the number of limit cycles of system (B.2). Applyingour algorithm NormalForm by taking k = drdθ = Õ i = ε i F i ( θ , r ) + O( ε ) . (B.3)Here we give only the expression of F ( θ , r ) , the explicit expres-sions of F i ( θ , r ) for i = , . . . , F ( θ , r ) = S (cid:16) (− C a , + C c , − Cc , ) r + (− C a , + C c , − c , ) r − Ca , r − a , (cid:17) + (cid:16) ( C b , − C d , + C e , − C b , + C d , − d , ) r + ( C b , + C e , − Cb , ) r + ( C b , + C e , − b , ) r + Ce , (cid:17) ith C = cos θ and S = sin θ . For the case F =
0, the algorithm
DSolutions returns a list [ z , { < z } , AveragedFunction f ( z ) = − πz (cid:16) ( b , + d , − e , ) z + ( b , − e , ) (cid:17) . Therefore f ( z ) can have at most one positive simple real root.From Theorem 2.3 it follows that the first order averaging providesthe existence of at most one limit cycle of system (B.2) and thisnumber can be reached since the coefficients in f ( z ) are indepen-dent constants.From now on, for each k = , . . . ,
5, we will perform the calcu-lation of the averaged function f k under the hypothesis f j ≡ j = , . . . , k − e , = d , + b , / e , = b , , and computing f weobtain f ( z ) = − πz (cid:16) A , z + A , z + A , (cid:17) , where A , = − a , d , + b , c , + c , d , , A , = − a , d , + a , b , − a , e , + b , c , + b , c , + b , + d , − e , , A , = a , b , − a , e , + c , e , + b , − e , . It is obvious that f ( z ) can have at most two positive simple realroots. Therefore, the second order averaging provides the existenceof at most two limit cycles of system (B.2) and this number can bereached.Letting e , = A , / + e , , e , = A , / + e , , a , = A , / d , + a , , and computing f we have f ( z ) = πz d , (cid:16) A , z + A , z + A , z + A , (cid:17) , where A , = d , ( b , + d , )( b , d , − c , ) , the expressionsof A , i for i = , , f ( z ) has at most three positive simplereal roots. Then the third order averaging provides the existenceof at most three limit cycles of system (B.2) and this number canbe reached.To consider the fourth order averaging theorem we have twocases from the expression of A , . CASE 1 : [5 b , + d , = CASE 2 : [3 b , d , − c , = CASE 1 , be-cause the maximum number of limit cycles in Lemma 4.2 underthis case can be reached up to k = f ≡ e , = − A , / d , + e , , e , = − A , / d , + e , , a , = − A , / d , + a , , b , = − d , / f we obtain f ( z ) = − πz d , (cid:16) A , z + A , z + A , z + A , z + A , (cid:17) , where A , = c , d , ( c , + d , ) , here again we omit theexpressions of A , i for i = , , ,
3, because they are too long.Hence f ( z ) has at most four positive simple roots. Then the fourthorder averaging provides the existence of at most four limit cyclesof system (B.2) and this number can be reached. Note that d , ,
0, so in order to let f ≡
0, we take c , = e , = A , / d , + e , , e , = A , / d , + e , , a , = A , / d , + a , , b , = − A , / d , ( c , + d , ) + b , . Computing f we obtain f ( z ) = − πz d , (cid:16) A , z + A , z + A , z + A , z + A , z + A , (cid:17) , where A , = d , , and A , = − d , (cid:16) a , + a , c , − b , d , + b , + b , e , + c , − e , − c , (cid:17) . We do not explicitly provide the expressions of A , i for i = , , ,
3, because they are very long. Therefore f ( z ) can have atmost five positive real roots. Then the fifth order averaging pro-vides the existence of at most five limit cycles of system (B.2) andthis number can be reached. C LIMIT CYCLES FOR CERTAINDIFFERENTIAL SYSTEMS WITH UNIFORMISOCHRONOUS CENTERS
This appendix is an overflow from Subsection 4.2. The followingcharacterization of planar quartic polynomial differential systemswith an isolated UIC at the origin is provided by Chavarriga, Garcíaand Giné [3], in 2001.
Theorem C.1.
A planar quartic differential system has a UIC atthe origin if and only if it can be written as Û x = − y + x (cid:16) ax + bxy + cx + dxy (cid:17) , Û y = x + y (cid:16) ax + bxy + cx + dxy (cid:17) , (C.1) where a , b , c , d ∈ R . A classification of the global phase portraits of the quartic dif-ferential systems of the form (C.1) is provided in [27].In this subsection, we study the limit cycles that bifurcate fromsome periodic orbits surrounding the origin of the Collins FirstForm and a subclass of system (C.1) (we take a = b = C.1 Bifurcation of limit cycles of Collins FirstForm
Consider the following perturbation of the Collins First Form: Û x = − y + x y + ε Õ i + j = a i , j x i y j , Û y = x + xy + ε Õ i + j = b i , j x i y j . (C.2) e study the number of limit cycles of system (C.2) by using thefirst order averaging method. Applying our algorithm Normal-Form by taking k = drdθ = F ( θ , r ) + εF ( θ , r ) + O( ε ) , (C.3)where F ( θ , r ) = r cos θ sin θ , and F ( θ , r ) = S (cid:16) r (cid:0) a , − a , + b , − b , (cid:1) C + r (cid:0) − a , + b , − b , (cid:1) C + (cid:0) r (− a , + a , − b , ) + r (− a , − a , + a , − b , − b , + b , ) (cid:1) C + (cid:16) r ( a , − b , ) + r ( a , − b , − b , + b , ) (cid:17) C + (cid:16) r a , + r ( a , + a , + b , ) + r ( a , + b , ) (cid:17) C + r b , + b , (cid:17) + r (cid:0) a , − a , − b , + b , (cid:1) C + r (cid:0) a , − a , + b , (cid:1) C + (cid:0) r (− a , + a , + b , − b , ) + r (− a , − a , + a , + b , + b , − b , ) (cid:1) C + (cid:0) r (− a , + a , − b , ) + r (− a , − a , + a , − b , ) (cid:1) C + (cid:0) r ( a , − b , ) + r ( a , + a , − b , − b , + b , ) + r ( a , − b , ) (cid:1) C + (cid:0) r a , + r ( a , + a , + b , ) + a , (cid:1) C + r b , + rb , with S = sin θ , C = cos θ . Using the algorithm DSolutions weobtain a list [ r ( θ , z ) , SIs , Y ( θ , z ) ], where r ( θ , z ) = z p z ( cos θ − ) + , S Is = (cid:26) < z , < z √− z + (cid:27) , Y ( θ , z ) = ( z ( cos θ − ) + ) / . To obtain the interval D in this case, we construct an equivalentsolution set SIs : = { < z , < z − z + } . Computing the solution ofthe set SIs , we have D = { < z < } . Now assuming z ∈ D andapplying the algorithm AveragedFunction for the system (C.2),we obtain f ( z ) = − πz (cid:16)(cid:0) ( b , − b , − b , ) z + (− a , − a , − b , − b , + a , + b , ) z + a , + b , − a , − b , (cid:1) + p − z × (cid:0) ( a , − b , ) z − a , − b , + a , + b , (cid:1)(cid:17) . By taking the transformation z = √ − s , 0 < s <
1, we then have z f ( z ) = π ( − s ) (cid:16) ( b , − b , − b , ) s + ( b , − a , − b , + b , ) s + ( a , − a , + a , + b , ) s + a , + a , + a , (cid:17) . It is not hard to check that f ( z ) can have 3 simple zeros in theinterval ( , ) . Thus system (C.2) can have up to 3 limit cycles byusing the first order averaging method from Theorem 2.3. There-fore, using our algorithmic approach we verified the first result in([34], Thm. 1.1). In the following, we provide a concrete example todemonstrate that there exist some systems in the form (C.2) havingexactly three limit cycles. Consider a family of systems Û x = − y + x y + ε (cid:16) (− − b , − b , ) x + ( + b , ) x + (− + b , ) xy (cid:17) , Û y = x + xy + ε (cid:16) ( + b , + b , ) y + b , x y + b , y (cid:17) . (C.4)Applying our algorithmic approach to system (C.4), we obtain thefollowing first order averaged function f ( z ) z = √ − s = ⇒ π √ − s √ + s (cid:0) ( s − )( s − )( s − ) (cid:1) , (C.5)where z and s are defined as before. Apparently, f ( z ) has exactlythree positive zeros, denoted by z = √ , z = √ , z = √ , corresponding to s = / s = /
3, and s = /
5, respectively,in z ∈ ( , ) . Then it follows from Theorem 2.3 that system (C.4)has exactly three limit cycles by the first order averaging method.We explicitly provide the expressions of these three limit cycles asfollows: r θ , √ ! = √ √ θ + , r θ , √ ! = √ √ θ + , r θ , √ ! = √ √
24 cos θ + . Remark
C.2 . (i) Usually for the nonlinear differential systems, thehigher the averaging order is, the more complex are the computa-tional operations to calculate the averaged functions. Even for thefirst order averaging, if one consider a perturbation of a given dif-ferential system inside the class of polynomial differential systemsof the same degree (see [23, 26, 34]), the calculations of the aver-aged function require powerful computational resources. In nextsubsection, we will try to analyze a class of quartic systems by thesecond order averaging method. (ii) To our best knowledge, the bi-furcation of medium amplitude limit cycles for the Collins SecondForm has not been studied by the averaging method. Here we at-tempt to do this, but our algorithm DSolutions can not return thedesired results. In fact, for the Collins Second Form, one can obtainthe unperturbed normal form: dr / dθ = Ar cos θ sin θ + r cos θ bythe algorithm NormalForm . However, determining the explicitlyfundamental solutions r ( θ , z ) and Y ( θ , z ) for this parametric differ-ential system dr / dθ seems to be very hard. .2 Bifurcation of limit cycles of system (C.1) | a = b = Consider system (C.1) with a = b = Û x = − y + x ( cx + dxy ) , Û y = x + y ( cx + dxy ) . (C.6)First, note that we may use a spatial scaling x → x / √ d , y → y / √ d ( d ,
0) to system (C.6) to obtain Û x = − y + x ( αx + xy ) , Û y = x + y ( αx + xy ) , (C.7)where α = c / d . We now consider the following class of perturba-tion of system (C.7): Û x = − y + x ( αx + xy ) + Õ k = ε k (cid:16) λ k x + Õ i + j = a k , i , j x i y j (cid:17) , Û y = x + y ( αx + xy ) + Õ k = ε k (cid:16) λ k y + Õ i + j = b k , i , j x i y j (cid:17) , (C.8)where λ k , a k , i , j ’s and a k , i , j ’s are real constants. Here we study thenumber of limit cycles of system (C.8) by the second order averag-ing method. We remark that special case of system (C.8) ( α = λ =
0) has been studied by the first order averaging method (see[49]), In fact, the system there can be written as the form (C.7) | α = by doing a linear change of variables.Now applying our algorithm NormalForm by taking k = drdθ = F ( θ , r ) + εF ( θ , r ) + ε F ( θ , r ) + O( ε ) , (C.9)where F ( θ , r ) = r cos θ (cid:0) ( α − ) cos θ + (cid:1) , and F ( θ , r ) = S (cid:16) r ( α − )( a , , − a , , + a , , + b , , − b , , ) C − r (cid:0) α ( a , , − a , , + b , , ) − a , , + a , , − a , , − b , , + b , , (cid:1) C − r (cid:0) a , , − a , , − b , , + b , , − b , , (cid:1) C + r (cid:0) αa , , − a , , + a , , − b , , (cid:1) C + r (cid:0) a , , − b , , + b , , (cid:1) C + r a , , C + r b , , (cid:17) + r ( α − )( a , , − a , , − b , , + b , , − b , , ) C − r ( α ( a , , − a , , − b , , + b , , ) − a , , + a , , + b , , − b , , + b , , ) C + r (cid:0) a , , − a , , + a , , + b , , − b , , (cid:1) C + r (cid:0) α ( a , , − b , , )− a , , + a , , + b , , − b , , (cid:1) C − r (cid:0) a , , − a , , + b , , − b , , (cid:1) C + r (cid:0) a , , − b , , (cid:1) C + r (cid:0) a , , + b , , (cid:1) C + rλ with S = sin θ , C = cos θ , the expression of F ( θ , r ) is quite long sowe omit it here. Using the algorithm DSolutions we obtain a list[ r ( θ , z ) , SIs , Y ( θ , z ) ], where r ( θ , z ) = z [ − ( α − ) z cos θ sin θ − ( α + ) z sin θ ] / , S Is = n < z , < z ( − ( α + ) z ) / , < z ( + ( α + ) z ) / o , Y ( θ , z ) = [ − ( α − ) z cos θ sin θ − ( α + ) z sin θ ] / . To obtain the interval D in this case, we construct an equivalentsolution set SIs : = { < z , < z −( α + ) z , < z + ( α + ) z } . By us-ing the Maple command SolveTools[SemiAlgebraic] to SIs , weobtain the interval D as follows: D = n < z , z < − ( α + ) / o , if α < − / , { < z } , if α = − / , n < z , z < ( α + ) / o , if α > − / . Now assuming the case α > − / Av-eragedFunction for the system (C.8) when k =
1, we find thatMaple may not evaluate the following kind of parametric integrals: I i , j = ∫ π (cid:18) cos i θ sin j θ − ( α − ) z cos θ sin θ − ( α + ) z sin θ (cid:19) dθ , i , j ∈ N . In view of this, from now on we restrict to the case α =
1. Comput-ing the first averaged function by the algorithm
AveragedFunc-tion and taking the transformation z = h − s ( + s ) i / for 0 < s < f ( z ) z = h − s ( + s ) i / = ⇒ ¯ f ( s ) = π / ( − s ) / ( + s ) / ( + s ) / (cid:0) N s + N s + N s + N s + N s + N s + N (cid:1) , (C.10)where N = a , , + a , , − b , , − b , , − b , , + λ , N = − a , , + b , , + a , , − b , , − b , , + λ , N = a , , − a , , − b , , + b , , − b , , + λ , N = − a , , + a , , + b , , − b , , + b , , + λ . As a result of the symmetry of coefficients of the function ¯ f ( s ) with respect to s , we know that if s , f ( s ) , so is1 / s . Hence, the fact that ¯ f ( s ) has at most three zeros in s ∈ ( , ) implies that there exist at most three zeros for f ( z ) in z ∈ ( , / √ ) .Moreover, the number three can be reached since the constants N , N , N and N are independent. It follows from Theorem 2.3 thatthe first order averaging provides the existence of at most threelimit cycles of system (C.8), and this number can be reached. Hence,with the nonzero constant λ the perturbed system (C.8) can pro-duce one more limit cycle than the case without it.In order to consider the second order averaging, we must usethe conditions of f ( z ) =
0. The following lemma follows fromequation (C.10).
Lemma C.3.
The first averaged function f ( z ) ≡ if and only if a , , = b , , , a , , = b , , , and b , , = λ = . Now we consider the second order averaging of system (C.8)with α =
1. First, we update the obtained dr / dθ in (C.9) by Lemma .3. We put the resulting expression of y ( θ , z ) as follows. y ( θ , z ) = z ( z sin θ − ) / (cid:16) m · ln ( − z sin θ ) + m · cos θ sin θ + m · cos θ sin θ + m · cos θ + m · cos θ + m · cos θ − b , , z cos θ + m · sin θ + m (cid:17) , (C.11)where m = (cid:0) − a , , + b , , (cid:1) z + (cid:0) − a , , + a , , + b , , − b , , (cid:1) z + ( a , , − a , , − a , , − b , , + b , , ) , m = − z (cid:16) (cid:0) a , , − a , , − a , , + b , , − b , , (cid:1) z + ( a , , − a , , − a , , − b , , + b , , ) (cid:17) , m = z ( b , , − a , , − b , , + a , , + a , , ) , m = − z ( b , , − a , , − b , , + a , , + a , , ) , m = z ( b , , − b , , ) , m = z (cid:16) (cid:0) a , , − a , , + b , , (cid:1) z + b , , − a , , − b , , + a , , + a , , (cid:17) , m = z (cid:16) ( a , , + a , , + a , , + b , , + b , , ) z − ( a , , + a , , − a , , − b , , + b , , ) z + ( a , , − a , , − a , , − b , , + b , , ) (cid:17) , m = z (cid:16) (cid:0) b , , + b , , (cid:1) z + (cid:0) − b , , − a , , + b , , − a , , + a , , (cid:1) z + ( b , , − b , , + a , , − a , , − a , , ) (cid:17) . Applying the algorithm
AveragedFunction to the updated nor-mal form dr / dθ when k =
2, we find that Maple was consuming toomuch of the CPU during a calculation, and can not able to allocateenough memory. This is mainly because the expressions involvedin the analysis are huge, and sometimes Maple can not evaluatethe integrals for certain complicated functions. For instance, in ourcase here, Maple can not evaluate the integrals in z ∈ ( , / √ ) of the following form: J i = ∫ π ln (cid:0) − z sin θ (cid:1) cos i θ z ( cos θ − ) + z sin θ − ! dθ , i ∈ N . To simplify our computation, we let the term ln (cid:0) − z sin θ (cid:1) in y ( θ , z ) be identically zero (i.e., m ≡ Corollary C.4.
The term ln (cid:0) − z sin θ (cid:1) in y ( θ , z ) vanishesif and only if a , , = b , , , a , , = b , , , and a , , = . Applying Lemma C.3 and Corollary C.4 to system (C.8), we ob-tain the following second order averaged function by the algorithm
AveragedFunction . f ( z ) = − π z √ − z (cid:16) ( h z + h z + h ) p − z − ( − z ) (cid:0) b , , z + h z + h (cid:1)(cid:17) , where h = a , , + a , , − b , , − b , , + b , , − λ , h = a , , − a , , − b , , + b , , − b , , , h = − a , , + a , , + b , , − b , , + b , , , h = − a , , + b , , − b , , . After making the transformation as before, f ( z ) becomes f ( z ) = π / ( − s ) / ( + s ) / ( + s ) / (cid:0) H s + H s + H s + H s + H s + H s + H (cid:1) , (C.12)where H = a , , + a , , − b , , − b , , − b , , + λ , H = − a , , + b , , + a , , − b , , − b , , + λ , H = a , , − a , , − b , , + b , , − b , , + λ , H = − a , , + a , , + b , , − b , , + b , , + λ . It follows from Theorem 2.3 and equation (C.12) that the averagingmethod up to second order provides the existence of at most threelimit cycles of system (C.8), and this number can be reached byusing similar arguments to the function f ( z ) .We summarize our results as follows based on the above analy-sis. Theorem C.5.
For sufficiently small parameter | ε | > , the per-turbed system (C.8) with α = has at most three limit cycles bifur-cating from periodic orbits of the unperturbed one by the first orderaveraging method, and by the second order averaging method underthe condition: a , , = b , , , a , , = b , , , and a , , = . In eachcase this number can be reached.. In eachcase this number can be reached.