LU分解
をテンプレートにして作成
home
>
サイトマップ
開始行:
RIGHT:寄稿:東條遼平
* 下三角行列と上三角行列との分解 [#b2d86abf]
次のようにある行列Aを下三角行列Lと上三角行列Uに分解するこ...
#ref(img173.png,nolink)
これは例えば行列Aが3×3の行列であれば次のようになります。
#ref(img174.png,nolink)
このように分解することのメリットは、連立1次方程式を解くと...
下三角行列と上三角行列に分解していれば後は簡単に解くこと...
#ref(img175.png,nolink)
のような連立1次方程式を解くことを考えてみます。LU分解を行...
#ref(img176.png,nolink)
とでき、ベクトルcを用いて、
#ref(img177.png,nolink)
とおくと、
#ref(img178.png,nolink)
となります。ここで例えば3次であれば、
#ref(img179.png,nolink)
のようになり、
#ref(img181.png,nolink)
を解くことと同じ事になります。するとベクトルcはすぐにわか...
今度は、
#ref(img177.png,nolink)
より、
#ref(img180.png,nolink)
を解けばベクトルxが求められます。これは先程と同じで、
#ref(img182.png,nolink)
のようになりますので、ベクトルxが簡単に求まります。
このように、行列を下三角行列と上三角行列に分解できれば連...
* 分解する方法 [#u5bc9d5a]
行列Aを下三角行列Lと上三角行列Uに分解すれば連立方程式が簡...
どうやったらそのように分解できるのでしょうか。とりあえず、
#ref(img183.png,nolink)
を計算してみる事にします。
#ref(img184.png,nolink)
どうやら対角成分より下と、それ以外で場合分けができそうで...
もちろん、3次のときにどうなったからといって、数学的にそ...
ある程度自明なところがあります。
一般的に行列のかけ算を考えると、
#ref(img185.png,nolink)
となりますが、ここで対角成分よりも下側とそれ以外で分ける...
#ref(img186.png,nolink)
であり、行列Lの対角成分が1であることを考えると、
#ref(img187.png,nolink)
とできるため、行列L、行列Uの成分は、
#ref(img188.png,nolink)
となります。
* 実装の流れ [#ze6d41bc]
以上で一般的にプログラムが出来そうですが、よく見ると下三...
上三角行列の下の部分も無駄です。 しかも行列Lの対角成分は...
問題ありません。そのため行列Lと行列Uの2つの行列をまとめ...
#ref(img189.png,nolink)
のように値を保持します。それでは3次のときの例で流れを見...
#ref(lu1.png,nolink)
今この行列は元の行列Aを表しています。
#ref(lu2.png,nolink)
このような順で計算すると、
#ref(lu3.png,nolink)
このようになり、この右下2×2の小行列を考えると先程の3×...
同じ状態になっています。そこでここまでを1セットと考え、...
#ref(lu4.png,nolink)
と計算すると、
#ref(lu5.png,nolink)
となってLU分解が完成します。
* 問題点 [#u4e0c255]
上のような流れでLU分解が進んでいきますが、それぞれの1セ...
行の交換をする必要がある場合があります。というのは上の流...
そのときの対角成分の左上の値でその列の残りの部分を割って...
そのため割る値がゼロであれば困りますし、ゼロでなくてもゼ...
誤差も大きくなります。そこで毎回1セットが始まるときに、...
最初の場合であれば、3行全ての行の内、1列目の要素の絶対...
交換し、2セット目であれば2行目以降の中から2列目の値の...
交換します。もちろん交換が行われれば、どことどこが交換さ...
わかるようにどこかにその情報を保持しておかねばなりません。
* ソースコード [#n60eeb58]
width = m->width;
height = m->height;
if((order = (int*)malloc(sizeof(int)*height)) == NULL){
fprintf(stderr, "Allocation Error\n");
exit(1);
}
for(i=0; i<height; i++) order[i] = i;
動的に整数の配列を作成し、orderにポインタを保持しておきま...
その後0,1,2・・・と初期化していますが、
これは行を交換したときにこの値も同じく交換しておくことで、
例えば2行目と3行目を交換すれば、0、2、1・・・となる...
d = i;
max = fabs(m->a[i*width + i]);
for(j=i+1; j<height; j++){
if(max < fabs(m->a[j*width + i])){
d = j;
max = fabs(m->a[j*width + i]);
}
}
if(d!=i){
temp = order[i];
order[i] = order[d];
order[d] = temp;
ExchangeRowvector(m, i, d);
}
if(!(m->a[i*width + i])) continue;
for文までは現在注目している列の中から絶対値が最大の行を探...
そして最大の行が現在の行でなければ現在の行と最大の行を入...
ExchangeRowvectorはi行目とd行目を交換する関数です。
交換しても対角成分がゼロであればその回はLU分解が終わって...
continueします。
肝心のLU分解の心臓部ですが、以上のような部分がごちゃごち...
上の肝じゃない部分を除くと次のようになります。
for(i=0; i<height-1; i++){
for(j=i+1; j<height; j++){
m->a[j*width + i] /= m->a[i*width + i];
for(k=i+1; k<width; k++){
m->a[j*width + k] -= m->a[j*width + i]*m->...
}
}
}
長々と書いて来ましたが、実際のLU分解の心臓部はたったのこ...
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
終了行:
RIGHT:寄稿:東條遼平
* 下三角行列と上三角行列との分解 [#b2d86abf]
次のようにある行列Aを下三角行列Lと上三角行列Uに分解するこ...
#ref(img173.png,nolink)
これは例えば行列Aが3×3の行列であれば次のようになります。
#ref(img174.png,nolink)
このように分解することのメリットは、連立1次方程式を解くと...
下三角行列と上三角行列に分解していれば後は簡単に解くこと...
#ref(img175.png,nolink)
のような連立1次方程式を解くことを考えてみます。LU分解を行...
#ref(img176.png,nolink)
とでき、ベクトルcを用いて、
#ref(img177.png,nolink)
とおくと、
#ref(img178.png,nolink)
となります。ここで例えば3次であれば、
#ref(img179.png,nolink)
のようになり、
#ref(img181.png,nolink)
を解くことと同じ事になります。するとベクトルcはすぐにわか...
今度は、
#ref(img177.png,nolink)
より、
#ref(img180.png,nolink)
を解けばベクトルxが求められます。これは先程と同じで、
#ref(img182.png,nolink)
のようになりますので、ベクトルxが簡単に求まります。
このように、行列を下三角行列と上三角行列に分解できれば連...
* 分解する方法 [#u5bc9d5a]
行列Aを下三角行列Lと上三角行列Uに分解すれば連立方程式が簡...
どうやったらそのように分解できるのでしょうか。とりあえず、
#ref(img183.png,nolink)
を計算してみる事にします。
#ref(img184.png,nolink)
どうやら対角成分より下と、それ以外で場合分けができそうで...
もちろん、3次のときにどうなったからといって、数学的にそ...
ある程度自明なところがあります。
一般的に行列のかけ算を考えると、
#ref(img185.png,nolink)
となりますが、ここで対角成分よりも下側とそれ以外で分ける...
#ref(img186.png,nolink)
であり、行列Lの対角成分が1であることを考えると、
#ref(img187.png,nolink)
とできるため、行列L、行列Uの成分は、
#ref(img188.png,nolink)
となります。
* 実装の流れ [#ze6d41bc]
以上で一般的にプログラムが出来そうですが、よく見ると下三...
上三角行列の下の部分も無駄です。 しかも行列Lの対角成分は...
問題ありません。そのため行列Lと行列Uの2つの行列をまとめ...
#ref(img189.png,nolink)
のように値を保持します。それでは3次のときの例で流れを見...
#ref(lu1.png,nolink)
今この行列は元の行列Aを表しています。
#ref(lu2.png,nolink)
このような順で計算すると、
#ref(lu3.png,nolink)
このようになり、この右下2×2の小行列を考えると先程の3×...
同じ状態になっています。そこでここまでを1セットと考え、...
#ref(lu4.png,nolink)
と計算すると、
#ref(lu5.png,nolink)
となってLU分解が完成します。
* 問題点 [#u4e0c255]
上のような流れでLU分解が進んでいきますが、それぞれの1セ...
行の交換をする必要がある場合があります。というのは上の流...
そのときの対角成分の左上の値でその列の残りの部分を割って...
そのため割る値がゼロであれば困りますし、ゼロでなくてもゼ...
誤差も大きくなります。そこで毎回1セットが始まるときに、...
最初の場合であれば、3行全ての行の内、1列目の要素の絶対...
交換し、2セット目であれば2行目以降の中から2列目の値の...
交換します。もちろん交換が行われれば、どことどこが交換さ...
わかるようにどこかにその情報を保持しておかねばなりません。
* ソースコード [#n60eeb58]
width = m->width;
height = m->height;
if((order = (int*)malloc(sizeof(int)*height)) == NULL){
fprintf(stderr, "Allocation Error\n");
exit(1);
}
for(i=0; i<height; i++) order[i] = i;
動的に整数の配列を作成し、orderにポインタを保持しておきま...
その後0,1,2・・・と初期化していますが、
これは行を交換したときにこの値も同じく交換しておくことで、
例えば2行目と3行目を交換すれば、0、2、1・・・となる...
d = i;
max = fabs(m->a[i*width + i]);
for(j=i+1; j<height; j++){
if(max < fabs(m->a[j*width + i])){
d = j;
max = fabs(m->a[j*width + i]);
}
}
if(d!=i){
temp = order[i];
order[i] = order[d];
order[d] = temp;
ExchangeRowvector(m, i, d);
}
if(!(m->a[i*width + i])) continue;
for文までは現在注目している列の中から絶対値が最大の行を探...
そして最大の行が現在の行でなければ現在の行と最大の行を入...
ExchangeRowvectorはi行目とd行目を交換する関数です。
交換しても対角成分がゼロであればその回はLU分解が終わって...
continueします。
肝心のLU分解の心臓部ですが、以上のような部分がごちゃごち...
上の肝じゃない部分を除くと次のようになります。
for(i=0; i<height-1; i++){
for(j=i+1; j<height; j++){
m->a[j*width + i] /= m->a[i*width + i];
for(k=i+1; k<width; k++){
m->a[j*width + k] -= m->a[j*width + i]*m->...
}
}
}
長々と書いて来ましたが、実際のLU分解の心臓部はたったのこ...
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
ページ名:
home
>
Modified by
物理のかぎプロジェクト
PukiWiki 1.4.5_1
Copyright © 2001-2005
PukiWiki Developers Team
. License is
GPL
.
Based on "PukiWiki" 1.3 by
yu-ji
Powered by PHP 5.3.29HTML convert time to 0.002 sec.