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

Help File for Hurwitz

FUNCTION: Hurwitz         --- decide whether or not a polynomial
                              has all its zeros strictly in the
                              left half plane.
FUNCTION: `Hurwitz/data`  --- Auxiliary routine for the above, 
                              may be called directly.
   
CALLING SEQUENCE:
  Hurwitz(p,z,'s','g');          --- usual call
  `Hurwitz/data`(p,z);           --- recommended if symbolic coefficients
                                     present, as type testing is not
                                     attempted.
   
PARAMETERS:
 -  p is a polynomial with complex coefficients, possibly with undefined
    parameters which evalc and hence this routine assumes to be real.
 -  z is the variable of the polynomial
 - 's' is the optional name in which the sequence of partial fractions
    of the Stieltjes continued fraction of (p - p*)/(p + p*) is returned.
 - 'g' is the optional name to return the gcd of p and its paraconjugate p*.
   
REFERENCE:   
 Complex Variables, by Norman Levinson and Raymond M. 
 Redheffer, Holden-Day 1970, pp 298-306.
   
SYNOPSIS:   
 The paraconjugate p* of p is defined as the polynomial whose roots
 are the roots of p reflected across the imaginary axis.
    
 A polynomial is a Hurwitz polynomial if all its roots are in the
 left half plane. 
     
`Hurwitz/data` computes the sequence of quotients in the Stieltjes continued
 fraction of the appropriate rational function to determine if p 
 is a Hurwitz polynomial.  It returns a list with two elements --- the
 first element is a list of the partial fractions in the Stieltjes continued
 fraction, and the second element is the gcd of p and p*.
   
 The first element of the returned sequence is special.  If it is
 of higher degree than 1 in z, p is not Hurwitz.  If it
 is of the form a + b*z, where Re(a) <> 0 or  b < 0, then p
 is not Hurwitz.  If each subsequent polynomial in the
 returned sequence is of the form a + b*z, where Re(a) = 0 and
 b > 0, then p is Hurwitz, and not otherwise. 
   
 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.
   
 The optional argument 's', if present, will contain on exit the
 partial fractions of the Stieltjes continued fraction so you can
 examine them yourself.  This is particularly useful if p has symbolic
 coefficients, and you can then decide the ranges of the coefficients
 that make p Hurwitz.
   
 The optional argument 'g', if present, will contain on exit the gcd
 of p and its paraconjugate.  The zeros of this gcd are precisely
 the zeros of p which are symmetrical under reflection across the 
 imaginary axis. 
   
WARNING: **********************************************************
         If the coefficients of p contain radicals or RootOfs, the
         code MAY work or it MAY not.  The problem is that radicals
         are converted to RootOfs to allow the GCD computation to
         proceed, and that if alpha is, say, RootOf(z^3 - z + 1,z)
         then evalc(conjugate(alpha)) is ANOTHER (or possibly the
         same) root of the same polynomial.  It is RECOMMENDED that
         symbols be used instead (e.g. use c or c + I*d instead of
         alpha) and then the resulting sequence should be tested at
         the numerical value of the parameter desired.
   
         This approach is problematical if division by the polynomial
         defining alpha occurs, but in that case the gcd should be
         nontrivial on specialization so one can still decide that the
         polynomial is not Hurwitz.
         ***********************************************************
   
EXAMPLES:   
> p1 := z^2 + z + 1;
                                       2
                                p1 := z  + z + 1
   
> Hurwitz(p1,z);
   
                                      true
   
> p2 := 3*z^3 + 2*z^2 + z + c;
                                    3      2
                           p2 := 3 z  + 2 z  + z + c
   
> Hurwitz(p2,z,'s2','g2');
                                      FAIL
   
> s2;   
                                  z            z (- 2 + 3 c)
                  [3/2 z, - 4 ---------, - 1/2 -------------]
                              - 2 + 3 c              c
   
> g2;   
                                       1
   
   
The elements of s2 are all positive if and only if 0 < c < 2/3, 
by inspection. Thus we can use the information returned even when 
the direct call to Hurwitz fails.
   
Separate calls to Hurwitz in the cases c = 0 and c = 2/3
give nontrivial gcd's between p2 and its paraconjugate, so we
conclude that the stability criteria are satisfied only as above.
   
   
> p3 := 4*z^4 + z^3 + z^2 + c;
                                     4    3    2
                            p3 := 4 z  + z  + z  + c
   
> Hurwitz(p3,z,'s3','g3');
                                      FAIL
   
> s3;   
                            [0, 4 z, z, - z/c, - z]
   
We see that the last term has coefficient -1, so in fact we can
say unequivocally that p3 is not Hurwitz, for any value of c.
   
   
> g3;   
                                       1
   
> p4 := z^5 + 5*z^4 + 4*z^3 + 3*z^2 + 2*z + c;
                           5      4      3      2
                    p4 := z  + 5 z  + 4 z  + 3 z  + 2 z + c
   
> Hurwitz(p4,z,'s4','g4');
                                      FAIL
   
> s4;   
                                                  2                      2
          25              z              (1 + 5 c)  z     (- 2 + 48 c + c ) z
 [1/5 z, ---- z, 289/5 -------, - 1/17 ---------------, - -------------------]
          17           1 + 5 c                       2        c (1 + 5 c)
                                       - 2 + 48 c + c
   
> g4;   
                                       1
   
Inspection of s4 leads us to believe that p4 is Hurwitz only if
c > -1/5, and -2 + 48*c + c^2 < 0, and c > 0.  This can be simplified
to the conditions 0 < c < -24 + 17*sqrt(2) = 0.04...
   
> p5 := p2 + I*d;
                                 3      2
                        p5 := 3 z  + 2 z  + z + c + I d
   
Evalc (and hence this routine) assumes that names have real values.
   
> Hurwitz(p5,z,'s5','g5');
                                      FAIL
   
> s5;   
                          I d             z
         [3/2 z, - 8 ------------ - 4 ---------,
                                2     - 2 + 3 c
                     (- 2 + 3 c)
   
                             2                                3
                  (- 2 + 3 c)  I d                 (- 2 + 3 c)  z
             ------------------------- - 1/2 -------------------------]
                       2      3      2                 2      3      2
             4 c - 12 c  + 9 c  - 8 d        4 c - 12 c  + 9 c  - 8 d
   
> g5;   
                                       1
    
   
The coefficients of s5 can be inspected according to the rules laid
down, but it is a tedious process.
   
> p6 := expand((x-1)*(x^2+2)*(x-c));
                    4    3        2            3    2
             p6 := x  - x  c + 2 x  - 2 x c - x  + x  c - 2 x + 2 c
   
> Hurwitz(p6,x,'s6','g6');
                                     false
   
> g6;   
                                      2
                                     x  + 2
   
> p7 := x + sqrt(2);
                                            1/2
                                 p7 := x + 2
   
> Hurwitz(p7,x);   
                                      FAIL
   
This failed because Maple was unable to immediately decide that
sqrt(2) is positive.  A simple check reveals that the test sqrt(2) > 0
fails, as perhaps it should (problems of this kind can be arbitrarily
difficult to decide, so Maple just does the simple ones automatically).
   
> p8 := x + sqrt(2 + sqrt(5));
                                            1/2 1/2
                            p8 := x + (2 + 5   )
   
> Hurwitz(p8,x);   
                                      FAIL
   
Likewise.  Examination of the sequence would allow the user to
decide if this was indeed Hurwitz or not.
   
> alias(alpha=RootOf(x^5 + x + 1,x));
                                   I, alpha
   
> p9 := x^3 + alpha*x^2 + (alpha^2-1)*x + 1;
                           3          2         2
                    p9 := x  + alpha x  + (alpha  - 1) x + 1
   
> Hurwitz(p9,x,'s9','g9');
Error, (in evala/indep/simple) reducible RootOf detected.  Substitutions are, {
alpha = RootOf(_Z**3-_Z**2+1), alpha = RootOf(_Z**2+_Z+1)}
   
This is the kind of message we get from evala(GCD(...)).  To deal with
it, make alpha one or the other of the suggested substitutions and re-
do the calculation.
   
> subs(alpha=RootOf(x^3-x+1,x),p9);
         3            3            2             3          2
        x  + RootOf(_Z  - _Z + 1) x  + (RootOf(_Z  - _Z + 1)  - 1) x + 1
   
> Hurwitz(",x,'s9','g9');
                                      FAIL
   
> s9;   
                                  2          3
                        x       %1  x     (%1  - %1 - 1) x
                     [----, ------------, ----------------]
                       %1     3                  %1
                            %1  - %1 - 1
   
                                     3
%1 :=                       RootOf(_Z  - _Z + 1)
   
> g9;   
                                       1
   
Examination of s9 with the various numerical values of the RootOf
will allow the user to decide what is happening.  A tricky point
is that some of the RootOf's have been conjugated, so for complex
values of the RootOf the above answer is almost certainly wrong.
   
> p10 := subs(alpha=c,p9);
                               3    2       2
                       p10 := x  + x  c + (c  - 1) x + 1
   
> Hurwitz(p10,x,'s10','g10');
                                      FAIL
   
> s10;   
                                 2         3
                                c  x     (c  - c - 1) x
                       [x/c, ----------, --------------]
                              3                 c
                             c  - c - 1
   
> g10;   
                                       1
   
   
Examination of the above (for real values of c) is a more
reliable way of telling whether or not the polynomial is Hurwitz.
   
> p11 := expand((c*z^2 + 1)*(z+1)*(z^2 + 2*z + 2));
                    5        4        3        2    3      2
          p11 := c z  + 3 c z  + 4 c z  + 2 c z  + z  + 3 z  + 4 z + 2
   
> Hurwitz(p11,z,'s11','g11');
                                      FAIL
   
> s11;   
                                       []
   
> g11;   
                                       2
                                    c z  + 1
   
"c" might be zero, so Hurwitz cannot tell if all the zeros are
in the left half plane or not.
   
   
This last example is a "forcing" of a very rare (perhaps impossible)
occurrence, to show what can happen.
   
> infolevel[Hurwitz] := 1;
                            infolevel[Hurwitz] := 1
   
We force the occurrence by saying what the Stieltjes continued
fraction will be (the empty list) and the gcd is (1).
   
> `Hurwitz/data`(p12,z) := [[],1];
                        Hurwitz/data(p12, z) := [[], 1]
   
> Hurwitz(p12,z,'s12','g12');
Hurwitz:   Continued fraction computed
Hurwitz:   Starting stability testing...
Hurwitz?:   Conditions vacuously satisfied --- check manually
                                      true
   
   
The above was insurance --- it should never happen.  However, if
it should turn out that the gcd is 1 while the sequence of partial
fractions is empty, then the conditions for being a Hurwitz polynomial
are vacuously satisfied.  A manual check is recommended, though this
warning will never be seen unless infolevel[Hurwitz] > 1.
   
   
SEE ALSO:  fsolve, the Jenkins-Traub package in the share library,
   and perhaps try the routine
   Roots := (p,z) -> evalf(Eigenvals(linalg[companion](p/lcoeff(p,z),z))):



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