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

第6回:連立一次方程式の解の計算1

1. 導入

今回から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では次のように演算子\を用いて求めることができます。

ただし、演算子\によって行われる計算のアルゴリズムは複雑であり、処理の内容が基本的にユーザーには分かりません。

参考URL:https://jp.mathworks.com/help/matlab/ref/mldivide.html

この授業では上記の演算子を使用して満足するのではなく、連立一次方程式の解を求めるための具体的なアルゴリズムを学習します。

特に、今回と次回の授業では直接法(direct method)に分類されるアルゴリズムによって連立一次方程式の厳密解を求め、まとめの回を挟んでその後の二回の授業では反復法(iterative method)に分類されるアルゴリズムによって連立一次方程式の近似解を求めます。

2. ガウスの消去法

ガウスの消去法(Gaussian elimination)は、連立一次方程式の解を求めるための直接法の有名なアルゴリズムです。

ガウスの消去法は、前進消去(forward elimination)後退代入(backward substitution)という二つの過程に分けられます。

2.1. 前進消去

連立一次方程式 $Ax=b$ を考えるとき、ガウスの消去法の前進消去では、拡大係数行列 $(A\,|\,b)$ を行列の行基本変形によって $\hat{A}$ が上三角行列であるような $(\hat{A}\,|\,\hat{b})$ に変形します($\hat{A}$ の対角成分を $1$ にする必要はありません)。

ここで、行列の行基本変形は

  1. 二つの行を入れ替える

  2. ある行を $0$ でないスカラー倍する

  3. ある行に別の行のスカラー倍を加える

の三つであり、ガウスの消去法の前進消去では行基本変形の1と3を用います。

特に、前進消去の過程で着目する対角成分に $0$ が現れない場合には、後ほど説明するピボット選択は必要なく行基本変形の3のみを用いればよいので、ここではまずそのような場合に限定して考えることにします。

先ほどの $A$ と $b$ を例として、前進消去を進めます。拡大係数行列を格納する変数を用意する方法もありますが、ここでは $A$ と $b$ を別に考えてそれぞれの成分を更新していきます。前進消去が終わった段階では、変数 $A$ と $b$ がそれぞれ上記の $\hat{A}$ と $\hat{b}$ の役割になります。

まず、$A$ の第 $1$ 行を用いて、第 $1$ 列の対角より下の成分が $0$ になるように行基本変形3を行います。$b$ も合わせて変形します。

次に、$A$ の第 $2$ 行を用いて、第 $2$ 列の対角より下の成分が $0$ になるように行基本変形3を行います。$b$ も合わせて変形します。

さらに、$A$ の第 $3$ 行を用いて、第 $3$ 列の対角より下の成分が $0$ になるように行基本変形3を行います。$b$ も合わせて変形します。

$A$ を上三角行列に変形することができており、これで前進消去は終了です。

演習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} $$

に対する変形の結果を表示してください。ただし、$A$ が $4$ 次以外の正方行列の場合にも正しく動作するコードにすること。

2.2. 後退代入

ガウス消去法の前進消去によって、元の $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$ を後ろから順に求めることができます。

演習2

演習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.3. ピボット選択

ガウスの消去法の前進消去において着目する対角成分に $0$ が現れる場合には、その列の対角より下の成分を $0$ にすることができません。そのような場合には、着目する対角成分が $0$ でなくなるように行基本変形1の行の入れ替えを行います。

特に、着目している対角成分より下の成分のうち、絶対値の最も大きな成分を利用することが望ましいです。これを、ガウスの消去法においてピボット選択(pivoting)を行うと言います。

演習3(オプション)

ピボット選択を行うガウスの消去法のアルゴリズムを実現するコードを書き、ピボット選択が必要な $A$ と $b$ の例を与えて $Ax=b$ の厳密解を求めてください。

第6回レポート課題

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