logo

Mueller Method 📂Numerical Analysis

Mueller Method

Method

Let f(α)=0f (\alpha) = 0. With the initial value x0,x1,x2x_{0} , x_{1} , x_{2} and wn:=f[xn,xn−1]+f[xn,xn−2]−f[xn−2,xn−1] w_{n} := f [x_{n} , x_{n-1} ] + f [ x_{n} , x_{n-2} ] - f [ x_{n-2} , x_{n-1} ] , define the sequence {xn}\left\{ x_{n} \right\} as xn+1:=xn−2f(xn)wn±wn2−4f(xn)f[xn,xn−1,xn−2] 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≈1.84p \approx 1.84 when n→∞n \to \infty.

  • Note, (wn±wn2−4f(xn)f[xn,xn−1,xn−2])∈C\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∣:=x2+y2|x + iy| := \sqrt{ x^2 + y^2}.

Description 1

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 (xi,f(xi))i=n−1n+1\left( x_{i} , f (x_{i}) \right)_{i = n-1}^{n+1}, the unique quadratic polynomial p2(x)=f(xn+1)+(x−xn+1)f[xn+1,xn]+(x−xn+1)(x−xn)f[xn+1,xn,xn−1] 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(xi)=f(xi)p ( x_{i} ) = f(x_{i} ) well. Let’s transform this into a quadratic equation for (x−xn+1)(x - x_{n+1}).

f[xn+1,xn,xn−1]=f[xn+1,xn−1]−f[xn−1,xn]xn+1−xn  ⟹  (xn+1−xn)f[xn+1,xn,xn−1]=f[xn+1,xn−1]−f[xn−1,xn]  ⟹  (x−xn)f[xn+1,xn,xn−1]=f[xn+1,xn−1]−f[xn−1,xn]+(x−xn+1)f[xn+1,xn,xn−1]  ⟹  f[xn+1,xn]+(x−xn)f[xn+1,xn,xn−1]=f[xn+1,xn]+f[xn+1,xn−1]−f[xn−1,xn]+(x−xn+1)f[xn+1,xn,xn−1]  ⟹  f[xn+1,xn]+(x−xn)f[xn+1,xn,xn−1]=wn+(x−xn+1)f[xn+1,xn,xn−1]  ⟹  (x−xn+1)f[xn+1,xn]+(x−xn)(x−xn+1)f[xn+1,xn,xn−1]=wn(x−xn+1)+(x−xn+1)2f[xn+1,xn,xn−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, p2(x)=f(xn+1)+(x−xn+1)f[xn+1,xn]+(x−xn+1)(x−xn)f[xn+1,xn,xn−1]=f(xn+1)+wn(x−xn+1)+(x−xn+1)2f[xn+1,xn,xn−1] \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=xn+1x = x_{n+1} gives p2(xn+1)=f(xn+1)+0+0p_{2} (x_{n+1}) = f(x_{n+1}) + 0 + 0, so finding (x−xn+1)(x - x_{n+1}) that makes it p2(x)=0p_{2} (x) = 0 is essentially the same as finding f(xn+1)=0f(x_{n+1} ) = 0. Since the original goal was to find xnx_{n} satisfying lim⁡n→∞f(xn)=0\displaystyle \lim_{n \to \infty} f ( x_{n} ) = 0, for p2(x)=f(xn+1)+wn(x−xn+1)+(x−xn+1)2f[xn+1,xn,xn−1] 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−x2)(x - x_{2}). The formula yields xn+1−xn=−wn±wn2−4f(xn)f[xn,xn−1,xn−2]2f[x2,x1,x0] 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 xn+1=xn−2f(xn)wn±wn2−4f(xn)f[xn,xn−1,xn−2] 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′f ' to exist, but if f’(α)≠0f’ ( \alpha ) \ne 0, then it’s wn→f′(α)w_{n} \to f ' ( \alpha ) when n→∞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 (x0,x1,x2)(x_{0} , x_{1} , x_{2}), but in practice, it can determine x2:=x0+x12\displaystyle x_{2} := {{x_{0} + x_{1}} \over {2}} with just x0,x1x_{0} , x_{1}, so it uses only two initial values. Even if (x0,x1)(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)=x2+9=f(x) = x^2 + 9 =, x=−3ix = -3 i, and one of the complex roots of g(x)=x2+x+1=0g(x) = x^2 + x + 1 = 0, x=−1+3i2\displaystyle x = - {{1 + \sqrt{3} i } \over {2}}.

20190321\_165937.png

  1. Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p73~74. ↩︎