第7-3回　反復法（jacobi法、gauss-seidel法）[python, julia, fortran]

• English (United States)
• 日本語

第7-3回では数値解析における逆行列の解法を説明したいと思います。

1．数値解析における逆行列の解法

そもそも，第7回では陰解法で解こうという趣旨で進めています。そして陰解法では次のように式変形できることが前回までの流れで分かったかと思います。
$$\bf A \bf x = \bf b$$

ここで\(\bf A,\bf b\)はすでに分かっており，\(\bf x\)を求めたいという状況です。

このような行列を解く方法として，「直接解法」，「定常反復法」，「非定常反復法」の3つが挙げられます。

①直接解法
有限回の演算で解が求まる方法です。第2-2回の「ガウスの消去法」や第3-2回の「TDMA法」は直接解法になります。
確実に答えがでるのですが，記憶容量が大量に必要であり演算に時間がかかります。
TDMA法のように要素の位置が決まっており，計算が簡略出来る場合には使えるかもしれませんが，この方法はあまり使われません。

②定常反復法
ある初期ベクトル\(\bf x^{(0)} \)から真の解\(\bf x^*\)に収束するよう反復計算を行う方法。その内，収束率が一定の方法。

③非定常反復法
定常反復法と同様に反復計算を行う方法。その内，収束率が一定でないの方法。

反復法は収束条件をみたさない場合があるが（数値解析では基本的に大丈夫），収束までの演算量は直接解法より少なく数値解析の主流はこちらです。

２．反復法を用いて逆行列を解く

ヤコビ法とガウス・サイデル法の比較

今回は定常反復法のヤコビ法とガウス・ザイデル法を紹介します。

最初は何をしているか意味がわかりにくいですが、一通り見てみてください。

前回の陰解法一次元拡散方程式を使います。離散化した形は下記の通りです。

$$\rho C \frac{\Delta T^n_i}{\Delta t}=\frac{1}{\Delta x^2}(\lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2}(\Delta T_{i}^n -\Delta T_{i-1}^n))+R^n_i$$

ただし、$$\Delta T^n_i =T^{n+1}_i-T^n_i, \ \ \ R^n_i=\frac{1}{\Delta x^2}(\lambda{i+1/2} (T_{i+1}^n-T_i^n)-\lambda_{i-1/2}(T_{i}^n -T_{i-1}^n)) $$

\(\Delta T_i^n\)を左辺に寄せる式変形をします．

$$ \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2} }{\Delta x^2} + \frac{ \lambda_{i-1/2} }{\Delta x^2} \right ) \Delta T^n_i =\frac{1}{\Delta x^2} \left ( \lambda_{i+1/2} \Delta T_{i+1}^n-\lambda_{i-1/2} \Delta T_{i-1}^n \right )+R^n_i$$

$$ \Delta T^n_i =\left [ \frac{1}{\Delta x^2} \left ( \lambda_{i+1/2} \Delta T_{i+1}^n-\lambda_{i-1/2} \Delta T_{i-1}^n \right )+R^n_i \right ] / \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2} }{\Delta x^2} + \frac{ \lambda_{i-1/2} }{\Delta x^2} \right ) $$

ここで反復数\(k\)を用いて，次のように書きます。

$$ \Delta T^{k+1}_i =\left [ \frac{1}{\Delta x^2} \left ( \lambda_{i+1/2} \Delta T_{i+1}^k-\lambda_{i-1/2} \Delta T_{i-1}^k \right )+R^n_i \right ] / \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2} }{\Delta x^2} + \frac{ \lambda_{i-1/2} }{\Delta x^2} \right ) $$

上式の計算を適当な初期値、例えば\( \Delta T_i=0, \ \ k=0 \)から始め、\(k\)を増やしていきます。十分変化しないくらい収束したら計算が終わりです。

混乱しないでほしいのですが、収束した結果は次の時間ステップの解です。

つまり、ある時間\(n\)の解が分かっており\(k=0\)から\(k=k^*\)まで計算し収束した場合、\(\Delta T_{i}^{n+1}=\Delta T_{i}^{k^*}\)とみなすわけです。

また、反復するべき値は\(\bf A \bf x =\bf b \)の求めたい解\(\bf x\)だけなので、上式の\(R^n_i\)は\(n\)のままです。

これがヤコビ法です。

証明について簡単に書いておきます。

$$\bf A \bf x=\bf b$$
$$⇔(\bf L+\bf D+\bf U) \bf x=\bf b$$
$$⇔\bf D \bf x=-(\bf L+\bf U)\bf x+\bf b$$
$$⇔\bf x=-\bf D^{-1}(\bf L+\bf U)\bf x+\bf D \bf b$$

次の通り\(\bf M,\bf N\)で置き換えます。また真の解に\(*\)を付けます。 $$⇔\bf x^{*}=\bf M \bf x^{*}+\bf N \ \ \ \ \ \ ①$$

ヤコビ法は次のように書けます。
$$⇔\bf x^{k+1}=\bf M \bf x^{k}+\bf N \ \ \ \ \ \ ②$$

(①-②)より
$$⇔\bf x^{*}-\bf x^{k+1}=\bf M (\bf x^{*}-\bf x^{k})$$
$$⇔\bf x^{*}-\bf x^{k+1}=\bf M ({\bf M \bf x^{*}+\bf N}-{\bf M \bf x^{k-1}+\bf N })$$
$$⇔\bf x^{*}-\bf x^{k+1}=\bf M ^2 (\bf x^{*}-\bf x^{k-1})$$
$$ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots $$
$$⇔\bf x^{*}-\bf x^{k+1}=\bf M ^{k+1} (\bf x^{*}-\bf x^{0})$$

両辺で2-ノルムを取ります。

2-ノルムは座標空間上における原点からの距離を一般化したものと考えてください。

具体例としてはこんな感じです。
$$\bf x = [1 \ \ 2] ⇔ ||\bf x||=\sqrt{1^2+2^2}=\sqrt{5}$$

では適用します。
$$⇔||\bf x^{*}-\bf x^{k+1}||=||\bf M ||^{k+1} ||\bf x^{*}-\bf x^{0}||$$

\(||\bf M ||<1\)であれば、\(k→\infty\)のとき\(||\bf M ||^{k+1}→0\)より
$$\bf x^{*}=\bf x^{k+1}$$

\(||\bf M ||<1\)は対角優位のとき成り立ちます。

深くは掘り下げませんが、拡散方程式の場合、次の時間におけるセル\(i\)の情報は現在の時間のセル\(i\)の情報が一番効いてきます。そのため自動的に対角優位の行列になるため、この反復法が使えます。

ヤコビ法は反復数を更新する際、\(k\)を用いて\(k+1\)のすべてを計算し、更新します。

しかし計算は順番に行うので、セル\(i+1\)の更新の際はセル\(i\)の計算はすでに終わっているわけです。

だから、セル\((i+1,k+1)\)の計算にセル\((i,k+1)\)を用います。
これがガウス・ザイデル法で、ヤコビ法より収束が速いです。

収束条件は同様の対角優位になります。

メリット、デメリットとしては以下のことがよく言われます。

ヤコビ法：並列化が簡単、収束性が悪い
ガウス・ザイデル法：実装が楽、収束性が良い

現在でもガウス・ザイデル法を応用したものが数値解析界隈では使用されています。

次回は二次元拡散方程式に対し、式を整理し実装できる形に持っていきます。

*コードは第7-5回にあります。

by hide