PHYS 2063 Wave Physics:
The Forced Mechanical Oscillator

Remove["Global`*"]

The Equation of Motion

Write the equation of motion in terms of x, k, m, b, ω, and F0, as we did in class

eq0=m x''[t]+b x'[t]+k x[t]==F0 Cos[ω t];

Make the standard substitutions

eq1=Simplify[eq0/.{b->2m r,k->ω0^2m,F0->f0 m},m>0]

Recovering the simple harmonic oscillator

The simple, unforced, undamped harmonic oscillator, with arbitrary initial conditions

xSimple=x[t]/.DSolve[{eq0/.{F0->0,b->0},x[0]==x0,x'[0]==v0},x[t],t][[1]]

Simplify[xSimple/.{k->ω0^2m},ω0>0]

Good. This is just the same simple solution we started working on in this class.

Solving the full equation of motion

We can let Mathematica handle this without using a complex solution. The result looks horribly messy, mostly because of the transient components of the solution, but look at the end of the expression for the important denominator that affects the amplitude.

xForced=x[t]/.DSolve[{eq1,x[0]==x0,x'[0]==v0},x[t],t]

Plotting some unforced solutions

Store two versions of the solution, one underdamped and one overdamped

rSubU=r->ω0/10;
xForcedU=xForced/.rSubU;
xForcedO=xForced/.r->1.25ω0;

Pick some values for k and m. We will always assume SI units, i.e. mass, kilograms, and seconds (so Newtons for force and Joules for energy)

ω0Sub={ω0->Sqrt[k/m]};
kmVals={k->500,m->5};

Now make some plots with the forcing turned off.

noForce={ω->0,f0->0};
x0PlotU=xForcedU/.ω0Sub/.kmVals/.noForce;
x0PlotO=xForcedO/.ω0Sub/.kmVals/.noForce;

ics={x0->5,v0->0};
Plot[{x0PlotU/.ics,x0PlotO/.ics},{t,0,5}]

Plotting the forced oscillator

Try experimenting with the driving frequency and the initial conditions

f0Val=100;
f0ωSubs={f0->f0Val,ω->0.5ω0};
xPlotU=xForcedU/.f0ωSubs/.ω0Sub/.kmVals;
xPlotO=xForcedO/.f0ωSubs/.ω0Sub/.kmVals;

ics={x0->5,v0->0};
Plot[{xPlotU/.ics,xPlotO/.ics},{t,0,5}]

The transients die away after some time. Replot with longer time base.

Plot[{xPlotU/.ics,xPlotO/.ics},{t,0,25}]

Confirm that the driving period looks right

N[(2 Pi/ω)/.f0ωSubs/.ω0Sub/.kmVals]

Using Manipulate to study the motion

xManipulate=xForcedU/.{f0->f0Val}/.ω0Sub/.kmVals/.ics;

ω0Val=ω0/.ω0Sub/.kmVals;
ωMlo=ω0Val/5;
ωMhi=3ω0Val;

Calculate the expected amplitude at resonance

B=(f0/ω0^2)/.{f0->f0Val}/.ω0Sub/.kmVals;
Q=(ω0 m/b)/.b->2 m r/.rSubU;

The driving wave is plotted with an amplitude equal to that expected at resonance

tMax=8;
Manipulate[
Show[
Plot[xManipulate/.ω->ωM,{t,0,tMax},PlotRange->{-5,5}],
Plot[Q B Cos[ωM t],{t,0,tMax},PlotStyle->Red]],
{ωM,ωMlo,ωMhi}]

Is the resonance at the right frequency?

ω0Val

Animating the motion

Set values for the animation; Use the decay lifetime to set a maximum time for running the animation

side=2B;
dist=10B;
valsAnimate={ω->5,x0->0,v0->0};
xminAnimate=-2B;xmaxAnimate=dist+Q B+side;
dxAnimate=xmaxAnimate-xminAnimate;
tmax=(10/r)/.rSubU/.ω0Sub/.kmVals;

Get  expressions for u0(t) and x(t). Specify the first part of the solution expression because it will go into a list.

u0=B Cos[ω t]/.valsAnimate;
xAnimate=xForcedU[[1]]/.{f0->f0Val}/.ω0Sub/.kmVals/.valsAnimate;

Define the points that are the polygon which represent the block and the mass. The driver amplitude B is used for the scale.

ptsB={{u0-side/4,0},{u0+side/4,0},{u0+side/4,side},{u0-side/4,side}};
uM=dist+xAnimate;
ptsM={{uM-side/2,0},{uM+side/2,0},{uM+side/2,side},{uM-side/2,side}};

Animate: Plot sets the axes and aspect ratio, then the block and the mass.

Animate[
Show[
Plot[0,{x,xminAnimate,xmaxAnimate},
PlotRange->{{xminAnimate,xmaxAnimate},{0,side}},
AspectRatio->side/dxAnimate,Axes->False],
Graphics[{Red,Polygon[ptsM/.t->time]}],
Graphics[{Black,Polygon[ptsB/.t->time]}]
],{time,0,tmax},AnimationRate->0.5,AnimationRunning->False]