2.5 Muller’s method#

The secant method uses a linear function that passes through \((p_{0},f(p_{0}))\) and \((p_{1},f(p_{1}))\) to find the next iterate \(p_{2}.\) Muller’s method takes three initial approximations, passes a parabola (quadratic polynomial) through \((p_{0},f(p_{0})),\) \((p_{1},f(p_{1}))\), \((p_{2},f(p_{2}))\), and uses one of the roots of the polynomial as the next iterate.

Let the quadratic polynomial be written in the following form

(15)#\[P(x)=a(x-p_{2})^{2}+b(x-p_{2})+c.\]

Solve the following equations for \(a,b,c\)

\[\begin{align*} P(p_{0})=f(p_{0}) & =a(p_{0}-p_{2})^{2}+b(p_{0}-p_{2})+c\\ P(p_{1})=f(p_{1}) & =a(p_{1}-p_{2})^{2}+b(p_{1}-p_{2})+c\\ P(p_{2})=f(p_{2}) & =c \end{align*}\]

to get

(16)#\[\begin{split}\begin{align} c & =f(p_{2}) \nonumber \\ b & = \frac{(p_0-p_2)(f(p_1)-f(p_2))}{(p_1-p_2)(p_0-p_1)}-\frac{(p_1-p_2)(f(p_0)-f(p_2))}{(p_0-p_2)(p_0-p_1)}\\ %b & =\frac{(p_{0}-p_{2})^{2}\left(f(p_{1})-f(p_{2})\right)-(p_{1}-p_{2})^{2}\left(f(p_{0})-f(p_{2})\right)}{(p_{0}-p_{2})(p_{1}-p_{2})(p_{0}-p_{1})}\\ a &= \frac{f(p_0)-f(p_2)}{(p_0-p_2)(p_0-p_1)}-\frac{f(p_1)-f(p_2)}{(p_1-p_2)(p_0-p_1)}.\nonumber %a & =\frac{(p_{1}-p_{2})\left(f(p_{0})-f(p_{2})\right)-(p_{0}-p_{2})\left(f(p_{1})-f(p_{2})\right)}{(p_{0}-p_{2})(p_{1}-p_{2})(p_{0}-p_{1})}. \end{align}\end{split}\]

Now that we have determined \(P(x)\), the next step is to solve \(P(x)=0\), and set the next iterate \(p_{3}\) to its solution. To this end, put \(w=x-p_{2}\) in (15) to rewrite the quadratic equation as

\[\begin{equation*} aw^{2}+bw+c=0. \end{equation*}\]

From the quadratic formula, we obtain the roots

(17)#\[\hat{w}=\hat{x}-p_{2}=\frac{-2c}{b\pm\sqrt{b^{2}-4ac}}.\]

Let \(\Delta=b^2-4ac\). We have two roots (which could be complex numbers), \(-2c/(b+\sqrt{\Delta})\) and \(-2c/(b-\sqrt{\Delta})\), and we need to pick one of them. We will pick the root that is closer to \(p_2\), in other words, the root that makes \(|\hat{x}-p_2|\) the smallest. (If the numbers are complex, the absolute value means the norm of the complex number.) Therefore we have

(18)#\[\begin{split}\hat{x}-p_2 = \begin{cases} \frac{-2c}{b+\sqrt{\Delta}} & \text{if } |b+\sqrt{\Delta}|>|b-\sqrt{\Delta}|\\ \frac{-2c}{b-\sqrt{\Delta}} & \text{if } |b+\sqrt{\Delta}|\leq |b-\sqrt{\Delta}| \end{cases}.\end{split}\]

The next iterate of Muller’s method, \(p_3\), is set to the value of \(\hat{x}\) obtained from the above calculation, that is,

\[\begin{equation*} p_3=\hat{x} = \begin{cases} p_2-\frac{2c}{b+\sqrt{\Delta}} & \text{if } |b+\sqrt{\Delta}|>|b-\sqrt{\Delta}|\\ p_2-\frac{2c}{b-\sqrt{\Delta}} & \text{if } |b+\sqrt{\Delta}|\leq |b-\sqrt{\Delta}| \end{cases}. \end{equation*}\]

Remark 5

  1. Muller’s method can find real as well as complex roots.

  2. The convergence of Muller’s method is superlinear, that is,

    \[\begin{equation*} \lim_{n\rightarrow\infty}\frac{|p-p_{n+1}|}{|p-p_{n}|^{\alpha}}=\left|\frac{f^{(3)}(p)}{6f'(p)}\right|^{\frac{\alpha-1}{2}} \end{equation*}\]

    where \(\alpha\approx1.84\), provided \(f\in C^{3}[a,b],p\in(a,b)\), and \(f'(p)\neq0\).

  3. Muller’s method converges for a variety of starting values even though pathological examples that do not yield convergence can be found (for example, when the three starting values fall on a line).

Python code for Muller’s method#

The following Python code takes initial guesses \(p_0,p_1,p_2\) (written as pzero, pone, ptwo in the code), computes the coefficients \(a,b,c\) from Equation (16), and sets the root \(p_3\) to \(p\). It then updates the three initial guesses as the last three iterates, and continues until the stopping criterion is satisfied.

We need to compute the square root, and the absolute value, of possibly complex numbers in Equations (17) and (18). The Python function for the square root of a possibly complex number \(z\) is \(\textbf{complex(z)}^{0.5}\), and its absolute value is np.abs(z).

import numpy as np
def muller(f, pzero, pone, ptwo, eps, N):
    n = 1
    p = 0
    while n <= N:
        c = f(ptwo)
        b1 = (pzero-ptwo) * (f(pone)-f(ptwo)) / ((pone-ptwo)*(pzero-pone))
        b2 = (pone-ptwo) * (f(pzero)-f(ptwo)) / ((pzero-ptwo)*(pzero-pone))
        b = b1 - b2
        a1 = (f(pzero)-f(ptwo)) / ((pzero-ptwo)*(pzero-pone))
        a2 = (f(pone)-f(ptwo)) / ((pone-ptwo)*(pzero-pone))
        a = a1 - a2
        d = (complex(b**2-4*a*c))**0.5
        if np.abs(b-d) < np.abs(b+d):
            inc = 2*c/(b+d)
        else:
            inc = 2*c/(b-d)
        p = ptwo - inc
        if np.isclose(f(p), 0) or np.abs(p-ptwo)<eps:
            print('p is ', p, ' and the iteration number is ', n)
            return
        pzero = pone
        pone = ptwo
        ptwo = p
        n += 1
    y = f(p)
    print('Method did not converge. The last iteration gives ', 
          p, ' with function value ', y)

The polynomial \(x^5+2x^3-5x-2\) has three real roots, and two complex roots that are conjugates. Let’s find them all, by experimenting with various initial guesses.

muller(lambda x: x**5+2*x**3-5*x-2, 0.5, 1.0, 1.5, 1e-5, 10)

p is  (1.3196411677283386+0j)  and the iteration number is  4

muller(lambda x: x**5+2*x**3-5*x-2, 0.5, 0, -0.1, 1e-5, 10)

p is  (-0.43641313299908585+0j)  and the iteration number is  5

muller(lambda x: x**5+2*x**3-5*x-2, 0, -0.1, -1, 1e-5, 10)

p is  (-1+0j)  and the iteration number is  1

muller(lambda x: x**5+2*x**3-5*x-2, 5, 10, 15, 1e-5, 20)

p is  (0.05838598289491982+1.8626227582154478j)  and the iteration number is  18