Dynamics of a Pendulum of Variable Length and Similar Problems

In this chapter we study three mechanical problems: dynamics of a pendulum of variable length, rotations of a pendulum with elliptically moving pivot and twirling of a hula-hoop presented in three subsequent sections. The dynamics of these mechanical systems is described by similar equations and is studied with the use of common methods. The material of the chapter is based on publications of the authors [1-7] with the renewed analytical and numerical results. The methodological peculiarity of this work is in the assumption of quasi-linearity of the systems which allows us to derive higher order approximations by the averaging method. All the approximate solutions are compared with the results of numerical simulation demonstrating good agreement. Supplementary, in Appendix (section 5) we briefly presented the method of averaging with higher order approximations which is used in sections 2, 3, and 4.


Introduction
In this chapter we study three mechanical problems: dynamics of a pendulum of variable length, rotations of a pendulum with elliptically moving pivot and twirling of a hula-hoop presented in three subsequent sections. The dynamics of these mechanical systems is described by similar equations and is studied with the use of common methods. The material of the chapter is based on publications of the authors [1][2][3][4][5][6][7] with the renewed analytical and numerical results. The methodological peculiarity of this work is in the assumption of quasi-linearity of the systems which allows us to derive higher order approximations by the averaging method. All the approximate solutions are compared with the results of numerical simulation demonstrating good agreement. Supplementary, in Appendix (section 5) we briefly presented the method of averaging with higher order approximations which is used in sections 2, 3, and 4.

Pendulum with periodically variable length
Oscillations of a pendulum with periodically variable length (PPVL) is the classical problem of mechanics. Usually, the PPVL is associated with a child's swing, see Fig. 1. Everyone can remember that to swing a swing one must crouch when passing through the middle vertical position and straighten up at the extreme positions, i.e. perform oscillations with a frequency which is approximately twice the natural frequency of the swing. Among previous works we cite [8][9][10][11][12][13][14][15] in which analytical and numerical results on dynamic behavior of the PPVL were presented.
The present section is devoted to the study of regular and chaotic motions of the PPVL. Asymptotic expressions for boundaries of instability domains near resonance frequencies are derived. Domains for oscillation, rotation, and oscillation-rotation motions in parameter space are found analytically and compared with numerical study. Chaotic motions of the pendulum depending on problem parameters are investigated numerically. Here we extend our results published in [1][2][3][4] in investigating dynamics of this rather simple but interesting mechanical system.

Main relations
Equation for motion of the swing can be derived with the use of angular momentum alteration theorem, see [8][9][10][11]. Taking into account also linear damping forces we obtain where m is the mass, l is the length, θ is the angle of the pendulum deviation from the vertical position, g is the acceleration due to gravity, and t is the time, Fig. 1. It is assumed that the length of the pendulum changes according to the periodic law l = l 0 + aϕ(Ωt) > 0, (2) where l 0 is the mean pendulum length, a and Ω are the amplitude and frequency of the excitation, ϕ(τ) is the smooth periodic function with period 2π and zero mean value.
We introduce new time τ = Ωt and three dimensionless parameters where Ω 0 = g l 0 is the eigenfrequency of the pendulum with constant length l = l 0 and zero damping. In this notations equation (1) takes the form where the upper dot denotes differentiation with respect to new time τ.
Stability and oscillations of the system governed by equation (4) will be studied in the following subsections via analytically under the assumption that the excitation amplitude ε and the damping coefficient ε are small. For rotational orbits we will also assume the smallness of the frequency ω which means high excitation frequency compared with the eigenfrequency Ω 0 .

Instability of the vertical position
It is convenient to change the variable by the substitution η = θ (1 + εϕ(τ)) . ( 5 ) Using this substitution in equation (4) and multiplying it by 1 + εϕ(τ) we obtain the equation for η asη This equation is useful for stability study of the vertical position of the pendulum as well as analysis of small oscillations.
Let us analyze the stability of the trivial solution η = 0 of the nonlinear equation (6). Its stability with respect to the variable η is equivalent to that of the equation (4) with respect to θ due to relation (5). According to Lyapunov's theorem on stability based on a linear approximation for a system with periodic coefficients the stability (instability) of the solution η = 0 of equation (6) is determined by the stability (instability) of the linearized equation This equation explicitly depends on three parameters: ε, β and ω. Expanding the ratio in (7) into Taylor series and keeping only first order terms with respect to ε and β we obtain This is a Hill's equation with damping with the periodic function −(φ(τ)+ω 2 ϕ(τ)).I ti s known that instability (i.e. parametric resonance) occurs near the frequencies ω = k/2, where k = 1, 2, . . .. Instability domains in the vicinity of these frequencies were obtained in [16,17] analytically. In three-dimensional space of the parameters ε, β and ω,thesedomainsare described by half-cones where r k = 3 4 a 2 k + b 2 k is expressed through the Fourier coefficients of the periodic function ϕ(τ) Inequalities (9) give us the first approximation of the instability domains of the vertical position of the swing. These inequalities were obtained in [1] using different variables.
Note that each k-th resonant domain in relations (9) depends only on k-th Fourier coefficients of the periodic excitation function. Particularly, for ϕ(τ)=cos(τ), k = 1 we obtain a 1 = 1, b 1 = 0, and r 1 = 3/4. Thus, the first instability domain takes the form  The boundary of the first instability domain (k = 1) is presented in Fig. 2 by the solid red line demonstrating a good agreement with the numerically obtained instability domain which is marked black. These boundaries are also drawn in Fig. 5 and 6 by solid white lines. It is easy to see from (8) that for the second resonance domain (k = 2, ω = 1) the excitation function −(φ(τ)+ω 2 ϕ(τ)) is zero for ϕ(τ)=cos(τ). This explains why the second resonance domain is empty, and the numerical results confirm this conclusion, see Fig. 2. Inside the instability domains (9) the vertical position η = 0 becomes unstable and motion of the system can be either regular (limit cycle, regular rotation) or chaotic.
We can express Poincaré variables q(τ) and ψ(τ) via ϑ andθ from (15) as q 2 = ϑ 2 +θ 2 /ω 2 and ψ = arctan −θ ϑω . We differentiate these equations with respect to time and substitute expressions forθ,θ,a n dϑ in terms of q and ψ obtained from (12), and (15). Then, in the resonant case we have equations for the slow amplitude q(τ) and phase shift where f (ϑ,θ, τ)=ε f 1 (ϑ,θ, τ)+ε 2 f 2 (ϑ,θ, τ)+o(ε 2 ). System (16)-(17) has small right hand sides because we assumed that ω − 1/2 = O(ε). As a result of averaging in the second approximation, see (121) in the Appendix, we get the system of averaged differential equationṡ where Q and Z are the averaged variables corresponding to q and ζ. This system gives steady solutions forQ = 0,Ż = 0. Thus, besides the trivial one Q = 0, in the first approximation we obtain from system (18)- (19) expressions for the averaged amplitude and phase shift as where j = ...,−1,0,1,2,... and "arctan" givesthe major function value lying between zero and π;s u b i n d e x" {1}" denotes the order of approximation with which the corresponding variable is obtained. Solution of system (16)- (17) in the first approximation is Solution of system (16)- (17) in the second approximation is the following, see (119) and (121),  (19) in the second order approximation. Substitution of these expressions into (15) yields the second order approximate solution of (4) in the following form

Regular rotations
We say that the system performs regular rotations if a nonzero average rotational velocity exists Ve l o c i t y b is a rational number because regular motions can be observed only in resonance with excitation. Motion with fractional average velocity such as |b| = 1/2 in Fig. 4 a) is usually called oscillation-rotation. Let us first study monotone rotations, where velocityθ has constant sign and integer average value b, see Fig. 4 b) and c).
In order to describe resonance rotations of the PPVL we use the method of averaging [9,18,19] which requires rewriting (4) in the Bogolubov's standard form. For that reason we assume  that ε, β and ω are small parameters, ε being of order ω 2 ,andβ of order ω 3 , which makes the system quasi-linear.
We introduce the vector of slow variables x and the fast time s = |b| τ,w h er ex 1 = θ − b τ is the phase mismatch, x 2 = dθ ds is the velocity, x 3 = 1 + ε cos s |b| is the excitation. From here the dot denotes derivative with respect to the new time s. Thus, equations (4) takes the standard

Dynamics of a Pendulum of Variable Length and Similar Problems
where it is assumed that With the method of averaging we can find the first, second and the following order approximations of equations (25).
Resonance rotation domains of PPVL for various |b| are presented in Fig. 5. We see that greater values of relative rotational velocities |b| are possible for higher excitation amplitudes ε. Numerically obtained rotational regimes are depicted in Fig. 5 by color points in parameter space (ω, ε) with β = 0.05 and initial conditions θ(0)=π,θ(0)=0.05. Domains of these points are well bounded below by analytically obtained curves for corresponding |b|.

Rotations with relative velocity |b| = 1
It is the third order approximation of averaged equation where regular rotations with |b| = 1 can be observed, see Fig. 4 b). In the third order approximation averaged equations take the following form, see (122),Ẋ where X 1 and X 2 are the averaged slow variables x 1 and x 2 . Auxiliary variable x 3 = 1 + ε cos(s/b) has unit average X 3 = 1 and is excluded from the consideration. Excluding variable X 2 from the steady state conditionsẊ 1 = 0andẊ 2 = 0 in (26) we obtain the equation for the averaged phase mismatch X 1 Thus, it is clear from (27)  Stability of the solutions obtained from (27) was studied in [2]. There was found the condition for asymptotic stability cos(X 1 ) > 0. Hence, if inequality (28) is strict, then there are asymptotically stable steady solutions and unstable solutions Thus, we conclude that if the parameters satisfy strict inequality (28) there are two stable

Rotations with relative velocity |b| = 2
Rotations with higher averaged velocities |b| = 2, . . . correspond to higher excitation amplitudes ε. That is why we consider the coefficient ω being of order ε,a n dβ being of order ε 3 . With this new ordering we obtain the sixth order approximation of the averaged equations for |b| = 2, see (128), which have steady state solutions determined by the following equation From equation (32) we get that the domain of rotations with |b| = 2 in the parameter space has the following boundary condition depicted in Fig. 5 with a bold solid line System (31) has similar structure to system (26). That is why stability condition for its steady state solutions appears to be the same: cos(X 1 ) > 0. Hence, if inequality (33) is strict, we find from (32) that there are asymptotically stable steady solutions and unstable steady solutions Thus, as in the previous case, if the parameters satisfy strict inequality (33) there are two stable

Basins of attractions and transitions to chaos
In order to determine domains of chaos we calculate maximal Lyapunov exponents presented in Fig. 6. We recall that positive Lyapunov exponents correspond to chaotic motions. Note that chaotic motion includes passing through the upper vertical position, i.e. irregular oscillations-rotations. This is usually called tumbling chaos. We have observed two types of transition to chaos. The first type is when the system goes through the cascade of period doubling (PD) bifurcations occurring within the instability domain of the vertical position when the excitation amplitude ε increases, for example at ω = 0.5 in Fig. 7(a). The second type is when chaos immediately appears after subcritical Andronov-Hopf (AH) bifurcation when the system enters the instability domain of the lower vertical position of PPVL, for   We can see the change of the system dynamics in its route to chaos along ω = 0.5 in the bifurcation diagram shown in Fig. 7(a), where red points denote rotations with mean angular velocity equal to one excitation frequency (|b| = 1) and green points denote those equal to two excitation frequencies (|b| = 2). The domain with the most complex regular dynamics is surrounded by the red rectangle, where the system can have coexisting oscillations, rotations and rotations-oscillations. Fig. 8 have been plotted using program Dynamics [21]. These basins track the changes of the system dynamics in its route to chaos along ω = 0.67. In Fig. 8(a) the oscillatory attractor (limit cycle) coexists with stationary attractor (lower vertical position of PPVL). In Fig. 8(b) we can see the first emergence of two rotational attractors with counterrotations. This picture is in a good agreement with condition (18)  small which means that at ω = 0.67 the transition to chaos through subcritical AH bifurcation is the most typical. In the middle of Fig. 8(d) the manifold of dark blue points reveals a typical strange attractor structure. The strange attractor inherits the basin of attraction from disappeared stationary attractor.

Elliptically excited pendulum
Elliptically excited pendulum (EEP) is a mathematical pendulum in the vertical plane whose pivot oscillates not only vertically but also horizontally with π/2 phase shift, so that the pivot has elliptical trajectory, see Fig. 9. EEP is a natural generalization of pendulum with vertically vibrating pivot that is one of the most studied classical systems with parametric excitation. It is often referred to simply as parametric pendulum, see e.g. [9,[22][23][24][25][26][27] and references therein. Stability and dynamics of EEP have been studied analytically and numerically in [28][29][30]. Approximate oscillatory and rotational solutions for EEP are the common examples in literature [31][32][33][34] on asymptotic methods. Sometimes EEP is presented in a slightly more general model of unbalanced rotor [31][32][33], where the phase shift between vertical and horizontal oscillations of the pivot can differ from π/2. EEP is also a special case of generally excited pendulum in [35]. The usual assumption for approximate solution in the literature is the smallness of dimensionless damping and pivot oscillation amplitudes in the EEP's equation of motion. We could find only one paper [36], where oscillations of EEP with high damping and yet small relative excitation were studied.
In this section we study rotations of EEP with not small excitation amplitudes and with both small and not small linear damping. Our analysis uses the exact solutions for EEP with the absence of gravity and with equal excitation amplitudes, when elliptical trajectory of the pivot

Main relations
Equation of EEP's motion can be derived with the use of angular momentum alteration theorem, see [30]  where l is the distance between the pivot and the concentrated mass m; c is the viscous damping coefficient; θ is the angle of the pendulum deviation from the vertical position; t is time; g is gravitational acceleration at the angle δ with respect to the negative direction of the axis y.
It is assumed that the pivot of the pendulum moves according to the periodic law where X, Y,andΩ are the amplitudes and frequency of the excitation.
We introduce new time τ = Ωt and the following dimensionless parameters With this notation equation (36) with substituted (37) in it takes the following form w h e r ew eu s et h ef o r m u l aY cos 1 Here the upper dot denotes differentiation with respect to the new time τ.

Exact rotational solution when ε = 0 and ω = 0
Conditions ε = ω = 0 mean that we find the mode of rotation for the circular excitation X = Y with absence of gravity g = 0. In this case, we call equation (39) the unperturbed equation 14 Will-be-set-by-IN-TECḦ where constants θ 0 are defined by the following equality provided that |β| ≤ μ.
To investigate the stability of these solutions we present the angle θ as θ = θ 0 − τ + η,where η = η(τ) is a small addition, and substitute it in equation (40). Then linearizing (40) and using equality (42), we obtain the linear equation According to the Lyapunov stability theorem based on the linear approximation, solution (41) is asymptotically stable if all eigenvalues of linearized equation (43) have negative real parts. Which happens when all coefficients in (43) are positive due to the Routh-Hurwitz conditions. From conditions (44), assumption μ > 0 in (38), and equality (42), it follows for β > 0 that the solutions are asymptotically stable, while the solutions are unstable, where k is any integer number. For negative damping, β < 0, both these solutions are unstable. From now on we will assume that the following conditions are satisfied which ensure the existence of stable rotational solution (45) as it is seen from (42) and (44). Indeed, in order to guarantee asymptotic stability β should be not only positive, but also strictly less than μ because of the second condition in (44), which can be transformed to inequality μ cos(θ 0 )= μ 2 − β 2 > 0 with the use of the positive root for μ cos(θ 0 ) from (42).

Approximate rotational solutions when ε ≈ 0 and ω ∼ √ ε
We assume that values of ε and ω 2 are small of the same order of smallness, i.e. ε ∼ ω 2 ≪ 1, so we can introduce new parameter w = ω 2 /ε. One can deduce from (38) and current assumptions that either gravity g is small or the frequency of excitation Ω is high with such damping c and mass m so that damping coefficient β ∼ 1. All small terms are in the right-hand side of equation (39). To solve equation (39) we assume that general solution of equation (39) has the form θ = −τ + θ 0 + εθ 1 + ε 2 θ 2 + ...
Next the general solution is substituted into equation (39), where sines are expanded into the Taylor series with respect to ε. By grouping together the terms with the same powers of ε and equating to zero, the set of differential equations is obtained We have already found solution (41) for equation (49) in the previous section. Here we consider the same stable regular rotations 1:1 (with the period equal to the period of excitation) whose zero approximation is given by (45). Hence, θ 0 is a constant. Thus, equations (50) and (51) can be written in the following waÿ where we denote μ sin(θ 0 )=β and μ cos(θ 0 )= μ 2 − β 2 with the use of relation (42) and the second condition in (44).

First order approximation
In consequence of conditions (47) non-homogeneous linear differential equation (52) can be presented in the following form where

lower index denotes harmonics number. Equation (54) has a unique periodic solution
where . Thus, the solution for (39) in the first approximation can be 83 Dynamics of a Pendulum of Variable Length and Similar Problems written as follows where constant θ 0 is defined in (45).

Second order approximation
Equation (54) takes the following form where coefficients in the right-hand side are the following Periodic solution for equation (57) has the form Thus, second order approximate solution can be shortly written in the following form where constant θ 0 is defined in (45), function θ 1 in (55), and function θ 2 in (59). In Fig. 10 it is shown how first and second order approximate solutions approach the numerical solution.
We will use this assumption in the next section to apply the classical averaging technique to the problem of small damping β ∼ √ ε.

Approximate rotational solutions when
One can see in (38) that assumptions ω ∼ β ∼ √ ε are valid for the high frequency of excitation Ω ∼ 1/ √ ε with other parameters being of order 1. Another option is small gravity g ∼ ε along with small ratio c/m ∼ √ ε.
After change of variable θ = −τ + √ εϑ equation (39) takes the following form with small right-hand side, where we denoteβ = β/ √ ε and as in the previous section w = ω 2 /ε. With zero right-hand side equation (61)θ + μϑ −β = 0 would describe harmonic oscillations aboutβ/μ value with frequency √ μ. After Taylor's expansion of sines in the right-hand side of (61) about ϑ = 0 we obtain the following equation which describes oscillator with both basic and parametric excitations. To solve equation (62) we use the method of averaging [9,18,19]. For that purpose we write (62) in the standard form of first order differential equations with small right-hand sides. First, we use Poincaré variables q and ψ defined via the following solution of generating systemθ In Poincaré variables equation (62) becomes a system of first order differential equationṡ where small function f (τ, q, ψ)= √ ε f 1 (τ, q, ψ)+ε f 2 (τ, q, ψ)+o(ε) is the right hand side of (61), where meaning that f (τ, q, ψ)=O( √ ε). Our next assumption is that which means that excitation frequency is close to the first resonant frequency of basic excitation component sin(τ − δ) and to the first resonant frequency of parametric excitation component cos(2τ) in equation (62). Thus, system (64) is transformed by ψ = ζ + τ to the standard forṁ with small right-hand side, where new slow variable ζ is often referred to as phase mismatch.
In the second approximation so called averaged equations can be obtained from the system of equations (67) and (68) as follows, see (121) in the Appendix, where Q and Z are the averaged variables corresponding to q and ζ.

First order approximation
Stationary solutions (Q = 0,Ż = 0) of (69)-(70) in the first approximation are the following which does not contain higher harmonics observed numerically. That is why we need to proceed to the second order approximation.

Second order approximation
In the second approximation averaged equations can be obtained as described in Appendix. Stationary solutions (Q = 0,Ż = 0) of (69), (70) in the second approximation can be found numerically or with the absence of gravity (ω = 0) analytically. Solution of system (67), (68) in the second approximation is the following, see (119), Substitution of these expressions into (63) yields the second order approximate solution of (61), which after changes of variable θ = −τ + √ εϑ and parameters w = ω 2 /ε,β = β/ √ ε results in the approximate solution of the original equation (39) Agreement of solution (75) with the numerical experiment is shown in Fig. 11. We see that the amplitude of angular velocity oscillations is much higher than that for not small β in Fig. 10.

Twirling of a hula-hoop
A hula-hoop is a popular toy -a thin hoop that is twirled around the waist, limbs or neck. In recent decades it is widely used as an implement for fitness and gymnastic performances . 2 To twirl a hula-hoop the waist of a gymnast carries out a periodic motion in the horizontal plane. For the sake of simplicity we consider the two-dimensional problem disregarding the vertical motion of the hula-hoop. We assume that the waist is a circle and its center moves along an elliptic trajectory close to a circle.
Previously considered was the simple case in which a hula-hoop is treated as a pendulum with the pivot oscillating along a line, see [37,38]. The stationary rotations of a hula-hoop  excited in two directions have been studied by an approximate method of separate motions in [33]. The similar problem of the spinner mounted loosely on a pivot with a prescribed bi-directional motion has been treated numerically and experimentally in [39].
Here we derive the exact solutions in the case of a circular trajectory of the waist center and approximate solutions in the case of an elliptic trajectory. We also check the condition of keeping contact with the waist during twirling. There is also a rolling resistance due to the waist deformation (right).

Main relations
We assume that the center O ′ of a gymnast's waist moves in time according to the elliptic law x = a sin ωt, y = b cos ωt with the amplitudes a, b and the excitation frequency ω > 0, Fig. 12.
The equations of motion in the waist-fixed coordinate system take the following form where θ is the rotation angle around center of mass C, I C = mR 2 is the central moment of inertia of the hula-hoop, ϕ is the angle between axis x and radius CO ′ , r is the radius of the waist, m and R are the mass and radius of the hula-hoop. Equation (76) describes change of angular momentum due to linear viscous damping with coefficient k, rolling drag (rolling resistance) with coefficient d, and the tangential friction force F T between the waist and the hoop. Equations (77) and (78) describe the motion of the hula-hoop in the longitudinal and transverse directions to the radius CO ′ ,whereN is the normal reaction force of the hula-hoop to the waist. Equations (77) and (78) contain additional inertial forces since the waist-fixed reference system is noninertial.
Assuming that slipping at the point of contact is absent we obtain the kinematic relation We exclude from equations (76) and (77) the force F T and with relation (79) obtain the equation of motionφ From equation (78) we find the normal force and imply the condition N > 0as which means that the hula-hoop during its motion keeps contact with the waist of the gymnast.
We introduce new time τ = ωt and non-dimensional parameters where γ and δ are the damping and rolling resistance coefficients, μ and ε are the excitation parameters. Relation between μ and ε determines the form of ellipse -the trajectory of the waist center. For ε = μ the trajectory is a line, and for ε = 0 it is a circle. Then equation (80) and inequality (81) take the form where the dot means differentiation with respect to the time τ.

Stability of the exact solutions
Let us investigate the stability of the obtained solutions. For this purpose we take the angle ϕ in the form ϕ = τ + ϕ 0 + η(τ) where η(τ) is a small quantity, and substitute it into equation (85). Taking linearization with respect to η and with the use of (87) we obtain a linear equation According to Lyapunov's theorem on the stability based on a linear approximation [17] solution (86), (88) is asymptotically stable if all the eigenvalues of linearized equation (89) have negative real parts. From Routh-Hurwitz criterion [17] we obtain the stability conditions as γ + 2δ > 0, μ (sin ϕ 0 + 2δ cos ϕ 0 ) < 0.

Approximate solutions
Let us find approximate solutions for the case of close but not equal amplitudes a ≈ b.F o r the sake of simplicity from now on we will keep δ = 0 and assume that a ≥|b| which means ε ≥ 0, μ ≥ 0. Taking ε as a small parameter we apply perturbation method assuming that the exact solution ϕ s (τ) of (83) can be expressed in a series ϕ s (τ)=τ + ϕ 0 + εϕ 1 (τ)+o(ε).
After substitution of series (94) in (83) and grouping the terms by equal powers of ε we derive the following chain of equations ε 0 : γ + μ cos(ϕ 0 )=0 (95) Taking solution of equation (95) for μ > 0 corresponding to the stable unperturbed solution (86), (91) we write equation (96) as It has a unique periodic solution where constants C and D are defined as follows (100)  We see that the approximate solution ϕ(τ)=τ + ϕ 0 + εϕ 1 (τ) with (97), (99), (100) differs from the exact solution (86), (91) of the unperturbed system by small vibrating terms of frequency 2, see (99). Note that the approximate solutions were obtained with the assumption that the excitation amplitudes and damping are not small.

Condition for hula-hoop's contact with the waist
The condition of twirling without losing contact (84) takes the following form Conditions (103), (104) imply restrictions to ε, i.e. how much the elliptic trajectory of the waist center differs from the circle.

Comparison with numerical simulations
In Fig. 13 the approximate analytical solution is presented and compared with the results of numerical simulation for the case when the excitation parameter μ and damping coefficient γ are not small.

Small excitation amplitudes and damping
It is interesting to consider the case when the excitation amplitudes and damping coefficient are small having the same order as ε. Then we introduce new parametersμ = μ/ε andγ = γ/ε and assume that the solution has the form where ρ is the angular velocity of rotation, and the functions ϕ 0 (τ), ϕ 1 (τ) are supposed to be bounded.

Clockwise rotation
For clockwise rotation ρ = 1, see Fig. 14 a), from equations (105), (106), and (107) we obtain in the first approximation the solution where function ϕ 1 (τ) is found up to the addition of a constant, which we have set to zero for determinacy. Thus, we let only ϕ 0 contain a constant term of the solution. The first expression in (109) is the special case of (99), (100) when μ = γ = 0.
To verify the stability conditions for solution (108) we use damped Mathieu-Hill equation (102) and for the case of small damping and excitation amplitudes get γ > 0, sin ϕ 0 < 0, see [17]. These conditions are similar with inequalities (90) derived for the undisturbed exact solution. Thus, the stable solution (108) with ϕ 0 = − arccos(−γ/μ)+2πn exists for For solution (108) condition (84) of keeping contact in the first approximation reads and holds true for sufficiently small ε.
Conditions (115) in physical variables take the form meaning that the trajectory of the waist should be sufficiently prolate. Coexisting clockwise and counterclockwise rotations are illustrated in Fig. 14.

Comparison with numerical simulations
In Fig. 15 the approximate analytical solutions for rotations in both directions are presented and compared with the results of numerical simulation for the case of small excitation parameters μ, ε and the damping coefficient γ.T h e v a l u e s o f μ, ε correspond to the dimensional parameters a = 15cm, b = 10cm, r = 10cm, R = 50cm.

Conclusions
In section 2 we showed that the pendulum with periodically varying length exhibits diversity of behavior types. We recognized that the analytical stability boundaries of the vertical position of the pendulum and the frequency-response curve for limit cycles are in a good agreement with the numerical results. The second resonance zone appeared to be empty. The stability conditions of limit cycles are derived based on direct use of Lyapunov's theorem on stability of periodic solutions. We found numerically regular rotation, oscillation, and rotation-oscillation regimes with various periods and mean angular velocities of the pendulum including high-speed rotations and rotations with fractional relative velocities (it is rotation-oscillation regime when the pendulum makes regular sequence of rotations in both directions). We derived analytically the conditions for existence of regular rotation and oscillation regimes which agree with the numerical results. Domains for chaotic motions are found and analyzed numerically in the parameter space via calculation of Lyapunov exponents and bifurcation diagrams. Basins of attractions of different regimes of the pendulum motion were plotted and analyzed.
In section 3 we studied the planar rotational motion of the pendulum with the pivot oscillating both vertically and horizontally when the trajectory of the pivot is an ellipse close to a circle. The analysis of motion was based on the exact rotational solutions in the case of circular pivot trajectory and zero gravity. The conditions for existence and stability of such solutions were derived. Assuming that the amplitudes of excitations are not small while the pivot trajectory has small ellipticity the approximate solutions were found both for large and small linear damping. Comparison between approximate and numerical solutions was made for different values of the damping parameter demonstrating good accuracy of the method involved.
Finally, in section 4 we assumed that the waist of a sportsman twirling a hula hoop is a circle and its center moves along an elliptic trajectory close to a circle. We studied the system with both small and not small linear viscous damping as well as with some rolling resistance. For the case of the circular trajectory, two families of the exact solutions were obtained, similar to those in section 3. Both of them correspond to twirling of the hula-hoop with a constant angular speed equal to the speed of the excitation. We showed that one family of the solutions is stable, while the other one is unstable. These exact solutions allowed us to obtain the approximate solutions for the case of an elliptic trajectory of the waist. An interesting effect of inverse twirling was described when the waist moves in opposite direction to the hula-hoop rotation. It is shown that the approximate analytical solutions agree with the results of numerical simulation.