5.2 連立一次方程式 (定常反復法)

n 元連立一次方程式 Ax = b に対して次のような形式の反復法を定常反復解法という. 上付きの (k) は k 回目の反復を表す.
\[
x^{(k+1)} = M x^{(k)} + N b
\] ここで, M と N は反復回数 k に依存しない行列である. M は反復行列と呼ばれる.

この反復が収束するための条件は反復行列 M のスペクトル半径 (行列の固有値の最大絶対値) \rho(M) が 1 より小さいことである.
\[
\rho(M) = max|\lambda_i| < 1 \space (\lambda_i は行列 M の固有値, i = 1, \dots, n) \] 代表的な定常反復解法としては, ヤコビ法, ガウス・ザイデル法, SOR (逐次過剰緩和) 法がある. SOR 法は実用的にもよく使用されている.

5.2.1 ヤコビ法

n 元連立一次方程式 Ax = b の第 i 番目の式
\[
\sum_{j=1}^n a_{ij} x_j = b_i
\] を \(x_i\) について解き, その値を新たな \(x_i\) とする反復を \(i = 1, \dots, n\) の全てについて行う.
\[
x_i^{(k+1)} = (b_i – \sum_{j \ne i, j=1}^n a_{ij} x_j^{(k)}) / a_{ii}
\] いま, A を左下三角行列 L (対角要素は 0), 対角行列 D, 右上三角行列 U (対角要素は 0)に分離して, A = L + D + U とする. そうすると上の反復は次のように表すことができる.
\[
x^{(k+1)} = D^{-1}(b – (L + U) x^{(k)}) = -D^{-1}(L + U) x^{(k)} + D^{-1}b
\] すなわち, 反復行列 M は次のように表される.
\[
M = -D^{-1}(L + U)
\] 次式が全ての i について成り立つときに行列 A は狭義の対角優位であると言う.
\[
|a_{ii}| \ge \sum_{j \ne i, j=1}^n |a_{ij}|
\] A が狭義の対角優位であるとき, ヤコビ法の反復行列 M のスペクトル半径 \(\rho\)(M) は 1 より小さい. すなわち, ヤコビ法は収束する.

5.2.2 ガウス・ザイデル法

ヤコビ法と似ているが, \(i = 1, \dots, n\) の全部の計算を待たずに \(x_i\) をすぐに更新して現在の近似解とする. すなわち, \(i = 1, \dots, n\) の順に更新する.
\[
x_i^{(k+1)} = (b_i – \sum_{j \lt i} a_{ij} x_j^{(k+1)} – \sum_{i \lt j} a_{ij} x_j^{(k)}) / a_{ii}
\] 行列表示では次のようになる.
\[
\begin{align}
x^{(k+1)} & = D^{-1}(b – Lx^{(k+1)} – Ux^{(k)}) \\
& = -(D + L)^{-1} Ux^{(k)} + (D + L)^{-1}b \\
\end{align}
\] 反復行列 M は次のように表される.
\[
M = -(D + L)^{-1}U
\]

A が狭義の対角優位であるとき, ガウス・ザイデル法の反復行列 M のスペクトル半径 \(\rho\)(M) は 1 より小さく, ガウス・ザイデル法は収束する.

この方法は, 式で書くとヤコビ法より面倒に思えるが, x をすぐに書き換えてもよいのでプログラムを書くのはこちらの方が楽である.

5.2.3 SOR (逐次過剰緩和) 法

ガウス・ザイデル法の反復における本来の更新量を \(\omega\) 倍して加速する方法である. \(\omega\) は加速パラメータと呼ばれる.

行列表示では次のように表される.
\[
x^{(k+1)} = (I + \omega D^{-1}L)^{-1} ((1 – \omega)I – \omega D^{-1}U) x^{(k)} + \omega(D + \omega L)^{-1}b
\] 反復行列 M は次のように表される.
\[
M = (I + \omega D^{-1}L)^{-1} ((1 – \omega)I – \omega D^{-1}U)
\] \(\omega\) の最適値を前もって知ることは難しく, 経験的に決められることが多い. 最適値を使えば SOR 法はガウス・ザイデル法よりかなり速く収束する.

A が狭義の対角優位行列であるとき, \(\omega\) の値が \(0 < \omega < 2\) であれば SOR 法は収束する.

5.2.3.1 SOR 法の ω の最適値

スペクトル半径 \(\rho\)(M) は \(\omega\) により変化する. 5.2.4.3 数値実験 (1) の例題において, 反復行列 M の固有値を求めてノルムが最大のものすなわちスペクトル半径 ρ(M) をプロットする. 横軸は ω, 縦軸はスペクトル半径 ρ(M) である.

尖った底になっている点は ω の最適値で 1.53 であった. ω が 2 を超えるとスペクトル半径が 1 より大きくなる (収束しなくなる).

ω の最適値は次式から求めることができる (文献[6]参照).
\[
\omega_{opt} = \frac{2}{1 + (1 – \mu_0^2)^{1/2}}
\] ここで \(\mu_0\) はヤコビ法の反復行列のスペクトル半径である. この場合, \(\mu_0\) = 0.951 であるから, 上式に代入すると \(\omega_{opt}\) = 1.53 となり実験結果と一致する.

なお, ω = 1, すなわちガウス・ザイデル法の場合は ρ(M) = 0.9045 であったが, これはヤコビ法の反復行列の ρ(M) (= \(\mu_0\)) = 0.951 の 2 乗に一致する. また, \(\omega_{opt}\) を超えると ρ(M) と ω の関係は直線になり, ρ(M) = ω – 1 になる.

5.2.3.2 適応的 SOR 法

SOR 法において ω の最適な値を前もって知ることは難しいが, 反復中に収束率を見積り ω の値を自動的に調節する適応的プログラムも研究されている. このように, ω の値を自動的に最適値に調節する機能を持った SOR 法を適応的 SOR 法 (Adaptive SOR method) という.

一例として, SOR 法の反復中に得られる情報を使ってスペクトル半径を見積もることにより ω の値を自動的に調節する方法が提案されている (詳細は文献[10]参照). この方法は係数行列が対称正定値行列であることなど適用できる問題に制限があるがうまく動作するようである.

5.2.4 定常反復法の比較

例としてポアソン方程式の境界値問題を解く.

5.2.4.1 ポアソン方程式の境界値問題の例

2 次元の正方形領域 (0 ≤ x ≤ 1, 0 ≤ y ≤ 1) におけるポアソン方程式を考える.
\[
-\Delta u = f
\] ここで, \(\Delta u = \nabla \cdot(\nabla u) = \partial^2u/\partial x^2 + \partial^2u/\partial y^2\) である.

境界条件は次のとおりである.
\[
\begin{align}
& u(x, 0) = 0, \space u(x, 1) = 0 \space (0 \le x \le 1) \\
& u(0, y) = 0, \space u(1, y) = 0 \space (0 \le y \le 1) \\
\end{align}
\] すなわち, 正方形領域の境界上では 0 とする.

5.2.4.2 5 点差分近似による離散化

偏微分方程式の微分を差分に置き換えて差分方程式を作り, それを解くことにより偏微分方程式の近似解を求める方法が差分法である.

2 階微分の差分近似は次のように表される.
\[
\begin{align}
& \partial^2 u/\partial x^2 = (u(x-h, y) – 2u(x, y) + u(x+h, y)) / h^2 \\
& \partial^2 u/\partial y^2 = (u(x, y-h) – 2u(x, y) + u(x, y+h)) / h^2 \\
\end{align}
\] ここで, h は小さな値の刻みであるが, 正方形の定義域の x, y 方向それぞれを N 等分して格子を作りその各辺を刻みとして使うことにすると, この場合 h = 1 / N になる. さらに, 格子点での値を \(u_{ij} = u(i/N, j/N) \space (i = 0, \dots, N, \space j = 0, \dots, N)\) と表すと, 格子点 (i, j) での差分近似は次のようになる.
\[
\begin{align}
& \partial^2 u/\partial x^2 = (u_{i-1,j} – 2u_{i,j} + u_{i+1,j}) / h^2 \\
& \partial^2 u/\partial y^2 = (u_{i,j-1} – 2u_{i,j} + u_{i,j+1}) / h^2 \\
\end{align}
\] これをポアソン方程式に代入すると次の近似式が得られる.
\[
u_{i,j} = (u_{i-1,j} + u_{i+1,j} + u{i,j-1} + u{i,j+1})/4
\] 図で表すと, 下のように前後左右の値の平均値で真ん中の値を近似していることがわかる. これは 5 点差分近似とよばれる.

ここで N = 10 とすると, 格子点は 11×11 = 121 点になる. そのうち境界上の 40 点 (一番外側の 4 辺上の点) は境界条件で与えられているから, 未知の \(u_{i,j}\) は内側の 81 点 (i = 1, …, 9, j = 1, …, 9) である. また, 上の近似式もこれらの 81 点について定義されるから 81 本ある. 従って, この近似式からなる 81 元連立一次方程式を解くことにより全ての \(u_{i,j}\) を求めることができる. すなわち, 差分法により偏微分方程式の境界値問題は連立一次方程式 (差分方程式) に帰着された.

こうして得られた差分方程式の係数行列は, 変数を \(u_{i,j}\) ((i,j) = (1,1), (2,1), …, (9,1), (1,2), (2,2), …, (9,9)) の順に並べると, 下に示すパターンの行列になる. 青で表されているのが非ゼロ要素 (この場合は 369 個で, 全体の 5.6% にあたる) で, その他の要素はすべて 0 の疎行列である. また, この場合, 正定値対称行列である.

差分法の場合, 精度を上げるために格子を細かくしていくと差分方程式は N の 2 乗に比例して大規模になっていく. この例では, N = 20 では 361 元, N = 100 では 9801 元というようにすぐに大きくなる.

5.2.4.3 数値実験 (1) 差分方程式を解く

ポアソン方程式の例題の差分方程式を解く. ここでは N = 10 (81 元連立一次方程式) とし, 解が (1, 1, …, 1) となるように右辺を設定した.

係数行列 A の対角要素は全て 4 であるから \(|a_ii|\) = 4 である. そして, 各行の対角要素以外の非ゼロ要素の値は -1 でその個数は最大で 4 個であるから, \(\sum |a_ij| \le 4\) である. 従って, この係数行列 A は狭義の対角優位である. すなわち, ヤコビ法, ガウス・ザイデル法, SOR 法 (0 < ω < 2) は収束する. ヤコビ法と ω = 1, 1.43 ～ 1.73 に変化させた SOR 法 (ω = 1 のときガウス・ザイデル法に一致) で計算してみる. 初期値を (0, 0, ..., 0) としたときの収束の様子を次に示す. 横軸は反復回数, 縦軸は残差ノルム (右辺ベクトルのノルムで正規化したもの) である.

図のようにヤコビ法よりガウス・ザイデル法の方が速いことがわかる.

SOR 法では, この場合は ω = 1 から ω が増加するに伴って収束が速くなっていくが, 最適値 (この場合は ω = 1.53)を超えると逆に遅くなっていくのがわかる. 最適値を超えると収束が不安定になる傾向がみえる.

5.2.4.4 数値実験 (2) 適応的 SOR 法により差分方程式を解く

上と同じポアソン方程式の例題に適応的 SOR 法を適用した数値実験結果を示す.

比較のため, ω = 1.53 (最適値) および 1.63 の場合も示した. この例では, 最適な ω を使用した場合と遜色ない結果を示した.