ガウスの消去法による3元連立1次方程式

[glafk]

コード:

#include<stdio.h>
#define N 3
int main(void)
{
    int h,i,j,k,x,y,moji;
    double a[N][N+1];
 
    printf("3元連立一次方程式を入力して下さい.\n");
    for(i=0;i<N;i++){
        moji='x';
        for(j=0;j<N+1;j++){
            if(j<N){
                printf("%d行目の式の%cの係数を入力:",i+1,moji);
                scanf("%lf",&a[i][j]);
                moji++;
            }
            else{
                printf("%d行目の式の右辺の数を入力:",i+1);
                scanf("%lf",&a[i][j]);
            }
        }
        putchar('\n');
    }
 
    printf("連立方程式\n");
    for(i=0;i<N;i++){
        moji='x';
        for(j=0;j<N+1;j++){
            if(j==N)
                printf(" = ");
            if(a[i][j]>0)
                printf("+");
            if(j<N){
                printf("%.0f%c",a[i][j],moji);
                moji++;
            }
            else
                printf("%.0f\n",a[i][j]);
        }
    }
    putchar('\n');
     
    //前進消去
    for(h=0;h<N;h++){
        x=a[h][h];
        for(i=0;i<N+1;i++)
            a[h][i]=a[h][i]/x;
        for(j=h+1;j<N;j++){
            y=a[j][h];
            for(k=0;k<N+1;k++)
                a[j][k]=a[j][k] - a[h][k] * y;
        }
    }
     //後退代入
    for(i=N-1;i>0;i--)
        for(j=i-1;j>=0;j--)
            a[j][N] = a[j][N] - a[j][i] * a[i][N];
             
    printf("解\n");
    moji='x';
    for(i=0;i<N;i++){
        printf("%c = %f\n",moji,a[i][N]);
        moji++;
    }
             
    return 0;   
}

見づらくて申し訳ありません、C言語を勉強しているのですが上手く解が出ません。
指摘お願い致します。

[みけCAT]

とりあえずピポット選択が無いですね。

コード:

1 2 0 4
5 6 7 8
9 10 11 12

を入力すると

コード:

3元連立一次方程式を入力して下さい.
1行目の式のxの係数を入力:1行目の式のyの係数を入力:1行目の式のzの係数を入力:1行目の式の右辺の数を入力:
2行目の式のxの係数を入力:2行目の式のyの係数を入力:2行目の式のzの係数を入力:2行目の式の右辺の数を入力:
3行目の式のxの係数を入力:3行目の式のyの係数を入力:3行目の式のzの係数を入力:3行目の式の右辺の数を入力:
連立方程式
+1x+2y0z = +4
+5x+6y+7z = +8
+9x+10y+11z = +12

解
x = -2.000000
y = 3.000000
z = -0.000000

と出力され、Octaveでの計算結果

コード:

octave:1> A = [1 2 0; 5 6 7; 9 10 11]
A =

    1    2    0
    5    6    7
    9   10   11

octave:2> B = [4; 8; 12]
B =

    4
    8
   12

octave:3> A^-1 * B
ans =

  -2.0000e+000
  3.0000e+000
  -1.7764e-015

octave:4>

と比較すると解は正しそうでした。
解が存在するが、このプログラムでは上手く解が出ないテストケース(入力)を教えていただけると助かります。
【編集】テストケースが見つかったので取り消し

[みけCAT]

連立方程式の表示部分は、

コード:

    printf("連立方程式\n");
    for(i=0;i<N;i++){
        moji='x';
        for(j=0;j<N+1;j++){
            if(j==N)
                printf(" = ");
            if(j>0 && j<N && a[i][j]>=0)
                printf("+");
            printf("%.0f%c",a[i][j],j<N?moji++:'\n');
        }
    }
    putchar('\n');

とするとより自然になると思います。

コード:

1 -2 0 4
5 6 -7 8
9 10 11 12

を入力すると、出力は

コード:

3元連立一次方程式を入力して下さい.
1行目の式のxの係数を入力:1行目の式のyの係数を入力:1行目の式のzの係数を入力:1行目の式の右辺の数を入力:
2行目の式のxの係数を入力:2行目の式のyの係数を入力:2行目の式のzの係数を入力:2行目の式の右辺の数を入力:
3行目の式のxの係数を入力:3行目の式のyの係数を入力:3行目の式のzの係数を入力:3行目の式の右辺の数を入力:
連立方程式
+1x-2y0z = +4
+5x+6y-7z = +8
+9x+10y+11z = +12

解
x = 2.385870
y = -0.807065
z = -0.130435

となり、Octaveでの計算結果

コード:

octave:1> [1 -2 0;5 6 -7;9 10 11]^-1 * [4; 8; 12]
ans =

   2.38710
  -0.80645
  -0.12903

octave:2>

と比較して大幅にずれることを確認しました。検証に入ります。

[みけCAT]

x,yはdouble型の値を格納することを意図した変数なのに、int型になっているため、
情報が欠落し、間違った解が出るようです。
x,yをdouble型にしたところ、出力が改善しました。

コード:

#include<stdio.h>
#define N 3
int main(void)
{
    int h,i,j,k,moji;
    double x,y,a[N][N+1];
 
    printf("3元連立一次方程式を入力して下さい.\n");
    for(i=0;i<N;i++){
        moji='x';
        for(j=0;j<N+1;j++){
            if(j<N){
                printf("%d行目の式の%cの係数を入力:",i+1,moji);
                scanf("%lf",&a[i][j]);
                moji++;
            }
            else{
                printf("%d行目の式の右辺の数を入力:",i+1);
                scanf("%lf",&a[i][j]);
            }
        }
        putchar('\n');
    }
 
    printf("連立方程式\n");
    for(i=0;i<N;i++){
        moji='x';
        for(j=0;j<N+1;j++){
            if(j==N)
                printf(" = ");
            if(a[i][j]>0)
                printf("+");
            if(j<N){
                printf("%.0f%c",a[i][j],moji);
                moji++;
            }
            else
                printf("%.0f\n",a[i][j]);
        }
    }
    putchar('\n');
     
    //前進消去
    for(h=0;h<N;h++){
#ifdef DO_PIVOT_SELECTION
        // ピポット選択
        double best=a[h][h]<0?-a[h][h]:a[h][h];
        int bestidx=h;
        for(j=h+1;j<N;j++){
            double cur=a[j][h]<0?-a[j][h]:a[j][h];
            if(cur>best){
                best=cur;
                bestidx=j;
            }
        }
        if(bestidx!=h){
            for(i=0;i<N+1;i++){
                double temp=a[h][i];
                a[h][i]=a[bestidx][i];
                a[bestidx][i]=temp;
            }
        }
#endif
        x=a[h][h];
        for(i=0;i<N+1;i++)
            a[h][i]=a[h][i]/x;
        for(j=h+1;j<N;j++){
            y=a[j][h];
            for(k=0;k<N+1;k++)
                a[j][k]=a[j][k] - a[h][k] * y;
        }
    }
     //後退代入
    for(i=N-1;i>0;i--)
        for(j=i-1;j>=0;j--)
            a[j][N] = a[j][N] - a[j][i] * a[i][N];
             
    printf("解\n");
    moji='x';
    for(i=0;i<N;i++){
        printf("%c = %f\n",moji,a[i][N]);
        moji++;
    }
             
    return 0;   
}

入力

コード:

1 -2 0 4
5 6 -7 8
9 10 11 12

出力

コード:

3元連立一次方程式を入力して下さい.
1行目の式のxの係数を入力:1行目の式のyの係数を入力:1行目の式のzの係数を入力:1行目の式の右辺の数を入力:
2行目の式のxの係数を入力:2行目の式のyの係数を入力:2行目の式のzの係数を入力:2行目の式の右辺の数を入力:
3行目の式のxの係数を入力:3行目の式のyの係数を入力:3行目の式のzの係数を入力:3行目の式の右辺の数を入力:
連立方程式
+1x-2y0z = +4
+5x+6y-7z = +8
+9x+10y+11z = +12

解
x = 2.387097
y = -0.806452
z = -0.129032

[glafk]

レスありがとうございます。
単純なところを見落としていました…ifの中身も含めご指摘ありがとうございました。