How to fix BodePlot that comes with Mathematica?

I am using Mathematica to go through the examples and exercises on the book Modern Control Systems, 6th edition by Dorf. On page 605, there is a table (Table 8.5) with the Bode plot for several transfer functions. In what follows there is a piece of code that attempts to build the very same table.

Here is the code:

With[{ τ1 = 20, τ2 = 2, τ3 = 0.4, τ4 = 0.05, τa = 10, τb = 1, k = 10},
Grid[
Partition[
Table[ BodePlot[ sys, PlotLabel->sys, GridLines -> Automatic], { sys,
{ k/(s τ1 + 1), (k(s τa + 1))/(s(s τ1 + 1)(s τ2 + 1)),
k/((s τ1 + 1)(s τ2 + 1)), k/s^2, k/((s τ1 + 1)(s τ2 + 1)(s τ3 + 1)),
k/(s^2 (s τ1 + 1)), k/s, (k(s τa + 1))/(s^2 (s τ1 + 1)),
k/(s(s τ1 + 1)), k/s^3, k/(s(s τ1 + 1)(s τ2 + 1)),
(k (s τa + 1))/s^3, (k (s τa + 1)(s τb + 1))/s^3,
(k (s τa + 1))/(s^2 (s τ1 + 1)(s τ2 + 1)),
(k (s τa + 1)(s τb + 1))/(s(s τ1 + 1)(s τ2 + 1)(s τ3 + 1)(s τ4 + 1)) }
}
],  2], Frame->All, Spacings->6] ]


All the transfer functions with 1/s^n ( n > 1 ) give the wrong result as far as the phase plot is concerned. Is there a simple way to fix this? Wolfram does not have a time line to go through the problem and sort it out.

• Could you explain why "the result" is wrong ??? I edited your question, you can see how it should be done clicking edit under your question. – Artes Jan 4 '13 at 13:13
• Dear Artes - Many thanks. – Ed Mendes Jan 4 '13 at 16:20
• Dear Artes - Many thanks. Wrong means that it does not follow the standard convention when plotting a Bode Diagram. Check all the most adopted textbooks on Control Systems and see the Bode Diagram for the transfer function 1/s^2. I bet that in all of them the phase is -180 and not +180 as Mathematica BodePlot shows you. Why don't they follow the convention? 1 pole at s=0 gives -90, two poles at s=0 gives -180, 3 poles at s=0 gives -270 ... – Ed Mendes Jan 4 '13 at 16:27
• @EdMendes If you want to write a comment to a selected user write @name (e.g. to me @Artes). Have you read Details and Options in the documentation pages of BodePlot ? It says : option specifications include: opt->val use val for both the magnitude and the phase plot. If you find a correct solution to the problem you can answer your own question as well. This is a welcome practice. – Artes Jan 5 '13 at 21:27

This is fixed for v10. In v10 there is also a PhaseRange option for BodePlot and NicholsPlot that can be used to override the default range (if you need to wrap it between $-\pi$ and $\pi$, etc).

With[{τ1 = 20, τ2 = 2, τ3 = 0.4, τ4 = 0.05, τa = 10, τb = 1, k = 10},
Grid[Partition[Table[BodePlot[sys, PlotLabel -> sys, GridLines -> Automatic,
PlotLayout -> "Phase"], {sys, {k/(s τ1 + 1),
(k (s τa + 1))/(s (s τ1 + 1) (s τ2 + 1)),
k/((s τ1 + 1) (s τ2 + 1)), k/s^2,
k/((s τ1 + 1) (s τ2 + 1) (s τ3 + 1)),
k/(s^2 (s τ1 + 1)),
k/s, (k (s τa + 1))/(s^2 (s τ1 + 1)),
k/(s (s τ1 + 1)), k/s^3,
k/(s (s τ1 + 1) (s τ2 + 1)), (k (s τa + 1))/
s^3, (k (s τa + 1) (s τb + 1))/
s^3, (k (s τa + 1))/(s^2 (s τ1 + 1) (s τ2 +
1)), (k (s τa + 1) (s τb + 1))/(s (s τ1 +
1) (s τ2 + 1) (s τ3 + 1) (s τ4 + 1))}}], 2],
Frame -> All, Spacings -> 6]]


• Many thanks for the reply. – Ed Mendes Oct 1 '14 at 23:26

To address the above comment by Ed

However if we change the example. Matlab - bode(tf(10*[10 1],[1 0 0 0])) - phase is negative (-270 to -180). Mathematica - See plots above - phase from +90 to +180. If instead of 10(10s+1)/sˆ3, one uses 10(10s+1)/(s+0.0001)ˆ3, the phase is negative

This below is a direct implementation of the phase plot part of Bode, using ArcTan. I used the following 2 transfer functions to compare with Mathematica BodePlot: 10(10s+1)/(s+0.0001)ˆ3 and 10(10s+1)/(s)ˆ3 and the result does show that there is a sudden phase change shift by 180 which does not show when using straight calculations using ArcTan to find the phase. Conclusion: There seems to be some convention used that causes this change as I would have expected it to match the ArcTan direct method.

Mathematica BodePlot phase diagram

Clear[s];
expr1 = (100 s + 10)/(s)^3;
expr2 = (100 s + 10)/(s + 0.0001)^3;

Grid[{{expr1, expr2},
BodePlot[TransferFunctionModel[#, s], GridLines -> Automatic,
ImageSize -> 300, PlotLayout -> "List",
FrameLabel -> {{{"magnitude (db)", None}, {None,"Bode plot"}},
{{"phase(deg)",None}, {"Frequency (rad/sec)", None}}},
ScalingFunctions -> {{"Log10", "dB"}, {"Log10", "Degree"}},
PlotRange -> {{{0.001, 10}, All}, {{0.001, 10}, All}}][[2]] & /@ {expr1, expr2}}]


Direct implementation of the phase plot using ArcTan

Clear[s];
ticks[min_, max_] := Table[{i, Superscript[10, i]}, {i, Ceiling[min], Floor[max], 1}];

makePhasePlot[expr_, s_Symbol] := Module[{ex, w, re, im, data},

ex = expr /. s -> (w I);
re = ComplexExpand[Re[ex]];
im = ComplexExpand[Im[ex]];

data =Table[{Log[10, w], 180/Pi ArcTan[re, im]}, {w, 0.001, 10, 0.001}];

ListPlot[data,
Joined -> True,
PlotRange -> All,
FrameTicks -> {{Automatic, Automatic}, {ticks, Automatic}},
AxesOrigin -> {0, 0},
Frame -> True,
ImageSize -> 300,
FrameLabel -> {{"angle(deg)", None}, {"Frequency (rad/sec)", None}},
Axes -> False]
];

expr1 = (100 s + 10)/s^3;
expr2 = (100 s + 10)/(s + 0.0001)^3;
Grid[{{expr1, expr2},makePhasePlot[#, s] & /@ {expr1, expr2}}, Spacings -> {3, 0}]


• Thanks Nasser. That is exactly the point what I was trying to get across. There is a conflict within BodePlot that other control toolboxes do not seem to have (at least not in this case). Thanks to Mr. Wizard and Mr. James Cunnane for showing that I should not have used the word decent. – Ed Mendes Jan 5 '13 at 12:30

No error seen in those phase plots.

For a transfer function of the form K/s^n, i.e. n poles at the origin, we expect constant phase of (-90 * n) degrees, plus or minus some integer multiple of 360 degrees - which is exactly what your Mathematica plots show.

• I disagree - although -180 e +180 are the same thing - Mathematica BodePlot does not follow the usual standard convention when plotting a Bode Diagram - If you have a pole at s=0 the phase is -90 (270 sounds really weird). If you have two poles at s=0, one would expect -90 + -90 = -180 and not +180 (this means that the phase margin is 0 and not 360). Any decent control toolbox returns the angle in the standard convention. – Ed Mendes Jan 4 '13 at 16:17
• @EdMendes Can you provide a reference for your assertion? – Mr.Wizard Jan 4 '13 at 16:37
• @EdMendes Contrariwise - see MATLAB Bode plot documentation; the first and last Bode plots all show phases greater than 180 degrees. Unless your comment is intended to suggest that MATLAB is not included in "any decent control toolbox". – James Cunnane Jan 4 '13 at 18:05
• Many thanks for pointing out the Matlab documentation page. Please note that I am not talking about "phases greater than 180". Let us use the Matlab example - H = tf([1 0.1 7.5],[1 0.12 9 0 0]); bode(H) - In my mac Pro it shows negative phases (Actually the phase starts from -180 (due to the poles at s=0) goes near to -45 and comes back to -180. – Ed Mendes Jan 4 '13 at 22:04
• Let us do the same thing with Mathematica. Thanks to you I can get the same result if I use BodePlot[TransferFunctionModel[(s^2 + 0.1*s + 7.5)/(s^4 + 0.12*s^3 + 9*s^2), s]]. That is great! Many thanks. However if we change the example. Matlab - bode(tf(10*[10 1],[1 0 0 0])) - phase is negative (-270 to -180). Mathematica - See plots above - phase from +90 to +180. If instead of 10(10s+1)/sˆ3, one uses 10(10s+1)/(s+0.0001)ˆ3, the phase is negative. – Ed Mendes Jan 4 '13 at 22:23