next up previous
Next: Source Code for Up: Stability of Polynomials Previous: Stability of Polynomials

Mathematical Background.

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,

and so

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 Hurwitz 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);

> 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)

                             I (56 + a)  b
                                        2    3        2
                - 12544 + 2688 a + 108 a  + a  + 216 b

                                       (56 + a)  lambda
                     - 1/6 ---------------------------------------], 1]
                                                   2    3        2
                           - 12544 + 2688 a + 108 a  + a  + 216 b
The 1 on the end is the GCD of p and its paraconjugate.


  1. A simple necessary condition for a real polynomial to be Hurwitz is that all its coefficients be positive. Find a degree three polynomial with positive coefficients that is not Hurwitz.
  2. For which values of c is Hurwitz?
  3. If is the truncation to degree n of the Taylor series for about x=0, show that is Hurwitz for k = 0, 1, 2, 3, 4, 5 but not for k = 6. In fact it is not Hurwitz for any k > 5 [].
  4. If is the Bessel polynomial of degree n, show that it is stable for n = 5, 10, 15. It can be shown [] that these polynomials are stable for all n.

next up previous
Next: Source Code for Up: Stability of Polynomials Previous: Stability of Polynomials

Robert Corless
Wed Feb 28 18:25:56 EST 1996