差分法

提供: testwiki
ナビゲーションに移動 検索に移動

テンプレート:Otheruses テンプレート:出典の明記 テンプレート:Differential equations テンプレート:計算物理学 数値解析における有限差分法(ゆうげんさぶんほう、テンプレート:Lang-en-short)あるいは単に差分法は、微分方程式を解くために微分有限差分近似(差分商)で置き換えて得られる差分方程式で近似するという離散化手法を用いる数値解法である。18世紀にオイラーが考案したと言われる[1]

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

精度と誤差

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

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

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

f(xi)=f(x0+ih).

一般に注目すべきはテンプレート:仮リンクで、典型的にはこれを O-記法で表す。局所打ち切り誤差は、各点における誤差について言うもので、真値 テンプレート:Math と近似値 テンプレート:Mvar との差

f(xi)f'i

である。この誤差の評価には、テイラー展開の剰余項を見るのが簡便である。式 テンプレート:Math に対するテイラー展開のラグランジュ型剰余項

Rn(x0+h)=f(n+1)(ξ)(n+1)!(h)n+1(x0<ξ<x0+h)

から、局所打ち切り誤差の支配項が求められる。例えば、一階差分近似 (テンプレート:Math) を考えれば

f(x0+ih)=f(x0)+f(x0)ih+f(ξ)2!(ih)2

である。この右辺は有限差分法で得られる近似値である。一方、0階差分近似(n=0)を考えれば

f(x0+ih)=f(x0)+f(x0)ih

よって、0階差分近似での支配的な誤差は

f(ξ)2!(ih)2

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

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

簡単な例

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

u(x)=3u(x)+2

これを解くには、差分商

u(x+h)u(x)hu(x)(h0)

を用いて

u(x+h)=u(x)+h(3u(x)+2)

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

例 熱伝導方程式

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

Ut=Uxx

左辺は時刻tによる微分、右辺は座標xによる2階微分である。また、境界条件および初期条件は以下とする:

U(0,t)=U(1,t)=0(境界条件)
U(x,0)=U0(x)(初期条件)

これを数値的に解く1つの方法は、すべての微分を差分で近似することである。空間の領域をメッシュx0,,xJで、時間の領域をメッシュt0,,tNで分割しよう。どちらの分割も等間隔とし、空間点の間隔をh、時刻の間隔をkとする。U(xj,tn)の数値的近似をujnで表す。

陽解法

時刻tnには前進差分を用い、空間点xjで2次微分に対して2次中央差分を用いれば、次の漸化式:

ujn+1ujnkΔt=uj+1n2ujn+uj1nh2

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

ujn+1の値は次のように得られる:

ujn+1=(12r)ujn+ruj1n+ruj+1n

ただしここでr=kΔt/h2拡散数と呼ばれる)である。

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

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

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

Δu=O(k)+O(h2)

陰解法

時刻tn+1に後退差分を用い、空間点xj で2階中央差分を用いれば、漸化式:

ujn+1ujnkΔt=uj+1n+12ujn+1+uj1n+1h2

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

線形方程式系:

(1+2r)ujn+1ruj1n+1ruj+1n+1=ujn

を解けば、ujn+1が得られる。この方法は常に数値的に安定で収束するが、時刻ごとに方程式系を解く必要があるため、陽解法よりも繁雑である。誤差は時間ステップ数と空間ステップ数の4乗とに比例する。

クランク・ニコルソン法

さいごに、時刻tn+1/2で中央差分を、空間点xjでの空間微分に2階中央差分を用いれば、漸化式:

2(ujn+1ujnkΔt)=uj+1n+12ujn+1+uj1n+1h2+uj+1n2ujn+uj1nh2

が得られる。これをクランク・ニコルソン法(Crank-Nicolson method)という。

線形方程式系:

(2+2r)ujn+1ruj1n+1ruj+1n+1=(22r)ujn+ruj1n+ruj+1n

を解けば、ujn+1が得られる。

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

Δu=O(k2)+O(h4)

しかし、境界付近では誤差はO(h4 ) でなくO(h2 ) となることが多い。

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

参考文献

テンプレート:Reflist

関連文献

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

関連項目

外部リンク

テンプレート:偏微分方程式の数値解法

テンプレート:Authority control

  1. テンプレート:Cite
  2. テンプレート:Cite book
  3. 高橋亮一、棚田芳弘:「差分法:数値シミュレーションの基礎」、培風館、ISBN 978-4-56303378-1 (1991年7月)
  4. テンプレート:Cite book
  5. 5.0 5.1 テンプレート:Cite book
  6. 6.0 6.1 テンプレート:Cite journal
  7. テンプレート:Cite book
  8. テンプレート:Cite book
  9. 矢嶋信男、野木達夫:「発展方程式の数値解析」、岩波書店(1977年9月26日)