I am trying to plot a phase angle vs. frequency curve numerically for the steady state condition of the following equation
x''[t] + 2 ζ x'[t] + x[t] == Cos[w t]
where w is the forcing frequency, and the initial conditions are x[0] == 0, x'[0] == 0. A similar plot is shown here:
Method 1 : Using Bode Plot from signals and systems
One more direct solution can be by making use of Mathematica's inbuilt BodePlot. You just need to find the Transfer function, which can be easily found in your equation by taking the Laplace Transform on both sides of the equation and is given as:
H[s_] := 1/(s^2 + 2 \[Zeta] s + 1)
Now, use the BodePlot function, just for Phase by this code:
BodePlot[1/(s^2 + 2 (#) s + 1) & /@ {0.005, 0.1, 0.38, 1},
PlotLayout -> "Phase",
PlotLegends ->
Placed[LineLegend[{0.005, 0.1, 0.38, 1}, LegendLabel -> "\[Zeta]",
LegendFunction -> (Framed[#, RoundingRadius -> 5] &),
LegendMargins -> 3], {0.8, 0.7}]]
Which gives you a similar answer as given by @Nasser. You just need to change your units to the ones described in the reference.
Method 2: Using Direct formula from the theory of vibrations
You can easily find the theory of forced harmonic vibrations in any of the books on mechanical vibrations, like Theory of Vibration with Applications by William Thomson and Theory of Vibration, An Introduction by Ahmed A. Shabana, which directly provides you with the formula for phase angle ($\psi$) in terms of damping factor ($\zeta$) and frequency ratio ($r$) as:
$$ \psi = \tan^{-1} \bigg( \frac{2 r \zeta}{1-r^2} \bigg) $$
Then it's really easy to plot the phase-frequency curve:
\[Psi][\[Zeta]_, r_] := ArcTan[1 - r^2, 2 r \[Zeta]];
Plot[Evaluate@
Table[\[Psi][\[Zeta], r], {\[Zeta], {0.005, 0.1, 0.38, 1}}], {r, 0,
4}, PlotRange -> {{0, 3}, {0, 3.5}}, AspectRatio -> 1/GoldenRatio,
PlotPoints -> 1000, PlotTheme -> {"Scientific", "NeonColor"},
FrameLabel -> {"Frequency ratio \[Rule]", "Phase angle \[Rule]"},
LabelStyle -> {Black, Bold, 15},
PlotLegends ->
Placed[SwatchLegend[Automatic, {0.005, 0.1, 0.38, 1},
LegendMarkers -> "Bubble", LegendMarkerSize -> 20,
LegendFunction -> (Framed[#, RoundingRadius -> 10,
Background -> GrayLevel[0.85]] &),
LegendLabel -> "\[Zeta]"], {0.9, 0.25}], ImageSize -> Large]
which gives you the desired plot:
x''[T] == (1/\[Epsilon])*((((1 + WE)*(1 + \[Epsilon]*x[T])^(-3*\[Kappa] -1))) - ((1 + \[Epsilon]* x[T])^-1) - (WE*(1 + \[Epsilon]*x[T])^-2) - (4*damping*\[Epsilon]^2*(x'[T])*(1 + (\[Epsilon]* x[T]))^-2) - ((3*\[Epsilon]^2*(x'[T]^2)*(1 + \[Epsilon]*x[T])^-1)/2) - \[Epsilon]^3*f*Sin[T* (w0 + \[Epsilon]^2 w)]) with the initial conditions x[0]=1, x'[0]=0
$\endgroup$
WE=0.2, \[Kappa]=1.4, \[Epsilon]=0.1, and damping=0.03
$\endgroup$
Code
phase[r_, z_] := Module[{t}, t = ArcTan[(2 z r)/(1 - r^2)];
If[t < 0, t = t + Pi];
180/Pi t];
plotOneZeta[z_, f_] := Module[{r, p1, p2},
p1 = Plot[f[r, z], {r, 0, 3}, PlotRange -> All,
PlotStyle -> Blue];
p2 = Graphics[Text[z, {1.1, f[1.1, z]}]];
Show[{p1, p2}]];
p = Map[plotOneZeta[#, phase] &, {.005, 0.1, 0.38, 1}];
Show[p, FrameLabel -> {{"Phase in degrees",
None}, {"r= ω/ωn",
"Phase vs. Frequency ratio for different ξ"}}, Frame -> True,
GridLines -> Automatic, GridLinesStyle -> Dashed, ImageSize -> 400,
AspectRatio -> 1]
Use degree for phase not radians. You can easily change that.
ps. I stole the code from this page
