Muller's method

Template:Short description Template:More footnotes Muller's method is a root-finding algorithm, a numerical method for solving equations of the form Template:Math. It was first presented by David E. Muller in 1956.

Animation illustrating Muller's method applied to the function Template:Math. At each iteration, a parabola interpolating three points is constructed, and one of its roots is used to generate the next approximation.

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 Template:Mvar 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.

Derivation

Muller's method uses three initial approximations of the root, $x_0,x_1$ and Template:Tmath, and determines the next approximation $x_3$ by considering the intersection of the x-axis with the parabola through Template:Tmath, $(x_1,f(x_1))$ and Template:Tmath.

Consider the quadratic polynomial $$P(x)=a(x-x_2{)}^2+b(x-x_2)+c,$$

that passes through Template:Tmath, $(x_1,f(x_1))$ and Template:Tmath. To simplify notation, define the differences

$$h_0=x_1-x_0,h_1=x_2-x_1$$

and $${\delta }_0=\frac{f(x_1)-f(x_0)}{h_0},{\delta }_1=\frac{f(x_2)-f(x_1)}{h_1}.$$

Substituting each of the three points Template:Tmath, $(x_1,f(x_1))$ and $(x_2,f(x_2))$ into $P(x)$ and solving simultaneously for $a,b$ and $c$ gives

$$a=\frac{{\delta }_1-{\delta }_0}{h_1+h_0},b=a{h}_1+{\delta }_1,c=f(x_2).$$

The quadratic formula is then applied to $P(x)$ to determine $x_3$ as

$$x_3-x_2=\frac{-2c}{b\pm \sqrt{b^2-4ac}}.$$

The sign preceding the radical term is chosen to match the sign of $b$ to ensure the next iterate is closest to Template:Tmath, giving

$$x_3=x_2-\frac{2c}{b+\mathrm{sign}(b)\sqrt{b^2-4ac}}.$$

Once $x_3$ 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 $\mathrm{sign}(b)$ with $\mathrm{sign}(-a)$: $$x_3=x_2-\frac{2c}{b-\mathrm{sign}(a)\sqrt{b^2-4ac}}.$$

Algorithm

The following pseudocode describes Muller's method for approximating a root of the equation Template:Math. It returns an estimate $x_3$ that satisfies Template:Tmath, where $\epsilon$ is a user-specified tolerance.

input: function Template:Mvar, initial approximations Template:Math, Template:Math, Template:Math, tolerance Template:Mvar, maximum iterations Template:Math
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

Example

Suppose we seek a solution of the polynomial $$f(x)=-x^3-x+7.$$

Applying Muller's method with initial guesses Template:Tmath, $x_1=2$ and $x_2=3$ and a tolerance of $10^{-9}$ gives the following approximations:

Iteration	$a$	$b$	$c$	Denominator	$x_n$	$f(x_n)$
$3$	$-6.0$	$-26.0$	$-23$	$-37.13552$	$1.7612940604716703$	$-2.2510\times 10^{-1}$
$4$	$-6.76129$	$-10.01078$	$-0.22510$	$-19.71273$	$1.7384555878656442$	$7.53555\times 10^{-3}$
$5$	$-6.49974$	$-10.03787$	$7.53555\times 10^{-3}$	$-20.08549$	$1.7392059353750402$	$-2.08961\times 10^{-5}$
$6$	$-5.23895$	$-10.07452$	$-2.08961\times 10^{-5}$	$-20.14903$	$1.7392038612200968$	$-3.42854\times 10^{-11}$

To the required tolerance, Muller's method generates the solution $1.7392038612200968$ on the sixth iteration.

By modifying the initial guesses, Muller's method can also generate complex solutions. Applying Muller's method again but with initial guesses Template:Tmath, $x_1=-2$ and $x_2=-3$ and a tolerance of $10^{-9}$ gives the following approximations:

Iteration	Denominator	$x_n$	$f(x_n)$
$3$	$-26-14.56021i$	$-0.8333333333333335-1.2133516482134197i$	$4.73148+1.95484i$
$4$	$6.38915-24.72107i$	$-0.777821059050328-1.6104880890500035i$	$2.19616+0.356473i$
$5$	$9.46045-19.15717i$	$-0.8389286926696546-1.8095901140685215i$	$-.187858-0.295347i$
$6$	$13.02224-18.42210i$	$-0.8699227902463772-1.8080759128505213i$	$-3.44378\times 10^{-3}+2.09461\times 10^{-3}i$
$7$	$13.08065-18.85860i$	$-0.8696017710815587-1.8079333552255168i$	$-1.67055\times 10^{-7}-2.34332\times 10^{-6}i$
$8$	$13.07447-18.86619i$	$-0.8696019306084998-1.8079332269621338i$	$-2.10054\times 10^{-12}+6.79900\times 10^{-13}i$

To the required tolerance, Muller's method generates the complex solution $-0.8696019306084998-1.8079332269621338i$ on the eighth iteration.

Speed of convergence

Suppose the real-valued function $f(x)$ has a root at Template:Tmath. Let $$e_n=x_n-\alpha$$ denote the error at the $n^{th}$ iteration of Muller's method. By definition, we have $$e_{n+1}=K{e}_n^p$$ where $p$ is the order of convergence. Therefore, $$e_n=K{e}_{n-1}^p$$ and $$e_{n-1}=K{e}_{n-2}^p$$ which gives $$e_{n-1}\sim e_n^{1/p}$$ and $$e_{n-2}\sim e_n^{1/p^2}.$$ Near the root Template:Tmath, Taylor's theorem gives $$f(x)=f^{\prime }(\alpha )(x-\alpha )+\frac{f^{″}(\alpha )}{2}(x-\alpha {)}^2+R(x)$$ where $R(x)$ is the remainder term, defined as $$R(x)=\frac{f^{‴}(\xi )}{6}(x-\alpha {)}^3$$ for some $\xi$ between $x$ and Template:Tmath. Using the three consecutive iterates of Mullers method, Template:Tmath, $x=x_{n-1}$ and $x=x_{n-2}$ evaluated at the root gives

$$R(\alpha )=O(e_n{e}_{n-1}{e}_{n-2})$$ or $$e_{n+1}=O(e_n{e}_{n-1}{e}_{n-2})$$ since the next iterate is obtained by setting the interpolating quadratic to zero. Asymptotically, $$e_{n+1}=K{e}_n{e}_{n-1}{e}_{n-2}.$$ Using $$e_{n-1}\sim e_n^{1/p}$$ and $$e_{n-2}\sim e_n^{1/p^2}$$ gives $$e_{n+1}=K{e}_n^{1+\frac{1}{p}+\frac{1}{p^2}}.$$

By definition, $$e_{n+1}=K{e}_n^p$$ so equating exponents gives $$p=1+\frac{1}{p}+\frac{1}{p^2}$$ which can be rearranged to give

$$p^3-p^2-p-1=0.$$

Solving this equation gives Template:Tmath, so the order of convergence of Muller's method is 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 $\alpha$ denotes a single root of $f(x)$ and $f^{\prime }(\alpha )\ne 0$ and $f(x)$ is three times continuously differentiable, if the initial guesses Template:Tmath, $x_1$ and $x_2$ are taken sufficiently close to Template:Tmath, then the iterates satisfy $$\underset{k\to \mathrm{\infty }}{lim}\frac{|x_k-\alpha |}{|x_{k-1}-\alpha {|}^{\mu }}={\left|\frac{f^{‴}(\alpha )}{6f^{\prime }(\alpha )}\right|}^{(\mu -1)/2}$$

where $\mu \approx 1.84$ is the positive solution of Template:Tmath, 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 Template:Math, Template:Math and Template:Math in each iteration. One can generalize this and fit a polynomial Template:Math of degree Template:Mvar to the last Template:Math points in the Template:Mvar-th iteration. Our parabola Template:Math is written as Template:Math in this notation. The degree Template:Mvar must be 1 or larger. The next approximation Template:Math is now one of the roots of the Template:Math, i.e. one of the solutions of Template:Math. Taking Template:Math we obtain the secant method whereas Template:Math gives Muller's method.

Muller calculated that the sequence Template:Math generated this way converges to the root ξ with an order Template:Math where Template:Math is the positive solution of Template:Tmath.

As Template:Mvar approaches infinity the positive solution for the equation approaches 2. The method is much more difficult though for Template:Math than it is for Template:Math or Template:Math 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 Template:Math to pick as the next approximation Template:Math for Template:Math.

These difficulties are overcome by Sidi's generalized secant method which also employs the polynomial Template:Math. Instead of trying to solve Template:Math, the next approximation Template:Math is calculated with the aid of the derivative of Template:Math at Template:Math in this method.

See also

• Halley's method, with cubic convergence
• Householder's method, includes Newton's, Halley's and higher-order convergence

References

• Atkinson, Kendall E. (1989). An Introduction to Numerical Analysis, 2nd edition, Section 2.4. John Wiley & Sons, New York. Template:ISBN.
• Burden, R. L. and Faires, J. D. Numerical Analysis, 4th edition, pages 77ff.
• Template:Cite book

Further reading

• A bracketing variant with global convergence: Template:Cite journal