掃出法で連立一次方程式の解を求める

質問

掃出法で連立一次方程式の解を求めるプログラムを作ってみたのですが、ポインタと浮動小数点のエラーが出てしまい、実行できません。どこが間違っているのかさえ分からず困っています。訂正箇所を教えてください。宜しくお願い致します。

#include<stdio.h>
#include<math.h>
#include <float.h>
#define N 3
#define EPSILON 1.0E-5
#define TRUE 1
#define FALSE 0
void sweep(int *flag);
void swap(float *wk1,float *wk2);
float a[N][N]={{ 2, 6, 3},
{-1, 5,-2},
{-2,-1, 6}};
float x[N],b[N]={6,3,14};
int flag;
void main()
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%10.4f",a[i][j]);
printf("%10.4f\n",b[i]);
}
flag=TRUE;
sweep(&flag);

if(flag==TRUE)
{
printf("連立方程式の解\n");
for(i=0;i=N;i++)
printf("x[%d]=%10.4f\n",i+1,x[i]);
}
else
printf("解なし\n");
}

void swap(float *wk1,float *wk2)
{
float w;
w=*wk1;
*wk1=*wk2;
*wk2=w;
}

void sweep(int *flag)
{
int i,j,k,ik;
float ak,aik;
for(k=0;k<N;k++)
{
ak=a[k][k];
if(fabs(ak)<=EPSILON)
{
ik=k+1;
while((ik<N)&&(fabs(a[ik][k])<EPSILON))
ik++;
if(ik<N)
{
for(j=k;j<N;j++)
swap(&a[k][j],&a[ik][j]);
swap(b[k],b[ik]);
ak=a[k][k];
}

else
{
printf("ピボットが零です\n");
*flag=FALSE;
goto end;
}
}

for(j=k;j<N;j++)
a[k][j]=a[k][j]/ak;
b[k]=b[k]/ak;

for(i=0;i<N;i++)
{
if(i!=k)
{
aik=a[i][k];
for(j=k;j<N;j++)
a[k][j]=a[i][j]-aik*a[k][j];
b[i]=b[i]-aik*b[k];
}
}
for(k=0;k<N;k++)
x[k]=b[k];
end:;
}
}

回答No.3

asuncion

> w = *wk1;
> *wk2 = w;

この２行の間に、
*wk1 = *wk2;
の行が必要だと思います。

お礼

ご回答をくださった方々へ
ありがとうございました。
解決することができました。

回答No.2

mac_res

とりあえず動くようにしました。
-- 8< ---- 8< ---- 8< ---- 8< ---- 8< ---- 8< ---- 8< ---- 8< ---- 8< ---- 8< --
#include<stdio.h>
#include<math.h>
#include <float.h>
#define N 3
#define EPSILON 1.0E-5
#define TRUE 1
#define FALSE 0
void sweep(int *flag);
void swap(float *wk1, float *wk2);
float a[N][N] = { {2, 6, 3},
{-1, 5, -2},
{-2, -1, 6}
};
float x[N], b[N] = { 6, 3, 14 };
int flag;

int
main(void)
{
int i, j;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++)
printf("%10.4f", a[i][j]);
printf("%10.4f\n", b[i]);
}
flag = TRUE;
sweep(&flag);

if (flag == TRUE) {
printf("連立方程式の解\n");
for (i = 0; i < N; i++)
printf("x[%d]=%10.4f\n", i + 1, x[i]);
}
else
printf("解なし\n");

return 0;
}

void
swap(float *wk1, float *wk2)
{
float w;
w = *wk1;
*wk2 = w;
}

void
sweep(int *flag)
{
int i, j, k, ik;
float ak, aik;
for (k = 0; k < N; k++) {
ak = a[k][k];
if (fabs(ak) <= EPSILON) {
ik = k + 1;
while ((ik < N) && (fabs(a[ik][k]) < EPSILON))
ik++;
if (ik < N) {
for (j = k; j < N; j++)
swap(&a[k][j], &a[ik][j]);
swap(&b[k], &b[ik]);
ak = a[k][k];
}

else {
printf("ピボットが零です\n");
*flag = FALSE;
goto end;
}
}

for (j = k; j < N; j++)
a[k][j] = a[k][j] / ak;
b[k] = b[k] / ak;

for (i = 0; i < N; i++) {
if (i != k) {
aik = a[i][k];
for (j = k; j < N; j++)
a[k][j] = a[i][j] - aik * a[k][j];
b[i] = b[i] - aik * b[k];
}
}
for (k = 0; k < N; k++)
x[k] = b[k];
end:;
}
}

補足

回答してくださったのを参考にプログラムを書き換え、実行できたのですが、出力された解が間違っています。どのようにすれば正しい解にたどりつけるのでしょうか？

回答No.1

jacta

環境が書かれていないのでよく分かりませんが...

> どこが間違っているのかさえ分からず困っています。訂正箇所を教えてください。

コンパイルエラーになるので、そのときに行番号が示されると思います。それとも、そんなエラーレポートもできないようなコンパイラを使っているのでしょうか？