The Vertically Driven Pendulum

Solving the Lagrangian

Setup the problem with Potential Energy (U) and Kinetic Energy (T)
The pendulum suport point moves as a function of time, yd[t].  The coordinates of the pendulum mass are xm and ym.

xm[t]=l Sin[phi[t]];
ym[t]=yd[t]-l Cos[phi[t]];
U=Simplify[m g ym[t]]
T = Simplify[(1/2) m (D[xm[t],t]^2+D[ym[t],t]^2)]

Solve the Euler-Lagrange Equation, with Damping Term

Fphi=Simplify[D[T-U,phi[t]]]-b D[phi[t],t]
pphi=Simplify[D[T-U,phi'[t]]]
dpphi=Simplify[D[pphi,t]]
eomd=Simplify[Fphi==dpphi,Assumptions->{m>0,l>0}]

Solve the Equations of Motion

List dependent variables

dvars={phi};

Set initial condisions and drive function

incon={phi[0]==phi0,phi'[0]==v0};
y1[t]=ad Sin[fd Sqrt[g/l] t];
eom=eomd/.{yd''[t]->D[D[y1[t],t],t]}

Define all parameters

pars={m->1,g->10,l->1,b->0.1,ad->0,fd->0.1,phi0->0.1,v0->0};
tmax=30 2π/Sqrt[g/l]/.pars;

Solve with NDSolve

sol=NDSolve[{eom,incon}//.pars,dvars,{t,0,tmax},MaxStepSize->N[1/(10 fd)/.pars]]

Plot the solution

plotpars={PerformanceGoal->"Quality",ImageSize->Full,AspectRatio->1/4};
Plot[Evaluate[phi[t]/.sol],{t,0,tmax},Evaluate[plotpars]]

Solve the Equations of Motion:
No dirve, big swing with damping

Try different values of b and note how the damping changes

pars={m->1,g->10,l->1,b->0.3,ad->0,fd->0.1,phi0->π/2,v0->0};
tmax=5 2π/Sqrt[g/l]/.pars;
sol=NDSolve[{eom,incon}//.pars,dvars,{t,0,tmax},MaxStepSize->N[1/(10 fd)/.pars]]
Plot[Evaluate[phi[t]/.sol],{t,0,tmax},PlotRange->Full,Evaluate[plotpars]]

Solve the Equations of Motion:
Parametric Resonance

Try different values of fd and note how the behavior changes at fd=2

pars={m->1,g->10,l->1,b->0.1,ad->0.1,fd->2,phi0->0.1,v0->0};
y1[t]/.pars
tmax=20 Sqrt[g/l]/.pars;
sol=NDSolve[{eom,incon}//.pars,dvars,{t,0,tmax},MaxStepSize->N[1/(10 fd)/.pars]]
Plot[{Evaluate[phi[t]/.sol],Evaluate[5y1[t]/.pars]},{t,0,tmax},PlotRange->Full,Evaluate[plotpars]]

Solve the Equations of Motion:
High Frequency Drive - strange

pars={m->1,g->10,l->1,b->0.1,ad->0.1,fd->14,phi0->π+0.1,v0->0};
tmax=5 Sqrt[g/l]/.pars;
sol=NDSolve[{eom,incon}//.pars,dvars,{t,0,tmax},MaxStepSize->N[1/(10 fd)/.pars]]
Plot[Evaluate[phi[t]/.sol],{t,0,tmax},Evaluate[plotpars]]

Solve the Equations of Motion:
Higher Frequency Drive - Kapitza

Try different values of fd and note how the stability and frequency of the π oscillations change

pars={m->1,g->10,l->1,b->0.1,ad->0.1,fd->16,phi0->π-0.1,v0->0};
tmax=5 Sqrt[g/l]/.pars;
sol=NDSolve[{eom,incon}//.pars,dvars,{t,0,tmax},MaxStepSize->N[1/(10 fd)/.pars]]
Plot[Evaluate[phi[t]/.sol],{t,0,tmax},Evaluate[plotpars]]