Durand–Kerner method: Difference between revisions

Revision as of 22:10, 23 April 2019

Latest revision as of 05:14, 4 August 2025

(38 intermediate revisions by 29 users not shown)

Line 1:

{{Short description|Root-finding algorithm for polynomials}}

In [[numerical analysis]], the '''Durand–Kerner method''', discovered by [[Karl Weierstrass]] in 1891 and rediscovered independently by Durand in 1960 and Kerner in 1966, is a [[root-finding algorithm]] for solving [[polynomial]] [[equation (mathematics)|equations]].<ref name="Petković">{{cite book|last1=Petković|first1=Miodrag|title=Iterative methods for simultaneous inclusion of polynomial zeros|date=1989|publisher=Springer|___location=Berlin [u.a.]|isbn=978-3-540-51485-5|pages=31–32}}</ref> In other words, the method can be used to solve numerically the equation

{{refstyle|date=November 2020}}

In [[numerical analysis]], the '''Weierstrass method''' or '''Durand–Kerner method''', discovered by [[Karl Weierstrass]] in 1891 and rediscovered independently by Durand in 1960 and Kerner in 1966, is a [[root-finding algorithm]] for solving [[polynomial]] [[equation (mathematics)|equations]].<ref name="Petković">{{cite book |last1=Petković |first1=Miodrag |title=Iterative methods for simultaneous inclusion of polynomial zeros |date=1989 |publisher=Springer |___location=Berlin [u.a.] |isbn=978-3-540-51485-5 |pages=31–32}}</ref> In other words, the method can be used to solve numerically the equation ''f''(''x''){{=}}0, where ''f'' is a given polynomial, which can be taken to be scaled so that the leading coefficient is&nbsp;1.

: ƒ(''x'') = 0

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

==Explanation==

TheThis explanation is forconsiders equations of [[Degree of a polynomial|degree]] four. It is easily generalized to other degrees.

Let the polynomial &fnof;''f'' be defined by

:&fnof; <math>f(''x'') = ''x''^{^4} + ''ax''^{^3} + ''bx''^{^2} + ''cx'' + ''d''</math>

for all ''x''. The known numbers ''a'', ''b'', ''c'', ''d'' are the [[coefficient]]s.

for all ''x''.

TheLet knownthe (potentially complex) numbers ''aP'', b''Q'', c''R'', d''S'' arebe the [[coefficient]]sroots of this polynomial ''f''. Then

: <math>x_1:=x_0-\frac{f(x_0x) = (x - P)}{(x_0x - Q)(x_0x - R)(x_0x - S)},</math>

Let the (complex) numbers ''P,Q,R,S'' be the roots of this polynomial &fnof;.

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

Then

:&fnof;(''x'') <math>P = (''x'' &minus;- ''P''\frac{f(x)}{(''x'' &minus;- ''Q'')(''x'' &minus;- ''R'')(''x'' &minus;- ''S'')}.</math>

So if used as a [[fixed point (mathematics)|fixed -point]] [[iteration]]

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

: <math>x_{k+1}x_1 :=x_k x_0 - \frac{f(x_kx_0)}{(x_kx_0 -q Q)(x_kx_0 -r R)(x_kx_0 -s S)},</math>

it is strongly stable in that every initial point ''x''₀ ≠ ''Q'', ''R'', ''S'' delivers after one iteration the root ''P'' = ''x''₁. Furthermore, 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

: <math>Px_{k+1} :=x x_k - \frac{f(xx_k)}{(xx_k -Q q)(xx_k -R r)(xx_k -S s)}.,</math>

So if used as a [[fixed point (mathematics)|fixed point]] [[iteration]]

:<math>x_1:=x_0-\frac{f(x_0)}{(x_0-Q)(x_0-R)(x_0-S)},</math>

it is strongly stable in that every initial point ''x₀'' ≠ ''Q,R,S''

delivers after one iteration the root ''P=x₁''.

Furthermore, 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

:<math>x_{k+1}:=x_k-\frac{f(x_k)}{(x_k-q)(x_k-r)(x_k-s)},</math>

since

: <math>P - \frac{f(P)}{(P - q)(P - r)(P - s)} = P - 0 = P.</math>

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

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.

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&nbsp;''i'')⁰ ;,

: ''q''₀ := (0.4 + 0.9&nbsp;''i'')¹ ;,

: ''r''₀ := (0.4 + 0.9&nbsp;''i'')² ;,

: ''s''₀ := (0.4 + 0.9&nbsp;''i'')³ ;.

There is nothing special about choosing 0.4&nbsp;+&nbsp;0.9&nbsp;''i'' except that it is neither a [[real number]] nor a [[root of unity]].

Make the substitutions for ''n'' = 1, 2, 3,&middot;&middot;&middot; ...:

|: <math> q_np_n = q_p_{n-1} - \frac{f(q_p_{n-1})}{ (q_p_{n-1} -p_n q_{n-1})(q_p_{n-1} - r_{n-1})(q_p_{n-1} - s_{n-1}) }; ,</math>

:{|

|: <math> r_nq_n = r_q_{n-1} - \frac{f(r_q_{n-1})}{ (r_q_{n-1} - p_n)(r_q_{n-1} -q_n r_{n-1})(r_q_{n-1} - s_{n-1}) }; ,</math>

|-

|: <math> p_nr_n = p_r_{n-1} - \frac{f(p_r_{n-1})}{ (p_r_{n-1} -q_{n-1} p_n)(p_r_{n-1} -r_{n-1} q_n)(p_r_{n-1} - s_{n-1}) }; ,</math>

|: <math> 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) }. </math>

|-

|<math> 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}) }; </math>

|-

|<math> 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}) }; </math>

|-

|<math> 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) }. </math>

|}

Re-iterate until the numbers ''p'', ''q'', ''r'', ''s'' essentially stop changing relative to the desired precision. They then have the values ''P'', ''Q'', ''R'', ''S'' in some order and in the chosen precision. So the problem is solved.

Re-iterate until the numbers ''p, q, r, s''

stop essentially changing in relative to the desired precision.

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 must be used, and that the roots are found simultaneously rather than one at a time.

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 variants are effective root-finding algorithms.

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 &fnof;''f''(''x''), e.g. the circle around the origin with radius <math>1 + \max\big(|a|, |b|, |c|, |d|\big)</math>, (where 1, ''a'', ''b'', ''c'', ''d'' are the coefficients of &fnof;''f''(''x''))

by some other procedure, even randomly, but in a way that

*they are inside some not-too-large circle containing also the roots of &fnof;(''x''), e.g. the circle around the origin with radius <math>1+\max(|a|,|b|,|c|,|d|)</math>, (where 1,''a,b,c,d'' are the coefficients of &fnof;(''x''))

and that

* they are not too close to each other,

which may increasingly become a concern as the degree of the polynomial increases.

as the degree of the polynomial increases.

If the coefficients are real and the polynomial has odd degree, then it must have at least one real root. To find this, use a real value of ''p''₀ as the initial guess and make ''q''₀ and ''r''₀, etc., [[complex conjugate]] pairs. Then the iteration will preserve these properties; that is, ''p''_''n'' will always be real, and ''q''_''n'' and ''r''_''n'', etc., will always be conjugate. In this way, the ''p''_''n'' will converge to a real root ''P''. Alternatively, make all of the initial guesses real; they will remain so.

== Example ==

This example is from the reference Jacoby (1992). The equation solved is {{nowrap|1=''x''³ − 3''x''² + 3''x'' − 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. Also notice that, in this example, the roots are used as soon as they are computed in each iteration. In other words, the computation of each second column uses the value of the previous computed columns.

::{|class="wikitable"

Line 157 ⟶ 125:

Those coefficients are, up to a sign, the [[elementary symmetric polynomial]]s <math>\alpha_1(\vec z),\dots,\alpha_n(\vec z)</math> of degrees ''1,...,n''.

To find all the roots of a given polynomial <math>f(X)=X^n+c_{n-1}X^{n-1}+\cdots+c_0</math> with coefficient vector <math>(c_{n-1},\dots,c_0)</math> simultaneously is now the same as to find a solution vector to the [[Vieta's_formulas|Vieta's system]]

:<math>\begin{matrix}

Line 187 ⟶ 155:

In the [[quotient ring]] (algebra) of [[residue class]]es modulo &fnof;(''X''), the multiplication by ''X'' defines an [[endomorphism]] that has the zeros of &fnof;(''X'') as [[eigenvalue]]s with the corresponding multiplicities. Choosing a basis, the multiplication operator is represented by its coefficient matrix ''A'', the [[companion matrix]] of &fnof;(''X'') for this basis.

Since every polynomial can be reduced modulo &fnof;(''X'') to a polynomial of degree ''n''&nbsp;&minus;&nbsp;1 or lower, the space of residue classes can be identified with the space of polynomials of degree bounded by ''n''&nbsp;&minus;&nbsp;1. A problem-specific basis can be taken from [[Lagrange interpolation]] as the set of ''n'' polynomials

A problem specific basis can be taken from [[Lagrange interpolation]] as the set of ''n'' polynomials

:<math>b_k(X)=\prod_{1\le j\le n,\;j\ne k}(X-z_j),\quad k=1,\dots,n,</math>

Line 212 ⟶ 179:

From the transposed matrix case of the [[Gershgorin circle theorem]] it follows that all eigenvalues of ''A'', that is, all roots of &fnof;(''X''), are contained in the union of the disks <math>D(a_{k,k},r_k)</math> with a radius <math>r_k=\sum_{j\ne k}\big|a_{j,k}\big|</math>.

Here one has <math>a_{k,k}=z_k+w_k</math>, so the centers are the next iterates of the Weierstrass iteration, and radii <math>r_k=(n-1)\left|w_k\right|</math> that are multiples of the Weierstrass updates. If the roots of &fnof;(''X'') are all well isolated (relative to the computational precision) and the points <math>z_1,\dots,z_n\in\mathbb C</math> are sufficidentlysufficiently 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 <math>TAT^{-1}</math> of ''A'' is as well a companion matrix of &fnof;(''X''). Choosing ''T'' as diagonal matrix leaves the structure of ''A'' invariant. The root close to <math>z_k</math> is contained in any isolated circle with center <math>z_k</math> regardless of ''T''. Choosing the optimal diagonal matrix ''T'' for every index results in better estimates (see ref. Petkovic et al. 1995).

Line 218 ⟶ 185:

==Convergence results==

The connection between the Taylor series expansion and Newton's method suggests that the distance from <math>z_k + w_k</math> to the corresponding root is of the order <math>O\big(|w_k|^2\big)</math>, 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 <math>\vec z = (z_1, \dots, z_n)</math> is close to some permutation of the vector of the roots of &fnof;''f''.

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

: <math>\max_{1 \le k \le n}\big |w_k\big| \le \frac1frac{1}{5n} \min_{1 \le j < k \le n}\big |z_k - z_j\big|,</math>

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

are disjoint and linear convergence with a contraction factor of ''1/2'' holds. Further, the inclusion disks can in this case be chosen as

: <math>\textstyle D\left(z_k + w_k, \frac14tfrac14 |w_k|\right),\qquadquad k = 1, \dots, n,</math>

each containing exactly one zero of &fnof;''f''.

==Failing general convergence==

The Weierstrass / Durand-Kerner method is not generally convergent: in other words, it is not true that for every polynomial, the set of initial vectors that eventually converges to roots is open and dense. In fact, there are open sets of polynomials that have open sets of initial vectors that converge to periodic cycles other than roots (see Reinke et al.)

==References==

{{Reflist}}

* {{cite conference|last=Weierstraß|first= Karl|authorlinkauthor-link=Karl Weierstraß|title=Neuer Beweis des Satzes, dass jede ganze rationale Function einer Veränderlichen dargestellt werden kann als ein Product aus linearen Functionen derselben Veränderlichen|booktitlebook-title=Sitzungsberichte der königlich preussischen Akademie der Wissenschaften zu Berlin|year=1891|url=http://bibliothek.bbaw.de/bibliothek-digital/digitalequellen/schriften/anzeige?band=10-sitz/1891-2&seite:int=00000565|access-date=2013-10-31|archive-date=2013-11-02|archive-url=https://web.archive.org/web/20131102093616/http://bibliothek.bbaw.de/bibliothek-digital/digitalequellen/schriften/anzeige?band=10-sitz%2F1891-2&seite%3Aint=00000565|url-status=dead}}

* {{cite conference|last=Durand|first=E.|booktitlebook-title=Solutions Numériques des Equations Algébriques |volume=1|editor=Masson |display-editors=etal |title=Equations du type ''F''(''x'')&nbsp;=&nbsp;0: Racines d'un polynome|year= 1960}}

* {{cite journal|last= Kerner|first= Immo O.|title=Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen|journal=Numerische Mathematik|volume=8|year= 1966|pagesissue= 290–2943|urlpages= http://www.springerlink.com/content/q5p055l61pm63206290–294|doi=10.1007/BF02162564 |s2cid= 115307022}}

* {{cite journal|author=Prešić, Marica|title=A convergence theorem for a method for simultaneous determination of all zeros of a polynomial |url=http://elib.mi.sanu.ac.rs/files/journals/publ/48/n042p159.pdf|journal=Publications de l'institutInstitut mathematiqueMathématique (Beograd)|series=Nouvelle (N.S.)Série|volume=28|pages=158–168 |year=1980 | issue=42}}

* {{cite journal|author=Petkovic, M.S., Carstensen, C. and Trajkovic, M.|title=Weierstrass formula and zero-finding methods|journal=Numerische Mathematik|volume=69|year=1995|issue=3|pages=353–372|urldoi=10.1007/s002110050097|citeseerx=10.1.1.53.7516|s2cid=18594004 }}

* 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-43–4, 1992, pp.&nbsp;83–105.

* {{cite book|last=Gourdon|first=Xavier|title=Combinatoire, Algorithmique et Geometrie des Polynomes|publisher=EcoleÉcole Polytechnique|___location= Paris|year=1996|url=http://algo.inria.fr/gourdon/thesis.html|access-date=2006-08-22|archive-url=https://web.archive.org/web/20061028060132/http://algo.inria.fr/gourdon/thesis.html|archive-date=2006-10-28|url-status=dead}}

* [[Victor Pan]] (May 2002): [https://web.archive.org/web/20060907205721/http://www.cs.gc.cuny.edu/tr/techreport.php?id=26 ''Univariate Polynomial Root-Finding with Lower Computational Precision and Higher Convergence Rates'']. Tech-Report, City University of New York

* {{cite journal|first= Arnold|last= Neumaier|title= Enclosing clusters of zeros of polynomials|journal= Journal of Computational and Applied Mathematics|volume= 156 |year=2003|issue= 2|url=httphttps://www.mat.univie.ac.at/~neum/papers.html#polzer|doi= 10.1016/S0377-0427(03)00380-7|pages= 389389–401|bibcode= 2003JCoAM.156..389N|doi-access= free}}

* Jan Verschelde, ''[httphttps://www2.math.uic.edu/~jan/mcs471f03/Project_Two/proj2/node2.html The method of Weierstrass (also known as the Durand–Kerner method)]'', 2003.

* Bernhard Reinke, Dierk Schleicher, and Michael Stoll, "[https://arxiv.org/abs/2004.04777 The Weierstrass root finder is not generally convergent]", 2020

** Bernhard Reinke, Dierk Schleicher and Michael Stoll: "The Weierstrass-Durand-Kerner root finder is not generally convergent", Math. Comp. vol.92 (2023), pp.839-866. DOI: https://doi.org/10.1090/mcom/3783 .

==External links==

* ''[httphttps://homesites.roadrunnergoogle.com/~jbmatthewssite/miscdrjohnbmatthews/groots.html?authuser=0 Ada Generic_Roots using the Durand–Kerner Method] [https://web.archive.org/web/20090701041632/http://home.roadrunner.com/~jbmatthews/misc/groots.html (archive)]'' &mdash; an [[open-source]] implementation in [[Ada programming language|Ada]]

* ''[http://sites.google.com/site/drjohnbmatthews/polyroots Polynomial Roots]'' &mdash; an [[Open-Source|open-source]] implementation in [[Java programming language|Java]]

* ''[https://web.archive.org/web/20060505043302/http://www.cpc.wmin.ac.uk/~spiesf/Solve/solve.html Roots Extraction from Polynomials : The Durand–Kerner Method]'' &mdash; contains a [[Java applet]] demonstration

Both variant are effective {{root-finding algorithms.}}

[[Category:Root-findingPolynomial factorization algorithms]]

{{DEFAULTSORT:Durand-Kerner method}}

[[Category:Root-finding algorithms]]