ガウスの消去法

ガウスの消去法(Gaussian elimination)と呼ばれる方法では， 前進消去と後退代入という2段階の計算により解を求める．

一般的な(中学校で習う)解き方の一つに加減法というものがある． たとえば，以下のような連立方程式で考える．

加減法では片方の式の両辺に実数をかけて，係数を合わせて，加減算する． ここでは下の式の両辺を2倍して，上の式から引くことにする．

あとは，の結果を1番目の式に代入することで，を得る．

これを行列で表すと，

に対して，2行目にをかけて引くことで，

という形にする．そして，一番下の行から未知数を解いて代入していく．

前半の処理では係数行列の対角要素を除いた左下部分がすべて0となる 以下に示すような上三角行列(upper triangle matrix)を求めている．

この処理のことを前進消去(forward elimination)と呼ぶ． 前進消去により，未知数はで求められ， このを行目に代入することでが得られ， を行目に代入することでが得られる． これをまで繰り返すことですべての解を得る． この後半部分の処理のことを後退代入(back substitution)と呼ぶ． 前進消去と後退代入により連立1次方程式の解を求める方法をガウスの消去法という．

ガウスの消去法の実装†

ガウスの消去法を元連立1次方程式に適用する． まず，に定数ベクトルを列の最後に追加したの行列を考える．

この行列は拡大行列と呼ばれる．

拡大行列の要素を

とする．ここで，C,C++で実装することを考えて要素番号を0始まりにしている． 今後，実装の説明では0で始まるインデックスを用いるので注意．

前進消去は次の式で表される．

ここで，

上式によりをで更新することで上三角行列を得る． この処理のC++のコードを以下に示す．

// 前進消去(forward elimination)
    //  - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする)
    for(int k = 0; k < n-1; ++k){
        double akk = A[k][k];
        for(int i = k+1; i < n; ++i){
            double aik = A[i][k];
            for(int j = k; j < n+1; ++j){ // 確認のため左下要素が0になるようにj=kとしたが，実際にはj=k+1でよい
                A[i][j] = A[i][j]-aik*(A[k][j]/akk);
            }
        }
    }

2次元配列A[n][n+1]には係数行列と定数ベクトルが格納されている．

次に後退代入の式を以下に示す．

後退代入のC++のコード例を以下に示す．

// 後退代入(back substitution)
    //  - x_nの解はb_n/a_nn，x_nをさらにn-1行の式に代入することでx_(n-1)を求める．
    //  - この作業を1行目まで続けることですべての解を得る．
    A[n-1][n] = A[n-1][n]/A[n-1][n-1];
    for(int i = n-2; i >= 0; --i){
        double ax = 0.0;
        for(int j = i+1; j < n; ++j){
            ax += A[i][j]*A[j][n];
        }
        A[i][n] = (A[i][n]-ax)/A[i][i];
    }

計算された結果は別の配列でなくA[i][n]に格納しているので注意．

ガウスの消去法の全体のコードは以下．

/*!
 * ガウスの消去法(ピボット交換なし)
 * @param[inout] A n×nの係数項とn×1の定数項(b)を併せたn×(n+1)の行列．n+1列目に解が入る．
 * @param[in] n n元連立一次方程式
 * @return 1:成功
 */
int GaussElimination(vector< vector<double> > &A, int n)
{
    // 前進消去(forward elimination)
    //  - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする)
    for(int k = 0; k < n-1; ++k){
        double akk = A[k][k];
        for(int i = k+1; i < n; ++i){
            double aik = A[i][k];
            for(int j = k; j < n+1; ++j){ // 確認のため左下要素が0になるようにj=kとしたが，実際にはj=k+1でよい
                A[i][j] = A[i][j]-aik*(A[k][j]/akk);
            }
        }
    }

    // 後退代入(back substitution)
    //  - x_nの解はb_n/a_nn，x_nをさらにn-1行の式に代入することでx_(n-1)を求める．
    //  - この作業を1行目まで続けることですべての解を得る．
    A[n-1][n] = A[n-1][n]/A[n-1][n-1];
    for(int i = n-2; i >= 0; --i){
        double ax = 0.0;
        for(int j = i+1; j < n; ++j){
            ax += A[i][j]*A[j][n];
        }
        A[i][n] = (A[i][n]-ax)/A[i][i];
    }

    return 1;
}

標準の2次元配列の代わりにSTLのvectorを使っている．