Part.21

熱伝導方程式の陰公式

(Part.20)で、1次元熱伝導方程式の数値解法について解説した。そこで用いた方法は陽公式と呼ばれ、安定性の条件

α = Dτ/h2 ≦1/2 (S)

を満たす必要がありました。この条件については、計算量の点で問題がある。例えば、時刻0から1まで計算を行うとする。時間方向の講師の総数は、 1/τ である。計算の手間を減らす為には τ をなるべく大きくとりたい。そこで、(S)から、

τ≒h2 /2

とすればよいだろう。一方、空間に関する分解能も高めたい。そこで、h を 1/10 倍小さくしたとする。すると、τ は上の関係より 1/100 倍小さくする必要がある。つまり、時間方向の格子点が100倍に増える。これでは、h を小さくしていったときに、あまりに格子点数が多くなってしまう。

この問題は、陽公式の安定性条件(S)から生じる。ここに紹介する陰公式には、(S)のような条件ない。

陰公式

(Part.20)の(1.5)式の右辺に現れる k を全て k+1 で置き換える。

(j = 1,2,...,N-1, k = 0,1,2,....)

この式をα = Dτ/h2 として整理すると、

となる。両端で 0 の境界条件では、任意の k に対して、

となる。上式を行列の形で書き直すと、以下の

に関する連立一次方程式が得られる。

このとき、

が既知であれば、この連立一次方程式を解くことによって、

を求める事ができる。この左辺の行列は、三重対角とよばれる特別な形をしており、効率よく解く事ができる。シミュレーションにおいては、陽公式では、漸化式を用いて

を求めたが、その部分が連立一次方程式を解く事に置き換わるだけである。

具体的には、(Part.20)のアルゴリズムの、

j := 1,2,....,N の順に
  new_uj := αuj-1 + (1 - 2α)uj + αuj-1
を繰り返す

の部分が連立一次方程式を解くアルゴリズムで置き換わる。

連立一次方程式の解法については、直接法と呼ばれる、行列を直接変形して求める方法(ガウスの消去法、LU分解法等)と、反復法と呼ばれる、適当な漸化式の反復計算から解の近似値を求める方法(ヤコビ法、ガウスーザイデル法、SOR法、CG法およびその亜種等)がある。

空間2次元以上については、反復法を用いることが多い。現象に関わる偏微分方程式をシミュレーションする場合、最終的には要素の殆どが0であるような超巨大な行列からなる連立一次方程式(未知数が数億とか数百億)を解く事となり、行列の性質をうまく使って反復回数を減らす等の工夫がなされる。現在も様々な計算方法が開発され続けている。

フーリエ分解の方法を使えば、上記陰公式は無条件安定、すなわち、どのような τ, h の値に対しても計算が安定に進む事がわかる。

課題1

フーリエ分解の方法を用いて、陰公式が無条件安定である事を示せ。

課題2

空間2次元熱伝導方程式についても同様に陰公式を得よ。その際出てくる連立一次方程式の左辺の行列は5重対角行列と呼ばれる特別な形をしている。

課題3

空間1次元熱伝導方程式を陰公式を用いてシミュレーションせよ。

課題4

課題2で得られた陰公式を、ヤコビ法(反復法の一種で最も単純)を用いて解く事によって、空間2次元熱伝導方程式のシミュレーションを完成せよ。