On the separation of motions in systems with a large fast excitation of general form

Eur. J. Mech. A/Solids 18 (1999) 527–538  Elsevier, Paris

On the separation of motions in systems with a large fast excitation of general form A. Fidlin LuK GS GmbH, Bühl, Germany (Received 27 November 1997; revised and accepted 31 July 1998) Abstract – In this study dynamic systems are considered, in which motion can be described through a system of second-order ordinary differential equations with the right sides depending both on the slow time t and on the fast time τ = ωt (ω  1 is a big dimensionsless parameter). It is assumed that the right sides are large (they have the magnitude order ω) and depend both on generalised coordinates and on the generalised velocities of the system. A motion separating procedure is developed for the systems described in twi ways. The procedure enables separate systems to derive for fast (oscillating) and slow components of the solution. Each of these separated systems is simpler than the original one. The equivalence of both procedures is shown. The first of them is based on the multiple scales method, the second one generalises the averaging method of Krylov–Bogoliubov–Mitropolskii. Motion of a linear oscillator excited through the large fast oscillations of the damping coefficient is analysed as an example of the established method usage. It is shown, that the excitation significantly changes the effective mass and in consequence the natural frequency of the original system. The analytic results are compared with numerical ones.  Elsevier, Paris oscillations theory / averaging / multiple scales

1. Introduction The separation of motions is one of the main ideas for asymptotic analysis of oscillating systems with small or big parameter. It is connected with the fact that solutions of many types of dynamic systems can be represented as a superposition of fast oscillations and slow evolution of the solution. The averaging method (Krylov and Bogoliubov, 1957; Bogoliubov and Mitropolskii, 1974) and the method of multiple scales (Nayfeh, 1973; Nayfeh and Mook, 1979), which differ significantly in form, are substantially very close. Together with the method of the direct separation of motions, which finds its origins in the works of Kapitsa (1951, 1954) and is more generally formulated in works of Blekhman (1988, 1994), these methods are the most effective practical realisations of this idea. Each of these methods has its own advantages in solving the special problems in the theory of oscillations. In this paper we study the response of non-linear systems to a strong excitation in general form, whose carrier frequency ω is much higher than natural frequencies of the system. Such systems are one of the main objects of interest for the method of the direct separatoin of motions, which however is proved only for special types of fast excitation. Method of multiple scales was also successfully used to analyse the response of some particular systems to these types of strong excitation (Nayfeh and Mook, 1979; Nayfeh and Nayfeh, 1995). We are going first to establish a formal procedure for the separation of motions in non-linear systems with strong fast excitation in its most general form based on the multiple scales method technique. We are then going to generalise the averaging method technique, which enables us to demonstrate the internal distinctions of the considered problem and to prove developed technique. A simple example of the use of this procedure can be found at the end of this paper.

528

A. Fidlin

2. Problem statement. A formal solution through the multiple scales technique Consider a system x •• = F (x, x • , t) + ω8(x, x • , t, ωt).

(1)

Here x is an n-dimensional vector of the generalised coordinates and x • is a vector of the generalised velocities. We take F and 8 to be n-dimensional vectors of forces. The components of F have to be bounded functions, satisfying Lipschitz conditions on the first and second arguments, and the components of 8 have to satisfy these conditions together with their first partial derivatives with respect to x and x • . We take also 8 to be 2π -periodic vector-function with respect to τ = ωt, and ω to be a big parameter. This system differs from that considered in the work of Blekhman (1994), because 8 does not depend only on the generalised coordinates x, but depends also on the generalised velocities x • , i.e., we are considering the large fast excitation in its most general form. In order to apply the multiple scales technique to (1), we have to convert to two independent variables t and τ , i.e., from the system of ordinary differential equations (1) to the following system with partial derivatives 







2 ∂ 2ϕ ∂ 2ϕ ∂ϕ ∂ϕ ∂ϕ ∂ϕ 2∂ ϕ + 2ω = F ϕ, + ω + ω , t + ω8 ϕ, + ω , t, τ . 2 2 ∂t ∂t∂τ ∂τ ∂t ∂τ ∂t ∂τ

(2)

The relationship between (1) and (2) is defined through the condition that, if ϕ(t, τ ) is a solution of (2), then x = ϕ(t, ωt) is a solution of (1). In other words, system (2) is more general then (1), and any solution of (2) taken at the straight line τ = ωt satisfies (1). We require ϕ(t, τ ) to be a 2π -periodic function of τ and try to find ϕ as a formal asymptotic expansion in terms of ε = 1/ω ϕ(t, τ ) = ψ0 (t, τ ) + εψ1 (t, τ ) + ε 2 ψ2 (t, τ ) + · · · . Substituting this expression in (2) and balancing the terms with equal orders of ε, we obtain ε −2 :

ε

−1

∂ 2 ψ0 = 0, ∂τ 2

(3)





∂ 2 ψ1 ∂ 2 ψ0 ∂ψ1 ∂ψ0 ∂ψ0 : + 2 = 8 ψ0 , + +ω , t, τ , ∂τ 2 ∂t∂τ ∂τ ∂t ∂τ 

(4) 

∂ 2 ψ2 ∂ 2 ψ1 ∂ 2 ψ0 ∂8 ∂8 ∂ψ1 ∂ψ2 ε : + 2 =F + . + ψ1 + • + ∂τ 2 ∂t∂τ ∂t 2 ∂x ∂x ∂t ∂τ 0

(5)

The general solution of (3) is ψ0 (t, τ ) = X(t) + A(t)τ.

(6)

According to the periodicity of ψ0 A(t) ≡ 0. Hence ψ0 = X(t) depends only on the slow time t. The objective of the following analysis is to discover the differential equations for this still unknown function, which do not contain the fast time τ .

Separation of motions in systems with a large fast excitation of general form The Eq. (4) takes the form



529



∂ 2 ψ1 ∂ψ1 = 8 X, X • + , t, τ . 2 ∂τ ∂τ

(7)

The main assumption of this investigation is that we take to know a general 2 π -periodic with respect to τ solution of the system of n differential equations of the first order ∂u = 8(X, u, t, τ ) ∂τ

(8)

depending on a vector of arbitrary constants C (here X and t are taken to be fixed) u = U (X, t, τ, C). We determine these constants to annihilate the average of

∂ψ1 ∂τ

(ψ1 has to be periodic)



U (X, t, τ, C) = X • .

Here and further, h i means the averaging of the expression in brackets, calculated during the time interval of periodicity of all the considered functions. Hence u = U (X, X •, t, τ ) and ψ1 =

Z 0

τ

(U − X • ) dτ + X1 (t).

The unknown function X1 (t) is a small (magnitude order ε) slow correction to the main slow part of the solution X(t). Denoting 9(t, τ ) =

Z

0

τ

(U − X • ) dτ

(9)

(9 is a periodic function of τ and depends on t both direct and indirect through X and X • ), we obtain ψ1 (t, τ ) = 9(t, τ ) + X1 (t).

(10)

The Eq. (5) takes now the following form ∂ 2 ψ2 ∂8 ∂ψ2 ∂8 ∂ψ1 ∂8 ∂ 2 ψ1 = + + F + ψ − 2 − X •• . 1 ∂τ 2 ∂x • ∂τ ∂x ∂x • ∂t ∂t∂τ

(11)

1 All the functions here should be calculated at the point (X, X • + ∂ψ , t, τ ). ∂τ System (10) is a system of inhomogeneous linear differential equations in respect to τ . All the coefficients 2 are periodic and the unknown function is ∂ψ . ∂τ As is known from the general theory of linear systems with periodic coefficients, for the existence of periodic solutions of (11) the projections of the inhomogeneous parts of the equations on the solutions of the system

530

A. Fidlin

conjugated to the homogeneous one must vanish. [A reference to the classical work of Poincaré (1957) seems here to be in order and to underscore the close relationship of Poincaré’s method to the method of multiple scales, which is in this case a procedure to find periodic in respect to τ solutions of (2) with a not isolated undisturbed solution (6).] In other words, we introduce a homogeneous system ∂w ∂8 = •w ∂τ ∂x and the conjugated one 

∂w∗ ∂8 =− ∂τ ∂x •

T

w∗ .

(12)

Let us assume we have got n independent periodic solutions of (12), the matrix of which we will design as W∗ . Then for the existence of a periodic solution of (11) is necessary that 



W∗T

∂8 ∂8 ∂ψ1 ∂ 2 ψ1 F+ ψ1 + • −2 − X •• ∂x ∂x ∂t ∂t∂τ

or

X •• = W∗T

−1





W∗T F +



=0

∂8 ∂8 ∂ψ1 ∂ 2 ψ1 ψ1 + • −2 ∂x ∂x ∂t ∂t∂τ



.

(13)

1 The Eq. (13) seems to contain through ψ1 and ∂ψ side by side with the already known functions W∗ and 8 ∂t and the desired variable X also the unknown function X1 (t) [see (10)]. Let us show that all the terms containing X1 (t) vanish. Really





W∗T





∂8 • ∂8 X = W∗T • X1• = ∂x • 1 ∂x



∂8 ∂x •

T

T

W∗

X1• = −



∂W∗ ∂τ

T

X1• .

But we have assumed, that W∗ (τ ) is periodic. Hence 

and therefore

∂W∗ ∂τ





=0 

W∗T

∂8 • X = 0. ∂x • 1

In order to calculate hW∗T ∂8 X i let us consider (8) and derive it partially in respect to X ∂x 1 ∂ ∂τ



∂u ∂X



=

∂8 ∂8 ∂u + • . ∂x ∂x ∂X

Multiplying the last equation with W∗T one obtains 

W∗T

∂8 ∂u ∂ W∗T = ∂x ∂τ ∂X



Separation of motions in systems with a large fast excitation of general form and because of the periodicity of W∗ , u and 

∂u ∂X



∂8 X1 = ∂x

W∗T

531





∂ ∂u W∗T ∂τ ∂X



X1 = 0.

(14)

So the function X1 does not occur in (13) and ψ1 can be replaced with the known function 9(t, τ ) (9). Let us notice in addition, that 



W∗T

∂8 ∂9 ∂ 29 − ∂x • ∂t ∂t∂τ







∂ ∂9 =− W∗T ∂τ ∂t



=0

because of the periodicity of W∗ , 9 and ∂9 . ∂t So we get the desired equation for X(t) ••

X =



−1 W∗T



W∗T



∂8 ∂ 29 F+ 9− ∂x ∂t∂τ



.

(15)

With this procedure we have formally separated the fast oscillating part of the solution (1), which is determined through (8), (9), called according to Blekhman (1994) the equations of fast notions, from the slow ‘evolutionary’ component of the solution. This was determined through the equations of slow motions (15). So, to find X(t) we have first of all to solve system (8), the order of which is two times lower than the order of the original system, then one has to solve the linear system (12) and to determine the matrix W∗ . At last it is possible to get the system (15), which differs from (1) because it does not contain the fast time τ = ωt. Equation (15) is not solved in respect to X •• . In fact [see (9)] ∂ 29 ∂U • ∂U ∂U •• X − X •• . = X + + ∂t∂τ ∂X ∂t ∂X • Substituting this relationship in (15), we obtain M(X, X • , t)X •• = F (X, X • , t) + V (X, X • , t) with

−1 W∗T



(16)



∂U M(X, X , t) = , ∂X •   

−1

−1 T   ∂8 ∂U • ∂U V (X, X • , t) = W∗T W∗T + W∗T W∗ F (X, U, t) − F (X, X • , t) . 9− X − ∂x ∂X ∂t •

W∗T

Function V (X, X • , t) can be called, according to Blekhman (1994), the ‘vibration force’. The matrix coefficient M(X, X • , t) can be interpreted as a matrix of efficient mass or inertia of the averaged system in respect to its slow motions. We should notice that this matrix depends on the solutions of the equations of fast motions, i.e., on the amplitude of fast excitation. This described procedure together with its validation is enough for the asymptotic analysis of (1), but we are going to investigate this system once more with the help of the averaging method in order to expose the internal essence of the transformations carried out. At the same time we are going to validate our procedure.

532

A. Fidlin

3. The generalised averaging method. Validation of the procedure As before, we are going to consider system (1). We also assume to have found a periodic solution for the equations of fast motions (7), satisfying together with their first partial derivatives with respect to X and X • Lipschitz conditions. By solving this equation we assume X, X • and t to be some constant parameters [Eq. (7) is a system with partial derivatives in respect to τ ]. Then we can use the following transformation x = X(t) + εψ1 (X, Y, t, τ ),

(17)

∂ψ1 (X, Y, t, τ ) x = Y (t) + . ∂τ •

Here X and Y are assumed to depend only on the slow time t but not on τ . The substitution similar to (17) was used by Bogoliubov and Mitropolskii (1974) to investigate motions of a pendulum with the oscillating suspension point, and by Blekhman and Malakhova (1974) to validate the method of the direct separation of motions for systems similar to (1), but with 8 depending only on x, t and τ , but not on X • . If one substitutes (17) into (1), one gets ∂ψ1 • ∂ψ1 • ∂ψ1 X +ε Y +ε = Y, ∂X ∂Y ∂t ∂ 2 ψ1 • ∂ 2 ψ1 • ∂ 2 ψ1 1 ∂ 2 ψ1 Y• + X + Y + + ∂X∂τ ∂Y ∂τ ∂t∂τ ε ∂τ 2     ∂ψ1 1 ∂ψ1 = F X + εψ1 , Y + , t + 8 X + εψ1 , Y + , t, τ . ∂τ ε ∂τ X• + ε

(18)

Expanding the right side of the last equation in terms of ε and balancing the terms of 1ε -order one gets the equation of fast motion 



∂ 2 ψ1 ∂ψ1 = 8 X, Y + , tτ . 2 ∂τ ∂τ

(19)

This equation is equivalent to (7) if X and Y do not depend on τ . If, as assumed, we known function ψ1 for any constant X, Y , t (accurate to an arbitrary additive function of the slow time t), we can rewrite system (18) in the following form X • = Y + ε O(1), 

Y• = E +

−1 

∂ 2 ψ1 ∂Y ∂τ

F+

∂8 ∂ 2 ψ1 ∂ 2 ψ1 ψ1 − Y− ∂X ∂X∂τ ∂t∂τ



+ ε O(1).

(20)

Here and further, O(1) is a value of the magnitude order 1 and E is a unit matrix. Converting to the fast time τ as an independent variable we get X 0 = εY + ε 2 O(1), 

∂ 2 ψ1 Y =ε E + ∂Y ∂τ 0

−1 

∂8 ∂ 2 ψ1 ∂ 2 ψ1 F+ ψ1 − Y− ∂X ∂X∂τ ∂t∂τ



+ ε 2 O(1).

(21)

It seems to be a system in standard form (Bogoliubov and Mitropolskii, 1974), to which one can apply the averaging method. But it is wrong. The averaging method does not mean simply the averaging of the right parts of a system in standard form in respect to τ . The averaging method is a perturbing change of variables

Separation of motions in systems with a large fast excitation of general form X = ξ + εu(ξ, η, t, τ ),

533 (22)

Y = η + εv(ξ, η, t, τ ). Here u, ν are 2π -periodic functions of τ , and ξ and η have to fulfill the equations ξ 0 = ε4(ξ, η, t),

(23)

η0 = εN(ξ, η, t)

the right sides of which do not depend on τ . But the substitution (22) stays in an evident contradiction to (17) and the described way to solve (19). Indeed, we have assumed that X and Y depend only on t. But now we suppose that they depend also on τ through the functions u and ν. It means that in reality not X and Y , but ξ and η are the slow variables, for which we are going to get the slow equations (23). In other words, if we are going to change variables according to (22), (23), we have to substitute the new variables not only in Eq. (21), but also in the equations for fast motions (19) or directly in the original Eq. (18). This leads to the following result: the equation of fast motions takes the form 

∂ 2 ψ1 ∂ψ1 = 8 ξ + εu, η + εν + , t, τ ∂τ 2 ∂τ



or after the expansion in terms of ε 



∂ 2 ψ1 ∂ψ1 ∂8 ∂8 = 8 ξ, η + , t, τ + ε u + ε • v. ∂τ 2 ∂τ ∂x ∂x So, if ψ1 fulfills (19) with constant X and Y , it does not fulfill it with constant ξ and η, or more exactly formulated it fulfills this equation accurately to the terms of magnitude order ε, which have to be added to the right side of (20) or (21) if one wants to get an approximate solution to the original system (1). Hence system (21) can be rewritten as following X 0 = εY + ε 2 O(1), 

∂ 2 ψ1 Y =ε E+ ∂Y ∂τ 0

−1 



∂8 ∂ 2 ψ1 ∂ 2 ψ1 ∂8 ∂8 F+ ψ1 − Y− + u + • v + ε 2 O(1). ∂X ∂X∂τ ∂t∂τ ∂x ∂x

Here X, Y , X 0 , Y 0 have to be expressed in terms of ξ , η, 4, N according to substitution (22), (23). After this, one will get the following system of differential equations with partial derivatives to find the unknown periodic in terms of τ functions u and v 4(ξ, η, t) +

∂u = η + ε O(1), ∂τ

∂ 2 ψ1 ∂ 2 ψ1 ∂v ∂8 ∂ 2 ψ1 ∂8 ∂8 N(ξ, η, t) + N+ =F + ψ1 − η− + u + • v + ε O(1). ∂η∂τ ∂τ ∂x ∂ξ ∂τ ∂t∂τ ∂x ∂x

(24)

From the first equation and the requirement of the periodicity of u(τ ) follows 4 = η,

u = A(t).

(25)

534

A. Fidlin

Here A(t) is an arbitrary function of the slow time t. The arbitrary function contained in ψ1 can be added to A(t) without any loss of generality. Then the second equation of (24) takes the form 

∂ν ∂ 29 ∂8 ∂8 ∂ 29 ∂ 29 = •v +F + 9− η+ N+ ∂τ ∂x ∂x ∂ξ ∂τ ∂η∂τ ∂t∂τ



−N +

∂8 A(t) + ε O(1). ∂x

(26)

Here as before, 9 is defined through (9). One should notice that the expression in brackets is nothing other than the derivative in respect to τ from the full partial derivative of 9(t, τ ) in respect to t (accounting that 9 depends on t both directly and indirectly through ξ and η). So one rewrite (26) ∂ν ∂8 ∂8 ∂ 29 ∂8 = •v +F + 9− −N + A(t) + ε O(1). ∂τ ∂x ∂x ∂t∂τ ∂x

(27)

System (27) is very similar to (11) and differs significantly from the conventional equations of the averaging ∂8 method through term ∂x • v. The necessary condition of the existence of the periodic solution of (27) under account of (14) has the following form 



W∗T

∂8 ∂ 29 F+ 9− −N ∂x ∂t∂τ



+ ε O(1) = 0

or together with (25) ξ

••

=



−1 W∗T





W∗T

∂8 ∂ 29 F+ 9− ∂x ∂t∂τ



+ ε O(1) = 0.

(28)

This equation is equivalent to the original system (1). Let us consider it together with the equation expressing the first approximation of the averaging method

ξ¯ •• = W∗T

−1





W∗T

∂8 ∂ 29 F+ 9− ∂x ∂t∂τ



(29)

which coincides with (15) and differs from the exact Eq. (28) through the omission of terms ε O(1). The justification of this jump from (28) to (29) with the help of Gronwalls lemma does not differ from the classical validation of the averaging method (Zhuravlev and Klimov, 1988) and ensures that kξ − ξ¯ k 6 C1 ε for the time interval τ 6 O(ε −1 ) or t 6 O(1), if all functions on the right side of (29) and W∗ confirm Lipschitz conditions and, besides that, if matrix hW∗T i−1 exists, i.e., the determinant of hW∗T i−1 is never equal to zero. Remark 1: If we have the initial conditions for the original system x|t =0 = x0 ,

x • |t =0 = v0

then the initial conditions for the averaged system (15) or (29) can be calculated as follows (we use the designations of (15)) X|t =0 = x0 ,



∂9 X |t =0 = v0 − . ∂τ t =0,τ =0 •

Separation of motions in systems with a large fast excitation of general form

535

Remark 2: If 8 does not depend on x • , then 9 also does not depend on X • and system (15) or (29) goes over into equations obtained both with the conventional averaging method and with the method of direct separation of motions (Blekhman, 1994) 



∂8 X = F+ 9 . ∂x ••

4. Example. Change in dynamic properties of a linear oscillator under the influence of large fast oscillations of ‘damping’ coefficient As an example of use of the described method we are going to analyse the following equation x •• + βx • + x = aωx • cos(ωt).

(30)

Let us assume that a and β both are magnitude order 1 and ω  1 is the big parameter. Then according to our designations F = −βx • − x,

8 = aωx • cos τ.

The equation of fast motions takes the form 



∂ 2 ψ1 ∂ψ1 = a X• + cos τ 2 ∂τ ∂τ and its periodic solution with the vanishing average is 



∂9 ea sin τ = − 1 X• , ∂τ I0 (a) 

1 9= I0 (a)

Z

τ

a sin ϑ

e 0

1 dϑ − τ + 2π I0 (a)

Here I0 (a) =

1 2π

Z

2π

Z

2π

e 0

ea sin ϑ dϑ

0

is the modified Bessels function of 0-order. The linear homogeneous system (12) has the form ∂W∗ = −aW∗ cos τ ∂τ and its periodic solution is W∗ = e−a sin τ ,

W∗ > 0.

Averaging of terms gives hW∗ i = I0 (a), hW∗ i = −XI0 (a) −

 a sin ϑ

βX • , I0 (a)

(31) •

dϑ − π X .

536

A. Fidlin 





∂ 29 1 W∗ = X •• − I0 (a) . ∂t∂τ I0 (a)

After substitution of these expressions into (15), one gets the following equation of slow motion 1 βX • X •• = −X − 2 I0 (a) I0 (a)

Figure 1. Comparison of analytic and numeric solutions for ω = 50.

Separation of motions in systems with a large fast excitation of general form

537

X •• + βX • + I02 (a)X = 0.

(32)

or

(32) is a very simple linear differential equation with constant coefficients. Its solution, together with the fast component (31) and formulas (17), gives us an approximate solution of (30). Let us notice some features of this solution.

Figure 2. Comparison of analytic and numeric solutions for ω = 10.

538

A. Fidlin

1) Big fast oscillations of the damping coefficient change the frequency of slow free oscillations of the system. This frequency can be controlled by the excitation amplitude. This result is similar to that reported by Nayfeh and Nayfeh (1995) for another particular system. Equation (16) shows that it is not an accident, but a characteristic property of systems with strong fast excitation in general form. In our case we can interpret the change of the ‘natural’ frequency of the averaged system as a consequence of the change of the effective mass of the system in respect to slow motions. 2) If there is no oscillation (a = 0), the averaged Eq. (32) goes over into the original one (30) because I0 (0) = 1. 3) For this problem and all problems we analyse here it is typical that the solution is a superposition of slow component and fast oscillations. These are small respectively to the generalised coordinate x, but their derivatives are not small respectively to the generalised velocities x • . A comparison of the analytical solution with a numerical one illustrates the last statement. It can be found in figures 1 and 2 for two different values of big parameter (ω = 50, figure 1; ω = 10, figure 2). In both cases the calculations were done for the following parameters: β = 0, a = 1, ⇒ I0 (1) ≈ 1,266. In the first case there is no visible difference between two predictions. For smaller values of ω the trajectories diverge already at early times. In figure 1 we can see that the oscillations period is decreased from 2π for the unperturbed system (a = 0) up to 2π/I0 (1) ≈ 4,963. One can also see that the solution for x contains almost only the slow component and the solution for x • is a superposition of a slow component and of big fast oscillations with slow variable amplitude. Acknowledgements The author wishes to thank I.I. Blekhman (Institute of Mechanical Engineering Problems of the Russian Academy of Sciences, St. Petersburg) and R. Fischer (LuK GS GmbH, Buehl, Germany) for their support. References Blekhman I.I., 1988. Synchronisation in Science and Technology. ASME Press, New York. Blekhman I.I., 1994. Vibration Mechanics. Nauka, Moscow (in Russian). Bogoliubov N.N., Mitropolskii Y.A., 1974. Asymptotic Methods in Non-linear Oscillation Theory. Nauka, Moscow (in Russian). Kapitsa P.L., 1951. Dynamic stability of a pendulum with oscillating suspension point. J. Exp. Theor. Phys. 21, 5. Kapitsa P.L., 1954. Pendulum with vibrating suspension. Achieve. Phys. Sci. 44, 1. Krylov N.M., Bogoliubov N.N., 1957. Introduction to Non-linear Mechanics. Princeton Univ. Press, Princeton, USA. Nayfeh A.H., 1973. Perturbation Methods. Wiley Interscience, New York. Nayfeh A.H., Mook D.T., 1979. Non-linear Oscillation. Wiley, New York. Nayfeh S.A., Nayfeh A.H., 1995. The response of nonlinear systems to modulated high frequency input. Nonlinear Dyn. 7, 310–315. Poincaré H., 1957. Les méthodes nouvelles de mèchanique céleste (3 vol.). Dover, New York. Zhuravlev V.F., Klimov D.M., 1988. Applied Methods in Oscillation Theory. Nauka, Moscow (in Russian).