| |
- 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 in the reals
which has complex conjugate zeroes usually fails (if not by hazard one
of the real zeroes is found). Laguerre's method finds one
zero of a quadratic in one step using the three values p(x), p'(x), p''(x)
(not the coefficients!) and is able to find complex
zeroes of real polynomials initialized with a real initial guess.
It is cubically convergent locally under the condition p'(z) ≠ 0).
Hence it can be the basis of finding all zeroes of a real polynomial.
We apply it here directly to a complex polynomial.
- The derivation of the method is a little bit dodgy and goes as follows:
For any polynomial p of degree n there holds
p'(x)/p(x) = ∑ i=1 to n 1/(x-zi)
with the n zeroes zi.
Now assume that x is nearer to z1 than to the other zeroes
and we want to approximate z1 first. Then we write
p'(x)/p(x) = 1/(x-z1) + (n-1)b(x) , b(x) = (1/(n-1)) ∑ i=2 to n 1/(x-zi) .
Then
(p'(x)/p(x))' = -1/(x-z1)2 + (n-1)b'(x).
Now we make the assumption that
b'(x) = -b(x)2
with a neglectable error. This might be the case
if the zeroes z3,...,zn are much farer away from x than z2.
Then, solving the two equations
- p'(x)/p(x) = 1/(x-z1) + (n-1)b(x)
- (p''(x)p(x) - p'(x)2)/p(x)2 = -1/(x-z1)2 -(n-1)b2(x)
for x-z1 we arrive at the Laguerre correction formula for x
x-z1 = (p(x)/p'(x))/( 1/n ± (n-1/n){1-(n/(n-1))p''(x)p(x)/(p'(x)2)}1/2 ) .
- If the evaluation of the radicand becomes questionable because of a too small |p'(x)|,
then we try a random perturbation of the current x, at most 10 times. We report
this as the number of restarts in the output.
- The assumption on b(x) made above is not as dodgy as it appears on the first look:
Write
b(x) = 1/(x-ξ(x)) with ξ(x) = x - 1/( (1/(n-1)) ∑i=2 to n1/(x-zi) ) .
Then
b'(x) = -b2(x)(1-ξ'(x)) .
The following table shows some values of ξ'(x)
Zeroes z(2) to z(6): 2,3,4,5,6
(2,0),(3,0),(4,0),(5,0),(6,0)
| x | abs((d/dx)xi(x) |
| 1.5 | 0.85300410 |
| 1.1 | 0.45623875 |
| 1.01 | 0.40843832 |
| 0.9 | 0.36043215 |
Zeroes z(2) to z(6): 4,5,6,7,8
| x | abs((d/dx)xi(x) |
| 1 | 9.59932804E-02 |
| 1.5 | 0.12425530 |
| 2 | 0.16858232 |
| 2.5 | 0.24556530 |
| 3 | 0.40364408 |
Hence we may hope the correction working well at least if x comes near to a zero separated
well by the others (in this case z(1)=1).
- We apply this correction iteratively and terminate if for two successive
approximations there holds
| xi+1-xi | <= C × ε × (| xi+1| +ε|
with the machine precision ε.
In order to cope with possible convergence problems we relax this by increasing C (originally set to n)
by 10 if the number of iterative steps becomes larger than a number maxcnt (originally set to 100)
and increase maxcnt by the factor 2, repeating this 8 times at most.
If termination did not occur then , we give up with a corresponding message.
- Having found a zero we use synthetic division in order to remove this from p
and diminish the degree. It is well known that this can introduce severe roundoff effects.
- Hence we try to find the zeroes in increasing absolute value, a heuristic first
introduced by Wilkinson. Mardens bound
B =2 max{|an-i/a0|1/(n-i)|}, i=0,..,n-1
is used for bounding all zeroes (a0 is the highest coefficient of the reflected
polynomial with the reciprocal zeroes) and
1/B = z0
is the first iterate.
- In addition, after acceptance of a zero, the try up to 3 steps of Newton's method (as long
as the absolute value of p decreases) in the
original, nondeflated polynomial in order to improve the precision.
- Although there is no proof of global convergence in practice this is found a
very reliable method for finding all zeroes.
- As a final error estimate we report the magnitude of the last Newton correction, multiplied by 2.
If there appears a -1, then this means that the correction failed because of a zero derivative.
|
| |
- The degree n.
- You have the choice: either you specify n zeroes and let the code compute the
coefficients and from these back the zeroes or you specify
- n+1 complex coefficients an, ... ,a0 (in this order)
- The 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.
an 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.
|