c02ag:: Zeros of Polynomials (NAG Toolbox)

The function generates a sequence of iterates

z_{1}, z_{2}, z_{3}, \dots

, such that

|P (z_{k + 1})| < |P (z_{k})|

and ensures that

z_{k + 1} + L (z_{k + 1})

‘roughly’ lies inside a circular region of radius

|F|

about

z_{k}

known to contain a zero of

P (z)

; that is,

|L (z_{k + 1})| \leq |F|

, where

F

denotes the Fejér bound (see Marden (1966)) at the point

z_{k}

. Following Smith (1967),

F

is taken to be

\min (B, 1.445 n R)

, where

B

is an upper bound for the magnitude of the smallest zero given by

B = 1.0001 \times \min (\sqrt{n} L (z_{k}), |r_{1}|, {|a_{n} / a_{0}|}^{1 / n}),

r_{1}

is the zero

X

of smaller magnitude of the quadratic equation

\frac{P^{''} (z_{k})}{2 n (n - 1)} X^{2} + \frac{P^{'} (z_{k})}{n} X + \frac{1}{2} P (z_{k}) = 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

|a_{0}| z^{n} + |a_{1}| z^{n - 1} + |a_{2}| z^{n - 2} + \dots + |a_{n - 1}| z - |a_{n}| = 0 .

Starting from the origin, successive iterates are generated according to the rule

z_{k + 1} = z_{k} + L (z_{k})

, for

k = 1, 2, 3, \dots

, and

L (z_{k})

is ‘adjusted’ so that

|P (z_{k + 1})| < |P (z_{k})|

and

|L (z_{k + 1})| \leq |F|

. The iterative procedure terminates if

P (z_{k + 1})

is smaller in absolute value than the bound on the rounding error in

P (z_{k + 1})

and the current iterate

z_{p} = z_{k + 1}

is taken to be a zero of

P (z)

(as is its conjugate

{\bar{z}}_{p}

z_{p}

is complex). The deflated polynomial

\tilde{P} (z) = P (z) / (z - z_{p})

of degree

n - 1

z_{p}

is real (

\tilde{P} (z) = P (z) / ((z - z_{p}) (z - {\bar{z}}_{p}))

of degree

n - 2

z_{p}

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.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

scal = 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 = b^{e_{\max} - p}

. You should note that no scaling is performed if the largest coefficient in magnitude exceeds

thresh

, even if

scal = true

. (

b

e_{\max}

and

p

are defined in Chapter X02.)

However, with

scal = true

, overflow may be encountered when the input coefficients

a_{0}, a_{1}, a_{2}, \dots, a_{n}

vary widely in magnitude, particularly on those machines for which

b^{(4 p)}

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^{(e_{\max} - 2 p)}

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)

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

d P (c z)

for some suitable values of

c

and

d

. For example, if the original polynomial was

P (z) = 2^{- 100} + 2^{100} z^{20}

, then choosing

c = 2^{- 10}

and

d = 2^{100}

, for instance, would yield the scaled polynomial

1 + z^{20}

, which is well-behaved relative to overflow and underflow and has zeros which are

2^{10}

times those of

P (z)

If the function fails with

ifail = 2

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

n_{R}

denote the number of roots found before the failure occurred. Then

z (1, n)

and

z (2, n)

contain the real and imaginary parts of the first root found,

z (1, n - 1)

and

z (2, n - 1)

contain the real and imaginary parts of the second root found,

\dots, z (1, n - n_{R} + 1)

and

z (2, n - n_{R} + 1)

contain the real and imaginary parts of the

n_{R}

th root found. After the failure has occurred, the remaining

2 \times (n - n_{R})

elements of z contain a large negative number (equal to

- 1 / (x02am () \times \sqrt{2})

Example

function c02ag_example


fprintf('c02ag example results\n\n');

a = [1;     2;     3;     4;     5;     6];
n = int64(5);

ex1(n,a);

ex2(n,a);



function ex1(n,a)

  [z, ifail] = c02ag(a, n);
  disp('Computed roots of polynomial:');
  cz = z(1,:) + i*z(2,:);
  disp(cz');

function ex2(n,a)

  [z, ifail] = c02ag(a, 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 j=1:2:n+1
    abar(j) = (1+epsbar)*abar(j);
  end
  for j=2:2:n+1
    abar(j) = (1-epsbar)*abar(j);
  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;

  cz = z(1,:) + i*z(2,:);
  czbar = zbar(1,:) + i*zbar(2,:);
  for j=1:n
    deltai = rmax;
    r1 = abs(cz(j));

    % Loop over all perturbed roots (stored in zbar)
    for k=1:n
      % Compare the current unperturbed root to all unmarked perturbed roots.
      if m(k) == 0
        r2 = abs(czbar(k));
        deltac = abs(r1-r2);
        
        if deltac < deltai
          deltai = deltac;
          jmin = k;
        end
      end
    end

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

    % Compute the relative error
    if r1 ~= 0
      r3 = abs(czbar(jmin));
      di = min(r1,r3);
      r(j) = max(deltai/max(di,deltai/rmax),x02aj);
    end
  end
  fprintf('\nComputed Roots of Polynomial   Error Estimate\n');
  fprintf('                              (machine dependent)\n');
  for j=1:n
    if (z(2,j)<0)
      fprintf('%9.4f - %7.4fi           %8.2e\n',z(1,j),-z(2,j),r(j));
    else
      fprintf('%9.4f + %7.4fi           %8.2e\n',z(1,j),z(2,j),r(j));
    end
  end

c02ag example results

Computed roots of polynomial:
  -1.4918 + 0.0000i
   0.5517 - 1.2533i
   0.5517 + 1.2533i
  -0.8058 - 1.2229i
  -0.8058 + 1.2229i


Computed Roots of Polynomial   Error Estimate
                              (machine dependent)
  -1.4918 +  0.0000i           1.04e-15
   0.5517 +  1.2533i           1.62e-16
   0.5517 -  1.2533i           1.62e-16
  -0.8058 +  1.2229i           1.52e-16
  -0.8058 -  1.2229i           1.52e-16

On entry,	$a (0) = 0.0$ ,
or	$n < 1$ .

NAG Toolbox: nag_zeros_poly_real (c02ag)

▸▿ Contents

Purpose

Syntax

Description