Derivation of the Routh array explained

The Routh array is a tabular method permitting one to establish the stability of a system using only the coefficients of the characteristic polynomial. Central to the field of control systems design, the Routh–Hurwitz theorem and Routh array emerge by using the Euclidean algorithm and Sturm's theorem in evaluating Cauchy indices.

The Cauchy index

Given the system:

\begin{align} f(x) & {} = a ^{n -1} _{1 x} + … + a _{n} & {} ( 1 ) \\ & {} = (x - r _{1)(x-r} _{2) … (x-r} _{n)} & {} ( 2 ) \\ \end{align} Assuming no roots of f(x) = 0 lie on the imaginary axis, and lettingN = The number of roots of f(x) = 0 with negative real parts, andP = The number of roots of f(x) = 0 with positive real partsthen we haveN + P = n ( 3 ) Expressing f(x) in polar form, we havef(x) = \rho(x)e ^{j\theta(x)} ( 4 ) where\rho(x) = \sqrt{ ak{Re} ^{2[f(x)]+ak{Im}} ^{2[f(x)]}} ( 5 ) and \theta(x) = \tan ^{-1} ( ak{Im}[f(x)]/ ak{Re}[f(x)] ) ( 6 ) from (2) note that\theta(x) = \theta _{r _{1}} (x) + \theta _{r _{2}} (x) + … + \theta _{r _{n}} (x) ( 7 ) where\theta _{r _{i}} (x) = \angle(x - r _{i)} ( 8 ) Now if the i^th root of f(x) = 0 has a positive real part, then (using the notation y=(RE[y],IM[y]))\begin{align} \theta _{r _{i}} (x) | _{x=-jinfty} & = \angle(x - r _{i)|} _{x=-jinfty} \\ & = \angle( 0- ak{Re}[r _{i],-infty-ak{Im}[r} _{i])} \\ & = \angle( - | ak{Re}[r _{i]|,-infty)} \\ & = \pi + \lim _{\phi} \tan ^{-1} \phi = (3 \pi)/(2) ( 9 )\\ \end{align} and\theta _{r _{i}} (x) | _{x=j0} = \angle( - | ak{Re}[r _{i]|,0)} = \pi - \tan ^{-1} 0= \pi ( 10 ) and\theta _{r _{i}} (x) | _{x=jinfty} = \angle( - | ak{Re}[r _{i]|,infty)} = \pi - \lim _{\phi} \tan ^{-1} \phi = (\pi)/(2) ( 11 ) Similarly, if the i^th root of f(x) =0 has a negative real part,\begin{align} \theta _{r _{i}} (x) | _{x=-jinfty} & = \angle(x - r _{i)|} _{x=-jinfty} \\ & = \angle( 0- ak{Re}[r _{i],-infty-ak{Im}[r} _{i])} \\ & = \angle(| ak{Re}[r _{i]|,-infty)} \\ & = 0 - \lim _{\phi} \tan ^{1} \phi =- (\pi)/(2) ( 12 )\\ \end{align} and\theta _{r _{i}} (x) | _{x=j0} = \angle(| ak{Re}[r _{i]|,0)} = \tan ^{-1} 0=0 ( 13 ) and\theta _{r _{i}} (x) | _{x=jinfty} = \angle(| ak{Re}[r _{i]|,infty)} = \lim _{\phi} \tan ^{-1} \phi = (\pi)/(2) ( 14 ) From (9) to (11) we find that \theta _{r _{i}} (x) | ^{x = j inf ty} _{x =- j inf ty} = - \pi when the i^th root of f(x) has a positive real part, and from (12) to (14) we find that \theta _{r _{i}} (x) | ^{x = j inf ty} _{x =- j inf ty} = \pi when the i^th root of f(x) has a negative real part. Thus,\theta(x) | ^{x = j inf ty} _{x =- j inf ty} = \angle(x - r _{1)|} ^{x = j inf ty} ^{x =- j inf ty} + \angle(x - r _{2)|} ^{x = j inf ty} ^{x =- j inf ty} + … + \angle(x - r _{n)|} ^{x = j inf ty} ^{x =- j inf ty} = \pi N - \pi P ( 15 ) So, if we define\Delta = (1)/(\pi) \theta(x) | ^{j inf ty} _{- j inf ty} ( 16 ) then we have the relationshipN - P = \Delta ( 17 ) and combining (3) and (17) gives usN = (n + \Delta)/(2) and P = (n - \Delta)/(2) ( 18 ) Therefore, given an equation of f(x) of degree n we need only evaluate this function \Delta to determine N , the number of roots with negative real parts and P , the number of roots with positive real parts.In accordance with (6) and Figure 1, the graph of \tan(\theta) vs \theta , varying x over an interval (a,b) where \theta _{a=\theta(x)|} _{x=ja} and \theta _{b=\theta(x)|} _{x=jb} are integer multiples of \pi , this variation causing the function \theta(x) to have increased by \pi , indicates that in the course of travelling from point a to point b, \theta has "jumped" from +inf ty to -inf ty one more time than it has jumped from -inf ty to +inf ty . Similarly, if we vary x over an interval (a,b) this variation causing \theta(x) to have decreased by \pi , where again \theta is a multiple of \pi at both x = ja and x = jb , implies that \tan \theta (x) = ak{Im}[f(x)]/ ak{Re}[f(x)] has jumped from -inf ty to +inf ty one more time than it has jumped from +inf ty to -inf ty as x was varied over the said interval.

Thus,

\theta(x) | ^{j inf ty} _{- j inf ty} is \pi times the difference between the number of points at which ak{Im}[f(x)]/ ak{Re}[f(x)] jumps from -inf ty to +inf ty and the number of points at which ak{Im}[f(x)]/ ak{Re}[f(x)] jumps from +inf ty to -inf ty as x ranges over the interval ( - j inf ty, + j inf ty) provided that at x = \pm j inf ty , \tan[\theta(x)] is defined.In the case where the starting point is on an incongruity (i.e. \theta _{a=\pi/2} \pm i\pi , i = 0, 1, 2, ...) the ending point will be on an incongruity as well, by equation (17) (since N is an integer and P is an integer, \Delta will be an integer). In this case, we can achieve this same index (difference in positive and negative jumps) by shifting the axes of the tangent function by \pi/ 2 , through adding \pi/ 2 to \theta . Thus, our index is now fully defined for any combination of coefficients in f(x) by evaluating \tan[\theta] = ak{Im}[f(x)]/ ak{Re}[f(x)] over the interval (a,b) = ( + j inf ty, - j inf ty) when our starting (and thus ending) point is not an incongruity, and by evaluating\tan[\theta'(x)] = \tan[\theta + \pi/ 2 ] = - \cot[\theta(x)] = - ak{Re}[f(x)]/ ak{Im}[f(x)] ( 19 ) over said interval when our starting point is at an incongruity.This difference, \Delta , of negative and positive jumping incongruities encountered while traversing x from - j inf ty to + j inf ty is called the Cauchy Index of the tangent of the phase angle, the phase angle being \theta(x) or \theta'(x) , depending as \theta _{a} is an integer multiple of \pi or not.

The Routh criterion

To derive Routh's criterion, first we'll use a different notation to differentiate between the even and odd terms of

f(x) :f(x) = a ^{n} _{0 x} + b ^{n -1} _{0 x} + a ^{n -2} _{1 x} + b ^{n -3} _{1 x} + … ( 20 ) Now we have: \begin{align} f(j\omega) & = a ^{n -1} _{0 (j\omega)} + a ^{n -2} _{1 (j\omega)} + b ^{n -3} _{1 (j\omega)} + … & {} ( 21 )\\ & = a ^{n -2} _{1 (j\omega)} + a ^{n -4} _{2 (j\omega)} + … & {} ( 22 )\\ & + b ^{n -1} _{0 (j\omega)} + b ^{n -3} _{1 (j\omega)} + b ^{n -5} _{2 (j\omega)} + … \\ \end{align} Therefore, if n is even, \begin{align} f(j\omega) & = ( -1 ) ^{n/2} [ a ^{n -2} _{1 \omega} + a ^{n -4} _{2 \omega} - … ] & {} ( 23 )\\ & + j( -1 ) ^{(n/2)-1} [ b ^{n -1} _{0 \omega} - b ^{n -3} _{1 \omega} + b ^{n -5} _{2 \omega} - … ] & {} \\ \end{align} and if n is odd:\begin{align} f(j\omega) & = j( -1 ) ^{(n-1)/2} [ a ^{n -2} _{1 \omega} + a ^{n -4} _{2 \omega} - … ] & {} ( 24 )\\ & + ( -1 ) ^{(n-1)/2} [ b ^{n -1} _{0 \omega} - b ^{n -3} _{1 \omega} + b ^{n -5} _{2 \omega} - … ] & {}\\ \end{align} Now observe that if n is an odd integer, then by (3) N + P is odd. If N + P is an odd integer, then N - P is odd as well. Similarly, this same argument shows that when n is even, N - P will be even. Equation (15) shows that if N - P is even, \theta is an integer multiple of \pi . Therefore, \tan(\theta) is defined for n even, and is thus the proper index to use when n is even, and similarly \tan(\theta') = \tan(\theta + \pi) = - \cot(\theta) is defined for n odd, making it the proper index in this latter case.

Thus, from (6) and (23), for

n even:\Delta = I ^{+inf ty} _{-inf ty} (- ak{Im)/([f(x)]}{ ak{Re}[f(x)]}= I) ^{+inf ty} ^{-inf ty} (b ^{n -1} _{0 \omega} - b ^{n -3} _{1 \omega} + …)/(a ^{n -2} _{1 \omega} + \ldots) ( 25 ) and from (19) and (24), for n odd:\Delta = I ^{+inf ty} _{-inf ty} (ak{Re)/([f(x)]}{ ak{Im}[f(x)]}= I) ^{+inf ty} ^{-inf ty} (b ^{n -1} _{0 \omega} - b ^{n -3} _{1 \omega} + \ldots)/(a ^{n -2} _{1 \omega} + \ldots) ( 26 ) Lo and behold we are evaluating the same Cauchy index for both:\Delta = I ^{+inf ty} _{-inf ty} (b ^{n -1} _{0 \omega} - b ^{n -3} _{1 \omega} + \ldots)/(a ^{n -2} _{1 \omega} + \ldots) ( 27 )

Sturm's theorem

Sturm gives us a method for evaluating

\Delta = I ^{+inf ty} _{-inf ty} (f _{2(x)})/(f _{1(x)}) . His theorem states as follows:

Given a sequence of polynomials

f _{1(x),f} _{2(x),} ... , f _{m(x)} where:

1) If

f _{k(x)} = 0 then f _{k-1} (x) ≠ 0 , f _{k+1} (x) ≠ 0 , and \operatorname{sign}[f _{k-1} (x)] = - \operatorname{sign}[f _{k+1} (x)]

2)

f _{m(x)} ≠ 0 for -inf ty < x < inf ty

and we define

V(x) as the number of changes of sign in the sequence f _{1(x),f} _{2(x),} ... , f _{m(x)} for a fixed value of x , then:\Delta = I ^{+inf ty} _{-inf ty} (f _{2(x)})/(f _{1(x)}) = V( -inf ty) - V( +inf ty) ( 28 ) A sequence satisfying these requirements is obtained using the Euclidean algorithm, which is as follows:

Starting with

f _{1(x)} and f _{2(x)} , and denoting the remainder of f _{1(x)/f} _{2(x)} by f _{3(x)} and similarly denoting the remainder of f _{2(x)/f} _{3(x)} by f _{4(x)} , and so on, we obtain the relationships:\begin{align} &f _{1(x)=} q _{1(x)f} _{2(x)} - f _{3(x)} ( 29 )\\ &f _{2(x)=} q _{2(x)f} _{3(x)} - f _{4(x)} \\ & \ldots \\ &f _{m-1} (x) = q _{m-1} (x)f _{m(x)} \\ \end{align} or in general f _{k-1} (x) = q _{k-1} (x)f _{k(x)} - f _{k+1} (x) where the last non-zero remainder, f _{m(x)} will therefore be the highest common factor of f _{1(x),f} _{2(x),} ... , f _{m-1} (x) . It can be observed that the sequence so constructed will satisfy the conditions of Sturm's theorem, and thus an algorithm for determining the stated index has been developed.

It is in applying Sturm's theorem (28) to (29), through the use of the Euclidean algorithm above that the Routh matrix is formed.

We get

f _{3(\omega)} = (a _{0})/(b _{0}) f _{2(\omega)} - f _{1(\omega)} ( 30 ) and identifying the coefficients of this remainder by c _{0} , - c _{1} , c _{2} , - c _{3} , and so forth, makes our formed remainder f _{3(\omega)} = c ^{n -2} _{0 \omega} - c ^{n -4} _{1 \omega} + c ^{n -6} _{2 \omega} - … ( 31 ) wherec _{0} = a _{1} - (a _{0})/(b _{0}) b _{1} = (b _{0a} _{1} - a _{0b} _{1})/(b _{0}) ; c _{1} = a _{2} - (a _{0})/(b _{0}) b _{2} = (b _{0a} _{2} - a _{0b} _{2})/(b _{0}) ;\ldots ( 32 ) Continuing with the Euclidean algorithm on these new coefficients gives usf _{4(\omega)} = (b _{0})/(c _{0}) f _{3(\omega)} - f _{2(\omega)} ( 33 ) where we again denote the coefficients of the remainder f _{4(\omega)} by d _{0} , - d _{1} , d _{2} , - d _{3} ,making our formed remainder f _{4(\omega)} = d ^{n -3} _{0 \omega} - d ^{n -5} _{1 \omega} + d ^{n -7} _{2 \omega} - … ( 34 ) and giving usd _{0} = b _{1} - (b _{0})/(c _{0}) c _{1} = (c _{0b} _{1} - b _{0c} _{1})/(c _{0}) ; d _{1} = b _{2} - (b _{0})/(c _{0}) c _{2} = (c _{0b} _{2} - b _{0c} _{2})/(c _{0}) ;\ldots ( 35 ) The rows of the Routh array are determined exactly by this algorithm when applied to the coefficients of (20). An observation worthy of note is that in the regular case the polynomials f _{1(\omega)} and f _{2(\omega)} have as the highest common factor f _{n+1} (\omega) and thus there will be n polynomials in the chain f _{1(x),f} _{2(x),} ... , f _{m(x)} .

Note now, that in determining the signs of the members of the sequence of polynomials

f _{1(x),f} _{2(x),} ... ,f _{m(x)} that at \omega = \pm inf ty the dominating power of \omega will be the first term of each of these polynomials, and thus only these coefficients corresponding to the highest powers of \omega in f _{1(x),f} _{2(x),} ... , and f _{m(x)} , which are a _{0} , b _{0} , c _{0} , d _{0} , ... determine the signs of f _{1(x)} , f _{2(x)} , ..., f _{m(x)} at \omega = \pm inf ty .

So we get

V( +inf ty) = V(a _{0,} b _{0,} c _{0,} d _{0,} ... ) that is, V( +inf ty) is the number of changes of sign in the sequence a ^{n} _{0inf ty} , b ^{n -1} _{0inf ty} , c ^{n -2} _{0inf ty} , ... which is the number of sign changes in the sequence a _{0} , b _{0} , c _{0} , d _{0} , ... and V( -inf ty) = V(a _{0,} - b _{0,} c _{0,} - d _{0,} ...) ; that is V( -inf ty) is the number of changes of sign in the sequence a ^{n} _{0 ( -inf ty)} , b ^{n -1} _{0 ( -inf ty)} , c ^{n -2} _{0 ( -inf ty)} , ... which is the number of sign changes in the sequence a _{0} , - b _{0} , c _{0} , - d _{0} , ...

Since our chain

a _{0} , b _{0} , c _{0} , d _{0} , ... will have n members it is clear that V( +inf ty) + V( -inf ty) = n since within V(a _{0,} b _{0,} c _{0,} d _{0,} ... ) if going from a _{0} to b _{0} a sign change has not occurred, within V(a _{0,} - b _{0,} c _{0,} - d _{0,} ... ) going from a _{0} to - b _{0} one has, and likewise for all n transitions (there will be no terms equal to zero) giving us n total sign changes.

As

\Delta = V( -inf ty) - V( +inf ty) and n = V( +inf ty) + V( -inf ty) , and from (18) P = (n - \Delta/ 2 ) , we have that P = V( +inf ty) = V(a _{0,} b _{0,} c _{0,} d _{0,} ... ) and have derived Routh's theorem -

The number of roots of a real polynomial

f(z) which lie in the right half plane ak{Re}(r _{i)} > 0 is equal to the number of changes of sign in the first column of the Routh scheme.

And for the stable case where

P = 0 then V(a _{0,} b _{0,} c _{0,} d _{0,} ... ) = 0 by which we have Routh's famous criterion:

In order for all the roots of the polynomial

f(z) to have negative real parts, it is necessary and sufficient that all of the elements in the first column of the Routh scheme be different from zero and of the same sign.

References

• Hurwitz, A., "On the Conditions under which an Equation has only Roots with Negative Real Parts", Rpt. in Selected Papers on Mathematical Trends in Control Theory, Ed. R. T. Ballman et al. New York: Dover 1964
• Routh, E. J., A Treatise on the Stability of a Given State of Motion. London: Macmillan, 1877. Rpt. in Stability of Motion, Ed. A. T. Fuller. London: Taylor & Francis, 1975
• Felix Gantmacher (J.L. Brenner translator) (1959) Applications of the Theory of Matrices, pp 177–80, New York: Interscience.