Hurwitz determinants

edit

They aren't built in yet, but you can implement them easily. For a single one:

def hurwitz_determinant(p, k):
    d = p.degree()
    a = list(reversed(p.list()))
    return matrix([[a[2*i-j+1] if 2*i-j+1 >= 0 and 2*i-j+1 <= d else 0 for j in range(k)] for i in range(k)]).det()

If you want all of them, then it's easier to construct the big Hurwitz matrix and take its leading principal minors:

def hurwitz_determinants(p):
    d = p.degree()
    a = list(reversed(p.list()))
    H = matrix([[a[2*i-j+1] if 2*i-j+1 >= 0 and 2*i-j+1 <= d else 0 for j in range(d)] for i in range(d)])
    return [H[0:k,0:k].det() for k in range(1,d+1)]

We can test it with a symbolic polynomial to check that we get the right formulas:

sage: K.<a0,a1,a2,a3,a4,a5> = PolynomialRing(QQ)
sage: R.<z> = PolynomialRing(K)
sage: p = a0*z^5 + a1*z^4 + a2*z^3 + a3*z^2 + a4*z + a5
sage: hurwitz_determinants(p)
[a1,
 a1*a2 - a0*a3,
 a1*a2*a3 - a0*a3^2 - a1^2*a4 + a0*a1*a5,
 a1*a2*a3*a4 - a0*a3^2*a4 - a1^2*a4^2 - a1*a2^2*a5 + a0*a2*a3*a5 + 2*a0*a1*a4*a5 - a0^2*a5^2,
 a1*a2*a3*a4*a5 - a0*a3^2*a4*a5 - a1^2*a4^2*a5 - a1*a2^2*a5^2 + a0*a2*a3*a5^2 + 2*a0*a1*a4*a5^2 - a0^2*a5^3]

The first three agree with Wikipedia, so probably it's correct.

We can also use it numerically:

sage: R.<z> = PolynomialRing(CC)
sage: p = 12*z^5 + 3*z^4 - 1*z^2 + 15
sage: hurwitz_determinants(p)
[3.00000000000000,
 12.0000000000000,
 528.000000000000,
 -32400.0000000000,
 -486000.000000000]

The stability criterion works:

sage: all(hd > 0 for hd in hurwitz_determinants(p))
False
sage: p.roots()
[(-1.08409079275927, 1),
 (-0.387674058792615 - 0.997697867303632*I, 1),
 (-0.387674058792615 + 0.997697867303632*I, 1),
 (0.804719455172253 - 0.599032525363998*I, 1),
 (0.804719455172253 + 0.599032525363998*I, 1)]

Note Routh improved on Hurwitz's criterion: Routh–Hurwitz stability criterion. You might also want to use the algorithms listed there.