掃き出し法による連立方程式の解法

[kuro-neko]

行列の掃き出し法で連立方程式を求めようとしていますが、
単位行列までは作れるのですが、解に相当するｘ[0][4]~x[3][4]
の値が入力値から変化していません。
プログラムを何度読み直しても、何処が違うのか分からず
困っています。

単純なミスの見落としだとは思うのですが、ご教授お願いします。

コード:

#include<stdio.h>

/* 掃き出し法により連立方程式を解く
w+4x-2y+4z=12
2w+3x+6y-2z=9
3w+7x+3y+2z=18
4w+x+7y-5z=3 */

int main(void){
	
	double x[4][5]={{1,4,-2,4,12},
                        {2,3,6,-2,9},
                        {3,7,3,2,18},
                        {4,1,7,-5,3}};
	
	int i,j,k;

/*三角行列を作る*/
	for(j=0;j<=2;j++){
		for(i=j;i<=2;i++){
			for(k=0;k<=4;k++){
			x[i+1][k]=x[i+1][k]
                        -(x[i+1][j]*x[j][k]/x[j][j]);		
			}
		}
	}

/*単位行列を作る*/
	
	for(k=0;k<=3;k++){	
		for(j=k;j<=2;j++){
			for(i=0;i<=4;i++){
				x[k][i]=x[k][i]
                        	-(x[j+1][i]*x[k][j+1]/x[j+1][j+1]);
			}
		}	
	}

	for(j=0;j<=3;j++){
		for(i=0;i<=4;i++){
			x[j][i]=x[j][i]/x[j][j];
		}
	}
/*画面出力*/

	for(i=0;i<=3;i++){
		printf("%3.0f,%3.0f,%3.0f,%3.0f,%3.0f\n"
                       ,x[i][0],x[i][1],x[i][2],x[i][3],x[i][4]);
	}

	printf("W=%3.0f\n",x[0][4]);
	printf("x=%3.0f\n",x[1][4]);
	printf("y=%3.0f\n",x[2][4]);
	printf("z=%3.0f\n",x[3][4]);
	printf("\n");
	return 0;
}

[salsaww]

そもそも、三角行列から作れていないで、まずはそこから作りましょう。



素直に、http://www.asahi-net.or.jp/~uc3k-ymd/Le ... epOut.html
などを見て、書き直した方がいいかと思います。

ijkなどの変数が何をしてるのか？ちょっと分り難いコードなのと、
微妙に数学の恒等式と、プログラミングにおける代入式を混同しているか、
ijkの変化によって、実際に対象がどの様に変化するのかが把握出来ていない様な気がします。

[kuro-neko]

計算されない理由は掃き出し法で０になった数値を使って計算していたからでした。
そのため各行の係数となる数字は別の変数に代入してループの外に出しました。
このことで正しい値が計算されるようになりました。

返信頂いた方どうもありがとうございました。
ijkの使い方が統一性なく解り辛かったですね。
ご教示いただいたHPも参考になりました。ありがとうございます。

以下のコードで正しく計算されました。

コード:

#include<stdio.h>
 
/* 掃き出し法により連立方程式を解く
w+4x-2y+4z=12
2w+3x+6y-2z=9
3w+7x+3y+2z=18
4w+x+7y-5z=3 */
 
int main(void){
    
    double x[4][5]={{1,4,-2,4,12},
                        {2,3,6,-2,9},
                        {3,7,3,2,18},
                        {4,1,7,-5,3}};
	double su1=0,su2=0;
    
    int i,j,k;
	
 
/*三角行列を作る*/
	for(k=0;k<=2;k++){
		for(j=k+1;j<=3;j++){
			su1=x[k][k];
			su2=x[j][k];
			for(i=0;i<=4;i++){
				x[j][i]=x[j][i]-x[k][i]*su2/su1;
			}
		}
	}
	
/*単位行列を作る*/
	for(k=0;k<=2;k++){
		for(j=k+1;j<=3;j++){
			su1=x[3-k][3-k];
			su2=x[3-j][3-k];
			for(i=0;i<=4;i++){
				x[3-j][i]=x[3-j][i]-x[3-k][i]*su2/su1;
			}
		}
	}


	su1=0;
    	for(j=0;j<=3;j++){
		su1=x[j][j];
        	for(i=0;i<=4;i++){
            		x[j][i]=x[j][i]/su1;
        	}
    	}

/*画面出力*/

printf("**********************計算結果***********************\n");
 
    for(i=0;i<=3;i++){
        printf("%3.0f,%3.0f,%3.0f,%3.0f,%3.0f\n"
                       ,x[i][0],x[i][1],x[i][2],x[i][3],x[i][4]);
    }
 
    printf("W=%3.0f\n",x[0][4]);
    printf("x=%3.0f\n",x[1][4]);
    printf("y=%3.0f\n",x[2][4]);
    printf("z=%3.0f\n",x[3][4]);
    printf("¥n");
    return 0;
}