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

Wed Feb 28 18:25:56 EST 1996