Center and isochronous center of a class of quasi-analytic switching systems
Feng Li Pei Yu Yirong Liu Yuanyuan Liu
























 









^aSchool of Science, Linyi University, Linyi, Shandong 276005,
P.R. China 
 
 
^bDepartment of Applied Mathematics, Western University, 

London, Ontario N6A 5B7, Canada 
 
^cSchool of Mathematics and Statistics, Central South University, 

Changsha, Hunan 410012, P.R. China 
 



In this paper, we study the integrability and linearization of a class
of quadratic quasi-analytic 
switching systems. We improve an existing
method to compute the focus values and periodic constants of
quasi-analytic switching systems. In particular, with our 
method, we demonstrate that the dynamical behaviors of quasi-analytic
switching systems are more complex than that of continuous quasi-analytic
systems, by showing the existence of six and seven limit 
cycles in the neighborhood of the origin and infinity, respectively, 
in a quadratic quasi-analytic switching system.
Moreover, explicit conditions are obtained for classifying the
centers and isochronous centers of the system.



Quasi-analytic switching systems, Lyapunov constant,
limit cycle, center, isochronous center.


 
MSC(2010) 34C07, 34C23













§ INTRODUCTION


The problem of characterizing the centers and isochronous centers of dynamical
systems has attracted the attention
of many researchers. So far, regarding the family of polynomial differential 
systems, a complete classification of centers and isochronous centers has 
only been solved for quadratic polynomial systems, or simply quadratic 
systems. Quadratic systems having a center were classified by 
Dulac <cit.>, Kapteyn <cit.>,
Bautin <cit.>, Żoła̧dek <cit.>, 
Yu and Han <cit.>, while quadratic 
systems having an isochronous center were characterized by 
Loud <cit.>. The centers of the cubic systems with homogeneous 
nonlinearities were studied in <cit.>,
and the isochronous centers for such cubic systems were further 
investigated by Pleshkan <cit.>. However it is still 
far away from obtaining a complete classification of the centers and 
isochronous centers for polynomial differential systems of degree three,
and it is extremely difficultly to study these problems when the degree 
of the systems is increased. For example, consider the following systems,

    ż=(λ+ i)z
    + ( zz̅)^d-5/2 (Az^4+jz̅^1-j+ Bz^3z̅^2 + Cz^2- jz̅^3+j+ Dz̅^5), 
      d = 2m + 1 ≥ 5,
    
          ż= iz+(zz̅)^d-4/2(Az^3z̅
    +Bz^2z̅^2+Cz̅^4),  d=2m ≥ 4,
    
          ż= (λ+i)z+(zz̅)^d-3/2(Az^3
    +Bz^2z̅+Czz̅^2+Dz̅^3),    d=2m+1≥3,
    
          ż= (λ+i)z+(zz̅)^d-2/2
    (Az^2+Bzz̅+Cz̅^2),  d=2m ≥ 2,

which have been investigated by Llibre and Valls
<cit.>,
and the conditions on centers and isochronous centers were obtained.
However, in these articles, the parameter d was restricted to
make the systems be polynomial systems.

Recently, switching systems have been widely used in modelling many 
practical problems in science and engineering. Theory suggests that 
switching systems can be considered as a uniform limit of continuous systems, 
and that the global dynamics of continuous models may be approximated 
by switching systems. In fact, the richness of dynamical behavior found in
switching systems covers almost all the phenomena discussed in general 
continuous systems, such as limit cycles, homoclinic and heteroclinic orbits, 
strange attractors.  For example, Leine and Nijmeijer <cit.>, 
and Zou and Kapper <cit.> considered non-smooth Hopf bifurcation 
in switching systems. Bifurcation of limit cycles from centers of 
discontinuous quadratic systems was studied by 
Chen and Du <cit.>. Limit cycles in a class of 
continuous and discontinuous cubic polynomial differential systems
were investigated by Llibre et al. <cit.>.
Bifurcation of limit cycles in discontinuous quadratic differential systems 
with two zones was considered in <cit.>.
The Melnikov function method has also been extended to study homoclinic 
bifurcation of non-smooth systems <cit.>. 
In addition, some general efficient methods have also been developed 
to study non-smooth 
systems. Among these methods, normal form computation for impact 
oscillators was given in <cit.>, and a general methodology
for reducing multidimensional flows to low dimensional maps in piecewise
nonlinear oscillators was proposed in <cit.>. 
Center and isochronous center conditions for switching systems associated 
with elementary singular points were discussed in <cit.>.

More recently, quasi-analytic systems have also been widely used in modelling
many practical problems. By “quasi-analytic”, we mean that 
the system may be analytic for some parameters but not for some other 
parameters. For example, an axis-symmetric quasi-analytical model was 
developed in order to simulate the behavior of a RFEC system during its 
operation <cit.>. A quasi-analytical model  for scattering 
infrared near-field microscopy has been designed for predicting and 
analyzing signals on layered samples <cit.>.
A simple quasi-analytical model was developed in <cit.> 
to study the response of ice-sheets to climate. 
On the other hand, a general 
type of quasi-analytic systems, described by 
 

    [ ẋ=δ x-y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2X_k(x,y),; ẏ=x+δ y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2Y_k(x,y), ]



 
where

    X_k(x,y)= ∑_α+β=kA_αβx^α y^β,     
    Y_k(x,y)= ∑_α+β=kB_αβx^α y^β,

has been studied by Liu  <cit.> and Liu et al. <cit.>. 
As special cases, quadratic quasi-analytic 
systems have been studied in <cit.> and 
cubic quasi-analytic systems in <cit.>.
In particular, generalized focal values and bifurcation of limit cycles 
for quadratic quasi-analytic systems were discussed in <cit.>. 
Here, a quadratic quasi-analytic system is defined by 
taking k=2 only in (<ref>). Similarly, cubit, etc. quasi-analytic 
systems can be defined.  

Similar to quasi-analytic continuous systems, in this paper, we propose 
to study the center and isochronous center conditions for the following class
of discontinuous planar systems: 


    [ ẋ=δ x-y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2F_k^+(x,y),; ẏ=x+δ y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2G_k^+(x,y), ]   (y>0), 
    [ ẋ= δ x-y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2F_k^-(x,y),;  ẏ=x+δ y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2G_k^-(x,y), ]   (y<0),


for (x,y) ∈ R^2, 
where the two subsystems (<ref>) 
and (<ref>) describe 
dynamics on the upper and lower half planes, called 
the upper and lower systems, respectively.  
For λ>0 (< 0), the linear terms in (<ref>) are the lowest 
(highest) order terms in these functions. Hence, when λ>0,
the origin of (<ref>) is a center or a focus. When λ<0, 
system (<ref>) has no real singular point in the equator of 
Poincaré compactification, but the point at infinity is a center 
or a focus. Therefore, it is necessary to determine, for λ≠ 0,
whether or not the origin (when λ >0) or infinity 
(when λ < 0) is a center (or a weak focus ).
As a continuous work of <cit.> and <cit.>, we 
generalize the study to consider center and isochronous center conditions 
of quasi-analytic switching systems and hope to promote the research 
in this direction. 

As compared with bifurcations, the center and isochronous center problems 
of switching systems have not received much attention. They should be 
considered carefully for switching systems because they are closely related 
to conditions on integrability and linearization.
Freire et al. <cit.> discussed the center problem in a 
piecewise linear system. But it is much hard to solve center or isochronous 
center problem for nonlinear or piecewise linear systems because  the 
classical methods for computing lyapunov constants and periodic
constants are no longer applicable. Thus, new techniques are needed to 
develop.

In this paper, we study quasi-analytic switching systems mainly from 
two aspects. First of all, we modify and improve existing methods to 
compute Lyapunov constants and periodic constants for quasi-analytic 
switching systems. Secondly, as an application, we study a quadratic 
quasi-analytic switching system and derive its center and isochronous center 
conditions.

The rest of the paper is organized as follows.  In Section 2, 
we present a method to compute the return map of system (<ref>). 
As an application, a class of quadratic quasi-analytic
switching systems is studied in Section 3, and the center and 
isochronous center are classified by using our method. Finally, 
our conclusion is drawn in Section 4.




§ LYAPUNOV CONSTANTS OF THE QUASI-ANALYTIC 
SWITCHING SYSTEM (<REF>)


Under the transformation of the polar coordinates,

    x=r^1/λcosθ,   y=r^1/λsinθ,

system (<ref>) becomes 


    ṙ   =λ r (δ+∑_k=1^∞φ_k+2(θ)r^k ),
    θ̇   =1+∑_k=1^∞ψ_k+2(θ)r^k,
         (r ≥ 0)

where φ_k(θ) and ψ_k(θ) are polynomial
functions in cosθ and sinθ, given in the form of

    φ_k(θ)=cosθ X_k-1(cosθ, sinθ)
    +sinθ Y_k-1(cosθ, sinθ),


    ψ_k(θ)= cosθ Y_k-1(cosθ, sinθ)
    -sinθ X_k-1(cosθ, sinθ).

Then, it follows from system (<ref>) that

    dr/dθ=λ  r δ+Σ_k=1^∞φ_k+2(θ)r^k/1+Σ_k=1^∞ψ_k+2(θ)r^k.

Obviously, the polar coordinate form of the quasi-analytic system 
(<ref>) differs from
 that of analytic systems by only a constant factor λ.
It is easy to see that (<ref>) is a special case of the
following equation,

    dr/dθ= r ∑_k=1^∞ R_k(θ)r^k,    (r ≥ 0).


By the method of small parameters of Poincaré,
the general solution of (<ref>) can be expressed as <cit.>

    r=r̃(θ,h)=∑_k=1^∞ v_k(θ)h^k,

where v_1(0)=1,  v_k(0)=0, ∀ k≥ 2.
Now, substituting the above solution r=r̃(θ,h)
into (<ref>) yields

    v'_1(θ)   =R_0(θ)v_1(θ),
    
    v'_2(θ)   =R_0(θ)v_2(θ)+R_1(θ)v_1(θ)^2,
       ⋮
    
    v'_m(θ)   =R_0(θ)Ω_1,m(θ)+R_1(θ)
    Ω_2,m(θ)+⋯+R_m-1(θ)Ω_m,m(θ).

Thus, we may solve v_k(θ) one by one to obtain

    v_1(θ)   =e^∮_0^ϑ R_0(φ)dφ,
    
    v_2(θ)   =2v_1(θ) ∮_0^ϑ R_1(φ)v_1(φ)dφ,
       ⋮
    
    v_m(θ)   =v_1(θ) ∮_0^ϑR_1(φ)
    Ω_2,m(φ)+⋯
    +R_m-1(φ)Ω_m,m(φ)/v_1(φ)  dφ.


Note that R_0(θ)=λδ for system (<ref>).
Further, we define the successive function as

    Δ(h)=r̃(2π,h)-h,

and thus the critical point being a center must satisfy Δ(h)=0, 
namely,

    r̃(2π,h)=h.

Many methods have been developed to compute the successive function
Δ (h) (e.g.,see <cit.>).

From the second equation of (<ref>), we can also obtain

    t=T(θ,h)=∫_0^θdϑ/1+∑_k=1^∞ψ_k+2(θ)r̃(θ,h)^k,

which implies that the critical point being an isochronous center 
should satisfy

    r̃(2π,h)= h    and   T(2π,h)= 2π.


However, the classical methods and formulas cannot
be directly applied to a non-analytic switching system
due to discontinuity. We need to modify the existing methods
to resolve this problem. 
Similar to the return map defined for analytic switching systems
<cit.>, the approach used in <cit.> (see Lemma 2.1) can be 
extended to define the return maps for the 
quasi-analytic switching system (<ref>). The basic idea is briefly
illustrated as follows, see Figure <ref>. First of all, we define 
the positive half-return map of the upper phase
of system (<ref>). Then, by a transformation y→ -y, 
the lower half phase could be transformed into the upper phase, 
as shown in Figure <ref>. Further, using a time reverse changing, the 
computation of this transformed half-return map of the lower phase is 
replaced by computing the positive half-return map of the following system,

    [ ẋ= δ x-y- ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2F_k^-(x,-y),;  ẏ=x+δ y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2G_k^-(x,-y), ]   (y>0),

which is shown in Figure <ref>. Therefore, we only need to compute 
the two positive half-return maps
for systems (<ref>) and (<ref>).



By defining the successive functions for systems (<ref>)
and (<ref>), respectively, as

    Δ_1(h)=r̃_1(π,h)-h
       and  Δ_2(h)=r̃_2(π,h)-h,

then we obtain the successive function for the switching
system (<ref>), defined as

    Δ(h)=Δ_1(h)-Δ_2(h)=r̃_1(π,h)-r̃_2(π,h).


Similarly, the period constants for systems (<ref>)
and (<ref>) can be defined as

    [ T_1(θ,h)= ∫_0^πdϑ/1+∑_k=1^∞ψ_2+k(θ)r̃_̃1̃^k(ϑ,h),; T_2(θ,h)= ∫_0^πdϑ/1+∑_k=1^∞ψ_2+k(θ)r̃_̃2̃^k(ϑ,h), ]

which in turn yield the period function for the switching
system (<ref>) in the form of

    T=T_1(π,h)+T_2(π,h)=2π+∑_k=1^∞ T_k h^k.


In particular, if the equations describing the lower half plane are given by

    [ ẋ= δ x -y+ ∑_k=2^∞
    (x^2+y^2)^(k-1)(λ-1)/2 F_k^-(x,y)=-y,;   ẏ=x+δ y+ ∑_k=2^∞ (x^2+y^2)^(k-1)
    (λ-1)/2 G_k^-(x,y)=x. ]

then we only need to compute Δ_1(h) and T_1(θ,h).

Based on the above results, we can define the focus values and
periodic constants for the quasi-analytic switching system (<ref>).


Δ(h) can be written as

    Δ(h)=∑_1^n [ u_1(π)-v_1(π) ] h^k=∑_k=1^∞ 
    V_k h^k,

where V_k is called the kth-order focus value  at the origin 
(or infinity) of the quasi-analytic switching system (<ref>).




T(h) can be expressed as

    T(h)=T_1(π,h)+T_2(π,h)=2π+∑_k=1^∞ T_k h^k,

where T_k is called the kth periodic constant at the origin (or infinity) 
of the quasi-analytic switching system (<ref>).


Having defined V_k and T_k, we now describe the steps in computing them.


  * Introduce the transformations: y→ -y and t→ -t
for the lower half plane.


  * By using the transformation of polar coordinates,

    x=r^1/λcosθ,   y=r^1/λsinθ,

in the systems (<ref>) and (<ref>), write the solutions 
for the systems (2.8) and (2.10) as

    r_1=r̃_1(θ,h)=Σ_k=1^∞ u_k(θ)h^k
       and  
    r_2=r̃_2(θ,h)=Σ_k=1^∞ v_k(θ)h^k,

respectively, satisfying u_1(0)=v_1(0)=1, u_k(0)=v_k(0)=0, ∀ k≥ 2.


  * Solve u_k(θ) and v_k(θ).


  * Compute the successive function for the switching system using the formula,

    Δ(h)=Δ_1(h)-Δ_2(h)=r̃_1(π,h)-r̃_2(π,h).



  * Compute the periodic constants for the switching system using the formula,

    T=T_1(π,h)+T_2(π,h).



Obviously, the symmetry principle for continuous systems is no longer
applicable for switching systems.
We need to redefine symmetry for switching systems
in order to derive the center conditions of switching systems.


If both systems (<ref>) and  (<ref>)
are symmetric with respect to the y-axis,
then system (<ref>) is said to be symmetric with respect to the y-axis.
Further, if the vector fields of systems (<ref>) and
(<ref>) satisfy

    F_k^+(x,y)=-F_k^-(x,-y)    and   G_k^+(x,y)= G_k^-(x,-y),

then system (<ref>) is said to be symmetric with respect to the x-axis.

So obviously, if system (<ref>) is symmetric with respect to 
the x-axis or the y-axis, then the origin of system (<ref>) 
is a center.


























§ A QUADRATIC QUASI-ANALYTIC SWITCHING SYSTEM


In this section, we consider a quadratic quasi-analytic switching system
to demonstrate the application of the formulae and results obtained in the
previous section. We will use our method to determine the center conditions
and isochronous center conditions of the system we will consider, given by 

    [ [ ẋ= δ x -y 
    +(x^2+y^2)^(λ-1)/2(a_20x^2+a_11xy+a_02y^2),; ẏ=x + δ y 
    +(x^2+y^2)^(λ-1)/2(b_20x^2+b_11xy+b_02y^2), ]                                                                                                             (y>0),;                                                                               [3.5ex]
    [  ẋ= δ x -y,; ẏ=x + δ y , ]                                                                                                             (y<0). ]

The case λ=1 (a polynomial system) has been studied 
in <cit.>, which becomes a special quadratic switching system.
It is shown in <cit.> that the highest order of focus values 
for this special system is 5, and 5 small-amplitude limit cycles 
are obtained.  
We want to extend the study to the case λ 1. 
However, when the lower system is not in a simple form, 
even for general quadratic switching systems,  
it is very difficult to determine the center conditions
and isochronous center conditions. Thus, in this paper we 
focus on the study of system (<ref>) for λ 1. 




 §.§ Center conditions and limit cycles for system (<ref>)


We first study the center conditions and bifurcation of limit cycles 
in system (<ref>). It has been recently noticed that Tian and Yu 
studied a quadratic switching Bautin system and obtained 10 
small-amplitude limit cycles <cit.>. Here, we want to show 
that system  (<ref>) with λ 1 can bifurcate 7 
limit cycles around the origin, 
two more than that of the system with λ =1. 

In order to consider the center and isochronous center conditions, 
and determine the number of limit cycles bifurcating in the small
neighborhood of the origin (or infinity), we need to compute the 
Lyapunov constants and periodic constants.  
With the aid of a computer algebra system –
Mathematica, we obtain the following Lyapunov constants of system 
(<ref>). 


The first three Lyapunov constants at the
origin are given by

    L_0     = 2 πδ, 
     
    L_1   =-2/3(a_11 + 2 b_02 + b_20)λ,
    
    L_2   =-π/8[ b_20 (a_20 + a_02) 
    + ( 2 a_20 + b_11) ( b_20 + b_02) ] λ .
 
For higher Lyapunov constants, we have two cases. 


 
Case (A) b_20≠ 0. For this case, L_3 is given by 

    L_3 =   -2/105λ{[ 6 a_20b_02 ( 2 a_20+ b_11) 
    + b_20 ( 3 a_20^2 + 4 b_02^2 ) ] (λ+6)  
                + 14 b_20[ 3 a_20 ( 2 a_20+b_11) 
    + 2 b_20 b_02] }.

Then, there are three sub-cases in computing L_i,   i≥ 4.


 
Case (A1)  
a_20 [  7 b_20 + b_02 (λ+6) ] ≠ 0, 
for which we have 

    L_4=    a_20 b_20λ π/1536 
    [ 7 b_20 + b_02 (λ+6)]^2[4 b_20+ b_02(λ+3) ] 
    [ 12 b_20 + b_02 ( λ + 9) ] 
       ×{ 3 a_20^2 ( λ+6)^2 
    - 28 b_20[ 7 b_20 + b_02 (λ+6) ] }.

(a)  If b_20 = -1/4b_02 (λ+3), then

    L_5   =- b_02 λ (λ+3) /3243240
    (λ-1)^2  f_1,
    
    L_6   =- π  b_02 λ^3 (λ+3)^2 /5308416 a_20
    (λ-1)^3  f_2,
    
    L_7   =- b_02 λ (λ+3)/6788231049600 a_20^2
    (λ-1)^4  f_3,
 
where

    f_1=    
    (3868 a_20^4 - 2909 a_20^2 b_02^2 - 1085 b_02^4 ) λ^5 
    + (34320 a_20^4 - 5962 a_20^2 b_02^2 - 
    1085 b_02^4 ) λ^4  
        + 2 (41274 a_20^4 + 3162 a_20^2 b_02^2 + 3563 b_02^4) λ^3 
    - 6 b_02^2 (2937 a_20^2 + 259 b_02^2) λ^2 
        + 243 b_02^2 (83 a_20^2 - 35 b_02^2) λ +5103 b_02^4 ,
    
    f_2=    [ (4 a_20^2 -7 b_02^2) λ^2 
    + 2 (24 a_20^2 - 7 b_02^2 ) λ + 3( 48 a_20^2 + 7 b_02^2) 
    ]
       ×[ 
    (196 a_20^4 + 25 a_20^2 b_02^2 + 4 b_02^4 ) λ^2 
    - 2 b_02^2 (53 a_20^2 + 4 b_02^2) λ+ 81 a_20^2 b_02^2 + 4 b_02^4 ],
    
    f_3=    56 b_02^8 (λ-1)^4 (λ+3)
    ( 1027505 λ^4 + 426036 λ^3 
    - 1498554 λ^2 
        +551124  λ +1476225 ) 
         +128 a_20^8 λ^5 ( 
    9756302 λ^4  
    + 203092731 λ^3 
    + 1398704409 λ^2  
        + 3748721013 λ +2969005833 ) 
        -4 a_20^6 b_02^2 λ^3 ( λ-1) 
    ( 488213167 λ^5 + 3344141799 λ^4 + 596192742 λ^3 
        
    - 12028807314 λ^2 
    + 57678458859 λ 
    + 102228640299 )


    +a_20^2 b_02^6 (λ- 1)^3 
    ( 24084923 λ^6
    - 1850127480 λ^5 
    - 6861739149 λ^4 
            
    + 1283638536 λ^3 
    + 11199899889 λ^2 
    - 5507145936 λ 
    -10721822175 
    ) 
        -3 a_20^4 b_02^4 (λ- 1)^2 λ 
    ( 192786817 λ^6
    + 2736149472 λ^5 
    + 4596112377 λ^4 
          
    - 12782671128 λ^3 
    - 13724853597 λ^2 
    - 2481265224 λ
    -22052498589).



(b)  If b_20 = -1/12b_02 (λ+9), then

    L_5   =- b_02λ(λ+9)/4169880(5λ+9)^3 f̃_1,
    
    L_6   =- π  b_02 λ^3(λ+9)/143327232a_20(5λ+9)^4 f̃_2,
    
    L_7   =- b_02 λ(λ+9) /183282238339200a_20^2
    (5  λ+9)^5 f̃_3,

where 

    f̃_1=   
    (43092 a_20^4 + 52305 a_20^2 b_02^2 + 8675 b_02^4 ) λ^6 
    
        + (386532 a_20^4 +597639 a_20^2 b_02^2 + 86980 b_02^4 ) λ^5 
        + 9 (45900 a_20^4 +56538 a_20^2 b_02^2 - 8377 b_02^4 ) λ^4 
        - 162 (5202 a_20^4 + 68961 a_20^2 b_02^2 + 11720 b_02^4) λ^3 
        - 729 b_02^2 (47283 a_20^2 + 6887 b_02^2) λ^2 
        - 177147 b_02^2 ( 159 a_20^2 + 28 b_02^2 ) λ 
      - 1594323 b_02^4 ,
    f̃_2=    [ (108 a_20^2 + 35 b_02^2 ) λ^2
    + 54 (24 a_20^2 + 7 b_02^2) λ 
    + 3888 a_20^2 + 567 b_02^2 ]
       ×[
    (11772 a_20^4 + 15255 a_20^2 b_02^2 + 2300 b_02^4) λ^4 
          
    + 24 (1035 a_20^4 + 2256 a_20^2 b_02^2 + 220 b_02^4) λ^3  
          
    - 18 (2034 a_20^4 + 13479 a_20^2 b_02^2 + 1036 b_02^4) λ^2 
          - 2592 b_02^2 (417 a_20^2 + 25 b_02^2 ) λ 
    - 12393 b_02^2 (81 a_20^2 +4 b_02^2) ],
    f̃_3=    279936 a_20^8 (λ-1) λ^5
    ( 38218052 λ^4 
    + 790049163  λ^3 
    +  5940016749  λ^2 
        
    + 19219084137  λ
    +22987749699 
    ) 
        +56 b_02^8 (λ+9) (5 λ+9)^4 
    ( 
     4362337 λ^5
    + 301281  λ^4 
    - 20129886  λ^3 
         
    - 60022134  λ^2 
    - 143311923  λ 
    -39858075  
    ) 
        +324 a_20^6 b_02^2 λ^3 (5 λ+9) 
    ( 
     12357727123 λ^6
    + 228212204300 λ^5 
                
     + 1303142758287 λ^4 
    +688198402656 λ^3 
    - 17536030592607 λ^2 
                
    - 57374137182060 λ 
    -55078300245699 
    ) 
        +3 a_20^2 b_02^6 (5 λ+9)^3 
    ( 
     5875349809 λ^7
    + 75871528713 λ^6 
    +173810927169 λ^5 
                 
    - 698576296239 λ^4 
    - 3526727842629 λ^3 
    - 8234719793325 λ^2


    - 10205502265773 λ 
    -2605402788525 ) 
        +27 a_20^4 b_02^4 λ (5 λ+ 9)^2 
    ( 
     15822388201 λ^7
    + 255144264153 λ^6 
                 
    + 1105548667713 λ^5
    - 1017555100407 λ^4 
    - 18518298666957 λ^3 
                 
    - 51274399423149 λ^2 
    - 74254559009373 λ 
    - 56058065936181).



(c)  If a_20^2 = 28 b_20 [7 b_20 + b_02(λ+6)]/3 (λ+6)^2, the following holds: 

    L_5=    128 b_20^5/715 (λ+6)^4 λ (λ-1) (2 λ-9) (5 λ+16),
    
    L_6=      0,
    
    L_7=     -128b_20^7/692835 (λ+6)^6 λ(λ-1)(2λ-9) 
       × (1587λ^3 -82636λ^2 -527988λ -838080).



Case (A2)  
a_20=0, under which L_3 is reduced to  

    L_3= - 8  b_20 b_02/105 λ[ 7 b_20 + b_02(λ+6 )],

and higher Lyapunov constants are  

    L_4=   -π b_02^3 b_11/3136 λ (λ-1) (5 λ+9) ,
    
    L_5=   -128 b_02^3/108153045 (λ+6) λ (λ-1) (2 λ-9) 
         ×[ 9 b_02^2 (5 λ+16)(λ+6)^2 
    +49 b_11^2 (27 λ+50) ],
    
    L_6=   -π b_02^3 b_11/59006976 (λ+6) λ (λ-1) 
    [ 3 b_02^2 (λ+6)^2 (562 λ^2 + 5085 λ +13365 ) 
        - b_11^2 ( 60074 λ^2 - 6615 λ -235935 ) ],
    
    L_7=   128 b_02^3 /231084662834025 (λ+6)^3λ ( λ-1) 
       ×[ 405 b_02^4 (2 λ -9)(λ+6)^4 
    (1587λ^3 -82636λ^2 -527988λ -838080) 
        + 441 b_02^2 b_11^2 (λ+6)^2  
    (2699488 λ^4 +1664883 λ^3 -12695418 λ^2 
        +64285812 λ +193185000 ) 
        + b_11^4 (17594187058 λ^4 + 24560280388 λ^3 
    -109958232048 λ^2 
            
    - 10609846128 λ +343997992800 ) ].


Case (A3)   
b_20=-1/7b_02(6 + λ), for which L_3 becomes 

    L_3= 2a_20^2 b_02/245 λ(λ+6)^2,
 
and higher Lyapunov constants are given by 

    L_4=   -π  b_02λ (λ-1) (9 + 5 λ) /43904 (6 + λ) [ a_20 b_02^2 λ^2 
    + 2 b_20^2 (20 a_20 + 7 b_11) λ
            +12 a_20 ( 17 b_02^2- 49 a_20^2 - 49 a_20 b_11 )
    + 21 b_11 (4 b_02^2 -7 a_20 b_11)  
    ],
    
    L_5=   18 a_20^4 b_02/105105 λ (2 λ + 9) ( 7 λ+12) (λ + 2) (λ+6) ,
    
    L_6=   π  a_20 b_02^3 λ (λ+6)/4956585984 [ 2 ( 34496 a_20^2 + 15375 b_02^2) λ^4 
            + 3 ( 29792 a_20^2 + 165045 b_02^2) λ^3  
    - 18( 279104 a_20^2 + 37917 b_02^2) λ^2 
            - 27 (130144 a_20^2 + 103149 b_02^2) λ  
    +1620 ( 1568 a_20^2 + 891 b_02^2)  ],
    
    L_7=   2a_20^2 b_02/245 λ(λ+6)^2.


Case (B)  b_20=0. For this case, we have 

    L_1=    -2(a_11 + 2 b_02)/3 λ, 
     
    L_2=    -π b_02(2a_20+ b_11)/8 λ,
    
    L_3=   -2 b_02(2a_20 + b_11)/315  
    (8 a_02 - 9 b_11) λ (λ + 6), 
      
    L_4=   -π b_02(2a_20+ b_11)/9216  
    (28 a_02^2 + 36 b_02^2 - 36 a_02 b_11+ 27 b_11^2)
     λ (λ+3) (λ+9),
    
    L_5=    b_02(2a_20+ b_11)/1216215 [ 27 b_11( 27 a_02^2 + 120 b_02^2 
    - 60 a_02 b_11+ 35 b_11^2) 
                    - 128 a_02 (10 a_02^2 + 27 b_02^2) ]
     λ (λ+2) (λ+12) (2 λ+9) ,
    
    L_6=   -π b_02(2a_20+ b_11)/28311552[ 16 (143 a_02^4 + 594 a_02^2 b_02^2 + 243 b_02^4) 
        - 288 a_02 b_11( 11 a_02^2 + 45 b_02^2) 
    + 1080 b_11^2 (3 a_02^2 + 7 b_02^2 
        - 45 b_11^3 ( 56 a_02- 27 b_11) ]
     λ (λ+3) (λ+6) (λ+15) (2 λ+3).

Note that in the above computations, 
L_k-1=0, k=1,2, ⋯,6 have been used in computing L_k. 


Now, by carefully analyzing the above Lyapunov constants, we obtain 
the following result. 

 
For system (<ref>), maximal six small-amplitude 
limit cycles can bifurcate from the origin and maximal 
seven small-amplitude 
limit cycles can exist in the neighborhood of infinity. 
Moreover, the first seven Lyapunov constants at the origin (or infinity) 
of system (<ref>) vanish if and only if one of the following 
conditions is satisfied: 

(i)  δ = a_11=b_20=b_02=0, 

(ii)  
δ = b_20= a_11 +2b_02 = b_11 +2a_20=0, 

(iii)  
δ = a_20=b_02= a_11 +b_20= b_11 +a_02 = 0, 

(iv)  
δ = λ-1 = a_02=a_20= b_02 +b_20 = a_11-b_20 = 0,


 
(v)  
{[             δ = λ-1 
= 2a_11b_20+3a_20^2-2b_20^2=2b_11+5a_20=0,; 8a_02b_20^2+a_20(8b_20^2-3a_20^2)=4b_02b_20
-3a_20^2+4b_20^2=0, ].


 
(vi)  
{[            δ = λ-9/2 = 12a_11b_20+27a_20^2 
-4b_20^2=4b_11+11a_20=0,; 32a_02b_20^2+3a_20(8b_20^2-9a_20^2)=24b_02b_20
-27a_20^2+16b_20^2=0. ]. 


 
First of all, note that δ = 0 is a necessary condition for all cases
in order to get limit cycles bifurcating from the origin (or infinity), 
under which L_0 = 0. 

We start from Case (B) in which b_20=0. 
It is easy to see that for this 
case all L_i,   i=2,3, …,6 contains a same factor 
b_02 ( 2 a_20 + b_11). Thus,  
L_i = 0,   i=2,3, … 6, if b_02 ( 2 a_20 + b_11) = 0, 
implying that the maximal limit cycles can be obtained is two.
When b_20 = b_02 =0, L_2 =0, and L_1 = 0 yields one solution: 
a_11 = 0, which gives the condition (i). 
If b_20 = 2 a_20 + b_11 =0, then L_2 = 0, and 
L_1 = 0 requires a_11 + 2 b_02 = 0, which 
yields the condition (ii). 

The Case (A3) in which b_20 = -1/7 b_02 (6 + λ) is 
simple since it is assumed b_20 0 and so b_02. 
It is also assumed a_20 0 for this case, yielding 
L_3  0, implying that three limit cycles can be obtained 
since one can choose appropriate values of a_11 and b_11 
to set L_1 = L_2 = 0. 

Next, consider Case (A2) in which a_20=0. 
Maximal 5 limit cycles may be obtained by choosing b_02 0, and 
setting  b_11 = 0 (so L_4 = 0), 
b_20 = 0 (so L_2 = L_3 =0), and 
a_11 + 2 b_02 = 0 (so L_1 = 0).  
If b_02 =0, then L_3 = L_4 = ⋯ = L_7 =0. 
Further, setting L_1=L_2 =0 we obtain a_11+ b_20 
= b_11 + a_02 =0, which is the condition (iii). 
Another possibility for this case is to set λ=1, giving 
L_4 = L_5 = L_6 = L_7 = 0. Further, letting 
L_1 = L_2 = L_3 = 0 yields 
a_02 = b_02 + b_20 = a_11 - b_20 = 0, leading to 
the condition (iv). 
 
For Case (A1)(c), it is noted that 
λ = 1 or λ = 9/2 yields 
L_5 = L_7 = 0 (L_6 is already equal to zero). 
L_1 = 0 gives b_02 = -1/2(a_11 + b_20) which
is substituted into a_20^2 = 28 b_20 [7 
b_20 + b_02 ( λ + 6) ] / 3 (λ+6)^2 to obtain 
2 a_11 b_20 + 3 a_20^2 - 2 b_20^2 = 0 for λ=1, 
and 12 a_11 b_20 + 27 a_20^2 - 4 b_20^2 = 0  
for λ = 9/2.  
Then, for case (v) for which λ = 1, 
using a_20^2 = 4/3 b_20(b_20+b_02) 
to simplify L_3 = 0 yields 2 b_11 + 5 a_20 = 0, and 
again using the expression of a_20^2 as well as 
b_11 = -5/2 a_20 to simplify L_2 = 0 we obtain 
8 a_02 + a_20( 8 b_20^2 - 3 a_20^2)=0. Finally, using 
a_11 = 2 b_20^2 - 3 a_20^2/ 2 b_20 to simplify 
L_1 = 0 results in 4 b_02 b_20 - 3 a_20^2 + 4 b_20^2=0. 
Summarizing the above results leads to the condition (v). 
Following the same procedure, we can obtain the condition (vi). 

Now we come to Cases (A1)(a) and (A1)(b). 
First note that for these two cases, b_02 0 since it is assumed
b_20 0. Therefore, instead of considering the equations 
L_5 = L_6 = L_7 = 0, 
we consider the polynomial equations f_1 = f_2 = f_3 = 0 
for Case (A1)(a) and f̃_1 = f̃_2 = f̃_3 = 0 
for Case (A1)(b). 
In order to obtain the maximal number of limit cycles, we need to find 
the conditions such that L_i=0, i=0,1, …, k-1, but L_k  0. 
Actually,  setting L_i=0,   i=0,1,2,3,4, we obtain 
two sets of solutions: 

    δ    = 0, 
     
    b_20   = -  b_02/4  (λ+3), 
     
    a_11   = b_02/4  (λ-5),
    
    a_02    =  a_20^2 λ - b_02^2 (λ-1) /6a_20,
    
    b_11    = -  a_20^2 (λ^2 + 21 λ + 6) 
    - b_02^2 (λ-1)(λ+3) /6 a_20 (λ-1),
 
for Case (A1)(a), and 

    δ    = 0, 
     
    b_20   =- b_02/12  (λ+9), 
     
    a_11    = b_02/12  (λ-15), 
    
    a_02   = -  9 a_20^2  λ  ( λ+13) 
                       + b_02^2 (λ-3) (5 λ+9) /18 a_20 (5 λ+9),
    
    b_11   =  a_20^2 (9 λ^2 - 45 λ + 162) 
    + b_02^2 (λ+9)(5 λ+9) /18 a_20 (5 λ+9),

for Case (A1)(b). 

Note that the functions f_1, f_2, f_3 and f̃_1, f̃_2, f̃_3 are homogeneous polynomials in a_20^2 and b_02^2. 
So we may introduce b_02^2 = k  a_20^2 (k>0) 
into these polynomials and note that f_2 and f̃_2 
have two factors, one is linear in k and one is quadratic in k, 
given as follows: 

    [                        f_2 = - a_20^6 f_2a f_2b =      - a_20^6 
    [ 7 (λ-1)(λ +3) k - 4 (λ + 6)^2  ];                                           [0.5ex] ×[ 
    4 (λ-1)^2 k^2 + (λ-1)(25 λ -81) k +196 λ^2 ]. ]
 
We first solve f_2a=0 to obtain 
k= 4 (λ + 6)^2/7 (λ-1)(λ +3) 
which is then substituted into f_1 and f_3 to yield 

    [                                                 f_1 = - 1296 a_20^4/7 (λ -1)(2 λ - 9) 
    (5 λ + 16) ( λ+3)^2,;                                                 f_3 =  -699840 a_20^8/49(λ -1)^2(2 λ - 9) 
    ( λ + 3)^3   C_3 ]
 
where 

    C_3 = (1587 λ^3 - 82636 λ^2 -  527988 λ - 838080).
 

Since (λ-1)(λ+3)  0, the only solution satisfies 
f_1 = 0 and f_3  0 is λ = -16/5. When 
λ = 9/2, L_i=0,  i=0,1, …, 7, but it is easy 
to verify that this is a special case of (vi). 

Moreover, for the solution λ = -16/5, 

    (J)= [ ∂ (L_5, L_6)/∂(k, λ)]  
    = 8384103915264/78125  a_20^10 0    for   a_20 0,
 
implying that seven limit cycles can bifurcate in the small neighborhood 
of infinity. 

For the second factor f_2b, we eliminate k from the two equations 
f_1 = f_2b = 0 to obtain the solution for k: 

    k = - 4 λ^2 (6337 λ^2+27872 λ+27783)/(λ-1)( 1721 λ^3 - 2945 λ^2 - 38097 λ-45927,
 
and a resultant equation,  

    [                                                     R_12 = λ (λ +1) (λ -1)
          (397378 λ^4+3696797 λ^3+12760835 λ^2;                                                                                           +19311435 λ+10762227) = 0 , ]
 
which has three real solutions: λ = -2.59473685 ⋯,  
-1.58363608 ⋯,  -1,  but none of them yields k>0. 
Hence, there are no solutions from f_2b=0 to generate seven limit cycles. 

 
When the conditions in (<ref>) are satisfied, similarly, we 
can use b_02^2 = k̃  a_20^2 to find that 

    [                                                                                  f̃_2=                                                                      a_20^6 f̃_2af̃_2b;                                                                                      =      a_20^6 [ 7 (λ+9) (5 λ+9) k̃ + 108 (λ+6)^2 ] [ 4 (5 λ+9)^2 (23 λ^2-30 λ-153) 
    k̃^2;                                                                                        +9 (5 λ L + 9) (339 λ^3 + 593 λ^2 
    - 6459 λ- 12393) k̃+ 108 λ^2 (λ- 1) (109 λ+ 339) ]. ]
 
Solving f̃_2a=0 to obtain 
k̃ = -108 (λ+6)^2/7 (λ+9) (5 λ+9), 
and then substituting it into f̃_1 and f̃_3 yields 

    [                                                              f̃_1 =              -34992 a_20^4/49 (λ -1)(2 λ - 9) 
    (5 λ + 16) ( λ+9)^2,;                                                             f̃_3 
    = -17061120 a_20^8/49(λ -1)^2(2 λ - 9) 
    (5 λ + 9)^2 ( λ + 9)^3   C_3, ]
  
where C_3 is given in (<ref>). 
Since (5 λ+9)(λ+9)  0, the only solution satisfies 
f_1 = 0 and f_3  0 is λ = -16/5. 
λ =1 and λ = 9/2 are not solutions 
since they yield k̃<0. 
For the solution λ = -16/5, we have 

    (J)= [ ∂(L_5,L_6)/∂(k̃, λ)] 
    = -  304277047914997479168/78125  a_20^10 0    for   a_20 0,
 
implying that seven limit cycles can bifurcate in the small neighborhood 
of infinity. 

For the second factor f̃_2b, similarly we eliminate k̃ 
from the two equations f̃_1=0 and f̃_2b=0 to
obtain the solution for k̃: 

    k̃= - 36 λ^2 (λ-1)
              (223 λ^4+8120 λ^3-127950 λ^2-778680 λ
               -741393)/(5 λ+9)
              /(6419 λ^6+31564 λ^5-1363197 λ^4
              -4806072 λ^3+20436705 λ^2 +88888428 λ+81310473),
 
and a resultant equation, 

    [                                                          R̃_12 
    = λ (λ -1) (λ -9) (5 λ +9)
    (8397602 λ^9+84616511 λ^8-1494342124 λ^7;                                                                                  -16706405616 λ^6 +16325397720 λ^5+375950483190 λ^4;                                                                            - 756410507892 λ^3 -10202463072792 λ^2
    -22941003813786 λ;                                                                                                                 -  15945065065773), ]
 
which has four real solutions: λ = -12.38286360 ⋯,  
-8.54108488 ⋯, 1,  9, 13.13951523 ⋯, 
but all of them yield k̃≤ 0.
Thus, there are no solutions from f̃_2b=0 to give 
seven limit cycles. 

Summarizing the above results obtained for Cases (A1)(a) and (A1)(b), 
we conclude that there exist two sets of infinite solutions 
such that system (<ref>) 
can have seven limit cycles bifurcating in the small neighborhood of 
infinity. 

Although we can not obtain seven limit cycles around the origin of 
system (<ref>), we may find an infinite number of 
solutions for six limit cycles which 
bifurcate in the small neighborhood of the origin, which is still better 
than the five limit cycles obtained in <cit.>. 
To find the solutions, it needs L_i,  i=0,1, ⋯, 5, but 
L_6  0. Thus, we only need to solve f_1 = 0 (or f̃_1 =0).  
Letting b_02^2 = k a_20^2, f_1 = 0 becomes 
f_1 = -  a_20^2 [ A_2 k^2 + A_1 k + A_0 ], where 

    [          A_2 = 7 (λ - 1)^2 (λ+3) (155 λ^2-243),; A_1 = λ (λ-1) (2909 λ^3+8871 λ^2+2547 λ+20169),;           A_0 = -4 λ^3 (967 λ^2+8580 λ +20637). ]
  
We want to find solutions satisfying k>0 and λ >0  (λ 1). 
It is easy to show that f_1 = 0 has a unique positive solution 
for k when λ≥ 9√(3/155), and 
does not have solutions when 1 < λ < 9√(3/155). 
When λ∈ (0,1), A_i < 0,   i=0,1,2, and 

    [                               Δ = A_1^2 - 4 A_2 A_0 = -81 λ^2 (1-λ)^3 (311721 λ^5+3409519 λ^4
    +14178654 λ^3;                                                                          +25596434 λ^2+14511609 λ-5022081), ]

which is positive for λ∈ (0, 0.23512585), leading to 
that f_1=0 has two solutions for k for each λ chosen 
from this interval.  
Hence, there exists an infinite number of solutions for 
k satisfying f_1 =0 when λ∈ (0, 0.23512585) 
∪ (9√(3/155), ∞).  
These solutions do not yield f_2 =0, since in the above we have already 
shown that f_1 = f_2 =0 does not have solutions satisfying 
k>0 and λ >0. 
This indicates that for Case (A1)(a) there exists an infinite 
number of solutions for the existence of six limit cycles around the 
origin of system (<ref>). Similarly, for Case (A1)(b),
we can prove that f_1 = 0 has two positive solutions for k when 
λ∈ (0, 1) and one positive solution when 
λ∈ (1,  5.132341426), implying that for Case (A1)(b)
there also exists an infinite number of solutions 
for the existence of six limit cycles around the
origin of system (<ref>). 

The proof is complete.
 


 
Note that the conditions (i)–(vi) given in Theorem <ref> yield 
L_i=0,  i=0,1, …, 7, implying that they are necessary conditions 
for the origin (or infinity) of system (<ref>) to be a center.  
In the following, we will show that these conditions are also sufficient 
for the origin (or infinity) of system (<ref>) to be a center.
We have the following theorem. 

 
The conditions (i)–(vi) given in Theorem <ref> 
are necessary and sufficient for the origin (or infinity) 
of system (<ref>) to be a center.
 

 
The necessity has been shown in the proof of Theorem <ref>. 
Hence, we only need to prove the sufficiency.
First note that for all the six cases, the lower-half plane is same 
(as δ = 0), described by

    [ ẋ=-y,;  ẏ=x, ]     (y<0).

This system has a first integral H_0(x,y) = x^2 + y^2, 
which is an even function of x (i.e., symmetric with the y-axis). 
Thus, in the following, for each case we only list the equations for
the upper-half plane.

When the condition (i) holds,
the equations for the upper-half plane of
system (<ref>) become

    [ ẋ=-y+ (a_20x^2+a_02y^2) (x^2+y^2)^(λ-1)/2 ,;            ẏ=x+ b_11xy  (x^2+y^2)^(λ-1)/2 , ]     (y>0),

which is symmetric with the y-axis, and so 
the origin (or infinity) is a center.

When the condition (ii) is satisfied, the equations for the
upper-half plane of system (<ref>) can be rewritten as

    [ ẋ=-y+ (a_20x^2-2b_02xy+a_02y^2) (x^2+y^2)^(λ-1)/2,;          ẏ=x- (2a_20x-b_02y)y  (x^2+y^2)^(λ-1)/2 , ]     (y>0),

which has an integrating factor
2/3λ(x^2+y^2)^(1-λ)/2.
Then, system (<ref>) becomes 

    [ ẋ=-2/3λ y(x^2+y^2)^(1-λ)/2
    +2/3λ(a_20x^2-2b_02xy+a_02y^2),;          ẏ=2/3λ x(x^2+y^2)^(1-λ)/2
    -2/3λ(2a_20x-b_02y)y ), ]

which has a first integral,

    [                                                                      H_1(x,y)= 6 λ (x^2 + y^2)^(3-λ)/2  
    + 2 λ (λ-3)  y 
    [ 3 (a_20x - b_02y)x  + a_02y^2 ) ]. ]
 
It is seen that H_1 is an even function of x when y=0, and so
the origin (or infinity) of system (<ref>) is 
a center (e.g., see Theorem 2.2 in <cit.>). 

When the condition (iii) holds, the equations for the upper-half
plane of system (<ref>) can be rewritten as

    [ ẋ=-y - (b_20x-a_02y) y  (x^2+y^2)^(λ-1)/2 ,; ẏ=x + (b_20x-a_02y) x   (x^2+y^2)^(λ-1)/2 , ]     (y>0).

It is easy to see that system (<ref>) has a first integral,

    H_2(x,y)=(x^2 + y^2)^λ/3,

which is an even function of x and so
the origin (or infinity) is a center.

When the condition (iv) is satisfied, the equations for the
upper-half plane of system (<ref>) become

    [                  ẋ=-y+b_20xy,; ẏ=x+(b_20x^2+b_11xy-b_20y^2), ]     (y>0).

It can be shown that system (<ref>) has a first integral,

    H_3(x,y)= (b_20x-1)[ b_20x+1+12(b_11-γ)y ]^α[ b_20x+1+12(b_11+γ)y ]^(1-α),

where

    α= 4b_20^2/γ(γ+b_11), 
      γ=√(b_11^2+8b_20^2).

H_3(x,y) is an even function of x when y=0 because
H_3(x,0) =b_20^2x^2-1,
and so the origin is a center <cit.>. 

When the condition (v) holds,
the upper-half plane has a first integral,

    H_4(x,y) = (2b_20x-a_20y-2)^2 [ 4(b_20x + 1)^2
    -(4a_20+12a_20 b_20x)y+(3a_20^2-8b_20^2)y^2  ],

which is an even function of x when y=0 since 
H_4(x,0)=16(1-b_20^2x^2)^2.
So the origin is a center <cit.>. 

Finally, when the condition (vi) is satisfied, the equations
for the upper-half plane of system (<ref>) become

    ẋ=    -y+ 1/96 b_20^2(x^2 + y^2)^7/4[ 96 a_20 b_20^2 x^2 - 8 b_20 ( 27 a_20^2 - 4 b_20^2) x y
       + 9 a_20 (9 a_20^2 - 8  b_20^2) y^2 ],
    ẏ=   x+1/24  b_20 (x^2+y^2)^7/4[ 24 b_20^2 x^2 - 66 a_20 b_20xy + (27 a_20^2 - 16 b_20^2 )y^2  
    ],

which has a first integral,

    H_5(x,y)=    9 a_20^2/64b_20^2 
    (9a_20^2+ 16b_20^2)(x^2+y^2)^3/2{ 4096 b_20^4 
        + (4 b_20 x -3 a_20 y)^4 
    [ 8 b_20^2 (x^2 - y^2) 
    - 3 a_20 (8 b_20 x - 3 a_20 y) y ] 
    (x^2 + y^2)^3/2
        + 128 b_20^2 y (4 b_20 x - 3 a_20 y)
    (24 a_20 b_20 x - 9 a_20^2 y + 16 b_20^2 y) 
    (x^2 + y^2)^3/4}.

Since 
H_5(x,0)=-288 a_20^2 b_20^2 (2 +b_20^2 x^6 
|x|^3)/(9 a_20^2 + 16b_20^2) |x|^3,
H_5(x,y) is an even function of x when y=0, 
and hence, the origin is a center <cit.>. 


Combining the results in Theorems <ref> and <ref>, we have 


For system (<ref>), the highest order of focus value is 7.





 §.§ Isochronous centers of system (<ref>)


Having established the center conditions in the previous section for
system (<ref>), we now discuss the isochronous center conditions for
this system.

First, by a direct computation, we can show that
under any of the conditions (ii)–(vi), 
no isochronous center can exist because these conditions cannot lead to
all periodic constants vanishing. When the condition (i) holds,
the periodic constants are obtained as

    [                                  τ_1= 2/3 (a_20 + 2 a_02-b_11),;                    τ_2= π/16a_02^2 [ 2a_20λ 
    + a_02 (3 λ+5) ] ,;                            [1.8ex]
    τ_3= - 2a_02^3/105  (λ^2-1),;                   [1.8ex]
    τ_4= 3 π a_02^4 /1024  (λ^2-1) (λ+3),;          τ_5= - 2 a_02^5/1001  (λ^2-1)  
    (14 λ^2 + 91λ + 139 ),; τ_6= 3 π a_02^6/4587521024  (λ^2-1) 
    (1190λ^2 + 8817λ +13898 ),;                                                               ⋮ ]

Therefore, we have the following theorem.


The origin (or infinity) of system (<ref>) is an isochronous center 
if and only if one of the following conditions holds:

    [                                                  (I)                    a_11=b_02=b_20=a_02= b_11-a_20=0;;                                                 (II) λ-1= a_11=b_02=b_20= b_11 +2a_02 
    = a_20 +4a_02 = 0;;                                                (III)   λ +1= a_11=b_02=b_20= b_11-3a_02 
    = a_20-a_02 = 0. ]





The necessity can be easily proved by setting 
the periodic constants τ_i = 0,   i=1,2, ⋯ 6. 
To prove the sufficiency, we consider the three 
systems under the three conditions (I), (II) and (III). 
First, consider condition (I). When this condition holds,  
the equations for the upper-half plane
of system (<ref>) can be written as

    [ ẋ=-y+a_20x^2(x^2+y^2)^(λ-1)/2,;   ẏ=x+a_20xy(x^2+y^2)^(λ-1)/2. ]

A simple calculation gives dθ/dt=1, and so the origin 
(or infinity) of system (<ref>) is an isochronous center. 

When the condition (II) is satisfied, system (<ref>) becomes

    [ [ ẋ=-y - a_02 (4 x^2-y^2),;             ẏ=x-2a_02xy, ]                                                 (y>0),;                     [        ẋ = -y,; [0.0ex]
    ẏ = x, ]                                                 (y<0). ]

This system has a transversal system

    [ [ ẋ =x- 4 a_02  xy   (1-a_02y),;    ẏ =y- a_02 y^2 (3-2a_02y), ]                                                           (y>0),;                                                   [ ẋ=x,; ẏ=y, ]                                                           (y<0). ]

So by Theorem 2.1 in <cit.>, the origin of system (<ref>)
is an isochronous center.

When the condition (III) is satisfied, introducing the transformation,

    u=x (x^2 + y^2)^-2/3,    v=y (x^2 + y^2)^-2/3,

into system (<ref>) yields

    [ [ u̇= - v - 1/3  
     a_02 (u^4 + 10 u^2 v^2 - 3v^4) ,;        v̇= u + 1/3  
    a_02 (5 u^2 - 7 v^2 )  u v , ]                                                                                                   (v>0),;                                                                                       [ u̇=-v,;  v̇=u, ]                                                                                                   (v<0). ]

The upper-half plane has an inverse integrating factor,

    V(u,v)=(u^2 + v^2)^2 f_7^5/6,

from which a first integral can be obtained as

    H(u,v)=u^2 + v^2/f_7^1/6 +16 a_02 (u^2 + v^2)
    ∫ u f_7^-5/6du,

where

    f_7=9 [ 1 + 2 a_02  v (u^2 - v^2) + a_02^2 v^2 (u^2 + v^2 )^2 ] .

Thus, this system has a transversal system,

    [ [        u̇=u (3-9a_02^2v^3 )  
    u^2+v^2/f_7^1/6H(u,v),; v̇=v (3+9 a_02u^2v-3a_02v^3 )
    u^2+v^2/f_7^1/6H(u,v), ]                                                                                                         (v>0),;                                                                                               [ u̇=u,; v̇=v, ]                                                                                                         (v<0), ]

implying that infinity of system (<ref>) is an isochronous center.





§ CONCLUSION


In this paper, quasi-analytic switching systems have been considered.
A modified and improved method for computing the return maps of
quasi-analytic switching systems is presented.
In particular, a quadratic quasi-analytic switching system is investigated
using this method. Center and isochronous center conditions are
explicitly derived. Compared to the special case λ=1 for which
five limit cycles are obtained around the origin <cit.>,
we have shown that there exists an infinite number of solutions of 
λ>0, λ 1
for the existence of six limit cycles around the 
origin, and two solutions for the existence of 
seven limit cycles for λ = -16/5,
which bifurcate in the small neighborhood of infinity of the system.
This shows that the dynamics of quasi-analytic switching systems
are more complex.





 
Acknowledgements  This research is partially supported by 
the National Nature Science
Foundation of China (NSFC No. 11371373 and No. 11601212)
and Applied Mathematics Enhancement Program of Linyi University,
and the Natural Science and Engineering Research Council of Canada
(NSERC No. R2686A02). 





















99

Amelikin1982
Amelikin B, Lukashivich H, Sadovski A.
Nonlinear Oscillations in Second Order Systems.
Minsk, BGY Lenin, B. I. Press, 1982 (in Russian).


 
Bautin1952
Bautin N N. 
On the number of limit cycles which appear with the 
variation of coefficients from an equilibrium position of focus or 
center type. 
Mat. Sbornik, 1952, 30: 181–196; 
Amer. Math. Soc. Transl., 1954, 100: 1–19


 
Chen2010
Chen X, Du Z.
Limit cycles bifurcate from centers of discontinuous quadratic systems.
Comp. Math. Appl., 2010, 59: 3836–3848


 
Du2005
Du Z, Zhang W.
Melnikov method for homoclinic bifurcation in nonlinear impact oscillators,
Comp. Math. Appl., 2005, 50: 445–458


 
Dulac1908
Dulac H. 
Dtermination et integration d'une certaine classe 
d'quations diffrentielle ayant par point singulier un centre. 
Bull. Sci. Math. Ser. (2), 1908, 32: 230–252


 
Freire2005
Freire E, Ponce E, Ros J.
The focus-center-limit cycle bifurcation in symmetric 3D piecewise 
linear systems.
SIAM. J. Appl. Math., 2005, 65: 1933–1951 


 
Friederiksson2000
Friederiksson M. H., Nordmark A. B.,
On normal form calculations in impact oscillators.
Proc. R. Soc. Lond. Ser. A, 2000, 456: 315–329 


 
GT2003
Gasull A, Torregrosa J.
Center-focus problem for discontinuous planar differential equations.
Int. J. Bifurcation and Chaos, 2003, 13: 1755–1765 


 
Hauer2012
Hauer B, Engelhardt A P, Taubner T.
Quasi-analytical model for scattering infrared near-field microscopy 
on layered systems.
Optics Expr., 2012, 12: 13173–13188 


 
Kapteyn1911
Kapteyn W. 
On the midpoints of integral curves of differential 
equations of the first degree.
Nederl. Akad. Wetensch. Verslag. Afd. Natuurk. Konikl.
Nederland, 1911, 2: 1446–1457 (in Dutch)


 
Kapteyn1912
Kapteyn W. 
Newinvestigations on the midpoints of integrals of 
differential equations of the first degree. 
Nederl. Akad. Wetensch. 
Verslag Afd. Natuurk., 1912, 20: 1354–1365 (in Dutch)


 
Kukuka2007
Kukuka P.
Melnikov method for discontinuous planar systems.
Nonlinear Anal. Ser. A, 2007, 66: 2698–2719


 
Leine2004
Leine R I, Nijmeijer H.
Dynamics and Bifurcations of Nonsmooth Mechanical Systems.
in: Lecture Notes in Applied and Computational Mechanics, vol. 18, 
Berlin: Springer-Verlag, 2004


 
Li2015
Li F, Yu P, Tian Y, Liu Y.
Center and isochronous center conditions for switching systems associated 
with elementary singular points.
Commun. Nonl. Sci. Numer. Simulat., 2015, 28: 81–97 


 
LL2012
Liu Y, Li J.
Center and isochronous center problems for quasi analytic systems.
Acta Math. Sinica, 2012, 24: 1569–1582 


 
Liu2002
Liu Y.
The generalized focal values and bifurcations of limit circles for 
quasi-quadratic system.
Acta. Math. Sinica, 2002, 45: 671–682 


 
LLH2009
Liu Y, Li J, Huang W.
Singular point vaules,center problem and bifurcations os limit circles 
of two dimensional differential autonomous systems.
Beijing, Science Press, 2009, 162–190 


 
Llibre2014
Llibre J, Lopes D, Moraes R.
Limit cycles for a class of continuous and discontinuous cubic 
polynomial differential systems.
Qual. Theory Dyn. Syst., 2014, 13: 129–148 


 
LM2014
Llibre J, Mereu C.
Limit cycles for discontinuous quadratic differential systems with two zones.
J. Math. Anal. Appl., 2014, 413: 763–775 


 
LV2009a
Llibre J, Valls C.
Classification of the centers, their cyclicity and isochronicity for 
a class of polynomial differential systems generalizing the linear 
systems with cubic homogeneous nonlinearities.
J. Differential Equations, 2009, 246: 2192–2204


 
LV2009b
Llibre J, Valls C.
Classification of the centers and isochronous centers for a class 
of quartic-like systems.
Nonlinear Anal., 2009, 71: 3119–3128


 
LV2009c
Llibre J, Valls C.
Classification of the centers, their cyclicity and isochronicity 
for the generalized quadratic polynomial differential systems.
J. Math. Anal. Appl., 2009, 357: 427–437


 
LV2009d
Llibre J, Valls C.
Classification of the centers, of their cyclicity and isochronicity 
for two classes of generalized quintic polynomial differential systems.
Nonlinear Differ. Eqn. Appl., 2009, 16: 657–679


 
Loud1999
Lloyd N G, Pearson J M. 
Bifurcation of limit cycles and integrability 
of planar dynamical systems in complex form. 
J. Phys. A: Math. Gen., 1999, 32: 1973–1984


 
Loud1964
Loud W S. 
Behavior of the period of solutions of certain plane 
autonomous systems near centres.
Contrib. Differential Equations, 1964, 3: 323–336


 
Pavlovskaia2007
Pavlovskaia E, Wiercigroch M.
Low-dimensional maps for piecewise smooth oscillators.
J. Sound Vib., 2007, 305: 750–771 


 
Pleshkan1969
Pleshkan I. 
A new method for investigating the isochronicity of a 
system of two differential equations. 
Differential Equations, 1969, 5: 796–802


 
Musolino2012
Musolino A, Rizzo R, Tripodi E.
A Quasi-Analytical Model for Remote Field Eddy Current Inspection.
Pro. Elec. Rese., 2012, 26: 237–249 


 
Oerlemans1999
Oerlemans J.
A quasi-analytical ice-sheet model for climate studies.
Nonlinear Pro. Geop., 1999, 10: 441–452 


 
Sabatini1992
Sabatini M.
Characterizing isochronous centers by Lie brackets.
Diff. Equa. Dyna. Syst., 1997, 5: 91–99 


 
Yu2015
Tian Y, Yu P.
Center conditions in a switching Bautin system.
J. Differential Equations, 2015, 259: 1203–1226 


 
Vulpe1988
Vulpe N I, Sibirskii K S. 
Centro-affine invariant conditions for 
the existence of a center of a differential system with cubic nonlinearities.
Dokl. Akad. Nauk SSSR, 1988, 301: 1297-1301. 
Translation in Soviet Math. Dokl., 1989, 38: 198–201 (in Russian)


 
P2005
Xiao P.
Critical point quantities and integrability conditions for
complex planar resonant polynomial differential systems.
Central South University, Ph.D. Thesis, 2005  


 
YuHan2012
Yu P, Han M. 
Four limit cycles from perturbing quadratic integrable systems by 
quadratic polynomials. 
Int. J. Bifurcation and Chaos, 2015, 22: 1250254 (28 pages)


 
Zoladek1994
Żoła̧dek H. 
Quadratic systems with center and their perturbations.
J. Differential Equations, 1994, 109: 223–273


 
Zou2006
Zou Y, Kpper T, Beyn W-J.
Generalized Hopf bifurcation for planar Filippov systems continuous 
at the origin.
J. Nonlinear Sci., 2006, 16: 159–177




