Dynamics of a parametrically excited double pendulum

PHYSICA ELSEVIER

Physica D 75 (1994) 541-558

Dynamics of a parametrically excited double pendulum A.C. S k e l d o n School of Mathematics, University of Bristol, Bristol, BS8 1TW, UK Received 27 October 1993; revised 28 February; accepted 23 March 1994 Communicated by A.V. Holden

Abstract The 2:2 mode interaction of a parametrically excited double pendulum is explored in the excitation frequency/excitation amplitude plane. To determine the bifurcation structure at small amplitudes of oscillation, the method of averaging combined with centre manifold reduction is used. The full equations are solved numerically to extend the bifurcation set to larger amplitudes of response. Numerical centre manifold reduction is employed to derive two maps, valid near two multiple bifurcation points which organise the dynamical phenomena of the mode interaction region. Iteration of these maps shows the existence of global bifurcations. These results are discussed in the light of numerical integrations which show that a whole range of interesting behaviour occurs including torus doubling, torus 'gluing' and chaos.

1. Introduction

An experimental study of a parametrically excited double pendulum was performed by Skeldon and Mullin [ 1 ]. The experiment was motivated by a Faraday crispation experiment [2], performed by Ciliberto and Gollub [3,4], in which a fluid in a cylindrical container is parametrically excited by moving the cylinder up and down sinusoidally. This system has been studied by a number of different authors using a variety of techniques [4-10]. The difficulties associated with studying this fluid dynamical problem are highlighted in [10] which includes a review of the different approaches. The fluid system and the pendulum system share several common features: for example, the primary instability is governed by Mathieu's equation, as shown for the fluid by Benjamin and Ursell [5]. The pendulum consisted of two rods, one suspended from the other, with the bottom rod hinged so that it could only swing in a plane orthogonal to the top rod (see Fig. 1 ). It was parametrically excited by moving the top pivot point up and down. The system has two normal modes of oscillation: the 0-mode oscillation (i.e. ~b = 0), where the whole pendulum moves as one about the top pivot point and the ~b-mode oscillation (0 = 0), where just the bottom rod moves about the hinge. For some values of the excitation frequency and excitation amplitude it was found that neither mode was excited: the whole system was simply translated up and down with the drive. In other parameter regimes, oscillations in either the pure 0 mode or the pure ~b mode were observed. 0167-2789/94/$07.00 ~) 1994 Elsevier Science B.V. All fights reserved SSDI 0167-2789(94)00099-C

542

A.C, Skeldon / Physica D 75 (1994) 541-558

A

I

Fig. 1. Geometry of the double pendulum.

However, in a region where the excitation frequency was close to twice the natural frequency of both modes, a range of complicated dynamical phenomena, including tori, torus doubling and chaos occurred. Numerical integrations showed that these same phenomena occur in the equations of motion. Furthermore, the use of numerical path-following methods showed the presence of two multiple bifurcation points in the mode interaction region. In this paper we expand on the results of Skeldon and Mullin and present an analytical and numerical study of the equations of motion. In Sections 2 and 3 we present a local analysis of the trivial solution through the use of the method of averaging [I 1 ]. We find that there exists a multiple bifurcation point, M22, where the trivial solution loses stability to both the 0 mode and the ~ modes of oscillation. The centre manifold theorem [ 12] is used at M22 to construct the local bifurcation set. This shows that, in the neighbourhood of the multiple bifurcation point, there can be a region of the control parameter plane where there is no stable steady-state solution of the averaged equations, and consequently, no stable periodic solution of the full equations. Within this region, there is either a large-amplitude solution which would not be captured by the local analysis or a non-periodic solution. In order to explore the solution structure of the mode-interaction region further, numerical path-following methods are used to construct a bifurcation set of the periodic solutions of the full equations. The results are discussed in Section 4. They show that a further multiple bifurcation point exists, DS~, when a pitchfork bifurcation interacts with a fold. In order to study both this second multiple bifurcation point and the mode interaction point in more detail, numerical centremanifold reduction is applied to the T-map of the original equations, where T is the period of the excitation. This produces two maps, one valid in the neighbourhood of M22 and the other valid in the neighbourhood of DSl. A study of these maps is presented in Section 5 and shows the existence of global bifurcations and thereby a mechanism for chaotic behaviour. In Section 6 the numerical and analytical results are compared with numerical integrations. The conclusions are presented in Section 7.

543

A.C. Skeldon / Physica D 75 (1994) 541-558 bp: bl"li ~.0

V

I1 i

075

0.5

0.25

0.0 0.5

,

,

i

i

0.75

t

i

i

[

"~.0

i

i

i

i

~

1.25 Excitation

,

,

i

i

,

1.5 frequenc~

~

,

,

i

,

1.75

,

,

i

i

2.0

,

r

,

i

i

,

,

i

i

2.25

I

2.5

~

Fig. 2. Computed bifurcation set showing the four largest regions of instability for the trivial solution of the linearised double pendulum equations. The heavy lines outline regions of instability to a 0-mode oscillation, whilst the lighter lines outline regions of instability to a ~b-mode oscillation.

2. The double pendulum The non-dimensionalised equations of motion are [ ( M + 3)L 2 + 3Lcos<~ + cos2~b]/~- (3L + 2cos~b) sin0~ + ill0 + [ ( M + 2)L + cos~b] (1 + Ag22cosg2z) sin0 = 0, + ½(3L + 2 cos~b)sin<~ 62 + fi2~ + cos0(1 + AI22 cosg2z)sin~b = 0,

(1)

where dots denote differentiation with respect to the nondimensional time r. In a typical experiment most of the damping is caused by friction at the pivot point; assuming a linear damping-law gives a good approximation. There are six control parameters: M, the mass ratio of the two rods; L, the length ratio; fil and ~2, the damping coefficients; A, the non-dimensional excitation amplitude and g2, the non-dimensional excitation frequency. We fix the first four of these control variables and consider the behaviour in the A, O-plane: the only parameters which can be varied easily in an experiment. The problem has two reflection symmetries, i.e. 0 --+ - 0 and ~b ---, -~b leave the equations invariant (Z2 • 2~2 symmetry). Characteristic of parametrically excited systems Eqs. (1) have a trivial solution, 0 = 6 = ~b = --- 0. This solution is not always stable: since the equations linearised about this solution are two uncoupled Mathieu equations, the instability regions form tongues directed down towards the excitation frequency axis. For illustration, the four largest tongues are shown in Fig. 2. Furthermore, if ~b = ~ = 0 then the equation of motion for a parametrically excited simple pendulum of natural frequency o~0 = ~ , where kl = ( ( M + 3)L 2 + 3L + 1) -1 and k2 = ( ( M + 2)L + 1), is obtained. Similarly if 0 = 6 = 0 then we have the equation of motion for a parametrically excited simple pendulum of natural frequency o9~ = 1. The parametrically excited simple pendulum itself has been extensively studied [ 14-18]. Period-doubling and chaos occur even in this problem, and hence each of the modes of the double pendulum individually exhibit these phenomena. Here we are only concerned with the mode interaction behaviour and our study

544

A.C. Skeldon / Physica D 75 (1994) 541-558

is restricted to a region of the control parameter plane close to M22. In Section 4 (and Fig. 5) it is shown that, for the parameter values we study, the complicated single pendulum behaviour occurs away from this region.

3. Local analysis of the trivial solution

3.1. Averaging of the equations The parametric excitation of the double pendulum makes the equations of motion non-autonomous which means, of course, that there are no steady-state solutions. In order to remove this underlying periodicity, the equations of motion are averaged and the bifurcation structure of the resulting, simpler, steady-state problem is examined. As we wish to study the 2:2 mode interaction region, the mass and length ratios are chosen so that the natural frequency of each mode is close to half the excitation frequency, that is, (klk2- (22/4) ---- ~ 1 and (1 - 0 2 / 4 ) = ff~'2, where e is a small parameter. The experimental pendulum is lightly damped and the mode interaction point M22 occurs at a low amplitude of the excitation. Furthermore, the swinging solutions have a small amplitude and we therefore scale the equations as follows:

A = cA',

~i = £ 6 [ , 0 = ~l/2yl,

0 : ¢~1/2y2,

q~ = ¢~1/2y3,

q~ : £1/2y 4.

The scaled equations are truncated at cubic order and the averaging theorem of Krylov and Bogliubov [ 1 1 ] is applied as described in Guckenheimer and Holmes [ 19 ] for forced systems. This results in the averaged equations, j3 = Ly + N(fi) + O ( ~ 2 ) ,

fi E [~4

(2)

and

-dl L =

N(~)

=

(bl + cl )

(bl - cl ) -dl

0 0

0 0

0

0

0

-d2

(b2 - c2)

0

(b2 + c2)

-dE

\

J

'

( - a l ( f i ~ + fi~)Y2- a2Y2Y3 - -2 - a3Y2fi~ - a4Ylfi3Y4 -2)- a ~ a2 - -2 al(fi~ + Y2 Yl 2YlY4 + a3YlY3 + aaY2Y3Y4 - a 5 ( Y ~ + Y~)Y4 -a6Y~Y4 -a7Y~Y4 + asYlY2Y3

as(Y~ + Y~)fi3a6fi~fi3 + a7fi~fi3 - asYlY2Y4 where

A.C. Skeldon / Physica D 75 (1994) 541-558

k~ = ( ( M + 3)L 2 + 3L + 1 k2 =

a5 = 8 ~ '

(

)

e 2k4 + klk3£-22 1612 '

a3 = 1 - ~

'(

a4 = - ~

6k4 - klk3 O2 2k4-klk3£2

2

)

bl = ½klkEA~2 Cl •

, k 3 = 3L + 2, k4 = klk2k3 - 1,

( M + 2)L + 1,

al = ~-~klk2, a2-

-1

545

£ - ~ ~)1,

a6 a7 -

3212

3k3~22 - 4

)'

320

°( k 3 ~ 2 + 4 )

as = 1 - ~

,bE = ½A~2, £ C2 -~ - ~ ~2 .

The full details of both this and the centre manifold reduction below are contained in Skeldon's thesis [20].

3.2. Stability of the trivial solution The averaged equations retain the symmetry of the full problem and therefore have the trivial solution ~ = 0. This is the only solution which has the full 7/2 G 7/2 symmetry of the equations. The eigenvalues of L show that this solution is stable for low amplitudes of the excitation frequency. However, when d ~ - ( b 2 - c 2) = 0 ,

(3)

there is a simple bifurcation to a 0-mode solution which breaks one of the 7/2 symmetries. Similarly, when d 2 - ( b 2 - c 2) = 0 ,

(4)

the trivial solution bifurcates to a ~b-mode solution, breaking the other 7/2 symmetry. In the full non-autonomous equations, Eqs. (3) and (4) describe curves of period-doubling bifurcations: at small amplitudes of the excitation they give good quantitative agreement with the numerical curves pdo and pd~ shown in Fig. 2. The mode interaction point, M22, occurs at the point of intersection of (3) and (4). The eigenvalues of L show that no Hopf bifurcation can occur from the trivial solution. This is also true for the damped, parametrically excited simple pendulum, where bifurcations from the trivial solution are sometimes mistakenly referred to as Hopf bifurcations. It can readily be shown that there are also no Hopf bifurcations from the pure mode solution branches.

3.3. Centre manifold reduction At a bifurcation point the centre manifold theorem [ 12] can be used to reduce a bifurcation problem to a simpler set of equations which locally encapsulates the bifurcation structure. The multiple bifurcation point M22 occurs at (12m, Am ) on the control parameter plane where ~'-2m and Am satisfy both Eqs. (3) and (4). At this point the averaged Eqs. (2) have two zero

546

A.C. Skeldon I Physica D 75 (1994) 541-558

eigenvalues and application of the centre manifold theorem results in two equations on the centre manifold. In turn, invoking the singularity theory of Golubistsky and Schaeffer [21 ], these centre manifold equations are Z2 G Z2 equivalent to the bifurcation problem, h ( u , v , 2 , Vn,'h,a) = [elu 3 + Fnuv 2 + e 2 ( / ~ - o-)u, nu2v + E3v3 + f 4 ~ ) ] ,

(5)

where =

m =

-

a2 r4~l

r2,

,

e3

=

sgn(-r4),

e4

=

sgn(-1),

n = -- ] ar--~ r3,

and

fl, =

k,k2

Am + ~

k,-~m +

'

almClm rl -dim '

1

r2 -

dl.,S2 m (blm - Clm) [a2m (blmCXm - b2~.Clm ) - a3m (blmC2m + b2,.clm) ] ,

1 b2m Qrn, a2 = 2 d2.,

[lbz,.Am c2., ( 1 ¼)] 1

r3 r4--

dl,,s2" (bl,, - c2,, ) [a6,, (bzmcl,, - bl,.Czm ) - aT,, (b2..Cl,. + bl..czm ) ] , d2.,

Here, ai,., bi,. and ci,. are the values which ai, bi and ci take at the mode interaction point (g2m,A m ). Furthermore, to first order, 2 and a are related to our control parameters by 2=fl2(Q-Qm)+a2(A-Am)'

° ' = ( f12

a 2( Q f l-l o) m ) ' o Q

Note that fixing a and varying 2 is qualitatively the same as fixing ~2 and varying A, and that u is related to the amplitude of the 0 mode and that v is related to the amplitude of the q~ mode. The general bifurcation problem (5) was considered by Golubitsky and Schaeffer. Note that, although the general form of this bifurcation problem arises from the Z2 G Z2 symmetry of the original problem, without the analysis we would have no expressions for the modal parameters and for the ei, and thus would not be able to give the particular unfolding for this case. For the numerical results contained in this paper, specific values of M and L and the damping coefficients, Oi were taken. There is a range of choices for which the scaling chosen in Section 3.1 is valid: we consider M = 7, L = 1.1 and Ol = 02 = 0.245. This gives 2oJ0 = 1.63 and 2o~ = 2 and that the mode interaction point occurs at (g-2m,Am) = (1.84,0.16). In turn, this gives el = + 1,¢2 = -1,¢3 = - 1 and e4 = - 1 , and the modal parameters as m = -10.52 and n = 9.78. This is within region VII of Golubitsky and Schaeffer's modal parameter plane classification and well away from the boundaries

547

A . C . Skeldon / Physica D 75 (1994) 5 4 1 - 5 5 8 0.17

~

.1675

"~/"

/'"

O.165

t" l

/,'"

1625

0.16.1575-

td

0.15 \\ ~.1525

\'~\

0.15 i J l , l , l l t l l , ~ j , , , , , i , , , , , , , , l l l ~ l 1.83 1,835 1.84

,,,~,,,i 1.845

1.85

Excitadon frequency, l~

Fig. 3. Bifurcation set for the neighbourhood of M22 resulting from the centre manifold analysis of the averaged equations. pdo, primary bifurcation curve to the 0-mode oscillation. pd~,, primary bifurcation curve to the ~b-mode oscillation. ..... se, secondary bifurcation curve from the 0-mode oscillation to the mixed mode. .............. s6, secondary bifurcation curve from the {,-mode oscillation to the mixed mode.

- _ _

of this region; for neighbouring values of M and L the qualitative bifurcation structure will be preserved. Analysis of (5) leads to the local bifurcation structure for the mode interaction point M22. There are three solution types: a trivial solution u = v = 0; the pure-mode solutions, u ~ 0, v = 0 and u = 0, v ~ 0; and a mixed-mode solution, u ~ 0, v ~ 0. In terms of the pendulum the pure u-mode solutions correspond to periodic swinging of the 0 mode and the pure v-mode solutions correspond to periodic swinging of the ~b mode. The mixed-mode solution corresponds to both the 0 and ~b mode swinging simultaneously and both with the same period. The curves of bifurcation points on the control parameter plane (the bifurcation set) are shown in Fig. 3. Figs. 4a,b and c are slices through this bifurcation set. These bifurcation diagrams show the amplitude of the 0-mode solution versus 12 in the plane of the paper, and the amplitude of the ~b-mode solution out of the plane of the page. The convention that stable solutions are shown with solid lines and unstable solutions with dotted lines is used. From Fig. 4a we see that for small values of the excitation frequency the trivial solution is stable. If the excitation amplitude is increased then this trivial solution becomes unstable at pdo where there is a supercritical bifurcation to a 0-mode solution: the whole pendulum starts to swing in the 0 direction. As A is increased further the amplitude of this solution increases. The ~b-mode solution bifurcates sub-critically from the trivial solution at pd~, but this solution is unstable at small amplitudes and so would not be observed in an experiment. As t2 is increased the two primary bifurcation points pdo and pd~ occur closer and closer together until when 12 = 12m they coincide (Fig. 4b). Increasing 12 further, as in Fig. 4c, interchanges the order of the two primary bifurcations so that the ~b mode bifurcates at the lower value of A. In order for this to occur, secondary bifurcations occur on the pure-mode solution branches at so and s~, producing a mixed-mode branch. At one end the mixed-mode branch has two positive eigenvalues and is therefore unstable, whilst at the other end it has two negative eigenvalues and is therefore stable. As Golubitsky and Schaeffer point out, this is consistent if somewhere along the mixed-mode branch a

548

A.C. Skeldon / Physica D 75 (1994) 541-558 /

\

(a) (b)

\

.4

Ic) Fig. 4. Bifurcation diagrams for (a) • < ~2m, (b) £2 = ~2m (c) £2 > ~2m. pdo and pde: mark the primary bifurcations to the 0 and ~ modes, respectively; M22 marks the 2:2 mode interaction point; so and s~ mark the secondary bifurcations from the pure 0 mode and pure ~ mode, respectively, to a mixed mode.

H o p f bifurcation occurs with its consequent double change in eigenvalues: a line of such bifurcations must emanate from the mode interaction point. In the full double pendulum equations this H o p f bifurcation corresponds to a transition from both the 0 mode and ~ mode swinging with the same period, to both swinging, but with slightly different periods. In Section 3.2 we stated that there can be no H o p f bifurcation from the trivial solution or from the pure-mode solution branches; there is no contradiction here as this H o p f bifurcation occurs on a mixed-mode solution branch. The H o p f bifurcation suggests the presence of more complicated dynamics in the full problem. In the averaged equations this H o p f bifurcation gives rise to a periodic orbit, whereas in the full problem it indicates the production of a torus. The original equations (1) were averaged to cubic order, but this is not sufficient to determine whether the H o p f bifurcation is super- or sub-critical or where it occurs along the mixed-mode branch. With the analysis of this section we have found the local periodic bifurcation structure of M22 and have shown that there exists a stable small amplitude mixed-mode solution. Furthermore this solution undergoes a H o p f bifurcation: if the locus of these bifurcations lies above the line pd~ in Fig. 3 (as we shall see in Section 5.1, this is the case for our problem), then there are no stable small amplitude periodic solutions. Instead there must be either a large amplitude periodic motion or some small amplitude non-periodic motion. In order to extend these local results to larger amplitude we turn to numerical methods.

A.C. Skeldon / Physica D 75 (1994) 541-558 o 4-

/,~

0 35-

.....

•.........

0.3~ i

549

".""

0.25- .. ~' .... . . . • ........ o,2-J o.15-

.

.

,,' ,.'/" ~

p4~

,.

"""

.,""'"

j/t~ / I s~

.s/" /

//

0.05 0.0 1.5

1.6

1.7

1.8

1.9

2.0

21

E×c~tat~onfrequency,

Fig. 5. Numerical bifurcation curves. - - p d e , primary bifurcation curve to the 0-mode oscillation. _ _ p d o , primary bifurcation curve to the O-mode oscillation. ..... se, sccondary bifurcation curve from the 0-mode oscillation to the mixcd mode. ..............s0, secondary bifurcation curve from the 0-mode oscillation to the mixed mode. ............ Is, limit points on the 0-mode branch. ...........I~, limit points on the 0-mode branch. ..... t0, symmetry breaking bifurcation on thc ~-mode branch.

4. Numerical bifurcation set

Numerical path-following methods were employed on the full nonlinear non-autonomous equations (1). As stated earlier, the periodic excitation means that all solutions are at least periodic. Numerically, it is therefore convenient to consider fixed points of the Poincard map constructed by integrating the equations through one cycle. Standard path-following ideas were applied to this map: an extended system was written for each type of bifurcation [22,23] and the Z2 ~)Z2 symmetry of the equations was used to reduce the dimension of the computational problem. Each extended system was solved for the relevant parameter values using Newton's method. The full details of the systems used are contained in Skeldon's thesis [20]. In Fig. 5 we show the lines which were calculated in the vicinity of the mode interaction point, M22. The lines pdo and pd¢~ are the curves of primary period-doubling bifurcations representing bifurcation from the trivial solution to the periodic 0-mode swinging solution and the periodic b-mode solution, respectively. These curves are the same as pdo and pd~ of Fig. 2. The mode interaction point M22 occurs at the intersection of these primary bifurcation curves, a t / 2 = 1.8472 and A = 0.1623. This agrees well with the values of 12 = 1.84 and A = 0.16 given by the analysis of Section 3.2. In good quantitative agreement with the centre manifold analysis, emerging from M22 are the secondary bifurcation lines so and s6. The two loci of bifurcation points lo and l~ are lines of limit points, the first on the 0-mode branch and the second on the ~b-mode branch. These join the primary bifurcation tongues at the point where there is a transition from a sub- to a supercritical primary bifurcation. The points of contact correspond to degenerate points on the primary bifurcation curves. Where the secondary bifurcation points s~ passes down and round the line of limit points l~ there is a second multiple

550

A.C. Skeldon / Physica D 75 (1994) 541-558

bifurcation point, DS1. This is a periodic double singular point. There is one further line shown on Fig. 5, t 0. This is a line of symmetry-breaking bifurcations which occurs on the pure q~-mode branch. In this case, the bifurcation does not break the Z2 • Z2 symmetry of the problem but breaks the temporal symmetry of the oscillation: the resulting solution is still a pure mode periodic q~-mode solution, but one in which the pendulum does not swing to the same height on both sides of its oscillation. This bifurcation line is the first in a sequence of bifurcations which leads to the single pendulum-like chaos involving only one of the modes as mentioned in Section 2. All the other bifurcations in the sequence lie above this line, as shown by Bryant and Miles I17]. Thus the mode interaction we consider occurs in a region of the parameter plane away from these complicated single-mode solutions. A similar line exists for symmetrybreaking of the pure 0 mode, but this line lies to the left of the region in parameter space shown in Fig. 5, and again well beyond the range of parameter space under discussion. The local analysis of Section 3 indicated that there is a region close to M22 where there may be no stable periodic solutions. These path-following curves show that the region where complicated behaviour is likely to occur is bounded by the lines So and so: outside this region local to M22 a stable periodic solution exists. On the boundary of this region occurs not only the mode interaction point M22, but also a second multiple bifurcation point, DS1. It is clear that if the dynamics of this region are to be understood, this second multiple bifurcation point must also be analysed. This cannot be done from the averaged equations as it does not occur local to the trivial solution, so we turn again to numerical methods.

5. Centre manifold maps at multiple bifurcation points The double pendulum equations (1) are two coupled second-order, non-autonomous equations. These can be re-written as four first order non-autonomous equations. In turn, by considering these equations at time T, where T is the period of the excitation, they may be considered as a four dimensional map. Centre manifold reduction may be performed on this map at bifurcation points. A method for performing this reduction numerically was devised by Skeldon [20] and Cliffe [24]. This method is not restricted to the study of bifurcation points which occur at small amplitudes or require an explicit expression for the solution at the multiple bifurcation point and therefore overcomes the restrictions imposed in Section 3 by the use of the method of averaging. We present the results obtained by applying this method to the two multiple bifurcation points M22 and DS1 discussed above. 5.1. Centre manifold map for the mode interaction point

At the mode interaction point there is an interaction between two period-doubling bifurcations. This means that the map has two Floquet multipliers of - 1 . Centre manifold reduction at this point therefore results in a two-dimensional map on the centre manifold. With the choice of M = 7, L = 1.1 and Jl = J2 = 0.24 used above, this map is

v where

( - 1 + Oz2,~q- fl2Q)V "1- CR2V "1- DV 3 + Hu4v + Iu2v 3 + J v 5 + h.o.t.

'

A.C. Skeldon / Physica D 75 (1994) 541-558

al = - 3 2 . 7 0 4 , fll = 19.662,

551

c~2 = - 3 . 3 0 2 0 , 1/2 = - 1 . 4 7 4 3 ,

A = 3.0998,

B = - 4 . 0 0 2 3 , C = 2.1711,

E = 554.01,

F = - 9 4 7 . 7 8 , G = 312.11,

H = 262.10,

I = -326.99,

D = -0.16985,

J = 0.22529,

.~ = A - Am and ~ = g2 - 12m. u is related to the motion of the 0 mode and v is related to the

motion of the ~b mode. Once again, the form of this map is known a priori from the symmetry of the equations and the nature of the two interacting bifurcations: the numerical method allows the computation of the coefficients. A similar bifurcation problem is discussed in Crawford, Knobloch and Riecke [ 10]. The period-doubling nature of the primary bifurcations means that local to M22 there are fixed points of the second iterate of this map, that is, and when Un+2 = Rn and Vn+2 = Vn. This gives the 772 G Z2 bifurcation problem, (a~ad- //~)u

+ A'U 3 + B'UV 2 + E ' u 5 + F'U3V 2 -q- G'uv 4 + h.o.t. = 0,

(oZ2a.q- / / ~ ) v

q- Ctu2v -~- D'V 3 + n ' u 4 v + I'u2v 3 + Jtv5 + h.o.t. = 0 ,

(7)

where the primed coefficients are the coefficients of the second iterate of the map. This bifurcation problem is equivalent to Eq. (5) discussed in Section 3 and results in qualitatively the same bifurcation set and diagrams presented in Figs. 3 and 4. Quantitatively, the numerical bifurcation results also agree well with the analytic results. In Section 3 we noted that a H o p f bifurcation necessarily occurs along the mixed-mode branch, but that neither its direction of bifurcation nor its position could be determined without higher-order terms. In this case the fifth-order terms have been computed and so this restriction no longer applies. From the results of Iooss [25,26] and of Zoladek [27] it is known firstly that a torus bifurcates from the map and secondly that it is unique. The locus of H o p f bifurcations, hi can be calculated from the map and is found to follow closely the line of secondary bifurcation points, s¢. Thus a stable pure 0-mode solution exists to the left of the line So and between So and hi a stable mixed mode exists (see Fig. 10). In order to find the stability of the solution produced by the H o p f bifurcation the expression given by Guckenheimer and Holmes [19] (p. 163) for the stability for a H o p f bifurcation of a map was used. Sufficiently close to the mode interaction point, the H o p f bifurcation is found to be supercritical and therefore produces a stable toms. In order to see how the torus develops as a parameter is changed, the map was iterated local to the mode interaction point. The results are illustrated in Figs. 6a,b and c. Note that the loops are not continuous but are a sequence of points. Only every second iterate is shown: the period-doubling nature of the primary bifurcation means that local to the origin, points in the first quadrant are mapped to the third quadrant and vice versa. Hence the attractor consists of two loops, one in the first quadrant and one in the third. In each case shown, a second toroidal attractor exists in the second and fourth quadrants related to the first by the 7/2 ~) 7/2 symmetry of the problem: which solution is observed depends on the initial conditions (we have not shown the transient part of the solutions). Each map represents a section through a torus of the flow. Fig. 6a shows the torus soon after it has appeared. In Fig. 6b the excitation amplitude has been decreased and the torus has grown. The torus continues to grow until it becomes homoclinic to the origin and is destroyed. Fig. 6c shows the torus just before it loses stability. Thereafter there are no

A.C. Skeldon / Physica D 75 (1994) 541-558

552

(a)

(b)

(c)

Fig. 6. Iteratesof the mode interactionmap for (a) ~ = 2.0x 10-4A" (c) ~ = 2.0 × 1 0 - 4 A" = -0.578 × 10 - 4 .

= - 0 . 5 5 × 10-4;

(b) ~ = 2.0× 1 0 - 4 . 4

=

-0.575× 10-4.

stable solutions to the map. This shows that not only is there a line of Hopf bifurcations emerging from M22, but also there is a line of homoclinic bifurcation points. 5.2. Centre manifold map for the periodic double singular point

Just as centre manifold reduction can be applied to the mode interaction point numerically, it can also be applied to the double singular point. The periodic double singular point occurs when the bifurcation from the pure ¢ mode to the mixed mode passes round the line of saddle-node bifurcations on the pure C-mode branch. The map of the full equations therefore has two unit eigenvalues and again the centre manifold theorem can be used to reduce the four-dimensional map to a two-dimensional map on the centre manifold. This point has only one Z2 symmetry: the Z2 symmetry associated with the ¢ mode has already been broken. The resulting map is

(:)

V + y A + t ~ + Du 2 + E v 2 + Fu2v + Gv 3 -F h.o.t

(8)

where A =3.20, E = -0.27, Oz = 71.45,

B =-6.66, F = 115.57, fl = - 3 3 . 2 4 ,

C=-94.53, G = -0.70, 7 = 4.92,

D=-2.17, ~ = 0.94.

Again .4 and ~ are measures of the distance from the multiple bifurcation point, u is related to the amplitude of the 0 mode and v is related to the amplitude of the ¢ mode. Steady state solutions of Eqs. (8) occur when u,+l = Un and Vn+l = Vn, that is u + a A u + f l ~ u + A u v + B u 3 + C u r 2 + h.o.t. = 0, V + ~.4+ ~ 2 + DU 2 + E v 2 + Fu2v + Gv 3 + h.o.t. = 0.

(9) (10)

This same steady-state problem arises in the study of a saddle-node/Hopf interaction, and in this context has been studied by a number of authors [19,28-30]. For the steady-state problem there are clearly two possible types of solution: u-mode solutions when u = 0, v # 0, and mixed-mode solutions when neither u or v are zero. In turn, these correspond to a periodic pure C-mode solution and a periodic mixed-mode solution, respectively. The corresponding bifurcation set and bifurcation diagrams are presented in Figs. 7 and 8a, b and c. In the discussion of these diagrams we shall again refer to the physical variables 0 and ¢ rather than u and v. In Fig. 8a it is seen that the saddle-node bifurcation point,/¢, gives rise to two pure

553

A.C. Skeldon / Physica D 75 (1994) 541-558 ,14335.14334.14333-

.14332-

\x xx

~ .14331 -

\

Nx \\\

....... 2:--.. "#g....'--.. ,14328-

•14327--t .14326.14325

5~}x

IIIIIIIIl~llllllllllllllllllLltbllllllllllli[I 1.8717

1,8718

1.8719

Exci at

on

~'1 1.872

1,8721

fr~quench,

Fig. 7. Bifurcation set local to the double singular point DS] constructed from the centre manifold map. .............. s~, secondary bifurcation curve from the 0-mode oscillation to the mixed mode. ............. l~, limit points on the t~-mode branch. ..... h2, line of Hopf bifurcations. . . . . . . . bet, line of heteroclinic bifurcations.

J-mode solutions, both unstable. At the point s~ the lower solution branch bifurcates to produce an unstable mixed mode. As Q is increased, s~ moves towards the saddle-node bifurcation point until when Q = Qd the two coincide at DSI, as shown in Fig. 8b. If ~ is increased further, s6 moves on to the upper pure J-mode branch, stabilising the ~ mode. In doing this the stability of the mixed-mode branch also changes. Again, for this to be consistent, a Hopf bifurcation must occur along the mixed-mode branch and a line of such bifurcations must emerge from the periodic double singular point. (This situation is analogous to that given in Golubitsky, Stewart and Schaeffer [29] (p. 428).) This toms is found to bifurcate supercritically, producing a stable toms. In order to discover how this toms changes with parameter, the map (8) is iterated local to the multiple bifurcation point. Fig. 9 shows a sequence of results for increasing .4. Fig. 9a, shows one toroidal attractor; there is a second attractor related to this by the Z2 symmetry. As the parameter is changed, the torus grows, as shown in 9b, until it touches the two fixed points, S + and S - , to form a heteroclinic cycle. This is shown in Fig. 9c. There is a line of such heteroclinic cycles het which emerges from the periodic double singular point.

6. Discussion of results From the numerical path-following of Section 4 it was found that the 2:2 mode interaction region is bounded by two lines of symmetry-breaking curves, so and s~ in Fig. 5. There are two multiple bifurcation points on the boundaries of this region; M22, the mode interaction point and DS1, the periodic double singular point. Both act as organising centres for the dynamics of the region. The previous section has shown that from M22 and D S I emerge lines of Hopf bifurcations and global bifurcations. The position of these lines is only known local to the multiple bifurcation

554

A.C. Skeldon / Physica D 75 (1994) 541-558

.....

- 7 7

..... /

1

/ /

x x

.

.

.

.

.

_

.

.

.

.

.

.

.

.

.

.

.

.

A

A

.

.

(b)

(a)

/ / /

l,

I k \

(c)

A

Fig. 8. Bifurcation diagram for (a) /2 < t2d, (b) 12 = t2d, (c) ~ > ~2d. l~ marks the limit point on the pure ~b mode; s¢ marks the symmetry breaking-bifurcation to the mixed mode; DSI marks the double singular point. S +

¢r

(a)

(b)

(c)

Fig. 9. Iterates of the double singular point map for (a) ~ = 2.0 × 10-4A "~ = = 2.0 x 10-4A"= - 0 . 3 0 0 × 10-4; (c) Q = 2.0 x 10-4A"--- -0.233 x 10 -4.

-0.317 × 10-4;

(b)

points, but they cannot stop in parameter space without reason: a possible picture of their extension is presented in Fig. 10. Between the lines ha and horn there should be a stable torus. The effect of the homoclinic bifurcation is to destroy this torus, thus between the lines horn and pd~ of Fig. 10 the previous sections leave us unable to determine the form of any stable solutions, even at small amplitude. Consideration of DSI showed that in this case, as well as the line of saddle-node bifurcations, l~ and the locus of symmetry-breaking bifurcations s~2, a line of supercritical Hopf bifurcations and a line of heteroclinic connections emerge from this point. The effect of this heteroclinic loop is also to terminate the torus. Hence to the left of the line het and the right of line s~ of Fig. 10 no stable solutions local to DSI a r e found. From the path-following we know that for this region the trivial solution, which is not local to DSI, is stable, and hence this is a candidate for an observable

A.C.

Skeldon

t~1 ?~o,,

o~gq

/,,,, / j : /

oi : !>~,.~"

018

i

/ ./

/.. ,,&'~

017

555

/ P h y s i c a D 75 ( 1 9 9 4 ) 5 4 1 - 5 5 8

,~ U,, ," •

/

016

/ ./.'// ./ /" ..// ./ .'// / ." / . /

/' / // "11 l~

•

/"

o~5

." i I

:. y

:: 1.8

1.825

1.85

1.875

1.9

1925

1.95

Excitation frequency. Q

Fig. 10. Composite bifurcation set. p d e , primary bifurcation curve to the O-mode oscillation. p d ~ , primary bifurcation curve to the ~b-mode oscillation. s e , secondary bifurcation curve from the 8-mode oscillation to the mixed mode. .............. s~, secondary bifurcation curve from the ~b-mode oscillation to the mixed mode. ............... 1~, limit points on the ~b-mode branch. ........... ht, line of Hopf bifurcations. ..... h:, line of Hopf bifurcations. _. _ h e t , line of heteroclinic bifurcations. h o r n , line of homoclinic bifurcations.

- _ _ .....

solution in the full equations. 6.1. N u m e r i c a l Integrations

An advantage of studying ordinary differential equations is that the theoretical results may readily be compared with numerical integrations. These were carried out using the Gear algorithm [ 31 ]. The resulting phase portrait is five-dimensional, including time, and is therefore difficult to visualise. In order to reduce the dimension we consider Poincar6 sections through the plane t = 2T, where T is the period of the excitation. The resulting four dimensional object is viewed on the 0 = ~ = 0 plane. The coordinates 0 and 6 are intimately related, as are ~b and ~ and so this gives a good qualitative representation of the motion on the attractors. There are, of course, still many spurious intersections and care must be taken in interpreting the results. Poincar6 sections were computed for a series of slices in parameter space, each slice defined by a constant excitation amplitude. For each slice, a stable periodic solution was first found outside the 2:2 mode interaction region. A small step in the excitation frequency was then taken and the equations integrated until the transients had decayed away. Typically the solutions converge in less than 100 cycles, but solutions were always integrated for 400 cycles, and 40,000 cycles if it was unclear whether the motion was still transient. Each converged solution was used as a starting solution for the next parameter step. The process was repeated for both increasing and decreasing excitation frequency to pick up any hysteresis. The resulting picture is extremely complicated: we have not sought to classify every piece of the mode interaction region as this in itself is a formidable task but we hope to convey some flavour of the dynamical phenomena that occur.

556

(a)

A.C. Skeldon / Physica D 75 (1994) 541-558

(b)

(c)

(d)

(e)

~~ ""~

(f)

--

(h) .

.

.

.

.

.

.

.

.

.

.

.

.

.

....

Fig. 11. Poincar6 sections produced by numerical integration of the full equations. (a) torus, A (b) symmetric glued torus, A = 0.15,12 = 1.8895; (b) asymmetric glued torus, A = 0.15,12 = glued toms, A = 0.15,I2 = 1.8830; (d) period-doubled asymmetric glued torus, A = 0.15,12 A = 0.15,t2 = 1.8840; (f) near heteroclinicity, A = 0.15,I2 = 1.8780; (g) near homoclinicity, A =

= 0.15,12 = 1.8910; 1.8860; (c) secondary = 1.8850; (e) chaos, 0.18,12 = 1.8720;

We consider a slice at constant excitation amplitude A = 0.15. At an excitation frequency of 12 = 1.91 the Poincar6 section is a single point on the 0 = 0 axis representing a pure ~b-mode solution. As ~2 is decreased to 1.895 this point moves off the 0 = 0 axis as the pure-mode solution bifurcates to a mixed-mode periodic solution. At 12 = 1.891 a torus is found, as predicted. Two different tori exist for the same parameter value and are related to each other by the reflection symmetry. The Poincar6 section for one of these tori is shown in Fig. 1 l a. As the excitation frequency is lowered, these two separate attractors grow until they touch each other through the unstable pure b-mode solution which also exists at these parameter values (Fig. l i b ) . Such a bifurcation is analogous to the gluing bifurcation discussed by Glendinning [32] and Gambaudo, Glendinning and Tresser [33], although in their cases they discuss the gluing of periodic orbits rather than tori. On further decreasing the excitation frequency this symmetric glued attractor first loses symmetry (Fig. 1 lc) and then undergoes a secondary gluing (Fig. 1 ld). Intermingled with this sequence are bands of other complicated attractors, for example the period-doubled asymmetric glued torus shown in Fig. 1 le and bands of chaos (Fig. 1 If). Decreasing the excitation frequency further produces ever more complicated attractors. Finally, all the structure collapses to give the simple loop shown in Fig. 1 If. This orbit is close to heteroclinicity: there is an unstable fixed point close to both the top and the bottom of the loop. This attractor is surprisingly similar to the heteroclinic connection found in the centre manifold map for the double singular point (Fig. 9c). For this amplitude of excitation no hysteresis was observed, however, as the excitation amplitude is raised, hysteresis does begin to appear. Many other types of stable solution are found for example doubling of the heteroclinic loop; chaos restricted to one quadrant of the projected phase space; chaos which visits all four of the quadrants of the projected phase space. Close to the symmetry-breaking line So, a torus is found to exist in a very small region of parameter space. As the excitation frequency is increased, this torus grows rapidly until it becomes homoclinic as predicted. An example of the motion close to homoclinicity is shown in Fig. 1 lh. If the excitation frequency is increased slightly, this homoclinic motion is replaced by larger scale heteroclinic-type orbits.

A.C. Skeldon / Physica D 75 (1994) 541-558

557

7. Conclusion A system of differential equations describing a parametrically excited double pendulum has been considered, with particular attention to the 2:2 mode interaction region. Three principal methods of investigating this problem have been used, each chosen to uncover different facets. The first of these was the method of averaging used in Section 3 to produce an autonomous set of equations to cubic order. These gave us expressions for the primary bifurcation curves and for the 2:2 mode interaction point M22. Centre manifold reduction was applied to these equations at M22 in order to construct the bifurcation set in the neighbourhood of M22. This approach enabled us to find explicit expressions for the coefficients which determine the particular unfolding of this point, and showed that both secondary and Hopf bifurcation lines emerge from M22. It showed that a cubic-order analysis was insufficient to determine even the local behaviour for this region. The method of averaging only allowed a restricted area of parameter space to be studied because of the assumption that the amplitudes of the excitation and of the response were small and that the an excitation frequency was close to 2. Hence, in order to extend the study to a wider region we turned to our second method of approach, numerical path-following. Good qualitative agreement was found for the bifurcation structure for the region around M22. We were able to show that the chaos which is observed in the parametrically excited simple pendulum occurs outside the region of the control parameter plane studied in this paper. Furthermore, a second multiple bifurcation point, a periodic double singular point, DS1, was found in the mode interaction region. Numerical centre manifold reduction was performed at both the multiple bifurcation points M22 and DSI. In both cases, this allowed the position of the lines of Hopf bifurcations emerging from the multiple bifurcation point to be found and these Hopf bifurcations were shown to produce stable tori. These tori were removed from the solution set by global bifurcations. One of the advantages of studying this supposedly simple problem is that the results may be compared with numerical integrations. In this way, evidence for both the homoclinic loop suggested by the analysis of M22 and the heteroclinic loop suggested by a similar study of DSI was found. Many other complicated attractors were found. These included high-order gluings of tori: phenomena which cannot occur in a two-dimensional map, and thus an indication that the local two-dimensional centre manifold is no longer a good approximation. We have adopted the approach of studying the multiple bifurcation points of the system, following the rationale that these points often act as organising centres. It is clear that, whilst the local methods used are valuable and have allowed us to show the existence of homoclinic and heteroclinic loops, in order to fully understand the behaviour away from the multiple bifurcation points, the organising nature of these global orbits must also be considered. Such behaviour is familiar for an orbit which is homoclinic to a saddle-focus steady-state solution [34] where the homoclinic orbit organises many period-doubling sequences. However, there has been little work on the understanding of the sequences of global bifurcations of periodic orbits.

Acknowledgements The author is grateful to Dr T. Mullin for originally suggesting the problem and for many useful discussions throughout the development of this project. She would also like to thank the SERC for

558

A.C. Skeldon / Physica D 75 (1994) 541-558

a Quota Award and the Leathersellers for a scholarship.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [I0] [ 11] [12] [ 13] [14] [15] [16] [17] [ 18 ] [19 ] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31 ] [32] [33] [34] [35]

A.C. Skeldon and T. Mullin, Phys. Lett. A 166 (1992) 224. M. Faraday, Philos. Trans. Roy. Soc. London 121 (1831) 319. S. Ciliberto and J.P. GoUub, Phys. Rev. Lett. 52 (1984) 922. S. Ciliberto and J.P. GoUub, J. Fluid. Mech. 158 (1985) 381. T.B. Benjamin and F. Ursell, Proc. Roy. Soc. London Ser. A 225 (1954) 505. E. Meron and I. Procaccia, Phys. Rev. A 34 (1986) 3221. Z.C. Feng and P.R. Sethna, J. Fluid Mech. 199 (1989) 495. M. Umeki and T. Kambe, J. Phys. Soc. Jpn. 58 (1989) 140. J.D. Crawford, E. Knobloch and H. Riecke, Phys. Lett. A 135 (1989) 20. J.D. Crawford, E. Knobloch and H. Riecke, Physica D 44 (1990) 340. N. Kryloff and N. Bogoliuboff, Introduction to non-linear mechanics (Princeton Univ. Press, Princeton, 1947). J. Carr, Appl. Math. Sci. 35 (1981). D.W. Jordan and P. Smith, Oxford Applied Maths and Computing Science Series (Oxford Univ. Press, Oxford, 1987). R.W. Leven and B.P. Koch, Phys. Lett. A 86 (1981) 71. B.P. Koch and R.W. Leven, Physica D 16 (1985) 1. R.W. Leven, B. Pompe, C. Wilke and B.P. Koch, Physica D 16 (1985) 371. P.J. Bryant and J.W. Miles, J. Austral. Math. Soc. Ser. B 30 (1990) 42. O. Chubykalo, private communication. J. Guckenheimer and P. Holmes, Appl. Math. Sci. 42 (1986). A.C. Skeldon, D.Phil. thesis, Oxford (1990). M. Golubitsky and D.G. Schaeffer, Appl. Math. Sci. 51 (1985 ). A.D. Jepson and A. Spence, SIAM J. Numer. Anal. 22 (1985) 347. A.D. Jepson and A. Spence, SIAM J. Numer. Anal. 22 (1985) 736. K.A. Cliffe, private communication. G. Iooss, C.R. Acad. Sci. Paris 296 (1983) 27. G. Iooss, C.R. Acad. Sci. Paris 296 (1983) 113. K. Zoladek, J. Diff. Eq. 67 (1987) 1. W.F. Langford and G. Iooss, in: Workshop on Bifurcation Problems and their Numerical Solution, H. Mittelman and H. Weber, eds. (Springer, Berlin, 1990). M. Golubitsky, I. Stewart and D.G. Schaeffer, Appl. Math. Sci. 69 (1990). V. Kirk, Phys. Lett. A 154 (1991) 243. C.W. Gear, Numerical initial value problems in ordinary differential equations (Prentice-Hall, Englewood Cliffs, NJ, 1971). P. Glendinning, Phys. Lett. A 103 (1984) 163. J.M. Gambaudo, P. Glendinning and C. Tresser, J. Phys. Lett. 46 (1985) 653. P. Glendinning and C. Sparrow, J. Slat. Phys. 35 (1984) 645. L.D. Landau and E.M. Lifschitz, Mechanics. Course of theoretical physics, Vol. 1 (Pergamon Press, Oxford, 1960).