2 Specification

3 Description

nag_zeros_real_poly (c02agc) attempts to find all the roots of the nth degree real polynomial equation

Pz=a0zn+a1zn-1+a2zn-2+⋯+an-1z+an=0.

The roots are located using a modified form of Laguerre's method, originally proposed by Smith (1967).

The method of Laguerre (see Wilkinson (1965)) can be described by the iterative scheme

Lzk=zk+1-zk=-nPzkP′zk±Hzk,

where Hzk=n-1n-1P′zk2-nPzkP′′zk and z0 is specified.

The sign in the denominator is chosen so that the modulus of the Laguerre step at zk, viz. Lzk, is as small as possible. The method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots.

The function generates a sequence of iterates z1,z2,z3,…, such that Pzk+1<Pzk and ensures that zk+1+Lzk+1 ‘roughly’ lies inside a circular region of radius F about zk known to contain a zero of Pz; that is, Lzk+1≤F, where F denotes the Fejér bound (see Marden (1966)) at the point zk. Following Smith (1967), F is taken to be minB,1.445nR, where B is an upper bound for the magnitude of the smallest zero given by

B=1.0001×minnLzk,r1,an/a01/n,

r1 is the zero X of smaller magnitude of the quadratic equation

P′′zk2nn-1X2+P′zknX+12Pzk=0

and the Cauchy lower bound R for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation

a0zn+a1zn-1+a2zn-2+⋯+an-1z-an=0.

Starting from the origin, successive iterates are generated according to the rule
zk+1=zk+Lzk, for
k=1,2,3,…, and
Lzk is ‘adjusted’ so that
Pzk+1<Pzk and
Lzk+1≤F. The iterative procedure terminates if
Pzk+1 is smaller in absolute value than the bound on the rounding error in
Pzk+1 and the current iterate
zp=zk+1 is taken to be a zero of Pz (as is its conjugate
z-p if zp is complex). The deflated polynomial
P~z=Pz/z-zp of degree n-1 if zp is real
(P~z=Pz/z-zpz-z-p of degree n-2 if zp is complex) is then formed, and the above procedure is repeated on the deflated polynomial until n<3, whereupon the remaining roots are obtained via the ‘standard’ closed formulae for a linear (n=1) or quadratic (n=2) equation.

On exit: the real and imaginary parts of the roots are stored in
z[i].re and z[i].im respectively, for i=0,1,…,n-1. Complex conjugate pairs of roots are stored in consecutive pairs of z; that is, z[i+1].re=z[i].re and z[i+1].im=-z[i].im.

6 Error Indicators and Warnings

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.

NE_POLY_NOT_CONV

The iterative procedure has failed to converge.
This error is very unlikely to occur. If it does, please contact NAG
immediately, as some basic assumption for the arithmetic has been violated.

NE_POLY_OVFLOW

nag_zeros_real_poly (c02agc) cannot evaluate Pz near some of its zeros without overflow. If this message occurs please contact NAG.

NE_POLY_UNFLOW

nag_zeros_real_poly (c02agc) cannot evaluate Pz near some of its zeros without underflow. If this message occurs please contact NAG.

NE_REAL_ARG_EQ

On entry, a[0]=0.0.
Constraint: a[0]≠0.0.

7 Accuracy

All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed.

8 Further Comments

If scale=Nag_TRUE, then a scaling factor for the coefficients is chosen as a power of the base b of the machine so that the largest coefficient in magnitude approaches thresh=bemax-p. You should note that no scaling is performed if the largest coefficient in magnitude exceeds thresh, even if scale=Nag_TRUE. (b, emax and p are defined in Chapter x02.)

However, with scale=Nag_TRUE, overflow may be encountered when the input coefficients a0,a1,a2,…,an vary widely in magnitude, particularly on those machines for which b4p overflows. In such cases, scale should be set to Nag_FALSE and the coefficients scaled so that the largest coefficient in magnitude does not exceed bemax-2p.

Even so, the scaling strategy used by nag_zeros_real_poly (c02agc) is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, you are recommended to scale the independent variable z so that the disparity between the largest and smallest coefficient in magnitude is reduced. That is, use the function to locate the zeros of the polynomial dPcz for some suitable values of c and d. For example, if the original polynomial was Pz=2-100+2100z20, then choosing c=2-10 and d=2100, for instance, would yield the scaled polynomial 1+z20, which is well-behaved relative to overflow and underflow and has zeros which are 210 times those of Pz.

If the function fails with fail.code=NE_POLY_NOT_CONV, NE_POLY_OVFLOW or NE_POLY_UNFLOW, then the real and imaginary parts of any roots obtained before the failure occurred are stored in z in the reverse order in which they were found.
Let nR denote the number of roots found before the failure occurred. Then z[n-2].re and z[n-2].im contain the real and imaginary parts of the first root found, z[n-3].re and z[n-3].im contain the real and imaginary parts of the second root found, …,z[n-nR-1].re and z[n-nR-1].im contain the real and imaginary parts of the nRth root found. After the failure has occurred, the remaining nR elements of z have identical real and imaginary parts, equal to the large relative number -1.0/nag_real_safe_small_number×2.