差分法

数値解析において、有限差分法( FDM ) は、導関数を有限差分で近似することで微分方程式を解く数値手法の一種です。空間領域と時間領域 (該当する場合) の両方が離散化、つまり有限数の区間に分割され、区間の終点における解の値は、有限差分と近傍点からの値を含む代数方程式を解くことで近似されます。

有限差分法は、非線形の場合もある常微分方程式(ODE) または偏微分方程式(PDE) を、行列代数技法で解くことができる線形方程式のシステムに変換します。現代のコンピュータはこれらの線形代数計算を効率的に実行でき、実装が比較的容易であることから、現代​​の数値解析では FDM が広く使用されています。^{[1]今日、FDM は}有限要素法 とともに、PDE の数値解を求める最も一般的なアプローチの 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の一次導関数の近似値を導出するプロセスです。 をh で割ると次のようになります。 を解くと、次のようになります。 $${\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x).}$$$${\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}.}$$

が十分に小さいと仮定すると、 fの一次導関数の近似は次のようになります。 ${\displaystyle R_{1}(x)}$$${\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 つの原因は、小数の量をコンピュータで丸めることによって精度が失われる丸め誤差と、元の微分方程式の正確な解と、完全な算術 (丸めなし) を前提とした正確な量との差である切り捨て 誤差または離散化誤差です。

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

一般的に関心の高い表現は、方法のローカルな切り捨て誤差です。通常、 Big-O 記法で表現されるローカルな切り捨て誤差は、方法の 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]}

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

同次ディリクレ境界条件を持つ1次元の正規化熱方程式を考える

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

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

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

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

時間における前方差分 と位置における空間微分に対する2次中心差分( FTCS )を使用すると、再帰方程式が得られます。 ${\displaystyle t_{n}}$${\displaystyle x_{j}}$

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

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

他の値から次のように 取得できます。${\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\Delta t/h^{2}.}$

したがって、この再帰関係と、時刻nにおける値がわかっている場合、時刻n +1における対応する値を取得できます。および は境界条件に置き換える必要があり、この例では両方とも 0 です。 ${\displaystyle u_{0}^{n}}$${\displaystyle u_{J}^{n}}$

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

時間における後方差分 と位置における空間微分に対する 2 次中心差分を使用すると(後方時間、中心空間法 "BTCS")、再帰方程式が得られます。 ${\displaystyle t_{n+1}}$${\displaystyle x_{j}}$

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

これは1次元熱方程式を解くための暗黙的な方法です。

線形方程式系を解くと 次の式が得られます。${\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}}$$

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

最後に、時間での中心差分と位置での空間微分に対する 2 次中心差分(「CTCS」) を使用すると、再帰方程式が得られます。 ${\displaystyle t_{n+1/2}}$${\displaystyle x_{j}}$

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k\Delta t}}={\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}$

1D では、ラプラス演算子は次のように近似されます 。 この近似は通常、次のステンシル と、対称三角行列を表すことによって表現されます。等距離グリッドの場合は、テプリッツ行列が得られます。 $${\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}}\,.}$$

上記の近似の一貫性は、 のような非常に規則的な関数に対しても示せます。 ${\displaystyle u\in C^{4}(\Omega )}$$${\displaystyle \Delta u-\Delta _{h}u={\mathcal {O}}(h^{2})\,.}$$

これを証明するには、離散ラプラス演算子に 3 次までの テイラー級数展開を代入する必要があります。

連続分数調波関数と同様に、有限差分近似の分数調波関数を定義することができる。${\displaystyle u_{h}}$ $${\displaystyle -\Delta _{h}u_{h}\leq 0\,.}$$

ポジティブタイプの一般的なステンシルは次のように 定義できる。$${\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}\,.}$$

が（離散）サブハーモニックである場合、次の平均値特性が成立します 。ここで、 近似はグリッドのポイント上で評価され、ステンシルは正のタイプであると想定されます。 ${\displaystyle u_{h}}$$${\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}}}\,,}$$

同様の平均値特性は連続的なケースにも当てはまります。

（離散）分数調和関数については、次が成り立ちます 。 ここで、 は連続領域 の離散化、 は境界です。 ${\displaystyle u_{h}}$$${\displaystyle \max _{\Omega _{h}}u_{h}\leq \max _{\partial \Omega _{h}}u_{h}\,,}$$${\displaystyle \Omega _{h},\partial \Omega _{h}}$${\displaystyle \Omega }$${\displaystyle \partial \Omega }$

同様の最大原理は連続的な場合にも当てはまります。

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

この方法は、微分演算子が部分和特性を示す有限差分に基づいています。通常、これらの演算子は、離散設定で部分積分を模倣するように設計された、慎重に選択された片側境界ステンシルを備えた内部の中央差分ステンシルを持つ微分行列で構成されます。SAT 手法を使用すると、PDE の境界条件が弱く課され、境界値は、正確に満たされるのではなく、目的の条件に向かって「引き寄せられ」ます。チューニング パラメータ (SAT 手法に固有のもの) が適切に選択されると、結果として得られる ODE システムは、連続 PDE と同様のエネルギー動作を示します。つまり、システムには非物理的なエネルギー増加はありません。これにより、4 次ルンゲ クッタ法などの、虚軸の一部を含む安定領域を持つ積分スキームが使用される場合、安定性が保証されます。このため、SAT 手法は、高次微分演算子が使用されると通常は安定しない注入法などとは対照的に、高次有限差分法に境界条件を課す魅力的な方法となります。

• KW Morton および DF Mayers、「偏微分方程式の数値解法入門」、ケンブリッジ大学出版局、2005 年。
• Autar KawとE. Eric Kalu、「Numerical Methods with Applications」（2008）[1]。第08.07章には、FDM（ODE用）の簡潔なエンジニアリング指向の紹介が含まれています。
• ジョン・ストリクヴェルダ (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 。