Mueller Method

Method

Let $f (\alpha) = 0$. With the initial value $x_{0} , x_{1} , x_{2}$ and $$ w_{n} := f [x_{n} , x_{n-1} ] + f [ x_{n} , x_{n-2} ] - f [ x_{n-2} , x_{n-1} ] $$ , define the sequence $\left\{ x_{n} \right\}$ as $$ x_{n+1} : = x_{n} - {{ 2 f ( x_{n} ) } \over { w_{n} \pm \sqrt{ w_{n}^{2} - 4 f (x_{n} ) f [ x_{n} , x_{n-1} , x_{n-2} ] } }} $$, which converges to $\alpha$ of order $p \approx 1.84$ when $n \to \infty$.

• Note, $\left( w_{n} \pm \sqrt{ w_{n}^{2} - 4 f (x_{n} ) f [ x_{n} , x_{n-1} , x_{n-2} ] } \right) \in \mathbb{C}$ is selected as the larger of $+$ and $-$.
• The modulus of a complex number is calculated as $|x + iy| := \sqrt{ x^2 + y^2}$.

Description ¹

The Muller method can be seen as an application of the secant method, still slower than the Newton-Raphson method, but it has the characteristic of being able to find complex roots. It solves using the root formula instead of derivatives or signs, so there’s no need to hang on to geometric meanings, and it covers complex numbers.

Example

Given $\left( x_{i} , f (x_{i}) \right)_{i = n-1}^{n+1}$, the unique quadratic polynomial $$ p_{2} (x) = f(x_{n+1} ) + (x - x_{n+1} ) f [ x_{n+1} , x_{n} ] + ( x - x_{n+1} ) ( x - x_{n} ) f [ x_{n+1} , x_{n} , x_{n-1} ] $$ satisfies $p ( x_{i} ) = f(x_{i} )$ well. Let’s transform this into a quadratic equation for $(x - x_{n+1})$.

$$ \begin{align*} & f [ x_{n+1} , x_{n} , x_{n-1} ] = { { f [x_{n+1} , x_{n-1} ] - f [ x_{n-1} , x_{n} ] } \over { x_{n+1} - x_{n} }} \\ \implies& ( x_{n+1} - x_{n} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ & = f [x_{n+1} , x_{n-1} ] - f [ x_{n-1} , x_{n} ] \\ \implies& ( x - x_{n} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ & = f [x_{n+1} , x_{n-1} ] - f [ x_{n-1} , x_{n} ] + ( x - x_{n+1} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ \implies& f [x_{n+1} , x_{n} ] + ( x - x_{n} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ & = f [x_{n+1} , x_{n} ] + f [x_{n+1} , x_{n-1} ] - f [ x_{n-1} , x_{n} ] + ( x - x_{n+1} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ \implies& f [x_{n+1} , x_{n} ] + ( x - x_{n} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ & = w_{n} + ( x - x_{n+1} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ \implies& ( x - x_{n+1} ) f [x_{n+1} , x_{n} ] + ( x - x_{n} ) ( x - x_{n+1} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ & = w_{n} ( x - x_{n+1} ) + ( x - x_{n+1} )^2 f [ x_{n+1} , x_{n} , x_{n-1} ] \end{align*} $$

Therefore, $$ \begin{align*} p_{2} (x) =& f(x_{n+1} ) + (x - x_{n+1} ) f [ x_{n+1} , x_{n} ] + ( x - x_{n+1} ) ( x - x_{n} ) f [ x_{n+1} , x_{n} , x_{n-1} ] \\ =& f(x_{n+1} ) + w_{n} (x - x_{n+1} ) + ( x - x_{n+1} )^2 f [ x_{n+1} , x_{n} , x_{n-1} ] \end{align*} $$

Substituting $x = x_{n+1}$ gives $p_{2} (x_{n+1}) = f(x_{n+1}) + 0 + 0$, so finding $(x - x_{n+1})$ that makes it $p_{2} (x) = 0$ is essentially the same as finding $f(x_{n+1} ) = 0$. Since the original goal was to find $x_{n}$ satisfying $\displaystyle \lim_{n \to \infty} f ( x_{n} ) = 0$, for $$ p_{2} (x) = f(x_{n+1} ) + w_{n} (x - x_{n+1} ) + ( x - x_{n+1} )^2 f [ x_{n+1} , x_{n} , x_{n-1} ] $$, use the quadratic formula for $(x - x_{2})$. The formula yields $$ x_{n+1} - x_{n} = {{ - w_{n} \pm \sqrt{ w_{n}^{2} - 4 f (x_{n} ) f [ x_{n} , x_{n-1} , x_{n-2} ] } } \over { 2 f [ x_{2} , x_{1} , x_{0} ] }} $$, which, because actual computation is done by computers, due to issues like significant figures, advises to use the reverse rationalized $$ x_{n+1} = x_{n} - {{ 2 f ( x_{n} ) } \over { w_{n} \pm \sqrt{ w_{n}^{2} - 4 f (x_{n} ) f [ x_{n} , x_{n-1} , x_{n-2} ] } }} $$ as the formula.

Meanwhile, it’s not necessary for $f '$ to exist, but if $f’ ( \alpha ) \ne 0$, then it’s $w_{n} \to f ' ( \alpha )$ when $n \to \infty$, hence it can be used to find derivatives as well.

Implementation

Below is the implementation of the Muller method in R. In theory, it starts with three points $(x_{0} , x_{1} , x_{2})$, but in practice, it can determine $\displaystyle x_{2} := {{x_{0} + x_{1}} \over {2}}$ with just $x_{0} , x_{1}$, so it uses only two initial values. Even if $(x_{0} , x_{1})$ is actually real, the data type must be complex.

dd<-function(f,X){
  if(length(X)==1) {return(f(X))}
  temp<-numeric(0)
  for(i in 1:length(X)) {temp[i]<-prod(X[i]-X[-i])}
  return(sum(f(X)/temp))
}
 
Muller<-function(f,x0,x1)
{
  x2<-(x0+x1)/2
  
  while(Mod(f(x2))>10^(-3)){
    w = dd(f,c(x2,x1)) + dd(f,c(x2,x0)) - dd(f,c(x0,x1))
    p = w + sqrt(w^2 - 4*f(x2)*dd(f,c(x2,x1,x0)))
    m = w - sqrt(w^2 - 4*f(x2)*dd(f,c(x2,x1,x0)))
    
    if(Mod(p)>Mod(m)) {w<-p} else {w<-m}
    
    x0<-x1
    x1<-x2
    x2 <- (x2 - (( 2*f(x2) ) / w))
  }
  return(x2)
}
 
f<-function(x) {x^2 + 9}
Muller(f,-6+0i,-5+0i)
 
g<-function(x) {x^2 + x + 1}
Muller(g,-1i,-2+0i)

The following is the result of executing the above code. It successfully found one of the imaginary roots of $f(x) = x^2 + 9 =$, $x = -3 i$, and one of the complex roots of $g(x) = x^2 + x + 1 = 0$, $\displaystyle x = - {{1 + \sqrt{3} i } \over {2}}$.