ガウスの消去法のプログラムがどうしてもうまく動きません。

質問

こんにちは。あまりにも困ってしまったので質問させていただきました。
よろしければご回答をよろしくお願い致します。

さて、今ガウスの消去法のプログラムを作っているのですが、
どうしてもどうしても、正しい解を得ることができません。
今日はほぼ徹夜でずっとパソコン画面とにらめっこしていたのですが、
何をしてもどこをいじってもさっぱり上手く行かず、正直嫌気が差し始めているところです（泣）

Cの詳しい知識などは皆無に近い人間ですが、
こんなド素人を助けていただけませんか？
↓が僕の書いたソースコードです。間違いだらけで非常に見苦しいと思いますがお許し下さい。

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

int main(){
float a[3][3] = {
{-2, 4000000, -6000000},
{-2, 0.03, -0.2},
{1, -0.2, -0.05}
};
float b[3] = {5000000, 0.1, 0.1};
float max[3];
float maxp;
float temp;
float x[3];
float m;
int pivot = 0;
int i,j,k;

for( k = 0; k < 3; k++ ){

/* scaling */
for( j = 0; j < 3; j++ ){
max[j] = 0;
for( i = 0; i < 3; i++ ){
if( max[j] < fabs(a[j][i]) ) max[j] = fabs(a[j][i]);
}
}

for( j = 0; j < 3; j++ ){
for( i = 0; i < 3; i++ ){
a[j][i] = a[j][i] / max[j];
}
b[j] = b[j]/max[j];
}

/* pivoting */
for( j = k; j < 3; j++ ){
if ( maxp < fabs(a[j][k]) ){
pivot = j;
maxp = fabs(a[j][k]);
}
}

for( i = k; i < 3; i++ ){
temp = a[k][i];
a[k][i] = a[pivot][i];
a[pivot][i] = temp;
}

temp = b[k];
b[k] = b[pivot];
b[pivot] = temp;

/* forword elimination */
for( j = k+1; j < 3; j++ ){
m = a[j][k] / a[k][k];
for( i = k; i < 3; i++ ){
a[j][i] -= a[k][i] * m;
}
b[j] -= b[k] * m;
}

/* backward substitution */
for( j = 1; j >= 0; j-- ){
for( i = j+1; i > 3; i++){
m += a[j][i] * b[i];
x[j] -= m;
x[j] = x[j] / a[j][j];
}
}
}

for( j = 0; j < 3; j++ ){
for( i = 0; i < 3; i++ ){
printf("%f ", a[j][i]);
}
printf("\n");
printf("%f", x[j]);
}

return(0);
}

ご回答お待ちしております。
改めて、お見苦しいソースコードだったとは思いますが、ご容赦下さい。

回答No.4

php504 · ベストアンサー

自分がやった方法を書きますね
float だと誤差が大きくてだめだったのでdouble にしました
double a[3][3] = {
{-2, 4000000, -6000000},
{-2, 0.03, -0.2},
{1, -0.2, -0.05}
};
double b[3] = {5000000, 0.1, 0.1};
/* scaling */
の処理は削除します
/* pivoting */
の後に次の処理を入れて a[k][k] が1になるようにします。
m = a[k][k];
for( i = k; i < 3; i++ ){
a[k][i] = a[k][i] / m;
}
b[k] = b[k]/m;
/* backward substitution */
は for( k = 0; k < 3; k++ )ループの外に出して次のように変更しました
/* backward substitution */
for( j = 2; j >= 0; j-- ){
x[j] = b[j];
for( i = j+1; i < 3; i++){
x[j] -= a[j][i] * x[i];
}
}
結果は
0.037865
-0.087719
-0.891813
という解になりました。

お礼

おお、わざわざソースコードまで書いていただいてどうもありがとうございます。
細かい部分を参考にさせていただき、何とかプログラムは完成まで持って行くことができました。
回答ありがとうございました。

回答No.3

php504

/* scaling */
の処理は不要では
各行の最大値で割っていますが1行1列目が 2/6000000 で0になってます。

お礼

なるほど。

スケーリング処理をしないと絶対値が安定せず、誤差が大きくなってしまう。
という話をきいていたのですが、無くても大丈夫なのですか……。

ご回答ありがとうございました！

回答No.2

taihei759

解決策　
１、ガウスの消去法は有名です。
　　数値計算の本に、サンプルがあります。
　　それを参考にしてください。
　　スキャナーで読み取り、e-typeでデジタル化し、
　　貼付けるのも手です。
２、多分に、コードはチェックされていると思います。
　　｝あるいは｛　の数と位置が合わないのではないでしょうか。
3、見えない所に、全角の部分はないでしょうか。再度打ち直すのも
　手です。私も、全角の部分があり、原因不明で、打ち直した事があり　ます。

お礼

ご回答ありがとうございます。

インターネットを見ていても色々なところにソースコードが載っていますし、
本当に有名な問題のようですねー。
実際に色々なソースコードを見ながらところどころを参考にして、
何とかプログラムを完成させることができました。

ありがとうございました。

回答No.1

buriburi3

/* backward substitution */
for( j = 1; j >= 0; j-- ){
for( i = j+1; i > 3; i++){
アルゴリズムは見てないけどココのループ条件は明らかに変。
何をどーやっても一度も実行されないでしょ？
初回 j=1の時、i=2から始まるけど i>3 を満たさないから一度も実行されない。
結果として
m += a[j][i] * b[i];
x[j] -= m;
x[j] = x[j] / a[j][j];
の部位は一度も実行されない。
初期化もしていないのでx[ ]の値は不定のまま。
毎回変な値が出るけど間違った答えが出ているのではなく何も計算されていないメモリー上のゴミが表示されているだけ。

お礼

おおー、なるほど。

ループの条件の細かいことを合わせるのがどうにも苦手なのですが、
たしかに言われてみればその通りですね。非常に参考にさせていただきました。

質問させていただいたソースコードは全て一度消して、
皆様の意見を参考に最初から組み直したところ、無事に動くプログラムを作ることができました。

ご回答ありがとうございました。