Class of numerical techniques
数値解析 において 、 有限差分法 （ FDM ）は、 導関数を 有限差分 で 近似することで 微分方程式 を解く数値解析手法の一種です 。空間領域と時間領域（該当する場合）の両方が 離散化 、つまり有限個の区間に分割され、区間の端点における解の値は、有限差分と近傍点の値を含む代数方程式を解くことで近似されます。
有限差分法は、 非線形である可能性のある 常微分方程式 （ODE）または 偏微分方程式 （PDE）を、 行列代数 技法によって解ける 線形方程式系 に変換します 。現代のコンピュータはこれらの 線形代数 計算を効率的に実行でき、実装の容易さも相まって、現代の数値解析においてFDMが広く使用されるようになりました。 [1]今日、FDMは 有限要素法 と並んで、PDEの数値解を求める最も一般的な手法の一つです 。 [1]
テイラー多項式から差商を導出する n 回微分可能な関数 の場合、 テイラーの定理 によりテイラー 級数 展開は次のように与えられる。 f ( x 0 + h ) = f ( x 0 ) + f ′ ( x 0 ) 1 ! h + f ( 2 ) ( x 0 ) 2 ! h 2 + ⋯ + f ( n ) ( x 0 ) n ! h n + R n ( x ) , {\displaystyle f(x_{0}+h)=f(x_{0})+{\frac {f'(x_{0})}{1!}}h+{\frac {f^{(2)}(x_{0})}{2!}}h^{2}+\cdots +{\frac {f^{(n)}(x_{0})}{n!}}h^{n}+R_{n}(x),}
ここで 、n ! は n の 階乗 を表し、 R n ( x ) は剰余項であり、 n 次テイラー多項式 と元の関数の差を表します。
以下は、テイラー多項式を切り捨てて
関数 f の1次導関数の近似値を導出する手順です。h で割ると次のように なります。 を 解くと ： f ( x 0 + h ) = f ( x 0 ) + f ′ ( x 0 ) h + R 1 ( x ) . {\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x).} f ( x 0 + h ) h = f ( x 0 ) h + f ′ ( x 0 ) + R 1 ( x ) h {\displaystyle {f(x_{0}+h) \over h}={f(x_{0}) \over h}+f'(x_{0})+{R_{1}(x) \over h}} f ′ ( x 0 ) {\displaystyle f'(x_{0})} f ′ ( x 0 ) = f ( x 0 + h ) − f ( x 0 ) h − R 1 ( x ) h . {\displaystyle f'(x_{0})={f(x_{0}+h)-f(x_{0}) \over h}-{R_{1}(x) \over h}.}
が十分に小さいと仮定すると、 f の 1 次導関数の近似は 次のようになります。 R 1 ( x ) {\displaystyle R_{1}(x)} f ′ ( x 0 ) ≈ f ( x 0 + h ) − f ( x 0 ) h . {\displaystyle f'(x_{0})\approx {f(x_{0}+h)-f(x_{0}) \over h}.}
これは、ゼロに向かう極限を除いて 、微分の定義に似ています
(このメソッドの名前はこれにちなんで付けられています)。 f ′ ( x 0 ) = lim h → 0 f ( x 0 + h ) − f ( x 0 ) h . {\displaystyle f'(x_{0})=\lim _{h\to 0}{\frac {f(x_{0}+h)-f(x_{0})}{h}}.}
正確さと秩序 ある手法の解における誤差は、近似解と正確な解析解との差として定義されます。有限差分法における誤差の2つの原因は 、小数点以下の数値をコンピュータで丸めることによる精度の低下 である丸め誤差 と、元の微分方程式の正確な解と、完全な算術演算（丸めなし）を仮定した場合の正確な数値との差である 切り捨て誤差 または 離散化誤差です。
有限差分法は、グリッド上で関数を離散化することに依存します。 有限差分法を用いて問題の解を近似するには、まず問題の領域を離散化する必要があります。これは通常、領域を均一なグリッドに分割することによって行われます（図を参照）。つまり、有限差分法は、多くの場合「時間ステップ」方式で、導関数の離散的な数値近似値の集合を生成します。
一般的に興味深い式は、 ある方法の 局所的打ち切り誤差です。通常、 Big-O 記法 を用いて表される局所的打ち切り誤差は、ある方法の 1 回の適用から生じる誤差を指します。つまり、 が 正確な値と数値近似値を参照する 場合、これは量です 。テイラー多項式の剰余項 は、 局所的打ち切り誤差 を解析するために使用できます。 に対するテイラー多項式の剰余の ラグランジュ形式 を用いると、 局所的打ち切り誤差の支配的な項を見つけることができます。例えば、 であることが分かっている最初の導関数の前進差分式を再び用いて、 代数的な操作を加えると、これは次の式に至り、さらに左側の量は有限差分法による近似値であり、右側の量は関心のある正確な量に剰余を加えたものであることに 留意 すると、明らかに剰余が局所的打ち切り誤差です。この例の最終的な式とその順序は次のようになります。 f ′ ( x i ) − f i ′ {\displaystyle f'(x_{i})-f'_{i}} f ′ ( x i ) {\displaystyle f'(x_{i})} f i ′ {\displaystyle f'_{i}} f ( x 0 + h ) {\displaystyle f(x_{0}+h)} R n ( x 0 + h ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ( h ) n + 1 , x 0 < ξ < x 0 + h , {\displaystyle R_{n}(x_{0}+h)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}(h)^{n+1}\,,\quad x_{0}<\xi <x_{0}+h,} f ( x i ) = f ( x 0 + i h ) {\displaystyle f(x_{i})=f(x_{0}+ih)} f ( x 0 + i h ) = f ( x 0 ) + f ′ ( x 0 ) i h + f ″ ( ξ ) 2 ! ( i h ) 2 , {\displaystyle f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}(ih)^{2},} f ( x 0 + i h ) − f ( x 0 ) i h = f ′ ( x 0 ) + f ″ ( ξ ) 2 ! i h , {\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+{\frac {f''(\xi )}{2!}}ih,} f ( x 0 + i h ) − f ( x 0 ) i h = f ′ ( x 0 ) + O ( h ) . {\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+O(h).}
この場合、局所的な打ち切り誤差はステップサイズに比例します。シミュレーションされるFDM解の品質と時間は、離散化方程式の選択とステップサイズ（時間ステップと空間ステップ）に依存します。ステップサイズが小さいほど、データ品質とシミュレーション時間が大幅に増加します。 [2] したがって、実用化には、データ品質とシミュレーション時間の間の適切なバランスが必要です。大きな時間ステップは、実際にはシミュレーション速度を向上させるのに役立ちます。しかし、時間ステップが大きすぎると、不安定性が生じ、データ品質に影響を与える可能性があります。 [3] [4]
フォン ・ノイマン 基準と クーラン・フリードリヒス・レヴィ 基準は、数値モデルの安定性を決定するためによく評価されます。 [3] [4] [5] [6]
例: 常微分方程式 たとえば、常微分方程式を考えてみましょう。 この方程式を解くオイラー法では、差分商を使用して微分 方程式 を 近似します。まずそれを に代入し、 次に少し代数を適用して (両辺に h を掛けて、 両辺に加算する)、次の式を得ます。 最後の方程式は差分方程式であり、この方程式を解くと、微分方程式の近似解が得られます。 u ′ ( x ) = 3 u ( x ) + 2. {\displaystyle u'(x)=3u(x)+2.} u ( x + h ) − u ( x ) h ≈ u ′ ( x ) {\displaystyle {\frac {u(x+h)-u(x)}{h}}\approx u'(x)} u ′ ( x ) {\displaystyle u'(x)} u ( x ) {\displaystyle u(x)} u ( x + h ) ≈ u ( x ) + h ( 3 u ( x ) + 2 ) . {\displaystyle u(x+h)\approx u(x)+h(3u(x)+2).}
例: 熱方程式 同次ディリクレ境界条件 のもとで、1次元の 正規 化熱方程式を考える。
{ U t = U x x U ( 0 , t ) = U ( 1 , t ) = 0 (boundary condition) U ( x , 0 ) = U 0 ( x ) (initial condition) {\displaystyle {\begin{cases}U_{t}=U_{xx}\\U(0,t)=U(1,t)=0&{\text{(boundary condition)}}\\U(x,0)=U_{0}(x)&{\text{(initial condition)}}\end{cases}}}
この方程式を数値的に解く一つの方法は、すべての導関数を有限差分で近似することです。まず、空間領域をメッシュ で 、時間領域をメッシュ で分割します 。空間と時間の両方で均一な分割を仮定すると、連続する2つの空間点間の差は h 、連続する2つの時間点間の差は k となります。点 x 0 , … , x J {\displaystyle x_{0},\dots ,x_{J}} t 0 , … , t N {\displaystyle t_{0},\dots ,t_{N}}
u ( x j , t n ) = u j n {\displaystyle u(x_{j},t_{n})=u_{j}^{n}}
の数値近似値を表す。 u ( x j , t n ) . {\displaystyle u(x_{j},t_{n}).}
明示的な方法 熱方程式の最も一般的な明示的方法の ステンシル 。 時間における 前方差分 と位置における空間微分に対する 2次 中心差分 ( FTCS )を使用すると 、次の再帰方程式が得られます。 t n {\displaystyle t_{n}} x j {\displaystyle x_{j}}
u j n + 1 − u j n k = u j + 1 n − 2 u j n + u j − 1 n h 2 . {\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}.}
これは 1次元 熱方程式を解くための 明示的な方法 です。
他の値からは次のように 取得できます。 u j n + 1 {\displaystyle u_{j}^{n+1}}
u j n + 1 = ( 1 − 2 r ) u j n + r u j − 1 n + r u j + 1 n {\displaystyle u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}
どこ r = k / h 2 . {\displaystyle r=k/h^{2}.}
したがって、この再帰関係と時刻n における値が分かれば 、時刻 n +1 における対応する値を取得できます。 は 境界条件に置き換える必要があり、この例では両方とも 0 です。 u 0 n {\displaystyle u_{0}^{n}} u J n {\displaystyle u_{J}^{n}}
この明示的解法は、 常に 数値的に安定 かつ 収束する ことが知られている。 [7] 数値誤差は時間ステップと空間ステップの2乗に比例する。 r ≤ 1 / 2 {\displaystyle r\leq 1/2} Δ u = O ( k ) + O ( h 2 ) {\displaystyle \Delta u=O(k)+O(h^{2})}
暗黙法 暗黙的メソッドのステンシル。 時間における 後方差分 と位置における空間微分の 2 次中心差分 (後方時間、中心空間法 "BTCS") を使用すると、次の再帰方程式が得られます。 t n + 1 {\displaystyle t_{n+1}} x j {\displaystyle x_{j}}
u j n + 1 − u j n k = u j + 1 n + 1 − 2 u j n + 1 + u j − 1 n + 1 h 2 . {\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}.}
これは 1次元 熱方程式を解く 暗黙的な方法 である。
線形方程式系を解くと 次の式が得られます。 u j n + 1 {\displaystyle u_{j}^{n+1}}
( 1 + 2 r ) u j n + 1 − r u j − 1 n + 1 − r u j + 1 n + 1 = u j n {\displaystyle (1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}}
この手法は常に 数値的に安定 かつ収束しますが、各時間ステップで数値方程式系を解く必要があるため、通常は陽解法よりも数値計算量が多くなります。誤差は時間ステップに対して線形、空間ステップに対しては二次関数です。 Δ u = O ( k ) + O ( h 2 ) . {\displaystyle \Delta u=O(k)+O(h^{2}).}
クランク・ニコルソン法 最後に、時間での中心差分 と位置での空間微分の 2 次中心差分 (「CTCS」) を使用して、再帰方程式を生成します。 t n + 1 / 2 {\displaystyle t_{n+1/2}} x j {\displaystyle x_{j}}
u j n + 1 − u j n k = 1 2 ( u j + 1 n + 1 − 2 u j n + 1 + u j − 1 n + 1 h 2 + u j + 1 n − 2 u j n + u j − 1 n h 2 ) . {\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {1}{2}}\left({\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}+{\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}\right).}
この式は クランク・ニコルソン法 として知られています。
クランク・ニコルソン ステンシル。 線形方程式系を解くと 次の式が得られます。 u j n + 1 {\displaystyle u_{j}^{n+1}}
( 2 + 2 r ) u j n + 1 − r u j − 1 n + 1 − r u j + 1 n + 1 = ( 2 − 2 r ) u j n + r u j − 1 n + r u j + 1 n {\displaystyle (2+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=(2-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}
この手法は常に 数値的に安定 かつ収束しますが、各時間ステップで数値方程式系を解く必要があるため、通常は数値計算量が多くなります。誤差は時間ステップと空間ステップの両方で2次関数となります。 Δ u = O ( k 2 ) + O ( h 2 ) . {\displaystyle \Delta u=O(k^{2})+O(h^{2}).}
比較 まとめると、通常、 クランク・ニコルソン法は 小さな時間ステップでは最も精度の高い法です。大きな時間ステップでは、計算負荷が低い暗黙法の方が適しています。陽的法は最も精度が低く、不安定になる可能性がありますが、実装が最も簡単で、計算負荷も最も低くなります。
ここに例を示します。下の図は、上記の方法によって得られた熱方程式の近似解を示しています。
U t = α U x x , α = 1 π 2 , {\displaystyle U_{t}=\alpha U_{xx},\quad \alpha ={\frac {1}{\pi ^{2}}},}
境界条件付き
U ( 0 , t ) = U ( 1 , t ) = 0. {\displaystyle U(0,t)=U(1,t)=0.}
正確な解決策は
U ( x , t ) = 1 π 2 e − t sin ( π x ) . {\displaystyle U(x,t)={\frac {1}{\pi ^{2}}}e^{-t}\sin(\pi x).}
例: ラプラス演算子 次元における （連続） ラプラス演算子 は で与えられる 。離散ラプラス演算子は 次元 に依存する 。 n {\displaystyle n} Δ u ( x ) = ∑ i = 1 n ∂ i 2 u ( x ) {\displaystyle \Delta u(x)=\sum _{i=1}^{n}\partial _{i}^{2}u(x)} Δ h u {\displaystyle \Delta _{h}u} n {\displaystyle n}
1次元では、ラプラス演算子は次のように近似されます。 この近似は通常、対称三角行列を表す次の ステンシル とで表されます
。等間隔グリッドの場合は、 テプリッツ行列 が得られます。 Δ u ( x ) = u ″ ( x ) ≈ u ( x − h ) − 2 u ( x ) + u ( x + h ) h 2 =: Δ h u ( x ) . {\displaystyle \Delta u(x)=u''(x)\approx {\frac {u(x-h)-2u(x)+u(x+h)}{h^{2}}}=:\Delta _{h}u(x)\,.} Δ h = 1 h 2 [ 1 − 2 1 ] {\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}1&-2&1\end{bmatrix}}}
2次元の場合、より一般的なn次元の場合のすべての特徴を示す。各2次偏微分は、1次元の場合と同様に近似する必要があり 、通常は次の ステンシルで与えられる。 Δ u ( x , y ) = u x x ( x , y ) + u y y ( x , y ) ≈ u ( x − h , y ) − 2 u ( x , y ) + u ( x + h , y ) h 2 + u ( x , y − h ) − 2 u ( x , y ) + u ( x , y + h ) h 2 = u ( x − h , y ) + u ( x + h , y ) − 4 u ( x , y ) + u ( x , y − h ) + u ( x , y + h ) h 2 =: Δ h u ( x , y ) , {\displaystyle {\begin{aligned}\Delta u(x,y)&=u_{xx}(x,y)+u_{yy}(x,y)\\&\approx {\frac {u(x-h,y)-2u(x,y)+u(x+h,y)}{h^{2}}}+{\frac {u(x,y-h)-2u(x,y)+u(x,y+h)}{h^{2}}}\\&={\frac {u(x-h,y)+u(x+h,y)-4u(x,y)+u(x,y-h)+u(x,y+h)}{h^{2}}}\\&=:\Delta _{h}u(x,y)\,,\end{aligned}}} Δ h = 1 h 2 [ 1 1 − 4 1 1 ] . {\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}&1\\1&-4&1\\&1\end{bmatrix}}\,.}
一貫性 上記の近似の整合性は、 のような非常に規則的な関数に対しても示せます 。その記述は u ∈ C 4 ( Ω ) {\displaystyle u\in C^{4}(\Omega )} Δ u − Δ h u = O ( h 2 ) . {\displaystyle \Delta u-\Delta _{h}u={\mathcal {O}}(h^{2})\,.}
これを証明するには、離散ラプラス演算子に 3 次までの テイラー級数 展開を代入する必要があります。
プロパティ
低調波 連続分数調波関数 と同様に、 有限差分近似の 分数調波関数 を定義することができる。 u h {\displaystyle u_{h}} − Δ h u h ≤ 0 . {\displaystyle -\Delta _{h}u_{h}\leq 0\,.}
平均値 ポジティブ型 の一般的な ステンシル は次のように 定義できる。 [ α N α W − α C α E α S ] , α i > 0 , α C = ∑ i ∈ { N , E , S , W } α i . {\displaystyle {\begin{bmatrix}&\alpha _{N}\\\alpha _{W}&-\alpha _{C}&\alpha _{E}\\&\alpha _{S}\end{bmatrix}}\,,\quad \alpha _{i}>0\,,\quad \alpha _{C}=\sum _{i\in \{N,E,S,W\}}\alpha _{i}\,.}
が （離散）サブハーモニックである場合、次の 平均値プロパティ が保持されます 。ここで、近似はグリッドのポイントで評価され、ステンシルは正のタイプであると想定されます。 u h {\displaystyle u_{h}} u h ( x C ) ≤ ∑ i ∈ { N , E , S , W } α i u h ( x i ) ∑ i ∈ { N , E , S , W } α i , {\displaystyle u_{h}(x_{C})\leq {\frac {\sum _{i\in \{N,E,S,W\}}\alpha _{i}u_{h}(x_{i})}{\sum _{i\in \{N,E,S,W\}}\alpha _{i}}}\,,}
同様の 平均値特性は 連続的なケースにも当てはまります。
最大原則 (離散)分数調和関数の場合、 次が成り立ちます
。 ここで、 は連続領域 の離散化 、 は境界 です 。 u h {\displaystyle u_{h}} max Ω h u h ≤ max ∂ Ω h u h , {\displaystyle \max _{\Omega _{h}}u_{h}\leq \max _{\partial \Omega _{h}}u_{h}\,,} Ω h , ∂ Ω h {\displaystyle \Omega _{h},\partial \Omega _{h}} Ω {\displaystyle \Omega } ∂ Ω {\displaystyle \partial \Omega }
同様の 最大原理 は連続的な場合にも当てはまります。
SBP-SAT法 SBP-SAT法（ 部分和 - 同時近似項 ）は、高次差分を用いて 適切に設定された線形 偏微分方程式を離散化し境界条件を課すための安定した正確な手法である。 [8] [9]
この手法は有限差分法に基づいており、微分演算子は 部分和 特性を示します。通常、これらの演算子は、内部に中央差分ステンシルを持つ微分行列と、離散設定における部分積分を模倣するように設計された慎重に選択された片側境界ステンシルで構成されます。SAT法を用いると、偏微分方程式の境界条件は弱く課され、境界値は所望の条件に厳密に適合するのではなく、望ましい条件に向かって「引き寄せられる」ようになります。SAT法に固有の調整パラメータが適切に選択されると、結果として得られる常微分方程式系は連続偏微分方程式と同様のエネルギー挙動を示します。つまり、系には非物理的なエネルギー増加はありません。これにより、4次 ルンゲ・クッタ法 など、虚軸の一部を含む安定領域を持つ積分法を用いる場合、安定性が保証されます。このため、SAT 手法は、高次微分演算子が使用されると通常は安定しない注入法などとは対照的に、高次有限差分法に境界条件を課す魅力的な方法となります。
参照
参考文献 ^ ab Christian Grossmann; Hans-G. Roos; Martin Stynes (2007). 偏微分方程式の数値的処理 . Springer Science & Business Media. p. 23. ISBN 978-3-540-71584-9 。 ^ Arieh Iserles (2008). 微分方程式の数値解析入門 . ケンブリッジ大学出版局. p. 23. ISBN 9780521734905 。 ^ ab Hoffman JD; Frankel S (2001). エンジニアと科学者のための数値解析法 . CRC Press, Boca Raton. ^ ab Jaluria Y; Atluri S (1994). 「計算熱伝達」. 計算力学 . 14 (5): 385– 386. Bibcode :1994CompM..14..385J. doi :10.1007/BF00377593. S2CID 119502676. ^ Majumdar P (2005). 熱および質量移動の計算手法 （第1版）. Taylor and Francis, New York. ^ Smith GD (1985). 偏微分方程式の数値解法：有限差分法 （第3版）. オックスフォード大学出版局. ^ Crank, J. 『拡散の数学 』第2版、オックスフォード、1975年、143ページ。 ^ Bo Strand (1994). 「d/dxの有限差分近似における部分和」. Journal of Computational Physics . 110 (1): 47– 67. Bibcode :1994JCoPh.110...47S. doi :10.1006/jcph.1994.1005. ^ Mark H. Carpenter、David I. Gottlieb、Saul S. Abarbanel (1994). 「双曲型システムを解く有限差分スキームの時間安定境界条件：方法論と高次コンパクトスキームへの応用」 Journal of Computational Physics . 111 (2): 220– 236. Bibcode :1994JCoPh.111..220C. doi :10.1006/jcph.1994.1057. hdl : 2060/19930013937 .
さらに読む KW Morton、DF Mayers著 『偏微分方程式の数値解法入門』 ケンブリッジ大学出版局、2005年。 Autar KawとE. Eric Kalu著 「Numerical Methods with Applications 」（2008年）[1]。第8章第7節には、工学的な観点から簡潔なFDM（常微分方程式用）の解説が掲載されている。 ジョン・ストリクヴェルダ (2004). 有限差分スキームと偏微分方程式 （第2版）. SIAM. ISBN 978-0-89871-639-9 。 スミス、GD（1985）、 偏微分方程式の数値解法：有限差分法、第3版 、オックスフォード大学出版局 ピーター・オルバー （2013年）偏微分方程式入門、シュプリンガー社。第5章 有限差分 。ISBN 978-3-319-02099-0 。 。 Randall J. LeVeque 、 「常微分方程式と偏微分方程式の有限差分法」 、SIAM、2007 年。 セルゲイ・レメシェフスキー、ピョートル・マトゥス、ドミトリー・ポリアコフ（編）：「正確な有限差分スキーム」、デ・グリュイテル（2016）。 DOI: https://doi.org/10.1515/9783110491326 。 ミハイル・シャシュコフ： 一般グリッドにおける保守的有限差分法 、CRC プレス、ISBN 0-8493-7375-1 (1996)。
