Bifurcations of Oscillatory and Rotational Solutions of Double Pendulum with Parametric Vertical Excitation

Abstract

This paper investigates dynamics of double pendulum subjected to vertical parametric kinematic excitation. It includes detailed bifurcation diagrams in two-parameter space (excitation’s frequency and amplitude) for both oscillations and rotations in the domain of periodic solutions.

1. Introduction

The construction of parametrically excited double pendulum is quite simple, but its dynamics are very complex. There is still lack of complete bifurcation diagrams of oscillatory and rotational motion. The influence of the parametric forcing on the dynamical system was investigated by many researchers. The different directions of excitation for the double pendulum are investigated in authors’ thesis [1, 2].

The chaotic motion is observed, as indicated in [3] or in [4]. The full bifurcation diagrams for horizontally, elliptically, and vertically excited single pendulum are presented by Horton et al. [5].

The detailed description of the parametrically excited double pendulum, with a special emphasis on 2 : 2 interaction mode, can be found in [6]. This work was further expanded [7], where authors show general explanation of the mode interaction. Here, the theory was proved by practical experiment.

The dynamical behaviour of the double pendulum, in the vicinity of a number of compound critical points, and the conditions of stability for the steady state solutions, in relation to the parameters of the systems, were explored by [8], resulting in closed form, achieved using bifurcation analysis for periodic and quasiperiodic solutions.

The approach towards numerical analysis of chaos in double pendulum is presented in [9]. It reveals regular behaviour of the pendula at the zero energy limit, which is transformed into quasiregular and globally chaotic motion, along with the energy increase. The bifurcation diagrams are not ordinary, due to energy application as the parameter, instead of one of the parameters occurring in the equations of motion.

Furthermore, an orthogonal double pendulum with a pivot point subjected to the excitation is also described in [10]. The used method of averaging allows investigating the resonant responses of the system, yielding to bifurcation analysis of the steady state constant solutions. It indicates the coexistence of planar and nonplanar harmonic motion for relatively large base excitations. Taking into consideration certain regions of parameter space, this motion becomes unstable and may bifurcate to quasiperiodic motion with changing amplitude with slow and fast frequency. Mentioned quasiperiodic motion is prone to lose stability and lead to chaotic, amplitude-modulated motion. Another approach towards investigation of the excited double pendulum is shown in [11]. It documents the parametric resonance causing the stable downward vertical equilibrium and the stabilization of the unstable upward equilibrium state by subjecting the whole system to the high-frequency motion. The issue of excited pendulum was also elaborated by Belyakov in [12]. He presented planar rotational motion of the mathematical pendulum with an elliptical excitation analysis, based on the exact solutions for the circular trajectory and no gravity. Derivation of existence and stability of these solutions allows finding approximate solutions for high and low linear damping.

In article [13], Thomsen examines the nonlinear dynamics of the elastically restrained double pendulum with nonconservative follower-type loading and linear damping for the presence of chaotic motion. It was revealed that static equilibrium solutions, stable periodic motion, and initially large changes in amplitude due to a destabilizing effect of both linear and nonlinear forces are all the result of a small, local nonlinear perturbation. A global numerical analysis proved that theoretical assumptions concerning occurrence of the chaotic motion in these regions are correct. Moreover, the combination of bifurcation cascade of static equilibrium points together with nonlinear destabilization may be treated as the cause of chaos trigger. The very similar system is presented in [14]; however, the investigation concerns destabilisation of the equilibrium position of the double pendulum with elastic end support and follower force loading. In Lyapunov criterion (for determining Lyapunov exponents from a time series, see [15]) this is a critical case, hence nonlinear analysis is performed, utilising centre manifold theory. The occurrence of heterochronic orbits suggested the possibility of divergent and periodic flutter motion; however, chaotic motion is also possible.

This paper is organized as follows. In Section 2, the considered model of double pendulum is presented. Section 3 is devoted to the numerical study of the bifurcation diagrams, performed in two-parameter space, namely, excitation frequency and amplitude, as well as in one-parameter space (i.e., amplitude) and basins of attractions. Finally, our results are summarised in Section 4.

2. Investigated System

In this paper, a system of double pendulum excited vertically is described. The system consists of planar double pendulum with nonelastic limbs. It is assumed that pendula are weightless with point masses 𝑚 on each end, respectively. The pivot point is vertically excited by parametric kinematic excitation. The excitation has amplitude 𝐴; excitation frequency is equal to 𝜔; 𝑙 denotes the length of each pendulum; 𝑐 is a damping coefficient on the pendula’ nodes.

In order to derive the equations of motion (4) second order Lagrange’s equations are used. For the system presented in Figure 1, the kinetic energy 𝑇 is formulated as follows: 1𝑇=2𝑚2𝐴²𝜔²cos²(𝜔𝑡)+𝑙2𝑙d𝜃₁d𝑡²+d𝜃₂()d𝑡−2𝐴𝜔cos𝜔𝑡sin𝜃₂+𝑙d𝜃₂+d𝑡d𝜃₁()d𝑡−4𝐴𝜔cos𝜔𝑡sin𝜃₁+2𝑙cos𝜃₁−𝜃₂d𝜃₂.d𝑡 (1) However, the potential energy 𝑉 is given by 𝑉=−𝑚𝑙𝑔−3+2cos𝜃₁+cos𝜃₂. (2) Finally, Rayleigh’s dissipation function 𝐷 takes the following form: 1𝐷=2𝑐d𝜃₁d𝑡²+d𝜃₁−d𝑡d𝜃₂d𝑡². (3)

Figure 1 

Scheme of the investigated system.

Putting (1), (2), and (3) into second order Lagrange equations, after some mathematical operations, one obtains equations of motion 𝑚𝑙²2d²𝜃₁d𝑡²+d²𝜃₂d𝑡²cos𝜃₁−𝜃₂+2𝑚𝑙sin𝜃₁𝑔+𝐴𝜔²sin𝜔𝑡+𝑐2d𝜃₁−d𝑡d𝜃₂d𝑡+𝑚𝑙²d𝜃₂d𝑡²sin𝜃₁−𝜃₂=0,𝑚𝑙²d²𝜃₁d𝑡²cos𝜃₁−𝜃₂+d²𝜃₂d𝑡²+𝑚𝑙sin𝜃₂𝑔+𝐴𝜔²sin𝜔𝑡+𝑐d𝜃₂−d𝑡d𝜃₁d𝑡−𝑚𝑙²d𝜃₁d𝑡²sin𝜃₁−𝜃₂=0. (4)

In the next part of this paper, dimensionless parameters are used: 𝜔²₀=𝑔/𝑙, 𝛾=𝜔/𝜔₀, 𝜌=𝑐/𝜔₀𝑚𝑙², 𝜏=𝜔₀𝑡, and 𝛽=𝐴/𝑙; symbols ̈‌ and ̇‌ denote d²/d𝜏² and d/d𝜏, respectively. These parameters yield to dimensionless equations as follows: 2𝜃̈₁+𝜃̈₂cos𝜃₁−𝜃₂+2sin𝜃₁1+𝛽𝛾²sin𝛾𝜏+𝜌2𝜃̇₁−𝜃̇₂+𝜃̇²₂sin𝜃₁−𝜃₂=0,𝜃̈₁cos𝜃₁−𝜃₂+𝜃̈₂+sin𝜃₂1+𝛽𝛾²sin𝛾𝜏+𝜌𝜃̇₂−𝜃̇₁−𝜃̇²₁sin𝜃₁−𝜃₂=0. (5)(6)

In the numerical calculations, we use the following value of dimensionless parameter 𝜌=0.04.

In order to calculate the natural frequencies of the system, equations (6) are linearised regardless of damping and excitation components. Consider 2𝜃̈₁+𝜃̈₂+2𝜃₁𝜃̈=0,₁+𝜃̈₂+𝜃₂=0. (7)

The natural frequencies calculated from (7) are 𝜔₁=2−√2, 𝜔₂=2+√2, which are 𝜔₁=0.765, 𝜔₂=1.848.

3. Results

In this section, results of the numerical analysis are presented, which include the bifurcation diagrams in two-parameter space (excitation frequency 𝛾 versus amplitude of excitation 𝛽), bifurcation diagrams in one-parameter space, and basins of attractions. The tool for numerical analysis is Auto 07P [16] and authors’ own programme using well-known 4th order Runge-Kutta method, with adaptive step size, from GNU GSL library version 1.15 [17]. The stability of periodic solutions is checked using Floquet multipliers [18]. In order to model physical behaviour of the system, the chosen range of parameters is as follows: 𝛾∈[0,5], 𝛽∈[0,3].

The system poses three main characteristic regimes of motion. The first one is when the excitation (parametrised by frequency 𝛾 and amplitude 𝛽), which is in vertical direction, is not big enough to start pendula motion. In such a case pendula are in stable fixed point 𝜃₁=𝜃₂=0. The other fixed points, (𝜃₁=0, 𝜃₂=𝜋), (𝜃₁=𝜋, 𝜃₂=0), and (𝜃₁=𝜃₂=𝜋), are unstable.

Provided that sufficient energy is delivered to the system via excitation, it is possible to observe oscillatory motion. Results for this type of motion are described in detail in Section 3.1.

Finally, providing higher energy to the system (especially for high frequency 𝛾), it is possible to obtain rotations. If this is the case, after transient time, the pendula rotate together (i.e., 𝜃₁=𝜃₂), just as in the oscillatory motion, as one pendulum. From the point of view of possible application of double pendulum as electric power generator, these types of solutions are particularly important for an engineer. Results for rotational motion are described in detail in Section 3.2.

This section is divided into three subsections: Section 3.1 oscillatory solutions, Section 3.2 rotations, and Section 3.3 basins of attractions. In Section 3.1, a case where both pendula oscillate is presented, which is typical for smaller values of excitation frequency 𝛾. In a case, where more energy is delivered to the system (i.e., for higher values of excitation frequency 𝛾 or excitation amplitude 𝛽), it is easier to set pendula in rotational motion, which is the topic of Section 3.2 rotations. Finally, the influence of the initial conditions of the system on the results is illustrated by means of basins of attractions in Section 3.3.

3.1. Oscillatory Solutions

In Figure 2(a), solutions for oscillatory motion of pendula are presented. The oscillatory periodic solutions are found usually for smaller energy delivered to the system (i.e., for smaller values of excitation frequency 𝛾 or excitation amplitude 𝛽). The system has enough energy to start its motion, however, not enough for the pendula to start rotation. In the investigated cases the pendula oscillate synchronised in phase (the lower pendulum follows the upper one). Numerical calculations find basically two locking ratios for oscillatory motion 1 : 1 : 1 and 2 : 1 : 1 (i.e., excitation period : upper pendulum period : lower pendulum period). Figure 2(a) focuses on 1 : 1 : 1 resonances and presents the borders of different solutions in two-dimensional space (𝛾, 𝛽). The black solid curve represents the beginning of motion when pendulum leaves the equilibrium state at stable fixed point (i.e., 𝜃₁=𝜃₂=𝜃̇₁=𝜃̇₂=0). It is worth mentioning that the tip of this curve corresponds to the value of natural frequency 𝜔₁=0.765 calculated in Section 2. The pendula oscillate synchronised in phase. Figure 2(b) shows the bifurcation scenario in one-parameter space (𝛽, max(𝜃₁)). The value of excitation frequency is set to 𝛾=0.68. The pendula starts motion for 𝛽=1.60, which corresponds to black curve from Figure 2(a). This can be mapped in Figure 2(b) as an intersection with 𝛽-axis. The grey dashed line presents unstable orbits, while the solid black one presents stable orbits. Three stable regions are encountered, which start in saddle-node (SN) and end in period doubling (PD) bifurcations. Their equivalents in two-parameter continuation can be found in Figure 2(a), marked with respective colors (thick lines mark SN bifurcations and solid lines mark PD bifurcations).

(a)

(b)

(c)

(d)

Figure 2 

Bifurcation diagrams of oscillations: the double pendulum (Figure 1) in two-parameter space (𝛾,𝛽). (a) Diagram representing 1 : 1 : 1 resonance. The thick lines apart from the black one represent saddle-node bifurcation, while the thin lines represent period doubling. (b) One-parameter space (𝛽, max𝜃1) continuation for 𝛾=0.68: thick line represents stable solution, while dashed line represents unstable solution. (c) Diagram representing 2 : 1 : 1 resonance. (d) One-parameter space continuation for 𝛾=1.2, red, 𝛾=1.4, blue, and 𝛾=1.7, black. Abbreviations of bifurcations: SN: saddle-node, PD: period doubling, and SB: symmetry breaking.

Figure 2(c) presents borders of oscillatory periodic solutions for 2 : 1 : 1 resonance. Here, five curves are found showing different types of bifurcations. The purple one represents the beginning of motion. Above this curve the stable fixed point (𝜃₁=𝜃̇₁=𝜃₂=𝜃̇₂=0) becomes unstable and pendula start to oscillate. For values of excitation frequency lower than 𝛾=1.535, the born solutions are unstable (subcritical bifurcations), while for higher frequencies solutions are stable (supercritical bifurcations). The black line corresponds to saddle-node bifurcations, where the unstable solutions coming from the purple curve can change stability properties. The next curve (blue) represents the symmetry breaking bifurcations; hence the motion of pendula becomes asymmetrical with respect to the equilibrium position. As it is well known [19–21], two mirror states are observed (one shifted in clockwise and second one shifted in counter-clockwise direction). The last one, brown line, also indicates a state when pendula leave the equilibrium and start period-2 oscillations. The red curve presents the beginning of motion.

To have a good understanding of what happened in the vicinity of resonance tongue (𝛾=1.535), we compute three one-parameter bifurcations diagrams, as shown in Figure 2(d). The parameters are 𝛾=1.2 (red line), 𝛾=1.4 (blue line), and 𝛾=1.7 (black line). Gray lines in Figure 2(d) correspond to positions of aforementioned one-parameter diagrams. As previously mentioned, the solid lines indicate stable solutions and dashed lines indicate unstable solutions. The red line (𝛾=1.2) for 𝛽=0.591 begins in the subcritical bifurcation, so the solutions along this curve are unstable. Hence, for 𝛽=0.043, we observe the saddle-node bifurcation, where curve turns but does not stabilize. Finally, for 𝛽=0.266, the symmetry of this unstable solution is broken. The blue line (𝛾=1.4) starts also in subcritical bifurcation; hence solutions are unstable. For 𝛽=0.022, we observe stabilization in the saddle-node bifurcation. Finally, for 𝛽=0.208 the symmetry breaking bifurcations occur. The last curve for 𝛾=1.7 (black line) shows different scenarios. The bifurcation where the motion originates from is supercritical, so from the beginning the solutions are stable. The change of stability properties occurs for 𝛽=0.42 in the symmetry breaking bifurcation. All the above described scenarios are different and can be treated as reference bifurcation sequences. Each of them corresponds to regions between cross sections of two-parameter bifurcation lines shown in Figure 2(c).

The curves obtained using the method of continuation are later compared with the results from numerical integration of the system equations of motion (5), (6). In Figure 3, the results are mapped with colored regions depending on the type of motion (oscillation, rotation) and order of periodicity. Additionally, curves from earlier figures are added in order to show the comparison between these 2 methods of system examining. The algorithm used in producing the color regions creates series of bifurcation diagrams along with excitation frequency 𝛾 for different values of excitation amplitude 𝛽. It analyses the obtained results, the order of periodicity, and character of motion (oscillations, rotations). In case of oscillations, it is obvious that the black curve from Figure 2, representing the beginning of motion, can be mapped by the green region of period-1 solution. The corresponding region for other curves presented in Figure 2 (saddle-node, period doubling pairs) can be found in Figure 3(b). It is especially clear that solid lines for period doubling borders match region calculated by numerical integration (green to red transition). For the brown curve from Figure 2(c), describing the 2 : 1 : 1 resonances, it is possible to find its equivalent in Figure 3(c). The red region indicating period-2 oscillation corresponds to the beginning of motion character of brown curve. The red curve from Figure 2(c) matches with period-2 region in Figure 3(c). This is a very clear tongue of resonance 2 : 1 : 1.

(a)

(b)

(c)

Figure 3 

Bifurcation diagrams of oscillations of the double pendulum (Figure 1) in two-parameter space (𝛾,𝛽) using 2 methods. (a) Curves from Figure 2(a) (1 : 1 : 1 resonance). (b) Detailed view in Figure 3(b). (c) Curves from Figure 2(c) (2 : 1 : 1 resonance). Initial conditions for numerical integration (𝜃1=𝜋, 𝜃̇1=−5, 𝜃2=𝜋−0.01, and 𝜃̇2=5). Abbreviations of bifurcations: SN: saddle-node, PD: period doubling, and SB: symmetry breaking.

3.2. Rotations

Figure 4 presents bifurcation curves for rotating motion of the pendulum. The rotating motion happens for bigger values of excitation frequency 𝛾. For most of the cases, both pendula are rotating. There exist several cases for rotating pendula. The 2 major ones are pendula rotating clockwise and pendula rotating counter-clockwise. The direction of rotation depends mainly on the initial conditions applied to the system. However, in periodic domain there exists another combination when the rotation of pendula is in the opposite direction. Borders of obtained solutions in two parameters space obtained by AUTO software are presented in Figure 4. The possible combinations (i.e., direction of oscillations) of rotations, computed using numerical integration, are presented in Figure 5. Let us focus on the borders of periodic solutions, using continuation method (Figure 4) and how they match the results from numerical integration. The cyan curve in Figure 4 represents 1 : 1 : 1 branching point. In that case, the results obtained using continuation and numerical integration converge particularly in the bottom region of the discussed curve. A significant number of period doublings (2 : 1 : 1 resonance) are found (3 cases). They are marked in Figures 4 and 5 as blue curves. The first one can be found in region 𝛾∈[2.5,5], 𝛽=0.75 (see also results obtained by numerical integration presented in Figure 5), its detailed analysis is shown in Figure 4(b). The next period doubling worth mentioning is located in the range 𝛾∈[2.5,5], 𝛽∈[1,2.5]. Once again the result from continuation method converges with the numerical integration. However, in this case only to some extent the 3rd 2 : 1 : 1 period doubling located left from the previously discussed curve has only partial coverage in numerical integration result.

(a)

(b)

(c)

Figure 4 

Bifurcation diagrams of rotations of the double pendulum. (a) Two-parameter bifurcation (excitation frequency 𝛾 versus excitation amplitude 𝛽) for rotational solutions. (b), (c) Detailed view of (a). Abbreviations of bifurcations: SN: saddle-node, PD: period doubling, BP: branching point, and SB: symmetry breaking.

Figure 5 

The rotation curves in two-parameter space (𝛾, 𝛽). Each color of background represents different period order and type of motion (i.e., oscillation, rotations). Initial conditions for numerical integration (𝜃1=𝜋, 𝜃̇1=−5, 𝜃2=𝜋−0.01, and 𝜃̇2=5). Abbreviations of bifurcations: SN: saddle-node, PD: period doubling, BP: branching point, and SB: symmetry breaking.

Two toruses are to be found in the range 𝛾∈[2.5,5], 𝛽=0.75, red (2 : 1 : 1) and green (4 : 1 : 1). As well as in previous cases, the numerical integration equivalent can be found for the above mentioned curves. Moving our attention to the left top corner of Figure 4(a) (detailed view, Figure 4(c)), one can observe saddle-node bifurcation 2 : 1 : 1 resonance and small period doubling region 8 : 4 : 4. These findings are confirmed also by numerical integration (see Figure 5).

3.3. Basins of Attraction

A phenomenon of coexistence of solutions is very often encountered in nonlinear dynamics. Therefore, in order to check if stable period solutions coexist, the system is checked using basins of attractions for different configuration of parameters. This section focuses on basins of attraction. The software used for that purpose is Dynamics 2 [22]. The chosen values of parameters are 𝛾=0.68, 𝛽=0.85 and 𝛾=0.75, 𝛽=1.5 which correspond to the regime of periodic oscillations. Because of the fact that the investigated systems have four-dimensional phase space, the obtained results are two-dimensional cross sections of four-dimensional phase space.

Figure 6(a) is in (𝜃₁, 𝜃₂) space, where 𝛾=0.68, 𝛽=0.85, 𝜃̇₁=0.5, 𝜃̇₂=0.5, 𝜃₁∈[−𝜋,𝜋], and 𝜃₂∈[−𝜋,𝜋]. One can observe 2 symmetric basins (marked as white and red, resp.). The blue area attracts to the stable solution (i.e., 𝜃₁=0, 𝜃̇₁=0, 𝜃₂=0, and 𝜃̇₂=0). The white and red basins attract the solution to periodic solution of period-1. In the case of Figure 6(b) the cross section through parameter space is done for different space variables. The basin is in (𝜃₂, 𝜃̇₂) space (position and velocity of the second pendulum), 𝛾=0.68, 𝛽=0.85, 𝜃₁=−𝜋/2, 𝜃̇₁=2,  𝜃₂∈[−𝜋,𝜋], and 𝜃̇₂∈[−5,5]. Here, once again one deals with symmetric basins. 3 regions are found: the cyan attracts to equilibrium (green dot), the purple basin attracts to period-1 oscillations, and white basin also attracts to period-1 oscillation. One can notice a central symmetry with (𝜃₂=0, 𝜃̇₂=0) as reflection point. Finally, the last basin presented in Figure 6(c) indicates 2 basins of attraction. The phase space is cut along (𝜃₂, 𝜃̇₂) plane, for 𝜃₁=−𝜋/2, 𝜃̇₁=2, 𝛾=0.75, and 𝛽=1.5. The detected basins of attractions attract to period-1 oscillation, which are central symmetrical with reflection point (0,0). The green dot attracts to white basin, while the blue one attracts to the purple basin.

(a)

(b)

(c)

Figure 6 

Basin of attractions in different planes and parameters: (a) plane (𝜃1,𝜃2), 𝛾=0.68, 𝛽=0.85, 𝜃̇1=0.5, and 𝜃̇2=0.5; (b) plane (𝜃2,𝜃̇2), 𝛾=0.68, 𝛽=0.85, 𝜃1=−𝜋/2, and 𝜃̇1=2; (c) plane (𝜃2,𝜃̇2), 𝛾=0.75, 𝛽=1.5, 𝜃1=−𝜋/2, and 𝜃̇1=2.

Looking at every obtained basin, one may draw a conclusion that the solution depends strongly on the initial conditions. While the basins for the investigated systems are mostly symmetrical and attract to the same type of pendula motion, there is a significant large basin that attracts to stable fixed point at pendula equilibrium.

4. Conclusions

This paper presents broad numerical analysis of the dynamics of double pendulum excited vertically. Bifurcation diagrams in two-parameter space (excitation frequency 𝛾 versus excitation amplitude 𝛽), as a continuation of periodic solutions, are shown. For oscillations of the pendula, two locking ratios are observed (i.e., 1 : 1 : 1 and 2 : 1 : 1). However, in case of rotational solutions, one obtains also different locking ratios, such as 4 : 1 : 1 or 8 : 4 : 4. The emerging of resonance zone happens about excitation frequency 𝛾=0.7. Additionally, an iterative numerical integration, in order to check the results from continuation, is performed. The result of this comparison is successful, because it is possible to observe similar regions using both methods. Finally, the basins of attractions give the answer on influence of the initial conditions on the solutions. Symmetrical basins are obtained for most of the cases with 2 periodic basins.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work has been supported by the Foundation for Polish Science, Team Programme, under Project TEAM/2010/5/5.