Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_zeros_poly_real (c02ag)

Purpose

nag_zeros_poly_real (c02ag) finds all the roots of a real polynomial equation, using a variant of Laguerre's method.

Syntax

[z, ifail] = c02ag(a, n, 'scal', scal)
[z, ifail] = nag_zeros_poly_real(a, n, 'scal', scal)

Description

nag_zeros_poly_real (c02ag) attempts to find all the roots of the n$n$th degree real polynomial equation
 P(z) = a0 zn + a1 zn − 1 + a2 zn − 2 + ⋯ + an − 1 z + an = 0 . $P(z) = a0 zn + a1 zn-1 + a2 zn-2 + ⋯ + an-1 z + 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
 L(zk) = zk + 1 − zk = ( − n P (zk) )/( P′ (zk) ± sqrt(H(zk)) ) , $L(zk) = zk+1 - zk = -n P (zk) P′ (zk) ± H(zk) ,$
where H(zk) = (n1)[(n1)(P(zk))2nP(zk)P(zk)]$H\left({z}_{k}\right)=\left(n-1\right)\left[\left(n-1\right){\left({P}^{\prime }\left({z}_{k}\right)\right)}^{2}-nP\left({z}_{k}\right){P}^{\prime \prime }\left({z}_{k}\right)\right]$ and z0${z}_{0}$ is specified.
The sign in the denominator is chosen so that the modulus of the Laguerre step at zk${z}_{k}$, viz. |L(zk)|$|L\left({z}_{k}\right)|$, 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,${z}_{1},{z}_{2},{z}_{3},\dots \text{}$, such that |P(zk + 1)| < |P(zk)|$|P\left({z}_{k+1}\right)|<|P\left({z}_{k}\right)|$ and ensures that zk + 1 + L(zk + 1)${z}_{k+1}+L\left({z}_{k+1}\right)$ ‘roughly’ lies inside a circular region of radius |F|$|F|$ about zk${z}_{k}$ known to contain a zero of P(z)$P\left(z\right)$; that is, |L(zk + 1)||F|$|L\left({z}_{k+1}\right)|\le |F|$, where F$F$ denotes the Fejér bound (see Marden (1966)) at the point zk${z}_{k}$. Following Smith (1967), F$F$ is taken to be min (B,1.445nR) $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(B,1.445nR\right)$, where B$B$ is an upper bound for the magnitude of the smallest zero given by
 B = 1.0001 × min ( sqrt(n) L (zk) ,|r1|,|an / a0|1 / n) , $B = 1.0001 × min( n L (zk) ,|r1|, | an / a0| 1/n ) ,$
r1${r}_{1}$ is the zero X$X$ of smaller magnitude of the quadratic equation
 ( P ′ ′ (zk) )/( 2 n (n − 1) ) X2 + ( P′ (zk) )/n X + (1/2) P(zk) = 0 $P′′ (zk) 2 n (n-1) X2 + P′ (zk) n X + 12 P(zk) = 0$
and the Cauchy lower bound R$R$ for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation
 |a0| zn + |a1| zn − 1 + |a2| zn − 2 + ⋯ + |an − 1| z − |an| = 0 . $|a0| zn + |a1| zn-1 + |a2| zn-2 +⋯+ |an-1| z - |an| = 0 .$
Starting from the origin, successive iterates are generated according to the rule zk + 1 = zk + L (zk) ${z}_{k+1}={z}_{k}+L\left({z}_{k}\right)$, for k = 1 , 2 , 3 ,$k=1,2,3,\dots \text{}$, and L(zk)$L\left({z}_{k}\right)$ is ‘adjusted’ so that |P(zk + 1)| < |P(zk)| $|P\left({z}_{k+1}\right)|<|P\left({z}_{k}\right)|$ and |L(zk + 1)| |F| $|L\left({z}_{k+1}\right)|\le |F|$. The iterative procedure terminates if P (zk + 1) $P\left({z}_{k+1}\right)$ is smaller in absolute value than the bound on the rounding error in P (zk + 1) $P\left({z}_{k+1}\right)$ and the current iterate zp = zk + 1 ${z}_{p}={z}_{k+1}$ is taken to be a zero of P(z)$P\left(z\right)$ (as is its conjugate zp ${\stackrel{-}{z}}_{p}$ if zp${z}_{p}$ is complex). The deflated polynomial (z) = P (z) / (zzp) $\stackrel{~}{P}\left(z\right)=P\left(z\right)/\left(z-{z}_{p}\right)$ of degree n1$n-1$ if zp${z}_{p}$ is real ( (z) = P (z) / ((zzp)(zzp)) $\stackrel{~}{P}\left(z\right)=P\left(z\right)/\left(\left(z-{z}_{p}\right)\left(z-{\stackrel{-}{z}}_{p}\right)\right)$ of degree n2$n-2$ if zp${z}_{p}$ is complex) is then formed, and the above procedure is repeated on the deflated polynomial until n < 3$n<3$, whereupon the remaining roots are obtained via the ‘standard’ closed formulae for a linear (n = 1$n=1$) or quadratic (n = 2$n=2$) equation.
Note that nag_zeros_quadratic_real (c02aj), nag_zeros_cubic_real (c02ak) and nag_zeros_quartic_real (c02al) can be used to obtain the roots of a quadratic, cubic (n = 3$n=3$) and quartic (n = 4$n=4$) polynomial, respectively.

References

Marden M (1966) Geometry of polynomials Mathematical Surveys 3 American Mathematical Society, Providence, RI
Smith B T (1967) ZERPOL: a zero finding algorithm for polynomials using Laguerre's method Technical Report Department of Computer Science, University of Toronto, Canada
Thompson K W (1991) Error analysis for polynomial solvers Fortran Journal (Volume 3) 3 10–13
Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

Parameters

Compulsory Input Parameters

1:     a(n + 1${\mathbf{n}}+1$) – double array
if a is declared with bounds (0 : n)$\left(0:{\mathbf{n}}\right)$, then a(i) ${\mathbf{a}}\left(\mathit{i}\right)$ must contain ai ${a}_{\mathit{i}}$ (i.e., the coefficient of zni ${z}^{n-\mathit{i}}$) , for i = 0,1,,n$\mathit{i}=0,1,\dots ,n$.
Constraint: a(1) 0.0 ${\mathbf{a}}\left(1\right)\ne 0.0$.
2:     n – int64int32nag_int scalar
n$n$, the degree of the polynomial.
Constraint: n1${\mathbf{n}}\ge 1$.

Optional Input Parameters

1:     scal – logical scalar
Indicates whether or not the polynomial is to be scaled. See Section [Further Comments] for advice on when it may be preferable to set scal = false${\mathbf{scal}}=\mathbf{false}$ and for a description of the scaling strategy.
Default: true$\mathbf{true}$

w

Output Parameters

1:     z(2$2$,n) – double array
The real and imaginary parts of the roots are stored in z(1,i) ${\mathbf{z}}\left(1,\mathit{i}\right)$ and z(2,i) ${\mathbf{z}}\left(2,\mathit{i}\right)$ respectively, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$. Complex conjugate pairs of roots are stored in consecutive pairs of elements of z; that is, z(1,i + 1) = z(1,i) ${\mathbf{z}}\left(1,\mathit{i}+1\right)={\mathbf{z}}\left(1,i\right)$ and z(2,i + 1) = z(2,i) ${\mathbf{z}}\left(2,i+1\right)=-{\mathbf{z}}\left(2,i\right)$.
2:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, a(0) = 0.0${\mathbf{a}}\left(0\right)=0.0$, or n < 1${\mathbf{n}}<1$.
ifail = 2${\mathbf{ifail}}=2$
The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG, as some basic assumption for the arithmetic has been violated. See also Section [Further Comments].
ifail = 3${\mathbf{ifail}}=3$
Either overflow or underflow prevents the evaluation of P(z)$P\left(z\right)$ near some of its zeros. This error is very unlikely to occur. If it does, please contact NAG. See also Section [Further Comments].

Accuracy

All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed. See also Section [Example].

If scal = true${\mathbf{scal}}=\mathbf{true}$, then a scaling factor for the coefficients is chosen as a power of the base b$b$ of the machine so that the largest coefficient in magnitude approaches thresh = bemaxp$\mathit{thresh}={b}^{{e}_{\mathrm{max}}-p}$. You should note that no scaling is performed if the largest coefficient in magnitude exceeds thresh$\mathit{thresh}$, even if scal = true${\mathbf{scal}}=\mathbf{true}$. (b$b$, emax${e}_{\mathrm{max}}$ and p$p$ are defined in Chapter X02.)
However, with scal = true${\mathbf{scal}}=\mathbf{true}$, overflow may be encountered when the input coefficients a0,a1,a2,,an${a}_{0},{a}_{1},{a}_{2},\dots ,{a}_{n}$ vary widely in magnitude, particularly on those machines for which b(4p) ${b}^{\left(4p\right)}$ overflows. In such cases, scal should be set to false and the coefficients scaled so that the largest coefficient in magnitude does not exceed b(emax2p)${b}^{\left({e}_{\mathrm{max}}-2p\right)}$.
Even so, the scaling strategy used by nag_zeros_poly_real (c02ag) is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, you are recommended to scale the independent variable (z)$\left(z\right)$ 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 dP(cz)$dP\left(cz\right)$ for some suitable values of c$c$ and d$d$. For example, if the original polynomial was P(z) = 2100 + 2100z20$P\left(z\right)={2}^{-100}+{2}^{100}{z}^{20}$, then choosing c = 210$c={2}^{-10}$ and d = 2100$d={2}^{100}$, for instance, would yield the scaled polynomial 1 + z20$1+{z}^{20}$, which is well-behaved relative to overflow and underflow and has zeros which are 210${2}^{10}$ times those of P(z)$P\left(z\right)$.
If the function fails with ${\mathbf{ifail}}={\mathbf{2}}$ or 3${\mathbf{3}}$, 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 ${n}_{R}$ denote the number of roots found before the failure occurred. Then z(1,n) ${\mathbf{z}}\left(1,n\right)$ and z(2,n) ${\mathbf{z}}\left(2,n\right)$ contain the real and imaginary parts of the first root found, z(1,n1) ${\mathbf{z}}\left(1,n-1\right)$ and z(2,n1) ${\mathbf{z}}\left(2,n-1\right)$ contain the real and imaginary parts of the second root found, , z(1, n nR + 1) $\dots ,{\mathbf{z}}\left(1,n-{n}_{R}+1\right)$ and z(2, n nR + 1 ) ${\mathbf{z}}\left(2,n-{n}_{R}+1\right)$ contain the real and imaginary parts of the nR ${n}_{R}$th root found. After the failure has occurred, the remaining 2 × (nnR) $2×\left(n-{n}_{R}\right)$ elements of z contain a large negative number (equal to 1 / (x02am() × sqrt(2)) $-1/\left(\mathbf{x02am}\left(\right)×\sqrt{2}\right)$).

Example

```function nag_zeros_poly_real_example
a = [1;
2;
3;
4;
5;
6];
n = 5;
[z, ifail] = nag_zeros_poly_real(a, int64(n));

% Estimate the accuracy of the computed roots using a perturbation analysis

% Form the coefficients of the perturbed polynomial
abar = a;
epsbar = 3*x02aj;
for i=1:2:n+1
abar(i) = (1+epsbar)*abar(i);
end
for i=2:2:n+1
abar(i) = (1-epsbar)*abar(i);
end

% Compute the roots of the perturbed polynomial
[zbar, ifail] = nag_zeros_poly_real(abar, int64(n));

% Perform error analysis

% Initialise markers to 0 (unmarked)
m =zeros(n, 1);
r =zeros(n, 1);
rmax = x02al;

for i=1:n
deltai = rmax;
r1 = abs(complex(z(1,i), z(2,i)));

% Loop over all perturbed roots (stored in zbar)
for j=1:n
% Compare the current unperturbed root to all unmarked perturbed roots.
if m(j) == 0
r2 = abs(complex(zbar(1,j), zbar(2,j)));
deltac = abs(r1-r2);

if deltac < deltai
deltai = deltac;
jmin = j;
end
end
end

% Mark the selected perturbed root
m(jmin) = 1;

% Compute the relative error
if r1 ~= 0
r3 = abs(complex(zbar(1,jmin),zbar(2,jmin)));
di = min(r1,r3);
r(i) = max(deltai/max(di,deltai/rmax),x02aj);
end

end

fprintf('\n   Computed Roots of Polynomial       Error Estimate\n');
fprintf('                                    (machine dependent)\n');
for i=1:n
fprintf(' z = %12.4e%+12.4ei       %9.1e\n', z(1,i), z(2,i), r(i));
end
```
```

Computed Roots of Polynomial       Error Estimate
(machine dependent)
z =  -1.4918e+00 +0.0000e+00i         1.2e-15
z =   5.5169e-01 +1.2533e+00i         1.1e-16
z =   5.5169e-01 -1.2533e+00i         1.1e-16
z =  -8.0579e-01 +1.2229e+00i         1.1e-16
z =  -8.0579e-01 -1.2229e+00i         1.1e-16

```
```function c02ag_example
a = [1;
2;
3;
4;
5;
6];
n = 5;
[z, ifail] = c02ag(a, int64(n));

% Estimate the accuracy of the computed roots using a perturbation analysis

% Form the coefficients of the perturbed polynomial
abar = a;
epsbar = 3*x02aj;
for i=1:2:n+1
abar(i) = (1+epsbar)*abar(i);
end
for i=2:2:n+1
abar(i) = (1-epsbar)*abar(i);
end

% Compute the roots of the perturbed polynomial
[zbar, ifail] = c02ag(abar, int64(n));

% Perform error analysis

% Initialise markers to 0 (unmarked)
m =zeros(n, 1);
r =zeros(n, 1);
rmax = x02al;

for i=1:n
deltai = rmax;
r1 = abs(complex(z(1,i), z(2,i)));

% Loop over all perturbed roots (stored in zbar)
for j=1:n
% Compare the current unperturbed root to all unmarked perturbed roots.
if m(j) == 0
r2 = abs(complex(zbar(1,j), zbar(2,j)));
deltac = abs(r1-r2);

if deltac < deltai
deltai = deltac;
jmin = j;
end
end
end

% Mark the selected perturbed root
m(jmin) = 1;

% Compute the relative error
if r1 ~= 0
r3 = abs(complex(zbar(1,jmin),zbar(2,jmin)));
di = min(r1,r3);
r(i) = max(deltai/max(di,deltai/rmax),x02aj);
end

end

fprintf('\n   Computed Roots of Polynomial       Error Estimate\n');
fprintf('                                    (machine dependent)\n');
for i=1:n
fprintf(' z = %12.4e%+12.4ei       %9.1e\n', z(1,i), z(2,i), r(i));
end
```
```

Computed Roots of Polynomial       Error Estimate
(machine dependent)
z =  -1.4918e+00 +0.0000e+00i         1.2e-15
z =   5.5169e-01 +1.2533e+00i         1.1e-16
z =   5.5169e-01 -1.2533e+00i         1.1e-16
z =  -8.0579e-01 +1.2229e+00i         1.1e-16
z =  -8.0579e-01 -1.2229e+00i         1.1e-16

```