ニュートン・ラフソン法 📂数値解析

ニュートン・ラフソン法

メソッド¹

$20180831\_192324.png$

$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$の近似解を求めるのに有用な方法である。

クアドラティック^quadraticに収束するという言葉は、収束率を語るときに$2$次で収束するという意味だが、適切な訳語がない。$1$次が「線形」と訳されるのに比べれば、単に$2$次で収束するという言い方はあまりにも漠然としていて使用例も少ないため、やむを得ず英語表現をそのまま使った。それでもあえて言及したのは、ほとんどの場合二分法に比べて収束する速度が速いことを数学的に示せるからである。

$20180831\_192719.png$ しかし、「$\alpha$の近傍」や「十分に近い初期値」という表現から察せられるように、実際には何も考えずに使える方法ではない。上のように運悪く初期値を間違えて与えて$0$で割る問題が起きたり、数列自体が解から遠いところへ行ってしまうこともある。そのため、コードを書くときも、特定の回数以内に解が求められなければ計算を諦めるなどの例外処理が必要である。

このような発散エラーは、数学的に収束性が証明されていても、実際に使うときには十分起こりうる問題となる。最も重要なユーザーが「$\alpha$の近傍」や「十分に近い初期値」を正確に把握できないかもしれないからである。二分法が遅かろうがどうであろうが必ず解を求めてくれるのと比べれば、明らかな短所だと言える。

証明

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

Part 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} $$

Part 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$が成り立つような$x_{0} \in I_{\delta}$が存在するように、十分小さい$\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$

Part 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$を取っている。

実装

$20180831\_193141.png$

以下は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)

Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p58. ↩︎

ニュートン・ラフソン法

メソッド1

説明

関連リンク

証明

実装

メソッド¹