{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2022年度プログラミング演習A・B\n",
"\n",
"# 第10回レポート課題の解説"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習1\n",
"\n",
"SOR法のアルゴリズムを実現するコードを書き、\n",
"\n",
"$$\n",
"A=\n",
"\\begin{pmatrix}\n",
"2 & -1 & 0 & 0\\\\\n",
"-1 & 2 & -1 & 0\\\\\n",
"0 & -1 & 2 & -1\\\\\n",
"0 & 0 & -1 & 2\n",
"\\end{pmatrix}\n",
",\\quad b=\n",
"\\begin{pmatrix}\n",
"1\\\\\n",
"1\\\\\n",
"1\\\\\n",
"1\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"に対して $Ax=b$ の近似解を求めてください。ただし、\n",
"\n",
"$$\n",
"x^{(0)}=\n",
"\\begin{pmatrix}\n",
"0\\\\\n",
"0\\\\\n",
"0\\\\\n",
"0\n",
"\\end{pmatrix}\n",
",\\quad EPS=10^{-5},\\quad K=100,\\quad w=1.5\n",
"$$\n",
"\n",
"とします。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"19回の反復で、次の近似解を求めることができた。\n",
"1.999996142430643\n",
"2.999994924546839\n",
"2.999998348368432\n",
"2.000001772195130\n"
]
}
],
"source": [
"A = [2,-1,0,0; -1,2,-1,0; 0,-1,2,-1; 0,0,-1,2];\n",
"b = [1;1;1;1];\n",
"EPS = 1E-5;\n",
"K = 100;\n",
"w = 1.5;\n",
"x = [0;0;0;0];\n",
"\n",
"n = size(A,1);\n",
"k = 0;\n",
"\n",
"while norm(b-A*x,Inf) >= EPS && k <= K\n",
" for i = 1:n\n",
" tmp = b(i);\n",
" for j = 1:i-1\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" for j = i+1:n\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" x(i) = (1-w)*x(i)+w*tmp/A(i,i); %wを用いてx(i)を更新する\n",
" end\n",
" k = k+1;\n",
"end\n",
"\n",
"if k > K\n",
" printf(\"%d回の反復では、近似解を求めることができなかった。\\n\",K)\n",
"else\n",
" printf(\"%d回の反復で、次の近似解を求めることができた。\\n\",k)\n",
" printf(\"%.15f\\n\",x)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習2\n",
"\n",
"演習1と同様に\n",
"\n",
"$$\n",
"A=\n",
"\\begin{pmatrix}\n",
"2 & -1 & 0 & 0\\\\\n",
"-1 & 2 & -1 & 0\\\\\n",
"0 & -1 & 2 & -1\\\\\n",
"0 & 0 & -1 & 2\n",
"\\end{pmatrix}\n",
",\\quad b=\n",
"\\begin{pmatrix}\n",
"1\\\\\n",
"1\\\\\n",
"1\\\\\n",
"1\n",
"\\end{pmatrix},\n",
"$$\n",
"\n",
"$$\n",
"x^{(0)}=\n",
"\\begin{pmatrix}\n",
"0\\\\\n",
"0\\\\\n",
"0\\\\\n",
"0\n",
"\\end{pmatrix}\n",
",\\quad EPS=10^{-5},\\quad K=100\n",
"$$\n",
"\n",
"に対してSOR法を適用する際に、\n",
"\n",
"$$\n",
"w=1,1.1,1.2,\\ldots,1.8,1.9\n",
"$$\n",
"\n",
"の中から最も収束の速い(反復回数の少ない)ものを求め、その $w$ の値と反復回数を表示してください。ただし、最適値の判定はプログラムの中で行うようにすること。"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"wの最適値は1.300000、反復回数は11\n"
]
}
],
"source": [
"A = [2,-1,0,0; -1,2,-1,0; 0,-1,2,-1; 0,0,-1,2];\n",
"b = [1;1;1;1];\n",
"EPS = 1E-5;\n",
"K = 100;\n",
"n = size(A,1);\n",
"\n",
"min_w = 0; %反復回数が最小値をとるようなwの候補\n",
"min_k = K+1; %反復回数の最小値の候補(初期値を十分大きな値にしておく)\n",
"\n",
"for w = 1:0.1:1.9 %wを変化させてSOR法を繰り返す\n",
" x = [0;0;0;0]; %変数xとkはSOR法を実行すると値が変わるため、for文の中で初期化する必要がある\n",
" k = 0;\n",
" while norm(b-A*x,Inf) >= EPS && k <= K\n",
" for i = 1:n\n",
" tmp = b(i);\n",
" for j = 1:i-1\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" for j = i+1:n\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" x(i) = (1-w)*x(i)+w*tmp/A(i,i);\n",
" end\n",
" k = k+1;\n",
" end\n",
" if k < min_k %反復回数が最小値の候補より小さければ、候補を更新する\n",
" min_w = w;\n",
" min_k = k;\n",
" end\n",
"end\n",
"\n",
"printf(\"wの最適値は%f、反復回数は%d\\n\",min_w,min_k)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習3(オプション)\n",
"\n",
"SOR法による収束の様子を可視化して確認します。\n",
"\n",
"$$\n",
"A=\n",
"\\begin{pmatrix}\n",
"2 & -1\\\\\n",
"1 & 2\n",
"\\end{pmatrix}\n",
",\\quad b=\n",
"\\begin{pmatrix}\n",
"1\\\\\n",
"3\n",
"\\end{pmatrix},\n",
"$$\n",
"\n",
"$$\n",
"x^{(0)}=\n",
"\\begin{pmatrix}\n",
"0\\\\\n",
"0\n",
"\\end{pmatrix}\n",
",\\quad EPS=10^{-2},\\quad K=100\n",
"$$\n",
"\n",
"に対してSOR法を適用する際のパラメータ $w$ の最適値(なるべく良い値)を求めた上で、連立一次方程式\n",
"\n",
"$$\n",
"\\left\\{\n",
"\\begin{matrix}\n",
"2x_1-x_2 = 1\\\\\n",
"x_1+2x_2 = 3\n",
"\\end{matrix}\n",
"\\right.\n",
"$$\n",
"\n",
"における各方程式が表す直線のグラフ(ただし、範囲は $0\\leq x_1\\leq 1.4$)と、求めた $w$ に対するSOR法により得られるベクトル列 $\\left\\{x^{(k)}\\right\\}$ の各項が表す点を、一つの図に描画してください。\n",
"\n",
"さらに、求めた $w$ の最適値および描画した図について自由に考察してください。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"wの最適値は0.910000、反復回数は3\n"
]
}
],
"source": [
"%wの最適値を求める部分は、演習2のコードを利用すればよい\n",
"A = [2,-1; 1,2];\n",
"b = [1;3];\n",
"EPS = 1E-2;\n",
"K = 100;\n",
"n = size(A,1);\n",
"\n",
"min_w = 0;\n",
"min_k = K+1;\n",
"\n",
"for w = 0.01:0.01:1.99 %wを0.01から1.99まで0.01刻みで変化させる\n",
" x = [0;0];\n",
" k = 0;\n",
" while norm(b-A*x,Inf) >= EPS && k <= K\n",
" for i = 1:n\n",
" tmp = b(i);\n",
" for j = 1:i-1\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" for j = i+1:n\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" x(i) = (1-w)*x(i)+w*tmp/A(i,i);\n",
" end\n",
" k = k+1;\n",
" end\n",
" if k < min_k\n",
" min_w = w;\n",
" min_k = k;\n",
" end\n",
"end\n",
"\n",
"printf(\"wの最適値は%f、反復回数は%d\\n\",min_w,min_k)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3回の反復で、次の近似解を求めることができた。\n",
"1.003797556409375\n",
"0.998616877139984\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"%上の行はコメントではなく、図をきれいにする(ベクター画像にする)ためのもの\n",
"\n",
"%図の描画は、前回の演習3のコードを利用すればよい\n",
"x_list = 0:0.01:1.4;\n",
"y1_list = 2*x_list-1; %一つ目の方程式に対応する直線の式\n",
"y2_list = -x_list/2+3/2; %二つ目の方程式に対応する直線の式\n",
"\n",
"plot(x_list,y1_list) %一つ目の直線の描画\n",
"grid on %グリッド線の表示\n",
"hold on %複数のplotを行うために必要\n",
"plot(x_list,y2_list) %二つ目の直線の描画\n",
"\n",
"w = min_w; %wを上で求めた最適値にする\n",
"x = [0;0];\n",
"plot(x(1),x(2)) %初期値ベクトルの点の描画\n",
"k = 0;\n",
"\n",
"while norm(b-A*x,Inf) >= EPS && k <= K\n",
" for i = 1:n\n",
" tmp = b(i);\n",
" for j = 1:i-1\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" for j = i+1:n\n",
" tmp = tmp-A(i,j)*x(j);\n",
" end\n",
" x(i) = (1-w)*x(i)+w*tmp/A(i,i);\n",
" end\n",
" plot(x(1),x(2)) %新たに計算したベクトルの点の描画\n",
" k = k+1;\n",
"end\n",
"\n",
"if k > K\n",
" printf(\"%d回の反復では、近似解を求めることができなかった。\\n\",K)\n",
"else\n",
" printf(\"%d回の反復で、次の近似解を求めることができた。\\n\",k)\n",
" printf(\"%.15f\\n\",x)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### (考察)\n",
"\n",
"SOR法では多くの問題に対して $1