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[[1]]; h1 > 0
h2 = cl[[1]] cl[[2]] - cl[[3]];

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]

enter image description here

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}]

enter image description here

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.

  • 2
    $\begingroup$ 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. $\endgroup$ – C. E. Apr 11 '15 at 19:08
  • $\begingroup$ @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? $\endgroup$ – 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"}]

Mathematica graphics

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"}]

Mathematica graphics

Note: In general, use


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.