3元連立方程式をガウスの消去法で解いてみた

タイトルまんまです（ガウスの消去法についてはwiki参照のこと）。逆行列を求めるのに似たような方法を使うのですが、微妙に忘れかけてたので復習を兼ねて解いてみた。本当は精度の問題とかもあって、もっとちゃんとしたアルゴリズムがあるらしいが、そんなのは知らん。あまりにひどいハードコーディングなので、明日はN元まで拡張するつもり。

#include <iostream>
#include <iomanip>

void PrintArray(const double *pArray, int col, int row)
{
    for(int i=0; i<col; ++i)
    {
        for(int j=0; j<row; ++j)
        {
            std::cout << std::setw(10) << pArray[i*row+j] << " ";
        }
        std::cout << std::endl;
    }
    std::cout << "------------------------------" << std::endl;

}

// 入力用データ
// x,y,z = 2, 3, 4
// 2x + 3y + 2z = 21 (1)
// 3x + 4y + 3z = 30 (2)
// 5x + 2y + 4z = 32 (3)
double mat[4*4] = 
{
    2, 3, 2, 21,
    3, 4, 3, 30,
    5, 2, 4, 32,
    1, 1, 1, 1,
};

int main(int argc, char **argv)
{
    //計算用配列
    double tmp[4*4] = {1};

    //まずは(1)式のxの係数を1にする
    for(int i=0; i<4; ++i)
    {
        tmp[i] = mat[i] / mat[0];
    }
    PrintArray(tmp, 4, 4);

    //(2),(3)式のxの係数を１に合わせる
    for(int i=0; i<4; ++i)
    {
        tmp[1*4+i] = mat[1*4+i] / mat[1*4];
        tmp[2*4+i] = mat[2*4+i] / mat[2*4];
    }
    PrintArray(tmp, 4, 4);

    //(1)式から(2),(3)式を引いたものをtmp[2,3][0-3]に格納
    for(int i=0; i<4; ++i)
    {
        tmp[1*4+i] = tmp[0*4+i] - tmp[1*4+i];
        tmp[2*4+i] = tmp[0*4+i] - tmp[2*4+i];
    }
    PrintArray(tmp, 4, 4);

    //次に(2)式のyの係数と(3)式のyの係数を合わせる
    double yk = tmp[1*4+1] / tmp[2*4+1];
    for(int i=1; i<4; ++i)
    {
        tmp[2*4+i] = tmp[2*4+i] * yk;
    }
    PrintArray(tmp, 4, 4);

    //(2)式から(3)式を引く
    for(int i=1; i<4; ++i)
    {
        tmp[2*4+i] = tmp[1*4+i] - tmp[2*4+i];
    }
    PrintArray(tmp, 4, 4);

    //z成分を求める
    double z = tmp[3*4+2] = tmp[2*4+3] / tmp[2*4+2];
    PrintArray(tmp, 4, 4);
    
    //y成分を求める
    double y = tmp[3*4+1] = (tmp[1*4+3] - (z*tmp[1*4+2])) / tmp[1*4+1];
    PrintArray(tmp, 4, 4);

    //xを求める
    double x = tmp[3*4+0] = (tmp[0*4+3] - (z*tmp[0*4+2] + y*tmp[0*4+1])) / tmp[0*4+0];
    PrintArray(tmp, 4, 4);

    std::cout <<"(x,y,z) = (" << x << ", " << y << ", " << z << ")" << std::endl;  

    return 0;
} 

どうでもいいがなんで2次元配列使わなかったんだろう俺…。あと所々でPrintArrayを呼び出してるのはデバッグのためです。あしからず。