## MC01TD

### Checking stability of a given real polynomial

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

```  To determine whether or not a given polynomial P(x) with real
coefficients is stable, either in the continuous-time or discrete-
time case.

A polynomial is said to be stable in the continuous-time case
if all its zeros lie in the left half-plane, and stable in the
discrete-time case if all its zeros lie inside the unit circle.

```
Specification
```      SUBROUTINE MC01TD( DICO, DP, P, STABLE, NZ, DWORK, IWARN, INFO )
C     .. Scalar Arguments ..
CHARACTER         DICO
LOGICAL           STABLE
INTEGER           DP, INFO, IWARN, NZ
C     .. Array Arguments ..
DOUBLE PRECISION  DWORK(*), P(*)

```
Arguments

Mode Parameters

```  DICO    CHARACTER*1
Indicates whether the stability test to be applied to
P(x) is in the continuous-time or discrete-time case as
follows:
= 'C':  Continuous-time case;
= 'D':  Discrete-time case.

```
Input/Output Parameters
```  DP      (input/output) INTEGER
On entry, the degree of the polynomial P(x).  DP >= 0.
On exit, if P(DP+1) = 0.0 on entry, then DP contains the
index of the highest power of x for which P(DP+1) <> 0.0.

P       (input) DOUBLE PRECISION array, dimension (DP+1)
This array must contain the coefficients of P(x) in
increasing powers of x.

STABLE  (output) LOGICAL
Contains the value .TRUE. if P(x) is stable and the value

NZ      (output) INTEGER
If INFO = 0, contains the number of unstable zeros - that
is, the number of zeros of P(x) in the right half-plane if
DICO = 'C' or the number of zeros of P(x) outside the unit

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (2*DP+2)
The leading (DP+1) elements of DWORK contain the Routh
coefficients, if DICO = 'C', or the constant terms of
the Schur-Cohn transforms, if DICO = 'D'.

```
Warning Indicator
```  IWARN   INTEGER
= 0:  no warning;
= k:  if the degree of the polynomial P(x) has been
reduced to (DB - k) because P(DB+1-j) = 0.0 on entry
for j = 0, 1,..., k-1 and P(DB+1-k) <> 0.0.

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  if on entry, P(x) is the zero polynomial;
= 2:  if the polynomial P(x) is most probably unstable,
although it may be stable with one or more zeros
very close to either the imaginary axis if
DICO = 'C' or the unit circle if DICO = 'D'.
The number of unstable zeros (NZ) is not determined.

```
Method
```  The stability of the real polynomial
2                DP
P(x) = p(0) + p(1) * x + p(2) * x  + ... + p(DP) x

is determined as follows.

In the continuous-time case (DICO = 'C') the Routh algorithm
(see [1]) is used. The routine computes the Routh coefficients and
if they are non-zero then the number of sign changes in the
sequence of the coefficients is equal to the number of zeros with
positive imaginary part.

In the discrete-time case (DICO = 'D') the Schur-Cohn
algorithm (see [2] and [3]) is applied to the reciprocal
polynomial
2               DP
Q(x) = p(DP) + p(DP-1) * x + p(DP-2) * x  + ... + p(0) x  .

The routine computes the constant terms of the Schur transforms
and if all of them are non-zero then the number of zeros of P(x)
with modulus greater than unity is obtained from the sequence of
constant terms.

```
References
```  [1] Gantmacher, F.R.
Applications of the Theory of Matrices.
Interscience Publishers, New York, 1959.

[2] Kucera, V.
Discrete Linear Control. The Algorithmic Approach.
John Wiley & Sons, Chichester, 1979.

[3] Henrici, P.
Applied and Computational Complex Analysis (Vol. 1).
John Wiley & Sons, New York, 1974.

```
Numerical Aspects
```  The algorithm used by the routine is numerically stable.

Note that if some of the Routh coefficients (DICO = 'C') or
some of the constant terms of the Schur-Cohn transforms (DICO =
'D') are small relative to EPS (the machine precision), then
the number of unstable zeros (and hence the value of STABLE) may
be incorrect.

```
```  None
```
Example

Program Text

```*     MC01TD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          DPMAX
PARAMETER        ( DPMAX = 10 )
*     .. Local Scalars ..
INTEGER          DP, DPP, I, INFO, IWARN, NZ
LOGICAL          STABLE
CHARACTER*1      DICO
*     .. Local Arrays ..
DOUBLE PRECISION DWORK(2*DPMAX+2), P(DPMAX+1)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
EXTERNAL         MC01TD
*     .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
READ ( NIN, FMT = * )
READ ( NIN, FMT = * ) DP, DICO
IF ( DP.LE.-1 .OR. DP.GT.DPMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) DP
ELSE
DPP = DP
READ ( NIN, FMT = * ) ( P(I), I = 1,DP+1 )
*        Determine whether or not the given polynomial P(x) is stable.
CALL MC01TD( DICO, DP, P, STABLE, NZ, DWORK, IWARN, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF ( IWARN.NE.0 ) THEN
WRITE ( NOUT, FMT = 99997 ) IWARN
WRITE ( NOUT, FMT = 99996 ) DPP, DP
END IF
IF ( STABLE ) THEN
WRITE ( NOUT, FMT = 99995 )
ELSE
WRITE ( NOUT, FMT = 99994 )
IF ( LSAME( DICO, 'D' ) ) THEN
WRITE ( NOUT, FMT = 99992 ) NZ
ELSE
WRITE ( NOUT, FMT = 99991 ) NZ
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' MC01TD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MC01TD = ',I2)
99997 FORMAT (' IWARN on exit from MC01TD = ',I2,/)
99996 FORMAT (' The degree of the polynomial P(x) has been reduced fro',
\$       'm ',I2,' to ',I2,/)
99995 FORMAT (' The polynomial P(x) is stable ')
99994 FORMAT (' The polynomial P(x) is unstable ')
99993 FORMAT (/' DP is out of range. ',/' DP = ',I5)
99992 FORMAT (/' The number of zeros of P(x) outside the unit ',
\$       'circle = ',I2)
99991 FORMAT (/' The number of zeros of P(x) in the right ',
\$       'half-plane = ',I2)
END
```
Program Data
``` MC01TD EXAMPLE PROGRAM DATA
4     C
2.0  0.0  1.0  -1.0  1.0
```
Program Results
``` MC01TD EXAMPLE PROGRAM RESULTS

The polynomial P(x) is unstable

The number of zeros of P(x) in the right half-plane =  2
```