# Analyze stability of equilibria using Routh-Hurwitz conditions

For an assignment, I need to analyze the stability of a system very close to equilibrium, using "Routh-Hurwitz conditions". I have already obtained the characteristic equation of my system, but I do not know how to proceed further. The question is worded quite vaguely, but I suspect obtaining a Routh array is the goal here.

My system:

Df1 = D[r[t], t] == (g*w) - m*r[t] - b*r[t]*v[t];
Df2 = D[l[t], t] == (p*b*r[t]*-v[t]) - (m*l[t]) - a*l[t];
Df3 = D[e[t], t] == (1 - p)*(b*r[t]*v[t]) + (a*l[t]) - d*e[t];
Df4 = D[v[t], t] == pi*e[t] - s*v[t];


State equilibrium conditions, substitute equation for Epsilon (to define near-equilibrium conditions):

eql = Solve[{Df1[[2]] == 0, Df2[[2]] == 0, Df3[[2]] == 0, Df4[[2]] == 0} /. v[t] -> 0, {r[t], l[t], e[t]}] // First
varsubs = Thread[{v[t], r[t], l[t], e[t]} -> ({0, r[t], l[t], e[t]} + \[Epsilon] {\[Delta]v[t], \[Delta]r[t], \[Delta]l[t], \[Delta]e[t]} /. eql)]

devl = Series[Df2[[2]] /. varsubs, {\[Epsilon], 0, 1}] // Normal
deve = Series[Df3[[2]] /. varsubs, {\[Epsilon], 0, 1}] // Normal
devv = Series[Df4[[2]] /. varsubs, {\[Epsilon], 0, 1}] // Normal


Generate matrix of substituted solutions:

(M = D[{devv, deve, devl}, {{\[Delta]v[t], \[Delta]e[t], \[Delta]l[t]}}]) //Simplify // MatrixForm


Get characteristic function from matrix:

P1 = CharacteristicPolynomial[M, t]


Output of aforementioned commands:

Can I analyze this polynomial directly in Mathematica, i.e. forming a Routh array? Some insights would be greatly appreciated.

Kind regards.

• Does it really need to be Routh-Hurwitz? I'll see if I can post my routine for this later... (should prolly also do Schur-Cohn and Jury at some point) – J. M. is away Oct 8 '18 at 13:47
• The general question is simply "mathematical analysis of stability near-equilibrium conditions", Routh-Hurwitz was added as a hint after pestering the TA, so any input you have would be great! – Matton968 Oct 8 '18 at 14:16
• @J.M.issomewhatokay. I'd be curious to see your routine if you have a chance! – Chris K Oct 8 '18 at 16:59

There are "modified Routh-Hurwitz criteria" that work directly on the Jacobian matrix that I prefer and don't seem very well known (Fuller 1968). Here's a function I wrote that uses them for up to the 3x3 case, which is based on eqns. 12.21-12.23 from Fuller.

RouthHurwitzCriteria[a_?MatrixQ]:=Module[{c3},
If[!SquareMatrixQ[a],Message[RouthHurwitzCriteria::nonsq];Return[a]];
Which[
Length[a]==1,
Return[Piecewise[{{True,a[[1,1]]<0},{False,a[[1,1]]>0}},Indeterminate]],
Length[a]==2,
Return[Piecewise[{
{True,Tr[a]<0&&Det[a]>0},{False,Tr[a]>0||Det[a]<0}},Indeterminate]],
Length[a]==3,
c3=Det[{
{a[[1,1]]+a[[2,2]],a[[2,3]],-a[[1,3]]},
{a[[3,2]],a[[1,1]]+a[[3,3]],a[[1,2]]},
{-a[[3,1]],a[[2,1]],a[[2,2]]+a[[3,3]]}}];
Return[Piecewise[{
{True,Tr[a]<0&&Det[a]<0&&c3<0},
{False,Tr[a]>0||Det[a]>0||c3>0}},Indeterminate]],
Length[a]>3,
Message[RouthHurwitzCriteria::toobig];Return[a]
]
];


A more direct way to find your Jacobian matrix may be

M = D[{Df2[[2]], Df3[[2]], Df4[[2]]}, {{l[t], e[t], v[t]}}]


Then

Simplify[RouthHurwitzCriteria[M /. eql]]


gives the stability criteria.

I'd be happy if someone could implement the results from Fuller for larger matrices, but it looks significantly uglier.

Fuller, A. T. 1968. Conditions for a matrix to have only characteristic roots with negative real parts. Journal of Mathematical Analysis and Applications 23: 71–98.

• Thank you so much, this looks very promising and pretty neat. Do you mind if I use your solution, with credit? – Matton968 Oct 8 '18 at 17:01
• No prob, that's the point of the site :) – Chris K Oct 8 '18 at 17:02