# Stability Boundaries via Routh-Hurwitz Criterion

I'm attempting to determine the stability boundaries of a 2nd order system via Routh-Hurwitz stability criterion. That is to say, I need to compute when a polynomial, which is in terms of variables "a" and "c", is equal to zero. This polynomial is very large. Where R is the system matrix,

R = GQ[7, X, Ginv];
cp = CharacteristicPolynomial[R, x];
cl = Reverse[CoefficientList[cp, x]];
h1 = cl[]; h1 > 0
h2 = cl[] cl[] - cl[];


Thus when h2>0 the system is stable. I first tried

RegionPlot[h2 > 0, {a, 0, 20}, {c, 0, 20}, Mesh -> All,
FrameLabel -> {"a", "c"}, PlotPoints -> 120] Attempt to numerically find boundaries failed miserably:

sol = Table[FindRoot[(h2 /. a -> n) == 0, {c, 0}], {n, 0, 3, 0.1}];
ListPlot[c /. sol, DataRange -> {0, 3}] FindRoot::lstol: "The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than \!$$TraditionalForm\MachinePrecision$$ digits of working precision to meet these tolerances. \!$$\*ButtonBox[\">>\", ButtonStyle->\"Link\", ButtonFrame->None, ButtonData:>\"paclet:ref/message/FindRoot/lstol\", ButtonNote -> \"FindRoot::lstol\"]$$"


I also tried

Solve[H2 == 0, c]
NSolve[H2 == 0, c]


but after quite some time of waiting, neither evaluated. Is it possible for Mathematica to determine these boundaries? RegionPlot provides a nice picture, but I'd like to have numerical values in order to know the exact curves of the boundaries.

• If you could provide sample GQ, Ginv, x and X that would be great. It would be nice to be able to recreate your graph or a similar one and then go from there. – C. E. Apr 11 '15 at 19:08
• @Pickett GQ[] is a Gaussian quadrature function I've written. I can include the matrix R. What is the preferred format of providing matrices? Can I upload it some way? – gKirkland Apr 12 '15 at 2:44

For simplicity, take as an example

h2 = c^2 + a^2;


Now, plot h == 0 instead of the region h > 0.

plt = ContourPlot[h2 == 1, {a, -1, 1}, {c, -1, 1}, FrameLabel -> {"a", "c"}]


and extract from it the Graphics elements that make up the curve.

pts = Cases[plt, GraphicsComplex[z_, __] :> z, Infinity]
(* {{{0.0350877, -0.999373}, {-2.22045*10^-16, -1.}, {-0.0357143, -0.999351}, ... *)


These are, I believe the boundary points you were seeking with FindRoot. The can be plotted with

ListPlot[pts, AspectRatio -> 1, AxesLabel -> {"a", "c"}] Your plot is, of course, much more complicated, with several disjoint regions. So, you will obtain multiple lists of points. For instance,

pltc = ContourPlot[Cos[c] + Cos[a] == 1/2, {a, 0, 4 Pi}, {c, 0, 4 Pi}];
ptsc = Cases[pltc, GraphicsComplex[z_, __] :> z, Infinity];
ListPlot[ptsc, AspectRatio -> 1, AxesLabel -> {"a", "c"}] Note: In general, use

plt//InputForm


to see what information is available inside a Graphics` object, but be prepared for a lot of output.