# 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. Jan 4 '13 at 13:13
• Dear Artes - Many thanks. 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 ... 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. 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. 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"}},
ScalingFunctions -> {{"Log10", "dB"}, {"Log10", "Degree"}},
PlotRange -> {{{0.001, 10}, All}, {{0.001, 10}, All}}][] & /@ {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. 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. Jan 4 '13 at 16:17
• @EdMendes Can you provide a reference for your assertion? 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". 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. 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. Jan 4 '13 at 22:23