Root-finding algorithm
From Freepedia
A root-finding algorithm is a numerical method or algorithm for finding a value x such that f(x) = 0, for a given function f. Such an x is called a root of the function f.
Finding a root of f − g is the same as solving the equation f(x) = g(x). Here, x is called the unknown in the equation. Conversely, any equation can take the canonical form f(x) = 0, so equation solving is the same thing as computing (or finding) a root of a function.
All numerical root-finding methods use iteration, producing a sequence of numbers (hopefully) converging towards a limit which is a root. The first values of this series are initial guesses. The method computes subsequent values based on the old ones and the function f.
The behaviour of root-finding algorithms is studied in numerical analysis.
Contents |
Specific algorithms
The simplest root-finding algorithm is the bisection method. It works when f is a continuous function and it requires previous knowledge of two initial guesses, a and b, such that f(a) and f(b) have opposite signs.
Newton's method, assumes the function f to have a derivative. Newton's method may not converge if you start too far away from a root. However, if it does converge, it is faster than the bisection method (convergence is quadratic). Newton's method is also important because it readily generalizes to higher-dimensional problems.
Replacing the derivative in Newton's method with a finite difference, we get the secant method. This method does not require the computation (nor the existence) of a derivative, but the price is slower convergence (the order is approximately 1.6).
The false position method, also called the regula falsi method, is like the bisection method. However, it does not cut the interval in two equal parts at every iteration, but it cuts the interval at the point given by the formula for the secant method. The false position method inherits the robustness of the bisection method and the superlinear convergence of the secant method.
The secant method also arises if one approximates the unknown function f by linear interpolation. When quadratic interpolation is used instead, one arrives at Muller's method. It converges faster than the secant method. A particular feature of this method is that the iterates xn may become complex.
This can be avoided by interpolating the inverse of f, resulting in the inverse quadratic interpolation method. Again, convergence is asymptotically faster than the secant method, but inverse quadratic interpolation often behaves poorly when the iterates are not close to the root.
Finally, Brent's method is a combination of the bisection method, the secant method and inverse quadratic interpolation. At every iteration, Brent's method decides which method out of these three is likely to do best, and proceeds by doing a step according to that method. This gives a robust and fast method, which therefore enjoys considerable popularity.
Finding roots of polynomials
Much attention has been given to the special case that the function f is in fact a polynomial. Of course, the method described in the previous section can be used. In particular, it is easy to find the derivative of a polynomial, so Newton's method is a viable candidate. But one can also choose a method that exploits the fact that f is a polynomial.
One possibility is to form the companion matrix of the polynomial. Since the eigenvalues of this matrix coincide with the roots of the polynomial, one can now use any eigenvalue algorithm to find the roots of the original polynomial.
Another possibility is Laguerre's method, which is a rather complicated method that converges very fast. In fact, it exhibits cubic convergence for simple roots, beating the quadratic convergence displayed by Newton's method.
If the polynomial has rational coefficients and only rational roots are of interest, then Ruffini's method can also be used.
In any case, it should be borne in mind that the problem of finding roots can be ill-conditioned, as the example of Wilkinson's polynomial shows.
How to solve an algebraic equation
A method is explained here for equations of degree four in order to simplify the notation. It is easily generalized to equations of other degrees.
An algebraic equation of degree four is of the form <math>\ f(x)=0 </math>, where <math>\ f</math> is a monic polynomial of degree four
- <math>\ f(x)=x^4+ax^3+bx^2+cx+d=(((x+a)x+b)x+c)x+d </math> for all <math>\ x</math>.
The known numbers <math>\ a,b,c,d</math> are the coefficients.
The product of four monic polynomials of degree one is a monic polynomial of degree four
- <math>\ f(x)=(x-P)(x-Q)(x-R)(x-S)=x^4+ax^3+bx^2+cx+d</math>
The equation <math>\ f(x)=0</math> has the four solutions <math>\ x=P,Q,R,S</math> because a product is zero whenever one of the factors is zero. <math>\ P,Q,R,S</math>are the unknown roots of the polynomial.
Note that <math>\ P=x-\frac{f(x)}{(x-Q)(x-R)(x-S)}</math> for all <math>\ x</math>
and that <math>\ P=P-\frac{f(P)}{(P-q)(P-r)(P-s)}</math> for all <math>\ q,r,s</math>.
This is the clue to the method.
Initialize four new variables <math>\ p,q,r,s</math>:
- <math>\ p:=(0.4+i0.9)^0</math>
- <math>\ q:=(0.4+i0.9)^1</math>
- <math>\ r:=(0.4+i0.9)^2</math>
- <math>\ s:=(0.4+i0.9)^3</math>
There is nothing special about choosing <math>\ 0.4+i0.9</math> except that it is neither a real number nor a root of unity.
Make the substitutions:
- <math>\ p:=p-\frac{f(p)}{(p-q)(p-r)(p-s)}</math>
- <math>\ q:=q-\frac{f(q)}{(q-p)(q-r)(q-s)}</math>
- <math>\ r:=r-\frac{f(r)}{(r-p)(r-q)(r-s)}</math>
- <math>\ s:=s-\frac{f(s)}{(s-p)(s-q)(s-r)}</math>
Re-iterate until things calm down. Then <math>\ p,q,r,s</math> have the values <math> \ P,Q,R,S</math> in some order. 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.
