Laguerre's method

In numerical analysis, Laguerre's method is a root-finding algorithm tailored to polynomials. In other words, Laguerre's method can be used to numerically solve the equation p(x) = 0 for a given polynomial p(x). One of the most useful properties of this method is that it is, from extensive empirical study, very close to being a "sure-fire" method, meaning that it is almost guaranteed to always converge to some root of the polynomial, no matter what initial guess is chosen. However, for computer computation, more efficient methods are known, with which it is guaranteed to find all roots (see Root-finding algorithm § Roots of polynomials) or all real roots (see Real-root isolation).

This method is named in honour of Edmond Laguerre, a French mathematician.

Definition

The algorithm of the Laguerre method to find one root of a polynomial p(x) of degree n is:

• Choose an initial guess x₀
• For k = 0, 1, 2, …
  • If ${\displaystyle p(x_{k})}$ is very small, exit the loop
  • Calculate ${\displaystyle G={\frac {p'(x_{k})}{p(x_{k})}}}$
  • Calculate ${\displaystyle H=G^{2}-{\frac {p''(x_{k})}{p(x_{k})}}}$
  • Calculate ${\displaystyle a={\frac {n}{G\pm {\sqrt {(n-1)(nH-G^{2})}}}}}$, where the sign is chosen to give the denominator with the larger absolute value, to avoid catastrophic cancellation.
  • Set ${\displaystyle x_{k+1}=x_{k}-a}$
• Repeat until a is small enough or if the maximum number of iterations has been reached.

If a root has been found, the corresponding linear factor can be removed from p. This deflation step reduces the degree of the polynomial by one, so that eventually, approximations for all roots of p can be found. Note however that deflation can lead to approximate factors that differ significantly from the corresponding exact factors. This error is least if the roots are found in the order of increasing magnitude.

Derivation

The fundamental theorem of algebra states that every nth degree polynomial ${\displaystyle p}$ can be written in the form

${\displaystyle p(x)=C(x-x_{1})(x-x_{2})\cdots (x-x_{n}),}$

such that ${\displaystyle x_{k}}$ where ${\displaystyle (k=1,2,...,n)}$ are the roots of the polynomial. If we take the natural logarithm of both sides, we find that

${\displaystyle \ln |p(x)|=\ln |C|+\ln |x-x_{1}|+\ln |x-x_{2}|+\cdots +\ln |x-x_{n}|.}$

Denote the derivative by

${\displaystyle G={\frac {d}{dx}}\ln |p(x)|={\frac {1}{x-x_{1}}}+{\frac {1}{x-x_{2}}}+\cdots +{\frac {1}{x-x_{n}}},}$

and the negated second derivative by

${\displaystyle H=-{\frac {d^{2}}{dx^{2}}}\ln |p(x)|={\frac {1}{(x-x_{1})^{2}}}+{\frac {1}{(x-x_{2})^{2}}}+\cdots +{\frac {1}{(x-x_{n})^{2}}}.}$

We then make what Acton calls a 'drastic set of assumptions', that the root we are looking for, say, ${\displaystyle x_{1}}$ is a certain distance away from our guess ${\displaystyle x}$, and all the other roots are clustered together some distance away. If we denote these distances by

${\displaystyle a=x-x_{1}\,}$

and

${\displaystyle b=x-x_{i},\quad i=2,3,\ldots ,n}$

then our equation for ${\displaystyle G}$ may be written as

${\displaystyle G={\frac {1}{a}}+{\frac {n-1}{b}}}$

and the expression for ${\displaystyle H}$ becomes

${\displaystyle H={\frac {1}{a^{2}}}+{\frac {n-1}{b^{2}}}.}$

Solving these equations for ${\displaystyle a}$, we find that

${\displaystyle a={\frac {n}{G\pm {\sqrt {(n-1)(nH-G^{2})}}}}}$,

where the square root of a complex number is chosen to produce larger absolute value of the denominator, or equivalently, to satisfy:

${\displaystyle \operatorname {Re} \,({\overline {G}}{\sqrt {(n-1)(nH-G^{2})}}\,)>0}$,

where Re denotes real part of a complex number, and G is the complex conjugate of G; or

${\displaystyle a={\frac {p(x)}{p'(x)}}\cdot \left({\frac {1}{n}}+{\frac {n-1}{n}}\,{\sqrt {1-{\frac {n}{n-1}}\,{\frac {p(x)p''(x)}{p'(x)^{2}}}}}\right)^{-1}}$,

where the square root of a complex number is chosen to have a non-negative real part.

For small values of p(x) this formula differs from the offset of the third order Halley's method by an error of O(p(x)³), so convergence close to a root will be cubic as well.

Note that, even if the 'drastic set of assumptions' does not work for some particular polynomial p(x), p(x) can be transformed into a related polynomial r for which the assumptions are correct, e.g. by shifting the origin towards a suitable complex number w, giving q(x) = p(x − w), to give distinct roots distinct magnitudes if necessary (which it will be if some roots are complex conjugates), and then getting r from q(x) by repeatedly applying the root squaring transformation used in Graeffe's method enough times to make the smaller roots significantly smaller than the largest root (and so, clustered in comparison); the approximate root from Graeffe's method can then be used to start the new iteration for Laguerre's method on r. An approximate root for p(x) may then be obtained straightforwardly from that for r.

If we make the stronger assumption that the terms in G corresponding to the roots x_i, i = 2, 3, …, n are negligibly small in comparison to the term corresponding to the root x₁, this leads to Newton's method.

Properties

Attraction zones of Laguerre's method for polynomial x^4 + 2*x^3 + 3*x^2 + 4*x + 1

If x is a simple root of the polynomial p(x), then Laguerre's method converges cubically whenever the initial guess x₀ is close enough to the root x. On the other hand, if x is a multiple root then the convergence is only linear. This is obtained with the penalty of calculating values for the polynomial and its first and second derivatives at each stage of the iteration.

A major advantage of Laguerre's method is that it is almost guaranteed to converge to some root of the polynomial no matter where the initial approximation is chosen. This is in contrast to other methods such as the Newton–Raphson method which may fail to converge for poorly chosen initial guesses. It may even converge to a complex root of the polynomial, because of the square root being taken in the calculation of a above may be of a negative number. This may be considered an advantage or a liability depending on the application to which the method is being used. Empirical evidence has shown that convergence failure is extremely rare, making this a good candidate for a general purpose polynomial root finding algorithm. However, given the fairly limited theoretical understanding of the algorithm, many numerical analysts are hesitant to use it as such, and prefer better understood methods such as the Jenkins–Traub algorithm, for which more solid theory has been developed. Nevertheless, the algorithm is fairly simple to use compared to these other "sure-fire" methods, easy enough to be used by hand or with the aid of a pocket calculator when an automatic computer is unavailable. The speed at which the method converges means that one is only very rarely required to compute more than a few iterations to get high accuracy.

References

• Acton, Forman S. (1970). Numerical Methods that Work . Harper & Row. ISBN   0-88385-450-3.
• Goedecker, S. (1994). "Remark on Algorithms to Find Roots of Polynomials". SIAM J. Sci. Comput. 15 (5): 1059–1063. doi:10.1137/0915064.
• Mekwi, Wankere R. (2001). "Iterative Methods for Roots of Polynomials". Master's Thesis, University of Oxford. Archived from the original on 2012-12-23.
• Pan, V. Y. (1997). "Solving a Polynomial Equation: Some History and Recent Progress". SIAM Rev. 39 (2): 187–220. doi:10.1137/S0036144595288554.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 9.5.3. Laguerre's Method". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN   978-0-521-88068-8.
• Ralston, Anthony; Rabinowitz, Philip (1978). A First Course in Numerical Analysis. McGraw-Hill. ISBN   0-07-051158-6.