Computing global orbits of the forced spherical pendulum

Physica D 165 (2002) 1–11

Computing global orbits of the forced spherical pendulum Thomas J. Bridges∗ , Kyriakos V. Georgiou Department of Mathematics and Statistics, University of Surrey, Guildford, Surrey GU2 7XH, UK Received 21 November 2000; received in revised form 3 December 2001; accepted 1 February 2002 Communicated by E. Ott

Abstract Orbits of the spherical pendulum with time-periodic forcing are considered. A numerical framework is developed, which allows orbits to explore the entire “globe”: the spherical pendulum is considered as an invariant manifold in an ambient six-dimensional Euclidean space. The numerical integrator is the second-order Störmer–Verlet method coupled with the Shake–Rattle algorithm. The algorithm preserves numerically the phase space of the sphere, which is a manifold, to machine accuracy. Poincaré sections, restricted to the configuration space, are used to illustrate the transition from oscillatory behavior to chaotic solutions, as the amplitude of the pivot motion is changed. The qualitative change in the Poincaré sections from regular to chaotic behavior is in excellent qualitative agreement with corresponding computations of the Lyapunov exponents (LEs). The LEs are also computed using a novel variant of the Shake–Rattle algorithm. The results show that irregular behavior can explore the entire sphere—even at low forcing amplitudes—and therefore local methods which parameterize only part of the sphere are inadequate in general, and may lead to spurious dynamics. The numerical framework provides a tool for detailed investigation of the symmetric chaos of the forced spherical pendulum. © 2002 Published by Elsevier Science B.V. Keywords: Euclidean space; Shake–Rattle algorithm; Chaotic; Forced spherical pendulum

1. Introduction The spherical pendulum is a model for a range of physical phenomena, and is a fundamental model for mechanical systems with a circular symmetry. The spherical pendulum is integrable, and a comprehensive treatment can be found in Cushman and Bates [4]. The forced spherical pendulum is not in general integrable, and it has been known for some time that the motion can be quite complicated, and its dynamics has been characterized as chaotic in a number of studies. The irregular motion of the spherical pendulum seems to have features which are distinctly more complicated than the motion of a forced planar pendulum. Because of the presence of a continuous symmetry, the forced spherical pendulum is also of interest as a model for the study of symmetric chaos. In the first systematic study of the dynamics of the forced spherical pendulum, Miles [9] showed that the weakly nonlinear motion could be highly irregular. The spherical pendulum in this case was a model for the sloshing of fluid in a partially filled cylindrical container, motivated by aerospace applications (cf. [11]). Revisiting this problem with the tools of modern dynamical systems, Miles [10] studied the weakly nonlinear normalized equations and showed ∗

Corresponding author. Tel.: +44-1483-259642; fax: +44-1483-876071. E-mail address: [email protected] (T.J. Bridges). 0167-2789/02/$ – see front matter © 2002 Published by Elsevier Science B.V. PII: S 0 1 6 7 - 2 7 8 9 ( 0 2 ) 0 0 3 8 1 - 0

2

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

that the irregular behavior was associated with a positive Lyapunov exponent (LE), a widely agreed signature of chaos. Experimental evidence [14] and further numerical studies of Miles’ model equations [15] support the evidence for complex chaotic behavior. However, all the above theoretical studies were based on approximate model equations which are valid only in a small neighborhood of the stable equilibrium position. Bryant [3] extended this work to finite-amplitude dynamics of the forced spherical pendulum, and identified parameter regions where the motion is chaotic. However, the sphere was approximated by the lower hemisphere, and therefore all orbits are restricted to lie below the equator. Aston [2] also studies finite-amplitude dynamics, using a combination of numerical and bifurcation techniques. The concept of blow-out bifurcation is used to give a more precise characterization of the chaos by finding a bifurcation of chaotic planar states to chaotic non-planar states. However, all the analyses were carried out on a single coordinate chart for the sphere, based on spherical coordinates, which is not regular at all points of the sphere. In this paper, we formulate and analyze the spherical pendulum as a constrained dynamical system. The ambient phase space is Euclidean six-dimensional space, with the spherical pendulum dynamics associated with a four-dimensional sub-manifold: the tangent bundle of the sphere. In this setting there is no restriction on orbits: they can freely explore any point on the “globe”. Let x0 (t) ∈ R3 be the specified position of the origin, and let q(t) = (q1 (t), q2 (t), q3 (t)) be coordinates for the position of the pendulum mass, relative to the position of the origin. As a constrained dynamical system, the governing equations for the spherical pendulum take the simple form q˙ =

1 p, m

p˙ = −mx¨ 0 (t) − mge2 − qλ(t),

(q, p) ∈ R3 × R3 ,

(1.1)

where m is the pendulum mass, g the gravitational constant and λ(t) a Lagrange multiplier associated with the constraint q · q = 2 . The difficulty with solving this equation is the required coupling with the sphere constraint, and the hidden constraint, p · q = 0. Therefore, the complete system is (1.1) with the tangent bundle constraints: q · q = 2 and p · q = 0. In this paper, we will concentrate on developing a numerical framework for solving these equations. The central difficulty is the preservation of the constraints numerically. Here we will use the Shake–Rattle algorithm which was originally proposed for molecular dynamics simulations [1], and then later found to have important geometric properties [7,8]. In particular, the method is symplectic, and can preserve the constraint manifold to machine accuracy. In Fig. 1, an example of the range of dynamics that can be computed with this numerical framework is shown. The figure shows the configuration space of the Poincaré section for the forced spherical pendulum at a fixed value

Fig. 1. Configuration space of the Poincar´e section for the forced spherical pendulum with m = 1, g = 1, = 1, forcing frequency ω = and forcing amplitude ε = 0.5 (details of forcing function and initial data are in Section 4).

√

2,

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

3

of the parameters. The position of the pendulum mass on the sphere at t = nT where T is the period of the periodic forcing and n a natural number is shown. It is clear that virtually the entire sphere is explored by this single orbit, and so parameterizing only part of the sphere or using a weakly nonlinear approximation would fail to capture the dynamics. Further results, illustrating the transition from regular to irregular behavior, can be found in Sections 4 and 5. The orbit in Fig. 1 is certainly irregular enough to suspect that it might be chaotic, and therefore it is natural to compute LEs. Here, we compute the LEs by combining linearized Shake–Rattle with the symplectic structure of the flow, to construct a relatively efficient algorithm for computing the largest LE. The results we present here are the first numerical calculations of LEs for the forced spherical pendulum, with no restriction on the amplitude or position of the orbit. Our results confirm that the flow in Fig. 1 has a positive LE, and we give other examples of parameter values with positive LEs. An outline of the paper is as follows. In Section 2, the precise governing equations are derived, as a constrained Hamiltonian system in R6 . In Section 3, the Shake–Rattle algorithm is applied to the forced spherical pendulum equations. The two nontrivial basic states of the unforced spherical pendulum are the conical pendulum solutions, and the planar—zero momentum—solutions. Both these classes of periodic solutions can be deduced purely from symmetry properties (cf. [12]). In Section 4, the effect of forcing on the conical pendulum orbits is considered. The main tool is the Poincaré section, which illustrates the qualitative change of the orbits, as the amplitude of the forcing is increased. In Sections 5 and 6, the planar and planar drift orbits are considered with forcing. Here the tools are Poincaré sections and computation of the largest LE. In Section 6, an algorithm for the LEs for a constrained system is developed. The appearance of positive LEs is in excellent qualitative agreement with the appearance of irregular orbits in the Poincaré section.

2. The forced spherical pendulum as a constrained dynamical system A schematic representation of the system is shown in Fig. 2. The components of the vector q(t) = (q1 (t), q2 (t), q3 (t))

Fig. 2. Schematic representation of the forced spherical pendulum, with the motion of the point of support prescribed.

4

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

represent the displacement of the mass m in the x, y and z directions, respectively, relative to the pivot point. The components of q are related by the equation of constraint q12 + q22 + q32 = 2 . The position of the base point s is specified, and can move in any direction. The position of s relative to the fixed origin o is given by some time-dependent vector x0 (t) = (x0 (t), y0 (t), z0 (t)). Typical examples for the motion of s are: (a) Purely horizontal periodic motion of amplitude ε and frequency ω, x0 (t) = ε sin (ωt),

y0 (t) = 0,

z0 (t) = 0.

(b) Purely vertical periodic motion of amplitude ε and frequency ω, x0 (t) = 0,

y0 (t) = ε sin (ωt),

z0 (t) = 0.

(c) Spherical quasiperiodic forcing of amplitude ε, x0 (t) = ε cos (ω1 t) sin (ω2 t),

y0 (t) = ε sin (ω1 t) sin (ω2 t),

z0 (t) = −ε cos (ω2 t),

where ω1 /ω2 is irrational. The absolute coordinates for the position of the mass m at time t are x(t) = (x(t), y(t), z(t)) with x(t) = x0 (t) + q1 (t),

y(t) = y0 (t) + q2 (t),

z(t) = z0 (t) + q3 (t).

˙ The velocity is u(t) = (u(t), v(t), w(t)) with u(t) = x(t). The Lagrangian for the system is then L(q, q˙ ) = 21 mu · u − V (q) − λ(t) 21 (q · q − 2 ), where λ is a Lagrange multiplier representative of the force required to keep the pendulum mass on the sphere. The function V (q) represents the potential energy (PE) V (q) = mg(y0 (t) + q2 ) = mgy0 (t) + mge2 · q, where g is the positive gravitational constant and e2 = (0, 1, 0). Substitution of the above quantities into the Lagrangian results in L(q, q˙ ) = 21 m(x˙ 0 (t) + q˙ ) · (x˙ 0 (t) + q˙ ) − mge2 · x0 (t) − mge2 · q − 21 λ(t)(q · q − 2 ). The first variation of L, with fixed endpoints leads to the governing equations m(x¨ 0 + q¨ + ge2 ) + λq = 0.

(2.1)

In terms of the position and linear momentum vector, p = (p1 , p2 , p3 ) = m˙q, the governing equations are q˙ =

1 p, m

p˙ = F (t) − qλ(t)

(2.2)

with F (t) = −mx¨ 0 (t) − mge2 . In the position–momentum formulation, the initial data is constrained to lie on the tangent bundle of the sphere: q · q = 2 and p · q = 0. The flow of Eq. (2.2) is symplectic in the sense that d dp ∧ dq = 0, dt where dq and dp satisfy the variational equations associated with (2.2).

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

5

3. Symplectic constraint-preserving numerics: Shake–Rattle The numerical integration of the constrained system is developed in this section. In addition to the constraint surface, it is desirable to preserve the symplectic structure of the system. Symplectic integrators are now well developed (cf. [16]), and have the important advantage of preserving energy— when the system is conservative—over very long time intervals, up to small fluctuations (cf. [13]), and generally have excellent stability properties. The Shake–Rattle algorithm which was developed for natural mechanical systems in molecular dynamics is designed precisely to preserve the constraints and the symplectic structure [1,7,8]. This algorithm is based on an explicit Störmer–Verlet time integrator for the unconstrained system, which is second-order accurate, coupled with Newton iteration to determine the Lagrange multipliers at each step. This algorithm leads to the following discrete approximation of the governing equation (2.2) and constraints. Shake–Rattle for the nonlinear system: pn+1/2 = pn + 21 h(F n − qn λn ), qn+1 = qn +

h p , m n+1/2

(3.1) (3.2)

pn+1 = pn+1/2 + 21 h(F n+1 − qn+1 µn ),

(3.3)

2 = qn+1 · qn+1 ,

(3.4)

0 = qn+1 · pn+1 .

(3.5)

The Lagrange multipliers λn and µn are determined by imposing the constraints (3.4) and (3.5), respectively. To determine λn , substitute (3.2) into (3.4), using (3.1). The result is the following quadratic equation for λn :
2
2 h(p1,n + 21 h(F1,n − q1,n λn )) h(p2,n + 21 h(F2,n − q2,n λn )) f (λn ) = q1,n + + q2,n + m m
2 h(p3,n + 21 h(F3,n − q3,n λn )) + q3,n + − 2 . (3.6) m This equation could be solved for λn using the quadratic formula, but it is simpler to use iteration based on Newton’s (0) method. Given an initial estimate for λn , denoted by λn —e.g. from the previous timestep—the Newton algorithm, within the nth timestep is (j +1) λn

=

(j ) λn

(j )

−

f (λn ) (j )

f (λn )

,

j = 0, 1, . . . .

(3.7)

The Lagrange multiplier µn is determined by imposing the constraint (3.5). Substitution of (3.3) into (3.5) results in a linear equation for µn , which is easily solved without any iteration, µn =

2(q1,n+1 p1,n+1/2 + q2,n+1 p2,n+1/2 + q3,n+1 p3,n+1/2 ) h(q1,n+1 q1,n+1 + q2,n+1 q2,n+1 + q3,n+1 q3,n+1 ) q1,n+1 F1,n+1 + q2,n+1 F2,n+1 + q3,n+1 F3,n+1 + . q1,n+1 q1,n+1 + q2,n+1 q2,n+1 + q3,n+1 q3,n+1

(3.8)

In summary, Newton iteration is used first to determine λn . With λn , pn+1 and qn+1 are determined using (3.1) and (3.2), respectively. Eq. (3.8) is then solved explicitly for µn . With µn substituted into (3.3), pn+1 can be determined.

6

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

If the Newton iteration is carried to convergence, the tangent bundle of the sphere is preserved at each step to machine accuracy.

4. The effect of forcing on conical pendulum orbits The conical pendulum orbits are one of the two basic nontrivial solutions of the unforced spherical pendulum, and are relative equilibria associated with the circular symmetry [12]. In the lower hemisphere, they can be characterized as the intersection of the plane q2 = a with − < a ≤ 0 and the sphere. Explicitly q1 (t) = r cos θ (t),

q2 (t) = a,

q3 (t) = r sin θ(t)

(4.1)

with r=

 2 − a 2 ,

a=−

mg , λ

θ (t) =

√ λt + θ0 .

If we assume that the momentum is specified (e.g. by specifying the initial data), then the value of the momentum determines the value of λ. The momentum is I = q˙3 q1 − q˙1 q3 =

1 p 3 q1 − p 1 q3 , m

(4.2)

√ and so, when evaluated on a conical pendulum solution, I = ( 2 − a 2 ) λ. 4.1. Poincaré section in the configuration space When the spherical pendulum is forced harmonically, the natural Poincaré section is the tangent bundle of the sphere at t = nT, where T is the period of the harmonic forcing. In this section, we illustrate the configuration space of this Poincaré section: for t = nT, we plot the location on the sphere of the position of the pendulum mass. The forcing function is the position of the pivot point of the pendulum and is specified to be x0 (t) = (x0 (t), y0 (t), z0 (t)) with x0 (t) = 0, y0 (t) = ε sin ωt, z0 (t) = 0 √ with ω = 2. The forcing frequency is equal to the frequency of the conical pendulum orbit, and the initial data is chosen to be on a conical pendulum state. Fig. 3 shows the resulting orbits for a range of values of ε, the forcing amplitude. In Fig. 3(a), the amplitude of forcing is zero and the resulting orbit is a conical pendulum solution. Since the forcing frequency is equal to the frequency of the conical pendulum orbit, the orbit in the Poincaré section is just a point (at the right edge of the sphere in the picture). However, when the forcing is turned on the orbit explores more of the globe with increasing complexity as ε increases. In each of the figures, the initial data is the same. When ε ≥ 0.5, the orbit becomes irregular and begins to explore the entire sphere. Note that an orbit which explores the entire sphere does not necessarily explore the entire phase space. Therefore, it is highly unlikely that such orbits are ergodic, or even closely approximated by an ergodic trajectory. The main point here is that these results show that a numerical algorithm needs the capability to explore the entire sphere—rather than approximate it by the lower hemisphere or other subset of the sphere.

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

7

Fig. 3. Configuration space of the Poincar´e sections for √ forced conical pendulum solutions, sampled once per√period of the forcing. Gravity points in the negative q2 direction, and m = 1, ω = 2, g = 1, = 1. The initial conditions are: q1 (t) = 21 3, q2 (t) = − 21 , q3 (t) = 0.0, √ p1 (t) = 0.0, p2 (t) = 0.0, p3 (t) = 3/2. The timestep is h = 0.001410439028, and the number of timesteps is 107 .

5. Effect of forcing on drift orbits: Poincaré sections The second basic non-trivial state of the unforced spherical pendulum is the family of planar states, which are obtained by setting the angular momentum to zero. In Fig. 4(a), an example of the unforced planar pendulum is shown. Because the forcing period differs from the natural period, the Poincaré section shows a “line” associated with this orbit. In the rest of the pictures in Fig. 4, non-zero momentum is added to the initial data, resulting in a quasiperiodic drift orbit when ε = 0, and increasingly complex dynamics when the forcing is turned on. The form of the forcing function is the same as in Section 4. The largest LE associated with these orbits will be computed in Section 6.

8

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

Fig. 4. Configuration space of the Poincar´e sections, sampled once per period of the forcing, for the forced drift orbits. Gravity point in the negative q2 direction. For (b–f), initial conditions are: q1 = 0.29, q2 = −0.9, q3 = 0.13, p1 = 0.28, p2 = 0.08, p3 = −1.0. For (a), the initial conditions are the same as in (b)–(f) except that p1 , p3 are set to zero, leading to zero momentum and a planar state. The parameters are m = 1, ω = 2, g = 1, l = 1, h = 0.0009973310014, and the number of timesteps is 107 .

6. Computing LEs for global orbits LEs are a qualitative tool for studying the complexity of orbits, particularly the exponential divergence of nearby orbits, and are obtained by integrating the linearized system. Since the phase space of the constrained system is the tangent bundle of the sphere, TS2 , we need to integrate the linearized discrete system on the tangent space of TS2 . A straightforward linearization of the discrete system (3.1)–(3.5) leads to the following system. Shake–Rattle for the linear system: dpn+1/2 = dpn − 21 h(dqn λn + qn dλn ), dqn+1 = dqn +

h , dp m n+1/2

(6.1) (6.2)

dpn+1 = dpn+1/2 − 21 h(dqn+1 µn + qn+1 dµn ),

(6.3)

0 = qn+1 · dqn+1 ,

(6.4)

0 = qn+1 · dpn+1 + dqn+1 · pn+1 .

(6.5)

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

9

The restriction from R6 to the tangent space of TS2 is obtained by linearization of the tangent bundle constraints, and the associated variations of the Lagrange multipliers dλn and dµn , which need to be computed at each step. They are computed using a linearized version of Shake–Rattle. Since dλn appears in (6.1), we substitute expression (6.2) into (6.4). Then, using (6.1), we obtain
 h(dp1,n − 21 h(dλn q1,n + λn dq1,n )) f (dλn ) = 2q1,n+1 dq1,n + m
 1 h(dp2,n − 2 h(dλn q2,n + λn dq2,n )) + 2q2,n+1 dq2,n + m
 h(dp3,n − 21 h(dλn q3,n + λn dq3,n )) + 2q3,n+1 dq3,n + . m This expression is linear in dλn , and so it is straightforward in principle to solve explicitly for dλn . However, the calculation is lengthy and therefore, we use one step of the Newton algorithm. To determine dµn we proceed in a similar way. Substitute the expressions for dpn+1 and dqn+1 into (6.5). This equation is linear in dµn and the expression is simple enough to solve explicitly for dµn . The resulting expression is then substituted into (6.3). In summary, given qn , pn , λn and µn , we first compute dλn using one step of the Newton iteration. Substitution of dλn into (6.1) then allows calculation of dpn+1/2 and dqn+1 . An expression for dµn is explicitly calculated as described above, and then dpn+1 is calculated using (6.3). If dλn is calculated so that (6.4) is satisfied to machine accuracy. One n-step of the above algorithm can be described abstractly as Y n+1 = An Y n ,

Y n = (dqn , dpn )

(6.6)

with Y 0 prescribed and of unit length. It follows from results of [7] that the mapping An is symplectic. An explicit expression for An can be given but will not be needed in the sequel. 6.1. An algorithm for computing the largest LE With the system in the abstract form (6.6), we are now in a position to apply standard results on the largest LE of a discrete dynamical system. Using Lemma 2.2 of Dieci and van Vleck [5], the largest LE is given by n−1

1 log|bj |. n→∞ n

Λ∞ = lim

(6.7)

j =0

An expression for bj will be given shortly. There are two ways this expression needs to be modified for the calculations. First, it is not necessary to compute this sum at every timestep. The timestep is typically very small because the Störmer–Verlet algorithm is explicit, but the sum in the definition of the LE may not change much in one step. Therefore, we typically take K steps before evaluating the sum. Secondly, to take into account that the discrete system comes from a flow, we need to add a scale factor (K (t)−1 where (t is the step size (see Eq. (16) of [6]), and so Λ∞ = lim ΛN , N→∞

ΛN =

N−1 1 1  log|bj |, K (t N j =0

where N indexes the steps at which the sum is evaluated.

(6.8)

10

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

The sequence bj , j = 1, 2, . . . , is determined as follows. Let zj = Y j  and express Y j = Qj zj , then Qj has unit length, and zj satisfies zn+1 = bn zn ,

bn = Qn+1 , An Qn ,

(6.9)

where ·, · is a standard Euclidean inner product on R6 . Given Y j , we normalize it to have unit length, denoted by Qj . Then Y j +1 = Aj Qj and   Y j +1 bj = Qj +1 , Aj Qj  = (6.10) , Aj Qj = Aj Qj , Y j +1  and so bj = Y j +1  when Y j is normalized to have unit length. To avoid summing the series in the expression for ΛN at each N -step, a recursion formula can be developed   N N−1   1 1 1 1  ΛN+1 = log bj = log bj + bN  , K (t N + 1 K (t N + 1 j =0

j =0

and so ΛN+1 =

N 1 1 ΛN + log bN . N +1 K (t N + 1

(6.11)

To illustrate the role of LEs, we compute them for the parameter values associated with Fig. 4 with ε, the amplitude of forcing, varying. The results are given in Table 1. In each case, the simulation was run for 107 timesteps. The results show excellent qualitative agreement with the Poincaré sections in Fig. 4. For 0 ≤ ε ≤ 0.125, the flow in Fig. 4 appears to be quasiperiodic—or in general regular—and the associated LEs are almost zero. However, it is apparent in Fig. 4 that a transition occurs between ε = 0.125 and 0.126. When ε = 0.126, the orbits are exploring most of the globe, and the largest LE reflects this complexity by becoming positive between ε = 0.125 and 0.126. Further increase in ε shows complex dynamics in the Poincaré sections, and the largest LE remains positive and increases with increasing ε. The positive LE is a good indication that the orbits in Fig. 4 are chaotic for ε > 0.126. A similar calculation at the parameter values associated with Fig. 1 showed that the largest LE is positive, with a value of Λ∞ ≈ 0.030014. Table 1 Computed LEs associated with the results in Fig. 4 as ε is varieda ε

Largest LE

0.0 0.1 0.12 0.125 0.126 0.13 0.2 0.3 0.4 0.5 0.6

1.4655 × 10−5 1.51907 × 10−5 1.54312 × 10−5 1.53176 × 10−5 0.0904216 0.0937221 0.124655 0.16336 0.184631 0.20834 0.226179

The parameters are ω = 2.0, m = 1.0, g = 1.0 and = 1.0, with initial data for the flow q1 = 0.29, q2 = −0.9, q3 = 0.13, p1 = 0.28, p2 = 0.08, p3 = −1.0. For the linearized system random initial data of unit length is used. a

T.J. Bridges, K.V. Georgiou / Physica D 165 (2002) 1–11

11

7. Concluding remarks One of the two main results of this paper is that global orbits of the—forced or unforced—spherical pendulum can be computed accurately by integrating in the ambient six-dimensional phase space, and constraining orbits to lie on the four-dimensional sub-manifold which defines the spherical pendulum. By using the Shake–Rattle algorithm, the tangent bundle is preserved to machine accuracy, and therefore numerical errors are restricted to location on the manifold. The second main result is that LEs associated with global orbits can be calculated using a variant of the Shake–Rattle algorithm for the linearized constrained system. The natural Poincaré section for orbits, when the forcing is T -periodic, is the tangent bundle of the sphere at each T . We showed that the configuration space of this Poincaré section provides interesting qualitative information about orbits, and a change from regular to chaotic behavior was confirmed by the calculation of the largest LE. Only the undamped system was considered, but it is straightforward to include damping in the numerical framework. Indeed, it is known that symplectic methods are also excellent for systems with uniform contraction, e.g. systems with constant linear damping, which is a standard model for damping in mechanical systems. The orbits of the forced system showed clear evidence of chaotic behavior. Moreover, it is clear from the results presented here that any general study of chaos in the forced spherical pendulum will require an algorithm which allows exploration of the entire sphere, as developed here. There are many open questions about the nature of the chaos in the forced spherical pendulum, since evidence suggests that it is more dramatic than in the forced planar pendulum. The numerical framework presented in this paper provides a robust tool for investigation of the qualitative properties of symmetric chaos in this system.

Acknowledgements Useful discussions with Philip J. Aston, Gianne Derks and Sebastian Reich are gratefully acknowledged. The authors are particularly grateful to Sebastian Reich for suggesting how Shake–Rattle could be modified to compute LEs. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

H.C. Anderson, Rattle: a ‘velocity’ version of the Shake algorithm for molecular dynamics calculations, J. Comp. Phys. 52 (1983) 24–34. P.J. Aston, Bifurcations of the horizontally forced spherical pendulum, Comp. Meth. Appl. Mech. Eng. 170 (1999) 343–353. P.J. Bryant, Breakdown to chaotic motion of a forced, damped, spherical pendulum, Physica D 64 (1993) 324–339. R.H. Cushman, L.M. Bates, Global Aspects of Classical Integrable Systems, Birkhauser, Basel, 1997. L. Dieci, E.S. van Vleck, Computation of a few Lyapunov exponents for continuous and discrete dynamical systems, Appl. Num. Math. 17 (1995) 275–291. K. Geist, U. Parlitz, W. Lauterborn, Comparison of different methods for computing Lyapunov exponents, Prog. Theoret. Phys. 83 (1990) 875–893. B. Leimkuhler, R.D. Skeel, Symplectic numerical integrators in constrained Hamiltonian dynamics, J. Comp. Phys. 112 (1994) 117–125. B. Leimkuhler, S. Reich, Geometric Integrators in Hamiltonian Mechanics, Cambridge University Press, Cambridge, 2002, in press. J.W. Miles, Stability of forced oscillations of a spherical pendulum, Quart. Appl. Math. 20 (1962) 21–32. J.W. Miles, Resonant motion of a spherical pendulum, Physica D 11 (1984) 209–323. J.W. Miles, Resonantly forced surface waves in a circular cylinder, J. Fluid Mech. 149 (1984) 15–31. J.A. Montaldi, R.M. Roberts, I.N. Stewart, Periodic solutions near equilibria of symmetric Hamiltonian systems, Phil. Trans. R. Soc. London A 325 (1988) 237–293. S. Reich, Backward error analysis for numerical integrators, SIAM J. Num. Anal. 36 (1999) 1549–1570. D.J. Tritton, Ordered and chaotic motion of a forced spherical pendulum, Eur. J. Phys. 7 (1986) 162–169. D.J. Tritton, M. Groves, Lyapunov exponents for the Miles’ spherical pendulum equations, Physica D 126 (1999) 83–98. J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994.