LU分解
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
これまで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 と 上三...
#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%);を算出することができ...
&ref(ls_lu.eq9.gif,nolink,70%);が求まったらそれを2番目の...
Uは上三角行列であるので，ガウスの消去法と同じく後退代入(b...
していくことで解&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)と...
* - 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...
得られた&ref(ls_lu.eq15.gif,nolink,70%);を2行目の式に代入...
#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...
#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 vec...
{
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;
}
}}
終了行:
これまで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 と 上三...
#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%);を算出することができ...
&ref(ls_lu.eq9.gif,nolink,70%);が求まったらそれを2番目の...
Uは上三角行列であるので，ガウスの消去法と同じく後退代入(b...
していくことで解&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)と...
* - 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...
得られた&ref(ls_lu.eq15.gif,nolink,70%);を2行目の式に代入...
#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...
#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 vec...
{
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;
}
}}
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
