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

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

#ref(ls_gauss_elimination.eq1.gif,nolink,70%)


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

#ref(ls_gauss_elimination.eq2.gif,nolink,70%)


あとは，&ref(ls_gauss_elimination.eq3.gif,nolink,70%);の結果を1番目の式に代入することで，&ref(ls_gauss_elimination.eq4.gif,nolink,70%);を得る．

これを行列で表すと，

#ref(ls_gauss_elimination.eq5.gif,nolink,70%)


に対して，2行目に&ref(ls_gauss_elimination.eq6.gif,nolink,70%);をかけて引くことで，

#ref(ls_gauss_elimination.eq7.gif,nolink,70%)


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

#ref(ls_gauss_elimination.eq8.gif,nolink,70%)


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

#ref(ls_gauss_elimination.eq10.gif,nolink,70%)


この処理のことを前進消去(forward elimination)と呼ぶ．
前進消去により，未知数&ref(ls_gauss_elimination.eq11.gif,nolink,70%);は&ref(ls_gauss_elimination.eq12.gif,nolink,70%);で求められ，
この&ref(ls_gauss_elimination.eq11.gif,nolink,70%);を&ref(ls_gauss_elimination.eq13.gif,nolink,70%);行目に代入することで&ref(ls_gauss_elimination.eq14.gif,nolink,70%);が得られ，
&ref(ls_gauss_elimination.eq15.gif,nolink,70%);を&ref(ls_gauss_elimination.eq16.gif,nolink,70%);行目に代入することで&ref(ls_gauss_elimination.eq17.gif,nolink,70%);が得られる．
これを&ref(ls_gauss_elimination.eq18.gif,nolink,70%);まで繰り返すことですべての解を得る．
この後半部分の処理のことを後退代入(back substitution)と呼ぶ．
前進消去と後退代入により連立1次方程式の解を求める方法をガウスの消去法という．

-ガウスの消去法の実装~

***ガウスの消去法の実装 [#e4c1bbf6]
ガウスの消去法を&ref(ls_gauss_elimination.eq19.gif,nolink,70%);元連立1次方程式に適用する．
まず，&ref(ls_gauss_elimination.eq9.gif,nolink,70%);に定数ベクトル&ref(ls_gauss_elimination.eq20.gif,nolink,70%);を列の最後に追加した&ref(ls_gauss_elimination.eq21.gif,nolink,70%);の行列を考える．
#ref(ls_gauss_elimination.eq22.gif,nolink,70%)

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

拡大行列の要素を
#ref(ls_gauss_elimination.eq23.gif,nolink,70%)

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

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

#ref(ls_gauss_elimination.eq24.gif,nolink,70%)


ここで，

#ref(ls_gauss_elimination.eq25.gif,nolink,70%)


上式により&ref(ls_gauss_elimination.eq26.gif,nolink,70%);を&ref(ls_gauss_elimination.eq27.gif,nolink,70%);で更新することで上三角行列を得る．
この処理のC++のコードを以下に示す．
#code(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]には係数行列と定数ベクトルが格納されている．

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

#ref(ls_gauss_elimination.eq28.gif,nolink,70%)


後退代入のC++のコード例を以下に示す．
#code(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]に格納しているので注意．

ガウスの消去法の全体のコードは以下．
#code(C){{
/*!
 * ガウスの消去法(ピボット交換なし)
 * @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を使っている．
トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS