課題・ガウスの消去法2（後退代入）

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

前回に引き続き、ガウスの消去法。
前進消去が終わって上三角行列になったので、ここからは答えを求める後退代入。
用語はあっているか知りませんけど。

現時点でここまで求めることが出来ています。


ここから3行目を利用して2行目と1行目を消していきます。

1. 2行目3列が0になるように3行目を定数倍して2行目から引く
   （2行目マイナス、3行目3列かける2行目）
2. 1行目3列が0になるように3行目を定数倍して2行目から引く
   （1行目マイナス、3行目3列かける1行目）
3. 1行目2列が0になるように2行目を定数倍して1行目から引く
   （1行目マイナス、2行目2列かける1行目）

上の式が、
1.

2行目3列が0になって、

2.

1行目3列が0になって、

3.

1行目2列が0になり、答えが求まります。


前進消去と似たような感じなので手早くプログラムを羅列していきます。

/*1*/
buf = a_matrix[1][2];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[2][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[2][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[2][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[2]* buf;

/*2*/
buf = a_matrix[0][2];
a_matrix[0][0] = a_matrix[0][0] - a_matrix[2][0] * buf;
a_matrix[0][1] = a_matrix[0][1] - a_matrix[2][1] * buf;
a_matrix[0][2] = a_matrix[0][2] - a_matrix[2][2] * buf;
b_matrix[0] = b_matrix[0] - b_matrix[2] * buf;

/*3*/
buf = a_matrix[0][1];
a_matrix[0][0] = a_matrix[0][0] - a_matrix[1][0] * buf;
a_matrix[0][1] = a_matrix[0][1] - a_matrix[1][1] * buf;
a_matrix[0][2] = a_matrix[0][2] - a_matrix[1][2] * buf;
b_matrix[0] = b_matrix[0] - b_matrix[1] * buf;

これを前までのプログラムの下に入れてあげると（すべてのforループの外）、

x1 = 3,
x2 = 2,
x3 = 1

という答えになります。

ここから、今まで書いたプログラムの中に上手く入るようにしていきます。
このくらい短いプログラムの省略だとまず恒常的に、一番回るのが遅いループを探します。

この場合、buf = a_matrix[1][2]の右側の配列番号が、手順3まで変化しません。
反対に左側の部分は、手順2で0に変化しています。
この一番外のループをiにします。

次に、2番目に遅く回るループを探します。
さっきの左側の部分は手順が次に行くまで変化しません。これをjとします。

残りは手順内で0~2と変化するのでこれをkとします。

するとこのプログラムは、

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;
}

となります。

前回を含めたプログラムリストに直すと、

#include

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

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

PRINT(a_matrix,b_matrix);

/*前進消去*/
for(i = 0;i < 3;i++){
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;
}
}

PRINT(a_matrix,b_matrix);

/*後退代入*/
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;
}

これで連立一次方程式のガウス消去法による解が求まりました。

しかしこれだけでは

のような連立一次方程式の解は得られません。
次はピボットについて書く予定です。

おつさまですた。