C言語関係掲示板

過去ログ

No.923 Gauss-Jordan法

[戻る] [ホームページ]
No.11686

Gauss-Jordan法
投稿者---senshou(2004/01/12 02:54:29)


題名のアルゴリズムによるプログラムを書いてみた(つもり)なのですが、正しい答えが出てくれません。
どこが悪いのかご指摘下さい。

/* Gauss-Jordan法 */
#include <stdio.h>

int main(void)
{
    /* この辺は適当 */
    int n = 3;
    double a[3][4] = { {2, 5, 7, 23}, {4, 13, 20, 58}, {8, 29, 50, 132} };

    int i, j, k;
    double p;   //ピボット
    
    /* ここから Gauss-Jordan法 */
    for (k = 0; k < n; k++) {
        p = a[k][k];        //ピボットの設定
        
        /* k行目をピボットで割り、頭の係数を1にする。*/
        for (j = 0; j < n + 1; j++)
            a[k][j] /= p;
            
        /* 係数が1となっているk行目の式をa[i][k]倍しi行目から、引く。ただしi行目がk行目のときはだめ */
        for (i = 0; i < n; i++)
            if (i != k)
                for (j = 0; j < n + 1; j++)
                    a[i][j] -= a[i][k] * a[k][j];
    }

    /* 解の表示 */
    for (i = 0; i < n; i++)
        printf("x[%d]=%f\n", i, a[i][n]);
    
    return(0);
}


No.11688

Re:Gauss-Jordan法
投稿者---かずま(2004/01/12 10:08:08)


    for (j = 0; j < n + 1; j++)
        a[i][j] -= a[i][k] * a[k][j];
-->
    for (p = a[i][k], j = 0; j < n + 1; j++)
        a[i][j] -= p * a[k][j];


No.11692

Re:Gauss-Jordan法
投稿者---senshou(2004/01/12 15:44:19)


かずまさんありがとうございます。
うまくいきました。
a[i][k]を使っている間にループでa[i][k]の値が書き換わっていたのがよくなかったのですね。
だから別のものに入れて使うか、j = k + 1 とすれば(j <= k は使わないから)うまくいくと。