logo

複数の点を使用した有限差分の導出 📂数値解析

複数の点を使用した有限差分の導出

定理

NN番の微分可能関数 f:RRf : \mathbb{R} \to \mathbb{R} が与えられているとする。ある点 tRt \in \mathbb{R} での d<Nd < N回目の微分の関数値 f(d)(t)f^{(d)} \left( t \right) は、十分に小さい h>0h > 0 に対して、基数がnNn \in \mathbb{N}有限集合 SZS \subset \mathbb{Z} と同じ数の点 {t+kh:kS}\left\{ t + k h : k \in S \right\} を使って、以下のように近似できる。 f(d)(t)kSckf(t+kh) f^{(d)} \left( t \right) \approx \sum_{k \in S} c_{k} f \left( t + kh \right) ここで、{ck1,,ckn}\left\{ c_{k_{1}} , \cdots , c_{k_{n}} \right\} は次のように求められる。 [ck1ckn]=[k10kn0k1N1knN1]1d!hd[δ0,dδN1,d] \begin{bmatrix} c_{k_{1}} \\ \vdots \\ c_{k_{n}} \end{bmatrix} = \begin{bmatrix} k_{1}^{0} & \cdots & k_{n}^{0} \\ \vdots & \ddots & \vdots \\ k_{1}^{N-1} & \cdots & k_{n}^{N-1} \end{bmatrix}^{-1} {{ d! } \over { h^{d} }} \begin{bmatrix} \delta_{0,d} \\ \vdots \\ \delta_{N-1,d} \end{bmatrix} ここで、点 {t+kh:kS}\left\{ t + k h : k \in S \right\}ステンシルポイントと呼ぶこともある。


  • X1X^{-1}逆行列である。
  • x!x!階乗である。
  • δij={1,if i=j0,if ij\delta_{ij} = \begin{cases} 1 & , \text{if } i = j \\ 0 & , \text{if } i \ne j \end{cases}クロネッカーのデルタである。

説明

f(t)f(t+h)f(t)h f ' (t) \approx {{ f(t+h) - f(t) } \over { h }} 離散化されたデータしかなく、実際の関数がないか、あるいは微分が困難な場合は、上記のように平均変化率を参照するしかない。しかし、密接に連結した二点の情報を使用するのではなく、次のように左右の二点を使用することもできる。 f(t)f(t+h)f(th)2h f ' (t) \approx {{ f(t+h) - f(t-h) } \over { 2h }} 一般的に、より多くの点を使用すると精度が上がり、データが悪ければ悪いほど、このような計算に依存する状況が生じるかもしれない。

証明 1

戦略:元の文書と同様に、例のように簡単に示されているだけである。テイラーの定理を通じて連立方程式を立てると終わりだ。


例えば、f(4)(t)f^{(4)} (t) を計算するために、次のように5つの点を使用すると考えてみよう。私たちの目標は、次のようにうまく近似する{c2,,c2}\left\{ c_{-2} , \cdots , c_{2} \right\} を見つけることである。 f(4)(t)c2f(t2h)+c1f(t1h)+c0f(t)+c1f(t+h)+ck2f(t+2h) f^{(4)} (t) \approx c_{-2} f (t - 2h) + c_{-1} f (t - 1h) + c_{0} f (t) + c_{1} f (t + h) + c_{k_{2}} f (t + 2h) ここで、各f(t+kh)f (t + kh)は、テイラーの定理に従って次のように展開される。 f(4)(t)=c2[f(t)+f(t)(2h)+f(t)2(2h)2+f(3)(t)6(2h)3+f(4)(t)24(2h)4]+c1[f(t)+f(t)(h)+f(t)2(h)2+f(3)(t)6(h)3+f(4)(t)24(h)4]+c0f(t)+c1[f(t)+f(t)h+f(t)2h2+f(3)(t)6h3+f(4)(t)24h4]+c2[f(t)+f(t)(2h)+f(t)2(2h)2+f(3)(t)6(2h)3+f(4)(t)24(2h)4]+O(h5)=[c2+c1+c0+c1+c2]f(t)+[2c2c1+0c0+c1+2c2]hf(t)+[4c2+c1+0c0+c1+4c2]h22f(t)+[8c2c1+0c0+c1+8c2]h36f(3)(t)+[16c2+c1+0c0+c1+16c2]h424f(4)(t)+O(h5) \begin{align*} & f^{(4)} (t) \\ =& c_{-2} \left[ f(t) + f ' (t) (-2h) + {{ f '' (t) } \over { 2 }} (-2h)^{2} + {{ f^{(3)} (t) } \over { 6 }} (-2h)^{3} + {{ f^{(4)} (t) } \over { 24 }} (-2h)^{4} \right] \\ & + c_{-1} \left[ f(t) + f ' (t) (-h) + {{ f '' (t) } \over { 2 }} (-h)^{2} + {{ f^{(3)} (t) } \over { 6 }} (-h)^{3} + {{ f^{(4)} (t) } \over { 24 }} (-h)^{4} \right] \\ & + c_{0} f (t) \\ & + c_{1} \left[ f(t) + f ' (t) h + {{ f '' (t) } \over { 2 }} h^{2} + {{ f^{(3)} (t) } \over { 6 }} h^{3} + {{ f^{(4)} (t) } \over { 24 }} h^{4} \right] \\ & + c_{2} \left[ f(t) + f ' (t) (2h) + {{ f '' (t) } \over { 2 }} (2h)^{2} + {{ f^{(3)} (t) } \over { 6 }} (2h)^{3} + {{ f^{(4)} (t) } \over { 24 }} (2h)^{4} \right] \\ & + O \left( h^{5} \right) \\ =& \left[ c_{-2} + c_{-1} + c_{0} + c_{1} + c_{2} \right] f(t) \\ & + \left[ -2 c_{-2} - c_{-1} + 0 c_{0} + c_{1} + 2 c_{2} \right] h f ' (t) \\ & + \left[ 4 c_{-2} + c_{-1} + 0 c_{0} + c_{1} + 4 c_{2} \right] {{ h^{2} } \over { 2 }} f '' (t) \\ & + \left[ -8 c_{-2} - c_{-1} + 0 c_{0} + c_{1} + 8 c_{2} \right] {{ h^{3} } \over { 6 }} f^{(3)} (t) \\ & + \left[ 16 c_{-2} + c_{-1} + 0 c_{0} + c_{1} + 16 c_{2} \right] {{ h^{4} } \over { 24 }} f^{(4)} (t) \\ & + O \left( h^{5} \right) \end{align*} いつものように、O(h5)O \left( h^{5} \right)は十分に小さくして捨てる。そしてf(t)f(t)f(t)f’(t)f(t)f^{''} (t)f(3)(t)f^{(3)} (t)f(4)(t)f^{(4)} (t) が何であれ、上の近似式が成り立つには、 [c2+c1+c0+c1+c2]=0[2c2c1+0c0+c1+2c2]=0[4c2+c1+0c0+c1+4c2]=0[8c2c1+0c0+c1+8c2]=0[16c2+c1+0c0+c1+16c2]h4f(4)(t)=24f(4)(t) \begin{align*} \left[ c_{-2} + c_{-1} + c_{0} + c_{1} + c_{2} \right] =& 0 \\ \left[ -2 c_{-2} - c_{-1} + 0 c_{0} + c_{1} + 2 c_{2} \right] =& 0 \\ \left[ 4 c_{-2} + c_{-1} + 0 c_{0} + c_{1} + 4 c_{2} \right] =& 0 \\ \left[ -8 c_{-2} - c_{-1} + 0 c_{0} + c_{1} + 8 c_{2} \right] =& 0 \\ \left[ 16 c_{-2} + c_{-1} + 0 c_{0} + c_{1} + 16 c_{2} \right] h^{4} f^{(4)} (t) =& 24 f^{(4)} (t) \end{align*} を満たせばよい。行列の形に簡単に変えてみると、 [111112101241014810181610116][c2c1c0c1c2]=[(2)0(1)0(0)0(1)0(2)0(2)1(1)1(0)1(1)1(2)1(2)2(1)2(0)2(1)2(2)2(2)3(1)3(0)3(1)3(2)3(2)4(1)4(0)4(1)4(2)4][c2c1c0c1c2]=1h4[00004!] \begin{align*} & \begin{bmatrix} 1 & 1 & 1 & 1 & 1 \\ -2 & -1 & 0 & 1 & 2 \\ 4 & 1 & 0 & 1 & 4 \\ -8 & -1 & 0 & 1 & 8 \\ 16 & 1 & 0 & 1 & 16 \end{bmatrix} \begin{bmatrix} c_{-2} \\ c_{-1} \\ c_{0} \\ c_{1} \\ c_{2} \end{bmatrix} \\ =& \begin{bmatrix} (-2)^{0} & (-1)^{0} & (0)^{0} & (1)^{0} & (2)^{0} \\ (-2)^{1} & (-1)^{1} & (0)^{1} & (1)^{1} & (2)^{1} \\ (-2)^{2} & (-1)^{2} & (0)^{2} & (1)^{2} & (2)^{2} \\ (-2)^{3} & (-1)^{3} & (0)^{3} & (1)^{3} & (2)^{3} \\ (-2)^{4} & (-1)^{4} & (0)^{4} & (1)^{4} & (2)^{4} \end{bmatrix} \begin{bmatrix} c_{-2} \\ c_{-1} \\ c_{0} \\ c_{1} \\ c_{2} \end{bmatrix} = {{ 1 } \over { h^{4} }} \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 4! \end{bmatrix} \end{align*} 最後の行で、左辺は定理で見た通りの形で、右辺では、私たちが興味を持つdd番目の行を除き、すべて00 になってしまったことを確認できる。左辺の行列を右辺に移すと、私たちが欲しかった近似式を得る。

例えば、{3,1,0,1}\left\{ -3, -1, 0, 1 \right\}の4つの点を使用して2次の微分を近似したい場合は、以下の行列方程式を使用すればよい。 [(3)0(1)0(0)0(1)0(3)1(1)1(0)1(1)1(3)2(1)2(0)2(1)2(3)3(1)3(0)3(1)3(3)4(1)4(0)4(1)4][c3c1c0c1]=[1111310191012710181101][c3c1c0c1]=1h2[002!0] \begin{align*} & \begin{bmatrix} (-3)^{0} & (-1)^{0} & (0)^{0} & (1)^{0} \\ (-3)^{1} & (-1)^{1} & (0)^{1} & (1)^{1} \\ (-3)^{2} & (-1)^{2} & (0)^{2} & (1)^{2} \\ (-3)^{3} & (-1)^{3} & (0)^{3} & (1)^{3} \\ (-3)^{4} & (-1)^{4} & (0)^{4} & (1)^{4} \end{bmatrix} \begin{bmatrix} c_{-3} \\ c_{-1} \\ c_{0} \\ c_{1} \end{bmatrix} \\ =& \begin{bmatrix} 1 & 1 & 1 & 1 \\ -3 & -1 & 0 & 1 \\ 9 & 1 & 0 & 1 \\ -27 & -1 & 0 & 1 \\ 81 & 1 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_{-3} \\ c_{-1} \\ c_{0} \\ c_{1} \end{bmatrix} = {{ 1 } \over { h^{2} }} \begin{bmatrix} 0 \\ 0 \\ 2! \\ 0 \end{bmatrix} \end{align*} そして、以下を得るように解く。直接計算はとても難しいので、コンピューターの力を借りよう。 f(t)[c3f(t3h)+c1f(th)+c0f(t)+c1f(t+h)]=1h2[0f(t3h)+f(th)2f(t)+f(t+h)] \begin{align*} f '' (t) \approx & \left[ c_{-3} f (t - 3h) + c_{-1} f (t - h) + c_{0} f(t) + c_{1} f (t + h) \right] \\ =& {{ 1 } \over { h^{2} }} \left[ 0 f (t - 3h) + f (t - h) - 2 f(t) + f (t + h) \right] \end{align*}

精度

じっくり見せていないが、その精度はO(hnd)O \left( h^{n - d} \right)とされていて、ステンシルポイントを多く使用するほど精確になり、微分の数を増やすほど、コストが高くなると言われている。

テーブル 2

関連項目


  1. https://web.media.mit.edu/~crtaylor/calculator.html ↩︎

  2. Fornberg. (1988). Generation of Finite Difference Formulas on Arbitrarily Spaced Grids: Ch4. ↩︎