2022年度プログラミング演習A・B

第10回:連立一次方程式の解の計算4

1. 前回の復習

前回の授業では、連立一次方程式 $Ax=b$ の近似解を求めるための反復法である、ヤコビ法とガウス・ザイデル法について扱いました。

具体的には、

$$ A=M+N $$

を満たす正則行列 $M$ と正方行列 $N$ を用いて初期値ベクトル $x^{(0)}$ と漸化式

$$ Mx^{(k)}=b-Nx^{(k-1)} $$

によってベクトル列 $\left\{x^{(k)}\right\}$ を生成することを考え、ヤコビ法においては

$$ M=D,\quad N=L+U, $$

ガウス・ザイデル法においては

$$ M=D+L,\quad N=U $$

としました。ただし、$D$ は $A$ の対角成分からなる行列、$L$ は $A$ の対角より下の成分からなる行列、$U$ は $A$ の対角より上の成分からなる行列であり、

$$ A=D+L+U $$

が成り立ちます。

今回の授業では、反復法の三つ目としてSOR法について説明します。

2. SOR法

2.1. SOR法の漸化式

SOR法(successive over-relaxation method、逐次過緩和法)においては、パラメータ $w\ (0<w<2)$ を用いて

$$ M=\frac{1}{w}D+L,\quad N=\frac{w-1}{w}D+U $$

とします。ここで、$M$ は全ての対角成分が $0$ でない下三角行列であるため、正則です。

このとき、前述の漸化式

$$ Mx^{(k)}=b-Nx^{(k-1)} $$

$$ \frac{1}{w}Dx^{(k)}=-\frac{w-1}{w}Dx^{(k-1)}+\left(b-Lx^{(k)}-Ux^{(k-1)}\right) $$

すなわち

$$ x^{(k)}=(1-w)x^{(k-1)}+wD^{-1}\left(b-Lx^{(k)}-Ux^{(k-1)}\right) $$

と変形され、成分ごとに書けば $x^{(k)}$ の $i$ 番目の成分 $x^{(k)}_i$ について

$$ x^{(k)}_i=(1-w)x^{(k-1)}_i+\frac{w}{a_{ii}}\left(b_i-\sum_{j=1}^{i-1}a_{ij}x^{(k)}_j-\sum_{j=i+1}^{n}a_{ij}x^{(k-1)}_j\right) $$

が得られます。

ガウス・ザイデル法と同様に、漸化式の右辺にも新たなベクトル $x^{(k)}$ の成分が含まれていますが、自然に $x^{(k)}_1,x^{(k)}_2,\ldots,x^{(k)}_n$ の順番で計算をする場合には $x^{(k)}_i$ の計算の時点で $x^{(k)}_j\;(j=1,\ldots,i-1)$ は既に求まっているため、問題なく計算できます。

2.2. SOR法の解釈

SOR法の漸化式のベクトル表現において

$$ \tilde{x}^{(k)}=D^{-1}\left(b-Lx^{(k)}-Ux^{(k-1)}\right) $$

とおくと、

$$ x^{(k)}=(1-w)x^{(k-1)}+w\tilde{x}^{(k)} $$

すなわち

$$ x^{(k)}=x^{(k-1)}+w\left(\tilde{x}^{(k)}-x^{(k-1)}\right) $$

が得られます。

ここで、$\tilde{x}^{(k)}$ はガウス・ザイデル法によって生成されるベクトルであり、$w=1$ のときにはSOR法はガウス・ザイデル法と一致することが分かります。

パラメータ $w$ はガウス・ザイデル法におけるベクトルの変化量を調節する役割を持ち、それによって収束を加速させることを狙います。多くの問題に対しては $1<w<2$ のときに収束が速くなりますが、$w$ の最適値を事前に知ることは難しく、問題ごとに良さそうな値を選択する必要があります。

なお、ヤコビ法およびガウス・ザイデル法と同様にSOR法も収束しない場合がありますが、収束の十分条件を表す定理として以下が知られています。

2.3. プログラムとしての実装

SOR法のアルゴリズムは、次のようにまとめられます。



赤字の箇所は敢えて曖昧に書いていますので、どうしたらよいか考えてください。

演習1

SOR法のアルゴリズムを実現するコードを書き、

$$ A= \begin{pmatrix} 2 & -1 & 0 & 0\\ -1 & 2 & -1 & 0\\ 0 & -1 & 2 & -1\\ 0 & 0 & -1 & 2 \end{pmatrix} ,\quad b= \begin{pmatrix} 1\\ 1\\ 1\\ 1 \end{pmatrix} $$

に対して $Ax=b$ の近似解を求めてください。ただし、

$$ x^{(0)}= \begin{pmatrix} 0\\ 0\\ 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-5},\quad K=100,\quad w=1.5 $$

とします。

演習2

演習1と同様に

$$ A= \begin{pmatrix} 2 & -1 & 0 & 0\\ -1 & 2 & -1 & 0\\ 0 & -1 & 2 & -1\\ 0 & 0 & -1 & 2 \end{pmatrix} ,\quad b= \begin{pmatrix} 1\\ 1\\ 1\\ 1 \end{pmatrix}, $$$$ x^{(0)}= \begin{pmatrix} 0\\ 0\\ 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-5},\quad K=100 $$

に対してSOR法を適用する際に、

$$ w=1,1.1,1.2,\ldots,1.8,1.9 $$

の中から最も収束の速い(反復回数の少ない)ものを求め、その $w$ の値と反復回数を表示してください。ただし、最適値の判定はプログラムの中で行うようにすること。

演習3(オプション)

SOR法による収束の様子を可視化して確認します。

$$ A= \begin{pmatrix} 2 & -1\\ 1 & 2 \end{pmatrix} ,\quad b= \begin{pmatrix} 1\\ 3 \end{pmatrix}, $$$$ x^{(0)}= \begin{pmatrix} 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-2},\quad K=100 $$

に対してSOR法を適用する際のパラメータ $w$ の最適値(なるべく良い値)を求めた上で、連立一次方程式

$$ \left\{ \begin{matrix} 2x_1-x_2 = 1\\ x_1+2x_2 = 3 \end{matrix} \right. $$

における各方程式が表す直線のグラフ(ただし、範囲は $0\leq x_1\leq 1.4$)と、求めた $w$ に対するSOR法により得られるベクトル列 $\left\{x^{(k)}\right\}$ の各項が表す点を、一つの図に描画してください。

さらに、求めた $w$ の最適値および描画した図について自由に考察してください。

(Markdownとして考察を書く)

第10回レポート課題

演習1~演習3に取り組んでください。