# Determine stability of a system using the Routh-Hurwitz stability criterion

The Routh-Hurwitz is a criteria which serves to prove or disprove the stability of an electric control system.

# Idea

Given a system which has an equation of the form P(s)/Q(s) where P(s) and Q(s) are polynomials of any degree, it is said to be stable if all the roots of the polynomial Q(s) are in the left half of the complex plane, which means the real part of the root is negative.

For polynomials of degree less than or equal to 2, it is trivial to determine the roots and notice whether or not the system is stable. For polynomials Q(s) of higher degrees, the Routh–Hurwitz theorem can be useful to determine whether the system is Hurwitz-stable or not.

# Process

One of the methods to apply this theorem for any degree is described on Wikipedia that consists of performing operations on rows of a table, then counting the number of sign changes in the first column of the table. It goes as follows:

• Fill the first two rows with the coefficients of the polynomial, in the following pattern: 1st row consists of the even coefficients s^n, s^n-2, ..., s^2, s^0; 2nd row consists of the odd coefficients s^n-1, s^n-3, ..., s^3, s^1;
• Keep calculating the following rows by calculating the determinant of the matrix (last 2 rows of the first column, last 2 rows of the next column) divided by (the value of the first column on the row above); it is detailed in this part of the Wikipedia article.

• If you encounter a row full of zeros, set it to the derivative of the polynomial that the row above it represents. Depending on the row it is on, the number on each cell represents a coefficient. It is illustrated in the example later.

• If you encounter a zero with other non-zero values, the system is not stable: you can stop there. (You could continue but you'll get the same result). A zero on the first column would mean a root of real value 0, which renders the system unstable.

• Otherwise, continue until the end. Then, if not all the values in the first column are of the same sign, the system is not stable.

# Examples

Here are some examples of that process:

# Polynomial: 1 3 6 12 8.00

s^4  1 6  8
s^3  3 12 0
s^2  2 8  0
s^1--0-0--0-- < Row of zero - replace it by:
s^1  4 0  0   < Derivative 2*s^2 + 8*s^0 = 4*s^1
(between each 2 columns the power of s is lowered by 2)
s^0  8 0  0   < Continue normally

# Sign Changes: 0 (equivalent to no root in the right half plane)
# Stable: Yes


Another one:

s^2   1  -2
s^1   1   0
s^0  -2   0

# Sign Changes: 1 (equivalent to 1 root in the right half plane)
# Stable: No


# Rules

• The program will keep reading polynomials from standard input until a blank line, and output 0/1 for stable/non-stable respectively for each polynomial on a new line.

• The program can read the polynomial coefficients in any form. It could be 1 3 6 12 8 for the first example (it suits this problem better, but it doesn't really matter in the end).

• It has to handle cases where 0 in the first column, or 0 in the following columns, or a row full of 0 is encountered properly. (Note: just in case, no need to see how many roots are in rhp/lhp/jw axis).

• It can short-circuit the process at any time if a condition is met, there is no need to show the RH table in the end.

• You do not have to use this method. Any method which will determine whether or not a system is Hurwitz-stable can do the job.

• It doesn't matter which language you chose to do it with, as long as it's not done in one single built-in function call. I did mine in Python, I'll post it soon.

• This is a code golf. The shortest code wins.

Here's a sample input/output (given in the above mentioned way) to test your program with:

#input:
1 12 6 10 10 6 14
1 1 -2
1 2 1 2
1 3 1 3 16
1 2 3 6 2 1
1 2 14 88 200 800
1 2 3 6 2 1
1 3 6 12 8
1 4 5 2 4 16 20 8

#output:
0
0
1
0
0
0
0
1
0


If you need more clarifications, let me know, I'll do my best to clarify.

• How are the rows calculated? Jun 11 '13 at 14:01
• It is written on the wikipedia page in details Jun 11 '13 at 14:43
• It's a great example of a Wikipedia page which is only of any use to anyone who already knows the topic well and just needs to be reminded on some detail. Since you've recently studied the subject and written a computer program to do it, I'm sure you can write a description which is actually implementable. Jun 11 '13 at 15:35
• Followup thought: why not simplify the question by just asking for a program which takes real polynomials as inputs and outputs whether or not they are Hurwitz-stable (with suitable constraints to disallow one-token Mathematica solutions)? That way people can approach the problem using this tabular approach, a principal-minor-of-Sylvester-matrix approach, factorisation to irreducible polynomials over the reals, ... Allowing a variety of approaches tends to give more interesting challenges. Jun 11 '13 at 15:43
• @PeterTaylor Yes, you are right. Sorry, it's the first time I post here. I'll try to fix it now :) Jun 11 '13 at 17:22

p(a:b:s)=(a,b):p s
p[a]=[(a,0)]
p[]=[]
s((a:w),m@(b:v))=a:s(m,m%zipWith(\c d->(b*c-a*d)/b)w v)
s((a:_),_)=[a]
m%(0:v)|all(==0)v=fst$foldr(\b(d,p)->(b*p:d,2+p))([],0)m m%v=v q(a:s)|all((>0).(a*))s="1\n" q _="0\n" main=interact$(>>=q.s.unzip.p.map read.words).lines


Running the sample inputs:

& runhaskell 11847-Stability.hs
1 12 6 10 10 6 14
1 1 -2
1 2 1 2
1 3 1 3 16
1 2 3 6 2 1
1 2 14 88 200 800
1 2 3 6 2 1
1 3 6 12 8
1 4 5 2 4 16 20 8
0
0
1
0
0
0
0
1
0


(The output is interleaved with the input if you type it one line at a time, as per spec. If you bulk paste, you see all the output after the input.)

# Python (645 bytes)

That's what I managed to do (and golf):

from decimal import*
R=xrange
def s(p):
d=len(p)
if 0 in p:return 0
elif sum(1 for c in p if c>0) not in (0,d):return 0
p=p+ if d%2 else p[:]
r=[[0 for j in R(len(p)//2)] for i in R(d)]
r[:2]=[[v for v in p[i::2]] for i in (0,1)]
for i in R(2,len(r)):
a,h,b = 1,0,0
for k in R(len(r[i])-1):
v=(r[i-1]*r[i-2][k+1]-r[i-2]*r[i-1][k+1])/r[i-1]
r[i][k],a,h,b=v,a and v==0,h or v==0,b or (h and v!=0)
if a:r[i]=[v*max(0,d-i-2*k) for k,v in enumerate(r[i-1])]
elif b:return False
return sum(1 for v in r if v>0) in (0,d)
while 1:
r=raw_input()
if not r:break
print 1 if s([Decimal(x) for x in r.split()]) else 0


It takes as input the coefficients of the polynomials in decreasing order (highest power of s to lowest 0).

• You should probably add that this is Python 2.x Jun 7 '14 at 9:24
• I haven't actually tested this, but it seems like R=xrange can drop the x and 1 if s([Decimal(x) for x in r.split()]) else 0 can be replaced by int(s([Decimal(x) for x in r.split()])). Jun 8 '14 at 7:30
• return False can probably be return 0 and can definitely be return 0>1. Jun 8 '14 at 7:31
• 0>1. That is plain evil. Oct 6 '14 at 8:56