ヤコビ法

提供: miniwiki
2018/8/19/ (日) 17:42時点におけるAdmin (トーク | 投稿記録)による版 (1版 をインポートしました)
(差分) ← 古い版 | 最新版 (差分) | 新しい版 → (差分)
移動先:案内検索

ヤコビ法とは[math]n[/math]元の連立一次方程式[math]A\vec{x}=\vec{b}[/math]反復法で解く手法の1つである。ドイツ数学者カール・グスタフ・ヤコブ・ヤコビの名前にちなむ。

[math]n[/math]正方行列[math]A[/math]は、上三角行列[math]U[/math]、下三角行列[math]L[/math]対角行列[math]D[/math]とすると、[math]A=L+D+U[/math]と書ける。このようにすると、まず以下のような変形ができる。

[math] \begin{array}{ccc} (L+D+U) \vec{x} &=& \vec{b} \\ D\vec{x} &=& \vec{b}-(L+U)\vec{x} \\ \end{array} [/math]

この式を満たす[math]\ x[/math]を求める。初期値[math]\vec{x}^{(0)}[/math]に対して、 [math]\ k[/math]回目の反復で得られた[math]x_1[/math]の値を[math]x_1^{(k)}[/math]と書くと、 以下のような反復法の漸化式ができる。

[math] D\vec{x}^{(k+1)} = \vec{b}-(L+U)\vec{x}^{(k)} [/math]

この式は以下のように変形できる。

[math] \vec{x}^{(k+1)} = D^{-1} \{ \vec{b}-(L+U)\vec{x}^{(k)} \} [/math]

もし、解が収束した場合、その場合は[math]x_1^{(k+1)}[/math][math]x_1^{(k)}[/math]は共通の値[math]x_1^{(*)}[/math]を持つことになる。このとき、

[math] \vec{x}^{(*)} = D^{-1} \{ \vec{b}-(L+U)\vec{x}^{(*)} \} [/math]

となり、変形していくと元の連立方程式の形に戻る。 したがって、ヤコビ法で解が収束した場合、その解は連立方程式の解となる。 また、その収束の十分条件は、係数行列の対角要素の絶対値が非対角要素の絶対値よりも相対的に大きい場合、すなわち対角優位な行列である場合に収束する。これはガウス=ザイデル法も同様である。

ヤコビ法の式はベクトル[math]\vec{x}[/math]の各成分ごとに次のような式で書くことができ、数値解析ではこの式が用いられる。

[math] x_{i}^{(k+1)} = \frac{1}{a_{ii}} \left( b_{i} - \sum_{j=1}^{i-1} a_{ij}x_{j}^{(k)} - \sum_{j=i+1}^{n} a_{ij}x_{j}^{(k)} \right) = \frac{1}{a_{ii}} \left( b_{i} - \sum_{j \neq i}^{n} a_{ij}x_{j}^{(k)} \right) [/math]

ガウス=ザイデル法とヤコビ法を加速する方法としてはSOR法が知られている。

具体例

3元の連立一次方程式、すなわち、

[math]\left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \end{array} \right)[/math]

を解くことを考える。[math]k[/math]回目の反復で得られた[math]x_1[/math]の値を[math]x_1^{(k)}[/math]と書く。 初期値[math]\vec{x}^{(0)}[/math]は、適当な値、例えばゼロベクトルでもかまわない。

[math]x_1^{(k+1)} = (b_1 - a_{12} x_2^{(k)} - a_{13} x_3^{(k)}) / a_{11}[/math]

[math]x_2^{(k+1)} = (b_2 - a_{21} x_1^{(k)} - a_{23} x_3^{(k)}) / a_{22}[/math]

[math]x_3^{(k+1)} = (b_3 - a_{31} x_1^{(k)} - a_{32} x_2^{(k)}) / a_{33}[/math]

という反復を繰り返していく。 ヤコビ法は、直列計算ではガウス=ザイデル法よりも遅いが、ガウス=ザイデル法と異なり各式が他の式に依存せず並列性があるため並列計算でも用いられる。

固有値問題

実対称行列の固有値および固有ベクトルを求める繰り返し計算手法においてもヤコビ法と呼ばれる解法がある。 [math]\ n[/math]次の実対称行列[math]\ A[/math]について次のように [math]\ G(p, q, \theta)[/math] による相似変換、すなわちギブンス回転を実行することにより、非対角要素 [math]\ a_{ij}(i \ne j)[/math] の最大値 [math]\ a_{pq}[/math] が0となるようにする。

[math]\ B=G^TAG[/math]

これによって行列 [math]\ B[/math] の各要素は次のようになる。但し、 [math]\ i,j \ne p,q [/math] である。

[math]\begin{cases} b_{pp} = a_{pp}\cos^2\theta + a_{qq}\sin^2\theta - 2a_{pq}\sin\theta\cos\theta ,\\ b_{qq} = a_{pp}\sin^2\theta + a_{qq}\cos^2\theta + 2a_{pq}\sin\theta\cos\theta ,\\ b_{pq} = b_{qp} = \frac{1}{2}(a_{pp}-a_{qq})\sin2\theta + a_{pq}\cos2\theta ,\\ b_{ij} = a_{ij} ,\\ b_{pj} = a_{pj}\cos\theta - a_{qj}\sin\theta ,\\ b_{qj} = a_{pj}\sin\theta + a_{qj}\cos\theta ,\\ b_{ip} = a_{ip}\cos\theta - a_{iq}\sin\theta ,\\ b_{iq} = a_{ip}\sin\theta + a_{iq}\cos\theta \end{cases}[/math]

ここで、[math]\ a_{pq} \ne 0[/math] のとき [math]\ b_{pq}=0[/math] となる [math]\ \theta[/math] は上式より

[math] \tan 2\theta = \frac{-2a_{pq}}{(a_{pp}-a_{qq})}[/math]

から求められることがわかる。 ギブンス回転をすべての非対角要素がほぼ0になるまで繰り返せば、実対称行列 [math]A\quad[/math] が対角化された形となるから、その対角要素が [math]A\quad[/math] の固有値となる。また、 [math]A\quad[/math] がk回変換された行列を [math]A_{k}\quad[/math] 、k回目のギブンス回転を表す直交行列を [math]G_{k}\quad[/math] と表せば、

[math] A_k = G_k^T A_{k-1} G_k = U_k^TAU_k[/math]
ここに、 [math] U_k = G_1 G_2 G_3 \cdots G_{k-1} G_k[/math]

となる。 [math]A_{k}\quad[/math] のすべての非対角要素がほぼ0となったとき、 [math]U_{k}\quad[/math] は固有ベクトルを並べた行列となっている。 なお、ギブンス回転の繰り返し過程において、一度は0になった要素がその後の変換により0でなくなることもあるが、変換の繰り返しによって非対角項は0に近づいてゆく。

関連項目