# On the Shape of Limit Cycles That Bifurcate from Isochronous Center

Authors: Guang Chen; Yuhai Wu
Journal: The Scientific World Journal
 (2014)

---

## Abstract

New idea and algorithm are proposed to compute asymptotic expression of limit cycles bifurcated from the isochronous center. Compared with known inverse integrating factor method, new algorithm to analytically computing shape of limit cycle proposed in this paper is simple and easy to apply. The applications of new algorithm to some examples are also given.


---

## Body

## 1. Introduction

Many physical, chemical, and biological systems show periodic activity. Mathematically, they can be modeled by limit cycles of vector field. For example, in [1], Van der Pol proved that a closed trajectory of a self-sustained oscillation occurring in a vacuum tube circuit was a limit cycle as defined by Poincaré. The study of limit cycles of real general planar vector field is closely related to Hilbert’s 16th Problem. As to the strongly nonlinear oscillation equation d
x
/
d
t
=
y, d
y
/
d
t
=
g
(
x
)
+
ɛ
f
(
x
,
y
), in [2], the first two order approximate expressions of limit cycles for small positive parameter ɛ were studied by the generalized KBM method, and, in [3], the shape of the limit cycles for moderately large positive parameter ɛ was plotted by using the perturbation-incremental method.In 1881–1886, Poincaré defined a center of planar vector field as an isolated singular point surrounded by a family of periodic orbits. Then one interesting problem is to ask whether limit cycles appear near the periodic orbits in the vicinity of the center as the planar vector field having a center is perturbed, and what are the shapes of these limit cycles if they exist? Literatures [4, 5] have applied the method of inverse integrating factor to analytically compute global shape of the limit cycles bifurcated from analytic isochronous center. The main idea of determining the shape of limit cycles of planar vector field (
P
(
x
,
y
)
,
Q
(
x
,
y
)
) in [4, 5] is to determine function V
(
x
,
y
)
=
∑
k
=
0
∞
‍
ε
k
V
k
(
x
,
y
) which satisfies the partial differential equation
(1)
P
∂
V
∂
x
+
Q
∂
V
∂
y
=
(
∂
P
∂
x
+
∂
Q
∂
y
)
V
,
and the limit cycles of planar vector field (
P
(
x
,
y
)
,
Q
(
x
,
y
)
) are implicitly determined by V
(
x
,
y
)
=
0. In other words, if one tries to find analytic expression of limit cycle, one should solve linear partial differential equations recursively. In this paper, a new idea and algorithm are developed to analytically compute the shape of the limit cycles bifurcated from the isochronous center. From Theorem  3.2 in [6], we know that any planar analytic system having isochronous center can be locally transformed into the above linear system x
˙
=
y, y
˙
=
-
x by analytic variable transformation and time scale. So without losing generality, we consider analytic expression of limit cycle of perturbed planar vector field x
˙
=
y
+
ɛ
U
(
x
,
y
,
ɛ
), y
˙
=
-
x
+
ɛ
V
(
x
,
y
,
ɛ
). The new algorithm proposed in the paper is based on the following lemmas.Lemma 1.
If planar analytic vector field(
P
(
x
,
y
)
,
Q
(
x
,
y
)
) has a limit cycle r
=
r
(
θ
) surrounding the origin O
(
0,0
), then r
(
θ
) is a periodic function with period 2
π, where (
x
,
y
)
=
(
r
cos
⁡
θ
,
r
sin
θ
).Proof.
From the periodicity of limit cycle and the property of polar coordinate system, we know that the conclusion of the lemma is true.Lemma 2 (see [7]).
Iff
(
θ
) is a C
2 periodic function with period 2
π, then
(2)
F
(
θ
)
=
∫
0
θ
f
(
ξ
)
d
ξ
=
g
θ
+
φ
(
θ
)
,
where g
=
(
1
/
2
π
)
∫
0
2
π
‍
f
(
ξ
)
d
ξ and φ
(
θ
) is a periodic function with period 2
π.
Further, ifF
(
θ
) is periodic function, then g
=
0.Proof.
Forf
(
ξ
) is a C
2 periodic function with period 2
π, so Fourier coefficients of functions f
(
ξ
) and f
′
(
ξ
) have the following relations:
(3)
a
0
′
=
1
π
∫
0
2
π
f
′
(
ξ
)
d
ξ
=
1
π
[
f
(
2
π
)
-
f
(
0
)
]
=
0
,
a
n
′
=
1
π
∫
0
2
π
f
′
(
ξ
)
cos
⁡
(
n
ξ
)
d
ξ
=
1
π
f
(
ξ
)
cos
⁡
(
n
ξ
)
|
0
2
π
+
n
π
∫
0
2
π
f
(
ξ
)
sin
(
n
ξ
)
d
ξ
=
n
b
n
,
n
=
1,2
,
…
,
b
n
′
=
1
π
∫
0
2
π
‍
f
′
(
ξ
)
sin
(
n
ξ
)
d
ξ
=
1
π
f
(
ξ
)
sin
(
n
ξ
)
|
0
2
π
-
n
π
∫
0
2
π
f
(
ξ
)
cos
⁡
(
n
ξ
)
d
ξ
=
-
n
a
n
,
n
=
1,2
,
…
.
So
(4)
|
a
n
cos
⁡
⁡
(
n
θ
)
+
b
n
sin
⁡
(
n
θ
)
|
≤
|
a
n
|
+
|
b
n
|
=
|
a
n
′
|
n
+
|
b
n
′
|
n
≤
1
2
(
a
n
′
2
+
b
n
′
2
)
+
1
n
2
.
By applying Bessel inequality, we get
(5)
∑
n
=
1
∞
‍
(
a
n
′
2
+
b
n
′
2
)
≤
1
π
∫
0
2
π
‍
f
′
2
(
ξ
)
d
ξ
.
By applying comparison test for convergence of series of functions, we get that Fourier series of f
(
ξ
) is uniformly convergent to f
(
ξ
) on [
0,2
π
]. Rewrite f
(
θ
) into the following Fourier series:
(6)
f
(
θ
)
=
a
0
2
+
∑
n
=
1
∞
‍
(
a
n
cos
⁡
⁡
(
n
θ
)
+
b
n
sin
⁡
(
n
θ
)
)
,
where
(7)
a
0
2
=
1
2
π
∫
0
2
π
‍
f
(
ξ
)
d
ξ
=
g
.
Integrating both sides of (6) with respect to variable θ from 0 to θ, we obtain (2) with
(8)
φ
(
θ
)
=
∑
n
=
1
∞
[
a
n
n
sin
⁡
(
n
θ
)
-
b
n
n
(
cos
⁡
⁡
(
n
θ
)
-
1
)
]
.
From the integration property of uniformly convergent series, we get that φ
(
θ
) is a periodic function with period 2
π.
IfF
(
θ
) is a periodic function with period 2
π, then we get g
·
(
2
π
-
0
)
+
φ
(
2
π
)
-
φ
(
0
)
=
0. From φ
(
2
π
)
=
φ
(
0
), we conclude that g
=
0.
The proof of the lemma is completed.The main goal of this paper is to develop a new approach for computing analytically the global shape of the bifurcated limit cycles from an isochronous center and the paper is organized as follows. In Section2, we develop a new algorithm to compute analytic expansion, up to arbitrary order of the parameter ɛ, of the limit cycles bifurcated from linear isochronous center. As applications, in Section 3, we compute the analytic expression of the unique limit cycle of the Van der Pol system x
′
=
y, y
′
=
-
x
+
ɛ
(
1
-
x
2
)
y up to order o
(
ε
7
). In Section 4, we study the analytic expression of the limit cycle bifurcated from a nonlinear isochronous center.


## 2. Asymptotic Expressions of Limit Cycles Bifurcated from the Center ofx
′
=
y, y
′
=
-
x

Consider the following planar system:(9)
x
′
=
y
+
∑
k
=
1
∞
ε
k
p
k
(
x
,
y
)
≡
y
+
ɛ
U
(
x
,
y
,
ɛ
)
,
y
′
=
-
x
+
∑
k
=
1
∞
ε
k
q
k
(
x
,
y
)
≡
-
x
+
ɛ
V
(
x
,
y
,
ɛ
)
,
where p
k
(
x
,
y
) and q
k
(
x
,
y
) are both analytic functions, p
k
(
0,0
)
=
q
k
(
0,0
)
=
0, k
=
1,2
,
…, and ɛ is a small real parameter. System (9) has an isochronous center at the origin when ɛ
=
0. As usual, the prime denotes derivative with respect to variable t. System (9) for ɛ
=
0 is called the unperturbed system, while system (9) for ɛ
≠
0 is called the perturbed one. Then the problem of studying shape of limit cycles bifurcated from isochronous center is to determine the number and analytic expansions of the families of limit cycles which emerge from the periodic orbits of the unperturbed system as the parameter ɛ is varied.The main idea of computing asymptotic expression of limit cycles of system (9) is the following.Firstly, we make a polar coordinates transformationx
=
r
cos
⁡
θ, y
=
r
sin
θ to system (9). By eliminating the variable t, we obtain
(10)
d
r
d
θ
=
ɛ
(
U
cos
⁡
⁡
θ
+
V
sin
⁡
θ
)
-
1
+
ɛ
(
V
cos
⁡
⁡
θ
-
U
sin
⁡
θ
)
/
r
.
By noticing that p
k
(
x
,
y
), q
k
(
x
,
y
) are both analytic functions, we rewrite system (10) into the following form:
(11)
d
r
d
θ
=
R
1
(
r
,
θ
)
ɛ
+
R
2
(
r
,
θ
)
ε
2
+
R
3
(
r
,
θ
)
ε
3
+
⋯
,
where R
1
(
r
,
θ
)
,
R
2
(
r
,
θ
)
,
… are analytic functions about r, cos
⁡
(
θ
), sin
(
θ
).Secondly, to obtain the polar coordinate form of limit cycles of perturbed system (9), we look for the following analytic solution to (11) as ɛ
≠
0,
(12)
r
(
θ
)
=
∑
k
=
0
∞
ε
k
r
k
(
θ
)
.From Lemma1, we know that if r
(
θ
) is limit cycle of system (10), then from the periodicity of r
(
θ
), we get that r
k
(
θ
) is periodic function with period 2
π, where k
=
0,1
,
2
,
….Thirdly, substituting (12) into (11) and considering the kth order terms of ɛ in the obtained system, then we can obtain a series of equations:
(13)
d
r
k
d
θ
=
f
~
k
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
,
k
=
0,1
,
2
,
…
.As to the formula off
~
k
(
θ
,
r
0
,
r
1
,
r
2
,
…
), we have the following lemma.Lemma 3.
Functionsf
~
k
(
θ
,
r
0
,
r
1
,
r
2
,
…
) obtained in (13) have the following properties:
(14)
f
~
0
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
=
0
,
f
~
1
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
=
f
1
(
r
0
,
θ
)
,
f
~
2
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
=
f
2
(
r
1
,
r
0
,
θ
)
,
⋮
f
~
k
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
=
f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
)
,
⋮
where f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
) is analytic function about r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
cos
⁡
(
θ
)
,
sin
(
θ
), k
=
0,1
,
2
,
….Proof.
According to (11) and (12), we get that
(15)
∑
k
=
0
∞
ε
k
d
r
k
(
θ
)
d
θ
=
R
1
(
r
,
θ
)
ɛ
+
R
2
(
r
,
θ
)
ε
2
+
R
3
(
r
,
θ
)
ε
3
+
⋯
.
Fork
=
0, it is easy to get that f
~
0
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
=
0.
Fork
=
1, the term on the right hand side of (15) contributing to ε
1 is ɛ
R
1
(
r
,
θ
). In detail, the constant in the term R
1
(
r
,
θ
) determines f
~
1
(
θ
,
r
0
,
r
1
,
r
2
,
…
). Noticing that r
(
θ
)
=
r
0
+
ɛ
r
1
(
θ
)
+
ε
2
r
2
(
θ
)
+
⋯, we get that f
~
1
(
θ
,
r
0
,
r
1
,
r
2
,
…
) only contains the term r
0. In other words, if function f
~
1
(
θ
,
r
0
,
r
1
,
r
2
,
…
) contains variable r
n
,
n
≥
1, then r
n appearing in the terms in the right hand side of (15) at least have term ε
n
+
1. This is contradiction, for f
~
1
(
θ
,
r
0
,
r
1
,
r
2
,
…
) corresponds to ε
1 term in the right hand side of (15). By using similar analysis, it can be shown that f
~
n
(
θ
,
r
0
,
r
1
,
r
2
,
…
) cannot contain term r
i
,
i
≥
n. Therefore, f
~
k
(
θ
,
r
0
,
r
1
,
r
2
,
…
)
=
f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
), k
=
0,1
,
2
,
….
The proof of the lemma is completed.To obtain analytic expression of limit cycler
=
r
(
θ
) of system (9), we need to determine r
k
(
θ
) in (12), k
=
0,1
,
2
,
…. From Lemma 3, we know that the determinations of r
k
(
θ
) in (12) are recursive.


### 2.1. Determination ofr
0 and the Poincaré-Melnikov Integral

Fromd
r
0
/
d
θ
=
0, we get r
0
(
θ
)
≡
r
0
(
constant
). To determine the constant r
0 in (12), the new approach we adopted is to utilize the expression of r
1
(
θ
).Fromd
r
1
/
d
θ
=
f
1
(
r
0
,
θ
), we obtain
(16)
r
1
(
θ
)
=
∫
0
θ
f
1
(
r
0
,
θ
)
d
θ
+
c
1
.By noting thatf
1
(
r
0
,
θ
) is a periodic function, according to Lemma 2, we know
(17)
∫
0
θ
f
1
(
r
0
,
θ
)
d
θ
=
g
1
θ
+
φ
1
(
r
0
,
θ
)
,
where g
1
=
(
1
/
2
π
)
∫
0
2
π
‍
f
1
(
r
0
,
θ
)
d
θ
≡
g
¯
1
(
r
0
), and φ
1
(
r
0
,
θ
) is periodic function with period 2
π.From Lemma1, we know that if r
=
r
(
θ
) is a limit cycle of system (9), then r
(
θ
) and r
1
(
θ
) in (16) are periodic functions, too.So from Lemma2, we know g
1
=
g
¯
1
(
r
0
)
=
0. By solving that algebraic equation, we can determine the value of constant r
0.Remark 4.
In fact, the functiong
¯
1
(
r
0
) is closely related to the first order Poincare-Melnikov integral of the perturbed system (9) near close orbit of unperturbed system (9)|
ɛ
=
0.
In detail,g
¯
1
(
r
0
)
=
-
(
1
/
2
π
r
0
)
∮
L
q
1
(
x
,
y
)
d
x
-
p
1
(
x
,
y
)
d
y, where close orbit L
:
x
2
+
y
2
=
r
0
2. From [8, 9], we know -
2
π
r
0
g
¯
1
(
r
0
) is the first order Poincare-Melnikov integral. So the zeros of g
¯
1
(
r
0
) are closely related to the number and position of limit cycles of the perturbed system (9).
It should be pointed out that the functiong
¯
1
(
r
0
) is also closely related to the first order averaging of 1-dimensional 2
π-periodic differential equation. First order (resp., second order) averaging method to study the existence and number of periodic orbits of planar differential equation is proposed in [10, 11]. The approach of high order averaging method is based on Brouwer degree theory (see [11] for more details).


### 2.2. Determination ofr
1
(
θ
)

Substitute the value ofr
0 into (16); we can obtain expression of φ
1
(
r
0
,
θ
). Thus we obtain
(18)
r
1
(
θ
)
=
φ
1
(
r
0
,
θ
)
+
c
1
.To determine the value ofc
1, new algorithm proposed in this paper needs the expression of r
2
(
θ
). From d
r
2
/
d
θ
=
f
2
(
r
1
,
r
0
,
θ
), we get
(19)
r
2
(
θ
)
=
∫
0
θ
f
2
(
r
1
,
r
0
,
θ
)
d
θ
+
c
2
.From Lemma2, we know
(20)
∫
0
θ
f
2
(
r
1
,
r
0
,
θ
)
d
θ
=
g
2
θ
+
φ
2
(
r
1
,
r
0
,
θ
)
,
where g
2
=
(
1
/
2
π
)
∫
0
2
π
‍
f
2
(
r
1
,
r
0
,
θ
)
d
θ
≡
g
¯
2
(
c
1
).From the fact thatr
2
(
θ
) is a periodic function and Lemma 2, we get g
¯
2
(
c
1
)
=
0.By solving the above algebraic equation, we determine the value ofc
1. Thus we have obtained r
1
(
θ
) by (18).


### 2.3.  Determination ofr
k
(
θ
)

Assuming that we have obtained the explicit expressions ofr
0
,
r
1
(
θ
)
,
…
,
r
k
-
1
(
θ
), now we start to determine r
k
(
θ
).Fromd
r
k
/
d
θ
=
f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
), we get
(21)
r
k
(
θ
)
=
∫
0
θ
f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
)
d
θ
+
c
k
.To determine the expression ofr
k
(
θ
) is to determine the value of c
k. According to the algorithm proposed in this paper, we resort to the expression of r
k
+
1
(
θ
). From d
r
k
+
1
/
d
θ
=
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
), we get
(22)
r
k
+
1
(
θ
)
=
∫
0
θ
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
)
d
θ
+
c
k
+
1
.From Lemma2, we know
(23)
∫
0
θ
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
)
d
θ
=
g
k
+
1
θ
+
φ
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
0
,
θ
)
,
where
(24)
g
k
+
1
=
1
2
π
∫
0
2
π
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
0
,
θ
)
d
θ
≡
g
¯
k
+
1
(
c
k
)
.Becauser
k
+
1
(
θ
) is a periodic function, from Lemma 2, we get that g
k
+
1
=
g
¯
k
+
1
(
c
k
)
=
0.By solving the algebraic equation, we obtain the value ofc
k; thus we determine r
k
(
θ
) by (20).Thus we can compute the shape of limit cycles of system (9) to any given order of ɛ explicitly and recursively.


## 2.1. Determination ofr
0 and the Poincaré-Melnikov Integral

Fromd
r
0
/
d
θ
=
0, we get r
0
(
θ
)
≡
r
0
(
constant
). To determine the constant r
0 in (12), the new approach we adopted is to utilize the expression of r
1
(
θ
).Fromd
r
1
/
d
θ
=
f
1
(
r
0
,
θ
), we obtain
(16)
r
1
(
θ
)
=
∫
0
θ
f
1
(
r
0
,
θ
)
d
θ
+
c
1
.By noting thatf
1
(
r
0
,
θ
) is a periodic function, according to Lemma 2, we know
(17)
∫
0
θ
f
1
(
r
0
,
θ
)
d
θ
=
g
1
θ
+
φ
1
(
r
0
,
θ
)
,
where g
1
=
(
1
/
2
π
)
∫
0
2
π
‍
f
1
(
r
0
,
θ
)
d
θ
≡
g
¯
1
(
r
0
), and φ
1
(
r
0
,
θ
) is periodic function with period 2
π.From Lemma1, we know that if r
=
r
(
θ
) is a limit cycle of system (9), then r
(
θ
) and r
1
(
θ
) in (16) are periodic functions, too.So from Lemma2, we know g
1
=
g
¯
1
(
r
0
)
=
0. By solving that algebraic equation, we can determine the value of constant r
0.Remark 4.
In fact, the functiong
¯
1
(
r
0
) is closely related to the first order Poincare-Melnikov integral of the perturbed system (9) near close orbit of unperturbed system (9)|
ɛ
=
0.
In detail,g
¯
1
(
r
0
)
=
-
(
1
/
2
π
r
0
)
∮
L
q
1
(
x
,
y
)
d
x
-
p
1
(
x
,
y
)
d
y, where close orbit L
:
x
2
+
y
2
=
r
0
2. From [8, 9], we know -
2
π
r
0
g
¯
1
(
r
0
) is the first order Poincare-Melnikov integral. So the zeros of g
¯
1
(
r
0
) are closely related to the number and position of limit cycles of the perturbed system (9).
It should be pointed out that the functiong
¯
1
(
r
0
) is also closely related to the first order averaging of 1-dimensional 2
π-periodic differential equation. First order (resp., second order) averaging method to study the existence and number of periodic orbits of planar differential equation is proposed in [10, 11]. The approach of high order averaging method is based on Brouwer degree theory (see [11] for more details).


## 2.2. Determination ofr
1
(
θ
)

Substitute the value ofr
0 into (16); we can obtain expression of φ
1
(
r
0
,
θ
). Thus we obtain
(18)
r
1
(
θ
)
=
φ
1
(
r
0
,
θ
)
+
c
1
.To determine the value ofc
1, new algorithm proposed in this paper needs the expression of r
2
(
θ
). From d
r
2
/
d
θ
=
f
2
(
r
1
,
r
0
,
θ
), we get
(19)
r
2
(
θ
)
=
∫
0
θ
f
2
(
r
1
,
r
0
,
θ
)
d
θ
+
c
2
.From Lemma2, we know
(20)
∫
0
θ
f
2
(
r
1
,
r
0
,
θ
)
d
θ
=
g
2
θ
+
φ
2
(
r
1
,
r
0
,
θ
)
,
where g
2
=
(
1
/
2
π
)
∫
0
2
π
‍
f
2
(
r
1
,
r
0
,
θ
)
d
θ
≡
g
¯
2
(
c
1
).From the fact thatr
2
(
θ
) is a periodic function and Lemma 2, we get g
¯
2
(
c
1
)
=
0.By solving the above algebraic equation, we determine the value ofc
1. Thus we have obtained r
1
(
θ
) by (18).


## 2.3.  Determination ofr
k
(
θ
)

Assuming that we have obtained the explicit expressions ofr
0
,
r
1
(
θ
)
,
…
,
r
k
-
1
(
θ
), now we start to determine r
k
(
θ
).Fromd
r
k
/
d
θ
=
f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
), we get
(21)
r
k
(
θ
)
=
∫
0
θ
f
k
(
r
k
-
1
,
r
k
-
2
,
…
,
r
1
,
r
0
,
θ
)
d
θ
+
c
k
.To determine the expression ofr
k
(
θ
) is to determine the value of c
k. According to the algorithm proposed in this paper, we resort to the expression of r
k
+
1
(
θ
). From d
r
k
+
1
/
d
θ
=
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
), we get
(22)
r
k
+
1
(
θ
)
=
∫
0
θ
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
)
d
θ
+
c
k
+
1
.From Lemma2, we know
(23)
∫
0
θ
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
)
d
θ
=
g
k
+
1
θ
+
φ
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
0
,
θ
)
,
where
(24)
g
k
+
1
=
1
2
π
∫
0
2
π
f
k
+
1
(
r
k
,
r
k
-
1
,
…
,
r
0
,
θ
)
d
θ
≡
g
¯
k
+
1
(
c
k
)
.Becauser
k
+
1
(
θ
) is a periodic function, from Lemma 2, we get that g
k
+
1
=
g
¯
k
+
1
(
c
k
)
=
0.By solving the algebraic equation, we obtain the value ofc
k; thus we determine r
k
(
θ
) by (20).Thus we can compute the shape of limit cycles of system (9) to any given order of ɛ explicitly and recursively.


## 3. The Shape of Limit Cycle of Van der Pol System

In this section we will apply the method just described in the above section to compute the analytic expansion of the unique limit cycle of the Van der Pol system(25)
x
′
=
y
,
y
′
=
-
x
+
ɛ
(
1
-
x
2
)
y
up to o
(
ε
7
).First we make a polar coordinates transformationx
=
r
cos
⁡
(
θ
)
y
=
r
sin
(
θ
) to system (25) and eliminate t; then we can obtain
(26)
d
r
d
θ
=
ɛ
r
(
1
-
r
2
cos
⁡
2
θ
)
sin
2
θ
-
1
+
ɛ
(
1
-
r
2
cos
⁡
2
θ
)
sin
θ
cos
⁡
θ
.Assumer
(
θ
)
=
∑
k
=
0
∞
ε
k
r
k
(
θ
) is the polar coordinates form of the limit cycles of (25) and substitute it into (26). By comparing first eight coefficients of terms ε
k
,
k
=
1,2
,
…
,
8 in both sides of the above equation, we get
(27)
d
r
0
d
θ
=
0
,
d
r
1
d
θ
=
-
r
0
(
1
-
r
0
2
cos
⁡
2
θ
)
sin
2
θ
≡
f
1
(
r
0
,
θ
)
,
d
r
2
d
θ
=
(
(
-
1
+
r
0
2
cos
⁡
2
θ
)
r
1
+
2
r
0
2
r
1
cos
⁡
2
θ
)
si
n
2
θ
+
(
-
1
+
r
0
2
cos
⁡
2
θ
)
r
0
×
si
n
2
θ
(
-
cos
⁡
θ
sin
θ
+
r
0
2
sin
θ
cos
⁡
3
θ
)
≡
f
2
(
r
1
,
r
0
,
θ
)
,
⋮
d
r
8
d
θ
≡
f
8
(
r
7
,
r
6
,
…
,
r
0
,
θ
)
.
Here for long expressions, the formula of f
k
(
r
k
-
1
,
…
,
r
1
,
r
0
,
θ
),  k
=
3,4
,
…
,
8 is omitted.Fromd
r
0
/
d
t
=
0, we get that r
0 is arbitrary constant.To determine the value ofr
0, we compute the following expression of r
1
(
θ
)
:
(28)
r
1
(
θ
)
=
∫
0
θ
f
1
(
r
0
,
θ
)
d
θ
+
c
1
=
g
1
θ
+
φ
1
(
r
0
,
θ
)
+
c
1
,
where
(29)
g
1
=
1
2
π
∫
0
2
π
f
1
(
r
0
,
θ
)
d
θ
=
-
r
0
(
1
-
r
0
2
cos
⁡
2
θ
)
sin
2
θ
d
θ
=
1
8
r
0
3
-
1
2
r
0
.
Because g
1
=
0, so we get r
0
=
2.Substituter
0
=
2 into (28); we obtain
(30)
r
1
(
θ
)
=
∫
0
θ
f
1
(
2
,
θ
)
d
θ
+
c
1
=
2
cos
⁡
θ
sin
θ
-
2
cos
⁡
3
θ
sin
θ
+
c
1
.To determine the value ofc
1, we compute the expression of r
2
(
θ
):
(31)
r
2
(
θ
)
=
∫
0
θ
f
2
(
r
1
,
r
0
,
θ
)
d
θ
+
c
2
=
g
2
θ
+
φ
2
(
r
1
,
r
0
)
+
c
2
,
where
(32)
g
2
=
1
2
π
∫
0
2
π
‍
f
2
(
2
cos
⁡
θ
sin
θ
-
2
cos
⁡
3
θ
sin
θ
+
c
1
,
2
,
θ
)
d
θ
=
c
1
.Becauseg
2
=
0, so we get c
1
=
0.Thus explicit expression ofr
1
(
θ
) is given by (30) with c
1
=
0.Substituter
1
(
θ
) and r
0
=
2 into (31); we obtain
(33)
r
2
(
θ
)
=
∫
0
θ
f
2
(
r
1
,
r
0
,
θ
)
d
θ
+
c
2
=
2
cos
⁡
2
θ
-
23
2
cos
⁡
4
θ
+
49
3
cos
⁡
6
θ
-
7
cos
⁡
8
θ
+
1
6
+
c
2
.
In a similar way, we determine the value of c
2 and obtain following results:
(34)
r
2
(
θ
)
=
2
cos
⁡
2
θ
-
23
2
cos
⁡
4
θ
+
49
3
cos
⁡
6
θ
-
7
cos
⁡
8
θ
+
17
96
,
r
3
(
θ
)
=
1
2304
(
(
9
θ
)
486
cos
⁡
θ
+
544
cos
⁡
(
3
θ
)
-
24
cos
⁡
(
5
θ
)
r
3
(
θ
)
=
g
h
h
h
h
-
621
cos
⁡
(
7
θ
)
-
297
cos
⁡
(
9
θ
)
)
sin
3
θ
,
r
4
(
θ
)
=
-
4433
6
cos
⁡
14
θ
+
715
4
cos
⁡
16
θ
-
1577
552960
r
4
(
θ
)
=
-
85021
1728
cos
⁡
6
θ
+
21505
18
cos
⁡
12
θ
+
17
96
cos
⁡
2
θ
r
4
(
θ
)
=
+
203747
576
cos
⁡
8
θ
-
28153
30
cos
⁡
10
θ
-
1045
1152
cos
⁡
4
θ
,
r
5
(
θ
)
=
-
4199
4096
cos
⁡
9
(
2
θ
)
sin
(
2
θ
)
+
1105
12288
cos
⁡
8
(
2
θ
)
r
5
(
θ
)
=
×
sin
(
2
θ
)
-
16237
46080
cos
⁡
6
(
2
θ
)
sin
(
2
θ
)
r
5
(
θ
)
=
+
412037
921600
sin
(
2
θ
)
cos
⁡
4
(
2
θ
)
r
5
(
θ
)
=
-
105841
518400
cos
⁡
2
(
2
θ
)
sin
(
2
θ
)
r
5
(
θ
)
=
+
623977
33177600
sin
(
2
θ
)
+
169
48
cos
⁡
7
(
2
θ
)
r
5
(
θ
)
=
×
sin
(
2
θ
)
-
18203
49152
cos
⁡
(
2
θ
)
sin
(
2
θ
)
r
5
(
θ
)
=
+
30479
13824
sin
(
2
θ
)
cos
⁡
3
(
2
θ
)
r
5
(
θ
)
=
-
797999
184320
cos
⁡
5
(
2
θ
)
sin
(
2
θ
)
,
r
6
(
θ
)
=
-
44679949
77414400
cos
⁡
7
(
2
θ
)
-
19738393
398131200
cos
⁡
3
(
2
θ
)
r
6
(
θ
)
=
+
16357283
55296000
cos
⁡
5
(
2
θ
)
-
2261
16384
cos
⁡
11
(
2
θ
)
r
6
(
θ
)
=
-
10462939
176947200
+
207757
442368
cos
⁡
9
(
2
θ
)
+
215852651
16588800
r
6
(
θ
)
=
×
cos
⁡
6
(
2
θ
)
-
52003
32768
cos
⁡
12
(
2
θ
)
+
5393
26542080
r
6
(
θ
)
=
×
cos
⁡
(
2
θ
)
+
166171249
132710400
cos
⁡
2
(
2
θ
)
r
6
(
θ
)
=
-
8240753
589824
cos
⁡
8
(
2
θ
)
+
122417
16384
cos
⁡
10
(
2
θ
)
r
6
(
θ
)
=
-
1624125463
265420800
cos
⁡
4
(
2
θ
)
,
r
7
(
θ
)
=
-
841356266407
2926264320000
cos
⁡
2
(
2
θ
)
sin
(
2
θ
)
r
7
(
θ
)
=
+
298012397
11796480
cos
⁡
9
(
2
θ
)
sin
(
2
θ
)
r
7
(
θ
)
=
+
3520666153
743178240
cos
⁡
8
(
2
θ
)
sin
(
2
θ
)
r
7
(
θ
)
=
-
14166301337
3612672000
cos
⁡
6
(
2
θ
)
sin
(
2
θ
)
r
7
(
θ
)
=
+
2104280130881
1300561920000
sin
(
2
θ
)
cos
⁡
4
(
2
θ
)
r
7
(
θ
)
=
+
334305
131072
cos
⁡
13
(
2
θ
)
sin
(
2
θ
)
+
260015
393216
r
7
(
θ
)
=
×
cos
⁡
12
(
2
θ
)
sin
(
2
θ
)
-
7451287
589824
cos
⁡
11
(
2
θ
)
sin
(
2
θ
)
r
7
(
θ
)
=
-
6238099
2211840
cos
⁡
10
(
2
θ
)
sin
⁡
(
2
θ
)
r
7
(
θ
)
=
+
283692942563
23410114560000
sin
⁡
(
2
θ
)
-
533135791
20643840
r
7
(
θ
)
=
×
cos
⁡
7
(
2
θ
)
sin
(
2
θ
)
+
8763093
26214400
cos
⁡
(
2
θ
)
sin
(
2
θ
)
r
7
(
θ
)
=
-
5835758083
1592524800
sin
(
2
θ
)
cos
⁡
3
(
2
θ
)
+
74189164357
5308416000
r
7
(
θ
)
=
×
cos
⁡
5
(
2
θ
)
sin
(
2
θ
)
.So the asymptotic expansion of limit cycle of system (25) for ɛ
>
0 and small is the following:
(35)
r
(
θ
)
=
∑
k
=
0
7
‍
ε
k
r
k
(
θ
)
+
o
(
ε
7
)
.The first seven terms in the above expansion ofr
(
θ
) are similar to the ones given in Section 3 of [4]. Here we present the expression of r
7
(
θ
) obtained in our method which was omitted in [4] for its long expression. By applying expansion (35), the shapes of limit cycles of Van der Pol system (25) for the values of ɛ
=
0, (
1
/
5
), (
1
/
2
), (
9
/
10
) are plotted in Figure 1. The periodic orbit x
2
+
y
2
=
4 of system (25) for ɛ
=
0 is drawn in solid line, the limit cycle of system (25) for ɛ
=
1
/
5 is drawn in dashed line, the limit cycle of system (25) for ɛ
=
1
/
2 is drawn in solid line, and the limit cycle of system (25) for ɛ
=
9
/
10 is drawn in dotted line.Figure 1
The periodic orbitx
2
+
y
2
=
4 and the analytical approximations of the limit cycle of Van der Pol system (25), for the values ɛ
=
0, (
1
/
5
), (
1
/
2
), (
9
/
10
).


## 4. The Shape of Limit Cycle of Perturbations of a System Having Nonlinear Isochronous Center

Consider the following perturbed system:(36)
x
′
=
-
y
+
2
x
y
-
2
y
3
+
ɛ
(
x
3
-
3
x
)
,
y
′
=
x
-
y
2
.
From [6], we know that as ɛ
=
0, nonlinear system (36) has isochronous center O
(
0,0
). To utilize new algorithm introduced in Section 2 to study the number and shape of limit cycles of perturbed system (36), we first apply the following analytic variable transformation:
(37)
x
=
u
+
v
2
,
y
=
v
,
t
=
-
τ
to system (36) and get
(38)
d
u
d
τ
=
v
+
ɛ
(
3
u
+
3
v
2
-
(
u
+
v
2
)
3
)
,
d
v
d
τ
=
-
u
.


### 4.1. The Shape of Limit Cycle of the Perturbed System (38)

In this subsection we start to compute the analytic expansion of the limit cycle of the perturbed system (38) to the second order of ɛ.First letu
=
r
cos
⁡
θ, v
=
r
sin
θ; then system (38) is transformed into the following polar coordinate form:
(39)
d
r
d
θ
=
ɛ
(
3
r
cos
⁡
θ
+
3
r
2
sin
2
θ
-
(
r
cos
⁡
θ
+
r
2
sin
2
θ
)
3
)
cos
⁡
θ
×
(
-
1
+
ɛ
(
(
r
cos
⁡
θ
+
r
2
sin
2
θ
)
2
)
sin
θ

-
3
cos
⁡
θ
-
3
r
sin
2
θ
+
(
cos
⁡
θ
+
r
sin
2
θ
)
ghhhhhhhhh
×
(
r
cos
⁡
θ
+
r
2
sin
2
θ
)
2
)
sin
θ
)
-
1
.
According to similar process in Section 3, we get
(40)
g
1
=
3
16
r
0
(
r
0
4
+
2
r
0
2
-
8
)
.By solvingg
1
=
0, we get the positive solution r
0
=
2.From(
d
(
-
2
π
r
0
g
1
)
/
d
r
0
)
|
r
0
=
2
=
-
9
2
π
<
0 and [9], we conclude that as ɛ
>
0, system (38) has a stable limit cycle, denoted by Γ
1, near the close curve u
2
+
v
2
=
2 on the phase plane.By applying the algorithm described in Section2 to system (39), we get the asymptotic expansion of stable limit cycle Γ
1 of system (38),
(41)
r
(
θ
)
=
r
0
+
ɛ
r
1
(
θ
)
+
ε
2
r
2
(
θ
)
+
O
(
ε
3
)
,
where r
0, r
1
(
θ
), r
2
(
θ
) are given in the following:
(42)
r
0
=
2
,
r
1
(
θ
)
=
1
35
sin
(
26
+
70
2
cos
⁡
5
θ
-
105
2
cos
⁡
3
θ
g
h
g
h
h
h
-
22
cos
⁡
2
θ
-
40
cos
⁡
6
θ
+
36
cos
⁡
4
θ
)
,
r
2
(
θ
)
=
-
914
35
cos
⁡
3
θ
+
14858
175
cos
⁡
5
θ
+
17408
105
cos
⁡
9
θ
-
7188
49
cos
⁡
7
θ
+
192
7
cos
⁡
13
θ
-
1168
11
cos
⁡
11
θ
-
208
49
2
cos
⁡
14
θ
-
258
35
2
cos
⁡
12
θ
+
12312
175
2
×
cos
⁡
10
θ
-
8367
70
2
cos
⁡
8
θ
+
1037
14
2
cos
⁡
6
θ
-
927
140
2
cos
⁡
4
θ
-
269
35
2
cos
⁡
2
θ
+
1130351
940800
2
.In Figure2 we illustrate the shape of the limit cycle Γ
1 of the system (38) by using formula (41) for the values ɛ
=
1
/
20. The periodic orbit u
2
+
v
2
=
2 of system (38) for ɛ
=
0 is drawn in solid line, and the limit cycle Γ
1 is drawn in dash line. We have also plotted the limit cycle Γ
1 for the value ɛ
=
1
/
20 by using the Runge-Kutta method in Figure 2. The close curve obtained numerically coincides with the one obtained analytically and we cannot distinguish between them with the eyes.Figure 2
The periodic orbitu
2
+
v
2
=
2 for ɛ
=
0 and the analytical and numerical approximations of the limit cycle for ɛ
=
1
/
20 of system (38).


### 4.2. The Shape of Limit Cycle of the Original Perturbed System (36)

In this subsection, we give the analytic expansion of the limit cycle of perturbed system (36) to the second order of ɛ.Rewrite the limit cycleΓ
1 of system (38) into the following parametric form:
(43)
u
=
r
(
θ
)
cos
⁡
θ
,
v
=
r
(
θ
)
sin
θ
,
where r
(
θ
) is given in (41).Thus from analytic transformation and time scale (37), corresponding to Γ
1, we obtain that limit cycle of the system (36) for ɛ
=
(
1
/
20
) is unstable and its parametric form is the following:
(44)
x
=
u
+
v
2
=
r
(
θ
)
cos
⁡
θ
+
(
r
(
θ
)
sin
θ
)
2
,
y
=
v
=
r
(
θ
)
sin
θ
.
The shape of limit cycle of the system (25) for ɛ
=
1
/
20 is plotted by using formula (44) in Figure 3.Figure 3
The periodic orbit forɛ
=
0 and the analytical and numerical approximations of the limit cycle for ɛ
=
1
/
20 of system (36).In Figure3, the periodic orbit (
x
-
y
2
)
2
+
y
2
=
2 of unperturbed system (36) for ɛ
=
0 is drawn in sold line, and the limit cycle of the perturbed system (36) for ɛ
=
1
/
20 is drawn in dash line.


## 4.1. The Shape of Limit Cycle of the Perturbed System (38)

In this subsection we start to compute the analytic expansion of the limit cycle of the perturbed system (38) to the second order of ɛ.First letu
=
r
cos
⁡
θ, v
=
r
sin
θ; then system (38) is transformed into the following polar coordinate form:
(39)
d
r
d
θ
=
ɛ
(
3
r
cos
⁡
θ
+
3
r
2
sin
2
θ
-
(
r
cos
⁡
θ
+
r
2
sin
2
θ
)
3
)
cos
⁡
θ
×
(
-
1
+
ɛ
(
(
r
cos
⁡
θ
+
r
2
sin
2
θ
)
2
)
sin
θ

-
3
cos
⁡
θ
-
3
r
sin
2
θ
+
(
cos
⁡
θ
+
r
sin
2
θ
)
ghhhhhhhhh
×
(
r
cos
⁡
θ
+
r
2
sin
2
θ
)
2
)
sin
θ
)
-
1
.
According to similar process in Section 3, we get
(40)
g
1
=
3
16
r
0
(
r
0
4
+
2
r
0
2
-
8
)
.By solvingg
1
=
0, we get the positive solution r
0
=
2.From(
d
(
-
2
π
r
0
g
1
)
/
d
r
0
)
|
r
0
=
2
=
-
9
2
π
<
0 and [9], we conclude that as ɛ
>
0, system (38) has a stable limit cycle, denoted by Γ
1, near the close curve u
2
+
v
2
=
2 on the phase plane.By applying the algorithm described in Section2 to system (39), we get the asymptotic expansion of stable limit cycle Γ
1 of system (38),
(41)
r
(
θ
)
=
r
0
+
ɛ
r
1
(
θ
)
+
ε
2
r
2
(
θ
)
+
O
(
ε
3
)
,
where r
0, r
1
(
θ
), r
2
(
θ
) are given in the following:
(42)
r
0
=
2
,
r
1
(
θ
)
=
1
35
sin
(
26
+
70
2
cos
⁡
5
θ
-
105
2
cos
⁡
3
θ
g
h
g
h
h
h
-
22
cos
⁡
2
θ
-
40
cos
⁡
6
θ
+
36
cos
⁡
4
θ
)
,
r
2
(
θ
)
=
-
914
35
cos
⁡
3
θ
+
14858
175
cos
⁡
5
θ
+
17408
105
cos
⁡
9
θ
-
7188
49
cos
⁡
7
θ
+
192
7
cos
⁡
13
θ
-
1168
11
cos
⁡
11
θ
-
208
49
2
cos
⁡
14
θ
-
258
35
2
cos
⁡
12
θ
+
12312
175
2
×
cos
⁡
10
θ
-
8367
70
2
cos
⁡
8
θ
+
1037
14
2
cos
⁡
6
θ
-
927
140
2
cos
⁡
4
θ
-
269
35
2
cos
⁡
2
θ
+
1130351
940800
2
.In Figure2 we illustrate the shape of the limit cycle Γ
1 of the system (38) by using formula (41) for the values ɛ
=
1
/
20. The periodic orbit u
2
+
v
2
=
2 of system (38) for ɛ
=
0 is drawn in solid line, and the limit cycle Γ
1 is drawn in dash line. We have also plotted the limit cycle Γ
1 for the value ɛ
=
1
/
20 by using the Runge-Kutta method in Figure 2. The close curve obtained numerically coincides with the one obtained analytically and we cannot distinguish between them with the eyes.Figure 2
The periodic orbitu
2
+
v
2
=
2 for ɛ
=
0 and the analytical and numerical approximations of the limit cycle for ɛ
=
1
/
20 of system (38).


## 4.2. The Shape of Limit Cycle of the Original Perturbed System (36)

In this subsection, we give the analytic expansion of the limit cycle of perturbed system (36) to the second order of ɛ.Rewrite the limit cycleΓ
1 of system (38) into the following parametric form:
(43)
u
=
r
(
θ
)
cos
⁡
θ
,
v
=
r
(
θ
)
sin
θ
,
where r
(
θ
) is given in (41).Thus from analytic transformation and time scale (37), corresponding to Γ
1, we obtain that limit cycle of the system (36) for ɛ
=
(
1
/
20
) is unstable and its parametric form is the following:
(44)
x
=
u
+
v
2
=
r
(
θ
)
cos
⁡
θ
+
(
r
(
θ
)
sin
θ
)
2
,
y
=
v
=
r
(
θ
)
sin
θ
.
The shape of limit cycle of the system (25) for ɛ
=
1
/
20 is plotted by using formula (44) in Figure 3.Figure 3
The periodic orbit forɛ
=
0 and the analytical and numerical approximations of the limit cycle for ɛ
=
1
/
20 of system (36).In Figure3, the periodic orbit (
x
-
y
2
)
2
+
y
2
=
2 of unperturbed system (36) for ɛ
=
0 is drawn in sold line, and the limit cycle of the perturbed system (36) for ɛ
=
1
/
20 is drawn in dash line.



---

*Source: 320406-2014-03-19.xml*