課題・ガウスの消去法3（ピボット）

ガウスの消去法1
ガウスの消去法2

引き続いてガウスの消去法。

前回はいったん答えが出たように思いましたが、このような

場合には正常な値が出ませんでした。

PRINTを前進消去の直後に入れると解るんですけど、
2行目以下の値がERRORになっています。

これは手計算してみると解るんですけど、
2行2列が0になってしまっています。

これによって2行目を2行2列で割って、1にするという計算が誤作動を起こします。
（定数割る0なので）

これを回避するために、ピボット（1にしたい対角要素）に合わせて
行を入れ替えるという方法があります。


↓



まず、このピボットの選択のためには今見ている行以下の行で、
見たい列の中の最大値がある行を現在の行と入れ替える方法を使います。

今の場合、見ている行は2行2列で、2行目以下（つまり2行目と3行目）で
2列目が最大になる行を選べばいいことになります。

ややこしいですけど2列目が一番大きい3行目を2行目と入れ替えるということです。


前進消去のループのしょっぱなには今見ている行列をbufに入れて、
bufで割ってしまっているのでピボットの選択はこれより前、
かつfor(i=0~)よりあとに入れます。


この式を例にとった場合、列の最大値をa_maxとでも置くと、

if(a_max > a_matrix[0][0]){
a_max = a_matrix[0][0];
}

if(a_max > a_matrix[1][0]){
a_max = a_matrix[1][0];
}
if(a_max > a_matrix[2][0]){
a_max = a_matrix[2][0];
}
//ここまで1行目(i=0)の場合

if(a_max > a_matrix[1][1]){
a_max = a_matrix[1][1];
}
if(a_max > a_matrix[2][1]){
a_max = a_matrix[2][1];
}
//ここまで2行目(i=1)の場合

if(a_max > a_matrix[2][2]){
a_max = a_matrix[2][2];
}
//ここまで3行目(i=2)の場合

次にこれを簡略化すると、jはこの時使用されていないのでjを使うと、

for(j = i;j < 3;j++){
if(a_max > a_matrix[j][i]){
a_max = a_matrix[j][i];
}
}

となります。これでその列の最大値が求まりましたが、
このままでは何行目にその最大値があったのかわかりません。
（最大値があった行と、もとの行を入れ替えるために必要）

なので適当な変数にその「最大値があった行」をいれておきます。
ここでは使われない(予定）のkをつかいます。
（できるだけ宣言する変数は少ない方がいいので、いままで宣言した変数の中で
その時値を変えても大丈夫な変数があればそれを使います。）

for(j = i;j < 3;j++){
if(a_max > a_matrix[j][i]){
a_max = a_matrix[j][i];
k = j;
}
}

これで最大値のある行が把握できました。

これからこのk行とi行とを入れ替えます。
入れ替えには1つ目を適当な変数に一度格納しておき、2つ目を1つ目に入れた後に
適当な変数から2つ目にいれるという方法を使います。

もとの行はi行、入れ替える行はk行なので

for(j = 0;j < 3;j++){
buf = a_matrix[i][j];
a_matrix[i][j] = a_matrix[k][j];
a_matrix[k][j] = buf;
}
buf = b_matrix[i];
b_matrix[i] = b_matrix[k];
b_matrix[k] = buf;

となります。

これらをすべてリストに直したものが以下にあります。
a_maxはbufで代用できるのでbufに変更、
for(i=0~)の最初に最大値の初期化としてbuf = 0;が入っています。

そして最大値buf（元a_max）は、マイナスであろうと0でなければいいので
絶対値をfabs()でとっています。fabs()はmath.hを#includeすれば使えます。

他にbuf(a_max)が0の時にERRORとして終了させると、
連立方程式が不成立の時に強制的に終了させることが出来ます。

以上、おつさまですた。

#include<stdio.h>
#include<math.h>

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{1,1,1},{1,1,3},{2,4,5}};
double b_matrix[3] = {6,8,19};
double buf;

PRINT(a_matrix,b_matrix);

//前進消去
for(i = 0;i < 3;i++){
buf = 0;
for(j = i;j < 3;j++){
if(buf < fabs(a_matrix[j][i])){
buf = a_matrix[j][i];
k = j;
}
}

for(j = 0;j < 3;j++){
buf = a_matrix[i][j];
a_matrix[i][j] = a_matrix[k][j];
a_matrix[k][j] = buf;
}
buf = b_matrix[i];
b_matrix[i] = b_matrix[k];
b_matrix[k] = buf;

buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

//後退代入
for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;

for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}