Localization of Liénard equation’s attractors and cycles^{1}

G.A.Leonov^{2} and E.A.Strigun leonov@math.spbu.ru yelena strigun@mail.ru

1 Introduction

There is a well known M.Cartwright theorem ([5], [24]) concerning the uniform by parameter μ boundedness of the y coordinate of solutions of Liénard equation

\frac{dy}{dt} = −μ[F (y) − E(t)] − x,  \frac{dx}{dt} = y.  (i)

Here F (y), E(t) satisfy a Lipschitz condition, parameter μ is positive. After M. Cartwright this problem was intensively investigated by other researchers [17], [18], [22], [29], [36].

A very important special case of system (i) is the Van der Pol equation

\frac{dy}{dt} = −μ(\frac{y^{3}}{3} − y) − x,  \frac{dx}{dt} = y.  (ii)

In [5] an estimation for sufficiently large parameter μ > 0 of the systems’s cycle was obtained.

At that time Van der Pol found a graphic solution for μ = 10. Then investigations of the equation for large values of the parameter were undertaken in [13], [16], [23].

Ph. le Colbeiller proved that the amplitude ρ(μ) (maximum absolute value of y) of the Van der Pol equation’s cycle is equal to 2 + o(1) for μ → ∞.

An estimation for sufficiently small μ > 0 was given by W. van Horssen [19]: ρ(μ) = 2 + (2/96)μ^{2} + O(μ^{3}).

An uniform estimation was obtained by P. Alsholm [1] : max |y| < 2.8025 for any μ ≠ 0. K. Odani improved this result to 2.5425 in [32] then to 2.3439 and 2.3233 in [33] and [34]. By computer experiment using the Runge-Kutta method he made the conjecture that 2 < ρ(μ) < 2, 0235 for any μ ≠ 0.

In this paper the method of constructions of comparison systems with discontinuous right-hand sides is applied to the estimation of Liénard system’s attractor. This method was introduced in [25]-[27] and it is based on comparison systems method of Chaplygin-Kamke [6], [20].

Further we deal with a special case of the Liénard system – Van der Pol equation(ii). For this equation a more complex comparison system is constructed. It uses the initial data previously founded from an application of a global theorem. By the computer assisted estimation we obtained the following localization: ρ(μ) < 2.015 for any μ > 0. It improves the conjecture of K. Odani.

The second important application of the method is an investigation of the equation (iii):

\frac{dy}{dt} = −μF (y) − x,  \frac{dx}{dt} = y + a.  (iii)

Here F (y) = \frac{y^{3}}{3} − y and μ, a are positive numbers.

An analysis of a solutions of the equation with respect of parameter a is one of the most interesting application of so called nonstandard analysis to the theory of differential equation. For equation (iii) there is revealed a phenomena of so called ”French ducks” [2], [37], [38], [9], [10], [30], [11], [3], [39]. Here is a numerical construction of the solution which makes the result of Cartier [4] more precise.

2 Main results

Let us consider the system

\frac{dy}{dt} = −μ[F (y) − E(t)] − x,  \frac{dx}{dt} = y,  (1)

where F (y), E(t) satisfy Lipschitz condition and μ is a positive number.

1Based on the lecture “Localization of Liénard attractor” hold at the DFG-Workshop “Dimensions of attractors”, St.

Petersburg, Russia, October 11-12, 1999

2Partly supported by the Russian Research Fund 98-01-01005

Further we will assume that for some positive numbers α and k the following inequality holds

\frac{F (y) − E(t)}{y} > \frac{αy − k sign y}{y} ,  ∀ t ∈ R^{1}, ∀ y ≠ 0.  (2)

Assumption (2) is usual and traditional for Liénard system (1) (see [7], [24], [35]).

First let us consider the case αμ ≥ 2. Here we consider the positive semi-trajectory of the system

\frac{dy}{dt} = −μαy − x + μk,  \frac{dx}{dt} = y  (3)

with initial data y(0) = 0, x(0) = −μk and the positive semi-trajectory of the system

\frac{dy}{dt} = −μαy − x − μk,  \frac{dx}{dt} = y  (4)

with the initial data y(0) = 0, x(0) = μk.

Solutions G_{1}(x, μk) and G_{2}(x, μk) correspond to the first-order equations

G \frac{dG}{dx} = −μαG − x + μk,  (5)

G \frac{dG}{dx} = −μαG − x − μk,  (6)

respectively. Further we deal with the set

Ω(α, k) = {x ∈ [−μk, μk], G_{2}(x, μk) ≤ y ≤ G_{1}(x, μk)}.

In the case of αμ < 2 we determine the set Ω(α, k) in a slightly different way. It is well-known [21] that the system

\frac{dy}{dt} = −μαy − x + μk sign y,  \frac{dx}{dt} = y  (7)

has a unique limit cycle x(t), y(t) with the initial data x(0) = ρ, y(0) = 0. Here

ρ = \frac{1 + exp (\frac{−λπ}{ω} )}{1 − exp (\frac{−λπ}{ω} )} μk,

where λ = αμ/2, ω = \sqrt{1 − α^{2}μ^{2}/4} .

Solutions of the discontinuous system (7) are defined in the sense of [12]. The solution G_{1}(x, ρ) of equation (5) corresponds to the part of the limit cycle which is situated in the phase half-space {x ≥ 0}. And the solution G_{2}(x, ρ) of equation (6) corresponds to the other part of the limit cycle of system (7) which is situated in the phase half-space{x ≤ 0}.

We define the set Ω(α, k) in the following way:

Ω(α, k) = {x ∈ [−ρ, ρ], G_{2}(x, ρ) ≤ y ≤ G_{1}(x, ρ)}.

Let us recall that a set

Φ(t, t_{0}) = ( \begin{matrix}x(t, t_{0}, Φ_{0}) \\ y(t, t_{0}, Φ_{0}) \end{matrix} ) , Φ(t_{0}, t_{0}) = Φ_{0},

is called global attracting set if for any ε > 0, x_{0} ∈ R^{1}, y_{0} ∈ R^{1} there exists T (ε, x_{0}, y_{0}) such that the solution (x(t, t_{0}, x_{0}, y_{0}), y(t, t_{0}, x_{0}, y_{0})) belongs to ε-neighborhood of the set Φ(t, t_{0}) for any t ≥ T (ε, x_{0}, y_{0}).

The set Φ(t, t_{0}) is called the minimal attracting set if there is no proper subset of Φ(t, t_{0}) which has the property of global attraction.

The minimal global attracting set Φ(t, t_{0}) is said to be the global attractor of system (1).

Theorem 2.1 The global attractor Φ(t, t_{0}) of system (1) belongs to Ω(α, k).

Proof. Let inequality (2) be fulfilled for some α > 0 and k = k_{0}. It obviously holds for all k ≥ k_{0} either. But then for any point (x_{0}, y_{0}) of the set R^{2} \ Ω(α, k_{0}) there exists a number k ≥ k_{0} such that (x_{0}, y_{0}) belongs to the boundary of Ω(α, k). Therefore we have a family of closed curves covering the set R^{2} \ Ω(α, k_{0}).

Our purpose is to demonstrate that all these curves are without contact everywhere except for points {y = 0, x ∈ R^{1}} and trajectories of system (1) go through these curves from outside to inside. For proving this fact the Chaplygin-Kamke comparison principle [6], [20], [28] and inequality (2) will be used:

\frac{dy}{dx} = \frac{−μ(F (y) − E(t)) − x}{y} < \frac{−μ(αy − k signy) − x}{x},

∀ x ∈ R^{1}, ∀ y ≠ 0, ∀ t ∈ R^{1}.

It follows from the comparison principle that for solutions G_{j}(x) corresponding to trajectories of the system

\frac{dy}{dt} = −μ(αy − k sign y) − x, \frac{dx}{dt} = y

and for the solution y(t), x(t) of system (1) the following inequality in a point t = t_{0}, x_{0} = x(t_{0}), y_{0} = y(t_{0}) = G_{j}(x_{0}) ≠ 0 is fulfilled:

\frac{dy}{dx} < \frac{dG_{j}}{dx} .

Hence curves y = G_{j}(x) are without contact regarding to a vector field of system (1). It follows from the fact that the set of curves are without contact almost everywhere that the theorem is true.

Let us consider the set

Ω_{0} = ⋂_{α,k} Ω(α, k),

where the intersection is taken for all parameters α, k satisfying condition (2).

Corollary 2.1 The global attractor of system (1) belongs to Ω_{0}.

The Van der Pol equation

Let us consider the Van der Pol equation, i.e. the case of E(t) = 0 and F (y) = \frac{y^{3}}{3} −y. Here condition (2) is satisfied if

k = \frac{2}{3} (α + 1)^{3/2} + ε,  (8)

where ε is any positive number. If α = 2/μ then the set Ω_{0} belongs to the strip

\{|x| ≤ \frac{2μ}{3} (\frac{2}{μ} + 1)^{3/2} , y ∈ R^{1}\} .  (9)

If αμ < 2, taking in account (8) we obtain the following inclusion:

Ω_{0} ⊂ \{|x| ≤ \frac{2μ}{3} min_{α∈(0,2/μ)} [(α + 1)^{3/2} \frac{1 + exp(−λπ/ω)}{1 − exp(−λπ/ω)}] , y ∈ R^{1}\} ,  (10)

where

τ(α) = \frac{1}{ω} arctg (\frac{ω}{λ} ) .

Now a more complex system of comparison will be constructed. Let us assume that for the function F (y) = \frac{y^{3}}{3} − y from the Van der Pol equation

\frac{dy}{dt} = −μF (y) − x,  y = \frac{dx}{dt},  (11)

and for some real numbers ν_{i} the following inequalities are fulfilled:

\frac{F (y)}{y} > \frac{ν_{i} sign y}{y} ,  ∀ t ∈ R^{1}, ∀ y ∈ Γ_{i},  (12)

where Γ_{i} are intervals of the form

Γ_{i} = \{ \begin{matrix}(γ_{i}, γ_{i+1}] & for & y > 0, \\ [−γ_{i+1}, −γ_{i}) & for & y < 0, \end{matrix}

where i = 0, …, N and N is some natural number. A sequence of γ_{i} contains numbers 0, 1 and defines a partition of the interval [0, γ_{N}] (γ_{N} ≥ 1) with some step S ≤ 1.

Let us construct such ν_{i} that inequalities (12) are valid. For |y| < γ_{i} = 1 the function F decreases and for |y| > γ_{i} = 1 it increases. So if we define ν_{i} as follows

ν_{i} = \{ \begin{matrix}F (γ_{i}) & for & γ_{i} ≤ 1, \\ F (γ_{i−1}) & for & γ_{i} > 1, \end{matrix}  (13)

then the function F will be bounded from above in the best way. For more clearness see Fig. 1, where the graph of the function F(y) is represented. Vertical lines correspond to the value of γ_{i} and horizontal ones correspond to appropriate ν_{i}.

Let us consider the equation

F \frac{dF}{dx} = −μf(F) − x,  (14)

where

f(y) = \{ \begin{matrix}−ν_{N} & for & y ∈ [−γ_{N}, −γ_{N−1}), \\ & ⋯ \\ −ν_{1} & for & y ∈ [−γ_{1}, 0), \\ ν_{1} & for & y ∈ (0, γ_{1}], \\ & ⋯ \\ ν_{N} & for & y ∈ (γ_{N−1}, γ_{N}]. \end{matrix}

for any ν_{i} < 0 and y.

We deal with solutions F_{1}(x) and F_{2}(x) of equation (14) which are situated in half-planes {F ≥ 0} and {F ≤ 0} respectively with initial data F_{1}(B_{0}) = 0. For parameter μ ≥ 1 the value B_{0} is taken from (9) as a maximum of |x| and for μ < 1 it is defined from (10). Let F_{2}(B_{1}) = 0, where B_{1} is such that F_{1}(x) > 0 for x ∈ (B_{0}, B_{1}). The value B_{2} is taken such that F_{2}(B_{2}) = 0 and F_{2}(x) < 0 for x ∈ (B_{1}, B_{2}) (see Fig. 2).

Let us define a set Ψ(ν_{i}) as follows

Ψ(ν_{i}) = \{ \begin{matrix}x ∈ [B_{2}, B_{1}], & F_{2}(x) ≤ y ≤ F_{1}(x), \\ x ∈ [B_{0}, B_{2}], & y ≥ 0 \end{matrix} \}.

Theorem 2.2 The global attractor Φ(t, t_{0}) of system (11) belongs to Ψ(ν_{i}).

Proof. It follows from Bendixon’s theorem [9] that all trajectories tend helically to the attracting set which in our case is a cycle around the equilibrium point (0,0). It follows from the choice of B_{0} and from theorem 1 that the cycle belongs to the strip B_{0} < x < −B_{0}. Trajectories wind on the cycle from the outside, hence |B_{2}| < |B_{0}|.

Let the Van der Pol equation’s cycle intersect the x-axis in points B_{0}^{1} and B_{1}^{1} (their existence is known from[9]) and B_{0}^{1} < B_{1}^{1}. In the same way from the choice of B_{0} and theorem 1 we can conclude that B_{0}^{1} > B_{0}. Trajectories wind on the cycle from the outside so B_{1}^{1} < B_{1} and further B_{2} < B_{0}^{1}. So the cycle of the Van der Pol equation belongs to the strip B_{2} < x < B_{1} (see Fig.3).

It is similar to the proof of theorem 1 that there is exists a family of closed curves covering the set R^{2} \ Ψ(ν_{i}).

Everywhere except for points {y = 0, x ∈ R^{1}} all these curves are without contact and trajectories of system (11) go through these curves from outside to inside(see Fig. 4).

Indeed, it follows from (12) that

\frac{dy}{dx} = \frac{−μF (y) − x}{y} < \frac{−μf(y) − x}{x}, ∀ x ∈ R^{1}, ∀ y ≠ 0, ∀ t ∈ R^{1}.

and it follows from the Chaplygin-Kamke comparison principle ([14–17]) that solutions G_{j}(x) corresponding to the trajectories of system

\frac{dy}{dt} = −μf(y) − x, \frac{dx}{dt} = y

and the solution (y(t), x(t)) of system (14) have such property in the point t = t_{0}, x_{0} = x(t_{0}), y_{0} = y(t_{0}) = G_{j}(x_{0}) ≠ 0:

\frac{dy}{dx} < \frac{dG_{j}}{dx} .

This means that all curves y = G_{j}(x) everywhere except for points y = 0, x ∈ R^{1} are without contact with respect to the vector field of system (11) and trajectories of system (11) go through these curves from outside to inside (see Fig.3). Thus the statement of theorem 2 follows.

Let Ψ_{0} = ∩Ψ(ν_{i}), where the intersection is taken for all parameters satisfying to condition (12).

Corollary 2.2 The global attractor of system (11) belongs to Ψ_{0}.

Algorithm for constructing Ψ_{0}

The solution of equation (14) is a solution of the corresponding discontinuous system. Each equation of the system is represented as follows:

x′′_{i}(t) + x_{i} = −μν_{i}, i = 1, …, M.

We describe the case of positive y only. If y < 0 it is easy to derive them similarly. The solution of each equation has the following representation:

x_{i}(t_{i}) = C_{i}^{1} cos(t_{i}) + C_{i}^{2} sin(t_{i}) − μν_{i}
y_{i}(t_{i}) = −C_{i}^{1} sin(t_{i}) + C_{i}^{2} cos(t_{i}),  (y_{i} = dx_{i}/dt), i = 1, …, M (19)

These are circles with centers in points (−μν_{i}, 0) and radii which are calculated from the formula

(x_{i} + μν_{i})^{2} + y_{i}^{2} = R_{i}^{2} = (C_{i}^{1})^{2} + (C_{i}^{2})^{2}.

We take the initial data B_{0} for constructing the solution of the first equation of the system (i = 1, y > 0) from formulas (9) and (10) as the maximal value of x.

The equations for points of intersection of solutions number i and i + 1 are

\begin{matrix}x_{i}(t_{i}) = & b_{i} = & x_{i+1}(t_{i}), \\ y_{i}(t_{i}) = & γ_{i} = & y_{i+1}(t_{i}). \end{matrix}  (20)

Here b_{i} = −μν_{i}−\sqrt{R_{i}^{2} − γ_{i}^{2}} are the x-coordinates of points of solutions’ intersections. Their y-coordinates obviously equal to γ_{i} (see Fig. 5). Taking into account representations of solutions (19) and these conditions we calculate coefficients C_{i+1}^{1,2} of the following solutions.

See Fig. 5. Here arcs D_{i−1}, D_{i}, D_{i+1} denote parts of solutions’ trajectories of equations with numbers i − 1, i and i + 1 respectively. The solution of equation (14) is the union of all arcs: ⋃_{i} D_{i}.

Obviously, it is possible to calculate values b_{i} only for R_{i} ≥ γ_{i}. If γ_{i} > 1 the radii R_{i} decrease with the increasing of index i and at the same time b_{i} increase. Hence there exist an index i such that the inequality will not be valid (see Fig. 6). The value of the i equals to mentioned earlier N exactly.

The construction solutions continues but we move along y-axis in inverse direction. Let us add on sequence γ such that γ_{j} for j > N will be decreasing sequence:

γ_{j} = γ_{2N−j}.

Values of b will be found from the expressions

b_{N} = −μν_{N} + \sqrt{R_{N}^{2} − γ_{N}^{2}},
b_{j} = −μν_{j} + \sqrt{R_{j}^{2} − γ_{j}^{2}} for j = N + 1, …, 2N − 1,

where ν_{j} are defined as

ν_{j} = ν_{2N−j}.

Thereby we have constructed the first part of the trajectory of the solution of equation (14) for y > 0. Let B_{1} the abscissa of the point of intersection of the trajectory and x-axis. The value x = B_{1} and y = 0 will be used again as the initial data for the construction of the trajectory for y < 0. The initial data for the following coil of the trajectory we find in the same way: it is the value of x = B_{2} - the abscissa of the point of intersection of the trajectory for y ≤ 0 and x-axis and y = 0 (see proof of Theorem 2 and Fig.2).

If |B_{0}| − |B_{2}| < 10^{−6} we denote B_{2} as B_{0} and continue construction of the trajectory of the solution of equation (14) using the described algorithm. Then the same constructions is made for the parameter μ from the interval (0; 30) with the step 0.01.

The value of ρ(μ) of the Van der Pol cycle is estimated by the following value: max |y(μ)| of the last coins of the constructed trajectories of the solutions.

For more precise estimation the same constructions is made with the decreased step S of the partitioning γ_{i}. The decreasing of the step S continues while there is an improvement of the estimation.

Let us consider an example:

According to figures 7, 8 you can see trajectories of solution of comparison systems with parameter μ_{1} = 0.1 and μ_{2} = 3.5 with the same accuracy. It is easy to understand the differences between them: for μ_{1} we need 11 winds and for μ_{2} - only 4.

An estimation of the Van der Pol cycle from inside

The following comparison system will be constructed to estimate the Van der Pol cycle from inside. Let instead of inequalities (12) for the function F (y) and for some real numbers ν_{i} the following inequalities hold

\frac{F (y)}{y} < \frac{ν_{i} sign y}{y} ,  ∀ t ∈ R^{1}, ∀ y ∈ Γ_{i},  (21)

where Γ_{i} are the same intervals.

Instead of f(y) from equation (14) we take the following function (for ∀ ν_{i} < 0, ∀ y):

f(y) = \{ \begin{matrix}−ν_{N} & for & y ∈ [−γ_{N}, −γ_{N−1}), \\ & ⋯ \\ −ν_{1} & for & y ∈ [−γ_{1}, 0), \\ ν_{1} & for & y ∈ (0, γ_{1}], \\ & ⋯ \\ ν_{N} & for & y ∈ (γ_{N−1}, γ_{N}], \end{matrix}

where γ_{i} is the same sequence (but it is not necessary that the sequence contains 1) and ν_{i} will be defined in a different way (see Fig. 9):

ν_{i} = \{ \begin{matrix}F (γ_{i}), & for & γ_{i} ≥ 1 \\ F (γ_{i+1}), & for & γ_{i} < 1, \end{matrix}  (22)

and the initial data for each value of parameter μ is the following point (B_{0} + δ,0), where B_{0} is the last initial data taken from the constructing of the trajectory of the solution for above estimation for the same value of parameter. And δ is such that the new constructing trajectory wind on the cycle from inside(see Fig.10).

Further the same simulation was made.

Look at Fig.11. Here are presented a graphics of (ρ_{1,2}(μ) − 2)) ∗ 10^{3} where ρ_{1,2} are estimations of the Van der Pol equation’s cycle from above and from below respectively. Parameter μ is taken from the interval (0, 36) to show the uniform estimation of the cycle.

There was a trembling of the graphics for μ < 0.1. It was the result from the lack of discharge of computer’s processor. A number in computer is presented as mantissa and degree. The discharge of mantissa is 16 and usually there is an approximation of numbers during a calculation. A library for more precise calculation was used. The result of modeling after the library linking is shown at fig. 11 exactly. Here is no the trembling but the strip between ρ_{1} and ρ_{2} is invisible except of intervals μ < 0.1 and μ > 35 where ρ_{2} − 2 becomes less than zero. It shows that the estimation is good.

So the following result is obtained: for any μ the maximum absolute value max |y(μ)| of the limiting cycle by y coordinate does not exceed 2.015. This is an improvement of Odani’s result.

3 An effect of French Ducks

Let us consider another example of comparisons systems method application. The following equation is the Van der Pol equation with a constant term:

\frac{dy}{dt} = −μF (y) − x,  \frac{dx}{dt} = y + a.  (23)

Here F (y) = \frac{y^{3}}{3} − y and μ, a are positive numbers.

Let us assume that for the function F (y) and for some real numbers ν_{i} the following inequalities are fulfilled:

\frac{F (y)}{y + a} > \frac{ν_{i} sign (y + a)}{y + a},  ∀ t ∈ R^{1}, ∀ y_{1} ∈ Γ_{i},  (24)

where Γ_{i} are intervals of the form

Γ_{i} = \{ \begin{matrix}(γ_{i}, γ_{i+1}] & for & y > −a, \\ [γ_{−i−1}, γ_{−i}) & for & y < −a. \end{matrix}

A sequence γ_{j} (j = −N, …, N, N is some natural number) can be constructed in the following way: The sequence contains numbers −1, 1, a and defines a partition of the interval [γ_{−N}, γ_{N}] with some step S ≤ 1 and

γ_{0} = −a, γ_{−i} = −γ_{i} + 2 ∗ a, i = 1, …, N.

For y ∈ (−1, 1) the function F decreases and for y ∉ (−1, 1) it increases. So if we define ν_{j}, for j = −N, …, N and j ≠ 0, as follows

for a < 1:

(for y > −a):

ν_{i} = \{ \begin{matrix}F (γ_{i−1}), & for & −a < y ≤ 1, \\ F (γ_{i}), & for & y > 1, \end{matrix}

(for y < −a):

ν_{−i} = \{ \begin{matrix}F (γ_{−i}), & for & −1 ≤ y < −a, \\ F (γ_{−i+1}), & for & y < −1, \end{matrix}  (25)

for a > 1:

(for y > −a):

ν_{i} = \{ \begin{matrix}F (γ_{i−1}), & for & −a < y ≤ −1, \\ F (γ_{i}), & for & −1 < y < 1, \\ F (γ_{i+1}), & for & y > 1, \end{matrix}

(for y < −a):

ν_{−i} = F (γ_{−i+1}).  (26)

i = 1, …, N then the function F will be bounded in the best way, i.e. inequalities (24) are valid. For more clearness see fig. 12,13 (for a < 1 and for a > 1), where the graph of the function F (y) is represented. Vertical lines correspond to the value of γ_{i}, and horizontal ones correspond to appropriate ν_{i}.

Let us consider the equation:

F \frac{dF}{dx} = −μf(F) − x,  (27)

where

f(y) = \{ \begin{matrix}ν_{−N} & for & y ∈ [γ_{−N}, γ_{−N+1}), \\ & ⋯ \\ ν_{−1} & for & y ∈ [γ_{−1}, γ_{0}), \\ ν_{1} & for & y ∈ (γ_{0}, γ_{1}], \\ & ⋯ \\ ν_{N} & for & y ∈ (γ_{N−1}, γ_{N}]. \end{matrix}

for any ν_{i} and y.

We deal with solutions F_{1}(x) and F_{2}(x) of equation (27) which are situated in half-planes {F ≥ −a} and {F ≤ −a} respectively with initial data F_{1}(B_{0}) = −a. For parameter μ ≥ 1 B_{0} is taken again from (9) as a maximum of |x| and for μ < 1 it is defined from (10). Let F_{2}(B_{1}) = −a, where B_{1} is such that F_{1}(x) > −a for x ∈ (B_{0}; B_{1}). The value B_{2} is taken such that F_{2}(B_{2}) = −a and F_{2}(x) < −a for x ∈ (B_{1}; B_{2}).

Let us define the set Ψ_{a}(ν_{i}) as

Ψ_{a}(ν_{i}) = \{ \begin{matrix}x ∈ [B_{2}, B_{1}], & F_{2}(x) ≤ y ≤ F_{1}(x), \\ x ∈ [B_{0}, B_{2}], & y ≥ 0. \end{matrix} \}

Theorem 3.1 The global attractor Φ(t, t_{0}) of system (23) belongs to Ψ_{a}(ν_{i}).

Let Ψ_{a0} = ∩Ψ_{a}(ν_{i}), where the intersection is taken for all parameters satisfying to condition (24).

Corollary 3.1 The global attractor of system (24) belongs to Ψ_{a0}.

Algorithm for constructing Ψ_{a0}

The algorithm is similar to the algorithm for constructing Ψ_{0}. The solution of equation (27) is a solution of the corresponding discontinuous system. Each equation of the system is represented as follows:

x′′_{i}(t) + x_{i} = −μν_{i}, i = 1, …, M, M > 0.

with additional condition

x′ = y + a

The solution of each equation has the following representation:

x_{i}(t_{i}) = C_{i}^{1} cos(t_{i}) + C_{i}^{2} sin(t_{i}) − μν_{i}
y_{i}(t_{i}) = −C_{i}^{1} sin(t_{i}) + C_{i}^{2} cos(t_{i}) − a, (y_{i} = dx_{i}/dt − a), i = 1, …, M (28)

These are circles with centers in the points (−μν_{i}, −a) and radii which are calculated from the formula

(x_{i} + μν_{i})^{2} + (y_{i} + a)^{2} = R_{i}^{2} = (C_{i}^{1})^{2} + (C_{i}^{2})^{2}.

We take the initial data B_{0} for constructing the solution of the first equation of the system (i = 1, y > −a) from formulas (9) and (10) like in the algorithm for constructing Ψ_{0}.

The equations for points of intersection of solutions number i and i + 1 are

\begin{matrix}x_{i}(t_{i}) = & b_{i} = & x_{i+1}(t_{i}), \\ y_{i}(t_{i}) = & γ_{i} = & y_{i+1}(t_{i}). \end{matrix}  (29)

Here b_{i} = −μν_{i} − \sqrt{R_{i}^{2} − γ_{i}^{2}} are x-coordinates of points of solutions’ intersections. Their y-coordinates obviously equal to γ_{i} (see Fig. 14). Taking in account representations of solutions (28) and these conditions we calculate coefficients C_{i+1}^{1,2} of the following solutions.

See Fig. 14. Here arcs D_{i−1}, D_{i}, D_{i+1} denote parts of solutions’ trajectories of equations with numbers i − 1, i and i + 1 respectively. The solution of equation (27) is the union of all arcs: ⋃_{i} D_{i}.

When R_{i} became less then γ_{i} ( index i is equal to mentioned earlier N) the construction of solution continues in inverse direction along the y-axis. Here we again add on sequence γ such that γ_{j} for j > N will be a decreasing sequence:

γ_{j} = γ_{2N−j}.

Values of b will be found from the expressions

b_{N} = −μν_{N} + \sqrt{R_{N}^{2} − γ_{N}^{2}},
b_{j} = −μν_{j} + \sqrt{R_{j}^{2} − γ_{j}^{2}} for j = N + 1, …, 2N − 1,

where ν_{j} are defined as

ν_{j} = ν_{2N−j}.

Thereby we have constructed the first part of the trajectory of the solution of equation (27) for y + a > 0. Let B_{1} the abscissa of the point of intersection of the trajectory and line y = −a. The value x = B_{1} and y = −a will be used again as the initial data for the construction of the trajectory for y < −a. The initial data for the following coil of the trajectory we find in the same way: the value x = B_{2} is the abscissa of the point of intersection of the trajectory for y ≤ −a and line y = −a.

If |B_{0}| − |B_{2}| > 10^{−6} we denote B_{2} as B_{0} and continue the construction of the trajectory of the solution of equation (27) using the described algorithm.

At last the Van der Pol cycle is estimated by the last coins of the constructed trajectory.

For more precise estimation the same constructions is made with the decreased step S of the partitioning γ_{i}. The decreasing of the step S continues while there is an improvement of the estimation.

You can see at figures 15-19 how the shape of the cycle changes with changing of the parameter a, μ = 0.01. a belongs to interval (a_{0}, a_{1}), a_{0} = 0.998740512, a_{1} = 0.998740513. There is a bifurcation depending on parameter a. For a < a_{0} the shape of the cycle changes very slowly. For a_{0} < a < a_{1} at first we see as down arc of the cycle goes up and so called ”duck with a head” appears (see Fig.17-18). Then the cycle decreases very fast. At last the cycle disappears and a point (−a, f(a)) becomes a stable stationary point. These figures corresponds to theoretical results (see [39], [4]). In [4] pictures only for a = a_{0} and for a = a_{1} are presented. Here we obtained more precise ones.

12

References

[1] P. Alsholm : Existence of limit cycles for generalized Liénard equation. J. Math. Anal. Appl., 171 (1992), no. 1, pp. 242–255.

[2] J.-L. Callot, F. Diener, M. Diener : Le probleme de la ”chasse au canard”. C.R.Acad. Sci. Peris, ser. 1,1978, 286: 22, p. 1059-1061.

[3] P. Cartier : Perturbations singulieres des équations différentielles ordinaires et analyse non-standard. (French)[Singular perturbations of ordinary differential equations and nonstandard analysis] Bourbaki Seminar, Vol. 1981/1982, pp. 21–44, Astérisque, 92–93, Soc. Math. France, Paris, 1982.

[4] P. Cartier : Singular perturbations of ordinary differential equations and non-standard analysis. (Russian) Translated from the French by A. K. Zvonkin, Uspekhi Mat. Nauk 39 (1984), no. 2(236), pp. 57-76.

[5] M. L. Cartwright : Van der Pol’s equation for relaxation oscillation. Contribution to the Theory of Non-linear Oscillations, vol. II, pp. 3–8. Princeton Univ. Press, Princeton, 1952.

[6] S. A. Chaplygin : A new method of an approximate integrating of differential equations. (Russian) ”Gostechizdat”, Moscow, 1950.

[7] L. Cezare : Asymptotic behavior and stability of solutions of ordinary differential equations. (Russian) Izdat. ”Mir”, Moscow, 1964, 477 pp.

[8] Ph. le Colbeiller : J. Inst. Elec. Eng., 79(1936), pp. 361–278.

[9] M. Diener : Nessie et les canards. Strasbourg: Publications IRMA, 1979.

[10] M. Diener : Etude générique descanars. Thése, Strasbourg, 1981.

[11] W. Eckhaus : Relaxation oscillations including a standard chase on French ducks. Asymptotic analysis, II, pp. 449–494, Lecture Notes in Math., 985, Springer, Berlin-New York, 1983.

[12] A. F. Filippov : Differential equations with discontinuous righthand sides. Translated from the Russian. Mathematics and its Applications (Soviet Series), 18. Kluwer Academic Publishers Group, Dordrecht, 1988. x+304 pp.

[13] D. A. Flanders and J. J. Stoker : The limit case of relaxation oscillation. Studies in Nonlinear Vibration Theory, pp. 50–64. Inst. for Math. and Mech., New York Univ., 1946.

[14] A. Kh. Gelig, G. A. Leonov, V. A. Yakubovich : The stability of nonlinear systems with a nonunique equilibrium state. (Russian) Izdat. ”Nauka”, Moscow, 1978, 400 pp.

[15] J. Grasman : Asymptotic Methods for Relaxation Oscillations and Applications. Applied Mathematical Sciences, Vol 63

[16] J. Haag : Examples concrets d’étude asymptotique d’oscillations de relaxation. (French) Ann. Ec. Norm.Sup., 61(1944), 73-117.

[17] O.Hájek : Cycle periods of plane dynamic systems. J. Math. Anal. Appl., vol. 9 (1964), pp. 75–84.

[18] J. K. Hale : Averaging methods for differential equations with retarded arguments and a small parameter. J. of Differential Equations, 2 (1966), no. 1, pp. 57–73.

[19] W. T. van Horssen : A perturbation method based on integrating factors. (English summary), SIAM J. Appl. Math., 58 (1999), no. 4, pp. 1427–1443 (electronic).

[20] E. Kamke : Zur Theorie der Systeme gewöhnlicher Differentialgleichungen II. Acta Math. Bd.58 (1932), pp. 57–85.

[21] M. V. Keldysh : About dumping with nonlinear characteristics. Collected works: Mech., (Russian) Izdat. ”Nauka”, Moscow, 1964, pp. 475–490.

[22] P. J. Lau : Bounds for the Length and periods of closed orbits of two-dimensional autonomous systems of differential equations. Differential Equation, 3 (1967), no. 3, pp. 330–342.

[23] J. P. LaSalle : Relaxation oscillations. Quart. Appl. Math., 7 (1949), pp. 1–20.

[24] S. Lefschetz : Differential equations: Geometric theory. 2nd Ed., Interscience,1963; reprint, Dover, NewYork, 1977.

[25] G. A. Leonov : Oscillations in systems with nonlinear damping. (Russian) Prikl. Mat. Mekh. 57 (1993), no. 5, pp. 183–184; translation in J. Appl. Math. Mech. 57 (1993), no. 5, pp. 945–946 (1994)

[26] G. A. Leonov : Localization of the attractors of the nonautonomous Liénard equation by the method of discontinuous comparison systems. (Russian) Prikl. Mat. Mekh. 60 (1996), no. 2, pp. 332–336; translation in J. Appl. Math. Mech. 60 (1996), no. 2, pp. 329–332 [27] G. A. Leonov : A lower bound on the number of cycles in two-dimensional dynamical systems. (Russian.

English, Russian summary) Vestnik S.-Peterburg. Univ. Mat. Mekh. Astronom., 1994, vyp. 1, pp. 42–46.

[28] G. A. Leonov, V. Reitmann, V. B. Smirnova : Non-local methods for pendulum-like feedback systems. (English.)Teubner-Texte zur Mathematik, Bd. 132 (1992).

[29] Y. S. Lim, L. F. Kazda : A study of second order nonlinear systems. J. Math. Anal. Appl., 8 (1964), pp. 423– 444.

[30] A. D. MacGillivray; B. Liu; N. D. Kazarinoff : Asymptotic analysis of the peeling-off point of a French duck. Methods Appl. Anal. 1 (1994), no. 4, pp. 488–509.

[31] R. E. O’Malley : Singular Perturbation Methods for Ordinary Differential Equations. Applied Mathematical Sciences, Vol 89.

[32] K. Odani : Existence of exactly N periodic solutions for Liénard systems. Funkcial. Ekvac., 39 (1996), no. 2, pp. 217–234.

[33] K. Odani : On the limit cycle of the van der Pol equation. ”Equadiff99 CD-ROM: Papers”, Z. Doslá, J. Kuben,J. Vosmanský, eds., Masaryk Univ., Czech, 1998, pp. 229–235.

[34] K. Odani : On the limit cycle of the van der Pol equation. preprint.

[35] V. A. Pliss : Non-local problems in the theory of oscillations. (Russian) Izdat. ”Nauka”, Moscow/Leningrad, 1964, 367 pp. Translatinon by Scripta Technica, Inc. Translation edited by Harry Herman Academic Press, New York-London 1966 xii+306 pp.

[36] G. E. H. Reuter : Boundedness theorems for non-linear differential equations of the second order. II. J. London Math. Soc. 27(1952), pp. 48–58.

[37] A. Troesch, E. Urlacher : Analyse non standard et l’équation de Van der Pol. Strasbourg: Publications IRMA, 1977.

[38] A. Troesch, E. Urlacher : Perturbations singuliéres et analyse nonstandard: C^{k}-convergense et crépitement des solutions. Strasbourg: Publications IRMA, 1977.

[39] A.K. Zvonkin and M.A. Shubin : Non-standard analysis and singular perturbations of ordinary differential equations. (Russian) Uspekhi Mat. Nauk 39 (1984), no. 2(236), pp. 77–127.

[40] B. van der Pol : Forced oscillations in a systems with nonlinear resistance. Philos. Mag. 2 (1926), pp. 978–992.