ミューラー法

メソッド

$f (\alpha) = 0$ としよう。初期値 $x_{0} , x_{1} , x_{2}$ と $$ w_{n} := f [x_{n} , x_{n-1} ] + f [ x_{n} , x_{n-2} ] - f [ x_{n-2} , x_{n-1} ] $$ によって、 $$ 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} ] } }} $$ のように定義された数列 $\left\{ x_{n} \right\}$ は $n \to \infty$ の時、$\alpha$ に $p \approx 1.84$ 次で収束する。

• ただし、$\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}$ は $+$ と $-$ の大きい方を選ぶ。
• 複素数の絶対値は $|x + iy| := \sqrt{ x^2 + y^2}$ のように計算される。

説明 ¹

ミュラーメソッドはセカントメソッドの応用と見ることができるが、まだニュートン-ラプソンメソッドより遅いが、複素数の解を見つけることができる特徴がある。微分や符号を使った解法ではなく、根の公式を使って解くため、幾何学的な意味にこだわる必要はなく、複素数をカバーする。

例

与えられた$\left( x_{i} , f (x_{i}) \right)_{i = n-1}^{n+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} ] $$ は一意に存在し、$p ( x_{i} ) = f(x_{i} )$ をよく満たしている。これを $(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*} $$

従って、 $$ \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*} $$

$x = x_{n+1}$ を代入すると、$p_{2} (x_{n+1}) = f(x_{n+1}) + 0 + 0$ となるため、$p_{2} (x) = 0$ となる $(x - x_{n+1})$ を見つければ、結局 $f(x_{n+1} ) = 0$ を見つけるのと同じだ。元々見つけたいものが $\displaystyle \lim_{n \to \infty} f ( x_{n} ) = 0$ を満たす $x_{n}$ だったので、 $$ 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} ] $$ における $(x - x_{2})$ に根の公式を使えばいい。根の公式で $$ 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} ] }} $$ を得て、実際の計算はコンピュータが行うため、有効数字などの問題で逆有理化を逆に使って $$ 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} ] } }} $$ を公式として使う。

一方で、$f '$が必ず存在する必要はないが、$f’ ( \alpha ) \ne 0$ならば、$n \to \infty$ の時、$w_{n} \to f ' ( \alpha )$であるため、微分係数を見つける用途で使うこともできる。

実装

以下はミュラーメソッドをRで実装したものである。理論上は3点 $(x_{0} , x_{1} , x_{2})$ を持って始めるが、実際には $x_{0} , x_{1}$ だけでも $\displaystyle x_{2} := {{x_{0} + x_{1}} \over {2}}$ を決定できるので、2つの初期値だけを受け取って使う。この時、$(x_{0} , x_{1})$ が実際に実数であっても、データ型は複素数でなければならない。

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)

以下は上記のコードを実行した結果である。$f(x) = x^2 + 9 =$の虚根の一つである$x = -3 i$と、$g(x) = x^2 + x + 1 = 0$の複素根の一つである$\displaystyle x = - {{1 + \sqrt{3} i } \over {2}}$をよく見つけ出したことが確認できる。