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

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

ls_gauss_elimination.eq1.gif

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

ls_gauss_elimination.eq2.gif

あとは,ls_gauss_elimination.eq3.gifの結果を1番目の式に代入することで,ls_gauss_elimination.eq4.gifを得る.

これを行列で表すと,

ls_gauss_elimination.eq5.gif

に対して,2行目にls_gauss_elimination.eq6.gifをかけて引くことで,

ls_gauss_elimination.eq7.gif

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

ls_gauss_elimination.eq8.gif

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

ls_gauss_elimination.eq10.gif

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

ガウスの消去法の実装

ガウスの消去法をls_gauss_elimination.eq19.gif元連立1次方程式に適用する. まず,ls_gauss_elimination.eq9.gifに定数ベクトルls_gauss_elimination.eq20.gifを列の最後に追加したls_gauss_elimination.eq21.gifの行列を考える.

ls_gauss_elimination.eq22.gif

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

拡大行列の要素を

ls_gauss_elimination.eq23.gif

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

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

ls_gauss_elimination.eq24.gif

ここで,

ls_gauss_elimination.eq25.gif

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
    // 前進消去(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]には係数行列と定数ベクトルが格納されている.

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

ls_gauss_elimination.eq28.gif

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
    // 後退代入(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]に格納しているので注意.

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
/*!
 * ガウスの消去法(ピボット交換なし)
 * @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を使っている.


添付ファイル: filels_gauss_elimination.eq15.gif 592件 [詳細] filels_gauss_elimination.eq20.gif 573件 [詳細] filels_gauss_elimination.eq19.gif 641件 [詳細] filels_gauss_elimination.eq9.gif 589件 [詳細] filels_gauss_elimination.eq21.gif 589件 [詳細] filels_gauss_elimination.eq11.gif 525件 [詳細] filels_gauss_elimination.eq12.gif 512件 [詳細] filels_gauss_elimination.eq22.gif 507件 [詳細] filels_gauss_elimination.eq14.gif 537件 [詳細] filels_gauss_elimination.eq24.gif 575件 [詳細] filels_gauss_elimination.eq3.gif 537件 [詳細] filels_gauss_elimination.eq13.gif 557件 [詳細] filels_gauss_elimination.eq4.gif 550件 [詳細] filels_gauss_elimination.eq5.gif 575件 [詳細] filels_gauss_elimination.eq18.gif 537件 [詳細] filels_gauss_elimination.eq25.gif 510件 [詳細] filels_gauss_elimination.eq23.gif 552件 [詳細] filels_gauss_elimination.eq2.gif 530件 [詳細] filels_gauss_elimination.eq8.gif 467件 [詳細] filels_gauss_elimination.eq6.gif 529件 [詳細] filels_gauss_elimination.eq1.gif 529件 [詳細] filels_gauss_elimination.eq10.gif 570件 [詳細] filels_gauss_elimination.eq27.gif 517件 [詳細] filels_gauss_elimination.eq26.gif 487件 [詳細] filels_gauss_elimination.eq17.gif 486件 [詳細] filels_gauss_elimination.eq7.gif 480件 [詳細] filels_gauss_elimination.eq16.gif 528件 [詳細] filels_gauss_elimination.eq28.gif 555件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-06-26 (火) 19:27:58 (3318d)