logo

セカント法 📂数値解析

セカント法

メソッド

20180903\_131323.png 20180903\_131333.png 20180903\_131343.png

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

α\alpha に十分近い初期値 x0,x1x_{0} , x_{1} に対して xn+1:=xn−f(xn)xn−xn−1f(xn)−f(xn−1) x_{n+1} := x_{n} - f ( x_{n} ) {{ x_{n} - x_{n-1} } \over { f ( x_{n} ) - f ( x_{n-1} ) }} と定義された数列 {xn}\left\{ x_{n} \right\} はn→∞n \to \infty の時 α\alpha に 1+52\displaystyle {{1 + \sqrt{5} } \over {2}} 次で収束する。

説明

黄金比

収束次数はかなり馴染み深いものがあるはずだ。それはまさに黄金比の 1+52=1.618⋯\displaystyle {{1 + \sqrt{5} } \over {2}} = 1.618 \cdots であり、数列の定義から3つの連続した項を使用して証明中にフィボナッチ数列が出現する。セカントメソッドの有用性を言及する前にかなり興味深い現象である。この次数は、セカント法がバイセクションメソッドよりも早いが、ニュートン・ラフソン法よりは遅いことを意味する。

ニュートン・ラフソン法がより早いにも関わらず、このような方法を考案した理由は、導関数を求めることが過度に難しい場合があるからである。さらに、計算量においては収束次数が低いため、ニュートン・ラフソン法に比べて計算回数は多いが、一回一回の計算コストの面ではセカントメソッドが良い場合がある。どの方法がより良いかは問題による。（しかし、正直に言ってニュートン法以外はほとんど使わない。）

初期値問題

20180903\_131820.png

ニュートン・ラフソン法と同様に、初期値が間違って与えられた場合、十分に α\alpha に近くなければ収束しないことがある。

証明 1

  • H{a,b,c,⋯ }\mathscr{H} \left\{ a,b,c, \cdots \right\} は a,b,c,⋯a,b,c, \cdots を含む最小の区間を表示する。

パート1. 差分表

f[x0,x1]:=f(x1)−f(x0)x1−x0f[x0,x1,x2]:=f[x1,x2]−f[x0,x1]x2−x0 \begin{align*} f [ x_{0} , x_{1} ] :=& {{ f ( x_{1} ) - f ( x_{0} ) } \over { x_{1} - x_{0} }} \\ f [ x_{0} , x_{1} , x_{2} ] :=& {{ f [ x_{1} , x_{2} ] - f [ x_{0} , x_{1} ] } \over { x_{2} - x_{0} }} \end{align*} を定義する。これらを差分表と呼び、f’(ξ)=f[x0,x1]f ’ ( \xi ) = f [ x_{0} , x_{1} ] と 12f’’(ζ)=f[x0,x1,x2]\displaystyle {{1} \over {2}} f ’’ ( \zeta ) = f [ x_{0} , x_{1} , x_{2} ] を満たす ξ∈H{x0,x1}ζ∈H{x0,x1,x2} \begin{align*} \xi \in& \mathscr{H} \left\{ x_{0} , x_{1} \right\} \\ \zeta \in& \mathscr{H} \left\{ x_{0} , x_{1} , x_{2} \right\} \end{align*} が存在する。

パート2. α−xn+1=−(α−xn−1)(α−xn)f′′(ζn)2f’(ξn)\displaystyle \alpha - x_{n+1} = - ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) {{ f ''( \zeta _{n} ) } \over { 2 f ’ ( \xi_{n} ) }}

xn+1=xn−f(xn)xn−xn−1f(xn)−f(xn−1)x_{n+1} = x_{n} - f ( x_{n} ) {{ x_{n} - x_{n-1} } \over { f ( x_{n} ) - f ( x_{n-1} ) }} の両辺に −1-1 を掛けて α\alpha を加えると α−xn+1=α−xn+f(xn)xn−xn−1f(xn)−f(xn−1) \alpha - x_{n+1} = \alpha - x_{n} + f ( x_{n} ) {{ x_{n} - x_{n-1} } \over { f ( x_{n} ) - f ( x_{n-1} ) }} であり、 α−xn+1=α−xn+xnf(xn)−xn−1f(xn)f(xn)−f(xn−1)=xnf(xn)−xn−1f(xn)+αf(xn)−xnf(xn)−αf(xn−1)+xnf(xn−1)f(xn)−f(xn−1)=1f(xn)−f(xn−1)[−xn−1f(xn)+αf(xn)−αf(xn−1)+xnf(xn−1)]=1f(xn)−f(xn−1)[(α−xn−1)f(xn)−(α−xn)f(xn−1)]=−(α−xn−1)(α−xn)f(xn)−f(xn−1)[f(xn)xn−α−f(xn−1)]=−(α−xn−1)(α−xn)f(xn)−f(xn−1)[f(xn)−f(α)xn−α−f(xn−1)−f(α)xn−1−α]=−(α−xn−1)(α−xn)1f(xn)−f(xn−1)[f[xn,α]−f[xn−1,α]]=−(α−xn−1)(α−xn)xn−xn−1f(xn)−f(xn−1)[f[xn,α]−f[xn−1,α]xn−xn−1]=−(α−xn−1)(α−xn)f[xn−1,xn,α]f[xn−1,xn] \begin{align*} \alpha - x_{n+1} =& \alpha - x_{n} + {{ x_{n} f ( x_{n} ) - x_{n-1} f ( x_{n} ) } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \\ =& {{ x_{n} f ( x_{n} ) - x_{n-1} f ( x_{n} ) + \alpha f(x_{n}) - x_{n} f(x_{n}) - \alpha f(x_{n-1}) + x_{n} f(x_{n-1}) } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \\ =& {{1 } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \left[ - x_{n-1} f ( x_{n} ) + \alpha f(x_{n}) - \alpha f(x_{n-1}) + x_{n} f(x_{n-1} ) \right] \\ =& {{1 } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \left[ ( \alpha - x_{n-1} ) f ( x_{n} ) - ( \alpha - x_{n} ) f(x_{n-1} ) \right] \\ =& - {{ ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \left[ {{ f ( x_{n} ) } \over { x_{n} - \alpha }} - {{ f(x_{n-1}) }} \right] \\ =& - {{ ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \left[ {{ f ( x_{n} ) - f(\alpha) } \over { x_{n} - \alpha }} - {{ f(x_{n-1} ) - f( \alpha) } \over { x_{n-1} - \alpha }} \right] \\ =& - ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) {{ 1 } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \left[ f[ x_{n} , \alpha ] - f[ x_{n-1} , \alpha ] \right] \\ =& - ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) {{ x_{n} - x_{n-1} } \over { f ( x_{n} ) - f ( x_{n-1} ) }} \left[ {{ f[ x_{n} , \alpha ] - f[ x_{n-1} , \alpha ] } \over { x_{n} - x_{n-1} }} \right] \\ =& - ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) {{ f [ x_{n-1} , x_{n} , \alpha ] } \over { f [ x_{n-1} , x_{n} ] }} \end{align*} である。整理すると α−xn+1=−(α−xn−1)(α−xn)f[xn−1,xn,α]f[xn−1,xn] \alpha - x_{n+1} = - ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) {{ f [ x_{n-1} , x_{n} , \alpha ] } \over { f [ x_{n-1} , x_{n} ] }} であり、差分表の性質により次のことが得られる。 α−xn+1=−(α−xn−1)(α−xn)f′′(ζn)2f’(ξn) \alpha - x_{n+1} = - ( \alpha - x_{n-1} ) ( \alpha - x_{n} ) {{ f ''( \zeta _{n} ) } \over { 2 f ’ ( \xi_{n} ) }}

パート3. 収束性

全ての与えられた仮定を満たす α\alpha の近傍で十分に小さい I:=[α−δ,α+δ]I : = [ \alpha - \delta , \alpha + \delta ] を考える。

M:=max⁡x∈I∣f′′(x)∣2min⁡x∈I∣f′(x)∣ M := {{ \max_{x \in I } | f '' (x) | } \over { 2 \min_{x \in I} | f '(x) | }} と定義すると、パート1で得られた結果により ∣α−xn+1∣≤∣α−xn∣∣α−xn−1∣M | \alpha - x_{n+1}| \le | \alpha - x_{n} | | \alpha - x_{n-1} | M 両辺に MM を掛けると ∣α−xn+1∣M≤(∣α−xn∣M)⋅(∣α−xn−1∣M) | \alpha - x_{n+1}| M \le ( | \alpha - x_{n} | M ) \cdot ( | \alpha - x_{n-1} | M ) 十分に小さい II を仮定したなら、 ε:=max⁡{(∣α−xn∣M),(∣α−xn−1∣M)}<1 \varepsilon := \max \left\{ ( | \alpha - x_{n} | M ) , ( | \alpha - x_{n-1} | M ) \right\} < 1 を設定することで ∣α−xn+1∣M≤ε2 | \alpha - x_{n+1}| M \le \varepsilon^2 ∣α−xn+3∣M≤(∣α−xn+2∣M)⋅(∣α−xn+1∣M)≤ε2⋅ε=ε3 | \alpha - x_{n+3}| M \le ( | \alpha - x_{n+2}| M ) \cdot ( | \alpha - x_{n+1}| M ) \le \varepsilon^2 \cdot \varepsilon = \varepsilon^3 ∣α−xn+4∣M≤(∣α−xn+3∣M)⋅(∣α−xn+2∣M)≤ε3⋅ε2=ε5 | \alpha - x_{n+4}| M \le ( | \alpha - x_{n+3}| M ) \cdot ( | \alpha - x_{n+2}| M ) \le \varepsilon^3 \cdot \varepsilon^2 = \varepsilon^5 このような方法で不等式 ∣α−xn+m+1∣M≤(∣α−xn+m∣M)⋅(∣α−xn+m−1∣M)≤εqm⋅εqm−1=εqm+1 | \alpha - x_{n+m+1}| M \le ( | \alpha - x_{n+m}| M ) \cdot ( | \alpha - x_{n+m-1}| M ) \le \varepsilon^{ q_{m} } \cdot \varepsilon^{ q_{m-1} } = \varepsilon^{ q_{m+1} } を得る。

フィボナッチ数列の一般項: 数列 {Fn}\left\{ F_{n} \right\} が Fn+1:=Fn+Fn−1F_{n+1} := F_{n} + F_{n-1} として定義されているとする。F0=F1=1F_{0} = F_{1} = 1 の場合、r0:=1+52\displaystyle r_{0} : = {{1 + \sqrt{5} } \over {2}} と r1:=1−52\displaystyle r_{1} : = {{1 - \sqrt{5} } \over {2}} に対して Fn=r0n+1−r1n+1r0−r1\displaystyle F_{n} = {{ {r_{0}}^{n+1} - {r_{1}}^{n+1} } \over { r_{0} - r_{1} }}

{qm}\left\{ q_{m} \right\} は他ならぬフィボナッチ数列で、十分に大きな mm に対して qm∼15(1.618)m+1 q_{m} \sim {{1} \over {\sqrt{5}}} ( 1.618 )^{m+1} である。したがって ∣α−xn+m∣≤1Mεqm | \alpha - x_{n+m}| \le {{1} \over {M}} \varepsilon^{ q_{m} } であり、n→∞n \to \infty の時 xn→αx_{n} \to \alpha であり、さらには 1+52∼1.618 {{1 + \sqrt{5} } \over {2}} \sim 1.618 次で収束することが確認できる。

■

実装

以下は Rで記述されたコードである。

itmax オプションは、最大反復回数であり、デフォルトでは 10001000 回反復しても十分に精度の高い解が得られない場合、計算を諦める。

Secant<-function(f,x0,x1,itmax=10^3,tol=10^(-8))
{
  if(f(x0)==f(x1)){stop('Wrong Initial Values')}
  
  x2 = x1 - f(x1) \ast ( ( x1  - x0 ) / ( f(x1) - f(x0) ) )
  for(i in 1:itmax)
  {
    x0 = x1
    x1 = x2
    x2 = x1 - f(x1) \ast ( ( x1  - x0 ) / ( f(x1) - f(x0) ) )
    if(abs(x2-x1)<tol) {return(x2)}
  }
  
  stop('Maybe Wrong Initial Point')
}
 
f<-function(x) {x^3 + 8}
Secant(f,-7,7)
 
g<-function(x) {x^6 - x - 1}
Secant(g,0,3)
 
h<-function(x) {exp(x)-1}
Secant(h,-2,-1)

上記のコードを実行した結果は以下の通りである。

20180903\_125940.png

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