ガウスの消去法の前進消去を拡張すると正方行列の逆行列を求めることができる. これは,ガウス・ジョルダン法(Gauss-Jordan method),もしくは,掃出し法(sweeping-out method)と呼ばれる.

まず,ls_gauss_jordan.eq1.gifの正方行列ls_gauss_jordan.eq2.gifを考える. ls_gauss_jordan.eq3.gifならばls_gauss_jordan.eq2.gifは正則行列(regular matrix)である (逆にls_gauss_jordan.eq4.gifなら特異行列(singular matrix)). ls_gauss_jordan.eq1.gifの単位行列をls_gauss_jordan.eq5.gifとすると,ls_gauss_jordan.eq2.gifが正則ならば,ls_gauss_jordan.eq6.gifを満足する行列ls_gauss_jordan.eq7.gifが1つだけ存在する. このls_gauss_jordan.eq7.gifを逆行列と呼びls_gauss_jordan.eq8.gifで表される.

さて,ガウスの消去法で対象とした連立1次方程式ls_gauss_jordan.eq9.gifにおいて,ls_gauss_jordan.eq10.gifls_gauss_jordan.eq7.gifに,ls_gauss_jordan.eq11.gifls_gauss_jordan.eq5.gifに置き換えることを考える. そうすると拡大行列は以下のようになる.

ls_gauss_jordan.eq12.gif

ls_gauss_jordan.eq6.gifであるので,両辺にls_gauss_jordan.eq8.gifをかけると

ls_gauss_jordan.eq13.gif

となる. つまり,拡大行列の左側半分が単位行列となるように変換することで,右側半分に逆行列ができる.

ls_gauss_jordan.eq14.gif

このとき,逆行列ls_gauss_jordan.eq8.gifは,

ls_gauss_jordan.eq15.gif

となる. 前進消去を拡張し,これらの処理を行うことで逆行列を求めるのが, ガウス・ジョルダン法である.

ガウス・ジョルダン法の実装

ガウスの消去法における前進消去では上三角行列にするために, ls_gauss_jordan.eq16.gifの条件を満たす要素のみ対象として処理していたが, ガウス・ジョルダン法ではすべての要素に対象を広げる. ガウス・ジョルダンでの前身消去の式は以下である.

ls_gauss_jordan.eq17.gif

ガウス・ジョルダン法で逆行列を求めるコード例を以下に示す.

  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
 35
/*!
 * ガウス・ジョルダン法(ピボット選択なし)
 * @param[inout] A n×2nの拡張行列
 * @param[in] n n元連立一次方程式
 * @return 1:成功
 */
int GaussJordan(vector< vector<double> > &A, int n)
{
    // 拡張行列の右半分を単位行列にする
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            A[i][j+n] = (i == j ? 1.0 : 0.0);
        }
    }
 
    // ガウス・ジョルダン法(Gauss-Jordan method)で逆行列計算
    for(int k = 0; k < n; ++k){
        double akk = A[k][k];
        // 対角要素を1にするために,k行目のすべての要素をa_kkで割る
        for(int j = 0; j < 2*n; ++j){
            A[k][j] /= akk;
        }
 
        // k列目の非対角要素を0にする
        for(int i = 0; i < n; ++i){
            if(i == k) continue;
            double aik = A[i][k];
            for(int j = 0; j < 2*n; ++j){
                A[i][j] -= A[k][j]*aik;
            }
        }
    }
 
    return 1;
}

添付ファイル: filels_gauss_jordan.eq9.gif 619件 [詳細] filels_gauss_jordan.eq15.gif 591件 [詳細] filels_gauss_jordan.eq6.gif 670件 [詳細] filels_gauss_jordan.eq7.gif 571件 [詳細] filels_gauss_jordan.eq5.gif 588件 [詳細] filels_gauss_jordan.eq17.gif 573件 [詳細] filels_gauss_jordan.eq14.gif 576件 [詳細] filels_gauss_jordan.eq8.gif 546件 [詳細] filels_gauss_jordan.eq12.gif 610件 [詳細] filels_gauss_jordan.eq4.gif 574件 [詳細] filels_gauss_jordan.eq1.gif 544件 [詳細] filels_gauss_jordan.eq13.gif 555件 [詳細] filels_gauss_jordan.eq11.gif 539件 [詳細] filels_gauss_jordan.eq16.gif 558件 [詳細] filels_gauss_jordan.eq10.gif 596件 [詳細] filels_gauss_jordan.eq3.gif 583件 [詳細] filels_gauss_jordan.eq2.gif 382件 [詳細]

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