2
$\begingroup$

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:

Output

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

Kind regards.

$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – J. M. is away Oct 8 '18 at 13:47
  • 1
    $\begingroup$ 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! $\endgroup$ – Matton968 Oct 8 '18 at 14:16
  • $\begingroup$ @J.M.issomewhatokay. I'd be curious to see your routine if you have a chance! $\endgroup$ – Chris K Oct 8 '18 at 16:59
3
$\begingroup$

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.

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

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.