差分法

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

差分法(FDM)は有限要素法(FEM)や境界要素法(BEM)などと並んで偏微分方程式の代表的な数値解析手法の1つである。

## 精度と誤差

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

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

一般に注目すべきはで、典型的にはこれをO-記法で表す。局所打ち切り誤差は、各点における誤差について言うもので、真値  と近似値  との差
: f'(x_i) - f'_i
である。この誤差の評価には、テイラー展開の剰余項を見るのが簡便である。式  に対するテイラー展開のラグランジュ型剰余項
: R_n(x_0 + h) = \frac{f^{(n+1)}(\xi)}{(n+1)!} (h)^{n+1}\quad (x_0 
から、局所打ち切り誤差の支配項が求められる。例えば、一階差分近似 () を考えれば
:  f(x_0 + i h) = f(x_0) + f'(x_0)i h + \frac{f*(\xi)}{2!} (i h)^{2}</math>
である。この右辺は有限差分法で得られる近似値である。一方、0階差分近似（n=0)を考えれば

<math> f(x_0 + i h) = f(x_0) + f'(x_0)i h</math>

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

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

## 簡単な例
最も簡単な例として、次の1階常微分方程式を考える：
: u'(x) = 3u(x) + 2 \, 
これを解くには、差分商
:\frac{u(x+h) - u(x)}{h} \approx u'(x)\quad (h \to 0)
を用いて
: u(x+h) = u(x) + h(3u(x)+2) \, 
と近似する。この方法を**[[オイラー法]]**という。この最後の方程式のように、微分方程式の微分を差分商に置き換えたものを、**差分方程式**（さぶんほうていしき、）と呼ぶ。

## 例 熱伝導方程式
偏微分方程式の例として、一様ディリクレ境界条件に従う1次元規格化熱伝導方程式を考える：
: U_t=U_{xx} \, 
左辺は時刻<math> t </math>による微分、右辺は座標 x による2階微分である。また、境界条件および初期条件は以下とする：
: U(0,t)=U(1,t)=0 \, （境界条件）
: U(x,0) =U_0(x) \, （初期条件）
これを数値的に解く1つの方法は、すべての微分を差分で近似することである。空間の領域をメッシュ x_0, \dots , x_J で、時間の領域をメッシュ t_0,\dots, t_N で分割しよう。どちらの分割も等間隔とし、空間点の間隔を h 、時刻の間隔を k とする。 U(x_j, t_n) の数値的近似を u_j^n で表す。

### 陽解法
時刻 t_n には前進差分を用い、空間点 x_j で2次微分に対して2次中央差分を用いれば、次の漸化式：
:\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} \,
が得られる。これを陽解法という。

 u_j^{n+1} の値は次のように得られる：
: u_{j}^{n+1} = (1-2r)u_{j}^{n} + ru_{j-1}^{n} + ru_{j+1}^{n}
ただしここで r=k \Delta t/h^2  （拡散数と呼ばれる）である。

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

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

誤差は時刻間隔 k の1乗と空間点間隔 h の2乗のオーダーである：
:\Delta u = O(k)+O(h^2) 

### 陰解法
時刻t_{n+1} に後退差分を用い、空間点 x_j  で2階中央差分を用いれば、漸化式：
:\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+2r)u_j^{n+1} - ru_{j-1}^{n+1} - ru_{j+1}^{n+1}= u_j^{n}
を解けば、 u_j^{n+1} が得られる。この方法は常に数値的に安定で収束するが、時刻ごとに方程式系を解く必要があるため、陽解法よりも繁雑である。誤差は時間ステップ数と空間ステップ数の4乗とに比例する。

### クランク・ニコルソン法
さいごに、時刻 t_{n+1/2} で中央差分を、空間点 x_j での空間微分に2階中央差分を用いれば、漸化式：
: 2\left(\frac{u_j^{n+1} - u_j^{n}}{k \Delta t }\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 
を解けば、 u_j^{n+1} が得られる。

この方法は常に数値的に安定で収束するが、各時刻で方程式系を解く必要があるので繁雑なことが多い。誤差は時間ステップ数の4乗と空間ステップ数の2乗とに比例する：
: \Delta u = O(k^2)+O(h^4)  \, 
しかし、境界付近では誤差はO(*h*4 ) でなくO(*h*2 ) となることが多い。

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

## 参考文献

## 関連文献
* 
* 水本久夫：「多様体上の差分法」、教育出版（シリーズ新しい応用の数学 2）(1973年11月10日)。※ リーマン面上の差分法など。

## 関連項目
* 有限体積法
* 有限要素法
* ドゥッチ数列

## 外部リンク
*  - スカラーペディア百科事典「差分法」の項目。