C202 Zeros of a Real Polynomial

Subroutine subprogram RMULLZ and DMULLZ compute the zeros of the polynomial

$P(z)\; =\; a$₀zⁿ+ a₁z^n-1+...+a_n-1z + a_n

of degree n with real coefficients $a$_k and $a$₀≠0 .

On computers other than CDC or Cray, only the double-precision version DMULLZ is available. On CDC and Cray computers, only the single-precision version RMULLZ is available.

Structure:

SUBROUTINE subprograms
User Entry Names : RMULLZ, DMULLZ
Files Referenced: Unit 6
External References: MTLMTR (N002), ABEND (Z035)

Usage:

For $t=R$ (type REAL), $t=D$ (type DOUBLE PRECISION),

    CALL tMULLZ(A,N,MAXIT,Z)

A

(type according to t) One-dimensional array of dimension (0:d), where $d\; \ge N$ , containing the coefficients $a$_k, (k = 0,1,...,n) .

N

( INTEGER) The degree n.

MAXIT

( INTEGER) The maximum number of iterations permitted.

Z

( COMPLEX for $t=R$ , COMPLEX*16 for $t=D$ ) One-dimensional array of length $\ge N$ . On exit, $Z(i)$ contains an approximation to the zero $z$_i , listed in roughly decreasing order of $|z$_i| .

Method:

The method of Muller (see Ref. 1) is used. This is based on iterated inverse quadratic interpolation followed by deflation to remove each zero as found.

Accuracy:

For well-conditioned polynomials (i.e. polynomials whose zeros are not unduly sensitive to small errors in the coefficients), the relative error of a computed zero of multiplicity m is of order $10-d/m$ where d is the machine precision expressed in decimal digits. For m>1, the m approximations to the single multiple zero are uniformly distributed on a small circle of radius of order $10-d/m$ around the exact zero. Therefore, if the polynomial is well-conditioned, the true value of the multiple zero will be close to the centre $(z$_k+1+...+z_k+m)/m of this circle.

Error handling:

Error C202.1: $a$₀= 0 .
Error C202.2: The number of iterations exceeds MAXIT.
In both cases, a message is written on Unit 6, unless subroutine MTLSET (N002) has been called. If the number of iterations exceeds MAXIT, those zeros which have not been found are set to $1020$ .

Notes:

For difficult cases which lead to too many iterations the following transformations may be applied, singly or together, to obtain a better-conditioned polynomial:

1. Reverse the order of the coefficients to obtain a polynomial whose zeros are $z$_i^-1 .
2. If the zeros $z$_i are clustered, or are too unsymmetrically positioned with respect to the origin, compute by synthetic division (see Ref. 3) the coefficients of the polynomial whose argument is $w=z-z$ , where $z=\; -a$₁/(n a₀) is the arithmetic mean of the zeros. The mean of the zeros of this new polynomial is situated at the origin, which is where the subprogram starts searching. Then, provided $|w$_i|<|z| for most i, $z$_i= w_i+z will be more accurate zeros.

References:

1. D.E. Muller, A method for solving algebraic equations using an automatic computer, MTAC (later renamed Math. Comp.) 10 (1956) 208--215.
2. J.W. Daniel, Correcting approximations to multiple roots of polynomials, Numer. Math. 9 (1966) 99--102.
3. F.B. Hildebrand, Introduction to numerical analysis, McGraw-Hill, New York (1956), Section 10.9.

C205