Chapter 11: Stability Analysis of Linear Systems

Lesson 3: Routh–Hurwitz Stability Criterion and Parameter Ranges

This lesson develops the Routh–Hurwitz stability criterion as an algorithmic alternative to explicitly computing poles (roots) of the characteristic polynomial. You will learn how to construct the Routh table, interpret sign changes to count right-half-plane poles, handle special (singular) cases, and derive admissible parameter ranges (e.g., gain intervals) that guarantee closed-loop stability. Computational implementations are provided in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

1. Prerequisites and Problem Setting

From Lesson 2, stability of a continuous-time LTI system can be determined from the pole locations of its transfer function. If the characteristic polynomial is \( p(s) = a_n s^n + a_{n-1} s^{n-1} + \cdots + a_1 s + a_0 \) with real coefficients and \( a_n \neq 0 \), then (for the standard stability region) the system is asymptotically stable if and only if all roots satisfy \( \Re(\lambda_i) < 0 \).

In many control engineering settings, the coefficients \( a_k \) depend on a design parameter (e.g., gain \( K \)). Direct root computation can be numerically sensitive or algebraically cumbersome. The Routh–Hurwitz criterion provides:

  • A root-counting rule for how many roots lie in the right half-plane (RHP), without solving for roots.
  • A direct pathway to parameter inequalities that enforce \( \Re(\lambda_i) < 0 \).

Throughout, assume \( a_n > 0 \) (if not, multiply \( p(s) \) by \( -1 \); this does not change root locations).

2. Routh Table Construction

The Routh table is a structured array built from the coefficients of \( p(s) \). Define rows labeled by descending powers of \( s \): \( s^n, s^{n-1}, \ldots, s^0 \).

The first two rows are formed directly from coefficients by separating even/odd index steps:

\[ \begin{array}{c|cccc} s^n & a_n & a_{n-2} & a_{n-4} & \cdots \\ s^{n-1} & a_{n-1} & a_{n-3} & a_{n-5} & \cdots \end{array} \]

Subsequent rows are computed by a determinant-like recursion. If we denote the entry in row \( i \), column \( j \) as \( r_{i,j} \), with row index increasing downward (\( i=0 \) corresponds to row \( s^n \)), then for \( i \ge 2 \):

\[ r_{i,j} \; ← \; \frac{ r_{i-1,0}\, r_{i-2,j+1} \;-\; r_{i-2,0}\, r_{i-1,j+1} }{r_{i-1,0}}, \quad j = 0,1,2,\ldots \]

Any “missing” coefficient in a row is treated as zero. The computation proceeds until the last row \( s^0 \).

flowchart TD
  A["Input coefficients: a_n ... a_0 (a_n > 0)"] --> B["Build first two rows: s^n and s^(n-1)"]
  B --> C["Compute next row using 2x2 formula (row-by-row)"]
  C --> D["Check special cases: zero leading entry or all-zero row"]
  D --> E["First column signs + sign changes"]
  E --> F["#RHP roots = number of sign changes"]
  F --> G["Stable iff sign changes = 0 and no jw-axis roots"]
        

3. Routh–Hurwitz Criterion

Let \( p(s) \) have real coefficients and no roots exactly on the imaginary axis (\( s = j\omega \)). Construct the Routh table.

Routh–Hurwitz Criterion (Routh form). The number of roots of \( p(s) \) in the open right half-plane \( \{s : \Re(s) > 0\} \) equals the number of sign changes in the first column of the Routh table (counted top to bottom).

Therefore, the polynomial is Hurwitz stable (all roots satisfy \( \Re(\lambda_i) < 0 \)) if and only if:

\[ \text{(i) no first-column sign changes} \quad \text{and} \quad \text{(ii) no imaginary-axis roots}. \]

In typical engineering usage, condition (ii) is addressed by the special-case handling in Section 4 (which detects and resolves borderline cases).

3.1 Proof Sketch (University-Level, Algorithmic Emphasis)

A complete proof is usually presented using the argument principle and the Cauchy index of a rational function built from the even/odd parts of \( p(s) \). We give a structured sketch focusing on the logical steps that justify the sign-change rule.

Define the even and odd parts of \( p(s) \): \( p(s) = p_e(s) + p_o(s) \), where \( p_e(-s)=p_e(s) \) and \( p_o(-s)=-p_o(s) \). Equivalently,

\[ p_e(s) = a_n s^n + a_{n-2} s^{n-2} + a_{n-4} s^{n-4} + \cdots, \quad p_o(s) = a_{n-1} s^{n-1} + a_{n-3} s^{n-3} + \cdots. \]

Consider the rational function \( F(s) = \dfrac{p_e(s)}{p_o(s)} \). Along the imaginary axis \( s=j\omega \), and under the assumption that \( p_o(j\omega)\neq 0 \) for all \( \omega \), the quantity \( F(j\omega) \) is well-defined.

Key idea. The change in argument of \( p(j\omega) \) as \( \omega \) runs from \( -\infty \) to \( +\infty \) is linked (by the argument principle) to the number of zeros in the RHP. This argument change can be related to the number of “jumps” in the real rational function \( F(j\omega) \) (its Cauchy index).

The Routh algorithm corresponds to computing a continued-fraction (Euclidean-algorithm) expansion of \( F(s) \). The successive leading coefficients in the first column of the Routh table determine the signs of the quotients in that continued fraction. Each sign change in the first column corresponds to a change in the Cauchy index, which corresponds to one root crossing into the RHP.

Concretely, one can show (under the “no imaginary-axis roots” assumption) that:

\[ N_{\text{RHP}} = V, \quad \text{where } V \text{ is the number of sign changes in the first column of the Routh table.} \]

This completes the proof sketch: the Routh table encodes the same root-counting information as the argument principle, without explicit root computation.

4. Special Cases and Their Remedies

Real-world characteristic polynomials can produce singularities in the Routh table. These cases must be handled correctly to avoid false stability conclusions.

4.1 Case A: First Column Entry Becomes Zero

If the first element of some row is exactly zero, the recursion divides by zero. A standard remedy is to replace the zero with a small positive number \( \varepsilon \) and proceed, then interpret sign changes in the limit \( \varepsilon → 0^+ \).

\[ r_{i,0} = 0 \;\Rightarrow\; r_{i,0} \; ← \; \varepsilon, \quad \varepsilon > 0,\quad \text{then analyze sign changes as } \varepsilon → 0^+. \]

Interpretation: a zero leading entry typically indicates a root on or symmetric with respect to the imaginary axis boundary (often signaling marginal stability conditions).

4.2 Case B: An Entire Row Becomes Zero

If a full row becomes zero, the polynomial has roots symmetrically located about the origin, often involving imaginary-axis roots. In this case, form the auxiliary polynomial from the row immediately above the zero row, using its coefficients as a polynomial in \( s \).

Suppose the zero row corresponds to power \( s^k \). Let the row above it (power \( s^{k+1} \)) contain coefficients \( b_0, b_1, \ldots \). Define:

\[ A(s) = b_0 s^{k+1} + b_1 s^{k-1} + b_2 s^{k-3} + \cdots \]

Replace the zero row with the coefficients of \( \dfrac{dA(s)}{ds} \) arranged in the same Routh format, then continue.

\[ \text{Zero row} \;\Rightarrow\; \text{replace by coefficients of } \frac{dA(s)}{ds}. \]

This procedure both detects and correctly counts imaginary-axis roots (and repeated roots on the axis).

5. Parameter Ranges via Routh Inequalities

Suppose \( p(s;K) \) depends on a real parameter \( K \) (gain, stiffness, damping, etc.). The Routh criterion turns stability into a set of inequalities: all first-column entries must be strictly positive (given \( a_n > 0 \) and no imaginary-axis roots).

If the first column is \( c_0(K), c_1(K), \ldots, c_n(K) \), stability requires: \( c_i(K) > 0 \) for all \( i \).

\[ \forall i \in \{0,1,\ldots,n\}:\quad c_i(K) > 0. \]

flowchart TD
  P["Given p(s;K)"] --> R["Build Routh table symbolically or numerically"]
  R --> C["Extract first-column entries c_i(K)"]
  C --> I["Solve inequalities: c_i(K) > 0 for all i"]
  I --> S["Intersect all solution sets -> stable K-range"]
  S --> X["Optional: verify by root check (numerical) on a test K"]
        

5.1 Worked Example (Cubic, Single Parameter)

Consider: \( p(s;K) = s^3 + 2 s^2 + (1+K)s + K \). Construct the Routh table:

\[ \begin{array}{c|cc} s^3 & 1 & 1+K \\ s^2 & 2 & K \\ s^1 & \dfrac{2(1+K)-1\cdot K}{2} & 0 \\ s^0 & K & \\ \end{array} \]

Simplify the \( s^1 \) first-column entry:

\[ \frac{2(1+K)-K}{2} = \frac{2+2K-K}{2} = \frac{2+K}{2}. \]

Stability requires positivity of the first column:

\[ 1 > 0,\quad 2 > 0,\quad \frac{2+K}{2} > 0,\quad K > 0 \;\Rightarrow\; K > 0. \]

Hence the system is asymptotically stable for \( K > 0 \) (and marginal at \( K=0 \)).

5.2 Practical Note on “Borderline” Parameters

Parameter values that make any first-column entry exactly zero correspond to imaginary-axis roots or repeated roots. In control engineering, such values are typically excluded from strict asymptotic stability: \( \Re(\lambda_i) < 0 \) requires strict inequalities in the Routh conditions.

6. Python Implementation

Recommended libraries: \( \texttt{control} \) (Control Systems Library), \( \texttt{scipy.signal} \) (transfer functions), \( \texttt{sympy} \) (symbolic inequalities), \( \texttt{numpy} \).


import numpy as np

def routh_table(coeffs, eps=1e-12):
    """
    Build the Routh table for a real-coefficient polynomial.

    coeffs: [a_n, a_{n-1}, ..., a_0] with a_n != 0
    returns: 2D numpy array representing the Routh table
    """
    coeffs = np.array(coeffs, dtype=float)
    n = len(coeffs) - 1  # degree
    rows = n + 1
    cols = int(np.ceil((n + 1) / 2.0))

    R = np.zeros((rows, cols), dtype=float)

    # first two rows
    R[0, :len(coeffs[0::2])] = coeffs[0::2]   # a_n, a_{n-2}, ...
    R[1, :len(coeffs[1::2])] = coeffs[1::2]   # a_{n-1}, a_{n-3}, ...

    for i in range(2, rows):
        # If leading element is zero, use epsilon method
        if abs(R[i-1, 0]) < eps:
            R[i-1, 0] = eps

        for j in range(0, cols - 1):
            a = R[i-2, 0]
            b = R[i-2, j+1]
            c = R[i-1, 0]
            d = R[i-1, j+1]
            R[i, j] = (c*b - a*d) / c

        # Special case: entire row becomes (near) zero
        if np.all(np.abs(R[i, :]) < eps):
            # Auxiliary polynomial from row i-1
            k = n - (i - 1)  # power associated with row i-1 is s^k
            # row i-1 corresponds to A(s) = sum row[i-1, m] * s^(k-2m)
            # Replace row i with derivative coefficients:
            for m in range(cols):
                power = k - 2*m
                if power < 0:
                    break
                # derivative: power * coeff * s^(power-1)
                R[i, m] = power * R[i-1, m]

    return R

def routh_sign_changes(R, eps=1e-12):
    first_col = R[:, 0].copy()
    # Treat near-zero as small positive to avoid spurious flips in counting
    first_col[np.abs(first_col) < eps] = eps
    signs = np.sign(first_col)
    # Count sign changes ignoring zeros
    changes = 0
    for k in range(len(signs)-1):
        if signs[k] * signs[k+1] < 0:
            changes += 1
    return changes

# Example: p(s;K) = s^3 + 2 s^2 + (1+K)s + K, test K=0.5
K = 0.5
coeffs = [1, 2, 1+K, K]
R = routh_table(coeffs)
print("Routh table:\n", R)
print("Sign changes in first column:", routh_sign_changes(R))

# Optional numerical verification using control (root computation)
# from control import tf
# den = [1, 2, 1+K, K]
# sys = tf([1], den)
# print("Poles:", sys.poles())
      

For symbolic parameter ranges, use \( \texttt{sympy} \) to construct the first-column expressions and solve inequalities. A common workflow is: (i) build the Routh table algebraically (symbolic), (ii) extract first column, (iii) solve all strict inequalities and intersect.


import sympy as sp

def routh_table_symbolic(coeffs):
    # coeffs: sympy list [a_n, ..., a_0]
    n = len(coeffs) - 1
    rows = n + 1
    cols = int(sp.ceiling((n + 1) / 2))

    R = sp.Matrix.zeros(rows, cols)
    R[0, :len(coeffs[0::2])] = sp.Matrix([coeffs[0::2]])
    R[1, :len(coeffs[1::2])] = sp.Matrix([coeffs[1::2]])

    for i in range(2, rows):
        for j in range(0, cols - 1):
            a = R[i-2, 0]
            b = R[i-2, j+1]
            c = R[i-1, 0]
            d = R[i-1, j+1]
            # keep as rational expression (assume c != 0 for now)
            R[i, j] = sp.simplify((c*b - a*d) / c)
    return R

K = sp.Symbol('K', real=True)
coeffs_sym = [1, 2, 1+K, K]  # s^3 + 2 s^2 + (1+K)s + K
R_sym = routh_table_symbolic(coeffs_sym)
first_col = [sp.simplify(R_sym[i, 0]) for i in range(R_sym.rows)]
print("First column:", first_col)

ineqs = [sp.StrictGreaterThan(expr, 0) for expr in first_col]
sol = sp.reduce_inequalities(ineqs, K)
print("Stable K-range:", sol)
      

7. C++ Implementation (From Scratch)

There is no universal “standard” system dynamics library in C++ comparable to MATLAB’s Control System Toolbox, so in control engineering practice you either: (i) implement the criterion directly (as below), (ii) use numeric root finding for verification (e.g., Eigen + companion matrix), or (iii) integrate with dedicated control toolkits when available.


#include <iostream>
#include <vector>
#include <cmath>

static int signum(double x, double eps=1e-12) {
    if (std::fabs(x) < eps) return 1; // treat near-zero as positive epsilon
    return (x > 0.0) ? 1 : -1;
}

std::vector<std::vector<double>> routhTable(const std::vector<double>& coeffs, double eps=1e-12) {
    int n = static_cast<int>(coeffs.size()) - 1; // degree
    int rows = n + 1;
    int cols = (n + 2) / 2; // ceil((n+1)/2)

    std::vector<std::vector<double>> R(rows, std::vector<double>(cols, 0.0));

    // first two rows
    int c0 = 0, c1 = 0;
    for (int i = 0; i <= n; i += 2) R[0][c0++] = coeffs[i];
    for (int i = 1; i <= n; i += 2) R[1][c1++] = coeffs[i];

    for (int i = 2; i < rows; ++i) {
        if (std::fabs(R[i-1][0]) < eps) R[i-1][0] = eps;

        for (int j = 0; j < cols - 1; ++j) {
            double a = R[i-2][0];
            double b = R[i-2][j+1];
            double c = R[i-1][0];
            double d = R[i-1][j+1];
            R[i][j] = (c*b - a*d) / c;
        }

        // Row-of-zeros handling (basic numeric version)
        bool allZero = true;
        for (int j = 0; j < cols; ++j) {
            if (std::fabs(R[i][j]) > eps) { allZero = false; break; }
        }
        if (allZero) {
            // Auxiliary polynomial from row i-1
            int k = n - (i - 1); // power for row i-1
            for (int m = 0; m < cols; ++m) {
                int power = k - 2*m;
                if (power < 0) break;
                R[i][m] = power * R[i-1][m];
            }
        }
    }
    return R;
}

int routhSignChanges(const std::vector<std::vector<double>>& R, double eps=1e-12) {
    int changes = 0;
    for (size_t i = 0; i + 1 < R.size(); ++i) {
        int s1 = signum(R[i][0], eps);
        int s2 = signum(R[i+1][0], eps);
        if (s1 * s2 < 0) changes++;
    }
    return changes;
}

int main() {
    // Example: p(s) = s^3 + 2 s^2 + 1.5 s + 0.5  (K=0.5 in the lesson example)
    std::vector<double> coeffs = {1.0, 2.0, 1.5, 0.5};

    auto R = routhTable(coeffs);
    std::cout << "Routh table:\n";
    for (const auto& row : R) {
        for (double v : row) std::cout << v << "\t";
        std::cout << "\n";
    }

    int V = routhSignChanges(R);
    std::cout << "Sign changes in first column: " << V << "\n";
    std::cout << "RHP roots (count): " << V << "\n";
    return 0;
}
      

8. Java Implementation

Java workflows often use \( \texttt{Apache Commons Math} \) for polynomial utilities, but the Routh table is straightforward to implement directly. The code below builds the table and counts sign changes.


import java.util.Arrays;

public class RouthHurwitz {

    static int signum(double x, double eps) {
        if (Math.abs(x) < eps) return 1;  // treat near-zero as +eps
        return (x > 0.0) ? 1 : -1;
    }

    static double[][] routhTable(double[] coeffs, double eps) {
        int n = coeffs.length - 1; // degree
        int rows = n + 1;
        int cols = (n + 2) / 2; // ceil((n+1)/2)

        double[][] R = new double[rows][cols];

        // first two rows
        int c0 = 0, c1 = 0;
        for (int i = 0; i <= n; i += 2) R[0][c0++] = coeffs[i];
        for (int i = 1; i <= n; i += 2) R[1][c1++] = coeffs[i];

        for (int i = 2; i < rows; i++) {
            if (Math.abs(R[i-1][0]) < eps) R[i-1][0] = eps;

            for (int j = 0; j < cols - 1; j++) {
                double a = R[i-2][0];
                double b = R[i-2][j+1];
                double c = R[i-1][0];
                double d = R[i-1][j+1];
                R[i][j] = (c*b - a*d) / c;
            }

            // row of zeros handling (basic numeric)
            boolean allZero = true;
            for (int j = 0; j < cols; j++) {
                if (Math.abs(R[i][j]) > eps) { allZero = false; break; }
            }
            if (allZero) {
                int k = n - (i - 1); // power of row i-1
                for (int m = 0; m < cols; m++) {
                    int power = k - 2*m;
                    if (power < 0) break;
                    R[i][m] = power * R[i-1][m];
                }
            }
        }

        return R;
    }

    static int signChanges(double[][] R, double eps) {
        int changes = 0;
        for (int i = 0; i + 1 < R.length; i++) {
            int s1 = signum(R[i][0], eps);
            int s2 = signum(R[i+1][0], eps);
            if (s1 * s2 < 0) changes++;
        }
        return changes;
    }

    public static void main(String[] args) {
        // Example polynomial: s^3 + 2 s^2 + 1.5 s + 0.5
        double[] coeffs = new double[]{1.0, 2.0, 1.5, 0.5};
        double eps = 1e-12;

        double[][] R = routhTable(coeffs, eps);
        System.out.println("Routh table:");
        for (double[] row : R) {
            System.out.println(Arrays.toString(row));
        }

        int V = signChanges(R, eps);
        System.out.println("Sign changes (first column): " + V);
        System.out.println("RHP roots (count): " + V);
    }
}
      

Optional verification: compute polynomial roots numerically using a companion-matrix eigenvalue method (requires a linear algebra package). In this course sequence, Routh is preferred because it avoids root solving and directly yields parameter inequalities.

9. MATLAB and Simulink

MATLAB’s Control System Toolbox supports pole computation (e.g., \( \texttt{pole} \)) and symbolic tools can solve inequalities. However, implementing the Routh table is instructive and enables direct parameter-range derivations.

9.1 MATLAB: Routh Table Function


function R = routh_table(coeffs)
%ROUTH_TABLE Build Routh table for polynomial with coefficients [a_n ... a_0]
    coeffs = coeffs(:).';
    n = length(coeffs) - 1;
    rows = n + 1;
    cols = ceil((n+1)/2);
    R = zeros(rows, cols);

    R(1,1:length(coeffs(1:2:end))) = coeffs(1:2:end);
    R(2,1:length(coeffs(2:2:end))) = coeffs(2:2:end);

    epsv = 1e-12;
    for i = 3:rows
        if abs(R(i-1,1)) < epsv
            R(i-1,1) = epsv;
        end
        for j = 1:cols-1
            a = R(i-2,1);
            b = R(i-2,j+1);
            c = R(i-1,1);
            d = R(i-1,j+1);
            R(i,j) = (c*b - a*d)/c;
        end

        % Row-of-zeros handling (basic numeric)
        if all(abs(R(i,:)) < epsv)
            k = n - (i-2); % power associated with row i-1
            for m = 1:cols
                power = k - 2*(m-1);
                if power < 0
                    break;
                end
                R(i,m) = power * R(i-1,m);
            end
        end
    end
end
      

9.2 MATLAB: Parameter Range Example with Symbolics


syms K real
% p(s;K) = s^3 + 2 s^2 + (1+K)s + K
% First-column conditions derived from the Routh table:
% 1 > 0, 2 > 0, (2+K)/2 > 0, K > 0
ineqs = [K > 0, (2+K)/2 > 0];
solK = solve(ineqs, K, 'ReturnConditions', true);
disp(solK.conditions)
      

9.3 Simulink: How Routh Fits into a Modeling Workflow

In Simulink, you typically build a block diagram model (plant + controller) and linearize it around an operating point. A standard workflow is:

  1. Build the closed-loop model using \( \texttt{Transfer Fcn} \) and \( \texttt{Sum} \) blocks (or state-space blocks later).
  2. Use MATLAB to extract the linear transfer function (or denominator polynomial) and apply the Routh test.

% Example outline (requires a Simulink model named 'myClosedLoopModel'):
% 1) Linearize model (tooling varies by MATLAB version/toolboxes)
% sys = linearize('myClosedLoopModel');

% 2) Convert to transfer function and extract denominator:
% [num, den] = tfdata(tf(sys), 'v');

% 3) Apply Routh:
% R = routh_table(den);
% firstCol = R(:,1);
% disp(firstCol);
      

Even when you verify stability by poles, the Routh inequalities remain valuable for tuning parameters without repeated root finding.

10. Wolfram Mathematica

Mathematica is well-suited for symbolic inequality solving (parameter ranges). Below is a compact Routh-table construction followed by a parameter range solve.


(* Build Routh table for polynomial coefficients {a_n, a_{n-1}, ..., a_0} *)
RouthTable[coeffs_List, eps_: 10^-12] := Module[
  {n, rows, cols, R, even, odd, i, j, a, b, c, d, allZero, k, power},
  n = Length[coeffs] - 1;
  rows = n + 1;
  cols = Ceiling[(n + 1)/2];
  R = ConstantArray[0, {rows, cols}];

  even = coeffs[[1 ;; ;; 2]];
  odd  = coeffs[[2 ;; ;; 2]];
  Do[R[[1, i]] = even[[i]], {i, 1, Length[even]}];
  Do[R[[2, i]] = odd[[i]],  {i, 1, Length[odd]}];

  For[i = 3, i <= rows, i++,
    If[Abs[R[[i - 1, 1]]] < eps, R[[i - 1, 1]] = eps];
    For[j = 1, j <= cols - 1, j++,
      a = R[[i - 2, 1]]; b = R[[i - 2, j + 1]];
      c = R[[i - 1, 1]]; d = R[[i - 1, j + 1]];
      R[[i, j]] = Simplify[(c*b - a*d)/c];
    ];

    allZero = And @@ (Abs[R[[i, #]]] < eps & /@ Range[cols]);
    If[allZero,
      k = n - (i - 2); (* power associated with row i-1 *)
      For[j = 1, j <= cols, j++,
        power = k - 2*(j - 1);
        If[power < 0, Break[]];
        R[[i, j]] = power * R[[i - 1, j]];
      ];
    ];
  ];
  R
];

(* Example: p(s;K) = s^3 + 2 s^2 + (1+K)s + K *)
Clear[K];
R = RouthTable[{1, 2, 1 + K, K}];
firstCol = R[[All, 1]] // Simplify
(* Solve strict inequalities for stability *)
Reduce[Thread[firstCol > 0] && K ∈ Reals, K]
      

11. Problems and Solutions

The problems below reinforce (i) building the Routh table, (ii) counting RHP roots, and (iii) extracting parameter ranges. All problems assume real coefficients and continuous-time stability region.


Problem 1 (Routh table for a cubic): For \( p(s)=s^3+4s^2+6s+3 \), build the Routh table and determine whether the polynomial is Hurwitz stable.

Solution:

Construct rows:

\[ \begin{array}{c|cc} s^3 & 1 & 6 \\ s^2 & 4 & 3 \\ s^1 & \dfrac{4\cdot 6 - 1\cdot 3}{4} & 0 \\ s^0 & 3 & \end{array} \]

Compute:

\[ \frac{24-3}{4} = \frac{21}{4} > 0. \]

First column is \( [1,4,21/4,3] \), all positive, hence no sign changes and (no singular case), so \( p(s) \) is Hurwitz stable.


Problem 2 (Counting RHP roots): Consider \( p(s)=s^4 + 2s^3 - s^2 - 2s + 1 \). Use the Routh criterion to find the number of roots with \( \Re(s) > 0 \).

Solution:

First two rows:

\[ \begin{array}{c|ccc} s^4 & 1 & -1 & 1 \\ s^3 & 2 & -2 & 0 \end{array} \]

Next row (\( s^2 \)):

\[ r_{2,0} = \frac{2(-1) - 1(-2)}{2} = \frac{-2+2}{2}=0, \quad r_{2,1} = \frac{2(1) - 1(0)}{2}=1. \]

The leading entry is zero; use the epsilon method: set \( r_{2,0} ← \varepsilon \), \( \varepsilon > 0 \). Then compute row \( s^1 \) first element:

\[ r_{3,0} = \frac{\varepsilon(-2) - 2(1)}{\varepsilon} = -2 - \frac{2}{\varepsilon} < 0. \]

First-column signs (top to bottom): \( +,+,+,-,+ \) (since last row is \( s^0=1 \)). Sign changes: from \( + \) to \( - \) (one), then \( - \) to \( + \) (second).

Therefore \( N_{\text{RHP}}=2 \).


Problem 3 (Parameter range, cubic): Find all \( K \in \mathbb{R} \) such that \( p(s;K) = s^3 + 5s^2 + (6+K)s + (2K+1) \) is Hurwitz stable.

Solution:

Routh table for a cubic:

\[ \begin{array}{c|cc} s^3 & 1 & 6+K \\ s^2 & 5 & 2K+1 \\ s^1 & \dfrac{5(6+K)-1(2K+1)}{5} & 0 \\ s^0 & 2K+1 & \end{array} \]

Simplify the \( s^1 \) first-column entry:

\[ \frac{30+5K-2K-1}{5} = \frac{29+3K}{5}. \]

Stability requires:

\[ 1 > 0,\quad 5 > 0,\quad \frac{29+3K}{5} > 0,\quad 2K+1 > 0. \]

Thus:

\[ 29+3K > 0 \Rightarrow K > -\frac{29}{3}, \quad 2K+1 > 0 \Rightarrow K > -\frac{1}{2}. \]

Intersecting yields \( K > -\dfrac{1}{2} \).


Problem 4 (Row of zeros and auxiliary polynomial): Consider \( p(s)=s^4 + 2s^2 + 1 \). Explain why a row of zeros occurs and determine the location of roots.

Solution:

First rows:

\[ \begin{array}{c|ccc} s^4 & 1 & 2 & 1 \\ s^3 & 0 & 0 & 0 \end{array} \]

The \( s^3 \) row is all zeros immediately, indicating symmetry and imaginary-axis structure. Build the auxiliary polynomial from the \( s^4 \) row:

\[ A(s) = s^4 + 2s^2 + 1 = (s^2+1)^2. \]

Hence roots are \( s=\pm j \), each with multiplicity 2. This is not asymptotically stable because roots lie on the imaginary axis.


Problem 5 (Design inequality interpretation): For a stable design you derived inequalities \( c_i(K) > 0 \). Explain why strict inequalities are required for asymptotic stability rather than \( c_i(K) \ge 0 \).

Solution:

If any first-column entry becomes zero, the Routh construction becomes singular, which corresponds to at least one root on the imaginary axis (or a root at the origin) under the Routh framework. Such roots imply the system is not asymptotically stable because asymptotic stability requires \( \Re(\lambda_i) < 0 \) strictly. Therefore stability regions derived from Routh are expressed as strict inequalities, and boundary parameter values must be analyzed separately (typically yielding marginal stability).

12. References

  1. Hurwitz, A. (1895). Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Teilen besitzt. Mathematische Annalen, 46, 273–284.
  2. Liénard, A., & Chipart, H. (1914). Sur le signe de la partie réelle des racines d’une équation algébrique. Bulletin de la Société Mathématique de France, 42, 1–24.
  3. Bistritz, Y. (2004). Zero location with respect to the imaginary axis of polynomials with real coefficients. Proceedings of the IEEE, 92(5), 825–853.
  4. Jury, E.I. (1964). On the stability of a general class of digital filters. IRE Transactions on Circuit Theory, 11(2), 185–195.
  5. Kharitonov, V.L. (1978). Asymptotic stability of an equilibrium position of a family of systems of linear differential equations. Differential Equations, 14, 1483–1485.
  6. Anderson, B.D.O. (1972). Stability of linear systems and the algebraic criterion of Hurwitz. International Journal of Control, 16(5), 1033–1043.
  7. Barmish, B.R. (1984). Invariance of the strict Hurwitz property under coefficient perturbations. IEEE Transactions on Automatic Control, 29(9), 845–846.
  8. Gantmacher, F.R. (1937). On the theory of matrices and stability criteria (early developments). Comptes Rendus (Doklady) de l’Académie des Sciences de l’URSS, various notes (1930s).

13. Summary

The Routh–Hurwitz criterion converts continuous-time stability of a characteristic polynomial into a table-based computation: the number of right-half-plane roots equals the number of sign changes in the first column. This avoids explicit pole solving and enables direct derivation of parameter ranges by enforcing strict positivity of first-column expressions. You also learned how to handle singular cases (zero leading element and row of zeros) that indicate marginal conditions or imaginary-axis roots.