第6回と第7回の授業では、
$$ A= \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} ,\quad x= \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{pmatrix} ,\quad b= \begin{pmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{pmatrix} $$として連立一次方程式
$$ Ax=b $$を与え、$A$ が正則（可逆）であると仮定した上で直接法（ガウスの消去法とLU分解）によってその厳密解（ただし、計算誤差は含まれうる）を求めることを考えました。
今回と次回の授業では、反復法に分類されるアルゴリズムによって $Ax=b$ の近似解を求めることを考えます。
正則行列 $M$ と正方行列 $N$ を用いて
$$ A=M+N $$と表せたとすると、$Ax=b$ は
$$ Mx=b-Nx $$と変形できます。初期値ベクトル $x^{(0)}$ が与えられるとき、$M$ が正則であることから漸化式
$$ Mx^{(k)}=b-Nx^{(k-1)} $$によってベクトル列 $\left\{x^{(k)}\right\}_{k=0}^{\infty}$ を定めることができます。このとき、もし $x^{(k)}\to\bar{x}$ と収束するならば、
$$ M\bar{x}=b-N\bar{x} $$より極限 $\bar{x}$ は $Ax=b$ の厳密解になることが分かります。実際には、漸化式による反復計算を有限回で打ち切るため、最後に求めた $x^{(k)}$ が近似解として得られることになります。
以下では、$M$ と $N$ の異なる決め方に基づく三つの反復法を扱います（そのうちの一つは次回の授業で扱います）。
準備として、$A$ の対角成分からなる行列を $D$、$A$ の対角より下の成分からなる行列を $L$、$A$ の対角より上の成分からなる行列を $U$ とし、
$$ A=D+L+U $$と表しておきます。具体的に成分で書けば、
$$ D= \begin{pmatrix} a_{11} & 0 & \cdots & 0\\ 0 & a_{22} & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & a_{nn} \end{pmatrix} ,\quad L= \begin{pmatrix} 0 & 0 & \cdots & 0\\ a_{21} & 0 & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ a_{n1} & \cdots & a_{n(n-1)} & 0 \end{pmatrix} ,\quad U= \begin{pmatrix} 0 & a_{12} & \cdots & a_{1n}\\ 0 & 0 & \ddots & \vdots\\ \vdots & \ddots & \ddots & a_{(n-1)n}\\ 0 & \cdots & 0 & 0 \end{pmatrix} $$です。
また、以下では $A$ の全ての対角成分が $0$ でないこと、すなわち
$$ \forall i\in\{1,\ldots,n\},\; a_{ii}\neq 0 $$を仮定します（そうでない場合は、仮定が満たされるように $A$ と $b$ に対して行の入れ替えを行えばよいです）。
ヤコビ法（Jacobi method）においては、
$$ M=D,\quad N=L+U $$とします。ここで、全ての対角成分が $0$ でない対角行列 $D$ は正則であり、
$$ D^{-1}= \begin{pmatrix} \frac{1}{a_{11}} & 0 & \cdots & 0\\ 0 & \frac{1}{a_{22}} & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & \frac{1}{a_{nn}} \end{pmatrix} $$です。
このとき、前述の漸化式
$$ Mx^{(k)}=b-Nx^{(k-1)} $$は
$$ x^{(k)}=D^{-1}\left(b-Lx^{(k-1)}-Ux^{(k-1)}\right) $$と変形され、成分ごとに書けば $x^{(k)}$ の $i$ 番目の成分 $x^{(k)}_i$ について
$$ x^{(k)}_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j=1}^{i-1}a_{ij}x^{(k-1)}_j-\sum_{j=i+1}^{n}a_{ij}x^{(k-1)}_j\right) $$が得られます。
この漸化式を用いれば、ベクトル $x^{(k-1)}$ の成分を使って新たなベクトル $x^{(k)}$ の成分を計算することができます。
ヤコビ法を含む反復法において、漸化式を用いた反復計算を終了する条件（収束判定基準）は様々に考えられます（参考リンク）。
ここでは、$Ax^{(k)}$ が $\infty$ ノルムの意味で $b$ とほぼ等しくなったとき、すなわち、微小正数 $EPS$ に対して
$$ \left\| b-Ax^{(k)}\right\|_{\infty} < EPS $$を満足したときに、「ベクトル列が収束した」とみなして反復を終了することにします。
さらに、ベクトル列 $\left\{x^{(k)}\right\}_{k=0}^{\infty}$ が収束しない場合もあるため、反復回数が事前に決めた上限 $K$ を超えたときには「ベクトル列が収束しなかった」とみなして反復を終了することにします。
ヤコビ法のアルゴリズムは、次のようにまとめられます。
入力：行列 $A$、ベクトル $b$、初期値ベクトル $x^{(0)}$、微小正数 $EPS$、反復回数の上限 $K$
出力：反復回数 $k$ と近似解 $x$、または「反復回数の上限を超えた」旨のメッセージ
Step 1：$A$、$b$、$EPS$、$K$ を初期化する。$x=x^{(0)}$ とする。
Step 2：$n=$（$A$の行数）とする。$k=0$ とする。$x\_new$ を $n$ 次元の縦ベクトルとして初期化する。
Step 3：$\left\| b-Ax\right\|_{\infty}\geq EPS$ かつ $k\leq K$ の間、次のStep 3-1～3-3を繰り返す。
Step 3-1：$i=1,\ldots,n$ について、次のStep 3-1-1～3-1-4を繰り返す。
Step 3-1-1：$tmp=b(i)$ とする。
Step 3-1-2：$j=1,\ldots,i-1$ について、$tmp=tmp-A(i,j)*x(j)$ を繰り返す。
Step 3-1-3：$j=i+1,\ldots,n$ について、$tmp=tmp-A(i,j)*x(j)$ を繰り返す。
Step 3-1-4：$x\_new(i)=tmp\,/\,A(i,i)$ とする。
Step 3-2：$x=x\_new$ とする。
Step 3-3：$k$ に $1$ を加算する。
Step 4：もし $k>K$ なら、「反復回数の上限を超えた」旨のメッセージを出力する。そうでなければ、反復回数 $k$ と近似解 $x$ を出力する。
ヤコビ法のアルゴリズムを実現するコードを書き、
$$ 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 $$とします。
%演習1のコード
A = [2,-1,0,0; -1,2,-1,0; 0,-1,2,-1; 0,0,-1,2];
b = [1;1;1;1];
EPS = 1E-5;
K = 100;
x = [0;0;0;0];
n = size(A,1);
k = 0;
x_new = zeros(n,1);
%ここに、不足しているコードを書く
%無限大ノルムは「norm(ベクトル,Inf)」を使えばよい
if k > K
printf("%d回の反復では、近似解を求めることができなかった。\n",K)
else
printf("%d回の反復で、次の近似解を求めることができた。\n",k)
printf("%.15f\n",x)
end
ガウス・ザイデル法（Gauss-Seidel method）においては、
$$ M=D+L,\quad N=U $$とします。ここで、全ての対角成分が $0$ でない下三角行列 $D+L$ は正則です（ただし、その逆行列はあまり簡単には表現できません）。
このとき、前述の漸化式
$$ Mx^{(k)}=b-Nx^{(k-1)} $$は
$$ x^{(k)}=D^{-1}\left(b-Lx^{(k)}-Ux^{(k-1)}\right) $$と変形され、成分ごとに書けば $x^{(k)}$ の $i$ 番目の成分 $x^{(k)}_i$ について
$$ x^{(k)}_i=\frac{1}{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)$ は既に求まっているため、問題なく計算できます。
むしろ、ヤコビ法と違って $x^{(k)}_i$ の計算の時点で $x^{(k-1)}_j\;(j=1,\ldots,i-1)$ の情報は必要ないため、プログラム上では変数 $x\_new$ を用意せずに変数 $x$ の成分を単に前から順に更新していけばよいことになります。
ガウス・ザイデル法のアルゴリズムは、次のようにまとめられます。
入力：行列 $A$、ベクトル $b$、初期値ベクトル $x^{(0)}$、微小正数 $EPS$、反復回数の上限 $K$
出力：反復回数 $k$ と近似解 $x$、または「反復回数の上限を超えた」旨のメッセージ
Step 1：$A$、$b$、$EPS$、$K$ を初期化する。$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$ を出力する。
赤字の箇所は敢えて曖昧に書いていますので、どうしたらよいか考えてください。
ガウス・ザイデル法のアルゴリズムを実現するコードを書き、
$$ 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 $$とします。
%演習2のコード
ヤコビ法とガウス・ザイデル法による収束の様子を可視化して確認します。
連立一次方程式
$$ \left\{ \begin{matrix} 2x_1-x_2 = 1\\ x_1+2x_2 = 3 \end{matrix} \right. $$における各方程式が表す直線のグラフ（ただし、範囲は $0\leq x_1\leq 1.4$）と、
$$ 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 $$に対してヤコビ法を適用して得られるベクトル列 $\left\{x^{(k)}\right\}$ の各項が表す点を、一つの図に描画してください。
同様の描画を、ガウス・ザイデル法についても行ってください（つまり、結果は二つの図になります）。
さらに、描画した二つの図について自由に考察してください。
%plot --format svg
%演習3のコード1（上の行は、図をきれいにするためにセルの一行目に書く（%が必要））
%ここで、ヤコビ法についての図を描画する
%グラフの描画方法は、第4回授業資料の真ん中辺りを参照
%複数のplotをするために「hold on」が必要
%点の描画は「plot(x(1),x(2))」のようにすればよい
%plot --format svg
%演習3のコード2（上の行は、図をきれいにするためにセルの一行目に書く（%が必要））
%ここで、ガウス・ザイデル法についての図を描画する
%グラフの描画方法は、第4回授業資料の真ん中辺りを参照
%複数のplotをするために「hold on」が必要
%点の描画は「plot(x(1),x(2))」のようにすればよい
（Markdownとして考察を書く）
以下が知られています。
行列 $A$ が狭義対角優位（strictly diagonally dominant）、すなわち $$\forall i\in\{1,\ldots,n\},\; |a_{ii}|>\sum_{j=1\\ j\neq i}^n|a_{ij}|$$ ならば、ヤコビ法およびガウス・ザイデル法は収束する。
行列 $A$ が正定値かつ対称ならば、ガウス・ザイデル法は収束する。
これらは、あくまで十分条件にすぎません。実際、演習1と演習2で考えている行列 $A$ は正定値かつ対称であっても狭義対角優位ではありませんが、ヤコビ法とガウス・ザイデル法のどちらも収束します。
演習1～演習3に取り組んでください。