logo

ニュートン-ラプソン法 📂数値解析

ニュートン-ラプソン法

メソッド 1

20180831_192324.png

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

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

説明

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

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

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

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

一緒に見る

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

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

証明

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

第1部. テイラー展開

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

第2部. 収束性

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

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

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

(1)(1)の両辺に絶対値をとると ∣α−xn+1∣≤M∣α−xn∣2 | \alpha - x_{n+1}| \le M | \alpha - x_{n} |^2 両辺にMMを掛けると M∣α−xn+1∣≤M2∣α−xn∣2 M | \alpha - x_{n+1}| \le M^2 | \alpha - x_{n} |^2 右辺を二乗でまとめると M∣α−xn+1∣≤(M∣α−xn∣)2 M | \alpha - x_{n+1}| \le ( M | \alpha - x_{n} | )^2 この過程をxn+1x_{n+1}からx0x_{0}が出るまで繰り返すと M∣α−xn+1∣≤(M∣α−xn∣)2≤(M∣α−xn−1∣)4≤⋯≤(M∣α−x0∣)2n+1 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}} 整理すると ∣α−xn∣≤1M(M∣α−x0∣)2n | \alpha - x_{n} | \le {{1 } \over {M}} ( M | \alpha - x_{0} | )^{2^{n}} M∣α−x0∣<1M | \alpha - x_{0} | < 1なので、n→∞n \to \inftyの時にxn→αx_{n} \to \alpha

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

(1)(1)により、α−xn+1(α−xn)2=−f’’(ξn)2f′(xn)\displaystyle {{ \alpha - x_{n+1} } \over { ( \alpha - x_{n})^2 }} = - {{f’’ (\xi_{n})} \over { 2 f '(x_{n}) }}であり、xn→αx_{n} \to \alphaの時にξn→α\xi_{n} \to \alphaなので lim⁡n→∞α−xn+1(α−xn)2=−lim⁡n→∞f′′(ξn)2f′(xn)=−f′′(α)2f′(α) \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 ) }}

■

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

実装

20180831_193141.png

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

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)

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