ガウスの消去法、後退代入について

解決済

質問者：neoon
質問日時：2009/06/22 11:12
回答数：1件

ガウスの消去法で、n変数n式の連立一次方程式を解くプログラムを
作っています。
ですが、前進消去およびピボット選択ははおそらくできたのですが
後退代入が悪いのか、答えが出ません。
もしかしたら前進消去部分も悪いのかもしれませんが；；

いろいろ調べたのですが
どうしたら直るのかわからなくなってしまいました。
どのようにして直せばいいのか教えてください。

プログラムは以下です。

#include<stdio.h>
#include<math.h>
#define SIZE 100/*上限MAXのマクロ定義*/
#define EPS (1.0e-5)

int main()
{
/*変数宣言*/
int n; /*連立方程式の式数*/
int i; /*1つ目のループカウンタ*/
int j; /*2つ目のループカウンタ*/
int k;/*3つ目のループカウンタ*/
int pivot; /*pivot選択を行う際に使用する変数*/
double max; /*pivot選択用の絶対値が最大の値*/
double matrix_a[SIZE][SIZE]; /*変数xijの係数aij*/
double matrix_b[SIZE]; /*式iの定数項bi*/
double a_hozon;/*matrix_aの値を一時保存する変数*/
double b_hozon;/*matrix_bの値を一時保存する変数*/
double r_hozon;/*前進消去法をする時使用する値を一時保存する変数*/
int l;
int m;

/*入力処理*/
/*未知数nの入力*/
printf("n変数n式の連立一次方程式を解きます。\n");
printf("式数n(1～100)を入力してください。 n=?\n");
scanf("%d",&n);
/*未知数nのチェック*/
if(n>SIZE||n<1)
{
printf("\n\nn=%dはプログラムの利用可能範囲外です。\n",n);
return -1;
}

printf("%d変数%d式の一次方程式を入力して下さい。\n",n,n);
printf("入力は、式毎に係数、定数項の順に行ってください。\n\n");

printf("連立一次方程式の拡大係数行列を表示します。\n");
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
scanf("%lf",&matrix_a[i][j]);

}
scanf("%lf",&matrix_b[i]);

}
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
printf("%6.2f",matrix_a[i][j]);
}
printf("|%6.2f\n",matrix_b[i]);
}

printf("計算します。\n");

/*pivot選択*/
for(k=0; k<n; k++)
{
max=0.0;
pivot=0;
for(i=k; i<n; i++)
{
if(fabs(matrix_a[i][k])>max)
{
max=fabs(matrix_a[i][k]);
pivot=i;
}
}
if(matrix_a[pivot][k]==0.0)
{
printf("正則性なし\n");
}
if(pivot!=k)
{
for(j=0; j<n; j++)
{
a_hozon = matrix_a[k][j];
matrix_a[k][j] = matrix_a[pivot][j];
matrix_a[pivot][j] = a_hozon;
}
b_hozon = matrix_b[k];
matrix_b[k] = matrix_b[pivot];
matrix_b[pivot] = b_hozon;
}

}

/*入れ換え後の行列を表示*/
for(l=0; l<n; l++)
{
for(m=0; m<n; m++)
{
printf("%6.2f",matrix_a[l][m]);
}
printf("|%6.2f\n",matrix_b[l]);
}

/*消去法*/

for(k=0; k<n; k++)
{
for(i=k+1; i<n-1; i++)
{
r_hozon = matrix_a[i][k]/matrix_a[k][k];
for(j=k; j<n; j++)
{
matrix_a[i][j] = matrix_a[i][j]-r_hozon*matrix_a[k][j];
}
matrix_b[i] = matrix_b[i]-r_hozon*matrix_b[k];
}
}

/*後退代入*/

for(i=n-1; i>=0; i--)
{
matrix_b[i]=0.0;
for(j=i+1; j<n; j++)
{
matrix_b[i]=matrix_b[i]-matrix_a[i][j]*matrix_b[j];
}
}
/*
for(i=n-1; i>0; i--)
{
temp = 0.0;
for(j=i+1; j<n; j++)
{

xn = b[n]/a[n][n];
xn-1 = (-a[n-1][n]*x[n]+b[n-1])/a[n-1][n-1];
xn-2 = (-a[n-2][n-1]*x[n-2]+a[n-2][n]*x[n]+b[n-2])/a[n-2][n-2];

xn-j = (-a[i][n-j]*x[i+1]+b[i])/a[i][i];

temp = matrix_a[][]

temp=temp+matrix_a[i][j]*matrix_a[j][n];
matrix_a[i][n]=matrix_a[i][n]-temp;
matrix_a[i][n]=matrix_a[i][n]/matrix_a[i][i];
}

}
*/

/*不定、不能の判別*/
if(matrix_a[k][k]<EPS)
{
printf("解は一意に求まらない。\n");
return -1;
}

/*出力処理*/
printf("連立一次方程式の解を表示します。\n");
for(i=0; i<n; i++)
{
printf("x%d = %6.2f\n",i,matrix_b[i]);
}

return 0;
}