前回の授業では、連立一次方程式 $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法について説明します。
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)$ は既に求まっているため、問題なく計算できます。
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法も収束しない場合がありますが、収束の十分条件を表す定理として以下が知られています。
SOR法のアルゴリズムは、次のようにまとめられます。
入力:行列 $A$、ベクトル $b$、初期値ベクトル $x^{(0)}$、微小正数 $EPS$、反復回数の上限 $K$、パラメータ $w$
出力:反復回数 $k$ と近似解 $x$、または「反復回数の上限を超えた」旨のメッセージ
Step 1:$A$、$b$、$EPS$、$K$、$w$ を初期化する。$x=x^{(0)}$ とする。
Step 2:$n=$($A$の行数)とする。$k=0$ とする。
Step 3:$\left\| b-Ax\right\|_{\infty}\geq EPS$ かつ $k\leq K$ の間、次のStep 3-1~3-2を繰り返す。
Step 3-1:$i=1,\ldots,n$ について、$x(i)$ の更新を繰り返す。
Step 3-2:$k$ に $1$ を加算する。
Step 4:もし $k>K$ なら、「反復回数の上限を超えた」旨のメッセージを出力する。そうでなければ、反復回数 $k$ と近似解 $x$ を出力する。
赤字の箇所は敢えて曖昧に書いていますので、どうしたらよいか考えてください。
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 $$とします。
%演習1のコード
%ガウス・ザイデル法のコード(前回の演習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$ の値と反復回数を表示してください。ただし、最適値の判定はプログラムの中で行うようにすること。
%演習2のコード
%for文でwの値を変化させ、変数min_wとmin_kを用意して反復回数kの最小値およびそのときのwの値を求めるとよい
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$ の最適値および描画した図について自由に考察してください。
%演習3のコード
(Markdownとして考察を書く)
演習1~演習3に取り組んでください。