Muller's method

Last updated September 22, 2023

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 is based on the secant method, which constructs at every iteration a line through two points on the graph of f. Instead, Muller's method uses three points, constructs the parabola through these three points, and takes the intersection of the x-axis with the parabola to be the next approximation.

Recurrence relation

Muller's method is a recursive method which generates an approximation of the root ξ of f at each iteration. Starting with the three initial values x₀, x₋₁ and x₋₂, the first iteration calculates the first approximation x₁, the second iteration calculates the second approximation x₂, the third iteration calculates the third approximation x₃, etc. Hence the k^th iteration generates approximation x_k. Each iteration takes as input the last three generated approximations and the value of f at these approximations. Hence the k^th iteration takes as input 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.

A parabola y_k(x) is constructed which goes 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)). When written in the Newton form, 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 next iterate x_k is now given as the solution closest to x_k-1 of the quadratic equation y_k(x) = 0. This yields the 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}]}}}}.

In this formula, the sign should be chosen such that the denominator is as large as possible in magnitude. We do not use the standard formula for solving quadratic equations because that may lead to loss of significance.

Note that x_k can be complex, even if the previous iterates were 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

The order of convergence of Muller's method is approximately 1.84. This can be compared with 1.62 for the secant method and 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$ .

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

<span class="mw-page-title-main">Numerical analysis</span> Study of algorithms using numerical approximation

Numerical analysis is the study of algorithms that use numerical approximation for the problems of mathematical analysis. It is the study of numerical methods that attempt at finding approximate solutions of problems rather than the exact ones. Numerical analysis finds application in all fields of engineering and the physical sciences, and in the 21st century also the life and social sciences, medicine, business and even the arts. Current growth in computing power has enabled the use of more complex numerical analysis, providing detailed and realistic mathematical models in science and engineering. Examples of numerical analysis include: ordinary differential equations as found in celestial mechanics, numerical linear algebra in data analysis, and stochastic differential equations and Markov chains for simulating living cells in medicine and biology.

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 single-variable function $f$ defined for a real variable $x$ , the function's derivative $f'$ , and an initial guess $x 0$ for a root of $f$ . If the function satisfies sufficient 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 mathematics and computing, a root-finding algorithm is an algorithm for finding zeros, also called "roots", of continuous functions. A zero of a function $f$ , from the real numbers to real numbers or from the complex numbers to the complex numbers, 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, expressed either as floating-point numbers or as small isolating intervals, or disks for complex roots (an interval or disk output being equivalent to an approximate output together with an error bound).

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. However, the secant method 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 number theory, the integer square root (isqrt) of a non-negative integer $n$ is the non-negative integer $m$ which is the greatest integer less than or equal to the square root of $n$ ,

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).

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-definite. 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 optimization, the line search strategy is one of two basic iterative approaches to find a local minimum $of an objective function . The other approach is trust region.$

Methods of computing square roots are numerical analysis algorithms for approximating the principal, or non-negative, square root of a real number. Arithmetically, it means given $, a procedure for finding a number which when multiplied by itself, yields; algebraically, it means a procedure for finding the non-negative root of the equation; geometrically, it means given two line segments, a procedure for constructing their geometric mean.$

In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a system of linear equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and is similar to the Jacobi method. Though it can be applied to any matrix with non-zero elements on the diagonals, convergence is only guaranteed if the matrix is either strictly diagonally dominant, or symmetric and positive definite. It was only mentioned in a private letter from Gauss to his student Gerling in 1823. A publication was not delivered before 1874 by Seidel.

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. Its early form was known to Seki Kōwa and was found for rectification of the circle, i.e. the calculation of π. It is most useful for accelerating the convergence of a sequence that is converging linearly.

Quasi-Newton methods are methods used to either find zeroes or local maxima and minima of functions, as an alternative to Newton's method. They can be used if the Jacobian or Hessian is unavailable or is too expensive to compute at every iteration. The "full" Newton's method requires the Jacobian in order to search for zeros, or the Hessian for finding extrema. Some iterative methods that reduce to Newton's method, such as SLSQP, may be considered quasi-Newtonian.

In numerical analysis, Halley's method is a root-finding algorithm used for functions of one real variable with a continuous second derivative. It is named after its inventor Edmond Halley.

<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
Newton	Newton's method
Quasi-Newton	Muller's method Secant method
Hybrid methods	Brent's method Ridders' method
Polynomial methods	Bairstow's method Jenkins–Traub method Laguerre's method