# Rotating solutions and stability of parametric pendulum by perturbation method

S Lenci, E Pavlovskaja, G Rega, M Wiercigroch

► To cite this version:

S Lenci, E Pavlovskaja, G Rega, M Wiercigroch. Rotating solutions and stability of parametric pendulum by perturbation method. Journal of Sound and Vibration, 2008, 310, pp.243 - 259. 10.1016/j.jsv.2007.07.069. hal-01403850

HAL Id: hal-01403850

https://hal.science/hal-01403850

Submitted on 28 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rotating solutions and stability of parametric pendulum
by perturbation method

S. Lenci<sup>a</sup>, E. Pavlovskaia<sup>b</sup>, G. Rega<sup>c</sup>, M. Wiercigroch<sup>b</sup>

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 first order, in detecting the period-doubling bifurcation where it loses stability, and an explanation is proposed.

## 1. Introduction

The mathematical pendulum is an arch type 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–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].

Over the last five decades various aspects of pendulum dynamics have been studied, ranging from simple harmonic and subharmonic motions [6,15–17] and complex dynamics [18,19] to control regular inverted position [19] and chaos [20], using analytical [6,15], numerical [10,17] and experimental methods [18,21,22]. Among various dynamical responses, rotating solutions under parametric excitation, i.e., *vertical motion* of the pivot, have recently attracted a great deal of attention, which is briefly summarized below.

<footnote><sup>a</sup>Department of Architecture, Buildings and Structures, Polytechnic University of Marche, Ancona, Italy</footnote>

<footnote><sup>b</sup>Centre for Applied Dynamics Research, Department of Engineering, University of Aberdeen, Scotland, UK</footnote>

<footnote><sup>c</sup>Department of Structural and Geotechnical Engineering, University of Rome “La Sapienza”, Rome, Italy</footnote>
Koch and Leven [12] applied the Melnikov's method to detect the saddle-node (SN) bifurcation where the harmonic and subharmonic rotating solutions are born for increasing excitation amplitude. They derived a closed form expression for the lower boundary of region in the excitation frequency-amplitude parameters space, $(\omega, p)$, in which the stable rotating solutions exist. They also applied an averaging method to study stability of the solutions originating from the SN bifurcation.

Capecchi and Bishop [23] studied rotating solutions analytically by the harmonic balance method, and obtained the asymptotic limit for $\omega \to \infty$. They also determined numerically the frequency and amplitude response curves, including the secondary stability region above the principal one. From these investigations it is clear that the stable region is between the SN bifurcation curve, from below, and the period-doubling (PD) bifurcation curve, from above.

Clifford and Bishop [24] investigated numerically rotating solutions with the aim of producing “an albeit partial classification of rotating orbits”. They introduced the term “tumbling chaos” to denote an erratic alternation of rotating clockwise, rotating anti-clockwise and oscillating motion. They also used the method of close returns to detect numerically certain rotating solutions which they classified according to the period motion and to the number of rotations in the period. More recently, Garira and Bishop [9] extended this classification identifying new types of rotating solutions, and classified them, together with the previously known ones, as purely rotating, oscillating rotating, straddling rotating and large amplitude rotating.

The studies of Szemplinska-Stupnicka and co-workers [10,17] are also entirely numerical and use a combination of basins of attraction, bifurcation diagrams and attractor-manifolds phase portraits. They are aimed at investigating global (i.e., homoclinic or heteroclinic) bifurcations responsible for onset of complex transient and/or steady-state dynamics and various other aspects including fractal basin boundaries and coexistence of rotating solutions with other (even non-rotating) attractors.

Most recently, Xu et al. [25] have performed extensive numerical simulations and produced various parameter space plots for different sets of initial conditions and damping. Their results allow follow the development of attractors in the excitation amplitude-frequency parameters space. They also report examples of various kind of rotating solutions.

Despite this great interest in pendulum dynamics, the rotating solutions have not been extensively investigated analytically, apart from the initial work by Koch and Leven [12] and a recent one by Xu and Wiercigroch [29]. This paper is specifically intended to fill this gap with the main aim to develop a robust approximate analytical solution. As exact solutions can be obtained only in the conservative case, the perturbation method will be applied by assuming that damping and excitation amplitudes are small quantities with respect to the Hamiltonian unperturbed case.

The current paper is focused on the development of solutions, which rotate with the same frequency as the excitation; however, the technique we adopt, can be also applied to subharmonic rotations. The objective of the work is to obtain formulas which—albeit comprehensive—are as simple as possible, in view of practical applications. The analysis is conducted in the phase space $(x, \dot{x})$, where the trajectory $\dot{x} = f(x)$ is unknown, instead of in the physical space $(t, x)$, where the time history $x(t)$ is unknown. In principle, the same perturbation scheme can be applied to the time history formulation; however, this would require the use of *incomplete elliptic integrals complicating the computation* [26].

The paper is organized as follows. In the next section, the dimensional and non-dimensional equations of motion are presented and transformed into an integral equation. The periodicity and synchronicity conditions are also formulated there. An approximate solution of the integral equation is then obtained in Section 3 by a perturbation method, while its stability is studied in Section 4. Section 5 gives comparisons with numerical results, while Section 6 contains some conclusions.

## 2. Problem formulation

The motion of a vertical parametric pendulum can be described by the following second-order ordinary differential equation [10]:

$$ml^2 \frac{d^2 x}{d\tau^2} + c \frac{dx}{d\tau} + ml \left( \frac{d^2 z}{d\tau^2} + g \right) \sin(x) = 0, \quad (1)$$
where *m* is the mass, *l* the distance between the centre of gravity and the pivot, *c* the viscous damping coefficient, *τ* the time, *g* the acceleration of gravity, *z*(*τ*) 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 *ω*₀ = √*g*/*l*, so Eq. (1) becomes

$$ \frac{d^2 x}{d\tau^2} + \frac{c}{ml^2} \frac{dx}{d\tau} + \left( \frac{1}{l} \frac{d^2 z}{d\tau^2} + \omega_0^2 \right) \sin(x) = 0. \quad (2) $$

Assuming $z(t) = -z \cos(\bar{\omega}t)$ we get

$$ \frac{d^2 x}{d\tau^2} + \frac{c}{ml^2} \frac{dx}{d\tau} + \left( \frac{z\bar{\omega}^2}{l} \cos(\bar{\omega}\tau) + \omega_0^2 \right) \sin(x) = 0 \quad (3) $$

while introducing the dimensionless time, damping coefficient, frequency and amplitude of the parametric excitation,

$$ t = \omega_0 \tau, \quad h = \frac{c}{\omega_0 ml^2}, \quad \omega = \frac{\bar{\omega}}{\omega_0}, \quad p = \frac{z\bar{\omega}^2}{l}, \quad (4) $$

respectively, gives the following dimensionless equation of motion [10]:

$$ \ddot{x} + h\dot{x} + [1 + p \cos(\omega t)] \sin(x) = 0, \quad (5) $$

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 *ω* very large (it is asymptotically exact for *ω* → ∞, 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

$$ \frac{d}{dt} \left[ \frac{\dot{x}^2}{2} - \cos(x) \right] = -h \dot{x}^2 - p \dot{x} \cos(\omega t) \sin(x). \quad (6) $$

For a given (initial) time $t_i$, let the (initial) position and velocity of the pendulum be $x_i = x(t_i)$ and $\dot{x}_i = \dot{x}(t_i)$, so that Eq. (6) can be expressed in the form

$$ \dot{x}^2(t) - 2 \cos[x(t)] - \dot{x}_i^2 + 2 \cos(x_i) = -2h \int_{t_i}^{t} \dot{x}(r) \frac{dx(r)}{dr} dr - 2p \int_{t_i}^{t} \cos(\omega r) \sin[x(r)] \frac{dx(r)}{dr} dr. \quad (7) $$

For rotating orbits, $\dot{x}(t)>0$, ∀*t*, so that the inverse function *t* = *t*(*x*) of *x*(*t*) exists. It is then possible to define unambiguously the function $\dot{x} = \dot{x}(t) = \dot{x}[t(x)] = f(x)$, which is central in the following development. We have

$$ t = t_i + \int_{x_i}^{x} \frac{ds}{f(s)}, \quad (8) $$

which explicitly provides *t* = *t*(*x*) when *f*(*x*) is known. For later use, we further note that $\dot{x}_i = f(x_i)$.

The previous considerations permit rewriting the differential Eq. (7) as the following integral equation:

$$ f(x)^2 - 2 \cos(x) - \dot{x}_i^2 + 2 \cos(x_i) = -2h \int_{x_i}^{x} f(s) ds - 2p \int_{x_i}^{x} \cos \left( \omega t_i + \omega \int_{x_i}^{z} \frac{ds}{f(s)} \right) \sin(z) dz, \quad (9) $$

where the unknown is *f*(*x*). Apart from the physical characteristics *h*, *p* and *ω*, the solution *f*(*x*) depends on the initial conditions $x_i$, $\dot{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,

$$ x_T = x(t_i + T) = x_i + 2\pi, \quad \dot{x}_T = \dot{x}(t_i + T) = \dot{x}_i, \quad (10) $$
where $T$ is the temporal period. From Eq. (8) we obtain the “synchronicity” equation

$$T = \int_{x_i}^{x_i+2\pi} \frac{dx}{f(x)}, \qquad (11)$$

while from Eq. (9) we obtain the “periodicity” equation

$$h \int_{x_i}^{x_i+2\pi} f(x) dx + p \int_{x_i}^{x_i+2\pi} \cos\left(\omega t_i + \omega \int_{x_i}^{x} \frac{ds}{f(s)}\right) \sin(x) dx = 0. \qquad (12)$$

**Remark.** Eq. (12) can be rewritten in the form

$$ha + p \cos(\omega t_i)b - p \sin(\omega t_i)c = 0, \qquad (13)$$

$$a = \int_{x_i}^{x_i+2\pi} f(x) dx, \quad b = \int_{x_i}^{x_i+2\pi} \cos\left(\omega \int_{x_i}^{x} \frac{ds}{f(s)}\right) \sin(x) dx, \quad c = \int_{x_i}^{x_i+2\pi} \sin\left(\omega \int_{x_i}^{x} \frac{ds}{f(s)}\right) \sin(x) dx. \qquad (14)$$

Let the angle $\alpha$ be defined by

$$\sin(\alpha) = \frac{b}{\sqrt{b^2 + c^2}}, \quad \cos(\alpha) = \frac{c}{\sqrt{b^2 + c^2}}, \qquad (15)$$

so that Eq. (13) becomes

$$ha + p\sqrt{b^2 + c^2} \sin(\alpha - \omega t_i) = 0 \quad \text{or} \quad \sin(\omega t_i - \alpha) = \frac{h}{p} \frac{a}{\sqrt{b^2 + c^2}}. \qquad (16)$$

Consequently, a solution exists if and only if

$$p \ge p_{cr} = h \frac{a}{\sqrt{b^2 + c^2}}. \qquad (17)$$

$p_{cr}$ represents a lower bound for the existence of periodic solutions, and it is attained for $\omega t_i = \alpha + \pi/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 = \epsilon h_1$ and $p = \epsilon p_1$ ($\epsilon$ 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 = 2\pi/\omega$, i.e., for rotating solutions having the same period as the excitation. The adopted procedure can be easily extended to subharmonic motions.

### 3. 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 = 2\pi$ at $t = t_i + T$. The governing equation, the synchronicity and periodicity conditions become, respectively,

$$f(x)^2 - 2 \cos(x) - \dot{x}_i^2 + 2 = -2\epsilon h_1 \int_0^x f(s) ds - 2\epsilon p_1 \int_0^x \cos\left(\psi_i + \omega \int_0^z \frac{ds}{f(s)}\right) \sin(z) dz, \qquad (18)$$

$$\frac{2\pi}{\omega} = \int_{0}^{2\pi} \frac{dx}{f(x)}, \qquad (19)$$
$$h_1 \int_0^{2\pi} f(x)dx + p_1 \int_0^{2\pi} \cos\left(\psi_i + \omega \int_0^x \frac{ds}{f(s)}\right) \sin(x)dx = 0, \quad (20)$$

where $\psi_i = \omega t_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:

$$\psi_i = \psi_0 + \epsilon\psi_1 + \dots, \quad \dot{x}_i = \dot{x}_0 + \epsilon\dot{x}_1 + \dots,$$

$$f(x) = f_0(x; \psi_0, \dot{x}_0) + \epsilon f_1(x; \psi_0, \psi_1, \dot{x}_0, \dot{x}_1) + \dots = f_0(x) + \epsilon f_1(x) + \dots \quad (21)$$

Inserting Eq. (21) in Eq. (18) and equating the coefficients of each $\epsilon$-power to zero we get

$$\epsilon^0: \quad f_0(x) = \sqrt{2 \cos x + \dot{x}_0^2 - 2}, \quad (22)$$

$$\epsilon^1: \quad f_0(x)f_1(x) = \dot{x}_0 \dot{x}_1 - h_1 \int_0^x f_0(s)ds - p_1 \int_0^x \cos\left(\psi_0 + \omega \int_0^z \frac{ds}{f_0(s)}\right) \sin(z)dz, \quad (23)$$

$$\epsilon^2: f_0(x)f_2(x) = \dot{x}_0 \dot{x}_2 + \frac{\dot{x}_1^2}{2} - \frac{f_1(x)^2}{2} - h_1 \int_0^x f_1(s)ds - p_1 \int_0^x \left(\psi_1 - \omega \int_0^z \frac{f_1(s)}{f_0(s)^2} ds\right) \sin\left(\psi_0 + \omega \int_0^z \frac{ds}{f_0(s)}\right) \sin(z)dz. \quad (24)$$

$f_0(x)$ is the classical solution for the unperturbed pendulum. The unknowns $\dot{x}_0, \dot{x}_1, \dots, \psi_0, \psi_1, \dots$, can be determined by the synchronicity

$$\epsilon^0: \quad \int_0^{2\pi} \frac{dx}{f_0(x)} = \frac{2\pi}{\omega}, \quad (25)$$

$$\epsilon^1: \quad \int_0^{2\pi} \frac{f_1(x)}{f_0(x)^2} dx = 0, \quad (26)$$

$$\epsilon^2: \quad \int_0^{2\pi} \frac{f_2(x)}{f_0(x)^2} dx = \int_0^{2\pi} \frac{f_1(x)^2}{f_0(x)^3} dx \quad (27)$$

and periodicity conditions

$$\epsilon^1: h_1 \int_0^{2\pi} f_0(x)dx + p_1 \int_0^{2\pi} \cos\left(\psi_0 + \omega \int_0^x \frac{ds}{f_0(s)}\right) \sin(x)dx = 0, \quad (28)$$

$$\epsilon^2: h_1 \int_0^{2\pi} f_1(x)dx + p_1 \int_0^{2\pi} \left(\psi_1 - \omega \int_0^x \frac{f_1(s)}{f_0(s)^2} ds\right) \sin\left(\psi_0 + \omega \int_0^x \frac{ds}{f_0(s)}\right) \sin(x)dx = 0. \quad (29)$$

It is straightforward to prove that Eqs. (28) and (29) guarantee periodicity of $f_1(x)$ and $f_2(x)$ respectively.

### 3.1. First-order terms

The solution is obtained in two steps. First, from the first synchronicity condition (25)

$$\frac{2\pi}{\omega} + \int_0^{2\pi} \frac{dx}{\sqrt{2 \cos(x) + \dot{x}_0^2 - 2}} = \frac{4}{\dot{x}_0} K\left(\frac{2}{\dot{x}_0}\right) \quad (30)$$

we obtain $\dot{x}_0 = \dot{x}_0(\omega)$, which is depicted in Fig. 1. $K(x) = \int_0^{\pi/2} [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 $\dot{x}_0 = \omega$ for $\omega \to \infty$ and satisfies the inequality $\dot{x}_0(\omega) > 2, \forall \omega$, which guarantees that the solution is rotating [5].
<figure><img src="image_2.png" /><figcaption>Fig. 1. Implicit function, $\dot{x}_0(\omega)$, determined from Eq. (30).</figcaption></figure>

In the second step, we determine the phase (or initial time) $\psi_0$ by the first periodicity condition (28), which can be rewritten as

$$0 = h_1 \int_0^{2\pi} f_0(x)dx + p_1 \cos(\psi_0) \int_0^{2\pi} \cos\left(\omega \int_0^x \frac{ds}{f_0(s)}\right) \sin(x)dx - p_1 \sin(\psi_0) \int_0^{2\pi} \sin\left(\omega \int_0^x \frac{ds}{f_0(s)}\right) \sin(x)dx \quad (31)$$

and is also the first-order expansion of Eq. (13).

As can be seen from Appendix A, the following relationships are satisfied:

$$\int_{0}^{2\pi} \cos\left(\omega \int_{0}^{x} \frac{ds}{f_{0}(s)}\right) \sin(x)dx = 0, \quad \int_{0}^{2\pi} \sin\left(\omega \int_{0}^{x} \frac{ds}{f_{0}(s)}\right) \sin(x)dx = 2, \quad \int_{0}^{\pi} \sin\left(\omega \int_{0}^{x} \frac{ds}{f_{0}(s)}\right) \sin(x)dx. \quad (32)$$

On the basis of Eq. (32), Eq. (31) becomes

$$h_1 \int_0^\pi f_0(x)dx = p_1 \sin(\psi_0) \int_0^\pi \sin\left(\omega \int_0^x \frac{ds}{f_0(s)}\right) \sin(x)dx \quad (33)$$

and $\psi_0 = \psi_0(\omega, h_1/p_1)$ can be determined by

$$\sin(\psi_0) = \frac{h_1}{p_1} \frac{\int_0^\pi f_0(x)dx}{\int_0^\pi \sin(\omega \int_0^x (ds/f_0(s))) \sin(x)dx} = \frac{h_1}{p_1} \gamma(\omega). \quad (34)$$

From Eq. (34) we see that there are two solutions for

$$p_1 > h_1 \gamma(\omega) = p_{1,SN}(\omega) \quad (35)$$

and no solutions for $p_1 < p_{1,SN}(\omega)$, so that $p_{1,SN}(\omega)$ 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(\psi_0) = 1$ and the phase is $\psi_0 = \pi/2$, while above it one solution has $\psi_0 < \pi/2$ and the other $\psi_0 > \pi/2$. Naturally, $p_{1,SN}$ is the first-order approximation of $p_{cr}$ as given in Eq. (17). The function $\gamma(\omega)$ is depicted in Fig. 2, along with the following asymptotic limits:

$$\gamma(\omega) \approx 2\omega, \quad \text{for} \quad \omega \to \infty; \quad \gamma(\omega) \approx 2/\omega, \quad \text{for} \quad \omega \to 0. \quad (36)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Function $\gamma(\omega)$, its Padé approximation $\gamma^*(\omega)$ and its asymptotic limits.</figcaption></figure>

The minimum of $\gamma(\omega)$ is attained approximately at $\omega = 1.06$ and is equal to 3.03821. The following Padé approximation of $\gamma(\omega)$ shown in Fig. 2 is robust:

$$ \gamma^*(\omega) = 2 \frac{1 + \omega^2 + \omega^4}{\omega + \omega^3}, \qquad (37) $$

so that it can be very useful in practical applications [27].

The function $\gamma(\omega)$ 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 = 0$, $\dot{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 rotations with $T = 2\pi n/\omega$, $n > 1$. The extension of the present approach to these solutions is straightforward and is not reported here.

## 3.2. Higher-order terms

Higher-order terms $\dot{x}_k$ and $\psi_k$, $k > 0$, can be obtained by the same steps as the first-order ones. First, we determine the velocity $\dot{x}_k$ from the synchronicity condition of order $k$, then the phase $\psi_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 $\dot{x}_1 = \dot{x}_1(\omega, h_1, p_1)$,

$$ \dot{x}_1 = h_1 \frac{\int_0^{2\pi} \frac{1}{f_0(x)^3} \left[ \int_0^x f_0(s) ds \right] dx}{\dot{x}_0 \int_0^{2\pi} (dx/f_0(x))^3} + p_1 \frac{\int_0^{2\pi} \left( \frac{1}{f_0(x)^3} \left[ \int_0^x \cos(\psi_0 + \omega \int_0^z (ds/f_0(s))) \sin(z) dz \right] dx \right)}{\dot{x}_0 \int_0^{2\pi} dx/f_0(x)^3}. \quad (38) $$
With this value of $\dot{x}_1$ it is possible to compute $f_1(x)$ by Eq. (23), so that from Eq. (29) we get $\psi_1 = \psi_1(\omega, h_1, p_1)$:

$$ \psi_1 = \frac{h_1}{p_1} \frac{\int_0^{2\pi} f_1(x) dx}{\int_0^{2\pi} \sin(\psi_0 + \omega \int_0^x (ds/f_0(s))) \sin(x) dx} + \omega \frac{\int_0^{2\pi} \left( \int_0^x (f_1(s)/f_0(s)^2) ds \right) \sin(\psi_0 + \omega \int_0^x (ds/f_0(s))) \sin(x) dx}{\int_0^{2\pi} \sin(\psi_0 + \omega \int_0^x (ds/f_0(s))) \sin(x) dx}. \quad (39) $$

## 4. Stability

To study the stability of the periodic solution we use the Poincaré section $\Sigma = \{(x, \dot{x}, t) \in S \times IR \times S$ such that $x = 0\}$. In fact, we take advantage of the space $2\pi$-periodicity of the pendulum which permits to identify the sections at $x = 2\pi n, \forall n$.

Let $(\dot{x}_i, t_i)$ be a point of $\Sigma$ and let $f(x; \dot{x}_i, t_i)$ be the solution of the integral Eq. (18) corresponding to $(\dot{x}_i, t_i)$. The image $(\dot{x}_{2\pi}, t_{2\pi})$ of this point on $\Sigma$ by the 2D Poincaré map $\Sigma \to \Sigma$ is given by

$$ \dot{x}_{2\pi} = f(2\pi; \dot{x}_i, t_i), \quad t_{2\pi} = t_i + \int_0^{2\pi} \frac{dx}{f(x)} - T. \quad (40) $$

The periodic solution corresponds to the fixed point of the map (40). In particular, from $\dot{x}_{2\pi} = \dot{x}_i$ we recover the periodicity condition, while from $t_{2\pi} = t_i$ we get the synchronicity condition. To study its stability we have to calculate the Jacobian

$$ \mathbf{J} = \begin{bmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{bmatrix} = \begin{bmatrix} \frac{\partial \dot{x}_{2\pi}}{\partial \dot{x}_i} & \frac{\partial \dot{x}_{2\pi}}{\partial t_i} \\ \frac{\partial t_{2\pi}}{\partial \dot{x}_i} & \frac{\partial t_{2\pi}}{\partial t_i} \end{bmatrix} \quad (41) $$

of the map and to compute the eigenvalues $\lambda_{1,2}$, which are the solutions of

$$ \lambda^2 - \lambda \operatorname{tr}(\mathbf{J}) + \det(\mathbf{J}) = 0. \quad (42) $$

The periodic solution is asymptotically stable as long as both $|\lambda_1| < 1$ and $|\lambda_2| < 1$. 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 $\mathbf{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

$$ \mathbf{J} = \begin{bmatrix} 1 - \varepsilon J_{11}^1 & \varepsilon J_{12}^1 \\ -J_{21}^0 + \varepsilon J_{21}^1 & 1 - \varepsilon J_{22}^1 \end{bmatrix} + \dots, \quad (43) $$

$$ \operatorname{tr}(\mathbf{J}) = 2 - \varepsilon(J_{11}^1 + J_{22}^1) + \dots = 2 - 2\varepsilon h_1 \left\{ \frac{\pi}{\omega} + \frac{2\omega}{\tan(\psi_0)} \int_0^\pi \frac{dx}{f_0(x)^3} \int_0^\pi f_0(x) dx \right\} + \dots, \quad (44) $$

$$ \det(\mathbf{J}) = 1 - \varepsilon(J_{11}^1 + J_{22}^1 - J_{21}^0 J_{21}^1) + \dots = 1 - 2\varepsilon h_1 \frac{\pi}{\omega} + \dots. \quad (45) $$

Eq. (44) can be further simplified by the relations:

$$ \int_0^\pi f_0(x) dx = 2 \dot{x}_0 E \left( \frac{2}{\dot{x}_0} \right). \quad (46) $$

$$ \int_0^\pi \frac{dx}{f_0(x)^3} = \frac{2}{\dot{x}_0(\dot{x}_0^2 - 4)} E\left(\frac{2}{\dot{x}_0}\right), \quad (47) $$

where $E(x) = \int_0^{\pi/2} [1 - x^2 \sin^2(s)]^{1/2} ds$ is the complete elliptic integral of the second kind [26].
The fact that $\det(J) < 1$ 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 $\epsilon p_{1,SN}$ and from above by $\epsilon p_{1,PD}$. Note that $\epsilon p_{1,SN} = 0.1\gamma(2) = 0.4168$ is just the value given by Eq. (35), while $\epsilon p_{1,PD} = \epsilon p_{1,SN}\delta(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 $\omega$ and $\epsilon h_1$. Before proceeding further, we note that in the unperturbed case $\epsilon = 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 $\lambda = 1$, one eigenvalue remains equal to 1 while the other becomes smaller than 1, and this can be easily achieved for small $\epsilon$. However, for the PD bifurcation, the eigenvalue has to reach, for growing $\epsilon$, the value -1. Since this is “far” from the starting point (for $\epsilon = 0$) $\lambda = 1$, it is expected that “large” values of $\epsilon$ would be required to achieve this bifurcation point. But in this $\epsilon$-range the perturbation analysis may fail, if we are out of the radius of convergence of the $\epsilon$-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 $\epsilon$-terms are likely required, but this is out of the scope of this work.

The SN bifurcation occurs when $\lambda = 1$. Substituting this value in Eq. (42) we obtain $\det(J) = 1 + \det(J)$. Using the expressions (44) and (45) we see that the SN occurs when $0 = \epsilon J_{21}^0 J_{12}^1 + \dots$. Since $J_{21}^0 > 0$ (see Appendix B), this is possible if and only if $J_{12}^1 = 0$, namely, by Eq. (B.9), if and only if $\cos(\psi_0) = 0$, which gives $\psi_0 = \pi/2$. This is exactly the same result as that obtained in Section 3.1 by a different reasoning, and thus it confirms the validity of the analysis.

The PD bifurcation occurs when $\lambda = -1$. Substituting this value in Eq. (42) we obtain $1 + \det(J) = 0$. Using the expressions (44) and (45) we see that PD occurs when $0 = 4 + \epsilon[-2(J_{11}^1 + J_{22}^1) + J_{21}^0 J_{12}^1] + \dots$,

<figure><img src="image_4.png" /><figcaption>Fig. 3. The eigenvalues for $\omega = 2$, $\epsilon h_1 = 0.1$ and for growing value of $\epsilon p_1$, $\lambda_{1,2}$ are the eigenvalues of the stable solution, $\mu_{1,2}$ the eigenvalues of the unstable solution: (a) large view and (b) zoom around the SN bifurcation.</figcaption></figure>
namely,

$$1 - \varepsilon h_1 \left\{ \frac{\pi}{\omega} + \frac{\omega}{\tan(\psi_0)} \int_0^\pi \frac{dx}{f_0(x)^3} \int_0^\pi f_0(x) dx \right\} + \dots = 0. \quad (48)$$

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 = \varepsilon J_{21}^0 J_{12}^1 + \dots$) has no $\varepsilon^0$ term, and therefore can be satisfied by small $\varepsilon$, while the equation for the PD bifurcation ($0 = 4 + \varepsilon[-2(J_{11}^1 + J_{22}^1) + J_{21}^0 J_{12}^1] + \dots$) has $\varepsilon^0$ term, and requires “large” values of $\varepsilon$ to be satisfied, which contradicts the philosophy of the perturbation analysis. Higher-order $\varepsilon$-terms have to be considered to improve the predictions on PD.

To determine the curve $p_{\text{PD}} = p_{\text{PD}}(\omega, h)$, we compute $\tan(\psi_0)$ from Eq. (48),

$$\tan(\psi_0) = \frac{\omega}{((1/\varepsilon h_1) - (\pi/\omega))} \int_0^\pi \frac{dx}{f_0(x)^3} \int_0^\pi f_0(x) dx, \quad (49)$$

so that

$$\frac{1}{\sin(\psi_0)} = \sqrt{1 + \frac{1}{\tan(\psi_0)^2}} = \sqrt{1 + \left[ \frac{(1/\varepsilon h_1) - (\pi/\omega)}{\omega \int_0^\pi dx/f_0(x)^3 \int_0^\pi f_0(x) dx} \right]^2} = \delta(\omega, \varepsilon h_1). \quad (50)$$

Note that if $\omega > \varepsilon h_1 \pi$ the right-hand side of Eq. (49) is positive, i.e., $\psi_0 < \pi/2$. As $\varepsilon h_1 = h$ is usually small, this is the common case. Substituting Eq. (50) in Eq. (34) we get

$$\varepsilon p_{1,\text{PD}}(\omega, \varepsilon h_1) = \varepsilon p_{1,\text{SN}}(\omega) \delta(\omega, \varepsilon h_1), \quad (51)$$

which provides the first-order approximation of $p_{\text{PD}}(\omega, h)$. The function $\delta(\omega, \varepsilon h_1)$ representing the ratio $p_{1,\text{PD}}/p_{1,\text{SN}}$ between the (first-order) critical thresholds for PD and SN bifurcations is depicted in Fig. 4 for various values of $\varepsilon h_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 $\omega$, 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

<figure><img src="image_5.png" /><figcaption>Fig. 4. The function $\delta(\omega, \varepsilon h_1)$ for various values of $\varepsilon h_1$.</figcaption></figure>
Poincarè 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.

## 5. 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 $\omega = 2$, $h = \varepsilon h_1 = 0.1$ for two different values of $p = \varepsilon p_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), \dot{x}(t)$] after eliminating the transient.

A very good correspondence is observed between the analytically predicted function $f_0(x) + \varepsilon f_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 $\varepsilon p_1$.

As expected, the two cases shown in Figs. 5a and b confirm that with increasing amplitude of the excitation, the $\varepsilon$ 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 $\omega = 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

<figure><img src="image_6.png" /><figcaption>Fig. 5. The analytical and numerical solutions for $\omega = 2$ and $h = \varepsilon h_1 = 0.1$ and (a) $p = \varepsilon p_1 = 0.42$, (b) $p = \varepsilon p_1 = 1.33$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. The numerical bifurcation diagram for $\omega = 2$. (a) $h = 0.01$, (b) $h = 0.05$, (c) $h = 0.1$ and (d) $h = 0.3$.</figcaption></figure>

<table><caption>Table 1
Comparison of analytical and numerical thresholds for $\omega = 2$</caption><thead><tr><th>$h = \varepsilon h_1$</th><th>$\varepsilon p_{1,SN}$ (theor.)</th><th>$p_{SN}$ (numer.)</th><th>$\varepsilon p_{1,PD}$ (theor.)</th><th>$p_{PD}$ (numer.)</th></tr></thead><tbody><tr><td>0.01</td><td>0.042</td><td>0.042</td><td>0.738</td><td>1.264</td></tr><tr><td>0.05</td><td>0.208</td><td>0.208</td><td>0.721</td><td>1.282</td></tr><tr><td>0.10</td><td>0.417</td><td>0.418</td><td>0.756</td><td>1.347</td></tr><tr><td>0.30</td><td>1.250</td><td>1.250</td><td>1.311</td><td>1.925</td></tr></tbody></table>

different results. In fact, while the theoretical $\varepsilon p_{1,SN}$ and numerical $p_{SN}$ thresholds coincide for all damping values, even quite large, those for PD are very different.

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 $\varepsilon$-series should be used.
<figure><img src="image_8.png" /><figcaption>Fig. 7. The numerical ($p_{SN}$ and $p_{PD}$) and theoretical ($\epsilon p_{1,SN}$ and $\epsilon p_{1,PD}$) critical thresholds for $h = \epsilon h_1 = 0.05$.</figcaption></figure>

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}/\epsilon p_{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 $\epsilon p_{1,PD}$ is to $p_{PD}$.

In the last numerical simulation we fix the damping to the intermediate value $h = \epsilon h_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 $\epsilon p_{1,SN}(\omega)$ of $p_{SN}(\omega)$ is always very accurate, while the first-order theoretical prediction $\epsilon p_{1,PD}(\omega)$ of $p_{PD}(\omega)$ 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.

## 6. 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 $\dot{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.
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.

## Acknowledgements

The support of the British Council, Italian MIUR and Italian CRUI, and by the Royal Society of London is gratefully acknowledged.

## Appendix A

In this appendix we prove Eq. (32). We start by noting that $f_0(x)$ is $2\pi$-periodic and even, so that $\int_0^{2\pi} f_0(s)ds = 2 \int_0^{\pi} f_0(s)ds$. Let $l(x) = \cos[\omega \int_0^x ds/f_0(s)] \sin(x)$. We have

$$ 
\begin{aligned}
l(x + 2\pi) &= \cos\left[\omega \int_{0}^{x+2\pi} \frac{ds}{f_0(s)}\right] \sin(x + 2\pi) \\
&= \cos\left[\omega \int_{0}^{x} \frac{ds}{f_0(s)} + \omega \int_{x}^{x+2\pi} \frac{ds}{f_0(s)}\right] \sin(x) \\
&= \text{by the periodicity of } f_0(s) \\
&= \cos\left[\omega \int_{0}^{x} \frac{ds}{f_0(s)} + \omega \int_{0}^{2\pi} \frac{ds}{f_0(s)}\right] \sin(x) \\
&= \text{by the synchronicity} \\
&= \cos\left[\omega \int_{0}^{x} \frac{ds}{f_0(s)} + \omega \frac{2\pi}{\omega}\right] \sin(x) \\
&= \cos\left[\omega \int_{0}^{x} \frac{ds}{f_0(s)}\right] \sin(x) = l(x),
\end{aligned}
\quad (A.1) $$

namely, $l(x)$ is $2\pi$-periodic. Furthermore, we note that

$$ 
\begin{aligned}
f_0(s) \text{ is even in } s &\Rightarrow \int_0^x \frac{ds}{f_0(s)} \text{ is odd in } x &\Rightarrow \cos\left[\omega \int_0^x \frac{ds}{f_0(s)}\right] \text{ is even in } x \\
&\Rightarrow \cos\left[\omega \int_0^x \frac{ds}{f_0(s)}\right] \sin(x) = l(x) \text{ is odd in } x &\Rightarrow \int_0^{2\pi} l(x)dx = 0.
\end{aligned}
\quad (A.2) $$

Analogously, it is possible to see that $m(x) = \sin[\omega \int_0^x ds/f_0(s)] \sin(x)$ is $2\pi$-periodic and even, so that $\int_0^{2\pi} m(x)dx = 2 \int_0^{\pi} m(x)dx$. We have thus proved Eq. (32).
Appendix B

In this appendix all components of Jacobian J, Eq. (41), are evaluated:

$$
J_{11} = \frac{\partial \dot{x}_{2\pi}}{\partial \dot{x}_i} = \frac{\partial f(2\pi; \dot{x}_i, t_i)}{\partial \dot{x}_i}.
$$

By taking the derivative of the integral equation with respect to dot{x}i and by defining g(x) = ∂f(x; x̂i, tᵢ) / ∂x̂i
we obtain the linear integral equation

$$
f(x)g(x) = \dot{x}_i - \varepsilon h_1 \int_0^x g(s)ds - \varepsilon p_1 \omega \int_0^x \sin \left( \psi_i + \omega \int_0^z \frac{ds}{f(s)} \right) \left( \int_0^z \frac{g(s)}{f(s)^2} ds \right) \sin(z)dz, \quad (B.1)
$$

permitting to determine g(x), from which we get J₁₁ = g(2π). The exact solution is not known, but according
to the previous section we assume g(x) = g₀(x) + εg₁(x) + ... and use a perturbation approach:

$$
\epsilon^0: \quad g_0(x) = \frac{\dot{x}_0}{f_0(x)} \quad \Rightarrow \quad g_0(2\pi) = 1, \qquad (B.2)
$$

$$
\epsilon^1: \quad f_0(x)g_1(x) = -f_1(x)g_0(x) + \dot{x}_1 - h_1 \int_0^x g_0(s) \, ds - p_1 \omega \int_0^x \sin\left(\psi_0 + \omega \int_0^z \frac{ds}{f_0(s)}\right) \left(\int_0^z \frac{g_0(s)}{f_0(s)^2} \, ds\right) \sin(z) \, dz \\
\Rightarrow \quad g_1(2\pi) = -h_1 \int_0^{2\pi} \frac{dx}{f_0(x)} - p_1 \omega \int_0^{2\pi} \sin\left(\psi_0 + \omega \int_0^x \frac{ds}{f_0(s)}\right) \left(\int_0^x \frac{ds}{f_0(s)^3}\right) \sin(x) \, dx. \quad (B.3)
$$

If we define

$$
\begin{align*}
J_{11}^1 &= h_1 \int_0^{2\pi} \frac{dx}{f_0(x)} + p_1 \omega \int_0^{2\pi} \sin\left(\psi_0 + \omega \int_0^x \frac{ds}{f_0(s)}\right) \left(\int_0^x \frac{ds}{f_0(s)^3}\right) \sin(x) \, dx \\
&= h_1 \frac{2\pi}{\omega} + p_1 \omega \sin(\psi_0) \int_0^{2\pi} \cos\left(\omega \int_0^x \frac{ds}{f_0(s)}\right) \left(\int_0^x \frac{ds}{f_0(s)^3}\right) \sin(x) \, dx \\
&\quad + p_1 \omega \cos(\psi_0) \int_0^{2\pi} \sin\left(\omega \int_0^x \frac{ds}{f_0(s)}\right) \left(\int_0^x \frac{ds}{f_0(s)^3}\right) \sin(x) \, dx, \tag{B.4}
\end{align*}
$$

we finally get

$$
\begin{align}
J_{11} &= 1 - \varepsilon J_{11}^1 + \dots, \tag{B.5} \\
J_{12} &= \frac{\partial \dot{x}_{2\pi}}{\partial t_i} = \frac{\partial f(2\pi; \dot{x}_i, t_i)}{\partial t_i}. \nonumber
\end{align}
$$

By taking the derivative of the integral equation with respect to tᵢ and by defining m(x) = (∂f(x; x̂ᵢ, tᵢ) / ∂tᵢ)
we obtain (remember that ψᵢ = ωtᵢ)

$$
f(x)m(x) = -\varepsilon h_1 \int_0^x m(s) \, ds + \varepsilon p_1 \omega \int_0^x \sin\left(\psi_i + \omega \int_0^z \frac{ds}{f(s)}\right) \left(1 - \int_0^z \frac{m(s)}{f(s)^2} \, ds\right) \sin(z) \, dz. \quad (B.6)
$$

Eq. (B.6) is the linear integral equation permitting to determine *m*(*x*), from which we get *J*₁₂ = *m*(2*π*). Using
the perturbation technique we assume *m*(*x*) = *m*₀(*x*) + *εm*₁(*x*) + ... and we get

$$
\begin{align*}
\epsilon^0: \quad m_0(x) &= 0, \tag{B.7} \\
\epsilon^1: \quad f_0(x)m_1(x) &= p_1 \omega \int_0^x \sin\left(\psi_0 + \omega \int_0^z \frac{ds}{f_0(s)}\right) \sin(z) \, dz \\
\Rightarrow \quad m_1(2\pi) &= \frac{p_1 \omega}{\dot{x}_0} \int_0^{2\pi} \sin\left(\psi_0 + \omega \int_0^x \frac{ds}{f_0(s)}\right) \sin(x) \, dx. \tag{B.8}
\end{align*}
$$
If we define (use is made of (33))

$$J_{12}^{1} = \frac{p_{1}\omega}{\dot{x}_{0}} \int_{0}^{2\pi} \sin\left(\psi_{0} + \omega \int_{0}^{x} \frac{ds}{f_{0}(s)}\right) \sin(x) \, dx = \dots = \frac{h_{1}\omega}{\dot{x}_{0}} \frac{1}{\tan(\psi_{0})} \int_{0}^{2\pi} f_{0}(x) \, dx, \quad (B.9)$$

we finally get

$$J_{12} = \varepsilon J_{12}^{1} + \dots \qquad (B.10)$$

$$J_{21} = \frac{\partial t_{2\pi}}{\partial \dot{x}_i} = \frac{\partial}{\partial \dot{x}_i} \left[ t_i + \int_0^{2\pi} \frac{dx}{f(x)} - T \right] = - \int_0^{2\pi} \frac{g(x)}{f(x)^2} dx$$

By the perturbation approach we get

$$J_{21} = -\dot{x}_0 \int_0^{2\pi} \frac{dx}{f_0(x)^3} + \varepsilon \int_0^{2\pi} \left[ 2\dot{x}_0 \frac{f_1(x)}{f_0(x)^4} - \frac{g_1(x)}{f_0(x)^2} \right] dx + \dots = -J_{21}^0 + \varepsilon J_{21}^1 + \dots \quad (B.11)$$

Note that $J_{21}^0 > 0$.

$$J_{22} = \frac{\partial t_{2\pi}}{\partial t_i} = \frac{\partial}{\partial t_i} \left[ t_i + \int_0^{2\pi} \frac{dx}{f(x)} - T \right] = 1 - \int_0^{2\pi} \frac{m(x)}{f(x)^2} dx.$$

By the perturbation approach we get

$$J_{22} = 1 - \varepsilon \int_{0}^{2\pi} \frac{m_{1}(x)}{f_{0}(x)^{2}} dx = 1 - \varepsilon p_{1}\omega \int_{0}^{2\pi} \frac{1}{f_{0}(x)^{3}} \left[ \int_{0}^{x} \sin\left(\psi_{0} + \omega \int_{0}^{z} \frac{ds}{f_{0}(s)}\right) \sin(z) dz \right] dx + \dots = 1 - \varepsilon J_{22}^{1} + \dots \quad (B.12)$$

## References

[1] F.M.A. Salam, S.S. Sastry, Dynamics of the forced Josephson junction circuit: the regions of chaos, *IEEE Transaction on Circuits Systems* CAS32 (1985) 784–796.

[2] G. Grüner, A. Zawadowski, P.M. Chaikin, Nonlinear conductivity and noise due to charge-density-wave depinning in NbSe$_3$, *Physical Review Letters* 46 (1981) 511–516.

[3] I.I. Blekhman, *Vibrational Mechanics*, World Scientific, Singapore, 2000.

[4] J. Guckenheimer, P. Holmes, *Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields*, Springer, New York, 1983.

[5] P. Hagedorn, *Nonlinear Oscillations*, Oxford University Press, New York, 1988.

[6] A.H. Nayfeh, *Introduction to Perturbation Techniques*, Wiley–Interscience, New York, 1981.

[7] E. Butikov, The rigid pendulum—an antique but evergreen physical model, *European Journal of Physics* 20 (1999) 429–441.

[8] R. van Dooren, Chaos in a pendulum with forced horizontal support motion: a tutorial, *Chaos, Solitons & Fractals* 7 (1996) 77–90.

[9] W. Garira, S.R. Bishop, Rotating solutions of the parametrically excited pendulum, *Journal of Sound and Vibration* 263 (2003) 233–239.

[10] W. Szemplinska-Stupnicka, E. Tyrkiel, A. Zubrzycki, The global bifurcations that lead to transient tumbling chaos in a parametrically driven pendulum, *International Journal of Bifurcation and Chaos* 10 (9) (2000) 2161–2175.

[11] S.R. Bishop, M.J. Clifford, The use of manifold tangencies to predict orbits, bifurcations and estimate escape in driven systems, *Chaos, Solitons & Fractals* 7 (1996) 1537–1553.

[12] B.P. Koch, R.W. Leven, Subharmonic and homoclinic bifurcations in a parametrically forced pendulum, *Physica D* 16 (1985) 1–13.

[13] P.B. Bornemann, U. Galvanetto, M.A. Crisfield, Some remarks on the numerical time integration of non-linear dynamical systems, *Journal of Sound and Vibration* 252 (2002) 935–944.

[14] E.M. Abdel-Rahman, A.H. Nayfeh, Z.N. Masoud, Dynamics and control of cranes: a review, *Journal of Vibrations and Control* 9 (2003) 863–908.

[15] W. Chester, The forced oscillations of a simple pendulum, *Journal of the Institute of Maths Applications* 15 (1975) 289–306.

[16] R.A. Struble, On the subharmonic oscillations of a pendulum, *ASME Journal on Applied Mechanics* 30 (1963) 301–303.

[17] W. Szemplinska-Stupnicka, E. Tyrkiel, The oscillation-rotation attractors in the forced pendulum and their peculiar properties, *International Journal of Bifurcation and Chaos* 12 (2002) 159–168.

[18] J.A. Blackburn, Experimental study of chaos in driven pendulum, *Physica D* 26 (1987) 385–395.

[19] P.L. Kapitsa, Dynamic stability of the pendulum when the point of suspension is oscillating, *Journal of Experimental and Theoretical Physics* 21 (1951) 588–597 (in Russian).
[20] J. Starrett, R. Tagg, Control of a chaotic parametrically driven pendulum, *Physical Review Letters* 74 (1995) 1974–1977.

[21] R.W. Leven, B. Pompe, C. Wilke, B.P. Koch, Experiments on periodic and chaotic motions of a parametrically forced pendulum, *Physica D* 16 (1985) 371–384.

[22] X. Xu, E. Pavlovskaia, M. Wiercigroch, F. Romeo, S. Lenci, Dynamic interactions between parametric pendulum and electro-dynamical shaker, *ZAMM—Zeitschrift für Angewandte Mathematik und Mechanik* 87 (2007) 172–186.

[23] Capecchi, D., Bishop, S., 1990. Periodic and non-periodic responses of a parametrically excited pendulum. Report 3/90 of the Dipartimento di Ingegneria delle Strutture, delle Acque e del Terreno, Università dell’Aquila, Italy.

[24] M.J. Clifford, S.R. Bishop, Rotating periodic orbits of the parametrically excited pendulum, *Physics Letters A* 201 (1995) 191–196.

[25] X. Xu, M. Wiercigroch, M.P. Cartmell, Rotating orbits of a parametrically-excited pendulum, *Chaos, Solitons & Fractals* 23 (2005) 1537–1548.

[26] M. Abramowitz, I. Stegun, *Handbook of Mathematical Functions*, Dover Publications, New York, 1970.

[27] G.A. Baker Jr., The theory and application of the Pade approximant method, in: K.A. Brueckner (Ed.), *Advances in Theoretical Physics*, Vol. 1, New York, Academic Press, 1965, pp. 1–58.

[28] S. Lenci, G. Rega, Higher-order Melnikov functions for single-dof. oscillators: theoretical treatment and applications, *Mathematical Problems in Engineering* 2 (2) (2004) 145–168.

[29] X. Xu, M. Wiercigroch, Approximate analytical solutions for oscillatory and rotational motion of a parametric pendulum, *Nonlinear Dynamics* 47 (2007) 311–320.