# Routh criterion calculator with symbolic variables

*Author: Rafael Ribeiro*


```python
from sympy import *
```


```python
def routh_criterion(coeffs):
    """
    Create a Routh-Hurwitz stability criterion matrix of a nth-degree polynomial 
    
    Argument:
    coeffs -- list of coefficients from the denominator of a transfer function from the highest to the 
              lowest power
              
    Return:
    
    routh_matrix -- Sympy matrix
    """
    
    first_line = []
    second_line = []

    for index in range(len(coeffs)):
        if (index%2 == 0) or (index==0):
            first_line.append(coeffs[index])
        else:
            second_line.append(coeffs[index])

    while(len(first_line)!=len(second_line)):
        second_line.append(0)

    routh_matrix = Matrix([first_line, second_line])

    line_n = []

    aux_matrix = [first_line, second_line]

    for line in range(len(coeffs)-2):

        for index in range(len(first_line)-1):

            line_n.append(
                simplify((routh_matrix[line+1,0]*routh_matrix[line,index+1] - routh_matrix[line+1,index+1]*routh_matrix[line,0] )\
                         /routh_matrix[line+1,0]))

        while(len(first_line)!=len(line_n)):

            line_n.append(0)

        aux_matrix.append(line_n)    

        routh_matrix = Matrix(aux_matrix)

        line_n = []

    return routh_matrix
```

Let's take 

$$s^{4} + 59s^{3} + 967s^{2} + (1000k + 3681)s + (50000k + 2772) $$ 

as an example. Note that the coefficients go from highest to lowest order.


```python
k = symbols('k')

coeffs = [1, 59, 967, 1000*k + 3681, 50000*k + 2772]
```


```python
routh_matrix = routh_criterion(coeffs)

routh_matrix
```




$\displaystyle \left[\begin{matrix}1 & 967 & 50000 k + 2772\\59 & 1000 k + 3681 & 0\\\frac{53372}{59} - \frac{1000 k}{59} & 50000 k + 2772 & 0\\\frac{250 \left(1000 k^{2} + 124359 k - 186813\right)}{250 k - 13343} & 0 & 0\\50000 k + 2772 & 0 & 0\end{matrix}\right]$




```python
routh_matrix[4,0]
```




$\displaystyle 50000 k + 2772$




```python
reduce_inequalities(routh_matrix[4,0]>0)
```




$\displaystyle - \frac{693}{12500} < k \wedge k < \infty$




```python
routh_matrix[3,0]
```




$\displaystyle \frac{250 \left(1000 k^{2} + 124359 k - 186813\right)}{250 k - 13343}$




```python
reduce_inequalities(routh_matrix[3,0]>0)
```




$\displaystyle \left(\frac{13343}{250} < k \wedge k < \infty\right) \vee \left(k < - \frac{124359}{2000} + \frac{1239 \sqrt{10561}}{2000} \wedge - \frac{1239 \sqrt{10561}}{2000} - \frac{124359}{2000} < k\right)$




```python

```

23ecb6ff589f6bdb6b76873f82d12396903132c5
ipynb
Jupyter Notebook
routh_criterion_calculator.ipynb
rmcnribeiro/routh-criterion-calculator
c0c1ccb6997aa3213eceed98ff5a9626a5fd35f1
MIT
routh_criterion_calculator.ipynb
rmcnribeiro/routh-criterion-calculator
c0c1ccb6997aa3213eceed98ff5a9626a5fd35f1
MIT
routh_criterion_calculator.ipynb
rmcnribeiro/routh-criterion-calculator
c0c1ccb6997aa3213eceed98ff5a9626a5fd35f1
MIT
Qwen/Qwen-72B
1. YES
2. YES