Durand–Kerner method

In numerical analysis, the Durand–Kerner method established 1960–66 and named after E. Durand and Immo Kerner, also called the method of Weierstrass, established 1859–91 and named after Karl Weierstrass, is a root-finding algorithm for solving polynomial equations. In other words, the method can be used to solve numerically the equation

ƒ(x) = 0

where ƒ is a given polynomial, which can be taken to be scaled so that the leading coefficient is 1.

Explanation

The explanation is for equations of degree four. It is easily generalized to other degrees.

Let the polynomial ƒ be defined by

ƒ(x) = x⁴ + ax³ + bx² + cx + d

for all x.

The known numbers a, b, c, d are the coefficients.

Let the (complex) numbers P,Q,R,S be the roots of this polynomial ƒ.

Then

ƒ(x) = (x − P)(x − Q)(x − R)(x − S)

for all x. One can isolate the value P from this equation,

$P=x-\frac{f(x)}{(x-Q)(x-R)(x-S)}.$

The substitution

$x:=x-\frac{f(x)}{(x-Q)(x-R)(x-S)}$

is a strongly stable fixed point iteration in that every initial point x ≠ Q,R,S delivers after one iteration the root P.

If one replaces the zeros Q, R and S by approximations q ≈ Q, r ≈ R, s ≈ S, such that q,r,s are not equal to P, then P is still a fixed point of the perturbed fixed point iteration since

$P-\frac{f(P)}{(P-q)(P-r)(P-s)} = P - 0 = P.$

Note that the denominator is still different from zero. This fixed point iteration is a contraction mapping for x around P.

The clue to the method now is to combine the fixed point iteration for P with similar iterations for Q,R,S into a simultaneous iteration for all roots.

Initialize p, q, r, s:

p₀ := (0.4 + 0.9 i)⁰ ;

q₀ := (0.4 + 0.9 i)¹ ;

r₀ := (0.4 + 0.9 i)² ;

s₀ := (0.4 + 0.9 i)³ ;

There is nothing special about choosing 0.4 + 0.9 i except that it is neither a real number nor a root of unity.

Make the substitutions for n = 1,2,3,···

$p_n = p_{n-1} - \frac{f(p_{n-1})}{ (p_{n-1}-q_{n-1})(p_{n-1}-r_{n-1})(p_{n-1}-s_{n-1}) };$

$q_n = q_{n-1} - \frac{f(q_{n-1})}{ (q_{n-1}-p_n)(q_{n-1}-r_{n-1})(q_{n-1}-s_{n-1}) };$

$r_n = r_{n-1} - \frac{f(r_{n-1})}{ (r_{n-1}-p_n)(r_{n-1}-q_n)(r_{n-1}-s_{n-1}) };$

$s_n = s_{n-1} - \frac{f(s_{n-1})}{ (s_{n-1}-p_n)(s_{n-1}-q_n)(s_{n-1}-r_n) }.$

Re-iterate until the numbers p, q, r, s stop essentially changing. Then they have the values P, Q, R, S in some order and in the chosen precision. So the problem is solved.

Note that you must use complex number arithmetic, and that the roots are found simultaneously rather than one at a time.

Variations

This iteration procedure, like the Gauss–Seidel method for linear equations, computes one number at a time based on the already computed numbers. A variant of this procedure, like the Jacobi method, computes a vector of root approximations at a time. Both variant are effective root-finding algorithms.

One could also choose the initial values for p,q,r,s by some other procedure, even randomly, but in a way that

• they are inside some not too large circle containing also the roots of ƒ(x), e.g. the circle around the origin with radius 1 + max( | a | , | b | , | c | , | d | ), (where 1,a,b,c,d are the coefficients of ƒ(x)) and that
• they are not too close to each other, which may increasingly become a concern as the degree of the polynomial increases.

Example

This example is from the reference 1992. The equation solved is x³ − 3x² +3x − 5 = 0 . The first 4 iterations move p, q, r seemingly chaotically, but then the roots are located to 1 decimal. After iteration number 5 we have 4 correct decimals, and the subsequent iteration number 6 confirms that the computed roots are fixed. This general behaviour is characteristic for the method.

it.-no.	p	q	r
0	+1.0000+0.0000i	+0.4000+0.9000i	−0.6500+0.7200i
1	+1.3608+2.0222i	−0.3658+2.4838i	−2.3858−0.0284i
2	+2.6597+2.7137i	+0.5977+0.8225i	−0.6320−1.6716i
3	+2.2704+0.3880i	+0.1312+1.3128i	+0.2821−1.5015i
4	+2.5428−0.0153i	+0.2044+1.3716i	+0.2056−1.3721i
5	+2.5874+0.0000i	+0.2063+1.3747i	+0.2063−1.3747i
6	+2.5874+0.0000i	+0.2063+1.3747i	+0.2063−1.3747i

Note that the equation has one real root and one pair of complex conjugate roots, and that the sum of the roots is 3.

Derivation of the method via Newton's method

For every n-tuple of complex numbers, there is exactly one monic polynomial of degree n that has them as its zeros (keeping multiplicities). This polynomial is given by multiplying all the corresponding linear factors, that is

$g_{\vec z}(X)=(X-z_1)\cdots(X-z_n).$

This polynomial has coefficients that depend on the prescribed zeros,

$g_{\vec z}(X)=X^n+g_{n-1}(\vec z)X^{n-1}+\cdots+g_0(\vec z).$

Those coefficients are, up to a sign, the elementary symmetric polynomials $\alpha_1(\vec z),\dots,\alpha_n(\vec z)$ of degrees 1,...,n.

To find all the roots of a given polynomial $f(X)=X^n+c_{n-1}X^{n-1}+\cdots+c_0$ with coefficient vector $(c_{n-1},\dots,c_0)$ simultaneously is now the same as to find a solution vector to the system

$\begin{matrix} c_0&=&g_0(\vec z)&=&(-1)^n\alpha_n(\vec z)&=&(-1)^nz_1\cdots z_n\\ c_1&=&g_1(\vec z)&=&(-1)^{n-1}\alpha_{n-1}(\vec z)\\ &\vdots&\\ c_{n-1}&=&g_{n-1}(\vec z)&=&-\alpha_1(\vec z)&=&-(z_1+z_2+\cdots+z_n). \end{matrix}$

The Durand–Kerner method is obtained as the multidimensional Newton's method applied to this system. It is algebraically more comfortable to treat those identities of coefficients as the identity of the corresponding polynomials, $g_{\vec z}(X)=f(X)$. In the Newton's method one looks, given some initial vector $\vec z$, for an increment vector $\vec w$ such that $g_{\vec z+\vec w}(X)=f(X)$ is satisfied up to second and higher order terms in the increment. For this one solves the identity

$f(X)-g_{\vec z}(X)=\sum_{k=1}^n\frac{\partial g_{\vec z}(X)}{\partial z_k}w_k=-\sum_{k=1}^n w_k\prod_{j\ne k}(X-z_j).$

If the numbers $z_1,\dots,z_n$ are pairwise different, then the polynomials in the terms of the right hand side form a basis of the n-dimensional space $\mathbb C[X]_{n-1}$ of polynomials with maximal degree n − 1. Thus a solution $\vec w$ to the increment equation exists in this case. The coordinates of the increment $\vec w$ are simply obtained by evaluating the increment equation

$-\sum_{k=1}^n w_k\prod_{j\ne k}(X-z_j)=f(X)-\prod_{j=1}^n(X-z_j)$

at the points X = z_k, which results in

$-w_k\prod_{j\ne k}(z_k-z_j)=-w_kg_{\vec z}'(z_k)=f(z_k)$, that is $w_k=-\frac{f(z_k)}{\prod_{j\ne k}(z_k-z_j)}.$

Root inclusion via Gerschgorin's circles

In the quotient ring (algebra) of residue classes modulo ƒ(X), the multiplication by X defines an endomorphism that has the zeros of ƒ(X) as eigenvalues with the corresponding multiplicities. Choosing a basis, the multiplication operator is represented by its coefficient matrix A, the companion matrix of ƒ(X) for this basis.

Since every polynomial can be reduced modulo ƒ(X) to a polynomial of degree n − 1 or lower, the space of residue classes can be identified with the space of polynomials of degree bounded by n − 1. A problem specific basis can be taken from Lagrange interpolation as the set of n polynomials

$b_k(X)=\prod_{1\le j\le n,\;j\ne k}(X-z_j),\quad k=1,\dots,n,$

where $z_1,\dots,z_n\in\mathbb C$ are pairwise different complex numbers. Note that the kernel functions for the Lagrange interpolation are $L_k(X)=\frac{b_k(X)}{b_k(z_k)}$.

For the multiplication operator applied to the basis polynomials one obtains from the Lagrange interpolation

$X\cdot b_k(X)\mod f(X)=X\cdot b_k(X)-f(X)$	$=\sum_{j=1}^n\Big(z_j\cdot b_k(z_j)-f(z_j)\Big)\cdot \frac{b_j(X)}{b_j(z_j)}$
	$=z_k\cdot b_k(X)+\sum_{j=1}^n w_j\cdot b_j(X)$,

where $w_j=-\frac{f(z_j)}{b_j(z_j)}$ are again the Weierstrass updates.

The companion matrix of ƒ(X) is therefore

$A = \mathrm{diag}(z_1,\dots,z_n) +\begin{pmatrix}1\\\vdots\\1\end{pmatrix}\cdot\left(w_1,\dots,w_n\right).$

From the transposed matrix case of the Gershgorin circle theorem it follows that all eigenvalues of A, that is, all roots of ƒ(X), are contained in the union of the disks D(a_k,k,r_k) with a radius $r_k=\sum_{j\ne k}\big|a_{j,k}\big|$.

Here one has a_k,k = z_k + w_k, so the centers are the next iterates of the Weierstrass iteration, and radii $r_k=(n-1)\left|w_k\right|$ that are multiples of the Weierstrass updates. If the roots of ƒ(X) are all well isolated (relative to the computational precision) and the points $z_1,\dots,z_n\in\mathbb C$ are sufficidently close approximations to these roots, then all the disks will become disjoint, so each one contains exactly one zero. The midpoints of the circles will be better approximations of the zeros.

Every conjugate matrix TAT ^{− 1} of A is as well a companion matrix of ƒ(X). Choosing T as diagonal matrix leaves the structure of A invariant. The root close to z_k is contained in any isolated circle with center z_k regardless of T. Choosing the optimal diagonal matrix T for every index results in better estimates (see ref. Petkovic et al. 1995).

Convergence results

The connection between the Taylor series expansion and Newton's method suggests that the distance from z_k + w_k to the corresponding root is of the order O( | w_k | ²), if the root is well isolated from nearby roots and the approximation is sufficiently close to the root. So after the approximation is close, Newton's method converges quadratically; that is: the error is squared with every step (which will greatly reduce the error once it is less than 1). In the case of the Durand–Kerner method, convergence is quadratic if the vector $\vec z=(z_1,\dots,z_n)$ is close to some permutation of the vector of the roots of ƒ.

For the conclusion of linear convergence there is a more specific result (see ref. Petkovic et al. 1995). If the initial vector $\vec z$ and its vector of Weierstrass updates $\vec w=(w_1,\dots,w_n)$ satisfies the inequality

$\max_{1\le k\le n}\big|w_k\big| \le \frac1{5n} \min_{1\le j<k\le n}\big|z_k-z_j\big|,$

then this inequality also holds for all iterates, all inclusion disks $\textstyle D\left(z_k+w_k,(n-1)|w_k|\right)$ are disjoint and linear convergence with a contraction factor of 1/2 holds. Further, the inclusion disks can in this case be chosen as

$\textstyle D\left(z_k+w_k,\frac14 |w_k|\right)\qquad k = 1,\dots, n,$

each containing exactly one zero of ƒ.

References

• Weierstraß, Karl (1891). "Neuer Beweis des Satzes, dass jede ganze rationale Function einer Veränderlichen dargestellt werden kann als ein Product aus linearen Functionen derselben Veränderlichen". Sitzungsberichte der königlich preussischen Akademie der Wissenschaften zu Berlin. 
• Durand, E. (1960). "Equations du type F(x) = 0: Racines d'un polynome". In Masson et al.. Solutions numeriques des equations algebriques, vol. 1. 
• Kerner, Immo O. (1966). "Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen". Numerische Mathematik 8: 290–294. doi:10.1007/BF02162564. http://www.springerlink.com/content/q5p055l61pm63206. 
• Prešić, Marica (1980). "A convergence theorem for a method for simultaneous determination of all zeros of a polynomial". Publications de l'institut mathematique (Beograd) (N.S.) 28 (42): 158–168. 
• Petkovic, M.S., Carstensen, C. and Trajkovic, M. (1995). "Weierstrass formula and zero-finding methods". Numerische Mathematik 69: 353–372. doi:10.1007/s002110050097. http://www.springerlink.com/content/x467nejrv3c8hq8j. 
• Bo Jacoby, Nulpunkter for polynomier, 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 Danish Society for Structural Science and Engineering) volume 63 no. 3-4, 1992, pp. 83–105.
• Gourdon, Xavier (1996). Combinatoire, Algorithmique et Geometrie des Polynomes. Paris: Ecole Polytechnique. http://algo.inria.fr/gourdon/thesis.html. 
• Victor Pan (May 2002): Univariate Polynomial Root-Finding with Lower Computational Precision and Higher Convergence Rates. Tech-Report, City University of New York
• Neumaier, Arnold (2003). "Enclosing clusters of zeros of polynomials". Journal of Computational and Applied Mathematics 156: 389. doi:10.1016/S0377-0427(03)00380-7. http://www.mat.univie.ac.at/~neum/papers.html#polzer. 
• Jan Verschelde, The method of Weierstrass (also known as the Durand-Kerner method), 2003.

External links

• Ada Generic_Roots using the Durand-Kerner Method — an open-source implementation in Ada

• Polynomial Roots — an open-source implementation in Java

• Roots Extraction from Polynomials : The Durand-Kerner Method — contains a Java applet demonstration