PHYS 2063 Wave Physics:
Nonlinear Dynamics and Chaos in the Damped Driven Pendulum

This discussion follows Taylor, Classical Mechanics, Chapter 12

Remove["Global`*"]

Set up the differential equation of motion

Start with Newton’s Second Law in terms of torques on the pendulum bob. Note that the moment of inertia is I=ml^2 and the damping force is -bv=-blϕ’.

eq=m l^2 ϕ''[t]==-m g l Sin[ϕ[t]]-b l^2ϕ'[t]+l F0 Cos[ω t]

Turn this into a form using Taylor’s parameterization. Note that a “natural” value for the acceleration of the bob is g=lω0^2, so γ plays the role of a dimensionless “strength” for the driving force.

eqreps={g->ω0^2l,b->2β m,F0->γ ω0^2l m};
eqDDP=Simplify[eq/.eqreps,l>0&&m>0]

Set up the standard equation that we’ll be solving. This picks a relatively large damping factor, but still sub critical, and a driving frequency that is 50% higher than the natural frequency of the pendulum. Notice that the driving period is defined to be one! We will use this fact to test later for the periodicity of the solutions.

eqn=((eqDDP/.{β->ω0/4})/.{ω0->1.5ω})/.ω->2 Pi;
ics={ϕ[0]==0,ϕ'[0]==0};

Checking the solution to the equation of motion

Solve, plot, and animate the solution with a small value of γ (Taylor Fig.12.2)

tMax=10;
sol=NDSolve[{eqn/.γ->0.2,ics},ϕ,{t,0,tMax}];

Plot[ϕ[t]/.sol,{t,0,6},PlotRange->{-0.5,0.5}]

len=1;
Bob=len Sin[ϕ[t]]/.sol;
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]

Confirm that the motion has a single period

NumberForm[ϕ[{4,5,6,7,8,9,10}]/.sol,5]

Now use the same equations in the nonlinear regime

Solve, plot, and animate the solution with γ=1.06 (Taylor Fig.12.4)

tMax=40;
sol=NDSolve[{eqn/.γ->1.06,ics},ϕ,{t,0,tMax}];

Plot[ϕ[t]/.sol,{t,0,15},PlotRange->{-Pi,5 Pi}]

len=1;
Bob=len Sin[ϕ[t]]/.sol;
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,15},AnimationRate->0.5,AnimationRunning->False]

Show that the motion still has a single period

NumberForm[ϕ[{34,35,36,37,38,39,40}]/.sol,5]

Next, drive it a little harder

Repeat for γ=1.073 and demonstrate frequncy doubling (Taylor Fig.12.5)

tMax=40;
sol=NDSolve[{eqn/.γ->1.073,ics},ϕ,{t,0,tMax}];

Plot[ϕ[t]/.sol,{t,0,30},PlotRange->{-3Pi,5 Pi}]

Plot[ϕ[t]/.sol,{t,20,30},PlotRange->{-8.8,-7.5}]

len=1;
Bob=len Sin[ϕ[t]]/.sol;
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,25},AnimationRate->0.5,AnimationRunning->False]

NumberForm[ϕ[{30,31,32,33,34,35}]/.sol,5]

And now, a little harder still

Repeat for γ=1.077 and demonstrate frequncy tripling (Taylor Fig.12.6)

tMax=40;
sol=NDSolve[{eqn/.γ->1.077,ics},ϕ,{t,0,tMax}];

Plot[ϕ[t]/.sol,{t,0,15},PlotRange->{-Pi,5 Pi}]

len=1;
Bob=len Sin[ϕ[t]]/.sol;
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,25},AnimationRate->0.5,AnimationRunning->False]

NumberForm[ϕ[{30,31,32,33,34,35,36,37,38}]/.sol,6]

And now, a little bit harder, showing non-periodic motion

Repeat for γ=1.105 and the onset to chaos (Taylor Fig.12.10, but I’m sticking with ϕ(0)=0.)

tMax=40;
sol=NDSolve[{eqn/.γ->1.105,ics},ϕ,{t,0,tMax}];

Plot[ϕ[t]/.sol,{t,0,30},PlotRange->{-1.2Pi, 1.2Pi}]

len=1;
Bob=len Sin[ϕ[t]]/.sol;
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,25},AnimationRate->0.5,AnimationRunning->False]

NumberForm[ϕ[{30,31,32,33,34,35,36,37,38,39}]/.sol,5]

Nonperiodic! But is this “chaos”? What do we mean by “chaos” anyway?

Summary of what we’ve seen so far

Construct a series of solutions for different values of γ

tMax=40;
eqnNonLin=eqn/.γ->1.06;
solnN1p06=NDSolve[{eqnNonLin,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
eqnNonLin=eqn/.γ->1.078;
solnN1p078=NDSolve[{eqnNonLin,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
eqnNonLin=eqn/.γ->1.081;
solnN1p081=NDSolve[{eqnNonLin,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];
eqnNonLin=eqn/.γ->1.0826;
solnN1p0826=NDSolve[{eqnNonLin,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];

From the plots, you see a progression from period=1, to period=2, to period=4, to period=8. This phenomenon is called a Period-Doubling Cascade.

Table[
{{Plot[ϕ[t]/.solnN1p06,{t,0,10},PlotStyle->Thick],
Plot[ϕ[t]/.solnN1p06,{t,28,40},PlotRange->{-2.5,-1.7}]},
{Plot[ϕ[t]/.solnN1p078,{t,0,10},PlotStyle->Thick],
Plot[ϕ[t]/.solnN1p078,{t,28,40},PlotRange->{-2.5,-1.7}]},
{Plot[ϕ[t]/.solnN1p081,{t,0,10},PlotStyle->Thick],
Plot[ϕ[t]/.solnN1p081,{t,28,40},PlotRange->{-2.5,-1.7}]},
{Plot[ϕ[t]/.solnN1p0826,{t,0,10},PlotStyle->Thick],
Plot[ϕ[t]/.solnN1p0826,{t,28,40},PlotRange->{-2.5,-1.7}]}}]

One more for now: γ=1.09

Pushing into the aperiodic regime again

tMax=40;eqnNonLin=eqn/.γ->1.09;

solnNonLin=NDSolve[{eqnNonLin,ϕ[0]==-Pi/2,ϕ'[0]==0},ϕ,{t,0,tMax}];

Plot[ϕ[t]/.solnNonLin,{t,0,tMax},PlotStyle->Thick]

The motion is no longer periodic!

NumberForm[
ϕ[{25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40}]/.solnNonLin,6]

The Bifurcation Diagram

Watch the different periods appear as you increase γ. Be careful of the numerical precision.

NumberForm[ϕ[{31,32,33,34,35,36,37,38}]/.solnN1p06,4]

NumberForm[ϕ[{31,32,33,34,35,36,37,38}]/.solnN1p078,4]

NumberForm[ϕ[{31,32,33,34,35,36,37,38}]/.solnN1p081,4]

NumberForm[ϕ[{31,32,33,34,35,36,37,38}]/.solnN1p0826,4]

See Taylor Figure 12.17 for the full bifurcation diagram

Quantifying Chaos: Sensitivity to Initial Conditions

Let’s go back and look at the “almost linear” case

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

Plot[ϕ0,{t,0,10}]

How different are the two solutions?

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

A log plot is a better way to see what is going on

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

In other words, the transients die away, and at some point we hit the limits of precision.

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

Now repeat for a larger value of γ

Choose γ=1.05. This is still motion with a single period, but very nonlinear nevertheless.

ϕ0=ϕ[t]/.NDSolve[{eqn/.γ->1.05,ϕ[0]==0.2,ϕ'[0]==0},ϕ,{t,0,tMax},
MaxSteps->∞,AccuracyGoal->30];
ϕ1=ϕ[t]/.NDSolve[{eqn/.γ->1.05,ϕ[0]==0.3,ϕ'[0]==0},ϕ,{t,0,tMax},
MaxSteps->∞,AccuracyGoal->30];

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

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

The difference rises at first, but then falls off, albeit more slowly than it did for γ=0.1

Now see what happens in the “chaotic” regime

Choose γ=1.105, and a very tiny difference in the initial conditions

ϕ0=ϕ[t]/.NDSolve[{eqn/.γ->1.105,ϕ[0]==0,ϕ'[0]==0},ϕ,{t,0,tMax},
MaxSteps->∞,AccuracyGoal->30];
ϕ1=ϕ[t]/.NDSolve[{eqn/.γ->1.105,ϕ[0]==0.0001,ϕ'[0]==0},ϕ,{t,0,tMax},
MaxSteps->∞,AccuracyGoal->30];

Plot the two motions. Even though the ICs are almost identical, the two solutions diverge after several periods and then never approach each other again. Sometimes this is referred to as the “butterfly effect”.

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

Plot the difference to demonstrate this

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

This is called “chaos”. There is an exponential sensitivity to initial conditions.