今回から4回分の授業(ただし、間にまとめの回を挟む)では、「連立一次方程式の解を求める」という問題を扱います。
変数の個数と式の本数がともに $n$ であるような連立一次方程式
$$ \left\{ \begin{matrix} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n & = & b_1\\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n & = & b_2\\ \vdots & \vdots & \vdots\\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n & = & b_n \end{matrix} \right. $$を考えます。このとき、
$$ 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$ が正則(可逆)である、すなわち逆行列 $A^{-1}$ が存在する場合には、連立一次方程式は唯一の解
$$ x=A^{-1}b $$をもちます。
例えば
$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} 11\\ -5\\ -5\\ 4 \end{pmatrix} $$のときには、行列式 $\det A\neq 0$ より $A$ は正則であり、$Ax=b$ の唯一の解をMATLABでは次のように演算子\
を用いて求めることができます。
A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1] %各行の成分をカンマで区切り、行をセミコロンで区切る
b = [11;-5;-5;4] %縦ベクトル
n = size(A,1) %Aの行数 size(A)(1)としても同じ
det(A) %Aの行列式の値
x = A\b %演算子「\」を用いてAx=bの解を求める
A*x-b %理論上はゼロベクトルになるが、実際には計算誤差の影響で少しずれる
A = 2 1 -2 3 2 0 4 1 3 0 5 2 1 -1 2 1 b = 11 -5 -5 4 n = 4 ans = -1.00000 x = 25.000 -14.000 -10.000 -15.000 ans = 0.0000e+00 3.5527e-15 0.0000e+00 -3.5527e-15
ただし、演算子\
によって行われる計算のアルゴリズムは複雑であり、処理の内容が基本的にユーザーには分かりません。
参考URL:https://jp.mathworks.com/help/matlab/ref/mldivide.html
この授業では上記の演算子を使用して満足するのではなく、連立一次方程式の解を求めるための具体的なアルゴリズムを学習します。
特に、今回と次回の授業では直接法(direct method)に分類されるアルゴリズムによって連立一次方程式の厳密解を求め、まとめの回を挟んでその後の二回の授業では反復法(iterative method)に分類されるアルゴリズムによって連立一次方程式の近似解を求めます。
ガウスの消去法(Gaussian elimination)は、連立一次方程式の解を求めるための直接法の有名なアルゴリズムです。
ガウスの消去法は、前進消去(forward elimination)と後退代入(backward substitution)という二つの過程に分けられます。
連立一次方程式 $Ax=b$ を考えるとき、ガウスの消去法の前進消去では、拡大係数行列 $(A\,|\,b)$ を行列の行基本変形によって $\hat{A}$ が上三角行列であるような $(\hat{A}\,|\,\hat{b})$ に変形します($\hat{A}$ の対角成分を $1$ にする必要はありません)。
ここで、行列の行基本変形は
二つの行を入れ替える
ある行を $0$ でないスカラー倍する
ある行に別の行のスカラー倍を加える
の三つであり、ガウスの消去法の前進消去では行基本変形の1と3を用います。
特に、前進消去の過程で着目する対角成分に $0$ が現れない場合には、後ほど説明するピボット選択は必要なく行基本変形の3のみを用いればよいので、ここではまずそのような場合に限定して考えることにします。
先ほどの $A$ と $b$ を例として、前進消去を進めます。拡大係数行列を格納する変数を用意する方法もありますが、ここでは $A$ と $b$ を別に考えてそれぞれの成分を更新していきます。前進消去が終わった段階では、変数 $A$ と $b$ がそれぞれ上記の $\hat{A}$ と $\hat{b}$ の役割になります。
A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1]
b = [11;-5;-5;4]
n = size(A,1)
A = 2 1 -2 3 2 0 4 1 3 0 5 2 1 -1 2 1 b = 11 -5 -5 4 n = 4
まず、$A$ の第 $1$ 行を用いて、第 $1$ 列の対角より下の成分が $0$ になるように行基本変形3を行います。$b$ も合わせて変形します。
alpha = -A(2,1)/A(1,1) %Aの(2,1)成分を0にするためのスカラー
A(2,:) = A(2,:)+alpha*A(1,:) %Aの第2行に第1行のalpha倍を加える
b(2) = b(2)+alpha*b(1) %bの第2成分に第1成分のalpha倍を加える
alpha = -A(3,1)/A(1,1) %Aの(3,1)成分を0にするためのスカラー
A(3,:) = A(3,:)+alpha*A(1,:) %Aの第3行に第1行のalpha倍を加える
b(3) = b(3)+alpha*b(1) %bの第3成分に第1成分のalpha倍を加える
alpha = -A(4,1)/A(1,1) %Aの(4,1)成分を0にするためのスカラー
A(4,:) = A(4,:)+alpha*A(1,:) %Aの第4行に第1行のalpha倍を加える
b(4) = b(4)+alpha*b(1) %bの第4成分に第1成分のalpha倍を加える
alpha = -1 A = 2 1 -2 3 0 -1 6 -2 3 0 5 2 1 -1 2 1 b = 11 -16 -5 4 alpha = -1.5000 A = 2.00000 1.00000 -2.00000 3.00000 0.00000 -1.00000 6.00000 -2.00000 0.00000 -1.50000 8.00000 -2.50000 1.00000 -1.00000 2.00000 1.00000 b = 11.0000 -16.0000 -21.5000 4.0000 alpha = -0.50000 A = 2.00000 1.00000 -2.00000 3.00000 0.00000 -1.00000 6.00000 -2.00000 0.00000 -1.50000 8.00000 -2.50000 0.00000 -1.50000 3.00000 -0.50000 b = 11.0000 -16.0000 -21.5000 -1.5000
次に、$A$ の第 $2$ 行を用いて、第 $2$ 列の対角より下の成分が $0$ になるように行基本変形3を行います。$b$ も合わせて変形します。
alpha = -A(3,2)/A(2,2) %Aの(3,2)成分を0にするためのスカラー
A(3,:) = A(3,:)+alpha*A(2,:) %Aの第3行に第2行のalpha倍を加える
b(3) = b(3)+alpha*b(2) %bの第3成分に第2成分のalpha倍を加える
alpha = -A(4,2)/A(2,2) %Aの(4,2)成分を0にするためのスカラー
A(4,:) = A(4,:)+alpha*A(2,:) %Aの第4行に第2行のalpha倍を加える
b(4) = b(4)+alpha*b(2) %bの第4成分に第2成分のalpha倍を加える
alpha = -1.5000 A = 2.00000 1.00000 -2.00000 3.00000 0.00000 -1.00000 6.00000 -2.00000 0.00000 0.00000 -1.00000 0.50000 0.00000 -1.50000 3.00000 -0.50000 b = 11.0000 -16.0000 2.5000 -1.5000 alpha = -1.5000 A = 2.00000 1.00000 -2.00000 3.00000 0.00000 -1.00000 6.00000 -2.00000 0.00000 0.00000 -1.00000 0.50000 0.00000 0.00000 -6.00000 2.50000 b = 11.0000 -16.0000 2.5000 22.5000
さらに、$A$ の第 $3$ 行を用いて、第 $3$ 列の対角より下の成分が $0$ になるように行基本変形3を行います。$b$ も合わせて変形します。
alpha = -A(4,3)/A(3,3) %Aの(4,3)成分を0にするためのスカラー
A(4,:) = A(4,:)+alpha*A(3,:) %Aの第4行に第3行のalpha倍を加える
b(4) = b(4)+alpha*b(3) %bの第4成分に第3成分のalpha倍を加える
alpha = -6 A = 2.00000 1.00000 -2.00000 3.00000 0.00000 -1.00000 6.00000 -2.00000 0.00000 0.00000 -1.00000 0.50000 0.00000 0.00000 0.00000 -0.50000 b = 11.0000 -16.0000 2.5000 7.5000
$A$ を上三角行列に変形することができており、これで前進消去は終了です。
ガウスの消去法における前進消去のアルゴリズムを実現するコードを書き、
$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} 11\\ -5\\ -5\\ 4 \end{pmatrix} $$に対する変形の結果を表示してください。ただし、$A$ が $4$ 次以外の正方行列の場合にも正しく動作するコードにすること。
%演習1のコード
A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1];
b = [11;-5;-5;4];
n = size(A,1);
%ここに、不足しているコードを書く
%for文のネスト(for文の中にfor文を入れる)を使うべき
%4ではなく変数nを使うことが重要
ガウス消去法の前進消去によって、元の $A$ と $b$ から次のような $\hat{A}$ と $\hat{b}$ が得られています。
$$ \hat{A}= \begin{pmatrix} \hat{a}_{11} & \hat{a}_{12} & \cdots & \hat{a}_{1(n-1)} & \hat{a}_{1n}\\ 0 & \hat{a}_{22} & \cdots & \hat{a}_{2(n-1)} & \hat{a}_{2n}\\ \vdots & \ddots & \ddots & \vdots & \vdots\\ 0 & \cdots & 0 & \hat{a}_{(n-1)(n-1)} & \hat{a}_{(n-1)n}\\ 0 & \cdots & 0 & 0 & \hat{a}_{nn} \end{pmatrix} ,\quad \hat{b}= \begin{pmatrix} \hat{b}_1\\ \hat{b}_2\\ \vdots\\ \hat{b}_{n-1}\\ \hat{b}_n \end{pmatrix} $$ここで、行列の行基本変形の性質から、変形前の $Ax=b$ と変形後の $\hat{A}x=\hat{b}$ の解は同じです。
$\hat{A}$ は上三角行列であるため、次のようにすれば解 $x$ の成分 $x_i$ を後ろから順に求めることができます。
$x_n=\hat{b}_n \bigm/ \hat{a}_{nn}$
$x_{n-1}=\left(\hat{b}_{n-1}-\hat{a}_{(n-1)n}x_n\right) \bigm/ \hat{a}_{(n-1)(n-1)}$
$\vdots$
$x_2=\left(\hat{b}_{2}-\hat{a}_{23}x_3-\cdots-\hat{a}_{2n}x_n\right) \bigm/ \hat{a}_{22}$
$x_1=\left(\hat{b}_{1}-\hat{a}_{12}x_2-\cdots-\hat{a}_{1n}x_n\right) \bigm/ \hat{a}_{11}$
演習1のコードにガウスの消去法における後退代入のアルゴリズムを実現するコードを追加し、
$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} 11\\ -5\\ -5\\ 4 \end{pmatrix} $$に対して $Ax=b$ の厳密解を求めてください。ただし、$A$ が $4$ 次以外の正方行列の場合にも正しく動作するコードにすること。
%演習2のコード
ガウスの消去法の前進消去において着目する対角成分に $0$ が現れる場合には、その列の対角より下の成分を $0$ にすることができません。そのような場合には、着目する対角成分が $0$ でなくなるように行基本変形1の行の入れ替えを行います。
特に、着目している対角成分より下の成分のうち、絶対値の最も大きな成分を利用することが望ましいです。これを、ガウスの消去法においてピボット選択(pivoting)を行うと言います。
ピボット選択を行うガウスの消去法のアルゴリズムを実現するコードを書き、ピボット選択が必要な $A$ と $b$ の例を与えて $Ax=b$ の厳密解を求めてください。
%演習3のコード
演習1~演習3に取り組んでください。