Stability of polynomials

Problems from the theory of ordinary differential equations, numerical mathematics and control theory lead to the investigation of stability of polynomials. The authors review some criteria and show possibilities  offered by Mathematica.

by János Tóth, László Szili and András Zachár

The polynomial p is said to be stable if all the roots of p lie in the open left half of the complex plane. In the complementary case the polynomial is said to be unstable. The problem emerged at the end of the last century in connection with the centrifugal regulator, theoretically treated by such scholars as Maxwell, Hermite and Liapunov. A practically useful approach to that special problem has been given by Vishnegradsky [7].

1 INTRODUCTION

We shall only consider polynomials with real coefficients, i.e. polynomials of the form

\(TraditionalForm
      \`p(z) := \(a\_0\) z\^N + \(a\_1\) z\^\(N - 1\) + ⋯ + 
          \(a\_\(N - 1\)\) z + 
          a\_N\ \ \ \ \((z\ \ ϵ\ \ ℂ)\)\)

where N is a fixed positive integer and (Subscript[a, n]) (n = 1, 2, ⋯, N) are real numbers.

An obvious approach to the problem is to calculate the roots of the polynomial and check if the stability condition holds or not. This is by far the least effective method—especially in the case of symbolical coefficients. Therefore the genuine problem is: How to check the stability condition only using the coefficients but without finding the roots to the polynomial?

It turns out that it is not a simple problem to find such procedures. Many people were involved in this area, the most well known ones being Routh, Stodola, Hurwitz, Liénard, Chippart and Mikhajlov. The theory of stable polynomials (together with historical remarks) has been outlined in [2,3,4].

Our aim in the present paper is twofold: First, we are going to present a short review of this area (as it seems to be much less known among both algebraists and analysts as it should be). Second, attached to each criteria we provide some Mathematica code to illustrate and use the given method.

Some of the given criteria can utilize the symbolic capabilities of  Mathematica; in those cases we formulated the code so as to make this possible.

No checking (as e.g. type check) is included into the codes in order to make them as transparent and simple as possible. Although we summarize our functions at the end into a single function with options, we think that further studies and experiments are needed until it can be considered as a 'product'.

Finally, we remark that there has been an earlier approach to the stability of polynomials using Mathematica by A. Urintsev, see in MathSource  under item number 0203–746.

2 STODOLA'S CRITERION

It is easy to find a necessary condition for a polynomial to be stable. According to the Stodola  condition [2, p. 488] the polynomial p can be stable only if all its coefficients are of the same sign.

The function below can be used to check this criterion:

\(Stodola[polyvalue_, variable_] 
      := 
\ \ \ \ \ \ \ Equal@@
        Sign[CoefficientList[polyvalue, variable]]\)

Let us remark that our aim is to check whether the sign of the coefficients of a polynomial are equal or not. The right hand side follows the syntax of the previous English sentence almost exactly.

Let us consider the following polynomials:

(Subscript[p, 1])[z_] := (z^3) + (z^2) + 4  z + 30(Subscript[p, 2])[z_] := 0  (z^3) + (z^2) + 5  z + 3(Subscript[p, 3])[z_, a_] := 3  (z^3) + (z^2) + 5  z + a

Let us check how the function Stodola[]works.

Stodola[(Subscript[p, 1])[z], z]

Our function considers the first nonzero coefficient to be the leading coefficient.

Stodola[(Subscript[p, 2])[z], z]

As it can be seen from the examples above and from some examples below, cases given by approximate numerical values are evaluated to give either False or True.  If the coefficients have exact numerical values (e.g. 1/2, (Sqrt[2]), (π/e), Log[2]) then previous versions (such as Version 2.2) may not yet give a False or True result.

polylist1 = 
        (Subscript[p, 3])[z, a] /. 
          a -> {-1, 0, (Sqrt[2]), Log[2], Log[2. ], Log[2, (1/2)]};TableForm[polylist1]

Map[Stodola[#, z] &, polylist1] // TableForm

If one also has parametric coefficients then a result of the type below will be given.

Stodola[(Subscript[p, 3])[z, a], z]

The criterion of Stodola is only necessary and not sufficient for the stability of a polynomial, as it is shown by the polynomial (Subscript[p, 1]).

Solve[(Subscript[p, 1])[z] == 0, z]

3 PLOTTING THE GRAPH OF A POLYNOMIAL AND SOLVING A POLYNOMIAL EQUATION

If one has an easy to use and efficient tool for drawing functions on a computer like Mathematica, one can try to represent the graph of the function. We can see below two examples.

(Subscript[p, 4])[z_] := 
      (z^5) + 31  (z^4) + 361  (z^3) + 1921  (z^2) + 4390  z + 2801Plot[(Subscript[p, 4])[x], {x, -10.1, -0.9}];

(Subscript[p, 5])[z_] := 4  (z^5) + 4  (z^4) + 8  (z^3) + 5  (z^2) + 3  z + 1Plot[(Subscript[p, 5])[x], {x, -2, 2}];

In the first case the polynomial has five real roots and it is a very lucky example because one can see every root and one can state easily that the polynomial (Subscript[p, 4]) is stable. The second polynomial has one real and four complex roots.

NSolve[(Subscript[p, 5])[z] == 0,  z]

But there are some problems, first, one may not get the exact values of the roots graphically, and, second, the Plot[] function does not show the complex roots (as for example in the second case).

Let us see what help can be obtained if we want to calculate the roots. The function Solve[] can be used to symbolically find the roots of many polynomial equations (even with symbolic coefficients), still its area of validity is rather limited. One can however solve the equation p(x) = 0 numerically with the NSolve[] function which gives a list of numerical approximations to the roots of a polynomial equation.

Here is a one-liner solving our problem based on NSolve[].

\(\(NegativeRealPartsQ[polyvalue_, variable_] 
      := 
	\ \ And@@
        Negative[Re[variable /. NSolve[polyvalue == 0, variable]]]
\)\)

Let us see some examples.

(Subscript[p, 6])[z_] := 3  (z^3) + 2  (z^2) + 5  z - 1(Subscript[p, 7])[z_] := (z^4) + 2  (z^3) + (z^2) + z + 1(Subscript[p, 8])[z_] := Evaluate[(∑_(k = 0)^4) ((z^k)/(k!))]\(p\_9[z_] := Evaluate[\(\(\ ∑\)\+\(k = 0\)\%5\) z\^k\/\(k!\)]\)NegativeRealPartsQ[#[z], z] & /@ 
      {(Subscript[p, 1]), (Subscript[p, 4]), (Subscript[p, 5]), (Subscript[p, 6]), (Subscript[p, 7]), (Subscript[p, 8]), (Subscript[p, 9])}

The above function should be applied carefully because if the degree of the polynomial is high the numerical method used by  NSolve[] may not provide a reliable result as a consequence of rounding errors. As an example we ask the reader to compare the known roots of

(Subscript[p, 100])[z_] := 
      Evaluate[Expand[(∏_(i = 1)^100)(z + N[(Sqrt[i])])]]

with the results obtained by

NSolve[(Subscript[p, 100])[z] == 0, z];

Increasing precision may lead to reliable results. (No wonder that neither Maple V, nor Derive nor MATLAB performs better at this point: the effect highly depends on the chosen machine precision of the given mathematical program package.) According to our experience one and a half time the degree of the polynomial as precision suffices. For example, the roots of the (approximately defined!) polynomial

(Subscript[q, 100])[z_] := 
      Evaluate[Expand[(∏_(i = 1)^100)(z + N[(Sqrt[i]), 150])]]

can be obtained by

NSolve[(Subscript[q, 100])[z] == 0, z, 150];

This rescue manoeuvre however makes the calculations essentially slower.

In the case of coefficients containing parameters the use of  NSolve[] is rather limited.

4 THE HERMITE–MIKHAJLOV CRITERION

This criterion of geometric character has been found by Hermite in 1851 and rediscovered and popularized among engineers by Mikhajlov in 1937.

Let us start from the fact that polynomial (1) is an analytical function on the complex plane. As a consequence of Cauchy's theorem it is completely determined by its values taken on on a single line of the complex plane. Let this line be the imaginary axis. Then the range (Subscript[Γ, p]) of the map

\(TraditionalForm
      \`\(ℝ ∋ 
        t ↦ p(i\ t) ∈ 
        ℂ\ \ \)\)

is the image of the imaginary axis. This curve is usually called as the amplitude-phase characteristic or Mikhajlov's hodograph of the polynomial p. Certain geometric properties of (Subscript[Γ, p]) are strongly conneceted with the stability of the underlying polynomial.

The following statement is said to be the amplitude-phase citerion of stability (see [4, p. 37]).

Theorem 1

The polynomial p is stable if and only if the curve (Subscript[Γ, p]) does not cross the origin and turns around the origine in positive direction (or: anticlockwise) at an angle of n π while its parameter t changes from -∞ to +∞.

Let us draw the amplitude-phase characteristic of the polynomial (Subscript[p, 5]).

ParametricPlot[{Re[(Subscript[p, 5])[I t]], Im[(Subscript[p, 5])[I t]]}, {t, -5, 5}, 
      PlotRange -> {-1.5, 1.5}, 
      Ticks -> {{1.5, 3}, {-1.5, 1.5}}, 
      Prolog -> 
        {Text["t=1.29", {4.5, -1}], PointSize[0.015], 
          Point[{Re[(Subscript[p, 5])[1.29  I]], Im[(Subscript[p, 5])[1.29  I]]}], 
          Point[{Re[(Subscript[p, 5])[-1.29 I]], Im[(Subscript[p, 5])[-1.29 I]]}]}];

The segment of the curve suggests that the criterion is fulfilled what really has been demonstrated in Section 3. There we also saw that the polynomial (Subscript[p, 9]) was unstable. The amplitude-phase diagram of that polynomial can be plotted in a similar way.

ParametricPlot[{Re[(Subscript[p, 9])[I t]], Im[(Subscript[p, 9])[I t]]}, {t, -3.7, 3.7}, 
      PlotRange -> {-1.5, 1.5}, 
      Ticks -> {{-0.5, 0.5}, {-1.5, 1.5}}, 
      Prolog -> 
        {Text["t=3.6", {1.5, 1.2}], Text["t=-3.6", {1.5, -1.2}], 
          PointSize[0.015], Point[{Re[(Subscript[p, 9])[3.6  I]], Im[(Subscript[p, 9])[3.6  I]]}], 
          Point[{Re[(Subscript[p, 9])[-3.6 I]], Im[(Subscript[p, 9])[-3.6 I]]}]}];

No one-liner is supplied for this method because the appropriate parameter domains can only be found by trial and error.

Theorem 1 can also be translated into the language of algebra. It is easy to see that the curve (Subscript[Γ, p]) is symmetrical to the real axis. A corllary of this fact is that one can work with polynomials of half the original degree.

Let us suppose that N > 2 and let us introduce two polynomials for positive arguments in the following way:

G ((t^2)) := Re (p (i t)),     H ((t^2)):=Im (p (i t))/
              t     (t ϵ ℝ)

The next statement is called Hermite's criterion or Mikhajlov's criterion (see e.g. [4, p. 48]).

Theorem 2

Let N > 2. Then the polynomial p in (1) is stable if and only if (Subscript[a, N]) (Subscript[a, N - 1]) > 0 and the roots (Subscript[g, 1]), (Subscript[g, 2, ]) ⋯,  and (Subscript[h, 1]), (Subscript[h, 2]), ⋯ of the polynomials G and H respectively are all positive and they alternately follow each other starting with (Subscript[g, 1]), i.e.

0 < (Subscript[g, 1]) < (Subscript[h, 1]) < (Subscript[g, 2]) < (Subscript[h, 2]) < ⋯

This theorem makes obviously possible to write a code which generally works faster then the one above. We are going however to state such theorems below which can serve as bases of much more efficient codes.

5 THE ROUTH–HURWITZ CRITERION

The best known criterion of stability is called the Routh–Hurwitz (sometime only Hurwitz) criterion, although the criterion given by Routh (1875) was different from the one below. (Later, we shall present another criterion which is much closer to Routh's original work.) The Routh–Hurwitz criterion has essentially been discovered by  Hurwitz  (1895), a student of Stodola, independently from Routh.

The Routh–Hurwitz matrix of polynomial (1) is defined as follows:

RH(p):=(({{(Subscript[a, 1]), (Subscript[a, 0]), 0, 0, 0, ⋯, 0, 0, 0, 0}, {(Subscript[a, 3]), (Subscript[a, 2]), (Subscript[a, 1]), (Subscript[a, 0]), 0, ⋯, 0, 0, 0, 0}, {⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮}, {0, 0, 0, 0, 0, 0, 0, (Subscript[a, N]), (Subscript[a, N - 1]), (Subscript[a, N - 2])}, {0, 0, 0, 0, 0, 0, 0, 0, 0, (Subscript[a, N])}}))

The elements of the above matrix can be written in a unified way with the following extension of the definition of the elements:

(Subscript[a, n]) := 0,    if    n < 0   or   n > N.

Now we have

RH(p)=\(\([a\_\(2  i - j\)]\)\_\(i = 1, 2, ⋯, N

            j = 1, 2, ⋯, N\)\)

Using the above notations the Routh–Hurwitz criterion (see [3, p. 559]) can be formulated as follows.

Theorem 3

Let us suppose that (Subscript[a, 0]) > 0.  Then the polynomial p in (1) is stable if and only if the chief minors of RH(p) are all positive, i.e.

(Subscript[Δ, n]):=det(({{(Subscript[a, 1]), (Subscript[a, 0]), ⋯, (Subscript[a, 2 - n])}, {(Subscript[a, 3]), (Subscript[a, 2]), ⋯, (Subscript[a, 4 - n])}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[a, 2  n - 1]), (Subscript[a, 2  n - 2]), ⋯, (Subscript[a, n])}}))>0     (n = 1, 2, ⋯, N) .

We shall give a function which tests the above criterion (without he assumption on positivity of (Subscript[a, 0])).

First we define a function which gives the Routh–Hurwitz matrix (3).

\(\(RH[polyvalue_, variable_] := 

      Module[{lc = 
Reverse[
CoefficientList[polyvalue, variable]], 

          deg = Exponent[polyvalue, variable], a}, 
a[k_ /; k < 0] := 0; 

        a[k_ /; 0 ≤ k ≤ deg] := 
          lc[\([k + 1]\)]\/Sign[First[lc]]; 
a[k_ /; k > deg] := 0; 

        Table[a[2  i - j], {i, deg}, {j, deg}]]		\)\)

The Routh–Hurwitz matrix of the polynomials

(Subscript[p, 10])[z_] := (z^3) + a (z^2) + b z + 1(Subscript[p, 11])[z_] := (z^4) + a (z^3) + b (z^2) + c z + 1

are as follows

TableForm[RH[(Subscript[p, 10])[z], z]]

TableForm[RH[(Subscript[p, 11])[z], z]]

The chief minors can be obtained in the following way.

\(Δ[polyvalue_, variable_, n_] := 

      Det[\(RH[polyvalue, variable]\)[\([Range[n], Range[n]]\)]]\)

Let us apply this function to the previous examples.

Δ[(Subscript[p, 10])[z], z, #] & /@ Range[3]

Δ[(Subscript[p, 11])[z], z, #] & /@ Range[4]

The following function tests the stability of the polynomial using the conditions of Theorem 3.

\(RouthHurwitz[polyvalue_, variable_] := 

      Module[
{deg = Exponent[polyvalue, variable], tf}, 

        tf = Array[Positive[Δ[polyvalue, variable, #]] &, deg, 
            1, And]; 
If[tf, "Stable", "Unstable", tf]]\)

If the coefficients are numerical values, then a well-defined answer will be obtained.

RouthHurwitz[#[z], z] & /@ 
      {(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 4]), (Subscript[p, 5]), (Subscript[p, 6]), (Subscript[p, 7]), (Subscript[p, 8]), (Subscript[p, 9])}

If there are symbolic numbers or parameters involved in the coefficients then the answer still contains relevant information.

RouthHurwitz[(Subscript[p, 10])[z], z]

RouthHurwitz[(Subscript[p, 11])[z], z]

In 1914 Liénard and Chippart discovered that the amount of calculation can significantly be reduced, namely, it is enough to calculate every second determinant. More precisely, the following statement (see e.g. [4, p. 96]) is true.

Theorem 4

Let us suppose that (Subscript[a, 0]) > 0.  Then the polynomial p in (1) is stable if and only if every coefficient of p and the chief minors

(Subscript[Δ, N - 1]), 
    (Subscript[Δ, N - 3]), (Subscript[Δ, N - 5]), 
    ⋯

are all positive.

The following function to test the Liénard–Chippart criterion is based upon the functions Stodola[]and Δ[]introduced above.

LienardChippart[polyvalue_, variable_] := 
      Module[{deg = Exponent[polyvalue, variable], tf}, 
        tf = Stodola[polyvalue, variable] && 
            And@@Positive[
                (Δ[polyvalue, variable, #1] &)/@
                  Range[deg - 1, 1, -2]]; 
        If[tf, "Stable", "Unstable", tf]]

Let us see some examples.

LienardChippart[#[z], z] & /@ {(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 4]), (Subscript[p, 5]), (Subscript[p, 6]), (Subscript[p, 7]), (Subscript[p, 8]), (Subscript[p, 9])}

LienardChippart[(Subscript[p, 10])[z], z]

LienardChippart[(Subscript[p, 11])[z], z]

6 EUCLIDEAN CRITERION

We obtain a necessary and sufficient condition for the stability problem using the principle of arguments and the theory of Cauchy indices of rational functions.

To the formulation of the final result let us consider polynomial (1) and introduce two polynomials in the following way.

(Subscript[f, 0]) (z):=(Subscript[a, 0]) (z^N)-(Subscript[a, 2]) (z^(N - 2))+ ⋯      (z ∈ ℂ),       (Subscript[f, 1]) (z):=(Subscript[a, 1]) (z^(N - 1))-(Subscript[a, 3]) (z^(N - 3))+ ⋯     (z ∈ ℂ).

From the long division (Euclidean) algorithm it follows that for the given polynomials (Subscript[f, 0]) and (Subscript[f, 1]) there exist a uniquely determined number m ∈ ℕ and nonzero polynomials (Subscript[f, 2]), (Subscript[f, 3]), ⋯, (Subscript[f, m]) and (Subscript[q, 2]), (Subscript[q, 3]), ⋯, (Subscript[q, m]) for which the following relations hold:

\(\(f\_0\) \((z)\) = 
        \(\(q\_1\) \((z)\)\) \(\(f\_1\) \((z)\)\) - \(f\_2\) \((z)\), 

      \(f\_1\) \((z)\) = 
        \(\(q\_2\) \((z)\)\) \(\(f\_2\) \((z)\)\) - 
          \(f\_3\) 
            \((z)\), 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ⋮\)\(\(f\_\(k - 1\)\) \((z)\) = 
            \(\(q\_k\) \((z)\)\) \(\(f\_k\) \((z)\)\) - 
              \(f\_\(k + 1\)\) 
                \((z)\), 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ⋮\) (Subscript[f, m - 1]) (z) = (Subscript[q, m]) (z) (Subscript[f, m]) (z),

where the degree of (Subscript[f, k]) is less than the degree of \(TraditionalForm
      \`\(f\_\(k - 1\)\ \ \(\((k = 2, 3, ⋯, m)\) . \)\ 
      \)\) The polynomial (-(Subscript[f, k + 1])) —resp. (Subscript[q, k])— is called the remainder—resp. quotient—in the k-th step   (k = 1, 2, ⋯, m) of the division algorithm (5).

The following statement is true (see [2, Theorem 6.7b]).

Theorem 5

The polynomial (1) is stable if and only if m = N in (5)  and the quotient polynomials (Subscript[q, k])   (k=1,2,⋯,m) are in the following form:

\(TraditionalForm
    \`\(q\_k\)(z) = 
      \(γ\_k\) 
        z\ \ \ \ \ \((z ∈ ℂ; \ k = 1, 2, 
          ⋯, m)\)\)

where

(Subscript[γ, k])>0     ( k = 1, 2, ⋯, m)

This criterion can be implemented in Mathematica as follows.

Euclidean[polyvalue_, variable_]:=
Module[{deg = Exponent[polyvalue, variable], op, f, tf},
\(op[{q_, d_, r_, v_}] := 

              {PolynomialQuotient[d, r, v], r, 

                \(-PolynomialRemainder[d, r, v]\), v}\);
f[k_] := 
              (∑_(i = 0)^(Quotient[deg - k, 2]))(((-1))^i) 
                  polyvalue[[-2 i - 1 - k]];
tf = And@@
                ((Length[#] == 2 && First[#] == 0 && Last[#] > 0) &/@
                    Together[
                      Simplify[
                        CoefficientList[
                          First[Transpose[
                              FixedPointList[op, 
                                {variable, f[0], f[1], variable}, 
                                SameTest -> (#2[[3]] == 0 &)]]], 
                          variable]]]);
If[tf, "Stable", "Unstable", tf]]

Euclidean[#[z], z] & /@ 
      {(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 4]), (Subscript[p, 5]), (Subscript[p, 6]), (Subscript[p, 7]), (Subscript[p, 8]), (Subscript[p, 9])}

Euclidean[(Subscript[p, 10])[z], z]

Euclidean[(Subscript[p, 11])[z], z]

We propose the reader to proceed with loading either the package Algebra`InequalitySolve` or the package Algebra`AlgebraicInequalities` and using them.

7 CONTINUED FRACTION CRITERION

Another necessary and sufficient condition for the stability problem may be given by means of continued fractions.

It is known that some real functions can be written in the following m-terminate continued fraction form:

\(b\_1\/\(1 + 
              \(\(b\_2\) z\)\/\(1 + \(\(b\_3\) z\)\/\(\ \ \ \ \ \ \ ⋮\ \ 

                    	1 + \(b\_m\) z\)\)\)\)

where (Subscript[b, k]) (k = 1, 2, , ⋯, m) are appropriate real numbers. It is obvious that there exist rational functions for which such a representation does not exist. Consider for example those which have the zero value at the point z = 0.

Now suppose that the rational function

\(TraditionalForm
    \`r(z) := 
      \(c\_0 + \(c\_1\) z + ⋯ + \(c\_s\) z\^s\)\/\(d\_0 + \(d\_1\) z + ⋯ + \(d\_t\) z\^t\)\)

has a representation of the form (8). Extend the definition of the coefficients (Subscript[c, k]) and (Subscript[d, k]) in (9) in the following way:

(Subscript[c, n]) := 0,    if n > s        and        (Subscript[d, n]) := 0,     if n > t .

It may be proved ([3, p. 547]) that in this case the numbers (Subscript[b, k]) in (8) satisfy the following recurrence relations:

((Subscript[c, n])^((-1))) := (Subscript[d, n]) ,      ((Subscript[c, n])^((0))) := 
            (Subscript[c, n])     (n = 0, 1, 2, ⋯),

(Subscript[b, k + 1]) := 
            (((Subscript[c, 0])^((k)))/((Subscript[c, 0])^((k - 1))))     (k = 0, 1, 2, 
                ⋯),

((Subscript[c, n])^((k + 1))) := 
            (Subscript[b, k + 1]) ((Subscript[c, n + 1])^((k - 1))) - 
              ((Subscript[c, n + 1])^((k)))      (k = 0, 1, 2, 
                    ⋯; n = 0, 1, 2, ⋯) 
                  .

The computation of (Subscript[b, k + 1])  (which in this manner is a variant of the original Routh algorithm) starts from the first two rows of ((Subscript[c, n])^((-1))) and ((Subscript[c, n])^((0))) (the coefficients of the denominator and that of the numerator) and then proceeds as indicated in the following table:

Table 1

From Theorems 12.7 b, 12.7 a and 12.6 c of [3] one obtains the continued fraction criterion:

Theorem 6

The polynomial

\(TraditionalForm
    \`p(z) := a\_0 + \(a\_1\) z + \(a\_2\) z\^2 + ⋯ + 
        \(a\_N\) z\^N\)

of degree N with real coefficients is stable if and only if the rational function

\(TraditionalForm
    \`\(a\_1 + \(a\_3\) z + \(a\_5\) z\^2 + ⋯\)\/\(a\_0 + \(a\_2\) z + \(a\_4\) z\^2 + ⋯\) = 
      \(∑\+\(k = 0\)\%\([\((N - 1)\)/2]\)\(a\_\(2  k + 1\)\) z\^k\)\/\(∑\+\(k = 0\)\%\([N/2]\)\(a\_\(2  k\)\) z\^k\)\)

(the Hurwitz alternant of the polynomial p) can be represented by such an N-terminating continued fraction of the form (8), in which every number (Subscript[b, k]) (k = 1, 2, ⋯, n) is positive.

A possible code of this criterion follows.

\(CFCriterion[polyvalue_, variable_] := 

      Module[
{cf = 
            CoefficientList[polyvalue, variable], 
\ \ \ \ \ \ \ deg, \ 
          cfde, cfnu, b, sQ, temp = {}, m = 1}, 
deg = Length[cf] - 1; 

        sQ = \(! FreeQ[NumberQ/@cf, False]\); 

        cfnu = Table[cf[\([k]\)], {k, 2, deg + 1, 2}]; 

        cfde = Table[cf[\([k]\)], {k, 1, deg + 1, 2}]; 

        While[\((\((First[cfnu] > 0\  ∧ \ m < deg + 1)\) 
              ∨ 
\ \ \ \ \ \((sQ ∧ m < deg + 1)\))\), 

          b[m] = First[cfnu]\/First[cfde]; 

          If[Length[cfnu] < Length[cfde], 
			cfnu = Join[cfnu, {0}]]; 

          temp = cfnu; 
cfnu = Drop[b[m]*cfde - cfnu, 1]; 
cfde = temp; 

          m = m + 1]; 

        If[m < deg, "Unstable", 
          If[\(! sQ\), 
            "Stable", 
\ \ \ \ \ Together[
              Simplify[Table[b[k], {k, 1, m - 1}]]]]]]\)

CFCriterion[#[z], z] & /@ 
        {(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 4]), (Subscript[p, 5]), (Subscript[p, 6]), (Subscript[p, 7]), (Subscript[p, 8]), (Subscript[p, 9])}

CFCriterion[(Subscript[p, 10])[z], z]

CFCriterion[(Subscript[p, 11])[z], z]

As it can be seen, in the case of paramteric coefficients those expressions are provided the positivity of which has to be checked.

8 A PACKAGE FOR SUMMARY

The code fragments above are summarized into a single function called StablePolynomialQ with the default value NegativeRealParts  of its option Method. Our aim was not to produce a final version of a package, as we think that further experiments (e.g. with automatic investigation of positivity with the new possibilities of version 3.0) are needed to produce a 'product' which can be used by everybody with confidence. We formulated the function in such a way that it gives an answer which is either True or False, or it gives an expression to be investigated further. This expression is given in different forms when different methods are used. The function itself is to be contained in a separate file attached to the electronic attachment of the journal. Below we show the results obtained when one investigates the same polynomials as above.

(Subscript[p, 3]) with different parameters

StablePolynomialQ[(Subscript[p, 3])[z, a], z];

StablePolynomialQ[(Subscript[p, 3])[z, a], z, Method -> NegativeRealParts];

StablePolynomialQ[(Subscript[p, 3])[z, a], z, Method -> Stodola]

StablePolynomialQ[#, z] & /@ polylist1

Results on polylist2

polylist2 = 
      {(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 4]), (Subscript[p, 5]), (Subscript[p, 6]), (Subscript[p, 7]), (Subscript[p, 8]), (Subscript[p, 9]), (Subscript[p, 10]), (Subscript[p, 11])};

Drop[StablePolynomialQ[#[z], z] & /@ polylist2, -2]

\(Drop[
      \(StablePolynomialQ[
#[z], z, Method -> NegativeRealParts] 
          &\)\ /@\ polylist2, \(-2\)]\)

StablePolynomialQ[#[z], z, Method -> Stodola] & /@ polylist2

StablePolynomialQ[#[z], z, Method -> RouthHurwitz] & /@ 
      polylist2

StablePolynomialQ[#[z], z, Method -> LienardChippart] & /@ 
      polylist2

StablePolynomialQ[#[z], z, Method -> Euclidean] & /@ 
      polylist2

StablePolynomialQ[#[z], z, Method -> ContinuedFraction] & /@ 
      polylist2

ACKNOWLEDGEMENT

The present work has partly been supported by the grants No. MKM 419/1995 and MKM 184/1996 of the Hungarian Ministry for Education and Culture. The authors benefited from the advice given by the referee.

REFERENCES

[1] Corless, R. M.: Simple engineering mathematics with Maple, MapleTech 3 (1) (1966) 2–7.

[2] Henrici, P.: Applied and Computational Complex Analysis. Vol. I. John  Wiley & Sons, New York, London, Sydney, Toronto (1974).

[3] Henrici, P.: Applied and Computational Complex Analysis. Vol. II. John  Wiley & Sons, New York, London, Sydney, Toronto (1977).

[4] Postnikov, M. M.: Stable polynomials. Nauka, Moscow (1981) (in Russian).

[5] Szili, L., Tóth, J.: Mathematics and Mathematica. ELTE Eötvös Kiadó, Budapest (1996) (in Hungarian).

[6] Wolfram, S.: The Mathematica Book. (Third edition) Wolfram Media–Cambridge University Press (1996).

[7] Pontryagin, L. S.: Ordinary differential equations, Nauka, Moscow, (1965) (in Russian)

ABOUT THE AUTHORS

János Tóth received his diploma in mathematics at the Loránd Eötvös University of Budapest in 1971. Currently he is head of the Department of Computer Science, Agricultural University, Gödöllõ, and participates in teaching at Loránd Eötvös University as well. His primary research interest is in deterministic and stochastic modelling of chemical reactions with emphasis on the reduction of the number of variables. He is the author of several books, including Mathematical Models of Chemical Reactions with P. Érdi copublished by Manchester University Press and Princeton University Press.

János Tóth
Department of Computer Science
Agricultural University
Páter K. u. 1.
H–2103 Gödöllõ, HUNGARY
tothj@mszi.gau.hu, jtoth@helka.iif.hu

László Szili received his diploma as a teacher of mathematics and physics at the Loránd Eötvös University of Budapest in 1978. Currently he is assistant professor at the Department of Numerical Analysis, Loránd Eötvös University, Budapest.  His primary research interest is in approximation theory on which topics he has published a dozen of papers. He is also interested in the theory, in the numerical methods and in the applications of differential equations.

László Szili
Department of Numerical Analysis
Loránd Eötvös University
Múzeum krt. 6–8.
H–1088 Budapest, HUNGARY
szili@ludens.elte.hu

András Zachár received his diploma in mechanical engineering at the Faculty for Agricultural Mechanical Engineering, Agricultural University, Gödöllõ. Currently he is a teaching assistant at the  Department of Computer Science, Agricultural University, Gödöllõ, and a PhD student in mathematical modelling at the Department of Physics and Process Control of the same university. His primary research interest is in modelling technical phenomena using ordinary and partial differential equations.

András Zachár
Department of Computer Science
Agricultural University
Páter K. u. 1.
H–2103 Gödöllõ, HUNGARY
azachar@mszi.gau.hu