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