Phys 238:  Driven Oscillations
Resonance

Clear["Global`*"]; DateString[]

SetDirectory[NotebookDirectory[]]; (* Find and save files alongside this notebook. *)

Differential Equation and Basic Parameters

The governing differential equation is:

    θ’’[t] = - ((ω0)^2) θ(t) - γ θ’(t) + α0 Sin[ωd t]
    
    where κ is the restoring torque, I is the moment of inertia, ω0 = Sqrt[ κ / I],  γ = drag, and α0 = τ0 / I, where τ0 is the amplitude of the driving torque, ωd is the frequency of the driving.   Related quantities are Q = ω0/γ and
    linear frequency f0 = ω0/2π  and fd = ωd/2 π   
    
     θ’’[t] = - ((ω0)^2) θ(t) -  (ω0/Q) θ’(t) + α0 Sin[ωd t]

Physical constants

ω0 = 4.4 (* rad/s, typical for torsional oscillator with 2 brass quadrants. *)
f0 = ω0/(2 π)
tmax1 = 60 (* Useful simulation limit *)

Do the numerical integration

damped = ParametricNDSolveValue[
{θ''[t] == -(2 π f0)^2 θ[t] - (2 π f0/Q)θ'[t] , θ[0] == 1, θ'[0] == 0},
  θ,
 {t, 0, tmax1}, {Q}]

Damped Oscillations -- the effect of Q

Manipulate[Plot[damped[Q][t], {t, 0, tmax1}, PlotRange ->{-1, 1}, PlotPoints -> 300, PlotLabel -> StringForm["Q = ``, (Subscript[t, 1/2]) = ``s", Q, NumberForm[2 Log[2] Q/ω0, {4, 2}]], LabelStyle->Larger, AxesLabel -> {"t (s)", "Angle (rad)"}, ImageSize->Scaled[0.8]],
{{Q, 20}, 1, 40, 1, Appearance->"Open"}]

Driven Oscillations -- initial transients

Set initial position to 0.  Set a longer integration time to see transients.

tmax2 = 90;

driven = ParametricNDSolveValue[
{θ''[t] == -(2 π f0)^2 θ[t] - (2 π f0/Q)θ'[t] + α0 Sin[2 π fd t], θ[0] == 0, θ'[0] == 0},
  θ,
 {t, 0, tmax2}, {Q, α0, fd}]

Manipulate[Plot[driven[Q, α0, fd][t], {t, 0, tmax2},PlotPoints -> 300,  PlotRange -> All, PlotLabel -> StringForm["Q = ``", Q], LabelStyle->Larger, AxesLabel -> {"t (s)", "Angle (rad)"}, ImageSize->Scaled[0.8]],
{{Q, 30}, 1, 40, 1, Appearance->"Labeled"},
{{α0, 5}, 1, 20, 1, Appearance -> "Labeled"},
{{fd,0.6},  0.2, 2, 0.02, Appearance -> "Open"}]

Steady State Response -- Resonance

resonance[α0_, f0_, f_, Q_] := ((α0)/(Sqrt[((((f0^2) - (f^2)))^2) + (((((f0 f)/Q)))^2)]))

Manipulate[Plot[resonance[α0, f0, fd, Q], {fd, 0, 2 f0},PlotPoints -> 300,  PlotRange -> All, PlotLabel -> StringForm["Q = ``", Q], LabelStyle->Larger, AxesLabel -> {"fd (Hz)", "Amplitude"}, ImageSize->Scaled[0.8]],
{{Q, 30}, 1, 40, 1, Appearance->"Labeled"},
{{α0, 5}, 1, 20, 1, Appearance -> "Labeled"}]

Plot[{resonance[1, f0, fd, 5], 
resonance[1, f0, fd, 10],
resonance[1, f0, fd, 15]},{fd, 0.4, 1.0},PlotPoints -> 300,  PlotRange -> All, PlotLegends -> {"Q = 5", "Q = 10", "Q = 15"},  LabelStyle->Larger, AxesLabel -> {"fd (Hz)", "Amplitude"}, ImageSize->Scaled[0.8]]

phase[f0_, f_, Q_] :=ArcTan[ (f0^2) - (f^2) , f0 f/Q]

Manipulate[Plot[{phase[f0, fd, Q], π}, {fd, 0, 2 f0},PlotPoints -> 300,  PlotRange -> {0, π}, PlotLabel -> StringForm["Q = ``", Q], PlotLegends -> {"Phase", "π"}, LabelStyle->Larger, AxesLabel -> {"fd (Hz)", "Phase"}, ImageSize->Scaled[0.8]],
{{Q, 5}, 1, 40, 1, Appearance->"Open"}]

Plot[{phase[f0, fd, 5], 
phase[f0, fd, 10],
phase[f0, fd, 15]},{fd, 0.4, 1.0},PlotPoints -> 300,  PlotRange -> All, PlotLegends -> {"Q = 5", "Q = 10", "Q = 15"},  LabelStyle->Larger, AxesLabel -> {"fd (Hz)", "Phase (radians)"}, ImageSize->Scaled[0.8]]