Weierstrass (Durand-Kerner) method

The Weierstraß (Durand-Kerner) method (W (DK) method) is an iterative method for the simultaneous determination of all zeros of a univariate polynomial . It is named after Karl Weierstraß , who developed it as part of a proof of the fundamental theorem of algebra between 1859 and 1891, and E. Durand and Immo Kerner, who converted it into a computer algorithm in 1960 and 1966, respectively.

The procedure

Let it be a standardized univariate polynomial with complex coefficients and with a leading coefficient 1. According to the fundamental theorem of algebra, this has exactly n zeros and can be broken down into linear factors, ${\displaystyle p(X)=X^{n}+p_{n-1}X^{n-1}+\cdots +p_{0}}$ ${\displaystyle \xi _{1},\dots ,\xi _{n}\in \mathbb {C} }$

${\displaystyle p(X)=(X-\xi _{1})\cdots (X-\xi _{n}).\,}$

Each of the zeros can be formally isolated from this formula, it applies

${\displaystyle \xi _{k}=X-{\frac {p(X)}{\prod _{j=1,\dots ,n,\;j\neq k}(X-\xi _{j})}}.}$

This formula can be understood as a trivial fixed point iteration, the iteration

${\displaystyle z_{k}^{(i+1)}:=z_{k}^{(i)}-{\frac {p(z_{k}^{(i)})}{\prod _{j=1,\dots ,n,\;j\neq k}(z_{k}^{(i)}-\xi _{j})}}}$

obviously becomes stationary in the zero after the first iteration step. ${\displaystyle \xi _{k}}$

If the other zeros in the iteration rule are replaced by good approximations, a fixed point remains and every iteration that starts near this zero also converges against it. The iteration of the W (DK) method results when approximation sequences are now determined simultaneously for all zeros by means of this iteration rule, and the respective determined approximation of a zero is immediately included in the determination of the next approximations of the other zeros. ${\displaystyle \xi _{k}}$

At the beginning of each iteration step there are n complex numbers that are different in pairs . For the first step, these numbers can be chosen randomly but differently in pairs, in later steps these numbers stand for approximations of the zeros of p (X) . ${\displaystyle z_{1}^{(i)},\dots ,z_{n}^{(i)}\in \mathbb {C} }$ ${\displaystyle z_{1}^{(0)},\dots ,z_{n}^{(0)}\in \mathbb {C} }$

The polynomial that has these complex numbers as zeros is assigned to the tuple . From this polynomial the derivative according to the indefinite X is determined. It apply ${\displaystyle {\vec {z}}^{(i)}=(z_{1}^{(i)},\dots ,z_{n}^{(i)})}$ ${\displaystyle g(X):=(X-z_{1}^{(i)})\cdot \ldots \cdot (X-z_{n}^{(i)})}$ ${\displaystyle g'(X)}$

${\displaystyle g(z_{k}^{(i)})=0}$ and ${\displaystyle g'(z_{k}^{(i)})=\prod _{j\neq k}(z_{k}^{(i)}-z_{j}^{(i)}).}$

Of p (X) and of the derivative which are Weierstrass corrections , k  = 1, ...,  n is determined and the corrected tuple obtained as a result of the next iteration step and input. ${\displaystyle g_{{\vec {z}}^{(i)}}'(X)}$ ${\displaystyle w_{k}^{(i)}=-{\frac {p(z_{k}^{(i)})}{g'(z_{k}^{(i)})}}}$ ${\displaystyle {\vec {z}}^{(i+1)}={\vec {z}}^{(i)}+{\vec {w}}^{(i)}=(z_{1}^{(i)}+w_{1}^{(i)},\dots ,z_{n}^{(i)}+w_{n}^{(i)})}$

The iteration can e.g. B. canceled if the correction falls below a previously defined return accuracy. (The calculation accuracy should be slightly higher than this return accuracy.)

Variants of the procedure

The derivation given above determines the new approximations simultaneously, in parallel, from the old approximations. In the simplest implementation of the method, however, the updates and thus the new approximations are determined one after the other, sequentially. Therefore the idea is to use every new approximation of a zero immediately instead of the old one. This difference between parallel and sequential updating corresponds to the analogous procedure in the Jacobi and Gauß-Seidel methods for the iterative solution of linear systems of equations.

The sequential variant converges faster in many cases, but this speed advantage is difficult to grasp theoretically. With high polynomial degrees, the parallel variant allows the use of methods of fast polynomial multiplication, such as Karatsuba or Schönhage-Strassen in the calculation of the update.

example

The cubic equation is to be solved . The complex parameter is selected as the start tuple . The following first iteration steps result for the parallel variant of the iteration ${\displaystyle x^{3}-3x^{2}+3x-5=0}$ ${\displaystyle (z_{1},z_{2},z_{3})=(a^{0},a^{1},a^{2})}$ ${\displaystyle a=0{,}4+0{,}9\cdot i}$

It no.	z 1	z 2	z 3
0	+1.000000 + 0.000000 i	+0.400000 + 0.900000 i	−0.650000 + 0.720000 i
1	+1.360773 + 2.022230 i	−1.398213 - 0.693566 i	+3.037440 - 1.328664 i
2	+0.980963 + 1.347463 i	−0.335252 - 0.644069 i	+2.354289 - 0.703394 i
3	+0.317181 + 0.936495 i	+0.490016 - 0.966141 i	+2.192804 + 0.029647 i
4th	+0.209016 + 1.572742 i	+0.041206 - 1.527519 i	+2.749778 - 0.045223 i
5	+0.212971 + 1.394827 i	+0.184678 - 1.384565 i	+2.602351 - 0.010262 i
6th	+0.206531 + 1.374879 i	+0.206001 - 1.374653 i	+2.587468 - 0.000226 i
7th	+0.206300 + 1.374730 i	+0.206299 - 1.374730 i	+2.587401 - 0.000000 i
8th	+0.206299 + 1.374730 i	+0.206299 - 1.374730 i	+2.587401 + 0.000000 i

and for the sequential variant of the iteration

It no.	z 1	z 2	z 3
0	+1.000000 + 0.000000 i	+0.400000 + 0.900000 i	−0.650000 + 0.720000 i
1	+1.360773 + 2.022230 i	−0.365804 + 2.483787 i	−2.385807 - 0.028361 i
2	+2.659661 + 2.713714 i	+0.597676 + 0.822483 i	−0.631985 - 1.671566 i
3	+2.270389 + 0.387972 i	+0.131179 + 1.312808 i	+0.282054 - 1.501550 i
4th	+2.542817 - 0.015337 i	+0.204444 + 1.371609 i	+0.205573 - 1.372072 i
5	+2.587418 - 0.000012 i	+0.206300 + 1.374733 i	+0.206299 - 1.374730 i
6th	+2.587401 - 0.000000 i	+0.206299 + 1.374730 i	+0.206299 - 1.374730 i
7th	+2.587401 - 0.000000 i	+0.206299 + 1.374730 i	+0.206299 - 1.374730 i

In the first 4 iterations of both variants, the triple is moved apparently chaotically, from the 5 step onwards, more and more decimal places remain constant, iteration 8 in the parallel and 7 in the sequential case confirm the correctness of iteration 7 or 6 within the specified accuracy. This general behavior is characteristic of this method. ${\displaystyle (z_{1},z_{2},z_{3})}$

As a 3rd degree equation with real coefficients, p (X) has a real zero and a conjugate pair of complex zeros. The approximations satisfy this pattern. According to Vieta's theorem , e.g. B. the sum of all zeros corresponds to the negative of the coefficient of the second highest degree, i.e. 3. This is also confirmed in the approximations.

Justification as Newton's method

The - at least local - convergence of the Weierstrasse iteration results from its interpretation as a multi-dimensional Newton method . The equation system for this results from the comparison of the coefficients of the same degree in the desired identity

${\displaystyle g({\vec {z}})(X)=\prod _{j=1}^{n}(X-z_{j})=p(X).}$

Since the polynomials are normalized on both sides (the highest degree coefficient is 1 ), the identity in degree n is trivial and there are n equations for the n unknowns.

In general, this identity is not fulfilled. The correction in each step of the Newton iteration results from the requirement that the identity ${\displaystyle {\vec {w}}=(w_{1},\dots ,w_{n})}$

${\displaystyle g({\vec {z}}+{\vec {w}})(X)=\prod _{j=1}^{n}(X-z_{j}-w_{j})=p(X)}$

in the first order is fulfilled. The linear equation results from the Taylor expansion of the first order ${\displaystyle {\vec {w}}}$ ${\displaystyle {\vec {w}}}$

${\displaystyle g({\vec {z}})(X)-\sum _{k=1}^{n}w_{k}\,\prod _{j\neq k}(X-z_{j})=p(X).}$

Each individual correction can be derived from it by inserting to ${\displaystyle w_{k}}$ ${\displaystyle X=z_{k}}$

${\displaystyle w_{k}=-{\frac {p(z_{k})-g({\vec {z}})(z_{k})}{\prod _{j\neq k}(z_{k}-z_{j})}}=-{\frac {p(z_{k})}{\partial _{X}g({\vec {z}})(z_{k})}}}$

can be obtained, which results in the Weierstraß correction given above.

A global proof of convergence for this method was already given in (K. Weierstrass 1891) as an alternative proof for the fundamental theorem of algebra .

Error estimation using Gerschgorin circles

An error estimate and an alternative derivation of the method is given in the article on Gershgorin's theorem . According to this, all zeros of the polynomial p (X) are contained in the union of the circular disks in each iteration step . If the circular disks are disjoint in pairs, each contains exactly one zero. (A. Neumaier 2003) shows the same statement for the somewhat smaller circular disks ${\displaystyle D{\big (}z_{k}+w_{k},(n-1)\cdot |w_{k}|{\big )}}$ ${\displaystyle D{\big (}z_{k}+n/2\cdot w_{k},n/2\cdot |w_{k}|{\big )}.}$

literature

• Karl Weierstrass: New proof of the theorem that every whole rational function of a variable can be represented as a product of linear functions of the same variable. In: Meeting reports of the Royal Prussian Academy of Sciences in Berlin. 1891. online
• E. Durand: Equations du type F (x) = 0: Racines d'un polynomials . In: Masson et al. (Ed.): Solutions numeriques des equations algebriques. Vol. 1. 1960,
• Immo O. Kerner: A total step procedure for calculating the zeros of polynomials . In: Numerical Mathematics . 8, 1966, pp. 290-294.
• Marica Prešić: A convergence theorem for a method for simultaneous determination of all zeros of a polynomial . In: Publications de l'institut mathematique (Beograd) (NS) . 28 (42), 1980, pp. 158-168.
• MS Petkovic, C. Carstensen, M. Trajkovic: Weierstrass formula and zero-finding methods . In: Numerical Mathematics . 69, 1995, pp. 353-372.
• Bo Jacoby: Zero point for polynomial. CAE-nyt (a periodical for Dansk CAE Gruppe [Danish CAE Group]), 1988.
• Agnethe Knudsen: Numeriske Metoder. (lecture notes), Københavns Teknikum.
• Bo Jacoby: Numerisk løsning af ligninger. Bygningsstatiske meddelelser (Published by the Danish Society for Structural Science and Engineering) 63, no. 3-4, 1992, pp. 83-105.
• Xavier Gourdon: Combinatoire, Algorithmique et Geometry des Polynomials . Ecole Polytechnique, Paris 1996 ( Postscript [accessed October 15, 2006]).
• Victor Pan: Univariate Polynomial Root-Finding with Lower Computational Precision and Higher Convergence Rates . Tech-Report, City University of New York, May 2002.
• Bini, DA; Gemignani, L .; Pan, VY: Inverse power and Durand-Kerner iterations for univariate polynomial root-finding. Comput. Math. Appl. 47 (2004), no. 2-3, 447-459
• Arnold Neumaier: Enclosing clusters of zeros of polynomials . In: Journal of Computational and Applied Mathematics . 156, 2003.

Web links

• Jan Verschelde, The method of Weierstrass (also known as the Durand-Kerner method) , 2003.
• Roots Extraction from Polynomials: The Durand-Kerner Method - with a Java applet to follow the iteration step by step