差分法

数値解析における有限差分法︵ゆうげんさぶんほう、英: finite-difference methods; FDM︶あるいは単に差分法は、微分方程式を解くために微分を有限差分近似︵差分商︶で置き換えて得られる差分方程式で近似するという離散化手法を用いる数値解法である。18世紀にオイラーが考案したと言われる。

今日ではFDMは偏微分方程式の数値解法として支配的な手法である。

精度と誤差

「有限差分係数」も参照

解の誤差とは、真の解析解と近似解との間の差として定義される。有限差分法における誤差の原因は丸め誤差および打ち切り誤差または離散化誤差である。

有限差分法は函数の定義域を格子に離散化することに基づく

問題に対する解の近似に有限差分法を用いるためには、まず初めに問題の領域を離散化しなければならない。これは普通は、その領域を一様な格子に分ければよい。これは有限差分法がしばしば﹁時間刻み﹂な仕方で微分に対する離散的な数値近似の集合を提供することを意味することに注意。

$f(x_{i})=f(x_{0}+ih)$$f(x_{i})=f(x_{0}+ih)$.

一般に注目すべきは局所打ち切り誤差で、典型的にはこれを O-記法で表す。局所打ち切り誤差は、各点における誤差について言うもので、真値 f'(x_i) と近似値 f'_i との差

$f'(x_{i})-f'_{i}$$f'(x_{i})-f'_{i}$
である。この誤差の評価には、テイラー展開の剰余項を見るのが簡便である。式 f(x₀ + h) に対するテイラー展開のラグランジュ型剰余項

${\displaystyle R_{n}(x_{0}+h)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}($h)^{n+1}\quad (x_{0}<\xi <x_{0}+h)} $R_{n}(x_{0}+h)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}(h)^{n+1}\quad (x_{0}<\xi <x_{0}+h)$
から、局所打ち切り誤差の支配項が求められる。例えば、一階差分近似 (n = 1) を考えれば

${\displaystyle f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}($ih)^{2}} $f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}(ih)^{2}$
である。この右辺は有限差分法で得られる近似値である。一方、0階差分近似︵n=0)を考えれば

$f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih$$f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih$
よって、0階差分近似での支配的な誤差は

${\displaystyle {\frac {f''(\xi )}{2!}}($ih)^{2}} ${\frac {f''(\xi )}{2!}}(ih)^{2}$
であり、この剰余項(n=1)が局所打ち切り誤差の支配項である。この場合、局所打ち切り誤差はほぼ刻み幅(h)の2乗に比例するということになる。有限差分法の近似解の精度と計算量は方程式の離散化の仕方や刻み幅の取り方に依存する。これらは刻み幅を小さくするにつれ著しく増加するから、実用上は必要な精度と計算時間を天秤にかけて十分合理的な条件で近似を行う。時間の刻み幅が大きければ多くの場合に計算速度は早くなるが、大きくしすぎると不安定性を生じ、データの精度に問題がでる。

数値モデルの安定性を決定するために、フォン・ノイマンの安定性解析を用いるのが普通である。

簡単な例

最も簡単な例として、次の1階常微分方程式を考える‥

${\displaystyle u'($x)=3u(x)+2\,} $u'(x)=3u(x)+2\,$
これを解くには、差分商

${\displaystyle {\frac {u(x+h)-u($x)}{h}}\approx u'(x)\quad (h\to 0)} ${\frac {u(x+h)-u(x)}{h}}\approx u'(x)\quad (h\to 0)$
を用いて

${\displaystyle u(x+h)=u($x)+h(3u(x)+2)\,} $u(x+h)=u(x)+h(3u(x)+2)\,$
と近似する。この方法をオイラー法という。この最後の方程式のように、微分方程式の微分を差分商に置き換えたものを、差分方程式︵さぶんほうていしき、difference equation︶と呼ぶ。

例 熱伝導方程式

偏微分方程式の例として、一様ディリクレ境界条件に従う1次元規格化熱伝導方程式を考える‥

$U_{t}=U_{xx}\,$$U_{t}=U_{xx}\,$
左辺は時刻$t$$t$による微分、右辺は座標$x$$x$による2階微分である。また、境界条件および初期条件は以下とする‥

$U(0,t)=U(1,t)=0\,$$U(0,t)=U(1,t)=0\,$︵境界条件︶

${\displaystyle U(x,0)=U_{0}($x)\,} $U(x,0)=U_{0}(x)\,$︵初期条件︶

これを数値的に解く1つの方法は、すべての微分を差分で近似することである。空間の領域をメッシュ$x_{0},\dots ,x_{J}$$x_{0},\dots ,x_{J}$で、時間の領域をメッシュ$t_{0},\dots ,t_{N}$$t_{0},\dots ,t_{N}$で分割しよう。どちらの分割も等間隔とし、空間点の間隔を$h$$h$、時刻の間隔を$k$$k$とする。$U(x_{j},t_{n})$$U(x_{j},t_{n})$の数値的近似を$u_{j}^{n}$$u_{j}^{n}$で表す。

陽解法

時刻$t_{n}$$t_{n}$には前進差分を用い、空間点$x_{j}$$x_{j}$で2次微分に対して2次中央差分を用いれば、次の漸化式‥

${\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}\,$${\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}\,$
が得られる。これを陽解法という。

$u_{j}^{n+1}$$u_{j}^{n+1}$の値は次のように得られる‥

$u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}$$u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}$
ただしここで$r=k/h^{2}$$r=k/h^{2}$ ︵拡散数と呼ばれる︶である。

ゆえに、時刻$t_{n}$$t_{n}$での値がわかれば、対応する時刻$t_{n+1}$$t_{n+1}$での値も漸化式を用いて求められる。$u_{0}^{n}$$u_{0}^{n}$と$u_{J}^{n}$$u_{J}^{n}$には境界条件︵この例ではどちらも0︶を適用する。

この陽解法は、$r\leq 1/2$$r\leq 1/2$であれば数値的に安定で収束することが知られている。

誤差は時刻間隔$k$$k$の1乗と空間点間隔$h$$h$の2乗のオーダーである‥

${\displaystyle \Delta u=O($k)+O(h^{2})\,} $\Delta u=O(k)+O(h^{2})\,$

陰解法

時刻$t_{n+1}$$t_{n+1}$に後退差分を用い、空間点$x_{j}$$x_{j}$ で2階中央差分を用いれば、漸化式‥

${\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}}}\,$${\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+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}$$(1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}$
を解けば、$u_{j}^{n+1}$$u_{j}^{n+1}$が得られる。この方法は常に数値的に安定で収束するが、時刻ごとに方程式系を解く必要があるため、陽解法よりも繁雑である。誤差は時間ステップ数と空間ステップ数の4乗とに比例する。

クランク・ニコルソン法

さいごに、時刻$t_{n+1/2}$$t_{n+1/2}$で中央差分を、空間点$x_{j}$$x_{j}$での空間微分に2階中央差分を用いれば、漸化式‥

$2\left({\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}\right)={\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}}}\,$$2\left({\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}\right)={\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}}}\,$
が得られる。これをクランク・ニコルソン法︵Crank-Nicolson method︶という。

線形方程式系‥

$(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+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}$
を解けば、$u_{j}^{n+1}$$u_{j}^{n+1}$が得られる。

この方法は常に数値的に安定で収束するが、各時刻で方程式系を解く必要があるので繁雑なことが多い。誤差は時間ステップ数の4乗と空間ステップ数の2乗とに比例する‥

$\Delta u=O(k^{2})+O(h^{4})\,$$\Delta u=O(k^{2})+O(h^{4})\,$
しかし、境界付近では誤差はO(h⁴ ) でなくO(h² ) となることが多い。

クランク・ニコルソン法は時間ステップ数が少なければたいてい最も正確な方法である。陽解法はそれより正確でなく不安定でもあるが、最も実行しやすく、繁雑さも最も少ない。陰解法は時間ステップ数が多い場合に最も優れている。

関連項目

• 有限体積法
• 有限要素法
• ドゥッチ数列

外部リンク

• Finite difference method （英語） - スカラーペディア百科事典「差分法」の項目。