Pendulum Examples: Nonlinear Phenomena for 29 Nov 2022

Remove["Global`*"]

First a little about solving Ordinary Differential Equations

Define a differential equation with initial conditions

eq=x''[t]+2β x'[t]+(Subscript[ω, 0])^2x[t]==0;
ic1=x[0]==A;
ic2=x'[0]==B;

Solve the differential equation analytically

sol=DSolve[{eq,ic1,ic2},x[t],t]

xA=x[t]/.sol[[1]]

Put in some values and plot the solution

vals={(Subscript[ω, 0])->2Pi,β->0.3,A->5,B->0};
xAplot=Re[Simplify[ComplexExpand[xA/.vals]]]

tMax=5;
Plot[xAplot,{t,0,tMax}]

You can also solve the equation analytically with the values already applied

eqVals=eq/.vals

ic1Vals=ic1/.vals;
ic2Vals=ic2/.vals;

sol=DSolve[{eqVals,ic1Vals,ic2Vals},x[t],t];
xAplot=x[t]/.sol[[1]]

Plot[xAplot,{t,0,tMax}]

Now solve the differential equation numerically

sol=NDSolve[{eqVals,ic1Vals,ic2Vals},x[t],{t,0,tMax}];
xNplot=x[t]/.sol[[1]]

Plot[xNplot,{t,0,tMax}]

Plot the difference between the two solutions

Plot[xAplot-xNplot,{t,0,tMax}]

Numerical solution of the free planar pendulum

Set up the fundamental differential equation

eq=θ''[t]+(Subscript[ω, 0])^2Sin[θ[t]]==0;
ic1=θ[0]==(Subscript[θ, 0]);
ic2=θ'[0]==0;

Set values for low and high initial values

valsLo={(Subscript[ω, 0])->2Pi,(Subscript[θ, 0])->0.1};
valsHi={(Subscript[ω, 0])->2Pi,(Subscript[θ, 0])->0.99Pi};

Now solve the differential equation numerically for each of the two different initial conditions

eqVals=eq/.valsLo;
ic1Vals=ic1/.valsLo;
ic2Vals=ic2/.valsLo;

solLo=NDSolve[{eqVals,ic1Vals,ic2Vals},θ[t],{t,0,tMax}];
θplotLo=θ[t]/.solLo[[1]];

eqVals=eq/.valsHi;
ic1Vals=ic1/.valsHi;
ic2Vals=ic2/.valsHi;

solHi=NDSolve[{eqVals,ic1Vals,ic2Vals},θ[t],{t,0,tMax}];
θplotHi=θ[t]/.solHi[[1]];

Plot the two solutions

Plot[θplotLo,{t,0,tMax}]

Plot[θplotHi,{t,0,tMax}]

Make an animation (switch below between solHi and solLo)

solAnim=solHi;

len=1;
rBob=0.05 len;
xBob=len Sin[θ[t]]/.solAnim;
yBob=-len Cos[θ[t]]/.solAnim;
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,rBob]}]],
{time,0,tMax},
AnimationRate->0.5,
AnimationRunning->False]

Now add linear drag to the pendulum

Add the linear drag term to the differential equation

eq=θ''[t]+2β θ'[t]+(Subscript[ω, 0])^2Sin[θ[t]]==0;

Solve and plot

eqVals=eq/.valsHi/.β->0.25;

solDrag=NDSolve[{eqVals,ic1Vals,ic2Vals},θ[t],{t,0,tMax}];
θplotDrag=θ[t]/.solDrag[[1]];

Plot[θplotDrag,{t,0,tMax}]

Animate

len=1;
rBob=0.05 len;
xBob=len Sin[θ[t]]/.solDrag;
yBob=-len Cos[θ[t]]/.solDrag;
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,rBob]}]],
{time,0,tMax},
AnimationRate->0.5,
AnimationRunning->False]

Chaos in the driven, damped pendulum

Add the driving term

eq=θ''[t]+2β θ'[t]+(Subscript[ω, 0])^2Sin[θ[t]]==γ (Subscript[ω, 0])^2  Cos[ω t];

Set both initial conditions to zero

ic1=θ[0]==0;
ic2=θ'[0]==0;

Use a standard set of parameters and prepare to change the driving term strength

eqBase=eq/.β->(Subscript[ω, 0])/4/.(Subscript[ω, 0])->1.5ω/.ω->2Pi

Solve, plot, and animate using a weak driving term

eqWeak=eqBase/.γ->0.1

tMax=10;
solWeak=NDSolve[{eqWeak,ic1,ic2},θ[t],{t,0,tMax}];

θplotWeak=θ[t]/.solWeak[[1]];
Plot[θplotWeak,{t,0,tMax}]

len=1;
rBob=0.05 len;
xBob=len Sin[θ[t]]/.solWeak;
yBob=-len Cos[θ[t]]/.solWeak;
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,rBob]}]],
{time,0,tMax},
AnimationRate->0.5,
AnimationRunning->False]

Now use a strong driving term

γVal=1.0;

eqStrong=eqBase/.γ->γVal;

tMax=40;
solStrong=NDSolve[{eqStrong,ic1,ic2},θ[t],{t,0,tMax}];

θplotStrong=θ[t]/.solStrong[[1]];
Plot[θplotStrong,{t,0,tMax}]

Focus in on long times, and check visually for periodicity

Plot[θplotStrong,{t,tMax/2,tMax},PlotRange->{-3,-2}]

For fun, show the animation

len=1;
rBob=0.05 len;
xBob=len Sin[θ[t]]/.solStrong;
yBob=-len Cos[θ[t]]/.solStrong;
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,rBob]}]],
{time,0,tMax},
AnimationRate->0.5,
AnimationRunning->False]

Now go back and increase the driving term,, and watch the transition through bifurcations and then chaos.
The values below are from Taylor Figure 12.8 and Figure 12.10

γ1=1.06;
γ2=1.078;
γ3=1.082;
γ4=1.0826;
γ5=1.105;