v. 5: updated 2019-12-03 to fix changes in Export of animated gifs since 20017
v. 4: updated 2019-11-21
v. 3: updated 2017-01-11

This notebook plots ϕ vs t for the driven damped pendulum.

This is a modification of a Mathematica notebook written by my son, Michael Swift, for a junior level mechanics course in the physics department at Rice University.  They use a slightly different parameters to define the system than “Chaos” by Alligood, Sauer and Yorke.   If ω0=1, then the three remaining parameters, c, ρ, and ω, are the same. 

Even though the driven damped pendulum can be written with 3 dimensionless parameters, this ODE has 4 parameters, c, ρ, ω, and ω0.

 If ω0=1, then the three remaining parameters, c, ρ, and ω, are the same as “Chaos” by Alligood, Sauer and Yorke.  This was the text for MAT 667 in Spring 2017.
 
 Usually, we will fix ω = 2π, and vary the parameters c, ρ, and ω0.  This has the advantage that the driving period is 1.  We also usually use ω0 = 1.5 ω.  This means that the driving frequency is slightly smaller than the natural frequency for small oscillations of the undamped system.

 This leaves us with two parameters, ρ and c to vary.  This is the strength of the driving term (ρ) and the friction constant (c).  We also have the initial conditions 

Use the Insert menu with Insert File Path... to set the place where animations are saved.

SetDirectory["/Users/jimswift/Dropbox/classes/MAT 665 2019.7/DrivenPendulumAnimations/"] (* for Jim's home *)
SetDirectory["C:\\Users\\jws8\\Dropbox\\classes\\MAT 665 2019.7\\DrivenPendulumAnimations\\"]  (* for Jim's office *)

First define a function which draws a pendulum of unit length at an angle phi, along with an arc representing the applied force.

Pendulum[phi_,F_]:=
Show[
ListPolarPlot[{{0,0},{phi-Pi/2,1}},PlotRange->{{-1.1,1.1},{-1.1,1.1}},Joined->True,AspectRatio->Automatic,ImageSize->{250}],
ListPolarPlot[{{phi-Pi/2,1}},PlotStyle->PointSize[.05]],
PolarPlot[0.3,{t,-Pi/2,(F-1)*Pi/2+10^-9},PlotStyle->{Thick,Red}]
]
Pendulum[0.1, .5]

Now define a function that solves the system, with the parameters c, ω0, etc. to be defined later.

runPendulum := Module[{runTime, Dummy},
{runTime, Dummy} = AbsoluteTiming[
grid2Pi[ϕMin_, ϕMax_] := (* puts grid lines at multiples of 2 Pi *)
Table[2 Pi i, {i, Floor[ϕMin/(2Pi)] , Ceiling[ϕMax/(2 Pi)] }];
If[nTransient > 0,
Do[sol=NDSolve[{ϕ''[t]+c * ω0*ϕ'[t]+ω0^2*Sin[ϕ[t]]==ρ*ω0^2*Cos[ω*t],ϕ[0]== ϕ0,ϕ'[0]==v0},ϕ,{t,0,1}];
{ϕ0, v0} = {ϕ[1],ϕ'[1]}/.sol[[1]](*replace ic with final state*)
,{Floor[nTransient]} (* do this nTransient times *)];
Print["After ", nTransient, " driving periods, start plotting with {ϕ0, v0} = ", {ϕ0, v0} ];
];
phivst = {};
Do[
sol= (* take the one and unique solution, Part 1 *)Part[NDSolve[{ϕ''[t]+c * ω0*ϕ'[t]+ω0^2*Sin[ϕ[t]]==ρ*ω0^2*Cos[ω*t],ϕ[0]== ϕ0,ϕ'[0]==v0},ϕ,{t,0,1}], 1];
phivst = Join[phivst, Table[{tt + tau, ϕ[tau]}/.sol, {tau, 0, 1-step, step}]];
{ϕ0, v0} = {ϕ[1],ϕ'[1]}/.sol[[1]] (*replace ic with final state*)
,{tt,0,Floor[tfin-1], 1}];
If[tfin - Floor[tfin] == 0,
AppendTo[phivst, {tfin, ϕ[1]/.sol}];
Print["Final {ϕ, ϕ'} is ", {ϕ[1], ϕ'[1]}/.sol],
(* tfin is not an integer.  Adjust step size *)
timeInterval = tfin - Floor[tfin];
newFrames = Round[timeInterval/step];
sol= (* take the one and unique solution, Part 1 *)Part[NDSolve[{ϕ''[t]+c * ω0*ϕ'[t]+ω0^2*Sin[ϕ[t]]==ρ*ω0^2*Cos[ω*t],ϕ[Floor[tfin]]== ϕ0,ϕ'[Floor[tfin]]==v0},ϕ,{t,Floor[tfin],tfin}], 1];
Do[tt=Floor[tfin] + i timeInterval/newFrames;
AppendTo[phivst, {tt, ϕ[tt]}/.sol],{i, 0, newFrames-1}];
Print["Final {ϕ, ϕ'} is ", {ϕ[tfin], ϕ'[tfin]}/.sol]
];
figRight =ListPlot[phivst, Joined -> True, InterpolationOrder -> 2(* -> 2 gives smoother curves but sometimes has gaps at edge of window *), GridLines ->  {None,grid2Pi}];
If[makeMovie, 
frames=Table[
GraphicsRow[
{Pendulum[phivst[[i, 2]],ρ*Cos[ω*phivst[[i, 1]] ] ],
Show[figRight,Graphics[{Red,PointSize[0.015], Point[phivst[[i]]] }]]}, ImageSize ->{560, 315}
],{i,1,Length[phivst]-1}]; (* t = 0 and t = tfin are same for periodic solutions, so don't include both frames *)
Export[filename, frames, "DisplayDurations" -> 3*step, AnimationRepetitions-> Infinity,ImageSize -> {560, 315}]
(* Each framed is displayed for 3*step seconds, so one driving period (Δt = 1) takes 3 seconds *)
];
];
Print["Took ", runTime, " seconds."];
Show[figRight]
]

Now define the parameters, and run the program.

An animated gif is produced if makeMovie is True.  The pendulum will be animated next to the graph of its behavior, and a red dot follows the pendulum’s current position on the graph.  The notebook will also show you how long it took to create the movie, which allows you to predict how long you’ll need to wait for a movie that has a lot of frames.

This setup is fairly flexible; you can change almost anything about the system and get an animation.  Simply run one of the cells below then change makeMove to True, and  the animation will appear in the directory specified in the SetDirectory command at the top.

This first cell shows no transients.

ρ=1.06;  (* amplitde of external force *)
ω=2*Pi;  (* The external force has period 1 *)
ω0=1.5*ω;  (* The natural frequency is 3/2 times the driving frequency *)
c = 1/2;  (* Damping constant *)
nTransient = 0; (* integer number of periods integrated before plotting *)
tfin=10; (* final time *)
ϕ0=- 2 Pi;  (* initial angle *)
v0=0;  (* inital velocity *)
makeMovie = False;
step = 1/30;  (* time step in frames of the movie, if made *)
If [makeMovie && tfin /step> 100, Print["You asked to make a movie with > 100 frames.  makeMovie set to false."];makeMovie = False];  
(* You can comment this out if your *really* want to make a long movie, but the export is long *)
filename = "period1.gif";  (* Can also use .mov or .avi or other video format *)
runPendulum

This second cell does 60 periods before plotting, and then does just one period  for the final time (tfin).  And it makes a movie.  The {phi, phi’} at the beginning and end of the plotting is the same to the accuracy shown.

ρ=1.06;  (* amplitde of external force *)
ω=2*Pi;  (* The external force has period 1 *)
ω0=1.5*ω;  (* The natural frequency is 3/2 times the driving frequency *)
c = 1/2;  (* Damping constant *)
nTransient = 60; (* integer number of periods integrated before plotting *)
tfin=1; (* final time *)
ϕ0=- 2 Pi;  (* initial angle *)
v0=0;  (* inital velocity *)
makeMovie = True;
step = 1/30;  (* time step in frames of the movie, if made *)
If [makeMovie && tfin /step> 100, Print["You asked to make a movie with > 100 frames.  makeMovie set to false."];makeMovie = False];  
(* You can comment this out if your *really* want to make a long movie, but the export is long *)
filename = "period1Asymmetric.gif";  (* Can also use .mov or .avi or other video format *)
runPendulum

The following is a symmetric oscillation.  This can be confirmed by using tfin = 0.5

ρ=.90;  (* amplitde of external force *)
ω=2*Pi;  (* The external force has period 1 *)
ω0=1.5*ω;  (* The natural frequency is 3/2 times the driving frequency *)
c = 1/2;  (* Damping constant *)
nTransient = 10; (* integer number of periods integrated before plotting *)
tfin=1; (* final time *)
ϕ0=0;  (* initial angle *)
v0=0;  (* inital velocity *)
makeMovie = False;
step = 1/30;  (* time step in frames of the movie, if made *)
If [makeMovie && tfin /step> 100, Print["You asked to make a movie with > 100 frames.  makeMovie set to false."];makeMovie = False];  
(* You can comment this out if your *really* want to make a long movie, but the export is long *)
filename = "period1Symmetric.gif";  (* Can also use .mov or .avi or other video format *)
runPendulum

The next movie shows a period 2.  with rho = 1.079 it didn’t settle down to period 2 with nTransient = 200.

ρ=1.078;  (* amplitde of external force *)
ω=2*Pi;  (* The external force has period 1 *)
ω0=1.5*ω;  (* The natural frequency is 3/2 times the driving frequency *)
c = 1/2;  (* Damping constant *)
nTransient = 200; (* integer number of periods integrated before plotting *)
tfin=2; (* final time *)
ϕ0=-.2;  (* initial angle *)
v0=18.5;  (* inital velocity *)
makeMovie = True;
step = 1/30;  (* time step in frames of the movie, if made *)
If [makeMovie && tfin /step> 100, Print["You asked to make a movie with > 100 frames.  makeMovie set to false."];makeMovie = False];  
(* You can comment this out if your *really* want to make a long movie, but the export is long *)
filename = "period2.gif";  (* Can also use .mov or .avi or other video format *)
runPendulum

(* Undriven undamped *)
filename = "undampedOverTop.gif";
ρ=0;
ω=2*Pi;
ω0=1.5*ω;  (* That is ω/ω0 = 2/3 *)
c = 0;
nTransient = 0;
tfin=1.406;  (* I found this by experimenting with tfin intil final phi = Pi *)
ϕ0=-Pi;
v0=0.1;
step = 1/30;
makeMovie = False;
runPendulum

(* Undriven undamped *)
filename = "undampedPeriodicSmallAmp.gif";
ρ=0;
ω=2*Pi;
ω0=1.5*ω;  (* That is ω/ω0 = 2/3 *)
c = 0;
nTransient = 0;
tfin=.6775;  (* I found this by experimenting with tfin = 5*0.67...  *)
ϕ0=0.5;
v0=0;
step = 1/30;
makeMovie = False;
runPendulum

(* Undriven undamped *)
filename = "undampedPeriodic.gif";
ρ=0;
ω=2*Pi;
ω0=1.5*ω;  (* That is ω/ω0 = 2/3 *)
c = 0;
nTransient = 0;
tfin=1.714;  (* I found this by experimenting with tfin = 5*1.7...  *)
ϕ0=3.0;
v0=0;
step = 1/30;
makeMovie = False;
runPendulum

makeMovie took 27.8 seconds with *.mov file and 30.0 sec with .gif

(* Undriven damped *)
filename = "undrivenDamped.gif";
ρ=0;
ω=2*Pi;
ω0=1.5*ω;
c = 1/10;
nTransient = 0;
tfin=7;
ϕ0=0;
v0=20.772;
step = 1/30;
makeMovie = False;
runPendulum

(* Period 3 attractor *)
ρ=1.077;
ω=2*Pi;
ω0=1.5*ω;
c = 1/2;
nTransient = 0;
tfin=15;
ϕ0=0;
v0=0;
step = 1/30;
filename = "dummy.gif";  (* I don't want to save this, but if I do by mistake write a gif file I don't want to over-write something good *)
makeMovie = False;
runPendulum

Now use that to get just the periodic solution.  I set tfin = 150 in the previous cell to get a good approximation to a point with period 3, and then used those final conditions as the intial condiditons in the next cell.

(* Period 3 attractor, from Figure 12.6 *)
filename = "period3.gif";
ρ=1.077;
ω=2*Pi;
ω0=1.5*ω;
c=1/2;
nTransient = 0;
tfin=3;
ϕ0=13.812249838984302`- 2 Pi;
v0=7.943343368297156`;
step = 1/30;
makeMovie = False;
runPendulum

(* Period 2 attractor *)
ρ=1.073;
ω=2*Pi;
ω0=1.5*ω;
c=1/2;
nTransient = 0;
tfin=50;
ϕ0=2 Pi;
v0=0;
step = 1/30;
filename = "dummy.gif";  (* I don't want to save this, but if I do by mistake write a gif file I don't want to over-write something good *)makeMovie = False;
runPendulum

(* Period 2 attractor *)
filename = "period2.gif";
ρ=1.073;
ω=2*Pi;
ω0=1.5*ω;
c=1/2;
nTransient = 0;
tfin=2;
ϕ0=-0.36061863303727104`;
v0=18.857553991046167`;
step = 1/30;
makeMovie = False;
runPendulum

(* Suspected Chaos *)
filename = "chaos.gif";
ρ=1.503;
ω=2*Pi;
ω0=1.5*ω;
c = 1/2;
nTransient = 0;
tfin=10;
ϕ0=-0.058410073786180294`;
v0=13.478925288721229`;
step = 1/30;
makeMovie = False;
runPendulum

(* Rolling motion, possibly chaotic *)
filename = "dummy.gif";
ρ=1.4;
ω=2*Pi;
ω0=1.5*ω;
c = 1/2;
nTransient = 0;
tfin=20;
ϕ0=-Pi/2 + 20 Pi ;
v0=0;
step = 1/30;
makeMovie = False;
runPendulum

(* Rolling motion*)
filename = "rolling.gif";
ρ=1.4;
ω=2*Pi;
ω0=1.5*ω;
c = 1/2;
nTransient = 0;
tfin=1;
ϕ0=-0.07228169156588723` + 2 Pi;
v0=8.805660683399568`;
step = 1/30;
makeMovie = False;
runPendulum