差分係数

数学では、微分を任意の精度に近似するために、有限差分を使用することができます。有限差分は、中心差分前方差分後方差分のいずれかです。

中心差分

この表には、中心差分の係数がいくつかの精度オーダーと均一なグリッド間隔で含まれています。[ 1 ]

デリバティブ 正確さ −5−4−3−2−1012345
1 2 −1/201/2
4 1/12−2/302/3−1/12
6 −1/603月20日−3/403/4−3/201/60
8 1/280−4/1051/5−4/504/5−1/54/105−1/280
2 2 1−21
4 −1/124/3−5/24/3−1/12
6 1/90−3/203/2−49/183/2−3/201/90
8 −1/5608/315−1/58/5−205/728/5−1/58/315−1/560
3 2 −1/210−11/2
4 1/8−18月13日0−13/81−1/8
6 −7/2403/10−169/12061/300−61/30169/120−3/107/240
4 2 1−46−41
4 −1/62−13/23月28日−13/22−1/6
6 7/240−2/5169/60−122/1591/8−122/15169/60−2/57/240
5 2 −1/22−5/205/2−21/2
4 1/6−3/213/3−29/606月29日−13/33/2−1/6
6 −13/28819/36−87/3213/2−323/480323/48−13/287/32−19/3613/288
6 2 1−615−2015−61
4 −1/43−1329−75/229−133−1/4
6 13/240−19/2487/16−39/2323/8−1023/20323/8−39/287/16−19/2413/240

例えば、2次の精度を持つ3次導関数は

f×012f×2+f×1f×+1+12f×+2h×3+h×2{\displaystyle f'''(x_{0})\approx {\frac {-{\frac {1}{2}}f(x_{-2})+f(x_{-1})-f(x_{+1})+{\frac {1}{2}}f(x_{+2})}{h_{x}^{3}}}+O\left(h_{x}^{2}\right),}

ここで、 は各有限差分区間間の均一なグリッド間隔を表し、 です。 h×{\displaystyle h_{x}}×n×0+nh×{\displaystyle x_{n}=x_{0}+nh_{x}}

精度 の 階微分には中心係数が存在する。これらは線形方程式系の解によって与えられる。 メートル{\displaystyle m}n{\displaystyle n}2p+12メートル+121+n{\displaystyle 2p+1=2\left\lfloor {\frac {m+1}{2}}\right\rfloor -1+n}1つのp1つのp+11つのp11つのp{\displaystyle a_{-p},a_{-p+1},...,a_{p-1},a_{p}}

1111pp+1p1pp2p+12p12p2p2pp+12pp12pp2p1つのp1つのp+11つのp+21つのp000メートル!0{\displaystyle {\begin{pmatrix}1&1&...&1&1\\-p&-p+1&...&p-1&p\\(-p)^{2}&(-p+1)^{2}&...&(p-1)^{2}&p^{2}\\...&...&...&...&...\\...&...&...&...&...\\...&...&...&...&...\\(-p)^{2p}&(-p+1)^{2p }&...&(p-1)^{2p}&p^{2p}\end{pmatrix}}{\begin{pmatrix}a_{-p}\\a_{-p+1}\\a_{-p+2}\\...\\...\\...\\a_{p}\end{pmatrix}}={\begin{pmatrix}0\\0\\0\\...\\m!\\...\\0\end{pmatrix}},}

ここで、右側の唯一のゼロでない値は- 行目にあります。 メートル+1{\displaystyle (m+1)}

1次元で任意の導関数と精度オーダーの有限差分係数を計算するオープンソース実装が利用可能です。[ 2 ] 左側の行列が転置されたヴァンデルモンド行列である場合、並べ替えると、係数は基本的に、点のウィンドウに - 次多項式を近似して導出することによって計算されることがわかります。したがって、係数は、多項式次数とウィンドウサイズがである完全に決定されたSavitzky–Golay フィルタの - 次導関数として計算することもできます。このために、オープンソース実装も利用できます。[ 3 ]係数の順序が異なる2つの定義が考えられます。離散畳み込みによるフィルタリング用フィルタ、または行列ベクトル積によるフィルタリング用フィルタです。上の表に示されている係数は後者の定義に対応しています。 JT{\displaystyle \mathbf {J} ^{T}}2p{\displaystyle 2p}2p+1{\displaystyle 2p+1}メートル{\displaystyle m}2p{\displaystyle 2p}2p+1{\displaystyle 2p+1}

ラグランジュ多項式理論は差分係数の明示的な公式を提供する。[ 4 ]最初の6つの導関数については次のようになる。

デリバティブ1つの0{\displaystyle a_{0}}1つのpp0{\displaystyle a_{p}(p\neq 0)}
10{\displaystyle 0}1!1p+1n!2pnp!n+p!{\displaystyle {\frac {1!(-1)^{p+1}(n!)^{2}}{p(n-p)!(n+p)!}}}
22Hn,2{\displaystyle -2H_{n,2}}2!(1)p+1(n!)2p2(np)!(n+p)!{\displaystyle {\frac {2!(-1)^{p+1}(n!)^{2}}{p^{2}(n-p)!(n+p)!}}}
30{\displaystyle 0}3!(1)p+1(n!)2p3(np)!(n+p)!(1p2Hn,2){\displaystyle {\frac {3!(-1)^{p+1}(n!)^{2}}{p^{3}(n-p)!(n+p)!}}(1-p^{2}H_{n,2})}
412(Hn,22Hn,4){\displaystyle 12(H_{n,2}^{2}-H_{n,4})}4!(1)p+1(n!)2p4(np)!(n+p)!(1p2Hn,2){\displaystyle {\frac {4!(-1)^{p+1}(n!)^{2}}{p^{4}(n-p)!(n+p)!}}(1-p^{2}H_{n,2})}
50{\displaystyle 0}5!(1)p+1(n!)2p5(np)!(n+p)!(1p2Hn,2+p42(Hn,22Hn,4)){\displaystyle {\frac {5!(-1)^{p+1}(n!)^{2}}{p^{5}(n-p)!(n+p)!}}\left(1-p^{2}H_{n,2}+{\frac {p^{4}}{2}}(H_{n,2}^{2}-H_{n,4})\right)}
6120Hn,23+360Hn,2Hn,4120Hn,6{\displaystyle -120H_{n,2}^{3}+360H_{n,2}H_{n,4}-120H_{n,6}}6!(1)p+1(n!)2p6(np)!(n+p)!(1p2Hn,2+p42(Hn,22Hn,4)){\displaystyle {\frac {6!(-1)^{p+1}(n!)^{2}}{p^{6}(n-p)!(n+p)!}}\left(1-p^{2}H_{n,2}+{\frac {p^{4}}{2}}(H_{n,2}^{2}-H_{n,4})\right)}

ここで、は一般化調和数です。 Hn,m{\displaystyle H_{n,m}}

前進有限差分

この表には、いくつかの精度のオーダーと均一なグリッド間隔での前方差分の係数が含まれています。 [ 1 ]

デリバティブ 正確さ 0 1 2 3 4 5 6 7 8
1 1−11       
2−3/22−1/2      
3−11/63−3/21/3     
4−25/124−34/3−1/4    
5−137/605−510月3日−5/41/5   
6−49/206−15/23月20日−15/46/5−1/6  
2 11−21      
22−54−1     
335/12−26/319/2−14/311月12日    
44月15日−77/6107/6−1361/12−5/6   
5203/45−87/5117/4−254/933/2−27/5137/180  
6469/90−223/10879/20−949/1841−201/101019/180−7/10 
3 1−13−31     
2−5/29−127−3/2    
3−17/471/4−59/249/2−41/47/4   
4−49/829−461/862−307/813−15/8  
5−967/120638/15−3929/40389/3−2545/24268/5−1849/12029/15 
6−801/80349/6−18353/1202391/10−1457/64891/30−561/8527/30−469/240
4 11−46−41    
23−1426−2411−2   
335/6−31137/2−242/3107/2−196月17日  
43月28日−111/2142−1219/6176−185/282/3−7/2 
51069/80−1316/1515289/60−2144/510993/24−4772/152803/20−536/15967/240

例えば、3次精度の1次導関数と2次精度の2次導関数は、

f(x0)116f(x0)+3f(x+1)32f(x+2)+13f(x+3)hx+O(hx3),{\displaystyle \displaystyle f'(x_{0})\approx \displaystyle {\frac {-{\frac {11}{6}}f(x_{0})+3f(x_{+1})-{\frac {3}{2}}f(x_{+2})+{\frac {1}{3}}f(x_{+3})}{h_{x}}}+O\left(h_{x}^{3}\right),}
f(x0)2f(x0)5f(x+1)+4f(x+2)f(x+3)hx2+O(hx2),{\displaystyle \displaystyle f''(x_{0})\approx \displaystyle {\frac {2f(x_{0})-5f(x_{+1})+4f(x_{+2})-f(x_{+3})}{h_{x}^{2}}}+O\left(h_{x}^{2}\right),}

対応する逆近似は次のように与えられる。

f(x0)116f(x0)3f(x1)+32f(x2)13f(x3)hx+O(hx3),{\displaystyle \displaystyle f'(x_{0})\approx \displaystyle {\frac {{\frac {11}{6}}f(x_{0})-3f(x_{-1})+{\frac {3}{2}}f(x_{-2})-{\frac {1}{3}}f(x_{-3})}{h_{x}}}+O\left(h_{x}^{3}\right),}
f(x0)2f(x0)5f(x1)+4f(x2)f(x3)hx2+O(hx2),{\displaystyle \displaystyle f''(x_{0})\approx \displaystyle {\frac {2f(x_{0})-5f(x_{-1})+4f(x_{-2})-f(x_{-3})}{h_{x}^{2}}}+O\left(h_{x}^{2}\right),}

後方差分

順方向近似の係数から逆方向近似の係数を得るには、前節の表に挙げたすべての奇数導関数に逆の符号を与え、偶数導関数には同じ符号を与える。次の表はこれを示している。[ 5 ]

デリバティブ 正確さ −8 −7 −6 −5 −4 −3 −2 −1 0
1 1       −11
2      1/2−23/2
3     −1/33/2−311月6日
2 1      1−21
2     −14−52
3 1     −13−31
2    3/2−712−95/2
4 1    1−46−41
2   −211−2426−143

任意のステンシルポイント

任意のステンシル点と、ステンシル点の数より1少ない次までの任意の微分に対して、差分係数は線形方程式を解くことによって得られる[ 6 ]。N{\displaystyle \displaystyle N}s{\displaystyle \displaystyle s}d<N{\displaystyle \displaystyle d<N}

(s10sN0s1N1sNN1)(a1aN)=d!(δ0,dδi,dδN1,d),{\displaystyle {\begin{pmatrix}s_{1}^{0}&\cdots &s_{N}^{0}\\\vdots &\ddots &\vdots \\s_{1}^{N-1}&\cdots &s_{N}^{N-1}\end{pmatrix}}{\begin{pmatrix}a_{1}\\\vdots \\a_{N}\end{pmatrix}}=d!{\begin{pmatrix}\delta _{0,d}\\\vdots \\\delta _{i,d}\\\vdots \\\delta _{N-1,d}\end{pmatrix}},}

ここで、 はクロネッカーのデルタであり、 の場合は 1 、それ以外の場合は 0 になります。 δi,j{\displaystyle \delta _{i,j}}i=j{\displaystyle i=j}

たとえば、 の場合、微分次数は次のようになります。 s=[3,2,1,0,1]{\displaystyle s=[-3,-2,-1,0,1]}d=4{\displaystyle d=4}

(a1a2a3a4a5)=(1111132101941012781018116101)1(000024)=(14641).{\displaystyle {\begin{pmatrix}a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\end{pmatrix}}={\begin{pmatrix}1&1&1&1&1\\-3&-2&-1&0&1\\9&4&1&0&1\\-27&-8&-1&0&1\\81&16&1&0&1\\\end{pmatrix}}^{-1}{\begin{pmatrix}0\\0\\0\\0\\24\end{pmatrix}}={\begin{pmatrix}1\\-4\\6\\-4\\1\end{pmatrix}}.}

近似の精度の順序は通常の形式(中心差分の場合はそれより良い形式)になります。 O(hx(Nd)){\displaystyle O\left(h_{x}^{(N-d)}\right)}

参照

参考文献

  1. ^ a b Fornberg, Bengt (1988)、「任意間隔グリッド上での有限差分式の生成」、Mathematics of Computation51 (184): 699– 706、doi : 10.1090/S0025-5718-1988-0935077-0ISSN  0025-5718
  2. ^ 「任意の次元数の有限差分数値導関数のためのPythonパッケージ」GitHub2021年10月14日。
  3. ^ "scipy.signal.savgol_filter" . Scipyオンラインドキュメント. 2008–2024.
  4. ^ 「有限差分係数」 . StackExchange . 2023年6月5日.
  5. ^ Taylor, Cameron (2019年12月12日). 「差分係数計算機」 . MIT.
  6. ^ 「差分係数計算機」