Muller's method

Last updated October 11, 2024

Muller's method is a root-finding algorithm, a numerical method for solving equations of the form f(x) = 0. It was first presented by David E. Muller in 1956.

Muller's method proceeds according to a third-order recurrence relation similar to the second-order recurrence relation of the secant method. Whereas the secant method proceeds by constructing a line through two points on the graph of f corresponding to the last two iterative approximations and then uses the line's root as the next approximation at every iteration, by contrast, Muller's method uses three points corresponding to the last three iterative approximations, constructs a parabola through these three points, and then uses a root of the parabola as the next approximation at every iteration.

Recurrence relation

Muller's method is a recursive method that generates a new approximation of a root ξ of f at each iteration using the three prior iterations. Starting with three initial values x₀, x₋₁ and x₋₂, the first iteration calculates an approximation x₁ using those three, the second iteration calculates an approximation x₂ using x₁, x₀ and x₋₁, the third iteration calculates an approximation x₃ using x₂, x₁ and x₀, and so on: the k^th iteration generates approximation x_k using x_k-1, x_k-2, and x_k-3. Each iteration takes as input the last three generated approximations and the value of f at these approximations: the values x_k-1, x_k-2 and x_k-3 and the function values f(x_k-1), f(x_k-2) and f(x_k-3). The approximation x_k is calculated as follows from those six values.

First, a parabola y_k(x) is constructed by interpolating through the three points (x_k-1, f(x_k-1)), (x_k-2, f(x_k-2)) and (x_k-3, f(x_k-3)) using a Newton polynomial. y_k(x) is

y_{k}(x)=f(x_{k-1})+(x-x_{k-1})f[x_{k-1},x_{k-2}]+(x-x_{k-1})(x-x_{k-2})f[x_{k-1},x_{k-2},x_{k-3}],\,

where f[x_k-1, x_k-2] and f[x_k-1, x_k-2, x_k-3] denote divided differences. This can be rewritten as

y_{k}(x)=f(x_{k-1})+w(x-x_{k-1})+f[x_{k-1},x_{k-2},x_{k-3}]\,(x-x_{k-1})^{2}\,

where

w=f[x_{k-1},x_{k-2}]+f[x_{k-1},x_{k-3}]-f[x_{k-2},x_{k-3}].\,

The iterate x_k is then given as the solution of the quadratic equation y_k(x) = 0 closest to x_k-1. Altogether, this implies the overall nonlinear third-order recurrence relation ^[1]

x_{k}=x_{k-1}-{\frac {2f(x_{k-1})}{w\pm {\sqrt {w^{2}-4f(x_{k-1})f[x_{k-1},x_{k-2},x_{k-3}]}}}},

where the sign of the square root should be chosen such that the total denominator is as large as possible in magnitude.

Note that x_k can be complex even when the previous iterates are all real. This is in contrast with other root-finding algorithms like the secant method, Sidi's generalized secant method or Newton's method, whose iterates will remain real if one starts with real numbers. Having complex iterates can be an advantage (if one is looking for complex roots) or a disadvantage (if it is known that all roots are real), depending on the problem.

Speed of convergence

For well-behaved functions, the order of convergence of Muller's method is approximately 1.84 and exactly the tribonacci constant. This can be compared with approximately 1.62, exactly the golden ratio, for the secant method and with exactly 2 for Newton's method. So, the secant method makes less progress per iteration than Muller's method and Newton's method makes more progress.

More precisely, if ξ denotes a single root of f (so f(ξ) = 0 and f'(ξ) ≠ 0), f is three times continuously differentiable, and the initial guesses x₀, x₁, and x₂ are taken sufficiently close to ξ, then the iterates satisfy

\lim _{k\to \infty }{\frac {|x_{k}-\xi |}{|x_{k-1}-\xi |^{\mu }}}=\left|{\frac {f'''(\xi )}{6f'(\xi )}}\right|^{(\mu -1)/2},

where μ ≈ 1.84 is the positive solution of $x^{3}-x^{2}-x-1=0$ , the defining equation for the tribonacci constant.

Generalizations and related methods

Muller's method fits a parabola, i.e. a second-order polynomial, to the last three obtained points f(x_k-1), f(x_k-2) and f(x_k-3) in each iteration. One can generalize this and fit a polynomial p_k,m(x) of degree m to the last m+1 points in the k^th iteration. Our parabola y_k is written as p_k,2 in this notation. The degree m must be 1 or larger. The next approximation x_k is now one of the roots of the p_k,m, i.e. one of the solutions of p_k,m(x)=0. Taking m=1 we obtain the secant method whereas m=2 gives Muller's method.

Muller calculated that the sequence {x_k} generated this way converges to the root ξ with an order μ_m where μ_m is the positive solution of $x^{m+1}-x^{m}-x^{m-1}-\dots -x-1=0$ .

The method is much more difficult though for m>2 than it is for m=1 or m=2 because it is much harder to determine the roots of a polynomial of degree 3 or higher. Another problem is that there seems no prescription of which of the roots of p_k,m to pick as the next approximation x_k for m>2.

These difficulties are overcome by Sidi's generalized secant method which also employs the polynomial p_k,m. Instead of trying to solve p_k,m(x)=0, the next approximation x_k is calculated with the aid of the derivative of p_k,m at x_k-1 in this method.

Computational example

Below, Muller's method is implemented in the Python programming language. It is then applied to find a root of the function $f (x) = x 2 - 612$ .

fromtypingimport*fromcmathimportsqrt# Use the complex sqrt as we may generate complex numbersNum=Union[float,complex]Func=Callable[[Num],Num]defdiv_diff(f:Func,xs:list[Num]):"""Calculate the divided difference f[x0, x1, ...]."""iflen(xs)==2:a,b=xsreturn(f(a)-f(b))/(a-b)else:return(div_diff(f,xs[1:])-div_diff(f,xs[0:-1]))/(xs[-1]-xs[0])defmullers_method(f:Func,xs:(Num,Num,Num),iterations:int)->float:"""Return the root calculated using Muller's method."""x0,x1,x2=xsfor_inrange(iterations):w=div_diff(f,(x2,x1))+div_diff(f,(x2,x0))-div_diff(f,(x2,x1))s_delta=sqrt(w**2-4*f(x2)*div_diff(f,(x2,x1,x0)))denoms=[w+s_delta,w-s_delta]# Take the higher-magnitude denominatorx3=x2-2*f(x2)/max(denoms,key=abs)# Advancex0,x1,x2=x1,x2,x3returnx3deff_example(x:Num)->Num:"""The example function. With a more expensive function, memoization of the last 4 points called may be useful."""returnx**2-612root=mullers_method(f_example,(10,20,30),5)print("Root: {}".format(root))# Root: (24.738633317099097+0j)

Related Research Articles

In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the i-th approximation is derived from the previous ones.

<span class="mw-page-title-main">Newton's method</span> Algorithm for finding zeros of functions

In numerical analysis, Newton's method, also known as the Newton–Raphson method, named after Isaac Newton and Joseph Raphson, is a root-finding algorithm which produces successively better approximations to the roots of a real-valued function. The most basic version starts with a real-valued function $f$ , its derivative $f'$ , and an initial guess $x 0$ for a root of $f$ . If $f$ satisfies certain assumptions and the initial guess is close, then

In mathematics, differential calculus is a subfield of calculus that studies the rates at which quantities change. It is one of the two traditional divisions of calculus, the other being integral calculus—the study of the area beneath a curve.

In numerical analysis, a root-finding algorithm is an algorithm for finding zeros, also called "roots", of continuous functions. A zero of a function $f$ is a number $x$ such that $f (x) = 0$ . As, generally, the zeros of a function cannot be computed exactly nor expressed in closed form, root-finding algorithms provide approximations to zeros. For functions from the real numbers to real numbers or from the complex numbers to the complex numbers, these are expressed either as floating-point numbers without error bounds or as floating-point values together with error bounds. The latter, approximations with error bounds, are equivalent to small isolating intervals for real roots or disks for complex roots.

In the mathematical subfield of numerical analysis, numerical stability is a generally desirable property of numerical algorithms. The precise definition of stability depends on the context. One is numerical linear algebra and the other is algorithms for solving ordinary and partial differential equations by discrete approximation.

In numerical analysis, the secant method is a root-finding algorithm that uses a succession of roots of secant lines to better approximate a root of a function f. The secant method can be thought of as a finite-difference approximation of Newton's method, so it is considered a quasi-Newton method. Historically, it is as an evolution of the method of false position, which predates Newton's method by over 3000 years.

In mathematics, the regula falsi, method of false position, or false position method is a very old method for solving an equation with one unknown; this method, in modified form, is still in use. In simple terms, the method is the trial and error technique of using test ("false") values for the variable and then adjusting the test value according to the outcome. This is sometimes also referred to as "guess and check". Versions of the method predate the advent of algebra and the use of equations.

In numerical analysis, inverse quadratic interpolation is a root-finding algorithm, meaning that it is an algorithm for solving equations of the form f(x) = 0. The idea is to use quadratic interpolation to approximate the inverse of f. This algorithm is rarely used on its own, but it is important because it forms part of the popular Brent's method.

In mathematics, the conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is positive-semidefinite. The conjugate gradient method is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation or other direct methods such as the Cholesky decomposition. Large sparse systems often arise when numerically solving partial differential equations or optimization problems.

In numerical optimization, the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an iterative method for solving unconstrained nonlinear optimization problems. Like the related Davidon–Fletcher–Powell method, BFGS determines the descent direction by preconditioning the gradient with curvature information. It does so by gradually improving an approximation to the Hessian matrix of the loss function, obtained only from gradient evaluations via a generalized secant method.

Methods of computing square roots are algorithms for approximating the non-negative square root $of a positive real number . Since all square roots of natural numbers, other than of perfect squares, are irrational, square roots can usually only be computed to some finite precision: these methods typically construct a series of increasingly accurate approximations.$

In numerical analysis, fixed-point iteration is a method of computing fixed points of a function.

In numerical analysis, Aitken's delta-squared process or Aitken extrapolation is a series acceleration method used for accelerating the rate of convergence of a sequence. It is named after Alexander Aitken, who introduced this method in 1926. It is most useful for accelerating the convergence of a sequence that is converging linearly. A precursor form was known to Seki Kōwa and applied to the rectification of the circle, i.e., to the calculation of π.

In numerical analysis, a quasi-Newton method is an iterative numerical method used either to find zeroes or to find local maxima and minima of functions via an iterative recurrence formula much like the one for Newton's method, except using approximations of the derivatives of the functions in place of exact derivatives. Newton's method requires the Jacobian matrix of all partial derivatives of a multivariate function when used to search for zeros or the Hessian matrix when used for finding extrema. Quasi-Newton methods, on the other hand, can be used when the Jacobian matrices or Hessian matrices are unavailable or are impractical to compute at every iteration.

In the mathematical field of numerical analysis, monotone cubic interpolation is a variant of cubic interpolation that preserves monotonicity of the data set being interpolated.

In numerical analysis, Halley's method is a root-finding algorithm used for functions of one real variable with a continuous second derivative. Edmond Halley was an English mathematician and astronomer who introduced the method now called by his name.

In mathematics, numerical analysis, and numerical partial differential equations, domain decomposition methods solve a boundary value problem by splitting it into smaller boundary value problems on subdomains and iterating to coordinate the solution between adjacent subdomains. A coarse problem with one or few unknowns per subdomain is used to further coordinate the solution between the subdomains globally. The problems on the subdomains are independent, which makes domain decomposition methods suitable for parallel computing. Domain decomposition methods are typically used as preconditioners for Krylov space iterative methods, such as the conjugate gradient method, GMRES, and LOBPCG.

<span class="mw-page-title-main">Fast inverse square root</span> Root-finding algorithm

Fast inverse square root, sometimes referred to as Fast InvSqrt or by the hexadecimal constant 0x5F3759DF, is an algorithm that estimates $, the reciprocal of the square root of a 32-bit floating-point number in IEEE 754 floating-point format. The algorithm is best known for its implementation in 1999 in Quake III Arena, a first-person shooter video game heavily based on 3D graphics. With subsequent hardware advancements, especially the x86 SSE instruction rsqrtss, this algorithm is not generally the best choice for modern computers, though it remains an interesting historical example.$

LOOP is a simple register language that precisely captures the primitive recursive functions. The language is derived from the counter-machine model. Like the counter machines the LOOP language comprises a set of one or more unbounded registers, each of which can hold a single non-negative integer. A few arithmetic instructions operate on the registers. The only control flow instruction is 'LOOP x DO...END'. It causes the instructions within its scope to be repeated x times.

Sidi's generalized secant method is a root-finding algorithm, that is, a numerical method for solving equations of the form $. The method was published by Avram Sidi.$

References

↑ Muller, David E. (1956). "A method for solving algebraic equations using an automatic computer". Mathematical Tables and Other Aids to Computation. 10 (56): 208–215. doi: 10.1090/S0025-5718-1956-0083822-0 . JSTOR 2001916. MR 0083822.

Atkinson, Kendall E. (1989). An Introduction to Numerical Analysis, 2nd edition, Section 2.4. John Wiley & Sons, New York. ISBN 0-471-50023-2.
Burden, R. L. and Faires, J. D. Numerical Analysis, 4th edition, pages 77ff.
Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 9.5.2. Muller's Method". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8.

v t e Root-finding algorithms
Bracketing (no derivative)	Bisection method Regula falsi ITP method
Householder	Newton's method Halley's method
Quasi-Newton	Broyden's method Secant method Newton–Krylov method Steffensen's method
Hybrid methods	Brent's method Ridders' method
Polynomial methods	Aberth method Bairstow's method Durand–Kerner method Graeffe's method Jenkins–Traub algorithm Lehmer–Schur algorithm Laguerre's method Splitting circle method
Other methods	Fixed-point iteration Inverse quadratic interpolation Muller's method Sidi's generalized secant method