これまでn元連立1次方程式
#ref(ls_lu.eq1.gif,nolink,70%)

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

n元連立1次方程式の係数行列を考える．
#ref(ls_lu.eq2.gif,nolink,70%)

これを以下の下三角行列(lower triangular matrix) L と 上三角行列 (upper triangular matrix) U に分解する．
#ref(ls_lu.eq3.gif,nolink,70%)

#ref(ls_lu.eq4.gif,nolink,70%)

ここで，&ref(ls_lu.eq5.gif,nolink,70%);である．
このように分解することをLU分解と呼ぶ．

LU分解した行列で連立1次方程式を解くことを考える．
まず，&ref(ls_lu.eq5.gif,nolink,70%);より，
#ref(ls_lu.eq6.gif,nolink,70%)

&ref(ls_lu.eq7.gif,nolink,70%);とすると，
#ref(ls_lu.eq8.gif,nolink,70%)

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

***LU分解の手順 [#zc7820d9]
係数行列AをLU分解するためには，以下の式を1行目から順番に適用していく．
#ref(ls_lu.eq11.gif,nolink,70%)

メモリを節約するために，LU分解した結果を以下のように1つの行列に格納する．
#ref(ls_lu.eq12.gif,nolink,70%)


***LU分解の実装 [#p0a9ff83]
LU分解をC++で実装した例を以下に示す．
#code(C){{
/*!
 * 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次方程式を解く手順 [#af94770d]
LU分解で連立1次方程式を解く方法は上でも述べたがここではより具体的な手順について説明する．
LU分解で連立1次方程式を解く手順は以下である．
+下三角行列Lと上三角行列Uを求める．
+前進代入により&ref(ls_lu.eq9.gif,nolink,70%);を求める．
+後退代入により&ref(ls_lu.eq10.gif,nolink,70%);を求める．
1はすでに述べているので2,3について説明する．

-前進代入
下三角行列Lと定数ベクトルbからyを求める．
#ref(ls_lu.eq13.gif,nolink,70%)
Lは下三角行列であるので，1行目は単純に&ref(ls_lu.eq14.gif,nolink,70%);で&ref(ls_lu.eq15.gif,nolink,70%);が算出される．
得られた&ref(ls_lu.eq15.gif,nolink,70%);を2行目の式に代入して同様に計算すると&ref(ls_lu.eq16.gif,nolink,70%);が得られる．これをn行目まで繰り返す．式で書くと以下となる．
#ref(ls_lu.eq17.gif,nolink,70%)
ここで，&ref(ls_lu.eq18.gif,nolink,70%);である．

-後退代入
前進代入で得られたyと上三角行列Uからxを求める．
#ref(ls_lu.eq19.gif,nolink,70%)
基本的にはガウスの消去法の後退代入と同じである．
ただし，各行は対角項で正規化されている(対角成分がすべて1)．
前進代入とは逆にn行目から順番に&ref(ls_lu.eq20.gif,nolink,70%);を求めていく．
#ref(ls_lu.eq21.gif,nolink,70%)
ここで，&ref(ls_lu.eq18.gif,nolink,70%);である．


***LU分解を用いた連立1次方程式の解法の実装 [#t6c3912c]
LU分解で連立1次方程式を解くコードをC++で実装した例を以下に示す．
#code(C){{
/*!
 * 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;
}
}}
トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS