**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.

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(*)

**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.

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 .FALSE. otherwise (see also NUMERICAL ASPECTS). 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 circle if DICO = 'D' (see also NUMERICAL ASPECTS).

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'.

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.

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.

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.

[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.

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

**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

MC01TD EXAMPLE PROGRAM DATA 4 C 2.0 0.0 1.0 -1.0 1.0

MC01TD EXAMPLE PROGRAM RESULTS The polynomial P(x) is unstable The number of zeros of P(x) in the right half-plane = 2