logo

ミューラー法 📂数値解析

ミューラー法

メソッド

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

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

説明 1

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

例

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

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

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

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

実装

以下はミュラーメソッドをRで実装したものである。理論上は3点 (x0,x1,x2)(x_{0} , x_{1} , x_{2}) を持って始めるが、実際には x0,x1x_{0} , x_{1} だけでも x2:=x0+x12\displaystyle x_{2} := {{x_{0} + x_{1}} \over {2}} を決定できるので、2つの初期値だけを受け取って使う。この時、(x0,x1)(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)=x2+9=f(x) = x^2 + 9 =の虚根の一つであるx=−3ix = -3 iと、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. ↩︎