Muller's method

Muller's method uses three initial approximations of the root, $x\_\{0\},\; x\_\{1\}$ and $x\_\{2\}$, and determines the next approximation $x\_\{3\}$ by considering the intersection of the x-axis with the parabola through $(x\_\{0\},\; f(x\_\{0\}))$, $(x\_\{1\},\; f(x\_\{1\}))$ and $(x\_\{2\},\; f(x\_\{2\}))$.

Consider the quadratic polynomial

{{NumBlk|:|$P(x)\; =\; a(x\; -\; x\_\{2\})^2\; +\; b(x\; -\; x\_\{2\})\; +\; c,$|{{EquationRef|1}}}}

that passes through $(x\_\{0\},\; f(x\_\{0\}))$, $(x\_\{1\},\; f(x\_\{1\}))$ and $(x\_\{2\},\; f(x\_\{2\}))$. Define the differences

h_{0} = x_{1} - x_{0}, \quad h_{1} = x_{2} - x_{1}

and

\delta_{0} = \frac{f(x_{1}) - f(x_{0})}{h_{0}}, \quad

\delta_{1} = \frac{f(x_{2}) - f(x_{1})}{h_{1}}.

Substituting each of the three points $(x\_\{0\},\; f(x\_\{0\}))$, $(x\_\{1\},\; f(x\_\{1\}))$ and $(x\_\{2\},\; f(x\_\{2\}))$ into equation ({{EquationNote|1}}) and solving simultaneously for $a,\; b$ and $c$ gives

a = \frac{\delta_{1} - \delta_{0}}{h_{1} + h_{0}}, \quad b = ah_{1} + \delta_{1}, \quad c = f(x_{2})

The quadratic formula is then applied to ({{EquationNote|1}}) 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 $x\_\{2\}$, giving

:$$

x_{3} = x_{2} - \frac{2c}{b + 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.

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'(ξ) ≠ 0), f is three times continuously differentiable, and the initial guesses x₀, x₁, and x₂ are taken sufficiently close to ξ, then the iterates satisfy

:$\backslash lim\_\{k\backslash to\backslash infty\}\; \backslash frac$

where μ ≈ 1.84 is the positive solution of $x^3\; -\; x^2\; -\; x\; -\; 1\; =\; 0$, 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_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\}\; -\; \backslash dots\; -\; x\; -\; 1\; =\; 0$.

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

Below, Muller's method is implemented in the Python programming language. It takes as parameters the three initial estimates of the root, as well as the desired decimals places of accuracy and the maximum number of iterations. The program is then applied to find a root of the function {{math|f(x) {{=}} x² − 612}}.

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

{{Reflist}}

• 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.
• {{Cite book | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press | publication-place=New York | isbn=978-0-521-88068-8 | chapter=Section 9.5.2. Muller's Method | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=466}}

• A bracketing variant with global convergence: {{cite journal |last1=Costabile |first1=F. |last2=Gualtieri |first2=M.I. |last3=Luceri |first3=R. |title=A modification of Muller's method |journal=Calcolo |date=March 2006 |volume=43 |issue=1 |pages=39–50 |doi=10.1007/s10092-006-0113-9|s2cid=124772103 }}

{{Root-finding algorithms}}

{{DEFAULTSORT:Mullers method}}

Category:Root-finding algorithms