差分法

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索

差分法(さぶんほう、difference method)とは、微分方程式を解く数値解析における離散化手法のひとつである。関数が2つの変数値に対してとる値の間の有限な差を差分(さぶん、difference)といい、この差分を変数値の差で割って得られる商を差分商(さぶんしょう、difference quotient)という。微分を差分商で近似することにより微分方程式を解くものである。18世紀にオイラーが考案したと言われる[1]

偏微分方程式の数値解法を特に有限差分法(ゆうげんさぶんほう、finite difference method, FDM)という。

概要[編集]

最も簡単な例として、次の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) \,

と近似する。この方法をオイラー法という。この最後の方程式のように、微分方程式の微分を差分商に置き換えたものを、差分方程式(さぶんほうていしき、difference equation)と呼ぶ。

変数軸に等間隔h の目盛り(メッシュ)をとり、それらに対する関数u の値で数列を作ることができる。例えば、 \dots, u_{n-1} = u(x-h), u_{n} = u(x), u_{n+1} = u(x+h),\dots とすれば、隣接項の差が差分となる。差分方程式は漸化式として表され、これを解けば微分方程式の近似解が得られる。

なお、離散的な問題にも差分方程式が必要とされるが、この場合には逆に差分の近似として微分を用いることもできる(作用素の項を参照)。

微分方程式の近似解と厳密解との誤差は、差分と微分との誤差(打ち切り誤差)で決まる。

差分の種類[編集]

次のように、本来の微分をとる点x あるいはn に対して、近似としてその次の点n + 1 とn との間で差分をとる方法を、前進(前方)差分forward difference)という:

 \Delta u(x) = u(x+h) - u(x) = u_{n+1} - u_{n}

同様に、n とその前の点n - 1 との間で差分をとる方法を、後退(後方)差分backward difference)という:

 \nabla u(x) = u(x) - u(x-h) = u_{n} - u_{n-1}

また、n に対して、n + 1 とn - 1 との間で差分をとる(前進差分と後退差分を平均する)方法を、中央(中心)差分central difference)という:

 \delta u(x) = u\left(x + \frac{1}{2}h\right) - u\left(x - \frac{1}{2}h\right) = \frac{\Delta u(x) + \nabla u(x)}{2} = \frac{u(x+h)-u(x-h)}{2} = \frac{u_{n+1} - u_{n-1}}{2}

2階微分の近似としては、2階差分2nd difference)を用いる。これは差分の差分に当たる:

 \delta^{2} u(x) = (u_{n+1} - u_n) - (u_n - u_{n-1}) = u_{n+1} - 2u_n + u_{n-1}

この例ではn を中央にしてn + 1 とn - 1 の間で差分をとっているので、2階中央(中心)差分という。差分商としてはh2 で割る。

作用素(演算子)[編集]

上記の各差分値をそれぞれΔu (x ) 、∇u (x ) 、δu (x ) 、δ2u (x) と書く。これらの関数の前につけた記号を差分作用素(差分演算子、difference operator)という。これらは微分作用素と同様、関数から関数への写像である(ただしこれらは同じ記号でもナブラ変分作用素ではないので、混同しないよう注意)。Δを前進差分演算子forward difference operator)、∇を後退差分演算子backward difference operator)、δを中心(中央)差分演算子central difference operator)、D を微分演算子(derivative operator)と呼ぶ。

差分のテイラー展開から、前進差分作用素Δは形式的には次のように、微分作用素(微分演算子) D によるテイラー展開として表され、差分作用素の有限部分が微分作用素であるとみることができる:

 \Delta = hD + \frac12 h^2D^2 + \frac1{3!} h^3D^3 + \cdots = \mathrm{e}^{hD} - 1

また逆に微分作用素D

 hD = \log(1+\Delta) = \Delta - \frac12 \Delta^2 + \frac13 \Delta^3 + \cdots \,

と書ける。

また後退差分∇と中央差分δに対してはそれぞれ

 hD = -\log(1-\nabla)
 hD = \, \operatorname{arcsinh} \left( \delta \right)

となる。

前進差分の誤差は(u が連続微分可能として):

 \frac{\Delta u(x)}{h} - u'(x) = O(h) \quad (h \to 0)

後退差分に対しても同じ式が成り立つ:

 \frac{\nabla u(x)}{h} - f'(x) = O(h)

中央差分はさらに正確な近似を与える。(u が2回連続微分可能として)誤差は間隔の2乗に比例する:

 \frac{\delta u(x)}{h} - u'(x) =  O(h^{2}) \!

例 熱伝導方程式[編集]

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

 U_t=U_{xx} \,

左辺は時間t による微分、右辺は座標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} =\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/h^2 拡散数と呼ばれる)である。

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

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

誤差は時間ステップ数と空間ステップ数の2乗とに比例する:

 \Delta u = O(k)+O(h^2) \,

陰解法[編集]

時刻t_{n+1} に後退差分を用い、空間点 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} \,

が得られる。これを陰解法という。

線形方程式系:

 (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}\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(h4 ) でなくO(h2 ) となることが多い。

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

参考文献[編集]

  1. ^ Joel H. Ferziger; Milovan Perić; 小林敏雄、谷口伸行、坪倉誠訳 『コンピュータによる流体力学』 シュプリンガー・フェアラーク東京、2003年、36頁。ISBN 4-431-70842-1 

関連項目[編集]

外部リンク[編集]