これまでn元連立1次方程式

ls_lu.eq1.gif

を直接法,反復法などで解くことだけを考えてきた. それでは実際の問題ではこのような線形システムはどのように扱われるだろうか. たとえば,制御の分野ではシステムの状態変化を捉えるために用いている. 係数行列Aがシステム内部を表し,右辺の定数ベクトルbで外乱などの状態を表す. システム内部が変わらず,外の状態が様々に変化したときの状態を知りたいとき, 係数行列Aは変化せず,bのみが変わるだけである. このように係数行列が固定で右辺の定数ベクトルだけが変化するということは, 物理学などの他の分野でも多くある. このとき,係数行列Aを解きやすい形に分解しておけば,計算量を大幅に減らすことができる. ここではそのような分解の一つであるLU分解について述べる.

n元連立1次方程式の係数行列を考える.

ls_lu.eq2.gif

これを以下の下三角行列(lower triangular matrix) L と 上三角行列 (upper triangular matrix) U に分解する.

ls_lu.eq3.gif
ls_lu.eq4.gif

ここで,ls_lu.eq5.gifである. このように分解することをLU分解と呼ぶ.

LU分解した行列で連立1次方程式を解くことを考える. まず,ls_lu.eq5.gifより,

ls_lu.eq6.gif

ls_lu.eq7.gifとすると,

ls_lu.eq8.gif

となる.Lは下三角行列であるので,1行目から順番に代入していくことで, 容易にls_lu.eq9.gifを算出することができる.この処理を前進代入(forward substitution)と呼ぶ. ls_lu.eq9.gifが求まったらそれを2番目の式に代入する. Uは上三角行列であるので,ガウスの消去法と同じく後退代入(backward substitution) していくことで解ls_lu.eq10.gifを求めることができる.

LU分解の手順

係数行列AをLU分解するためには,以下の式を1行目から順番に適用していく.

ls_lu.eq11.gif

メモリを節約するために,LU分解した結果を以下のように1つの行列に格納する.

ls_lu.eq12.gif

LU分解の実装

LU分解をC++で実装した例を以下に示す.

  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
/*!
 * LU分解(ピボット交換なし)
 *  - 行列A(n×n)を下三角行列(L:Lower triangular matrix)と上三角行列(U:Upper triangular matrix)に分解する
 *  - L: i >= j,  U: i < j の要素が非ゼロでUの対角成分は1
 *  - LとUを一つの行列にまとめた形で結果を返す
 * @param[inout] A n×nの係数行列.LU分解した結果を格納する.
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int LUDecomp(vector< vector<double> > &A, int n)
{
    if(n <= 0) return 0;
 
    for(int i = 0; i < n; ++i){
        // l_ijの計算(i >= j)
        for(int j = 0; j <= i; ++j){
            double lu = A[i][j];
            for(int k = 0; k < j; ++k){
                lu -= A[i][k]*A[k][j];    // l_ik * u_kj
            }
            A[i][j] = lu;
        }
 
        // u_ijの計算(i < j)
        for(int j = i+1; j < n; ++j){
            double lu = A[i][j];
            for(int k = 0; k < i; ++k){
                lu -= A[i][k]*A[k][j];    // l_ik * u_kj
            }
            A[i][j] = lu/A[i][i];
        }
    }
 
    return 1;
}

LU分解で連立1次方程式を解く手順

LU分解で連立1次方程式を解く方法は上でも述べたがここではより具体的な手順について説明する. LU分解で連立1次方程式を解く手順は以下である.

  1. 下三角行列Lと上三角行列Uを求める.
  2. 前進代入によりls_lu.eq9.gifを求める.
  3. 後退代入によりls_lu.eq10.gifを求める. 1はすでに述べているので2,3について説明する.
  • 前進代入 下三角行列Lと定数ベクトルbからyを求める.
    ls_lu.eq13.gif
    Lは下三角行列であるので,1行目は単純にls_lu.eq14.gifls_lu.eq15.gifが算出される. 得られたls_lu.eq15.gifを2行目の式に代入して同様に計算するとls_lu.eq16.gifが得られる.これをn行目まで繰り返す.式で書くと以下となる.
    ls_lu.eq17.gif
    ここで,ls_lu.eq18.gifである.
  • 後退代入 前進代入で得られたyと上三角行列Uからxを求める.
    ls_lu.eq19.gif
    基本的にはガウスの消去法の後退代入と同じである. ただし,各行は対角項で正規化されている(対角成分がすべて1). 前進代入とは逆にn行目から順番にls_lu.eq20.gifを求めていく.
    ls_lu.eq21.gif
    ここで,ls_lu.eq18.gifである.

LU分解を用いた連立1次方程式の解法の実装

LU分解で連立1次方程式を解くコードをC++で実装した例を以下に示す.

  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
/*!
 * LU分解した行列a(n×n)から前進代入・後退代入によりA・x=bを解く
 * @param[in] A LU分解された行列
 * @param[in] b 右辺ベクトル
 * @param[out] x 結果ベクトル
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int LUSolver(const vector< vector<double> > &A, const vector<double> &b, vector<double> &x, int n)
{
    if(n <= 0) return 0;
 
    // 前進代入(forward substitution)
    //  LY=bからYを計算
    for(int i = 0; i < n; ++i){
        double bly = b[i];
        for(int j = 0; j < i; ++j){
            bly -= A[i][j]*x[j];
        }
        x[i] = bly/A[i][i];
    }
 
    // 後退代入(back substitution)
    //  UX=YからXを計算
    for(int i = n-1; i >= 0; --i){
        double yux = x[i];
        for(int j = i+1; j < n; ++j){
            yux -= A[i][j]*x[j];
        }
        x[i] = yux;
    }
 
    return 1;
}

添付ファイル: filels_lu.eq18.gif 950件 [詳細] filels_lu.eq8.gif 887件 [詳細] filels_lu.eq4.gif 853件 [詳細] filels_lu.eq14.gif 863件 [詳細] filels_lu.eq10.gif 903件 [詳細] filels_lu.eq17.gif 799件 [詳細] filels_lu.eq12.gif 854件 [詳細] filels_lu.eq21.gif 814件 [詳細] filels_lu.eq3.gif 845件 [詳細] filels_lu.eq20.gif 859件 [詳細] filels_lu.eq11.gif 757件 [詳細] filels_lu.eq1.gif 888件 [詳細] filels_lu.eq15.gif 879件 [詳細] filels_lu.eq5.gif 908件 [詳細] filels_lu.eq6.gif 902件 [詳細] filels_lu.eq19.gif 734件 [詳細] filels_lu.eq13.gif 852件 [詳細] filels_lu.eq7.gif 832件 [詳細] filels_lu.eq16.gif 887件 [詳細] filels_lu.eq2.gif 915件 [詳細] filels_lu.eq9.gif 921件 [詳細]

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