On the approximation of the limit cycles function

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 New


Introduction
We consider the planar autonomous system 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.
EJQTDE, 2007 No. 28, p. 1 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 dx dt = P (x, y, a), dy dt = Q(x, y, a), 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 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: 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.

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.
EJQTDE, 2007 No. 28, p. 2 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 consisting of limit cycles of system (1.1).

Remark 2.1
The invariance of M implies that the functional m satisfies the relation ∂m ∂x (x, y)P (x, y, m (x, y)) + ∂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 for a ∈ R, n ∈ N .Using polar coordinates x = r cos ϕ, y = r sin ϕ system (2.3) reads (2.4) It is easy to verify that the circle 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 EJQTDE, 2007 No. 28, p. 3 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 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)) Introducing the displacement function δ by 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 dr dϕ = r(r 2 − a), r(0) = r 0 > 0.
As solution we get The following relations can be verified Thus, K a is a simple limit cycle for system (2.3) with n = 1.

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 Hence, if we have for some x ∈ dom l, it holds 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 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 Analogously, after k−times differentiating the displacement function δ we obtain If we assume x (x, l(x)) = 0 we get under the validity of (3.2) 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.
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.

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) EJQTDE, 2007 No. 28, p. 6 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 = T (x) 2π τ = µτ with µ = 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 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 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 (µ k i , a k i ) defined by , k = 0, 1, ...
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

Example
We consider the Liénard system 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.
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.
EJQTDE, 2007 No. 28, p. 9 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.

Figure 1 :Figure 2 :
Figure 1: Graph of the function a = l 9 (x) respectively, as k tends to infinity.The entries of the matrix J k can be calculated by solving the initial value problem