Driven Pendulum

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

A driven pendulum with linear damping satisfies an equation of motion of the form

(((∂^2) θ)/(∂(t^2)))+γ ((∂θ)/(∂t))+ ((Subscript[ω, 0])^2) Sin[θ]  == f Cos[ω t]

where f and ω are the amplitude and frequency of the driving force and γ is the damping coefficient.  The natural frequency for free small-amplitude oscillations is (Subscript[ω, 0]) = (Sqrt[g/ℓ]), but it is convenient to express frequency in units of (Subscript[ω, 0]) by setting  (Subscript[ω, 0]) -> 1 and rescaling the time unit accordingly.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

$DefaultFont={"Times", 12}; 
$TextStyle={FontFamily->"Times", FontSize->12,FontSlant->"Italic"};

Needs["Utilities`Notation`"];
Needs["Utilities`FilterOptions`"]

Symbolize[(Subscript[ω, 0])]; Symbolize[(Subscript[θ, 0])];Symbolize[\(θ\& . \_0\)];

Implementation

In this section we construct functions which define and solve the pendulum equation given its parameters and initial conditions.

\(\(\(pendulumEq::usage\  = \ \*""pendulumEq[γ,f,ω] describes the motion of a driven pendulum with linear damping, \n (θ^″) + γ (θ^′)+Sin[θ]==f Cos[ω t], \n where f and ω are the driving force and frequency and γ is the damping coefficient."";\)\(
\)
    \)\)
pendulumEq[γ_, 
        f_, ω_] :=  (θ^″)[
            t]  +  γ (θ^′)[t]  +  
          Sin[θ[t]] ==  f Cos[ω t]

pendulumSol::usage = "pendulumSol[γ,f,ω,(Subscript[θ, 0]),\!\(
\*TagBox[\(θ\& . \_0\),
NotationBoxTag,
Editable->True]\),tmin,tmax,opts] uses NDSolve to produce a numerical solution to the pendulum equation given initial conditions θ[0]==(Subscript[θ, 0]) and (θ^′)[0]==\!\(
\*TagBox[
TagBox[\(θ\& . \_0\),
NotationBoxTag,
Editable->True],
NotationBoxTag,
Editable->True]\) for {t,tmin,tmax}.  See pendulumEq for other parameter definitions.  The optional parameters opts are passed to NDSolve.";
		
\(pendulumSol[γ_, 
        f_, ω_, \(θ\_0\) : _, \(θ\& . \_0\) : _, tmin_, 
        tmax_, opts___]\  := \(NDSolve[{pendulumEq[γ, 
              f, ω], θ[
                0] == θ\_0, \(θ\^′\)[
                0] == θ\& . \_0}, θ, {t, tmin, tmax}, 
          opts]\)[[1]]\)

pendulumPlot::usage = "pendulumPlot[γ,f,ω,(Subscript[θ, 0]),\!\(
\*TagBox[\(θ\& . \_0\),
NotationBoxTag,
Editable->True]\),tmin,tmax,opts] produces a phase diagram for the pendulum equation.  The optional parameters opts can regulate either NDSolve or ParametricPlot as appropriate.  See pendulumEq and pendulumSol for parameter definitions.";

\(pendulumPlot[γ_, 
        f_, ω_, \(θ\_0\) : _, \(θ\& . \_0\) : _, tmin_, 
        tmax_, opts___]\  := \ 
      With[{optPlot\  = \ FilterOptions[ParametricPlot, opts]\ , 
          optND\  = 
            FilterOptions[NDSolve, opts]}, 
ParametricPlot[
          Evaluate[{θ[t], \(θ\^′\)[t]} /. 
              pendulumSol[γ, 
                f, ω, θ\_0, θ\& . \_0, tmin, tmax, 
                optND]], {t, tmin, tmax}, \ optPlot, 
          AspectRatio -> Automatic, Frame -> True, 
          FrameLabel -> {"position", \ "velocity"}, 
          PlotLabel -> "Pendulum Equation"]]\)

Examples

Free swing without damping

For small amplitudes we find the classic circular phase diagram characteristic of simple harmonic motion.  For larger amplitudes the nonlinearity of the restoring force distorts the trajectory and increases the period.

plots  =  
        Table[pendulumPlot[0, 0, 0, (Subscript[θ, 0]), 0, 0, 15, 
            DisplayFunction -> 
              Identity], {(Subscript[θ, 0]), (π/8), ((7  π)/8), (π/8)}];

Show[plots, DisplayFunction -> $DisplayFunction];

If we initiate the swing at the bottom with finite velocity, we find that when the total energy is negative the pendulum swings back and forth over a limited range, but for positive energy it orbits the pivot with monotonically increasing angle and a limited band of velocities.  The critical condition in which the bob just reaches the top and stops is represented by the separatrix which divides the phase space into these two distinct regions.  The inverted equilibrium is unstable, of course.

\(\(plots\  = \ 
        Table[pendulumPlot[0, 0, 0, 0, θ\& . \_0, \(-5\), 5, 
            PlotRange -> {{\(-2\) π, 2  π}, Automatic}, 
            DisplayFunction -> Identity], {θ\& . \_0, \(-3\), 3, 
            0.5}];\)\)

Show[plots, DisplayFunction -> $DisplayFunction];

Driven pendulum with damping

For a wide range of initial conditions, the motion eventually settles upon a limiting orbit that depends only on the driving force.

attractor  = 
        pendulumPlot[0.1, 0.1, 0.5, 0, 0, 160, 200, MaxSteps -> 2000, 
          PlotPoints -> 100, PlotRange -> All, 
          PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];

plot1  = 
        pendulumPlot[0.1, 0.1, 0.5, 0, 0, 0, 200, MaxSteps -> 2000, 
          PlotPoints -> 100, PlotRange -> All, 
          DisplayFunction -> Identity];
Show[{plot1, attractor}, 
        DisplayFunction -> $DisplayFunction];

plot2  = 
      pendulumPlot[0.1, 0.1, 0.5, 1, 1, 0, 200, MaxSteps -> 2000, 
        PlotPoints -> 100, 
        PlotRange -> {{-0.2, 0.2}, {-0.2, 0.2}}, 
        DisplayFunction -> Identity]; 
    Show[{plot2, attractor}, 
      DisplayFunction -> $DisplayFunction];

(((∂^2) θ)/(∂(t^2)))+γ ((∂θ)/(∂t))+  ((Subscript[ω, 0])^2)( 1  +  f Cos[ω t] ) 
            Sin[θ] == 0