| |
- Newton's method finds the zero of an affine linear function
in one step and this in connection with Taylor's formula is the basis of
the proof of its local quadratic convergence in the general case.
Newton's method, applied to a real polynomial p in the reals
which has complex conjugate zeroes fails (if not by hazard one
of the real zeroes is found) but cannot find any complex zeros. Bairstow's method overcomes this by
computing real quadratic factors of the polynomial which of course
will combine a complex conjugate pair of zeros if such are present but
otherwise may also combine two nearby real zeros.
- We assume here for simplicity that
p(x) = ∑ i=0 to n aixn-i
with a0an ≠ 0
- The unknowns in one main step of this method are the coefficients
r, q of a quadratic
x2 - r x - q .
Given some arbitrary r, q the factorization
p(x) = w(x)(x2 - r x - q) + f1 x + f2
will have nonzero coefficients f1, f2.
And, of course, f1, f2 will be functions of
r, q.
- The task makes sense for n >= 3 only, of course, for
n=1, 2 the zeros can be computed directly in a numerically safe way.
- If we let
w(x) = ∑ i=0 to n-2 bixn-2-i
then the following recursion computes the bi, f1, f2:
- b0 = a0 , b1 = b0 r + a1
- bi = bi-2 q + bi-1 r + ai, i=2,...,n-2
- f1 = bn-3 q + bn-2 r + an-1
- f2 = bn-2 q + an
- Hence this boils down to the nonlinear system in two unknowns
- f1(q,r) = 0
- f2(q,r) = 0
- We solve this system using the damped Newton's method (Bairstow himself used the full step
method only), since this one has a much better convergence for bad initial guesses, and indeed, a good
initial guess will not be handy here. The required Jacobian can easily be obtained by differentiating
the recursion above w.r.t. r, q.
- Having found a factor of p the process is repeated using w and decreasing the degree by 2
in this way until n=1 or n=2 is obtained.
- It is well known that this so called deflation may introduce large errors in the computed factors
and hence the zeros. Therefore we try to find factors with zeros near the origin first. To this end
we first compute a lower bound 1/B for the absolute values of the zeros of p
(0 is excluded by construction)
using Marden's bound
B=2 max{ |an-i/an|1/i: i=1,...,n }
for the reflected polynomial. At z=1/B we evaluate the second order Taylor expansion of p
using Horner's extended scheme and take this as the first guess for the factor after normalization.
Should p'' be very small there, then we take
r = q = 0
as initial guess.
This is our type of startvalue indicated as either 1 or 0 in the output.
If this fails within our tests, then we try a second initial value, namely
r=2z, q=-2z2
corresponding to a complex conjugate zero near zero, indicated as type 2.
- Having found a quadratic factor, we don't use it directly for deflation but add (the same) Newton's
process with the original nondeflated polynomial in order to weaken the effects of deflation even more.
This however might also generate a failure (false convergence to a factor already found earlier
for example), hence we give you the choice to suppress this feature.
- It may be that the damped Newton's method runs into trouble because of an almost
singular Jacobian. Then we perform a small random perturbation of the current values of r, q
and restart. This is repeated up to ten times, afterwards we give up with a corresponding message.
- We try to solve the nonlinear system with a very high relative precision in
r, q requiring
||(r,q)k+1-(r,q)k|| <= C ε (||(r,q)k+1||+ε)
k being the iteration counter, with the machine precision ε and
C initially set ot n within maxsteps=100 steps. If this fails, then we increase
C by 10 and maxsteps by 2 , this up to 8 times. If this does not help,
we give up with a corresponding message.
- With these features the method is more robust than the original one, but
nevertheless not nearly so as Hirano's.
There occurs a problem for n odd. There will be at least
one single real root which cannot be incorporated into such a quadratic and depending
on the quality of the automatic initial guess Newton's method might hang up without
finding a solution. (Try x2n+1-1!)
Therefore we provide the possibility to compute in this case one real root separately
before starting Bairstows method. In this case we use the Brent-Decker method for this root and
deflate it as the first. A similar problem occurs if for example there
is a real root near to a complex conjugate pair and all other real zeros are ''far away''.
|
| |
- The degree n, which is restricted by 50.
- You have the choice: either you specify n complex zeroes and let the code compute the
coefficients and from these back the zeroes (don't forget to use pairs of complex conjugate zeroes here, the coefficients
are assumed as real!)
- or you specify
- n+1 real coefficients a0, ... ,an (in this order)
- Complex numbers must be written as pairs in parentheses (re,im)
denoting re + sqrt(-1)*im.
- In order to avoid unnecessary preparative work we require
an*a0 ≠ 0.
a0 is the coefficient of zn.
- You may require a contour plot of |p(z)|. Then you must decide on a bounding box
given by a parameter b, which defines the size of the area to be plotted:
|p| will be shown over
[remin-b,remax+b]x [immin-b,immax+b],
where re/immin, re/immax are the algebraically smallest
resp. largest real and imaginary part of a zero.
Sometimes we recompute the area internally in order to obtain a useful picture.
- You might wish to avoid trouble for odd n by eliminating one real root
beforehand.
- You might want to suppress the purging of factors in the original polynomial.
|
| |
- You get a list of messages reporting possible trouble (or success).
- You get a list of the zeroes obtained so far, with the number of Newton steps and
the number of restarts tried for the corresponding factor, the successful initial guess
(as ''starttype''), and as an error estimation two times the absolute value of Newton's correction computed
from the original polynomial with the final value, in complex arithmetic.
- If you specified the zeroes originally,
this input.
- If you requested this, a contourplot of |p| around the zeroes.
This is adapted internally to get a reasonable representation (for example if all
zeroes are real but you did choose a large b.)
|