<< Graphics`Graphics`
<< LinearAlgebra`MatrixManipulation`
dB[z_] = 20 Log[10, Abs[z]] ; 
    dg[z_] = Im[Log[z]]/Degree ;

gr = (s + 1)/(α s + 1) 
          10/(s + 10);
gp = 1/m/s^2;
go = gr gp

para = {α ->  0.4, m ->  1};

example 6.16   magnetic levitation system      eigenvalues   p.124

chr = Numerator[Factor[1 + go]] /. para

pole = Solve[chr == 0, s]

λ = s /. pole[[3]];
Print["ωn=", 
      Abs[λ], "  ζ=", -Re[λ]/
        Abs[λ]]

example 6.17   Routh-Hurwitz criterion  p.126

chrg = Numerator[Factor[1 + kg go]] /. para

h3 = {{5, 10 kg, 0}, {0.4, 10, 10 kg}, {0, 5, 
            10 kg}};
MatrixForm[h3]
h2 = TakeMatrix[h3, {1, 1}, {2, 2}]; MatrixForm[h2]

Det[h3]

    Det[h2]

example 6.18   root locus p.127

loc[k_] := {pole = Solve[chrg == 0 /. kg ->  k, s]; 
        polelist = Table[{Re[s], Im[s]} /. pole[[i]], {i, 4}]; 
        ListPlot[polelist, PlotRange ->  {{-2, 0.5}, {-5, 5}}, 
          PlotStyle ->  {Hue[0.9], PointSize[0.03]}, 
          PlotLabel ->  "k="[k]]}
loc[6.25]

zuall = Table[loc[k], {k, 0, 7}]

Show[zuall]

example 6.19    Nyquist plot of Go(s)   p. 130

goω = go /. para /. 
          s ->  I ω;
zu2 = 
      ParametricPlot[{Re[goω], Im[goω]}, {ω, 0.1, 
          10}]

goω625 = 6.25 goω;
zu3 = 
      ParametricPlot[{Re[goω625], Im[goω625]}, {ω, 0.1, 
          10}]

Show[zu2, zu3, PlotRange ->  {{-2, 0.15}, {-1, 0.15}}, 
      AspectRatio ->  Automatic]

example 6.20   Bode plot of Go(s)  p. 132

zu4 = 
      LogLinearPlot[{dB[goω], dg[goω] 60/180}, {ω, 0.1, 
          10}, PlotRange ->  {{0.1, 10}, {-60, 60}}, 
        GridLines ->  Automatic]

example 6.21    closed loop   p.133

gc = go/(1 + go);

gcω = gc /. para /. 
          s ->  I ω;
LogLinearPlot[{dB[gcω], dg[gcω] 60/180 - 20}, {ω, 
        0.1, 100}, GridLines ->  Automatic, 
      PlotRange -> {{0.1, 100}, {-80, 40}}, 
      PlotPoints ->  1000]

unbalance

gu = m ϵ ω^2  gp/(1 + go);

guω = gu /. para /. ϵ ->  1 /. 
          s ->  I ω;
Plot[Abs[guω], {ω, 0, 5}, PlotRange ->  {0, 5}, 
      GridLines ->  Automatic]

step response and setady-stae error

xt = Re[InverseLaplaceTransform[gc/s /. para, s, t]];

Plot[xt, {t, 0, 20}, PlotRange ->  {0, 2}, 
      GridLines ->  Automatic]

intergal function for step disturbance

xds = 1/s gp/(1 + (gr + ki/s) gp) /. para

xdt0 = 
        Re[InverseLaplaceTransform[xds /. ki ->  0, s, 
            t]];

Plot[xdt0, {t, 0, 40}, PlotRange ->  {-1, 2}, 
      GridLines ->  Automatic]

xdt1 = 
        Re[InverseLaplaceTransform[xds /. ki ->  0.1, s, 
            t]];
xdt3 = 
        Re[InverseLaplaceTransform[xds /. ki ->  0.3, s, 
            t]];

Plot[{xdt0, xdt1, xdt3}, {t, 0, 40}, PlotRange ->  {-1, 2}, 
      GridLines ->  Automatic]