How Can I Plot the Phase-Frequency Curve Using Mathematica?

Question

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:

Answer 1

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:

Answer 2

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