logo

数値解析学における階差段 📂数値解析

数値解析学における階差段

定義

関数 f:R→Rf : \mathbb{R} \to \mathbb{R} に対して異なる x1,⋯ ,xnx_{1} , \cdots , x_{n} における 区分差分divided Differenceは以下のように定義される。

f[x0]:=f(x0)f[x0,x1]:=f(x1)−f(x0)x1−x0f[x0,x1,x2]:=f[x1,x2]−f[x0,x1]x2−x0f[x0,⋯ ,xn]:=f[x1,⋯ ,xn]−f[x0,⋯ ,xn−1]xn−x0 \begin{align*} f[x_{0}] :=& f( x_{0} ) \\ 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} }} \\ f [ x_{0} , \cdots , x_{n} ] :=& {{ f [ x_{1} , \cdots , x_{n} ] - f [ x_{0} , \cdots , x_{n-1} ] } \over { x_{n} - x_{0} }} \end{align*}

定理

  • [1]: f[x0,⋯ ,xn]f [ x_{0} , \cdots , x_{n} ] は x0,⋯ ,xnx_{0} , \cdots , x_{n} の順序に関わらず常に同じである。
  • [2]: f[x0,⋯ ,xn]=∑i=0nf(xi)∏j∈{0,⋯ ,n}∖{i}(xi−xj)f [ x_{0} , \cdots , x_{n} ] = \sum_{i=0}^{n} {{ f( x_{i} ) } \over { \prod_{j \in \left\{ 0 , \cdots , n \right\} \setminus \left\{ i \right\} } ( x_{i} - x_{j} ) }}
  • [3]: 次を満たす ξ\xi が H{x0,x1}\mathscr{H} \left\{ x_{0} , x_{1} \right\} に存在する。f’(ξ)=f[x0,x1]f’ ( \xi ) = f [ x_0 , x_{1} ]
  • [4]: 次を満たす ξ\xi が H{x0,⋯ ,xn}\mathscr{H} \left\{ x_{0} , \cdots , x_{n} \right\} に存在する。1n!f(n)(ξ)=f[x0,⋯ ,xn]{{1} \over {n!}} f^{(n)} ( \xi ) =f [ x_{0} , \cdots , x_{n} ]
  • [3]’: f[xi,xi]=f′(xi)f [ x_{i} , x_{i} ] = f ' ( x_{i} )
  • [4]’: f[xi,⋯ ,xi⏟n+1]=1n!f(n)(xi)f [ \underbrace{ x_{i} , \cdots , x_{i} }_{ n+1 } ] = {{1} \over {n!}} f^{(n)} ( x_{i} )
  • H{a,b,c,⋯ }\mathscr{H} \left\{ a,b,c, \cdots \right\} は a,b,c,⋯a,b,c, \cdots を含む最小の区間を表す。

説明

区分差分は、数値解析の様々な理論を表現する際に空間を節約するのに大きな助けとなる。

定理 [2]は実際の計算において特に便利で、 f[x0,x1,x2]=f(x0)(x0−x1)(x0−x2)+f(x1)(x1−x2)(x1−x0)+f(x2)(x2−x0)(x2−x1) f [ x_{0} , x_{1} , x_{2} ] = {{ f( x_{0} ) } \over { ( x_{0} - x_{1} ) ( x_{0} - x_{2} ) }} + {{ f( x_{1} ) } \over { ( x_{1} - x_{2} ) ( x_{1} - x_{0} ) }} + {{ f( x_{2} ) } \over { ( x_{2} - x_{0} ) ( x_{2} - x_{1} ) }} のように表現されるため、実際の定義と同じように反復計算を省くことができる。

[3]は実質平均値の定理であり、微分した際に「瞬間変化率」の概念を考えたことがあれば、なぜ区分差分が微分を代替し、Divided Differenceという名前を持つのか理解できるだろう。定理 [4]はこのように一般化可能である。証明は区分差分の変数が無限に近づく、つまり極限を取ることを考えると、[3]’, [4]‘のようにnn階の微分係数と見ることができる。このような概念的アプローチは、区分差分の定義から来る限界を克服し、数値解析の理論をさらに豊かにする。

例

簡単にf(x):=x2+1f(x) := x^2 + 1 に対していくつかの区分差分を計算してみよう。

一点に対しては、 f[3]=f(3)=10 f[3] = f(3) = 10

二点に対しては、 f[0,5]=f(0)−f(5)0−5=1−26−5=5 f[0,5] = {{ f(0) - f(5) } \over { 0 - 5}} = {{1 - 26 } \over { -5 }} = 5

三点に対しては、 f[2,3,−1]=f(2)−f(3)2−3−f(3)−f(−1)3−(−1)2−(−1)=5−10−1−10−243=5−23=1 \begin{align*} f[2,3,-1] =& {{ \displaystyle {{ f(2) - f( 3) } \over { 2 - 3 }} - {{ f(3) - f(-1) } \over { 3 - (-1) }} } \over { 2 - (-1) }} \\ =& {{ \displaystyle {{ 5 - 10 } \over { -1 }} - {{ 10 - 2 } \over { 4 }} } \over { 3}} \\ =& {{ 5 -2 } \over {3}} \\ =& 1 \end{align*} 区分差分はもちろん有限ベクトルに対して常に計算できるが、通常は上述のようにn=2n=2 まで使用され、それ以外ではほとんど使用されない。

実装

以下は区分差分を R コードで実装したもので、関数 ff とベクトル x=(x0,⋯ ,xn)\mathbf{x} = ( x_{0} , \cdots , x_{n} ) を受け取り、f[x0,⋯ ,xn]f [ x_{0} , \cdots , x_{n} ] を返す。

f<-function(x) {x^2 + 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))
}
dd(f,c(3))
dd(f,c(0,5))
dd(f,c(2,3,-1))

上記コードを実行した結果は以下の通りであり、f(x)=x2+1f(x) = x^2 +1 に対する手計算と正確に一致することが確認できる。

20190319\_161740.png

証明

[2]

戦略: ラグランジュの公式とニュートンの区分差分公式の次数を比較する。

Ψn(x):=(x−x0)⋯(x−xn) \Psi_{n} (x) := (x-x_{0}) \cdots (x - x_{n}) Ψn\Psi_{n} をxx に対して微分し、x=xjx = x_{j} を代入すると Ψ’n(xj):=(x−x0)⋯(x−xj−1)(x−xj+1)⋯(x−xn) \Psi’_{n} (x_{j}) := (x-x_{0}) \cdots (x - x_{j-1})(x - x_{j+1}) \cdots (x - x_{n})

  ⟹  Ψ(x)n(x−xj)Ψ’n(xj)≡(x−x0)⋯(x−xj−1)(x−xj+1)⋯(x−xn)(xj−x0)⋯(xj−xj−1)(xj−xj+1)⋯(xj−xn)=lj(x) \implies {{ \Psi (x)_{n} } \over { (x - x_{j}) \Psi’_{n}(x_{j}) }} \equiv {{ (x - x_{0} ) \cdots (x - x_{j-1} ) (x - x_{j+1} ) \cdots (x - x_{n} ) } \over { (x_{j} - x_{0} ) \cdots (x_{j} - x_{j-1} ) (x_{j} - x_{j+1} ) \cdots (x_{j} - x_{n} ) }} = l_{j} (x) ラグランジュの公式によれば、ff のポリノミアルの補間は pn(x)=∑j=0nΨn(x)(x−xj)Ψ’n(xj)⋅f(xj) p_{n} (x) = \sum_{j=0}^{n} {{ \Psi_{n} (x) } \over { (x - x_{j}) \Psi’_{n} (x_{j}) }} \cdot f( x_{j}) したがって、pnp_{n} の最高次項の係数は ∑j=0nf(xj)Ψ’n(xj)\displaystyle \sum_{j=0}^{n} {{ f( x_{j}) } \over { \Psi’_{n} (x_{j}) }} であり、ニュートンの区分差分公式でのpnp_{n} の最高次項の係数は f[x0,⋯ ,xn]f [x_{0} , \cdots , x_{n} ] であるため、 f[x0,⋯ ,xn]=∑i=0nf(xi)∏j∈{0,⋯ ,n}∖{i}(xi−xj) f [ x_{0} , \cdots , x_{n} ] = \sum_{i=0}^{n} {{ f( x_{i} ) } \over { \displaystyle \prod_{j \in \left\{ 0 , \cdots , n \right\} \setminus \left\{ i \right\} } ( x_{i} - x_{j} ) }}

■

[1]

定理 [2]により、[x0,⋯ ,xn][x_{0} , \cdots , x_{n}] 内の順序を変えることは∑i=0nf(xi)∏j∈{0,⋯ ,n}∖{i}(xi−xj)\displaystyle \sum_{i=0}^{n} {{ f( x_{i} ) } \over { \displaystyle \prod_{j \in \left\{ 0 , \cdots , n \right\} \setminus \left\{ i \right\} } ( x_{i} - x_{j} ) }} を計算する際に加算の順序を変えることに過ぎず、したがって常に同じである。

■

[4]

戦略: ラグランジュの公式とニュートンの区分差分公式の次数を比較する。

pn(x)=pn−1(x)+C(x)p_{n}(x) = p_{n-1}(x) + C(x) とすると deg⁡C=n\deg C = n であり、i=0,⋯ ,(n−1)i= 0, \cdots, (n-1) に対して pn(xi)=f(xi)=pn−1(xi)p_{n}(x_{i}) = f(x_{i}) = p_{n-1}(x_{i}) であるため、ある定数 an≠0a_{n} \ne 0 に対してC(x)=an(x−x0)⋯(x−xn−1)C (x) = a_{n} (x - x_{0}) \cdots (x - x_{n-1}) のように表される。C(x)C(x) にx=xnx=x_{n} を代入してみると C(xn)=pn(xn)−pn−1(xn) C(x_{n} ) = p_{n} (x_{n}) - p_{n-1}(x_{n})

C(xn)=an(xn−x0)⋯(xn−xn) C (x_{n}) = a_{n} (x_{n} - x_{0}) \cdots (x_{n} - x_{n}) pnp_{n} はff のポリノミアルの補間であるため、pn(xn)=f(xn)p_{n} (x_{n}) = f(x_{n}) である an(xn−x0)⋯(xn−xn−1)=f(xn)−pn−1(xn) a_{n} (x_{n} - x_{0}) \cdots (x_{n} - x_{n-1}) = f (x_{n}) - p_{n-1}(x_{n})

an=f(xn)−pn−1(xn−1)(xn−x0)⋯(xn−xn−1) a_{n} = {{ f (x_{n}) - p_{n-1}(x_{n-1}) } \over { (x_{n} - x_{0}) \cdots (x_{n} - x_{n-1}) }}

ポリノミアル補間と実際の関数との誤差: nn回微分可能なf:R→Rf : \mathbb{R} \to \mathbb{R} とあるξ∈H{x0,⋯ ,xn−1}\xi \in \mathscr{H} \left\{ x_{0} , \cdots , x_{n-1} \right\} に対して、ff のポリノミアル補間 pn−1p_{n-1} はあるxn∈Rx_{n} \in \mathbb{R} に対して f(xn)−pn−1(xn)=(xn−x0)⋯(xn−xn−1)n!f(n)(ξ) f(x_{n}) - p_{n-1} (x_{n}) = {{ (x_{n} - x_{0}) \cdots (x_{n} - x_{n-1}) } \over { n! }} f^{(n)} ( \xi )

ポリノミアル補間と実際の関数との誤差の式に従って、 an=f(n)(ξ)n! a_{n} = {{ f^{(n)} ( \xi ) } \over { n! }} 一方でニュートンの区分差分公式でのpnp_{n} の最高次項の係数はan=f[x0,⋯ ,xn]a_{n} = f [x_{0} , \cdots , x_{n} ] であるため、 1n!f(n)(ξ)=f[x0,⋯ ,xn] {{1} \over {n!}} f^{(n)} ( \xi ) =f [ x_{0} , \cdots , x_{n} ]

■

[3]

定理 [4]により自明である。

■