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.
Muller's method uses three initial approximations of the root, and , and determines the next approximation by considering the intersection of the x-axis with the parabola through , and .
Consider the quadratic polynomial
that passes through , and . To simplify notation, define the differences
and
Substituting each of the three points , and into and solving simultaneously for and gives
The quadratic formula is then applied to to determine as
The sign preceding the radical term is chosen to match the sign of to ensure the next iterate is closest to , giving
Once is determined, the process is repeated. Note that due to the radical expression in the denominator, iterates 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.
The method can easily be adopted for scenarios, such as collision detection, where not the closest root is of interest, but the smaller one by replacing with :
The following pseudocode describes Muller's method for approximating a root of the equation f(x) = 0. It returns an estimate that satisfies , where is a user-specified tolerance.
input: function f, initial approximations x0, x1, x2, tolerance õ, maximum iterations Nmax for n = 0 to Nmax do h0 â x1 â x0 h1 â x2 â x1 ô0 â (f(x1) â f(x0)) / h0 ô1 â (f(x2) â f(x1)) / h1 a â (ô1 â ô0) / (h1 + h0) b â a÷h1 + ô1 c â f(x2) if b âÂÂ¥ 0 then denominator â b + sqrt(bò â 4÷a÷c) else denominator â b â sqrt(bò â 4÷a÷c) end if x3 â x2 â 2÷c / denominator if |x3 â x2| < õ then return x3 end if x0 â x1 x1 â x2 x2 â x3 end for return x3
For well-behaved functions, the order of convergence of Muller's method is approximately 1.839 and exactly the tribonacci constant. This can be compared with approximately 1.618, 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<nowiki>'</nowiki>(þ) â 0), f is three times continuously differentiable, and the initial guesses x<sub>0</sub>, x<sub>1</sub>, and x<sub>2</sub> are taken sufficiently close to þ, then the iterates satisfy
where ü â 1.84 is the positive solution of , the defining equation for the tribonacci constant.
Muller's method fits a parabola, i.e. a second-order polynomial, to the last three obtained points f(x<sub>k-1</sub>), f(x<sub>k-2</sub>) and f(x<sub>k-3</sub>) in each iteration. One can generalize this and fit a polynomial p<sub>k,m</sub>(x) of degree m to the last m+1 points in the k<sup>th</sup> iteration. Our parabola y<sub>k</sub> is written as p<sub>k,2</sub> in this notation. The degree m must be 1 or larger. The next approximation x<sub>k</sub> is now one of the roots of the p<sub>k,m</sub>, i.e. one of the solutions of p<sub>k,m</sub>(x)=0. Taking m=1 we obtain the secant method whereas m=2 gives Muller's method.
Muller calculated that the sequence {x<sub>k</sub>} generated this way converges to the root þ with an order ü<sub>m</sub> where ü<sub>m</sub> is the positive solution of .
As m approaches infinity the positive solution for the equation approaches 2. 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<sub>k,m</sub> to pick as the next approximation x<sub>k</sub> for m>2.
These difficulties are overcome by Sidi's generalized secant method which also employs the polynomial p<sub>k,m</sub>. Instead of trying to solve p<sub>k,m</sub>(x)=0, the next approximation x<sub>k</sub> is calculated with the aid of the derivative of p<sub>k,m</sub> at x<sub>k-1</sub> in this method.