Electronic Journal of Qualitative Theory of Differential Equations 2007, No. 28, 1-11; http://www.math.u-szeged.hu/ejqtde/

On the approximation of the limit cycles function

Leonid Cherkas Belorussian State University of Informatics and Radioelectronics Brovka Street 6, 220127 Minsk, Belarus Alexander Grin

Grodno State University, Ozheshko Street 22, 230023 Grodno, Belarus

Klaus R. Schneider Weierstrass Institute for Applied Analysis and Stochastics Mohrenstraße 39 10117 Berlin, Germany

Abstract

We consider planar vector fields depending on a real parameter. It is assumed that this vector field has a family of limit cycles which can be described by means of the limit cycles function l. We prove a relationship between the multiplicity of a limit cycle of this family and the order of a zero of the limit cycles function. Moreover, we present a procedure to approximate l(x), which is based on the Newton scheme applied to the Poincaré function and represents a continuation method. Finally, we demonstrate the effectiveness of the proposed procedure by means of a Liénard system.

Key words: family of limit cycles, multiple limit cycle, Liénard system AMS(MOS) Subject Classification: 34C05 34C07 65L99

1 Introduction

We consider the planar autonomous system

\frac{dx}{dt} = P (x, y, a), \frac{dy}{dt} = Q(x, y, a)  (1.1)

depending on the scalar parameter a for (x, y) ∈ Ω ⊂ R^{2}. Under certain conditions, the phase portrait of system (1.1) in Ω is determined by the so-called singular trajectories, namely equilibria, separatrices and limit cycles of (1.1) in Ω (see, e.g., [1, 11]). The most difficult problem in studying these singular trajectories is to localize multiple limit cycles and to estimate the number of limit cycles. This problem is still unsolved even in the case of polynomial systems, and it represents the second part of the famous 16-th problem of D. Hilbert [9, 12]. As already D. Hilbert indicated, the investigation of the dependence of limit cycles on parameters should play a fundamental role in solving the posed problem.

Our main goal is to exploit the parameter dependence in (1.1) for detecting multiple limit cycles. If we consider the parameter a as an additional state variable, then we get from (1.1) the extended system

\frac{dx}{dt} = P (x, y, a), \frac{dy}{dt} = Q(x, y, a), \frac{da}{dt} = 0.  (1.2)

We assume that this system has an invariant manifold of the form

a = m(x, y)  (1.3)

consisting of limit cycles of system (1.1), and that there exists a smooth segment

S := {(x, y) ∈ R^{2} : y = s(x), x_{0} ≤ x ≤ x_{1}}

in the phase plane such that any limit cycle of the family (1.3) intersects S transversally and there are no two distinct limit cycles of this family intersecting S in the same point. Then, the limit cycles of the family (1.3) can be uniquely characterized by the x-coordinate of their intersection point (x, s(x)) with the segment S. We denote by L(x) the limit cycle of the family (1.3) intersecting S in the point (x, s(x)). Now we introduce the function l: [x_{0}, x_{1}] → R by

l(x) := m(x, s(x)),  (1.4)

which associates L(x) with the corresponding parameter a. This function is called the limit cycles function. Under some conditions it coincides with the Andronov-Hopf function [2, 4, 7, 8]. In this note we will show that the order of a zero of the derivative of l(x) is related to the multiplicity of the corresponding limit cycle L(x). Moreover, we present a procedure to approximate l(x), which is based on the Newton scheme applied to the Poincaré function and represents a continuation method. Finally, we demonstrate the effectiveness of the proposed procedure by means of a Liénard system. For all calculations we used the software MATHEMATICA 5.2.

2 Prelimineries

Let Ω be some simply connected region in R^{2}, and let I ⊂ R be some interval. We consider system (1.1) in Ω for a ∈ I under the following smoothness assumption.

(A_{1}). P and Q are n-times continuously differentiable with respect to x and y, and continuously differentiable with respect to a in Ω × I with n ≥ 1.

System (1.1) defines the planar vector field f := (P, Q) on Ω. We denote f as smooth if assumption (A_{1}) is satisfied.

An isolated periodic solution of (1.1) which is no stationary solution is called a limit cycle.

Our key assumptions with respect to the extended system (1.2) are the following ones.

(A_{2}). There is a smooth functional m: dom m ⊂ Ω → I such that system (1.2) has an invariant manifold of the form

M = {(x, y, a) ∈ R^{3} : a = m(x, y), (x, y) ∈ dom m},  (2.1)

consisting of limit cycles of system (1.1).

Remark 2.1 The invariance of M implies that the functional m satisfies the relation

\frac{∂m}{∂x}(x, y)P (x, y, m (x, y)) + \frac{∂m}{∂y}(x, y)Q(x, y, m (x, y)) ≡ 0 for (x, y) ∈ dom m. (2.2)

(A_{3}). There is a smooth segment S := {(x, y) ∈ Ω : y = s(x), x_{0} ≤ x ≤ x_{1}} intersecting transversally all limit cycles belonging to M. There are no two different limit cycles of M intersecting S in the same point.

Definition 2.1 The function l: [x_{0}, x_{1}] → I defined by (1.4) is called the limit cycles function.

We note that l is a smooth function. In case that the manifold M is connected with the bifurcation of a limit cycle from an equilibrium point and that the parameter a rotates the vector field f, the limit cycles function l coincides with the Andronov-Hopf function (see [2, 4, 7, 8]).

To illustrate the introduced concepts we consider the system

\begin{matrix}\frac{dx}{dt} = −y + x(x^{2} + y^{2} − a)^{n} ≡ P (x, y, a), \\ \frac{dy}{dt} = x + y(x^{2} + y^{2} − a)^{n} ≡ Q(x, y, a) \end{matrix}  (2.3)

for a ∈ R, n ∈ N. Using polar coordinates x = r cos φ, y = r sin φ system (2.3) reads

\begin{matrix}\frac{dr}{dt} = r(r^{2} − a)^{n}, \\ \frac{dφ}{dt} = 1. \end{matrix}  (2.4)

It is easy to verify that the circle

K_{a} := {(x, y) ∈ R^{2} : x^{2} + y^{2} = a, a > 0}

is a limit cycle of system (2.3). Thus, we have m(x, y) := x^{2} + y^{2} with dom m := R^{2}\{(0, 0)} for any n, and the invariant manifold M is connected with the bifurcation of the limit cycle K_{a} from the origin when a passes the value zero for increasing a. It is also obvious that the positive x−axis intersects the limit cycle K_{a} transversally. Hence, we have s(x) ≡ 0 and the limit cycles function l reads l(x) ≡ x^{2} and is defined for x > 0. From (2.3) we get the relation

P (x, y, a) \frac{∂Q}{∂a}(x, y, a) − Q(x, y, a) \frac{∂P}{∂a} (x, y, a) = n(x^{2} + y^{2})(x^{2} + y^{2} − a)^{n−1}

that implies that system (2.3) represent for n = 1 a rotated vector field, and thus the limit cycles function l coincides with the Andronov-Hopf function.

Since for x ∈ [x_{0}, x_{1}] the limit cycle L(x) of system (1.1) for a = l(x) intersects the segment S transversally, the Poincaré function π(x, a) (first return function) is defined for x ∈ [x_{0}, x_{1}] and a near l(x). From the property that a fixed point of the Poincaré function corresponds to a periodic solution of (1.1), we have according to the definition of l

π(x, l(x)) ≡ x for x ∈ [x_{0}, x_{1}].  (2.5)

Introducing the displacement function δ by

δ(x, a) := π(x, a) − x,  (2.6)

then an isolated root of the displacement function corresponds to a limit cycle of (1.1).

Definition 2.2 A limit cycle is called a limit cycle of multiplicity k if the corresponding isolated root of the displacement function δ has multiplicity k.

To determine the multiplicity of the limit cycle K_{a} of system (2.3) for n = 1 we consider the corresponding equation in polar coordinates and study the initial value problem

\frac{dr}{dφ} = r(r^{2} − a), r(0) = r_{0} > 0.

As solution we get

r(φ, r_{0}, a) = r_{0}\sqrt{\frac{a}{r_{0}^{2} − (r_{0}^{2} − a)^{2φa}}} .

Thus, we have

δ(r_{0}, a) = r(2π, r_{0}, a) − r_{0} = r_{0} (\sqrt{\frac{a}{r_{0}^{2} − (r_{0}^{2} − a)^{4πa}}} − 1) .

The following relations can be verified

δ(\sqrt{a}, a) = 0, δ_{r_{0}}(\sqrt{a}, a) = −1.

Thus, K_{a} is a simple limit cycle for system (2.3) with n = 1.

3 Multiple limit cycles

In this section we derive an result about the zeros of the derivative of the limit cycles function l and the existence of multiple limit cycles. From (2.5) and (2.6) we get

δ(x, l(x)) = π(x, l(x)) − x ≡ 0 for x ∈ dom l.

Differentiating this relation we obtain

π_{x}(x, l(x)) − 1 + π_{a}(x, l(x)) l′(x) ≡ 0.  (3.1)

Hence, if we have

π_{a}(x̃, l(x̃)) ≠ 0  (3.2)

for some x̃ ∈ dom l, it holds

l′(x̃) = π_{a}^{−1}(x̃, l(x̃)) (1 − π_{x}(x̃, l(x̃))) = −π_{a}^{−1}(x̃, l(x̃)) δ_{x}(x̃, l(x̃)).  (3.3)

Thus, under the condition (3.2), the property that for a = l(x̃) the limit cycle L(x̃) of system (2.3) represents a simple limit cycle is equivalent to the condition l′(x̃) ≠ 0. Further, differentiating relation (3.1) we obtain

π_{xx}(x, l(x)) + 2π_{xa}(x, l(x)) l′(x) + π_{aa}(x, l(x)) (l′(x))^{2} + π_{a}(x, l(x)) l′′(x) ≡ 0.

If we assume δ_{x}(x̃, l(x̃)) = 0, which implies l′(x̃) = 0, then we get from (3.3) under the condition (3.2) the relation

l′′(x̃) = −π_{a}^{−1}(x̃, l(x̃))δ_{xx}(x̃, l(x̃)).

Analogously, after k−times differentiating the displacement function δ we obtain

π_{x}^{(k)}(x, l(x)) + f_{1}(x) l′(x) + f_{2}(x) l′′(x) + ... + π_{a}(x, l(x)) l^{(k)}(x) ≡ 0.

If we assume

δ_{x}(x̃, l(x̃)) = δ_{xx}(x̃, l(x̃)) = ... = δ_{x}^{(k−1)}(x̃, l(x̃)) = 0

we get under the validity of (3.2)

l′(x̃) = l′′(x̃) = ... = l^{(k−1)}(x̃) = 0, l^{(k)}(x̃) = −π_{a}^{−1}(x̃, l(x̃))δ_{x}^{(k)}(x̃, l(x̃)).

Hence, we have the result

Theorem 3.1 Suppose that the assumptions (A_{1}) − (A_{3}) hold and that for some x̃ ∈ dom l the relation (3.2) is valid. Then a necessary and sufficient condition for L(x̃) to be a limit cycle of multiplicity k, k ≤ n, is that x̃ is a root of the function l′(x) of order k − 1.

To illustrate that theorem we consider the system

\begin{matrix}\frac{dx}{dt} = −y + x((x^{2} + y^{2} − 1)^{2} − a) ≡ P (x, y, a), \\ \frac{dy}{dt} = x + y((x^{2} + y^{2} − 1)^{2} − a) ≡ Q(x, y, a) \end{matrix}  (3.4)

for a ∈ R. Using polar coordinates x = rcos φ, y = rsin φ, system (3.4) reads

\begin{matrix}\frac{dr}{dt} = r((r^{2} − 1)^{2} − a), \\ \frac{dφ}{dt} = 1. \end{matrix}  (3.5)

The extended system to (3.4) has for a ≥ 0 the invariant manifold a = (x^{2} +y^{2} −1)^{2} consisting of limit cycles of system (3.4). The positive x-axis intersects these limit cycles transversally. Thus, the limit cycle function l reads l(x) := (x^{2} −1)^{2}. It has a simple minimum at x = 1, where l takes the value a = l(1) = 0. The corresponding limit cycle L(1) reads x^{2} + y^{2} = 1. To determine its multiplicity we consider the initial value problem

\frac{dr}{dφ} = r(r^{2} − 1)^{2}, r(0) = r_{0}

that has the implicit solution

\frac{1}{2(r_{0}^{2} − 1)} − \frac{1}{2(r^{2}(φ, r_{0}) − 1)} + ln (\frac{r(φ, r_{0})}{\sqrt{r^{2}(φ, r_{0}) − 1}}) − ln (\frac{r_{0}}{\sqrt{r_{0}^{2} − 1}}) = φ.

Using this relation we get

δ(1) = 0, δ_{r_{0}}(1) = 0, δ_{r_{0}r_{0}}(1) ≠ 0.

This can be verified e.g. by means of the software MATHEMATICA 5.2. Hence, the limit cycle x^{2} + y^{2} = 1 of system (3.4) with a = 0 has the multiplicity 2.

4 A numerical procedure to approximate the limit cycles function

In what follows we describe a numerical method to approximate to given x ∈ dom l the value of l(x). For this purpose we denote by (x̄(t, x, y, a), ȳ(t, x, y, a)) the solution of (1.1) satisfying x̄(0, x, y, a) = x, ȳ(0, x, y, a) = y. Under the hypotheses (A_{1}) − (A_{3}), to x ∈ dom l there is a unique parameter value a = l(x) and a unique primitive period T = T (x) > 0 such that the orbit of system (1.1) starting at the point (x, s(x)) is the limit cycle L(x) with primitive period T (x). Hence, the system of equations

x̄(T, x, s(x), a) − x = 0, ȳ(T, x, s(x), a) = s(x)  (4.1)

can be used to determine a and T as functions of x.

Our basic idea is to apply a continuation method to compute l(x) for x ∈ dom l. We assume that for some point x = x_{0} the value a_{0} = l(x_{0}) as well as the limit cycle L(x_{0}) is known and use these facts for our continuation method. For the following we assume that x_{0} and the corresponding value l(x_{0}) is connected with the appearance of a simple limit cycle from the equilibrium point (x_{e}, y_{e}) of system (1.1) (Andronov-Hopf bifurcation) when a passes the value l(x_{0}). From the theory of Andronov-Hopf bifurcation we know that the primitive period T_{0} of the bifurcating limit cycle satisfies T_{0} = 2π/b_{0}, where b_{0} is the imaginary part of the critical eigenvalues of the Jacobian of the right hand side of (1.1) for a = a_{0} = l(x_{0}) at the equilibrium (x_{e}, y_{e}). Moreover, we can use a segment of the straight line y = y_{e} as our segment S. In the next step we replace x_{0} by x_{0} + h, where h is a small positive number and apply Newton’s iteration scheme to determine a and T from system (4.1), where we use a_{0} and T_{0} = 2π/b_{0} as initial guess. In the practical realization we use a slight modification of this procedure. By means of the transformation t = \frac{T (x)}{2π} τ = μτ with μ = \frac{T (x)}{2π} we introduce a new time τ such that the primitive period of the limit cycle L(x) is always 2π, but then we have to determine the additional parameter μ. Using the new time τ, we get from (1.1) the system

\frac{dx}{dτ} = μP (x, y, a), \frac{dy}{dτ} = μQ(x, y, a).  (4.2)

If we denote by (x̃(τ, x, y_{e}, μ, a), ỹ(τ, x, y_{e}, μ, a)) the solution of system (4.2) satisfying

x̃(0, x, y_{e}, μ, a) = x, ỹ(0, x, y_{e}, μ, a) = y_{e},

then the system of equations which analogously to (4.1) determines the parameters a and μ has the form

\begin{aligned}φ_{1}(x, a, μ) & ≡ & x̃(2π, x, y_{e}, μ, a) − x = 0, \\ φ_{2}(x, a, μ) & ≡ & ỹ(2π, x, y_{e}, μ, a) = y_{e}. \end{aligned}  (4.3)

Suppose we have determined to the sequence x_{1}, ..., x_{i−1} the values μ_{1}^{∗}, ..., μ_{i−1}^{∗}, a_{1}^{∗}, ..., a_{i−1}^{∗} approximating the corresponding values μ(x_{1}) = T (x)/2π, ..., μ(x_{i−1}) = T (x_{i−1})/2π, and a(x_{1}) = l(x_{1}), ..., a(x_{i−1}) = l(x_{i−1}). In order to determine to x_{i} the corresponding approximating values (μ_{i}^{∗}, a_{i}^{∗}) we apply Newton’s method to system (4.3) yielding the sequence (μ_{i}^{k}, a_{i}^{k}) defined by

\binom{a_{i}^{k+1}}{μ_{i}^{k+1}} = \binom{a_{i}^{k}}{μ_{i}^{k}} − J_{k}^{−1}(x_{i})\binom{φ_{1}(x_{i}, μ_{i}^{k}, a_{i}^{k})}{φ_{2}(x_{i}, μ_{i}^{k}, a_{i}^{k},)}, k = 0, 1, ...

where a_{i}^{0} = a_{i−1}^{∗}, μ_{i}^{0} = μ_{i−1}^{∗},

J_{k}(x_{i}) = ( \begin{matrix}\frac{∂x̃}{∂a}(2π, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) & \frac{∂x̃}{∂μ}(2π, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) \\ \frac{∂ỹ}{∂a}(2π, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) & \frac{∂ỹ}{∂μ}(2π, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) \end{matrix} ) .

Remark 4.1 Under the assumption that the Jacobian matrix J_{k}(x_{i}) is invertible and that the difference |x_{i} − x_{i−1}| is sufficiently small for any i, the sequences {μ_{i}^{k}}, {a_{i}^{k}}, converge to μ_{i}^{∗}, a_{i}^{∗}, respectively, as k tends to infinity.

The entries of the matrix J_{k} can be calculated by solving the initial value problem

\frac{dx}{dτ} = μP (x, y, a), \frac{dy}{dτ} = μQ(x, y, a),

\frac{d(\frac{∂x}{∂a})}{dτ} = μ (\frac{∂P}{∂a} + \frac{∂P}{∂x} \frac{∂x}{∂a} + \frac{∂P}{∂y} \frac{∂y}{∂a}) ,

\frac{d(\frac{∂y}{∂a})}{dτ} = μ (\frac{∂Q}{∂a} + \frac{∂Q}{∂x} \frac{∂x}{∂a} + \frac{∂Q}{∂y} \frac{∂y}{∂a}) ,

\frac{d(\frac{∂x}{∂μ})}{dτ} = P + μ (\frac{∂P}{∂x} \frac{∂x}{∂μ} + \frac{∂P}{∂y} \frac{∂y}{∂μ}) ,

\frac{d(\frac{∂y}{∂μ})}{dτ} = Q + μ (\frac{∂Q}{∂x} \frac{∂x}{∂μ} + \frac{∂Q}{∂y} \frac{∂y}{∂μ}) ,

x(0, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) = x_{i}, y(0, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) = y_{e},

\frac{∂x}{∂μ}(0, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) = \frac{∂y}{∂μ}(0, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) = 0,

\frac{∂x}{∂a}(0, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) = \frac{∂y}{∂a}(0, x_{i}, y_{e}, μ_{i}^{k}, a_{i}^{k}) = 0.

5 Example

We consider the Liénard system

\begin{matrix}\frac{dx}{dt} = y + ax − \frac{1024x^{11} − 2816x^{9} + 2816x^{7} − 12325x^{5} + 220x^{3}}{1024} ≡ P (x, y, a), \\ \frac{dy}{dt} = −x ≡ Q(x, y, a) \end{matrix} (5.1)

depending on the real parameter a. This system has the following properties.

(i). It has a unique equilibrium point in the finite part of the phase plane, namely the origin, which does not depend on the parameter a.

(ii). Let α(a) ± iβ(a) be the eigenvalues of the linearization of system (5.1) for small a at the origin. We have

α(0) = 0, α′(0) = 1/2, β(0) > 0.

The corresponding Lyapunov number α_{1} (see, e.g., [6, 10]) is negative.

(iii). By (5.1) we have

Q(x, y)P′_{a}(x, y, a) ≡ −x^{2},

that is, the parameter a rotates the vector field f corresponding to (5.1) (see [6, 10]).

(iv) Any trajectory of (5.1) intersects the positive x-axis transversally.

From (ii) we get that Hopf bifurcation takes place when the parameter a crosses zero, where exactly one limit cycle Γ(a) bifurcates from the origin for increasing a. This limit cycle is simple and orbitally asymptotically stable. Moreover, by property (iv), for small a this limit cycle intersect the positive x-axis transversally. From the theory of rotated vector fields [6, 10] we get that Γ(a_{1}) and Γ(a_{2}) do not intersect for small a_{1} and a_{2} if a_{1} ≠ a_{2}. Thus, the family of bifurcating limit cycles {Γ(a)} can be characterized also by the x-coordinates of their intersection point with the positive x-axis, that means by the limit cycle function l, which coincides in the present case with the Andronov-Hopf function. In a given small neighborhood of the origin the number of limit cycles can change by the occurrence of multiple limit cycles. In our case, we restrict to the neighborhood x^{2} + y^{2} ≤ 1.2 and use the limit cycle function l to investigate the appearance of multiple limit cycles. By applying the numerical procedure described in the section before we compute l(x) for x = 0.01 + j 0.025, j = 1, .., 44. Then we approximate the obtained set of points by some polynomial l_{N} of order N. For this purpose we apply the method of least squares. We choose N in such a way that number of extrema of the polynomial l_{N} in the interval [0.01, 1.085] does not change if we increase N. By this way, we have got the result

l_{9}(x) =  0.154224x^{2} + 0.124769x^{3} − 1.63087x^{4}
+ 3.22805x^{5} − 5.36163x^{6} + 8.75029x^{7} − 7.9267x^{8} + 2.6726x^{9}.

By means of Sturm’s chains we can prove that the polynomial l_{9} has in the interval [0.01, 1.085] exactly four extrema: the points A(0.4522, 0.0121) and C(0.8726, 0.0111) are minima, the points B(0.6952, 0.0101) and D(0.9758, 0.0105) are maxima. The graph of the function a = l_{9}(x) is represented in Fig. 1. From these investigations we can conclude that system (5.1) has at least five limit cycles for 0.0105003 ≤ a ≤ 0.0111532.

Figure 1: Graph of the function a = l_{9}(x)

Figure 2: Ovals defined by Ψ(x, y) = 0

To check the quality of the obtained numerical results for a = 0.0107 we apply the method of Dulac-Cherkas function Ψ(x, y) as described in [3, 4, 5, 7]. We construct the function Ψ(x, y) in the form of a polynomial of 100-th degree. The curve Ψ(x, y) = 0 defines five nested ovals surrounding the unique finite equilibrium at the origin (Fig. 2) which are intersected transversally by the trajectories of system (5.1) for a = 0.0107. Using the techniques from [7] we can establish that system (5.1) for a = 0.0107 has exactly five limit cycles. We note that it is possible to construct a Dulac-Cherkas function which depends on the parameter a and to prove that system (5.1) has for any a not more than five limit cycles.

6 Acknowledgement

The second author acknowledgements the financial support by DAAD and the hospitality of the Weierstrass Institute for Applied Analysis and Stochastics in Berlin.

References

[1] A. A. Andronov, E. A. Leontovich, I. I. Gordon, A. G. Maier, Qualitative theory of second order dynamical systems, John Wiley and Sons, New York, 1973.

[2] L. A. Cherkas, The question on the analyticity of a manifold determining limit cycles, Differential Equations 18 (1982), 839-845.

[3] L. A. Cherkas, Dulac function for polynomial autonomous systems on a plane, Differential Equations 33 (1997), 692-701.

[4] L. A. Cherkas, J. C. Artes, J. Llibre, Quadratic systems with limit cycles of normal size, Bul. Acad. Stiinte Rep. Mold., Matematica. 41 (2003), 31-46.

[5] L. A. Cherkas, A. A. Grin, Algebraic aspects of finding a Dulac function for polynomial autonomous systems on the plane, Differential Equations 37 (2001), 411-417.

[6] F. Dumortier, J. Llibre, J. C. Artes, Qualitative theory of planar differential systems, Universitext. Springer, 2006.

[7] A. A. Grin, L. A. Cherkas, Dulac function for Liénard systems (in Russian), Transactions Inst. Math. Belorus. Nat. Acad. Sci., Minsk, 4 (2000), 29-38.

[8] A. A. Grin, L. A. Cherkas, Extrema of the Andronov-Hopf function of a polynomial Liénard system, Differential Equations 41 (2005), 50-60.

[9] D. Hilbert, Mathematical Problems. Reprinted english translation: Bull. Amer. Math. Soc. 37 (2000), 407-436.

[10] L. Perko, Differential equations and dynamical systems, Texts Appl. Math. 7, Springer, 2001.

[11] L.P. Shilnikov, A.L. Shilnikov, D.V. Turaev, L.O. Chua, Methods of qualitative theory in nonlinear dynamics. Part I, Part II, World Scientific Series on Nonlinear Science. Series A. 4,5. Singapore, 1998, 2001.

[12] S. Smale, Mathematical problems for the next century, Mathematical Intelligencer 20 (1998), 7-15.

(Received April 11, 2007)