# 4.5: Harmonically-driven, linearly-damped, plane pendulum

- Page ID
- 9583

The harmonically-driven, linearly-damped, plane pendulum illustrates many of the phenomena exhibited by non-linear systems as they evolve from ordered to chaotic motion. It illustrates the remarkable fact that determinism does not imply either regular behavior or predictability. The well-known, harmonically-driven linearly-damped pendulum provides an ideal basis for an introduction to non-linear dynamics^{1}.

Consider a harmonically-driven linearly-damped plane pendulum of moment of inertia \(I\) and mass \(m\) in a gravitational field that is driven by a torque due to a force \(F(t)=F_{D}\cos \omega t\) acting at a moment arm \(L\). The damping term is \(b\) and the angular displacement of the pendulum, relative to the vertical, is \(\theta\). The equation of motion of the harmonically-driven linearly-damped simple pendulum can be written as

\[I \ddot{\theta}+b\dot{\theta}+mgL\sin \theta =LF_{D}\cos \omega t \label{4.28}\]

Note that the sinusoidal restoring force for the plane pendulum is non-linear for large angles \(\theta\). The natural period of the free pendulum is

\[\omega _{0}=\sqrt{\frac{mgL}{I}}\]

A dimensionless parameter \(\gamma\), which is called the **drive strength,** is defined by \[\gamma \equiv \frac{F_{D}}{mg}\]

The equation of motion \ref{4.28} can be generalized by introducing dimensionless units for both time \(\tilde{t}\) and relative drive frequency \( \tilde{\omega}\) defined by

\[\tilde{t}\equiv \omega _{0}t\hspace{1in}\tilde{\omega}\equiv \frac{\omega }{ \omega _{0}}\]

In addition, define the inverse damping factor \(Q\) as

\[Q\equiv \frac{\omega _{0}I}{b}\]

These definitions allow Equation \ref{4.28} to be written in the dimensionless form \[\frac{d^{2}\theta }{d\tilde{t}^{2}}+\frac{1}{Q}\frac{d\theta }{d\tilde{t}} +\sin \theta =\gamma \cos \tilde{\omega}\tilde{t} \label{4.33}\]

The behavior of the angle \(\theta\) for the driven damped plane pendulum depends on the drive strength \(\gamma\) and the damping factor \(Q\). Consider the case where Equation \ref{4.33} is evaluated assuming that the damping coefficient \(Q=2\), and that the relative angular frequency \(\tilde{\omega}= \frac{2}{3},\) which is close to resonance where chaotic phenomena are manifest. The Runge-Kutta method is used to solve this non-linear equation of motion.

## Close to Linearity

For drive strength \(\gamma =0.2\) the amplitude is sufficiently small that \( \sin \theta \simeq \theta ,\) superposition applies, and the solution is identical to that for the driven linearly-damped linear oscillator. As shown in Figure \(\PageIndex{1}\), once the transient solution dies away, the steady-state solution asymptotically approaches one attractor that has an amplitude of \( \pm 0.3\) radians and a phase shift \(\delta\) with respect to the driving force. The abscissa is given in units of the dimensionless time \(\tilde{t} =\omega _{0}t\). The transient solution depends on the initial conditions and dies away after about \(5\) periods, whereas the steady-state solution is independent of the initial conditions and has a state-space diagram that has an elliptical shape, characteristic of the harmonic oscillator. For all initial conditions, the time dependence and state space diagram for steady-state motion approaches a unique solution, called an "**attractor**", that is, the pendulum oscillates sinusoidally with a given amplitude at the frequency of the driving force and with a constant phase shift \(\delta\), i.e.

\[\theta (t)=A\cos (\omega t-\delta ).\]

This solution is identical to that for the harmonically-driven, linearly-damped, linear oscillator discussed in chapter \(3.6.\)

## Weak nonlinearity

Figure \(\PageIndex{1}\) shows that for drive strength \(\gamma =0.9\), after the transient solution dies away, the steady-state solution settles down to one attractor that oscillates at the drive frequency with an amplitude of slightly more than \(\frac{\pi }{2}\) radians for which the small angle approximation fails. The distortion due to the non-linearity is exhibited by the non-elliptical shape of the state-space diagram.

The observed behavior can be calculated using the successive approximation method discussed in chapter \(4.2\). That is, close to small angles the sine function can be approximated by replacing

\[\sin \theta \approx \theta -\frac{1}{6}\theta ^{3}\]

in Equation \ref{4.33} to give

\[\ddot{\theta}+\frac{1}{Q}\dot{\theta}+\omega _{0}^{2}\left( \theta -\frac{1}{ 6}\theta ^{3}\right) =\gamma \cos \tilde{\omega}\tilde{t} \label{4.35}\]

As a first approximation assume that

\[\theta (\tilde{t})\approx A\cos (\tilde{\omega}\tilde{t}-\delta )\]

then the small \(\theta^{3}\) term in Equation \ref{4.35} contributes a term proportional to \(\cos ^{3}(\tilde{\omega}\tilde{t}-\delta )\). But

\[\cos ^{3}(\tilde{\omega}\tilde{t}-\delta )=\frac{1}{4}\left( \cos 3(\tilde{ \omega}\tilde{t}-\delta )+3\cos (\tilde{\omega}\tilde{t}-\delta )\right)\]

That is, the nonlinearity introduces a small term proportional to \(\cos 3(\omega t-\delta )\). Since the right-hand side of Equation \ref{4.35} is a function of only \(\cos \omega t,\) then the terms in \(\theta ,\dot{\theta},\) and \(\ddot{\theta}\) on the left hand side must contain the third harmonic \( \cos 3(\omega t-\delta )\) term. Thus a better approximation to the solution is of the form

\[\theta (\tilde{t})=A\left[ \cos (\tilde{\omega}\tilde{t}-\delta )+\varepsilon \cos 3(\tilde{\omega}\tilde{t}-\delta )\right]\]

where the admixture coefficient \(\varepsilon <1\). This successive approximation method can be repeated to add additional terms proportional to \(\cos n(\omega t-\delta )\) where \(n\) is an integer with \(n\geq 3\). Thus the nonlinearity introduces progressively weaker \(n\)-fold harmonics to the solution. This successive approximation approach is viable only when the admixture coefficient \(\varepsilon <1.\) Note that these harmonics are integer multiples of \(\omega\), thus the steady-state response is identical for each full period even though the state space contours deviate from an elliptical shape.

## Onset of complication

Figure \(\PageIndex{1}\) shows that for \(\gamma =1.05\) the drive strength is sufficiently strong to cause the transient solution for the pendulum to rotate through two complete cycles before settling down to a single steady-state attractor solution at the drive frequency. However, this attractor solution is shifted two complete rotations relative to the initial condition. The state space diagram clearly shows the rolling motion of the transient solution for the first two periods prior to the system settling down to a single steady-state attractor. The successive approximation approach completely fails at this coupling strength since \(\theta\) oscillates through large values that are multiples of \(\pi .\)

Figure \(\PageIndex{1}\) shows that for drive strength \(\gamma =1.078\) the motion evolves to a much more complicated periodic motion with a period that is three times the period of the driving force. Moreover the amplitude exceeds \( 2\pi\) corresponding to the pendulum oscillating over top dead center with the centroid of the motion offset by \(3\pi\) from the initial condition. Both the state-space diagram, and the time dependence of the motion, illustrate the complexity of this motion which depends sensitively on the magnitude of the drive strength \(\gamma ,\) in addition to the initial conditions, \((\theta (0),\omega (0))\) and damping factor \(Q\) as is shown in Figure \(\PageIndex{2}\)

## Period doubling and bifurcation

For drive strength \(\gamma =1.078,\) with the initial condition \(\left( \theta (0),\omega \left( 0\right) \right) =\left( 0,0\right) ,\) the system exhibits a regular motion with a period that is three times the drive period. In contrast, if the initial condition is \([\theta (0)=-\frac{\pi }{2} ,\omega \left( 0\right) =0]\) then, as shown in Figure \(\PageIndex{2}\), the steady-state solution has the drive frequency with no offset in \(\theta\), that is, it exhibits period-one oscillation. This appearance of two separate and very different attractors for \(\gamma =1.078,\) using different initial conditions, is called **bifurcation**.

An additional feature of the system response for \(\gamma =1.078\) is that changing the initial conditions to \([\theta (0)=-\frac{\pi }{2},\omega \left( 0\right) =0]\) shows that the amplitude of the even and odd periods of oscillation differ slightly in shape and amplitude, that is, the system really has period-two oscillation. This period-two motion, i.e. **period doubling**, is clearly illustrated by the state space diagram in that, although the motion still is dominated by period-one oscillations, the even and odd cycles are slightly displaced. Thus, for different initial conditions, the system for \(\gamma =1.078\) bifurcates into either of two attractors that have very different waveforms, one of which exhibits period doubling.

The period doubling exhibited for \(\gamma =1.078,\) is followed by a second period doubling when \(\gamma =1.081\) as shown in Figure \(\PageIndex{2}\). With increase in drive strength this period doubling keeps increasing in binary multiples to period \(8\), \(16\), \(32\), \(64\) etc. Numerically it is found that the threshold for period doubling is \(\gamma _{1}=1.0663,\) from two to four occurs at \(\gamma _{2}=1.0793\) etc. Feigenbaum showed that this cascade increases with increase in drive strength according to the relation that obeys

\[(\gamma _{n+1}-\gamma _{n})\simeq \frac{1}{\delta }(\gamma _{n}-\gamma _{n-1})\]

where \(\delta =4.6692016\), \(\delta\) is called a Feigenbaum number. As \( n\rightarrow \infty \ \)this cascading sequence goes to a limit \(\gamma _{c}\) where \[\gamma _{c}=1.0829\]

## Rolling motion

It was shown that for \(\gamma >1.05\) the transient solution causes the pendulum to have angle excursions exceeding \(2\pi\), that is, the system rolls over top dead center. For drive strengths in the range \(1.3<\gamma <1.4,\) the steady-state solution for the system undergoes continuous rolling motion as illustrated in Figure \(\PageIndex{3}\). The time dependence for the angle exhibits a periodic oscillatory motion superimposed upon a monotonic rolling motion, whereas the time dependence of the angular frequency \(\omega =\frac{ d\theta }{dt}\) is periodic. The state space plots for rolling motion corresponds to a chain of loops with a spacing of \(2\pi\) between each loop. The state space diagram for rolling motion is more compactly presented if the origin is shifted by \(2\pi\) per revolution to keep the plot within bounds as illustrated in Figure \(\PageIndex{3c}\).

## Onset of chaos

When the drive strength is increased to \(\gamma =1.105,\) then the system does not approach a unique attractor as illustrated by Figure \(\PageIndex{4 left}\) which shows state space orbits for cycles \(25-200\). Note that these orbits do not repeat implying the onset of chaos. For drive strengths greater than \( \gamma _{c}=1.0829\) the driven damped plane pendulum starts to exhibit chaotic behavior. The onset of chaotic motion is illustrated by making a \(3\)-dimensional plot which combines the time coordinate with the state-space coordinates as illustrated in Figure \(\PageIndex{4 right}\). This plot shows \(16\) trajectories starting at different initial values in the range \(-0.15<\theta <0.15\) for \(\gamma =1.168\). Some solutions are erratic in that, while trying to oscillate at the drive frequency, they never settle down to a steady periodic motion which is characteristic of chaotic motion. Figure \(\PageIndex{4 right}\) illustrates the considerable sensitivity of the motion to the initial conditions. That is, this deterministic system can exhibit either order, or chaos, dependent on miniscule differences in initial conditions.

^{1}A similar approach is used by the book *"Chaotic Dynamics"* by Baker and Gollub[Bak96].