CITE THIS NOTEBOOK: The Routh-Hurwitz stability criterion by Paul Abbott. Wolfram Community MAR 4 2023.
The Routh–Hurwitz stability criterion gives a necessary and sufficient stability condition for a control system. A polynomial
p
is stable—that is, all its roots have strictly negative real part—if and only if all the principal minors
Δ
n
of the Hurwitz matrix corresponding to
p
are positive. An application to optimising the ground-state cooling of a micromechanical oscillator, which involves Hurwitz Integrals, is presented.
TraditionalForm
input and output is used throughout so that the computational expressions most closely correspond to standard typeset notation.

Hurwitz Polynomials and Integrals

A Hurwitz polynomial is a polynomial with real positive coefficients and roots that are either negative or pairwise conjugate with negative real parts.
The integral
1
2π
∞
∫
-∞
(z)
(z)(-z)
z
, where
(z)
is Hurwitz can be computed using the residue theorem: if
g
and
h
are analytic at
z
i
and
g(
z
i
)≠0
,
h(
z
i
)=0
, and
′
h
(
z
i
)≠0
, then
f(z)=
g(z)
h(z)
has a simple pole at
z
i
, the residue of
f
at
z
i
is
Res(f,
z
i
)=
g(
z
i
)
′
h
(
z
i
)
, and
1
2π
∫
γ
g(z)
h(z)
z=
∑
i
g(
z
i
)
′
h
(
z
i
)
, where γ is a simple closed curve enclosing the
z
i
.
Hence, the integral reduces to the sum
∑
i
(
z
i
)
′

(
z
i
)(-
z
i
)
over the roots
z
i
of
(z)
and one can use
RootSum
to compute the integral.

Definitions

Polynomials

Define

n
(z)
using dynamic programming:
/:

n_
:=/:

n
=zEvaluate
n
∑
i=1
p
i
2(n-i)
z

Compute

4
(z)
:

4
(z)
p
1
6
z
+
p
2
4
z
+
p
3
2
z
+
p
4
Define

n
(z)
using dynamic programming:
/:

n_
:=/:

n
=zEvaluate
n
∑
i=0
q
i
n-i
z

Compute

4
(z)
:

4
(z)
q
0
4
z
+
q
1
3
z
+
q
2
2
z
+
q
3
z+
q
4

Stability

The best known criterion of stability is the Routh–Hurwitz criterion.
Compute the Hurwitz (or Routh–Hurwitz) matrix for

4
(z)
:
mat4=
[◼]
HurwitzMatrix
[

4
[z],z]
q
1
q
3
0
0
q
0
q
2
q
4
0
0
q
1
q
3
0
0
q
0
q
2
q
4
A polynomial
p
is stable if and only if the principal minors
Δ
n
of the Hurwitz matrix of
p
are all positive.
Compute the principal minors,
Δ
n
, using
Det
:
Δ
n_
:=ℳDet[Take[ℳ,n,n]]
Determine
Δ
n
for

4
(z)
:
Table[
Δ
n
[mat4],{n,4}]

q
1
,
q
1
q
2
-
q
0
q
3
,-
q
4
2
q
1
+
q
2
q
3
q
1
-
q
0
2
q
3
,-
q
0
q
4
2
q
3
+
q
1
q
2
q
4
q
3
-
2
q
1
2
q
4

Determine the conditions under which

4
(z)
is stable:
Reduce[And@@Thread[Join[%,Variables[mat4]]>0]]
q
0
>0∧
q
1
>0∧
q
2
>0∧0<
q
3
<
q
1
q
2
q
0
∧0<
q
4
<
q
1
q
2
q
3
-
q
0
2
q
3
2
q
1
With positive coefficients, only every second determinant is required (Liénard and Chippart criterion):
StabilityConditions[p_]:=zModulerh=
[◼]
HurwitzMatrix
[p[z],z],n=Exponent[p[z],z],Reduce[And@@Thread[Join[Table[
Δ
k
[rh],{k,n-1,1,-2}],Variables[rh]]>0],BacksubstitutionTrue]
Alternative determination of the conditions under which

4
(z)
is stable:
StabilityConditions[

4
][z]
q
0
>0∧
q
1
>0∧
q
2
>0∧0<
q
3
<
q
1
q
2
q
0
∧0<
q
4
<
q
1
q
2
q
3
-
q
0
2
q
3
2
q
1

Computation

Computation via RootSum

Compute the integral
ℐ
n
:=
1
2π
∞
∫
-∞

n
(z)

n
(z)

n
(-z)
z
using
RootSum
:
ℐ
n_
:=
ℐ
n
=RootSumzEvaluate@

n
(z),zEvaluate

n
(z)

n
(-z)
′

n
(z)
//Factor
Determine
ℐ
4
:
ℐ
4
p
4
q
3
2
q
0
-
p
4
q
1
q
2
q
0
+
p
3
q
1
q
4
q
0
-
p
2
q
3
q
4
q
0
-
p
1
q
1
2
q
4
+
p
1
q
2
q
3
q
4
2
q
0
q
4

q
4
2
q
1
-
q
2
q
3
q
1
+
q
0
2
q
3


Computation via Determinant of Hurwitz matrix

Observe that the denominator is proportional to the determinant of the Hurwitz matrix for

4
(z)
:
-2
q
0
q
1
q
3
0
0
q
0
q
2
q
4
0
0
q
1
q
3
0
0
q
0
q
2
q
4
//Factor
2
q
0
q
4

q
4
2
q
1
-
q
2
q
3
q
1
+
q
0
2
q
3

Also observe that the numerator is the determinant of matrix obtained from the Hurwitz matrix by inserting the coefficients of

4
(z)
into the first row:
p
1
p
2
p
3
p
4
q
0
q
2
q
4
0
0
q
1
q
3
0
0
q
0
q
2
q
4
p
4
q
3
2
q
0
-
p
4
q
1
q
2
q
0
+
p
3
q
1
q
4
q
0
-
p
2
q
3
q
4
q
0
-
p
1
q
1
2
q
4
+
p
1
q
2
q
3
q
4
Check that
ℐ
4
is recovered:
-
1
2
q
0
p
1
p
2
p
3
p
4
q
0
q
2
q
4
0
0
q
1
q
3
0
0
q
0
q
2
q
4
q
1
q
3
0
0
q
0
q
2
q
4
0
0
q
1
q
3
0
0
q
0
q
2
q
4

ℐ
4
//Factor
True
Finding this elegant and simple result, one would expect that it is well known; E.g., see Puri and Weygandt (1964) and, even earlier,
Kac (1952)
:
Let
P(p)/Q(p)
be the Laplace transform of
f(t)
, where
Q(p)=
a
0
n
p
+⋯+
a
n
and
P(p)
is a polynomial of degree
≤n-1
; let
P(p)P(-p)
g
0
2n−2
p
+⋯+
g
n-2
2
p
+
g
n-1
. The author proves that
J=
∞
∫
0
2
f(t)
t
n-1
(-1)
G/(2
a
0
D)
, where
D
a
1
a
3
a
5
⋯
0
a
0
a
2
a
4
⋯
0
0
a
1
a
3
⋯
0
0
a
0
a
2
⋯
0
⋮
⋮
⋮
⋱
⋮
0
⋯
⋯
a
n-2
a
n
and
G
is obtained by replacing the first row of
D
by the row
g
0
,⋯,
g
n-1
. A convenient method for the numerical computation of
G
is also given.

General Case

Compute
1
2π
∞
∫
-∞
(z)
(z)(-z)
z
for polynomials
(z)
and
(z)
, with
deg((z))≤n-2
, using
Det
:
ℐ[p_,q_]:=zModule{n=Exponent[q(z),z],ps,pmat,qmat,c},​​c=Coefficient[q(z),
n
z
];​​ps=PadLeft[CoefficientList[p(z),
2
z
]//Reverse,n];​​qmat=
[◼]
HurwitzMatrix
[q(z),z];​​pmat=ReplacePart[qmat,1ps];​​
n+1
(-1)
2c
pmat
qmat
//Factor

Application: Micromechanical Oscillator

Genes et al. (2008) provides a general framework to describe cooling of a micromechanical oscillator to its quantum ground state by means of radiation-pressure coupling with a driven optical cavity.
Consider equations (13)-(17):
along with equations (21)-(28):
The integrals in equation (13) involve Hurwitz polynomials so are straightforward to compute.

Implement equations (15)-(17)

Define the effective susceptibility, equation (17):
Define the (simplified) thermal noise spectrum, equation (15) simplified in equation (21):
Define the radiation pressure noise spectrum, equation (16):

Integral computation via Hurwitz Polynomials

A necessary condition for a polynomial to be Hurwitz is that it has positive real coefficients:
The constant coefficient in (z) (the first coefficient above) leads to the second stability condition:
Define the two stability conditions:
Define equations (23)-(25):