📂数値解析

ニュートン-ラプソン法

メソッド ¹

$f,f’,f’’$が$\alpha$の近くで連続であり、$f(\alpha) = 0, f '(\alpha) \ne 0$とする。

$\alpha$に十分近い初期値$x_{0}$に対して $$ x_{n+1} := x_{n} - {{ f ( x_{n} ) } \over { f ' ( x_{n} ) }} $$として定義された数列$\left\{ x_{n} \right\}$は、$n \to \infty$の時に$\alpha$に対してクアドラティックに収束する。

説明

ニュートン-ラフソン法は、ただのニュートン法とも呼ばれる。微分可能性や連続性などの条件があるが、それでも簡単で収束速度が速いため、方程式$f(x) = 0$の近似解を見つけるのに有用な方法である。

クアドラティックに収束するとは、収束率を語る時に$2$次で収束するという意味だが、適切な翻訳がない。「線形」に収束すると翻訳される$1$次と比べると、ただの$2$次で収束するという言葉はあまりにもあいまいで使用例が少ないため、やむを得ず英語表現をそのまま使った。それにもかかわらず、主に二分探索法より収束速度が速いことを数学的に示すことができるため、ここで言及された。

しかし、「$\alpha$の近く」とか「十分に近い初期値」という表現から分かるように、実際には考えなしに使える方法ではない。$0$で割る問題が発生するなど不運にも初期値を間違って与えたり、数列自体が解から遠く離れてしまう可能性もある。そのため、解を特定の回数内に求められなかった場合など、計算を諦めるなどの例外処理が必要だ。

このような発散エラーは、数学的に収束性が証明されていても、実際に使う場合は起こりうる問題であり、特にユーザーが「$\alpha$の近く」や「十分に近い初期値」を正確に把握できないことが原因である。二分探索法が遅くても確実に解を求めてくれるのと比較すると、明らかな欠点だと言える。

一緒に見る

• 高次元に対して一般化されたニュートン-ラフソン法
• ニュートン-フーリエ法: $f,f’,f’’$が閉区間$[a,b]$で連続であり、$f$が$[a,b]$で増加または減少して解を持つならば、$\displaystyle x_{n+1} = x_{n} - {{ f ( x_{n} ) } \over { f ' ( x_{n} ) }}$として定義された数列$\left\{ x_{n} \right\}$は、$n \to \infty$の時に$\alpha$にクアドラティックに収束する。

もう少し強い条件があれば、ニュートン-フーリエ法を考えることができ、収束性が完全に保証されつつクアドラティックに収束する。しかし、この方法は本質的にニュートン-ラフソン法の「$\alpha$の近く」に閉区間$[a,b]$を入れたものに過ぎない。したがって、条件をすべてチェックすれば、ニュートン-ラフソン法と完全に同じになり、別のコードを書く必要はない。

証明

戦略: テイラー展開を通じて$2$次の項を引き出した後、再帰的な不等式を作って収束性を示す。

第1部. テイラー展開

$f$を$x_{n}$でテイラー展開すると、$x$と$x_{n}$の間の$\xi_{n}$に対して $$ f(x) = f(x_{n} ) + (x - x_{n}) f '(x_{n}) + {{(x - x_{n})^2 } \over {2}} f '' (\xi_{n}) $$ $f( \alpha ) = 0$なので、$x = \alpha$を代入すると $$ 0 = f(x_{n} ) + ( \alpha - x_{n}) f '(x_{n}) + {{( \alpha - x_{n})^2 } \over {2}} f '' (\xi_{n}) $$ 両辺を$f ' (x_{n})$で割ると $$ 0 = {{f(x_{n} )} \over {f ' (x_{n})}} + \alpha - x_{n} + {{( \alpha - x_{n})^2 } \over {2}} {{f’’ (\xi_{n})} \over { f '(x_{n}) }} $$ よって$\displaystyle {{ f ( x_{n} ) } \over { f ' ( x_{n} ) }} - x_{n} = - x_{n+1}$だから $$ 0 = \alpha - x_{n+1} + {{( \alpha - x_{n})^2 } \over {2}} {{f’’ (\xi_{n})} \over { f '(x_{n}) }} $$ 項を整理すると $$ \begin{equation} \displaystyle \alpha - x_{n+1} = - {{f’’ (\xi_{n})} \over { 2 f '(x_{n}) }} ( \alpha - x_{n})^2 \end{equation} $$

---

第2部. 収束性

$f '$は$\alpha$の近くで連続であり、$f ' (\alpha) \ne 0$があるので、すべての$x \in I_{\delta}$に対して$f ' (x) \ne 0$を満たす$I_{\delta} : = [ \alpha - \delta , \alpha + \delta ]$が存在する。

$\displaystyle M( \delta ) := {{ \max_{x \in I_{\delta} } | f '' (x) | } \over { 2 \min_{x \in I_{\delta} } | f '(x) | }}$と置くと、$\delta$を小さくするたびに、$\displaystyle \max_{x \in I_{\delta} } | f '' (x) |$は増加せず、$\displaystyle \min_{x \in I_{\delta} } | f '(x) |$は減少しないので、$M ( \delta )$全体は増加しない。したがって、$M(\delta) | \alpha - x_{0} | < 1$が成り立つような十分に小さな$\delta = \delta_{0}$を選べる。まさにこの$x_{0}$が$\alpha$に十分近い初期値である。

この$\delta_{0}$について、今$M:= M(\delta_{0})$としよう。

$(1)$の両辺に絶対値をとると $$ | \alpha - x_{n+1}| \le M | \alpha - x_{n} |^2 $$ 両辺に$M$を掛けると $$ M | \alpha - x_{n+1}| \le M^2 | \alpha - x_{n} |^2 $$ 右辺を二乗でまとめると $$ M | \alpha - x_{n+1}| \le ( M | \alpha - x_{n} | )^2 $$ この過程を$x_{n+1}$から$x_{0}$が出るまで繰り返すと $$ M | \alpha - x_{n+1}| \le ( M | \alpha - x_{n} | )^2 \le ( M | \alpha - x_{n-1} | )^4 \le \cdots \le ( M | \alpha - x_{0} | )^{2^{n+1}} $$ 整理すると $$ | \alpha - x_{n} | \le {{1 } \over {M}} ( M | \alpha - x_{0} | )^{2^{n}} $$ $M | \alpha - x_{0} | < 1$なので、$n \to \infty$の時に$x_{n} \to \alpha$

---

第3部. クアドラティックに収束

$(1)$により、$\displaystyle {{ \alpha - x_{n+1} } \over { ( \alpha - x_{n})^2 }} = - {{f’’ (\xi_{n})} \over { 2 f '(x_{n}) }}$であり、$x_{n} \to \alpha$の時に$\xi_{n} \to \alpha$なので $$ \lim_{n \to \infty} {{\alpha - x_{n+1} } \over { ( \alpha - x_{n} )^2 }} = - \lim_{n \to \infty} {{f '' ( \xi_{n} )} \over { 2 f ' ( x_{n} ) }} = -{{f '' (\alpha)} \over { 2 f ' ( \alpha ) }} $$

■

説明で収束性が常に保証されるわけではないと話したが、実際に十分に小さな$I$を選ぶことが重要である。

実装

以下はRで書かれたコードだ。dfを導関数として入れれば、導関数を直接計算し、省略すれば微分係数を別に求めて使う。itmaxオプションは最大で繰り返す回数で、デフォルトで$1000$回繰り返しても十分な精度の解が得られなければ、計算を諦める。

d<-function(f,x,tol=10^(-8))
{
  h<-1
  d1<-1
  d2<-0
  
  while(abs(d1-d2)>tol)
  {
    d1<-d2
    d2<-(f(x+h)-f(x))/h
    h<-h/2
  }
  return(d2)
}
 
Newton<-function(f,x0,df=FALSE,itmax=10^3,tol=10^(-8))
{
  if(f(x0)==0){return(x0)}
  denom=0
  
  x1 = x0 - f(x0)/d(f,x0)
  
  for(i in 1:itmax)
  {
    if(is.function(df)){
      denom=df(x1)
    }else{
      denom=d(f,x1)
    }
    if(denom==0){stop('Zero Derivative\')}
    
    x2 = x1 - f(x1)/denom
    if(abs(x2-x1)<tol){
      return(x2)
    }else{
      x1 = x2
    }
  }
  
  stop('Maybe Wrong Initial Point')
}
 
f<-function(x) {x^3 + 8}
Newton(f,7)
 
g<-function(x) {x^6 - x - 1}
Newton(g,3)
 
h<-function(x) {exp(x)-1}
Newton(h,-2)