ガウスの消去法のプログラム

質問

ガウスの消去法（部分ピボット選択）のプログラムを組んでみたつもりなのですが上手くいきません。

間違いだらけだと思いますがどうかアドバイスをして頂けませんでしょうか？

#include <stdio.h>

double main(void){
int i,j,k,N,M,m;
double A[N+1][N+1][N+1],B[j],S,X[N+1];
printf("次数の入力。\n");
scanf("%d",&N);

for(k=1;k<=N;k++){
for(i=1;i<=N;i++){
for(j=1;j<=N+1;j++){
A[k][i][j]=0;
}
}
}

for(i=1;i<=N;i++){
for(j=1;j<=N+1;j++){
printf("係数の入力.\n A[1][%d][%d]?\n",i,j);
scanf("%f",&A[1][i][j]);
}
}

for(k=2;k<=N;k++){
if(A[k-1][k-1][k-1]=0){
for(M=k;M<=N;M++){
if(A[k-1][M][k-1]!=0){
for(j=k-1;j<=N+1;j++){
B[j]=A[k-1][k-1][j];
A[k-1][k-1][j]=A[k-1][M][j];
A[k-1][M][j]=B[j];
}
goto abc;
}
else {printf("解は無い\n");}
}
}

abc:
for(i=k;i<=N;i++){
for(j=k;j<=N+1;j++){
A[k][i][j]=A[k-1][i][j]-(A[k-1][i][k-1]/A[k-1][k-1][k-1])*A[k-1][k-1][j];
}
}
X[N]=A[N][N][N+1]/A[N][N][N];
printf("解Ｘ(Ｎ)は %f 。\n",X[N]);

for(k=N-1;k>=1;k--){
S=0;
for(m=N;m>=k+1;m--){
S+=A[k][k][m]*X[m];
}
X[k]=(A[k][k][N+1]-S)/A[k][k][k];
printf("解Ｘ(%d)は %f 。\n",k,X[k]);}
}
}

回答No.1

kary · ベストアンサー

ガウス消去法には3元配列は必要なく(Cで3元配列は使えないと思うのですが)2元配列でできます。概要は、例えば3元の連立方程式を解く場合、配列はX[3][4]を確保し、未知数の係数3×3を0～2列に、既知数の定数を3列に代入します。ピボット選択法ということですが、説明を簡単にするために、最大値の選択をしないとすると、

1.　0行目の全要素X[0][0～3]をX[0][0]で割る。
2.　1行目に対してX[1][i]-=X[0][i]*X[1][0](i=0～3)を行う。2行目も同様にする。これにより、0列目はX[0][0]=1.0、X[1～2][0]=0.0となる。
これを1ステップとして、次のステップは、
1.　1行目の全要素X[1][1～2]をX[1][1]で割る。(X[1][0]は0.0なので割らなくても良い)
...
これを2行目まで繰り返すと、X[0～2][0～2]は対角要素のみが1.0で、その他は0.0の行列となり、X[0～2][3]が答えとなる。
このようにすると、X[3][4]の配列のみで計算ができます(ループカウンタ等は必要ですが)。ピボット選択法の場合は、1.の前に行の入れ替えをします。
ガウス消去法についての知識がおありな様でしたので、とりあえず概要のみの回答とさせて頂きました。(こんなことぐらい知ってるよ、という内容で参考にならないかも知れませんが...)

お礼

こんなに早く回答してくださって有難うございました。

>こんなことぐらい知ってるよ、という内容で参考にならな>いかも知れませんが・・・

知らない事ばかりです、とても参考になりました。
Ｃ言語に関する知識なんかは本当に皆無に近い状態でしたのでご指摘を受けなければ、解らないままになっていたと思います