Routh Hurwitz

[KennyDang]

Does anyone have a program that can do Routh Hurwitz matrix?

[toshk]

http://www.hpmuseum.org/forum/thread-5528.html

[Han]

Let \( p(z) = a_0 z^n + a_1 z^{n-1} + \dotsm + a_{n-1} z + a_n \) then using HM({ \( a_0, a_1, \dotsm, a_{n-1}, a_n \) }) would produce

\[ \begin{pmatrix}
a_1 & a_3 & a_5 & \dots & \dots & \dots & 0 & 0 & 0 \\
a_0 & a_2 & a_4 & & & & \vdots & \vdots & \vdots \\
0 & a_1 & a_3 & & & & \vdots & \vdots & \vdots \\
\vdots & a_0 & a_2 & \ddots & & & 0 & \vdots & \vdots \\
\vdots & 0 & a_1 & & \ddots & & a_n & \vdots & \vdots \\
\vdots & \vdots & a_0 & & & \ddots & a_{n-1} & 0 & \vdots \\
\vdots & \vdots & 0 & & & & a_{n-2} & a_n & \vdots \\
\vdots & \vdots & \vdots & & & & a_{n-3} & a_{n-1} & 0 \\
0 & 0 & 0 & \dots & \dots & \dots & a_{n-4} & a_{n-2} & a_n
\end{pmatrix} \]

The source code for the HM() program:

Code:

// Compute Hurwitz matrix of p(z)
// a[0]z^n + a[1]z^(n-1) +...+ a[n-1]z + a[n]
// by Han

// input ai := { a[0], a[1], ..., a[n] }
EXPORT HM(ai)
BEGIN
  local n:=size(ai)-1;
  local m:=makemat(0,n,n);
  local j,k,r0,r1,r;

  for j from 1 to n do
    k:=min(j-1,n-j);
    r0:=j-k; r:=j+k+1;
    k:=min(j,n-j);
    r1:=j+k;
    for k from r0 to r1 do
      m[k,j]:=ai[r];
      r:=r-1;
    end; 
  end;
  return(m);
END;

[Han]

And here's a CAS version of the program above, with an extra modification that allows one to compute either the Hurwitz matrix, or the submatrix corresponding to the minors.

Usage: hurwitz({a0,a1,a2,...,an}, s) where s is the size of the submatrix (leading principal minor). If \( s=n \), then the entire matrix is produced. Otherwise, \( s \)-th minor is returned.

Code:

#cas
hurwitz(ai,s):=
BEGIN
  local n, m;
  local j,k,r0,r1,r;

  n:=size(ai)-1;
  if (s<1) or (s>n) then
    return("Bad Argument Value");
  end;
  m:=makemat(0,s,s);

  for j from 1 to s do
    k:=min(j-1,n-j);
    r0:=j-k; r:=j+k+1;
    k:=min(j,n-j);
    r1:=min(s,j+k);
    for k from r0 to r1 do
      m[k,j]:=ai[r];
      r:=r-1;
    end; 
  end;
  return(m);
END;
#end

[Brad Barton]

Thanks for posting both solutions Han. This is very instructive.

Brad