Taylor Chapter 12: The Driven Damped Pendulum

This notebook takes the place of in-person instruction for the  Nov 28 and Nov 30, 2023 class in PHYS 3101 Analytical Mechanics. It follows Chapter 12 in Taylor “Classical Mechanics.” It goes with homework exercises that should be done in Mathematica.

You should follow along with the text in Taylor, but here are the basics of what’s going on with this chapter. Most of the differential equations we study are linear. Among other things, this means that the phenomena they describe are insensitive to the initial conditions. That is, if you make a small change in the initial conditions, you get a small change in the behavior after long times. This is not true, in general, for nonlinear differential equations. In fact, many other types of behavior that you’ve come to take for granted are not true for phenomena governed by nonlinear equations.

This will all be demonstrated using a particular system, namely the damped, driven pendulum. This problem has just one degree of freedom, and has a nonlinearity (in the sine function) that you are very familiar with. We’ll start with the linearized problem, which you get with no driving term and only small angles, and show that the motion is just what you expect. We’ll then add driving, and watch the system behave in peculiar ways, at first with a period equal to the driving term, but then seeing the period double, triple, quadruple, and so forth, until we get to a point where there is no more periodic motion. This is “chaos” and we’ll quantify that by studying sensitivity to the initial conditions.

It is best for you to go through this notebook step by step. Feel free to modify things as you go along. I think that working on it with others is a smart way to learn this, you can all discuss what’s going on together at each step.

Let’s get started.

Getting started

Clear the memory and set up some things we’ll be using over and over again

Remove["Global`*"]

The fundamental differential equation is given by Taylor Equation (12.11)

eqnBase=ϕ''[t]+2β ϕ'[t]+ω0^2Sin[ϕ[t]]==γ ω0^2Cos[ω t];

The examples in Taylor make use of specific fundamental and driving periods and damping parameter

eqnExamples=ReplaceAll[ReplaceAll[ReplaceAll[eqnBase,
β->ω0/4],ω0->1.5ω],ω->2Pi];

We also define some parameters that will be useful for later.

tMax=40;len=1;

Sanity Check: No Damping and No Driving

Find and plot the solutions for small and large initial angles, with the pendulum starting from rest. Set the natural frequency ω0 equal to 2π so that the period for small oscillations equals one.

eqnSanity=ReplaceAll[ReplaceAll[ReplaceAll[eqnBase,
β->0],{γ->0,ω0->2Pi}],ω->2Pi];

sol=NDSolve[{eqnSanity,ϕ[0]==Pi/30,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕFreeLow=ϕ[t]/.sol;

Plot[ϕFreeLow,{t,0,tMax/10}]

This looks right. The motion looks sinusoidal with period equal to unity and the amplitude looks like π/30.

Now repeat for a large amplitude.

sol=NDSolve[{eqnSanity,ϕ[0]==0.99Pi,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕFreeHigh=ϕ[t]/.sol;

Plot[ϕFreeHigh,{t,0,tMax/10}]

This also looks right. The motion isn’t sinusoidal anymore, and the period is much longer than unity, but the amplitude is close to π

Animations are useful tools for visualization

Build the animation setup for the pendulum, and run the two sanity checks

First run it for the small starting angle solution

xBob=len Sin[ϕFreeLow];yBob=-len Cos[ϕFreeLow];pBob={xBob[[1]],yBob[[1]]};

Animate[
Show[
Plot[0,{x,-len,len},AspectRatio->1],
Graphics[{Thick,Line[{{0,0},pBob/.t->time}]}],
Graphics[{Black,Disk[pBob/.t->time,len/20]}]
],{time,0,5},AnimationRate->0.5,AnimationRunning->False]

To run the large angle solution, just replace the bob positions and rerun the animation

xBob=len Sin[ϕFreeHigh];yBob=-len Cos[ϕFreeHigh];pBob={xBob[[1]],yBob[[1]]};

One more sanity check: Add some weak damping

Use a damping parameter that is less than one, but not too small

eqnDamped=ReplaceAll[ReplaceAll[ReplaceAll[eqnBase,
β->0.2],{γ->0,ω0->2Pi}],ω->2Pi];

Find the solution starting from rest a ϕ=π/2.

sol=NDSolve[{eqnDamped,ϕ[0]==Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕDamped=ϕ[t]/.sol;

Plot and animate

Plot[ϕDamped,{t,0,10},PlotStyle->Thick]

xBob=len Sin[ϕDamped];yBob=-len Cos[ϕDamped];pBob={xBob[[1]],yBob[[1]]};

(Once again, go back to the animation slide to run it.)

Now drive it with a weak amplitude (F0/mg=0.2)

Set the parameters equal to those used to make Figure 12.2 in Taylor

sol=NDSolve[{eqnExamples/.γ->0.2,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p2=ϕ[t]/.sol;

Plot the solution to confirm that it looks like Figure 12.2

Plot[ϕTaylorFig12p2,{t,0,6},PlotStyle->Thick]

The motion looks to be sinusoidal, with a period given by the driving frequency. This is just what you’d expect from a linear (that is, a simple harmonic) oscillator.

Animate the result. It’s nice to see the bob start from “zero” and then rise up due to the forcing

xBob=len Sin[ϕTaylorFig12p2];yBob=-len Cos[ϕTaylorFig12p2];pBob={xBob[[1]],yBob[[1]]};

(Once again, go back to the animation slide to run it.)

Now a larger amplitude (F0/mg=0.9)

Set the parameters equal to those used to make Figure 12.3 in Taylor

sol=NDSolve[{eqnExamples/.γ->0.9,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p3=ϕ[t]/.sol;

Plot[ϕTaylorFig12p3,{t,0,6},PlotStyle->Thick]

Following is Figure 12.3(b). I just eyeballed the parameters of the sine curve to look like the textbook. The point is, of course, that the true motion is no longer a sine curve. That is, the nonlinear character of the forced damped pendulum is starting to show up.

Show[
Plot[ϕTaylorFig12p3,{t,4.9,6.1},PlotStyle->Thick],
Plot[2.75Sin[2 Pi t-0.175],{t,4.9,6.1},PlotStyle->{Red,Thick,Dashed}]
]

Nevertheless, the motion still looks like it has the period of the driving frequency. (We’ll check this in a moment.)

The animation also seems to show oddness that we can blame on nonlinearity

xBob=len Sin[ϕTaylorFig12p3];yBob=-len Cos[ϕTaylorFig12p3];pBob={xBob[[1]],yBob[[1]]};

(Once again, go back to the animation slide to run it.)

Confirm that the motion at long times is periodic

Just because the motion is nonlinear doesn’t mean that it isn’t periodic. We can check by looking at the value of ϕ at different times, well into the future after the transients have died out, that are separated by the driving period.

ϕvals=(ϕTaylorFig12p3/.t->{31,32,33,34,35,36,37,38,39,40})[[1]]

The values look to be the same. That is, for times separated by one unit of the driving period, the value of ϕ(t) is the same. Just to check, look at the differences between successive elements in the list.

Differences[ϕvals]

Well, they’re not exactly zero, but they are very small. We can blame this on the limits of numerical precision for the solution of the differential equation

Now, an even larger driving amplitude γ = F0/mg=1.06 > 1

This is Taylor Figure 12.4. The driving amplitude is larger than one, and the behavior is drastically different

sol=NDSolve[{eqnExamples/.γ->1.06,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p4=ϕ[t]/.sol;

Plot[ϕTaylorFig12p4,{t,0,15},PlotStyle->Thick]

The pendulum flips over a few times before finally settling into something that looks periodic. (We’ll check that shortly.) The animation is cool. I’m re-copying the Animation code here, so that it goes to longer times (and so that you don’t have to flip all the back to that earlier slide.)

xBob=len Sin[ϕTaylorFig12p4];yBob=-len Cos[ϕTaylorFig12p4];pBob={xBob[[1]],yBob[[1]]};

Animate[Show[Plot[0,{x,-len,len},AspectRatio->1],
Graphics[{Thick,Line[{{0,0},pBob/.t->time}]}],
Graphics[{Black,Disk[pBob/.t->time,len/20]}]
],{time,0,15},AnimationRate->0.5,AnimationRunning->False]

Crazy motion, but it still has period (pretty close to) one

Once again, look at the values at times separated by one period and check thir differences.

ϕvals=(ϕTaylorFig12p4/.t->{31,32,33,34,35,36,37,38,39,40})[[1]]

Differences[ϕvals]

The deviations from unity are not quite as small, but still very small. It might be worth going even further out in time to see if they get closer, but I’ll leave that up to you.

In fact, if you read Footnote 8 on page 468 of Taylor, you’ll see that he says you need to go up to period 46 before things equilibrate to give the same value for another significant digit. If you want to check this, you’ll need to change the value of tMax on the second slide.

Now, use γ = F0/mg=1.073

This is Taylor Figure 12.5. Something very interesting happens here. The motion undergoes “bifurcation.”

sol=NDSolve[{eqnExamples/.γ->1.073,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p5=ϕ[t]/.sol;

We need to go pretty far out in time before we see something that even looks periodic

Plot[ϕTaylorFig12p5,{t,0,30},PlotStyle->Thick]

However, the period has in fact doubled!! You can see this easily if you vertically expand the later times, as in Taylor Figure 12.5(b). (In the next slide, we will check numerically that the period is two.)

Plot[ϕTaylorFig12p5,{t,20,30},PlotStyle->Thick,PlotRange->{-8.7,-7.5}]

Watch the animation and see if you can spot the double period motion. I’m reproducing the animation code so that it will run out to the full extent of the times for which we solved the differential equation.

xBob=len Sin[ϕ[t]]/.sol;yBob=-len Cos[ϕ[t]]/.sol;pBob={xBob[[1]],yBob[[1]]};

Animate[Show[Plot[0,{x,-len,len},AspectRatio->1],
Graphics[{Thick,Line[{{0,0},pBob/.t->time}]}],
Graphics[{Black,Disk[pBob/.t->time,len/20]}]
],{time,0,tMax},
AnimationRate->0.5,AnimationRunning->False,AppearanceElements->All]

Numerical demonstration that the period is two

Exactly the same thing we’ve done before, except...

ϕvals=(ϕTaylorFig12p5/.t->{31,32,33,34,35,36,37,38,39,40})[[1]]

... you see that the repeat is every other value. That is, the period is two, not one. You can see this because the differences toggle back and forth between plus and minus the same number.

Differences[ϕvals]

A slight increase to γ = F0/mg=1.077

This is Taylor Figure 12.6. The period triples!

sol=NDSolve[{eqnExamples/.γ->1.077,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p6=ϕ[t]/.sol;

Plot[ϕTaylorFig12p6,{t,0,15},PlotStyle->Thick]

Also, it seems to have settled into its periodic state after only a short while. Check that the period is in fact three.

ϕvals=(ϕTaylorFig12p6/.t->{30,31,32,33,34,35,36,37,38})[[1]]

You know enough now to define the bob position and run the animation if you like.

You have a homework problem to reproduce Taylor Figure 12.7. This figure shows that for exactly the same equations that created Figure 12.6, except for ϕ(0)=-π/2 instead of ϕ(0)=0, the motion is completely different.

Taylor’s “Period Doubling Cascade”

This is Taylor Figure 12.8. It is the same Damped DPP we have been running, but with a different initial position. Four different driving strengths are used, and the motions show the period going from one, to two, to four, to eight.

sol=NDSolve[{eqnExamples/.γ->1.06,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p8γ1p06=ϕ[t]/.sol;

sol=NDSolve[{eqnExamples/.γ->1.078,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p8γ1p078=ϕ[t]/.sol;

sol=NDSolve[{eqnExamples/.γ->1.081,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p8γ1p081=ϕ[t]/.sol;

sol=NDSolve[{eqnExamples/.γ->1.0826,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p8γ1p0826=ϕ[t]/.sol;

Grid[
{{Plot[ϕTaylorFig12p8γ1p06,{t,0,10},PlotStyle->Thick],
Plot[ϕTaylorFig12p8γ1p06,{t,28,40},PlotRange->{-2.5,-1.7}]},
{Plot[ϕTaylorFig12p8γ1p078,{t,0,10},PlotStyle->Thick],
Plot[ϕTaylorFig12p8γ1p078,{t,28,40},PlotRange->{-2.5,-1.7}]},
{Plot[ϕTaylorFig12p8γ1p081,{t,0,10},PlotStyle->Thick],
Plot[ϕTaylorFig12p8γ1p081,{t,28,40},PlotRange->{-2.5,-1.7}]},
{Plot[ϕTaylorFig12p8γ1p0826,{t,0,10},PlotStyle->Thick],
Plot[ϕTaylorFig12p8γ1p0826,{t,28,40},PlotRange->{-2.5,-1.7}]}}]

Check to find out if the motion is still periodic

I’m going to output the numbers this time with a fixed number of significant digits.

NumberForm[ϕTaylorFig12p8γ1p06/.t->{31,32,33,34,35,36,37,38},4]

NumberForm[ϕTaylorFig12p8γ1p078/.t->{31,32,33,34,35,36,37,38},4]

NumberForm[ϕTaylorFig12p8γ1p081/.t->{31,32,33,34,35,36,37,38},4]

NumberForm[ϕTaylorFig12p8γ1p0826/.t->{25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40},4]

This looks right. By the time we get to the period-eight result, though, the precision looks a little dicey. For example “-0.08289” isn’t exactly the same as “-0.08278”, but this is probably due to the limits of numerical precision. It’s worth checking out what happens if we solved the equation for a much longer time.

One more: γ=F0/mg=1.09

You should read in Taylor about the “Feigenbaum number” and universailty in nonlinear systems. The period-doubling cascade we just went through used increasing values of the driving amplitude, but those increases were getting very tiny. In fact, they converge on the value 1.10829. For higher driving strengths, the motion is “chaotic”. We’ll be more careful about what we mean by chaos in a moment, but let’s look at the motion for one more driving strength, just above the critical value.

sol=NDSolve[{eqnExamples/.γ->1.09,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕγ1p09=ϕ[t]/.sol;

Plot[ϕγ1p09,{t,0,tMax},PlotStyle->Thick]

There doesn’t seem to be anything that looks periodic here. Let’s check the numbers.

NumberForm[ϕγ1p09/.t->{25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40},4]

None of those numbers repeat. It looks like the motion isn’t periodic at all. In fact, it is “chaotic”. Let’s see how we can quantify that statement.

Quantifying chaos

Chaos is characterized by the sensitivity to initial conditions. We will see for the DDP examples we ran above, that a tiny change in the initial conditions for a non-chaotic system lead, at least eventually, to a small change in the motion at some time. This is not true for a chaotic system. A tiny change in the initial conditions, in this case, leads to a difference in the motion which grows forever.

To illustrate this, let’s go back to an “almost linear” case and see what happens when we make a small change in the initial conditions. Use a small driving amplitude with our standard example equation, and in one case let the pendulum start at ϕ=0 and in the other case ϕ=0.1.

ϕ0=ϕ[t]/.NDSolve[{eqnExamples/.γ->0.1,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕ1=ϕ[t]/.NDSolve[{eqnExamples/.γ->0.1,ϕ[0]==0.1,ϕ'[0]==0},ϕ,{t,0,tMax}];

Plot[{ϕ0,ϕ1},{t,0,5}]

It sure looks like that, after just a couple of periods, the motion is the same for the two cases. In fact, let’s plot the difference as a function of time.

Δϕ=ϕ0-ϕ1;
Plot[Δϕ,{t,0,5}]

Yes. the difference in the two motions goes to zero after several periods.

A better way to plot the difference

See Taylor Figure 12.11.

It will turn out to be more useful to plot the difference on a log plot. This is a problem if the difference is negative, so take the absolute value. This leads to other screwiness, namely the “cusp” behavior when the difference goes negative, but don’t worry about that. It will not be an issue when we get into the serious nonlinear regime.

LogPlot[Abs[Δϕ],{t,0,10}]

We may be getting into trouble with precision after t=7 or so, but it seems clear that the differences are falling exponentially. In fact, for a linear system with damping coefficient β, we expect the difference to fall like Exp[-βt]. Let’s compare this with the curve above, using the same value of β that we’ve been using or our examples.

LogPlot[{Abs[Δϕ],Exp[-β t]/.β->1.5 2 Pi/4},{t,0,10}]

It looks like the same is true here. The differences fall like the damping coefficient, just as you’d expect for a linear system.

Now try γ=F0/mg=1.07

Now move towards Taylor Figure 12.12, for γ=1.07 with the same initial conditions as in Figure 12.11. This is very nonlinear, but before the period doubles.

ϕ0=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.07,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕ1=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.07,ϕ[0]==0.1,ϕ'[0]==0},ϕ,{t,0,tMax}];

Plot[{ϕ0,ϕ1},{t,0,tMax}]

This is not going to reproduce Figure 12.12! Maybe there is an additional sensitivity due to numerical solutions on my computer with my version of Mathematica, compared to what Taylor used to make his figure.

Let’s try again with a smaller difference in initial angles.

ϕ0=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.07,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕ1=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.07,ϕ[0]==0.05,ϕ'[0]==0},ϕ,{t,0,tMax}];

Plot[{ϕ0,ϕ1},{t,0,tMax}]

That’s better. Now plot the difference.

Δϕ=ϕ0-ϕ1;
LogPlot[Abs[Δϕ],{t,0,tMax}]

Well, that still looks weird, still nothing like Figure 12.12. Again I blame this on the precision. We could play around with the relevant parameters in NDSolve, or I can just find a different example that shows what Figure 12.12 is trying to show. I found the following, just by poking around. Note that γ=1.05 is also before bifurcation, just not as strongly nonlinear as γ=1.07.

ϕ0=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.05,ϕ[0]==0.2,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕ1=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.05,ϕ[0]==0.3,ϕ'[0]==0},ϕ,{t,0,tMax}];

Plot[{ϕ0,ϕ1},{t,0,tMax}]

Δϕ=ϕ0-ϕ1;
LogPlot[{Abs[Δϕ],Exp[-β t]/.β->1.5 2 Pi/4},{t,0,20},
PlotRange->{10^(-9),10}]

This makes our point. In this very nonlinear, but still periodic, regime for γ, the small change in initial conditions still leads to a small difference in the change of behavior after long times.. However, this time, the falloff is much slower than it is for the almost-linear case, namely Taylor Figure 12.11.

Now go into the “chaotic” region

Our task now is to reproduceTaylor Figure 12.13, for γ=1.105, which is well into the chaotic regime. The motion is shown in Taylor Figure 12.10, so let’s plot that first.

sol=NDSolve[{eqnExamples/.γ->1.105,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕTaylorFig12p11=ϕ[t]/.sol;

Plot[ϕTaylorFig12p11,{t,0,30}]

Well it certainly looks “chaotic”, that is, no sign of any periodic motion, but this does not reproduce the shape of Figure 12.11 very well. I’m again happy to blame this on the idiosyncrasies of numerical solutions for these nonlinear differential equations, but the essence featueres are the same.

In fact, the extreme sensitivity to initial conditions is apparent without having to plot the logarithm of the different. Following are two plots with an extremely small difference in the initial conditions, namely one part in 10,000.

ϕ0=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.105,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax}];
ϕ1=ϕ[t]/.NDSolve[{eqnExamples/.γ->1.105,ϕ[0]==0.0001,ϕ'[0]==0},ϕ,{t,0,tMax}];

Plot[{ϕ0,ϕ1},{t,0,30}]

The two motions look to be the same for the first six or so periods. But then the motions clearly diverge. This is chaos. A very tiny change in the initial conditions leads to a huge variation in the behavior after  long times. This is popularized as the “butterfly effect.” Chaos was first observed in weather pattern simulations, and people like to say that if a butterfly flaps its wings in Africa and you don’t take it into account (as a small change in the initial condition) then the weather in the United States will be very different than what you predict.

Plot the difference

Plotting the (logarithm) of the difference is a good way to see the “butterfly effect” quantitatively.

Δϕ=ϕ0-ϕ1;

LogPlot[Abs[Δϕ],{t,0,15}]

This is, effectively, Figure 12.13 in Taylor, although the differences in numerical precision between this version and Taylor’s are apparent. It grows a bit faster than the plot in Taylor, and levels off sooner. The leveling off, as explained in Taylor, is easy to understand, because the biggest possible difference in ϕ is 2π.

The difference in the two solutions grows exponentially in time. (This is clear because the the log of the difference increases linearly up through six or seven periods.) In other words, even an infinitesimal change in the initial conditions will yield a finite change in the motion. This is the signature of chaos. This is of course entirely different from the case in Figure 12.11 where the difference falls exponentially with time. If we write |Δϕ|=K Exp[λt], where λ is called the Liapunov Coefficient, and K is a positive constant, then we say that chaos is characterized by a positive Liapunov coefficient.

There is plenty more that can be learned from the Damped Driven Pendulum, but we’ll stop here. Your homework assignment will help you appreciate what’s been covered up to here. Of course, feel free to dive into the rest of chapter 12 if you like!