ガウス・ジョルダン法
ガウス・ジョルダン法はガウスの消去法同様広く利用されており、連立方程式を手動で解く場合にもよく利用されている。
ガウス・ジョルダン法はしばしば逆行列(元の行列にかけると単位行列になる様な行列。単位行列は対角成分が全て1で後は0の正方行列)を求める際にも用いられる。
今回はコンピュータプログラムへの実装も視野にいれつつ、より詳細にかつ線型代数学の手法をもちいてガウス・ジョルダン法について述べる。
ガウス・ジョルダン法は別名掃出し法とも呼ばれ拡大係数行列の係数行列部分の対角成分以外を0に掃き出すためこの呼び名で呼ばれることがある。
対角成分以外を0にする作業を行列の対角化といい、ガウス・ジョルダン法の特徴はこの対角化にあるといえる。
対角化
対角化は行列の対角成分以外を0にすることであるが、対角化された行列を対角行列と呼ぶ。
対角行列は左上から右下へ行列を対角線上に値を持つものを言う。右上から左下の場合は対角行列とは呼ばない
\[
\left[
\begin{array}{}
a_{11} & 0 & 0\\
0 &a_{22} & 0\\
0 & 0 &a_{33}\\
\end{array}
\right]
\]
対角行列
\[
\left[
\begin{array}{}
0 & 0 & a_{13}\\
0 &a_{22} & 0\\
a_{31} & 0 & 0\\
\end{array}
\right]
\]
対角行列ではない
ガウス・ジョルダン法で連立方程式を解く場合、最終的に拡大係数行列の係数行列部分での対角成分の値を全て1にする。
\[
\left[
\begin{array}{ccc|c}
1 & 0 & 0 &b_{1}\\
0 & 1 & 0 &b_{2}\\
0 & 0 & 1 &b_{3}\\
\end{array}
\right]
\]
拡大係数行列の係数行列の対角化は行基本変形を使って行う。
ガウス・ジョルダン法による連立方程式の解法
以下の様なn元1次連立方程式を解くことを考える
\[
\left\{
\begin{array}{1}
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 \\
a_{n1}x_{1} + a_{n2}x_{2} + \cdots = a_{nn}x_{n} =b_{n}
\end{array}
\right.
\]
まず以下の様に拡大係数行列変換する。そして対角化をすることで連立方程式の解を得る。
\[
\left[
\begin{array}{cccc|c}
a_{11} &a_{12} & \cdots &a_{1n} &b_{1}\\
a_{21} &a_{22} & \cdots &a_{2n} &b_{2}\\
\vdots & \vdots & \ddots & \vdots & \vdots \\
a_{n1} &a_{n2} & \cdots &a_{nn} &b_{n}\\
\end{array}
\right]
\]
まず左上端の対角成分\(a_{11}\)を1にしたいので、1行目に対して\(a_{11}\)で割る。
\[
\left[
\begin{array}{cccc|c}
1 &a_{12}' & \cdots &a_{1n}' &b_{1}'\\
a_{21} &a_{22} & \cdots &a_{2n} &b_{2}\\
\vdots & \vdots & \ddots & \vdots & \vdots \\
a_{n1} &a_{n2} & \cdots &a_{nn} &b_{n}\\
\end{array}
\right]
\]
次に対角成分の下の列成分を0にしたいので1行目を\(-a_{21}\)倍して2行目に加える。
すると2行目の行の先頭が0になる
\[
\left[
\begin{array}{cccc|c}
1 &a_{12}' & \cdots &a_{1n}' &b_{1}'\\
0 &a_{22}' & \cdots &a_{2n}' &b_{2}'\\
\vdots & \vdots & \ddots & \vdots & \vdots \\
a_{n1} &a_{n2} & \cdots &a_{nn} &b_{n}\\
\end{array}
\right]
\]
これを列全体に大して行い、1列目の対角成分以外を全て0にする。
\[
\left[
\begin{array}{cccc|c}
1 &a_{12}' & \cdots &a_{1n}' &b_{1}'\\
0 &a_{22}' & \cdots &a_{2n}' &b_{2}'\\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 &a_{n2}' & \cdots &a_{nn}' &b_{n}'\\
\end{array}
\right]
\]
これで1行目に対する操作は終了である。
次は2列目に対して対角成分\(a_{22}’\)に着目して1列目と同様自身で割って1にして、
2行目以外の全ての行の列成分に対してその負数をかけてその行に加える。
\[
\left[
\begin{array}{cccc|c}
1 &0 & \cdots &a_{1n}' &b_{1}'\\
0 &1 & \cdots &a_{2n}' &b_{2}'\\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots &a_{nn}' &b_{n}'\\
\end{array}
\right]
\]
この操作を全ての対角成分で行い、最終的に以下の形になる
\[
\left[
\begin{array}{cccc|c}
1 &0 & \cdots & 0 &b_{1}'\\
0 &1 & \cdots & 0 &b_{2}'\\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 1 &b_{n}'\\
\end{array}
\right]
\]
この時の定数項ベクトル成分\(b_{1}’ ,b_{2}’, \dots b_{n}’\)が連立方程式の解になる。
具体例として以下の4元1次連立方程式についてガウス・ジョルダン法で解く
\[
\left\{
\begin{array}{1}
x_{1} + x_{2} + x_{3} + x_{4} =10\\
2x_{1} + x_{2} + 3x_{3} + 2x_{4} =21\\
x_{1} + 3x_{2} + 2x_{3} + x_{4} =17\\
3x_{1} + 2x_{2} + x_{3} + x_{4} =14\\
\end{array}
\right.
\]
まず拡大係数係数行列に変換する
\[
\left[
\begin{array}{cccc|c}
1 & 1 & 1 & 1 & 10\\
2 & 1 & 3 & 2 & 21\\
1 & 3 & 2 & 1 & 17\\
3 & 2 & 1 & 1 & 14\\
\end{array}
\right]
\]
- 1行目/1
- 2行目−1行目×2
- 3行目−1行目
- 4行目−1行目×3
\[
\left[
\begin{array}{cccc|c}
1 & 1 & 1 & 1 & 10\\
0 & -1& 1 & 0 & 1\\
0 & 2 & 1 & 0 & 7\\
0 & -1&-2 &-2 & -16\\
\end{array}
\right]
\]
- 1行目−2行目×−1
- 2行目/−1
- 3行目−2行目×−2
- 4行目−1行目×−1
\[
\left[
\begin{array}{cccc|c}
1 & 0 & 2 & 1 & 11\\
0 & 1 & -1 & 0 & -1\\
0 & 0 & 3 & 0 & 9\\
0 & 0 &-3 &-2 & -17\\
\end{array}
\right]
\]
- 1行目−3行目×2/3
- 2行目−3行目×−1/3
- 3行目/3
- 4行目−1行目
\[
\left[
\begin{array}{cccc|c}
1 & 0 & 0 & 1 & 5\\
0 & 1 & 0 & 0 & 2\\
0 & 0 & 1 & 0 & 3\\
0 & 0 & 0 &-2 & -8\\
\end{array}
\right]
\]
- 1行目−4行目×−1/2
- 2行目−4行目×0
- 3行目−4行目×0
- 4行目/−2
\[
\left[
\begin{array}{cccc|c}
1 & 0 & 0 & 0 & 1\\
0 & 1 & 0 & 0 & 2\\
0 & 0 & 1 & 0 & 3\\
0 & 0 & 0 & 1 & 4\\
\end{array}
\right]
\]
\[
\therefore x_1=1,x_2=2,x_3=3,x_4=4
\]
以上がガウス・ジョルダン法である。
ピボット選択
プログラムの実装においておおよその計算結果は同じだが、係数行列の対角項の絶対値が大きい時、計算誤差が小さくなることが知られている。
そのため対角項の絶対値が大きくなる様に行の並び替えを行う。この操作をピボット選択という。
コード
for(int i=0;i<n;i++) {
double inv=1.0/a[i][i];
for(int k=i;k<m;k++) {
a[i][k]*=inv;
}
for(int j=0;j<n;j++) {
if(i==j)continue;
double alpha=a[j][i];
for(int k=i;k<m;k++) {
a[j][k]-=alpha*a[i][k];
}
}
}
サンプルコード
ピボット選択付きのガウス・ジョルダン法
//ピボット選択
for(int i=0;i<n;i++){
int diag=i;
double max=Math.abs(a[i][i]);
for(int j=i+1;j<n;j++) {
if(i==j)continue;
if(max<Math.abs(a[j][i])) {
max=Math.abs(a[j][i]);
diag=j;
}
}
if(i==diag)continue;
for(int j=0;j<m;j++) {
double work=a[i][j];
a[i][j]=a[diag][j];
a[diag][j]=work;
}
}
サンプルコード(ガウス・ジョルダン法)