MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Bode Plots in Mathematica

  • To: mathgroup at smc.vnet.net
  • Subject: [mg57232] Bode Plots in Mathematica
  • From: "David Park" <djmp at earthlink.net>
  • Date: Sat, 21 May 2005 02:40:16 -0400 (EDT)
  • Sender: owner-wri-mathgroup at wolfram.com

Following the discussion on the recent Bode Plots thread, along with DrBob, and learning from the advice of Mariusz Jankowski and others, I wrote a DrawGraphics BodePlot routine.

As I mentioned in a previous email, LogLinearPlot (and also LogLogPlot) has a serious flaw in that it sample points on a linear scale rather than a log scale. If the plot has a x scale that spans more than one or two decades, this results in severe underrepresentation of low values. In several of the examples below it is necessary to go to almost 10000 plot points to obtain a smooth curve at low values, when using LogLinearPlot. This is not only grossly inefficient, but results in over plotting at the high x values and resultant ragged lines. The routine below uses parameter transformation and the regular Plot routine to obtain an even sampling of points, and thus behaves as expected.

The db and phase functions are plotted full scale in the plot. The actual plot coordinates are the power domain for frequency s, and 0 to 100 for the y values. The db and phase curves are scaled to fit into the y range. NMinimize and NMaximize are used to obtain the minimum and maximum values of the curves. I found that the Automatic method in these routines did not always work, but the "DifferentialEvolution" method seem to work for all the cases considered. CustomTicks is used to make db and phase scales on the left and right and also the x scale. CustomGridLines is used to put a light log scale grid on the s axis and a linear grid for the y axis. The y grid is keyed to the underlying y coordinates and will not generally match either the db or phase ticks. The routine also uses the LinearScale routine from Graphics`Graphics`.

I have added the routine to the Examples section of the DrawGraphics Help, under Log Plots & Bode Plots so those who download the 20 May 2005 version will have the code available there. (DrawGraphics can be obtained from my web site below.)

Here is the code and the examples DrBob was using.

Needs["DrawGraphics`DrawingMaster`"]
Needs["Graphics`Graphics`"]

BodePlot::usage = 
    "BodePlot[h, {mindecade,maxdecade}, phasemode:0, opts] will make a Bode \
plot of the pure function h. The log frequence range will be from mindecade \
to maxdecade. The optional argument, phasemode, specifies the range for the \
phase function: -1 gives {-360,0}, 0 gives {-180, 180} and 1 gives {0,360}. \
Options may be passed to the plotting routine.";

BodePlot[hfunc_, {mindecade_, maxdecade_}, phasemode_:0, opts___?OptionQ] := 
    Module[{a, b, c, d, phasemin, phasemax, dbmin, dbmax, h, logstep, 
        logvalue, db, phase, plotdb, plotphase, nformat, work, diff, 
        dbticksiterator, phaseticksiterator},
      
      nformat :=
        Switch[#,
            1, 1,
            10, 10,
            1/10, 0.1,
            _, DisplayForm[SuperscriptBox[10, Log[10, #]]]] &;
      
      logstep[f_][z_] := f[10^z];
      logvalue[expr_] := Log[10, expr];
      db[h_][s_] := 20 Log[10, Abs[hfunc[I s]]];
      phase[h_][s_] :=
        Module[{arg = (Arg[hfunc[I s]] 180)/Pi},
          Switch[phasemode,
            -1, If[arg >= 0, arg - 360, arg],
            1, If[arg <= 0, arg + 360, arg],
            _, arg]];
      
      dbmin = 
        First[NMinimize[{logstep[db[h]][s], mindecade <= s <= maxdecade}, {s},
             Method -> "DifferentialEvolution"]];
      dbmax = 
        First[NMaximize[{logstep[db[h]][s], mindecade <= s <= maxdecade}, {s},
             Method -> "DifferentialEvolution"]]; 
      phasemin = 
        First[NMinimize[{logstep[phase[h]][s], 
              mindecade <= s <= maxdecade}, {s}, 
            Method -> "DifferentialEvolution"]]; 
      phasemax = 
        First[NMaximize[{logstep[phase[h]][s], 
              mindecade <= s <= maxdecade}, {s}, 
            Method -> "DifferentialEvolution"]];
      {a, b} = {a, b} /. 
          First@Solve[{dbmin a + b == 0, dbmax a + b == 100}]; {c, 
          d} = {c, d} /. 
          First@Solve[{phasemin c + d == 0, phasemax c + d == 100}];
      plotdb[h_][s_] := a db[h][s] + b;
      plotphase[h_][s_] := c phase[h][s] + d;
      
      work = LinearScale[dbmin, dbmax] /. {a_, b_} -> a;
      diff = Part[work, 2] - Part[work, 1];
      dbticksiterator = {First[work] - diff, Last[work] + diff, diff, 5};
      work = LinearScale[phasemin, phasemax] /. {a_, b_} -> a;
      diff = Part[work, 2] - Part[work, 1];
      phaseticksiterator = {First[work] - diff, Last[work] + diff, diff, 5};
      
      xticks[format_] := 
        CustomTicks[logvalue, {mindecade, maxdecade, {1}, Range[2, 9]}, 
          CTNumberFunction :> format];
      dbticks := CustomTicks[a #1 + b &, dbticksiterator]; 
      phaseticks := CustomTicks[c #1 + d &, phaseticksiterator]; 
      xgrids := 
        CustomGridLines[
          logvalue, {mindecade, maxdecade, Range[10]}, {Gainsboro}];
      ygrids := CustomGridLines[Identity, {0, 100, 10}, {Gainsboro}];
      
      Draw2D[{Draw[logstep[plotdb[h]][s], {s, mindecade, maxdecade}], Red, 
          Draw[logstep[plotphase[h]][s], {s, mindecade, maxdecade}] /. 
            SplitLineOnVertical[90]},
        opts,
        Frame -> True, 
        FrameTicks -> {xticks[nformat], dbticks, xticks["" &], phaseticks}, 
        GridLines -> {xgrids, ygrids}, 
        FrameLabel -> {s, "db", None, StyleForm["phase", FontColor -> Red]}, 
        PlotRange -> {-2, 102},
        PlotLabel -> SequenceForm["Bode Plot for ", hfunc[s]],
        Background -> Linen,
        ImageSize -> 500]];

Examples:

h[s_] = 100/(s + 20); 
BodePlot[h, {1, 3}]; 

h[s_] = (100*s + 100)/(s^2 + 110*s + 1000); 
BodePlot[h, {-2, 4}, 0, Background -> ColorMix[Mint, White][0.6]]; 

h[s_] = (10*(s + 10))/(s^2 + 3*s); 
BodePlot[h, {-1, 3}]; 

In the following plot we have used phasemode = 1 to force the phase function to lie between 0° and 360°. With phasemode = 0 the phase curve would have been discontinuous.

h[s_] = -((100*s)/(s^3 + 12*s^2 + 21*s + 10)); 
BodePlot[h, {-2, 3}, 1]; 

h[s_] = (30*(s + 10))/(s^2 + 3*s + 50); 
BodePlot[h, {-1, 3}]; 

h[s_] = (4*(s^2 + s + 25))/(s^3 + 100*s^2); 
BodePlot[h, {-1, 3}]; 

David Park
djmp at earthlink.net
http://home.earthlink.net/~djmp/ 




  • Prev by Date: Re: Nestwhile
  • Next by Date: Re: how to map function with over 1 arguments to list
  • Previous by thread: Re: Re: Clearing function definitions by argument type?
  • Next by thread: Re: Bode Plots in Mathematica