Rotating solutions and stability of parametric pendulum by perturbation method

Rotating solutions of a parametrically driven pendulum are studied via a perturbation method by assuming the undamped unforced Hamiltonian case as basic solution, and damping and excitation as small perturbations. The existence and stability of the harmonic solution are determined analytically. First-order terms are mainly considered, but the extensions to higher-order terms, as well as to subharmonic rotations, are straightforward. Setting up the analysis in the phase space instead of in the physical space has facilitated development of simple but comprehensive formulas. Comparison with numerical simulations shows that the analytical results are very effective in accurate predictions, even for fairly large amplitudes, and in detecting the saddle–node bifurcation where the rotating solution is born. A limited accuracy is observed, to ﬁrst order, in detecting the period–doubling bifurcation where it loses stability, and an explanation is proposed.


Introduction
The mathematical pendulum is an archetype for studying nonlinear dynamics and complex phenomena, in various fields of science and engineering (e.g., electrical circuits [1] and charge density waves [2]). In spite of the fact that the first scientific interest in pendula dates back centuries ago with seminal contributions from Galileo, Huygens and Foucault, pendula dynamics is still actively investigated, e.g., Refs. [3][4][5][6][7][8]. The aim of the current studies is to discover new motions [9], unveil new phenomena [10], test out various analytical [11,12] and numerical [13] techniques, and exploit the developed models and techniques for practical applications, e.g., Ref. [14].
where m is the mass, l the distance between the centre of gravity and the pivot, c the viscous damping coefficient, t the time, g the acceleration of gravity, z(t) the vertical displacement of the pivot and x the angle measured from the downward vertical position. The natural frequency of the small oscillations around the unperturbed equilibrium point x ¼ 0 is o 0 ¼ ffiffiffiffiffiffi ffi g=l p , so Eq. (1) becomes Assuming zðtÞ ¼ Àz cosðōtÞ we get while introducing the dimensionless time, damping coefficient, frequency and amplitude of the parametric excitation, respectively, gives the following dimensionless equation of motion [10]: where the dot represents differentiation with respect to t.
In this work we are interested in the rotating solutions of Eq. (3) (or (5)), which cannot be determined exactly, so that we have to resort to approximate solutions. These can be obtained by applying the perturbation method directly to the governing Eq. (3), after an appropriate choice of the smallness parameter. This has been done in Section 4 of Ref. [29], where an approximate solution is obtained for o very large (it is asymptotically exact for o-N, see their Fig. 10), without determining the stability.
Herein, those results are improved by following a different approach which allows for a more comprehensive analysis. To this aim, before applying the perturbation method, we reformulate problem (5) as an integral equation, which can be obtained by multiplying Eq. (5) by _ x and rearranging it in a classical way d dt For a given (initial) time t i , let the (initial) position and velocity of the pendulum be x i ¼ x(t i ) and _ x i ¼ _ xðt i Þ, so that Eq. (6) can be expressed in the form cosðorÞ sin½xðrÞ dxðrÞ dr dr.
For rotating orbits, _ xðtÞ40, 8t, so that the inverse function t ¼ t(x) of x(t) exists. It is then possible to define unambiguously the function _ x ¼ _ xðtÞ ¼ _ x½tðxÞ ¼ f ðxÞ, which is central in the following development. We have which explicitly provides t ¼ t(x) when f(x) is known. For later use, we further note that _ The previous considerations permit rewriting the differential Eq. (7) as the following integral equation: where the unknown is f(x). Apart from the physical characteristics h, p and o, the solution f(x) depends on the initial conditions x i , _ x i and t i . In the following, it is shown that only two of them are independent, while the third one can be assumed equal to 0 without loss of generality.
We are interested in investigating rotating solutions periodic in space and time, namely, 3 where T is the temporal period. From Eq. (8) we obtain the ''synchronicity'' equation while from Eq. (9) we obtain the ''periodicity'' equation Remark. Eq. (12) can be rewritten in the form Let the angle a be defined by so that Eq. (13) becomes Consequently, a solution exists if and only if p cr represents a lower bound for the existence of periodic solutions, and it is attained for ot i ¼ a+p/2. Indeed, above this curve two solutions exist, while none below. This implies that p cr is the critical threshold corresponding to the birth of the rotating periodic solution through a SN bifurcation [4].
Eqs. (11) and (12) allow to determine two independent unknown initial conditions. Actually, two different strategies can be pursued to develop analytical solutions based on spatial or temporal periodicity. In this work, use is made of spatial periodicity assuming x i ¼ 0. However, the problem cannot be solved in closed form, and perturbation analysis is performed assuming that the damping and the excitation are small, h ¼ eh 1 and p ¼ ep 1 (e is a non-dimensional small parameter) [6]. This means that the considered case is close to the Hamiltonian problem, for which the solutions are well known [5].
In the following we look for harmonic response, T ¼ 2p/o, i.e., for rotating solutions having the same period as the excitation. The adopted procedure can be easily extended to subharmonic motions.

Perturbation analysis
We assume x i ¼ 0, namely, we look for the periodic solution starting from x ¼ 0 at time t ¼ t i and ending at x ¼ 2p at t ¼ t i +T. The governing equation, the synchronicity and periodicity conditions become, respectively, where c i ¼ ot i can also be interpreted as the phase of the excitation. According to the perturbation approach, a solution of Eqs. (18)- (20) has the following generic form: Inserting Eq. (21) in Eq. (18) and equating the coefficients of each e-power to zero we get . . . ; can be determined by the synchronicity 0 : 1 : 2 : and periodicity conditions It is straightforward to prove that Eqs. (28) and (29) guarantee periodicity of f 1 (x) and f 2 (x) respectively.

First-order terms
The solution is obtained in two steps. First, from the first synchronicity condition (25) 2p , which is depicted in Fig. 1. KðxÞ ¼ R p=2 0 ½1 À x 2 sin 2 ðsÞ À1=2 ds denotes the complete elliptic integral of the first kind [26]. The function tends to the asymptotic limit _ x 0 ¼ o for o-N and satisfies the inequality _ x 0 ðoÞ42, 8o, which guarantees that the solution is rotating [5].
In the second step, we determine the phase (or initial time) c 0 by the first periodicity condition (28), which can be rewritten as and is also the first-order expansion of Eq. (13). As can be seen from Appendix A, the following relationships are satisfied: On the basis of Eq. (32), Eq. (31) becomes and c 0 ¼ c 0 (o,h 1 /p 1 ) can be determined by From Eq. (34) we see that there are two solutions for and no solutions for p 1 op 1,SN (o), so that p 1,SN (o) represents the critical threshold where a rotating solution is born via a SN bifurcation. Just above this threshold, there are two rotating solutions, one stable and another unstable (see Section 4 for more details). It is worth to note that at the bifurcation point we have sin(c 0 ) ¼ 1 and the phase is c 0 ¼ p/2, while above it one solution has c 0 op/2 and the other c 0 4p/2. Naturally, p 1,SN is the first-order approximation of p cr as given in Eq. (17). The function g(o) is depicted in Fig. 2, along with the following asymptotic limits: so that it can be very useful in practical applications [27]. The function g(o) is the same as the one presented by Koch and Leven [12, Eq. (2.16)], where this result is obtained by a Melnikov analysis for subharmonic orbits. Comparing that approach with the present one, the following points can be drawn: Both are perturbation techniques based on the same choice of the small quantities. Melnikov analysis provides only the critical threshold, while here we also obtain the initial conditions x i and t i (to be used in numerical simulations) as well as the entire solution (by inverting Eq. (8)).
The present approach allows for easier extension to higher-order terms (see Section 3.2 and, for example, Ref. [28]).
As the Melnikov approach does not provide information on the stability of rotating solutions, Koch and Leven added to the Melnikov analysis an approach suggested by Guckeneihmer and Holmes [4, Section 4.7] using action-angle coordinates and averaging. However, while obtaining the SN bifurcation curve which triggers the rotating solutions, they did not mention the upper bound determined by the PD curve, where the rotating solutions lose stability. This is now appropriately discussed in Section 4.
In Ref. [12] the Melnikov analysis is reported also for other subharmonic The extension of the present approach to these solutions is straightforward and is not reported here.

Higher-order terms
Higher-order terms _ x k and c k , k40, can be obtained by the same steps as the first-order ones. First, we determine the velocity _ x k from the synchronicity condition of order k, then the phase c k from the periodicity condition of order k+1. These equations are linear in their single unknown function, and their solutions are simple, even if they involve the computation of complex integrals. For example, from Eq. (26) we get _ With this value of _ x 1 it is possible to compute f 1 (x) by Eq. (23), so that from Eq. (29) we get c 1 ¼ c 1 (o,h 1 ,p 1 ):

Stability
To study the stability of the periodic solution we use the Poincare`section S ¼ fðx; _ x; tÞ 2 S Â IR Â S such that x ¼ 0}. In fact, we take advantage of the space 2p-periodicity of the pendulum which permits to identify the sections at x ¼ 2pn, 8n.
Let ð _ x i ; t i Þ be a point of S and let f ðx; _ x i ; t i Þ be the solution of the integral Eq. (18) corresponding to ð _ x i ; t i Þ. The image ð _ x 2p ; t 2p Þ of this point on S by the 2D Poincare`map S-S is given by The periodic solution corresponds to the fixed point of the map (40). In particular, from _ x 2p ¼ _ x i we recover the periodicity condition, while from t 2p ¼ t i we get the synchronicity condition. To study its stability we have to calculate the Jacobian of the map and to compute the eigenvalues l 1,2 , which are the solutions of l 2 À ltrðJÞ þ detðJÞ ¼ 0.
The periodic solution is asymptotically stable as long as both jl 1 jo1 and jl 2 jo1. The loss of stability can occur by a SN bifurcation when one eigenvalue is equal to 1, or by a period-doubling (PD) bifurcation when one eigenvalue is equal to À1. The Hopf bifurcation corresponding to the crossing of the unit circle in a complex point is not possible in the present case, as it is shown later.
Next each component of J is computed and details are given in Appendix B. After straightforward simplifications involving the use of Eq. (34), we obtain simple explicit formulae for the Jacobian, its trace and determinant Eq. (44) can be further simplified by the relations: where EðxÞ ¼ R p=2 0 ½1 À x 2 sin 2 ðsÞ 1=2 ds is the complete elliptic integral of the second kind [26]. 8 The fact that det (J)o1 proves that there are no Hopf bifurcations. Eq. (42) with the expressions (44) and (45) allows to determine the eigenvalues of the stable and unstable periodic solutions. A typical example of how they vary under increasing excitation amplitude is shown in Fig. 3, where both solutions are born from the SN bifurcation. The eigenvalues of the saddle always remain real, and one is greater than 1, one is between 0 and 1. The eigenvalues of the node are initially real and smaller than 1. They rapidly become complex (the node becomes a focus). Next change occurs at the point where they return to be real and negative, but still greater than À1. Slightly further, one of the eigenvalues becomes smaller than À1, and the node loses stability through a PD bifurcation. The stability region is bounded from below by ep 1,SN and from above by ep 1,PD . Note that ep 1,SN ¼ 0.1g(2) ¼ 0.4168 is just the value given by Eq. (35), while ep 1,PD ¼ ep 1,SN d(2,0.1) ¼ 0.4168*1.814 ¼ 0.7563 (see forthcoming Eq. (51)).
In the following we will generalize the previous result by systematically determining the stability region. Indeed, we will compute the SN and PD critical thresholds for every value of o and eh 1 . Before proceeding further, we note that in the unperturbed case e ¼ 0, i.e., no dissipation and excitation, the eigenvalues are coincident and equal to 1. If the perturbation is small, the eigenvalues vary but remain close to 1. This has strong implications in terms of reliability of the perturbation analysis. In fact, for the SN bifurcation occurring at l ¼ 1, one eigenvalue remains equal to 1 while the other becomes smaller than 1, and this can be easily achieved for small e. However, for the PD bifurcation, the eigenvalue has to reach, for growing e, the value À1. Since this is ''far'' from the starting point (for e ¼ 0) l ¼ 1, it is expected that ''large'' values of e would be required to achieve this bifurcation point. But in this e-range the perturbation analysis may fail, if we are out of the radius of convergence of the e-series, or it may require several terms of the series to provide adequate results. In other words, we expect that the first-order perturbation analysis developed in this section will provide good results as far as the SN bifurcation is concerned and poor results for the PD bifurcation. In the latter case, higher-order e-terms are likely required, but this is out of the scope of this work.
The namely, The different reliability of the quantitative predictions of SN and PD discussed above is evident also from a mathematical point of view. In fact, the equation for the SN bifurcation ð0 ¼ J 0 21 J 1 12 þ Á Á ÁÞ has no e 0 term, and therefore can be satisfied by small e, while the equation for the PD bifurcation ð0 ¼ 4 þ ½À2ðJ 1 11 þ J 1 22 Þ þ J 0 21 J 1 12 þ . . .Þ has e 0 term, and requires ''large'' values of e to be satisfied, which contradicts the philosophy of the perturbation analysis. Higher-order e-terms have to be considered to improve the predictions on PD.
To determine the curve p PD ¼ p PD (o,h), we compute tan(c 0 ) from Eq. (48), so that Note that if o4eh 1 p the right-hand side of Eq. (49) is positive, i.e., c 0 op/2. As eh 1 ¼ h is usually small, this is the common case. Substituting Eq. (50) in Eq. (34) we get which provides the first-order approximation of p PD (o, h). The function d(o, eh 1 ) representing the ratio p 1,PD /p 1,SN between the (first-order) critical thresholds for PD and SN bifurcations is depicted in Fig. 4 for various values of eh 1 . Fig. 4 shows that for small values of the excitation frequency the PD occurs just above the SN, and that the stable region vanishes. For large values of o, on the other hand, the stable region increases in magnitude and tends to become infinite, the smaller the damping, the larger the distance between the SN and PD. Thus, for high frequencies it is relatively easy to have a harmonic stable rotating solution, while for small frequencies it is practically impossible.
Remark. In the previous developments, we assumed x i ¼ 0 and took advantage of the space periodicity of the system. It is also possible to assume t i ¼ 0 and use the temporal periodicity, i.e., the standard stroboscopic Poincare`map. This approach has been performed using the perturbation technique, too. Although it requires more laborious computations, it provides the same results as those previously obtained.

Comparison with numerical results
The results analytically obtained in the previous Sections 3 and 4 are here compared with some numerical simulations. We start by comparing in Fig. 5 the analytically and numerically obtained functions f(x) in the case o ¼ 2, h ¼ eh 1 ¼ 0.1 for two different values of p ¼ ep 1 , where one is located just after the numerically observed SN bifurcation (see Fig. 6c), the other one is located just before the numerically observed PD bifurcation (Fig. 6c). These examples are chosen at the boundaries of the actual stability region of the rotating solution to analyse the possible difference in the quality of the obtained prediction as the value of p increases. The numerical solution is obtained by integrating the differential Eq. (5) using standard Runge-Kutta techniques, and the results are parametrically plotted as the curve ½xðtÞ; _ xðtÞ after eliminating the transient. A very good correspondence is observed between the analytically predicted function f 0 (x)+ef 1 (x) and the numerical function f(x) in both cases. These examples demonstrate that the perturbation analysis is very effective in determining the rotating solution for a large range of parameters, and thus we believe that it is overall reliable. Indeed, the agreement in Fig. 5b is somehow surprising, because here the excitation amplitude is not small and one could not expect good correspondence up to this ''high'' value of ep 1 .
As expected, the two cases shown in Figs. 5a and b confirm that with increasing amplitude of the excitation, the e term plays more and more significant role and the difference between the functions f(x) and f 0 (x) increases.
To check if the effectiveness of the perturbation approach holds also with respect to the stability of the rotating solutions, four numerical bifurcation diagrams, calculated for o ¼ 2 and different values of damping, are reported in Fig. 6. In all cases, the rotating solutions are born by the SN bifurcation, and have a certain range of stability, which is lost through the PD bifurcation. The PD bifurcation is followed by classical PD cascade ending in tumbling chaos after a boundary crisis (including Fig. 6d, where the route to chaos is out of the picture) [10]. Thus, we conclude that the described bifurcation scenario is damping independent.
Qualitatively, we find a good agreement between theoretical and numerical analysis. In particular, the range of stability, which is of principal interest for the present work, is confirmed to be confined by a SN from below and by a PD from above, irrespective of the damping. For a quantitative comparison, on the other hand, we refer to Table 1, which reports the analytical and numerical SN and PD thresholds, and shows that we have 2π 0 x Though the latter discrepancy is certainly large, the overall comparison indeed agrees with the expectations, according to which the first-order perturbation analysis works well for the SN and poorly for the PD (see above). To improve theoretical results in the latter case, higher-order terms in the e-series should be used.  However, it is worth to note that the first-order perturbation analysis correctly captures the fact that the critical threshold for SN linearly depends on the damping, while the critical threshold for PD varies scarcely with h. In particular, the ratio p PD /ep 1,PD is practically constant for damping up to h ¼ 0.1, showing that the theoretical analysis describes the qualitative dependence on h correctly, even if it fails in quantitative predictions. For larger damping (h ¼ 0.3), that ratio reduces, and the approximation is slightly more accurate, according to the fact that the smaller the difference between p SN and p PD , the closer the ep 1,PD is to p PD .
In the last numerical simulation we fix the damping to the intermediate value h ¼ eh 1 ¼ 0.05 and vary the excitation frequency to determine the numerical stability region in the excitation frequency-amplitude parameters space, which is reported in Fig. 7 together with its theoretical counterpart.
From a qualitative point of view, the main conclusion which can be drawn from Fig. 7 is that the theoretical analysis correctly describes the principal features of the stability region, in particular its vanishing for small frequencies and its rapid growing in magnitude for large frequencies. From a quantitative point of view, on the other hand, we can say that the theoretical prediction ep 1,SN (o) of p SN (o) is always very accurate, while the first-order theoretical prediction ep 1,PD (o) of p PD (o) is accurate for small frequencies but fails as the frequency increases. Overall, the PD threshold is underestimated. In this respect, since it marks the start of the transition route to possibly unwanted tumbling chaos, it turns out to be conservative.
The general conclusions previously drawn also hold for different damping values, as confirmed by other not reported here numerical simulations.

Conclusions
In this paper the analytical approximation of the rotating solutions of the parametrically excited pendulum has been obtained by the perturbation method. Existence and stability of this solution have been determined. The perturbation formulation has been implemented by working with the trajectory _ x ¼ f ðxÞ instead of the time history x(t), which allows to obtain a comprehensive but simple picture of the dynamic responses. An explicit formula for the first-order approximation of f(x) has been derived. The higher-order terms can be obtained by solving a sequence of linear equations whose coefficients contain some integrals to be computed numerically. Relative to existing work (e.g., Refs. [12,29]), the present approach has two main contributions, namely (i) it provides the initial conditions as part of the rotating solutions and (ii) it affords a stability analysis, up to first order, that, in addition to obtaining a lower bound (SN bifurcation) that triggers the basic period-T rotating solutions in the excitation parameter plane, also provides an upper bound (PD bifurcation) through which these solutions lose stability.

13
To check the robustness of the analytical results, comparisons with numerical simulations of trajectories, of bifurcation diagrams and of the stability region chart have been made. The main conclusions are the following: There is an excellent agreement between analytically and numerically determined trajectories, even for large amplitudes.
The perturbation method is very effective in predicting, both qualitatively and quantitatively, the SN bifurcation, namely, the lower critical threshold where the rotating solutions are born.
The qualitative properties of the PD bifurcation bounding the stability region from above are captured by the perturbation method.
The first-order terms used in this work are not able to provide accurate quantitative predictions of the PD bifurcations. An explanation of this fact is proposed.
In this last respect, the use of higher-order terms in the asymptotic framework is suggested, and this constitutes the most natural development of the present work. It can also be continued in some other directions. For example, considering other excitations (horizontal and/or rotational motion of the pivot, instead of vertical) deserves an interest in practical applications and should provide further insight into the complex dynamical behaviour of this archetypal mechanical system.