有限差分法

数値解析において、有限差分法（FDM ）は、微分方程式を有限差分で近似することにより解く数値手法の一種です。空間領域と時間領域（該当する場合）の両方が離散化され、有限個の区間に分割されます。そして、区間の端点における解の値は、有限差分と近傍の点の値を含む代数方程式を解くことによって近似されます。

有限差分法は、非線形である可能性のある常微分方程式（ODE）または偏微分方程式（PDE）を、行列代数技術で解くことができる線形方程式系に変換します。現代のコンピュータはこれらの線形代数計算を効率的に実行でき、その実装の容易さと相まって、現代の数値解析において有限差分法が広く使用されるようになりました。^{[ 1 ]}今日、有限差分法は、有限要素法 と並んで、偏微分方程式の数値解法の最も一般的なアプローチの1つです。^{[ 1 ]}

n回微分可能な関数に対して、テイラーの定理により、テイラー級数展開は次のように与えられる。 $${\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の一次導関数の近似値を導出する手順は次のとおりです。まず、テイラー多項式と剰余項を切り捨てます。 $${\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x).}$$h で割ると次のようになります。 $${\displaystyle {f(x_{0}+h) \over h}={f(x_{0}) \over h}+f'(x_{0})+{R_{1}(x) \over h}}$$ 解決する${\displaystyle f'(x_{0})}$: $${\displaystyle f'(x_{0})={f(x_{0}+h)-f(x_{0}) \over h}-{R_{1}(x) \over h}.}$$

と仮定すると${\displaystyle R_{1}(x)}$が十分に小さい場合、 fの 1 階微分の近似値は次のようになります。 $${\displaystyle f'(x_{0})\approx {f(x_{0}+h)-f(x_{0}) \over h}.}$$

これは導関数の定義と似ており、導関数の定義は次のとおりです。 $${\displaystyle f'(x_{0})=\lim _{h\to 0}{\frac {f(x_{0}+h)-f(x_{0})}{h}}.}$$ ただし、ゼロへの制限を除く（このメソッドはこれにちなんで名付けられている）。

手法の解における誤差は、近似解と厳密な解析解との差として定義されます。有限差分法における誤差の2つの原因は、1つは丸め誤差（コンピュータによる小数値の丸め処理による精度の低下）、もう1つは切り捨て誤差または離散化誤差（元の微分方程式の厳密解と、完全な算術演算（丸め誤差なし）を仮定した場合の厳密な値との差）です。

有限差分法を用いて問題の解を近似するには、まず問題領域を離散化する必要があります。これは通常、領域を均一なグリッドに分割することによって行われます（図を参照）。つまり、有限差分法は、多くの場合「時間ステップ」方式で、導関数の離散的な数値近似値のセットを生成します。

一般的に関心を集める表現として、メソッドの局所打ち切り誤差があります。通常、ビッグオー記法を用いて表現される局所打ち切り誤差は、メソッドを1回適用した際の誤差を指します。つまり、それは次の量です。${\displaystyle f'(x_{i})-f'_{i}}$もし${\displaystyle f'(x_{i})}$正確な値と${\displaystyle f'_{i}}$数値近似に対して。テイラー多項式の剰余項は、局所的な打ち切り誤差を解析するために使用できます。テイラー多項式の剰余のラグランジュ形式を使用して、${\displaystyle f(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,}$$ 局所打ち切り誤差の支配的な項を発見することができる。例えば、再び一次導関数の前方差分公式を使用し、${\displaystyle f(x_{i})=f(x_{0}+ih)}$、 $${\displaystyle f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}(ih)^{2},}$$ そして、いくつかの代数的な操作によって、これは次のようになる。 $${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+{\frac {f''(\xi )}{2!}}ih,}$$ さらに、左辺の量は有限差分法による近似値であり、右辺の量は目的とする正確な量に剰余を加えたものであることを考慮すると、明らかにその剰余は局所的な打ち切り誤差である。この例とその次数を最終的に表すと次のようになる。 $${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+O(h).}$$

この場合、局所的な打ち切り誤差はステップサイズに比例します。シミュレーションされた FDM 解の品質と所要時間は、離散化方程式の選択とステップサイズ (時間ステップと空間ステップ) に依存します。ステップサイズが小さいほど、データ品質とシミュレーション所要時間は大幅に増加します。^{[ 2 ]}したがって、実用上、データ品質とシミュレーション所要時間の適切なバランスが必要です。大きな時間ステップは、実際にはシミュレーション速度を上げるのに役立ちます。ただし、時間ステップが大きすぎると不安定性が発生し、データ品質に影響を与える可能性があります。^{[ 3 ] [ 4 ]}

数値モデルの安定性を判断するために、フォン・ノイマン基準とクーラン・フリードリヒス・レヴィ基準がよく評価される。 ^{[ 3 ] [ 4 ] [ 5 ] [ 6 ]}

例えば、常微分方程式を考えてみましょう。 $${\displaystyle u'(x)=3u(x)+2.}$$ この方程式を解くためのオイラー法は、有限差分商を使用します。 $${\displaystyle {\frac {u(x+h)-u(x)}{h}}\approx u'(x)}$$ 微分方程式を近似するには、まずそれを代入します。${\displaystyle u'(x)}$次に少し代数計算を行う（両辺にhを掛け、さらに${\displaystyle u(x)}$両側に） $${\displaystyle u(x+h)\approx u(x)+h(3u(x)+2).}$$ 最後の式は差分方程式であり、この式を解くことで微分方程式の近似解が得られる。

一次元における正規化された熱方程式を、同次ディリクレ境界条件の下で考察する。

$${\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}}}$$

この方程式を数値的に解く方法の一つは、すべての導関数を有限差分で近似することです。まず、メッシュを使用して空間領域を分割します。${\displaystyle x_{0},\dots ,x_{J}}$そして時間とともにメッシュを使用する${\displaystyle t_{0},\dots ,t_{N}}$空間的にも時間的にも均一な分割を仮定すると、連続する2つの空間点間の差はh 、連続する2つの時間点間の差はkとなる。

$${\displaystyle u(x_{j},t_{n})=u_{j}^{n}}$$

は、数値近似を表します。${\displaystyle u(x_{j},t_{n}).}$

明示的方法

時点のフォワード差分 を使用する${\displaystyle t_{n}}$そして、位置における空間微分に対する2次中心差分${\displaystyle x_{j}}$（FTCS）は漸化式を与える。

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}.}$$

これは、一次元熱方程式を解くための明示的な方法です。

入手できる${\displaystyle u_{j}^{n+1}}$他の値から次のようにして取得します。

$${\displaystyle u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}$$

どこ${\displaystyle r=k/h^{2}.}$

したがって、この漸化式を用い、時刻nにおける値が分かれば、時刻n + 1における対応する値を求めることができる。${\displaystyle u_{0}^{n}}$そして${\displaystyle u_{J}^{n}}$は境界条件に置き換える必要があります。この例では、両方とも0です。

この明示的な方法は、次の場合に数値的に安定かつ収束することが知られています。${\displaystyle r\leq 1/2}$^{[ 7 ]}数値誤差は時間ステップと空間ステップの二乗に比例する 。$${\displaystyle \Delta u=O(k)+O(h^{2})}$$

暗黙法

時刻における後方差分 を使用する${\displaystyle t_{n+1}}$そして、位置における空間微分に対する2次中心差分${\displaystyle x_{j}}$（逆時間中心空間法「BTCS」）は、以下の漸化式を与える。

$${\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}}}.}$$

これは、一次元熱方程式を解くための陰解法です。

入手できる${\displaystyle u_{j}^{n+1}}$連立一次方程式を解くことから：

$${\displaystyle (1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}}$$

この手法は常に数値的に安定かつ収束するが、各時間ステップで数値方程式系を解く必要があるため、通常は陽解法よりも計算負荷が高い。誤差は時間ステップに対して線形、空間ステップに対して二次関数的である。 $${\displaystyle \Delta u=O(k)+O(h^{2}).}$$

クランク・ニコルソン法

最後に、時間における中心差分を使用して${\displaystyle t_{n+1/2}}$そして、位置における空間微分に対する2次中心差分${\displaystyle x_{j}}$（「CTCS」）は、次の漸化式を与える。

$${\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).}$$

この式はクランク・ニコルソン法として知られています。

入手できる${\displaystyle u_{j}^{n+1}}$連立一次方程式を解くことから：

$${\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乗に比例する。 $${\displaystyle \Delta u=O(k^{2})+O(h^{2}).}$$

比較

要約すると、通常、クランク・ニコルソン法は小さな時間ステップにおいて最も精度の高い手法です。大きな時間ステップでは、計算負荷が少ないため、陰解法の方が適しています。陽解法は精度が最も低く不安定になる場合もありますが、実装が最も容易で、数値計算負荷も最も低くなっています。

以下に例を示します。下の図は、上記の方法で熱方程式を近似して得られた解を示しています。

$${\displaystyle U_{t}=\alpha U_{xx},\quad \alpha ={\frac {1}{\pi ^{2}}},}$$

境界条件付き

$${\displaystyle U(0,t)=U(1,t)=0.}$$

正確な解は

$${\displaystyle U(x,t)={\frac {1}{\pi ^{2}}}e^{-t}\sin(\pi x).}$$

有限差分法の比較

明示的な方法（不安定）

暗黙的方法（安定）

クランク・ニコルソン法（安定）

（連続）ラプラス演算子${\displaystyle n}$-次元は次のように与えられます${\displaystyle \Delta u(x)=\sum _{i=1}^{n}\partial _{i}^{2}u(x)}$離散ラプラス演算子${\displaystyle \Delta _{h}u}$次元によって異なる${\displaystyle n}$。

1次元では、ラプラス演算子は次のように近似される。 $${\displaystyle \Delta u(x)=u''(x)\approx {\frac {u(x-h)-2u(x)+u(x+h)}{h^{2}}}=:\Delta _{h}u(x)\,.}$$ この近似は通常、次のステンシルで表現されます。$${\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}1&-2&1\end{bmatrix}}}$$ これは対称な三重対角行列を表します。等間隔のグリッドの場合、トープレッツ行列が得られます。

2次元の場合、より一般的なn次元の場合のすべての特徴が明らかになります。各2階偏微分は、1次元の場合と同様に近似する必要があります。 $${\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}}}$$ これは通常、次のステンシルで示されます$${\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}&1\\1&-4&1\\&1\end{bmatrix}}\,.}$$

SBP-SAT（部分和 - 同時近似項）法は、高次の有限差分を用いて適切に定義された線形偏微分方程式を離散化し、境界条件を課すための安定かつ正確な手法である。 ^{[ 8 ] [ 9 ]}

この手法は、微分演算子が部分和特性を示す有限差分法に基づいています。通常、これらの演算子は、内部に中心差分ステンシルを持ち、離散設定での部分積分を模倣するように設計された、慎重に選択された片側境界ステンシルを持つ微分行列で構成されます。SAT手法を使用すると、偏微分方程式の境界条件は弱く課され、境界値は厳密に満たされるのではなく、望ましい条件に向かって「引き寄せられる」ことになります。SAT手法に固有の調整パラメータが適切に選択されていれば、結果として得られる常微分方程式系は、連続偏微分方程式と同様のエネルギー挙動を示し、つまり、非物理的なエネルギー増加は発生しません。これにより、4次ルンゲ・クッタ法のように、虚軸の一部を含む安定領域を持つ積分スキームが使用される場合、安定性が保証されます。このため、SAT 手法は、例えば高次の微分演算子を使用すると通常は安定しない注入法とは対照的に、高次の有限差分法に境界条件を課す魅力的な方法となります。SAT 項の加法的な定式化により、境界条件を個別の補正項として扱うことで、GPU やその他の最新の高性能アーキテクチャ上でこの手法を効率的に実装できます。これは、境界データを直接注入して境界条件を課す有限差分法とは対照的です。^{[ 10 ]}

• KW Morton および DF Mayers、『偏微分方程式の数値解法入門』、ケンブリッジ大学出版局、2005 年。
• Autar Kaw および E. Eric Kalu、『数値解析とその応用』（2008年）第8.07章には、FDM（常微分方程式用）に関する工学的な入門が簡潔に記述されています。
• ジョン・ストライクワーダ (2004).有限差分法と偏微分方程式(第 2 版). SIAM. ISBN 978-0-89871-639-9。
• Smith, 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 Press、ISBN 0-8493-7375-1（1996年）。