A complete discussion of the mathematical theory behind the approach taken in this book can be found for example in . The following abridged version should give enough information to use the given routine reliably. This routine is available from the Maple share library.
A polynomial is a Hurwitz polynomial if all its roots are in the left half plane. The paraconjugate of p is defined as the polynomial whose roots are the roots of p reflected across the imaginary axis. Algebraically, . So if , then .
The auxiliary Maple subroutine
`Hurwitz/data` computes the sequence of partial
quotients in the Stieltjes continued
fraction of the rational function
to determine if for , which will happen if and only if p is Hurwitz. The subroutine returns a list with two elements --- the first element of the list is another list, of the partial quotients in the Stieltjes continued fraction, and the second element of the first list is the gcd of p and . If this gcd is nontrivial, then there are pairs of roots symmetric on reflection through the imaginary axis, or else purely imaginary roots, and in either case p is not Hurwitz.
We then compute the partial quotients of the Stieltjes continued fraction for essentially by the Euclidean algorithm. We put where q is the quotient on polynomial division of by , and put and repeat the process with the resulting numerator and denominator. We terminate the process when we encounter a zero remainder.
For example, consider . We see and since the coefficients of p are real,
Then the first step of the algorithm gives and .
Replacing h by and starting again, we get (which is our first partial quotient) and , so . Repeating the process twice more gives quotients and . Since not all of these are positive, the original polynomial is not Hurwitz.
The Maple routine
`Hurwitz/data` does all this calculation
automatically and returns a list containing the GCD of p and
and a sub-list of the partial quotients.
The first element of the returned list of partial quotients
is special. If it is
of higher degree than 1 in , p is not Hurwitz. If it
is of the form , where or b < 0, then p
is not Hurwitz. If each subsequent polynomial in the
returned sequence is of the form , where and
b > 0, then p is Hurwitz, and not otherwise.
See  for details.
The Maple routine
will return true if it can use these rules to
decide if p is Hurwitz, false if it can decide that p is not
Hurwitz, and FAIL otherwise, in which case it is up to the user to
decide when p is Hurwitz, given the partial quotients. Of course,
if the polynomial has symbolic coefficients then there will usually
be no way that this routine can tell automatically when the
polynomial is Hurwitz, and thus the described mechanism for returning
the partial quotients to the user for inspection is required.
Consider the matrix
where c = a + ib. The characteristic polynomial of this matrix is
We can easily tell directly when the roots of this polynomial are in the left half plane as it is linear in c and hence a and b. Putting gives us the parametric equations and describing the location in parameter space where roots of p cross the imaginary axis. See Figure 5.1.
Figure 1: Values of a and b where the roots of have purely imaginary roots.
This parametric curve is symmetric and has a loop. Inspection shows that all roots are in the left half plane precisely when c = a+ib has values inside the loop.
We now show how our implementation of the Hurwitz criterion is used to determine the same information. We load in the Hurwitz package from the Maple share library and call the routine according to the instructions in the help file. The immediate result is the answer FAIL --- the Maple routine cannot decide by itself if the polynomial is Hurwitz or not, as this depends on the value of the parameter c, which Maple knows nothing about. This is the expected behaviour, and we then examine the gcd of p and , which turns out to be 1, and the sequence of partial quotients computed, which turns out to be
Now we see that the first partial quotient has positive real part (1/6) for the coefficient of and zero real part for the constant term (in fact zero imaginary part, also). The second partial quotient will have positive real part if and only if a > -56. The third partial quotient will have positive real part if and only if
(one sees this by isolating and factoring the right-hand side). Plotting of this region gives exactly the interior of the loop shown earlier.
The Maple calculation is as below.
> p := lambda^3 + 6*lambda^2 + 10*lambda + 4 - a - I*b; 3 2 p := lambda + 6 lambda + 10 lambda + 4 - a - I b > with(share): > readlib(Hurwitz): > Hurwitz(p,lambda); FAIL > Stieltjes_continued_fraction_and_GCD := `Hurwitz/data`(p,lambda); Stieltjes_continued_fraction_and_GCD := [ I b lambda [1/6 lambda, 216 --------- + 36 ------, 2 56 + a (56 + a) 2 I (56 + a) b --------------------------------------- 2 3 2 - 12544 + 2688 a + 108 a + a + 216 b 3 (56 + a) lambda - 1/6 ---------------------------------------], 1] 2 3 2 - 12544 + 2688 a + 108 a + a + 216 bThe 1 on the end is the GCD of p and its paraconjugate.