LU分解
•連立一次方程式を解く手段
•通常、連立方程式を解く場合は、変数を減らして
いくという方針で解いていくが、LU分解もその１
つ
•Gaussの消去法とLU分解の比較
•計算機で実装した場合、Gaussの消去法はLU分解に比
べて計算速度がかなり「遅い」
•逆行列を求める方法も「遅い」
•2
連立一次方程式の行列表現
•一般に連立一次方程式は、
•これを行列表現すると、
•3
A0,0x0+A0,1x1・・・A0,n-1xn-1 = b0
A1,0x0+A1,1x1・・・A1,n-1xn-1 = b1
An-1,0x0+An-1,1x1・・・An-1,n-1xn-1 = bn-1
A0,0 A0,1 ・
・A0,n-1
A1,0 A1,1 ・
・A1,n-1
： ： ： ：
An-1,0 An-1,1 ・
・An-1,n-1
x0
x1
：
xn-1
b0
b1
：
bn-1
=Ax=b
上三角行列，下三角行列
•三角行列であれば解を求めることが容易。
•LU分解では，行列Aを下三角行列Lと上三角行列U
の二つに分解する。
•4
A0,0 A0,1 ・
・A0,n-1
A1,0 A1,1 ・
・A1,n-1
： ： ： ：
An-1,0 An-1,1 ・
・An-1,n-1
1 0 ・
・
・0 0
l1,0 1・
・
・0 0
： ： ： ：
ln-2,0 ln-2,1 ・
・
・1 0
ln-1,0 ln-1,1 ・
・
・ln-1,n-2 1
u0,0 u0,1 ・
・
・u0,n-2 u0,n-1
0 u1,1 ・
・
・u1,n-2 u1,n-1
： ： ： ：
0 0 ・
・
・un-2,n-2 un-2,n-1
0 0 ・
・
・0 un-1,n-1
A=
=
=LU
LU分解で連立一次方程式を解く(1)
•一次方程式 Ax = b を解く。
•行列Aを下三角行列Lと上三角行列Uの二つにLU分
解する。
A = LU
•行列AがLU分解されると求める方程式は以下のよ
うになる。
Ax = LUx = b
•5
LU分解で連立一次方程式を解く(2)
•Ux=c とおいて、まず、
Lc = b
を解く
•c を求めたら、
Ux = c
を解くことで、ｘが求まる
•6
•LU分解する。
•7
A0,0 A0,1 A0,2
A1,0 A1,1 A1,2
A2,0 A2,1 A2,2
=
１0 0
ｌ
1,0
１0
ｌ
2,0
ｌ2,1 １
u0,0 u0,1 u0,2
0 u1,1 u1,2
0 0 u2,2
A00=u00
A01=u01
A02=u02
A10=l10*u00
A11=l10*u01+u11
A12=l10*u02+u12
A20=l20*u00
A21=l20*u01+l21*u11
A22=l20*u02+l21*u12+u22
u00=A00
u01=A01
u02=A02
l10=A10/u00
u11=A11-l10*u01
u12=A12-l10*u02
l20=A20/u00
l21=(A21-l20*u01)/u11
u22=A22-l20*u02-l21*u12
u00=A00
u01=A01
u02=A02
l10=A10/u00
l20=A20/u00
u11=A11-l10*u01
u12=A12-l10*u02
l21=(A21-l20*u01)/u11
u22=A22-l20*u02-l21*u12
A0,0 A0,1 ・
・A0,n-1
A1,0 A1,1 ・
・A1,n-1
： ： ： ：
An-1,0 An-1,1 ・
・An-1,n-1
1 0 ・
・
・
0 0
l1,0 1・
・
・
0 0
： ： ： ：
ln-2,0 ln-2,1 ・
・
・
1 0
ln-1,0 ln-1,1 ・
・
・
ln-1,n-2 1
u0,0 u0,1 ・
・
・
u0,n-2 u0,n-1
0 u1,1 ・
・
・
u1,n-2 u1,n-1
： ： ： ：
0 0 ・
・
・
un-2,n-2 un-2,n-1
0 0 ・
・
・
0 un-1,n-1
u0,k = A0,k where(0≦k≦n-1)
li,0 = Ai,0/u0,0 where(1≦i≦n-1)
u1,k = A1,k-u0,k・l1,0 where(1≦k≦n-1)
0
0
1
1
2
2
3
3
一般化すると、
•8
1 0 ・
・
・
0 0
l1,0 1・
・
・
0 0
ln-2,0 ln-2,1 ・
・
・
1 0
ln-1,0 ln-1,1 ・
・
・
ln-1,n-2 1
c0
c1
cn-2
cn-1
b0
b1
bn-2
bn-1
=
これを c について解くと、
c0= b0
c1= b1-l1,0・c0
c2= b2-l2,0・c0-l2,1・c1
ck= bk-∑lk,i・ciwhere(0≦k≦n-1)
i=0
k-1
(6) のLc = b
•9
xn-1 = cn-1/un-1,n-1
xn-2 = (cn-2-un-2,n-1・xn-1)/un-2,n-2
xn-3 = (cn-3-un-3,n-2・xn-2-un-3,n-1・xn-1)/un-3,n-3
xn-k = (cn-k-∑un-k,n-i・xn-i)/un-k,n-k where(1≦k≦n)
i=1
k-1
u0,0 u0,1 ・
・
・
u0,n-2 u0,n-1
0 u1,1 ・
・
・
u1,n-2 u1,n-1
0 0 ・
・
・
un-2,n-2 un-2,n-1
0 0 ・
・
・
0 un-1,n-1
c0
c1
cn-2
cn-1
x0
x1
xn-2
xn-1
=
最後に (7) のUx = c
これを x について解くと、
•10
プログラム例
入力行列は既にファイル（c:\lu_test.txt）
にあるものとする
3 2 6 1 17
2 4 1 6 23
5 4 1 3 23
3 2 5 6 26
c:\lu_data.txtの中身
3x0+ 2x1+ 6x2+ x3= 17
2x0+ 4x1+ x2+ 6x3= 23
5x0+ 4x1+ x2+3x2= 23
3x0 + 2x1+ 5x2+ 6x3= 26
これは、以下の方程式を表しているものとする。
•11
#include "stdio.h"
#pragma warning(disable:4996)
int main()
{const int M = 4;
double a[M][M];
double b[M];
double c[M];
double l[M][M];
double u[M][M];
double x[M];
int i,j,k;
FILE *fp;
char line[1000];
int ch;
fp = fopen("c:\\lu_data.txt", "r"); /* データの読み込み */
for(i = 0; i < M; i++){
for(j = 0; j < M; j++){
fscanf_s(fp,"%lf",&a[i][j]);
}
fscanf_s(fp,"%lf",&b[i]);
}
fclose(fp);
for(i = 0; i < M; i++){ /* L行列、U行列を1と0で初期化 */
for(j = 0; j < M; j++){
u[i][j] = 0;
if(i == j)
l[i][j] = 1;
else
l[i][j] = 0;
}
}
for(i = 0; i < M; i++){
for(j = i; j < M; j++){ /* U行列の生成 */
u[i][j] = a[i][j];
for(k = 0; k < i; k++){
u[i][j] -= u[k][j] * l[i][k];
}
}
for(j = i+1; j < M; j++){ /* L行列の生成 */
l[j][i] = a[j][i];
for(k = 0; k < i; k++){
l[j][i] -= u[k][i] * l[j][k];
}
l[j][i] /= u[i][i];
}
}
for(i = 0; i < M; i++){ /* c行列の生成 */
c[i] = b[i];
for(j = 0; j < i; j++){
c[i] -= l[i][j] * c[j];
}
}
for(i = M-1; i >= 0; i--){ /* x行列の生成 */
x[i] = c[i];
for(j = M-1; j > i; j--){
x[i] -= u[i][j] * x[j];
}
x[i] /= u[i][i];
}
printf("入力行列\n"); /* 入力行列の出力 */
for(i = 0; i < M; i++){
for(j = 0; j < M; j++){
printf("%10.5lf ",a[i][j]);
}
printf("%10.5lf\n",b[i]);
}
printf("\nL行列\n"); /* L行列の出力 */
for(i = 0; i < M; i++){
for(j = 0; j < M; j++){
printf("%10.5lf ",l[i][j]);
}
printf("\n");
}
printf("\nU行列\n"); /* U行列の出力 */
for(i = 0; i < M; i++){
for(j = 0; j < M; j++){
printf("%10.5lf ",u[i][j]);
}
printf("\n");
}
printf("\n");
for(i = 0; i < M; i++){ /* 解の出力 */
printf("x%d = %10.5lf\n",i,x[i]);
}
printf( "Enter キーを1,2回押してください. プログラムを終了します\n");
ch = getchar();
ch = getchar();
return 0;
} •12
実行結果
•13
LU分解できない場合
•Ai,iやuk,k等の対角上の要素が0であれば、0徐算が
発生しLU分解できない。
•14
0 3 4 1 20
3 2 2 2 19
5 6 2 3 27
4 8 5 2 14
例
ピボッティング(1)
•これを避けるために、0である要素を含む行と0で
ない要素を含む行を入れ替える（ピボッティン
グ）。
•15
3 2 2 2 19
5 6 2 3 27
4 8 5 2 14
03 4 1 20
03 4 1 20
3 2 2 2 19
5 6 2 3 27
4 8 5 2 14
例
ピボッティング(2)
•0でなくても、0に近い小さな値で割ったときは丸
め誤差が大きくなる。
•この誤差を小さくするために、0でなくても、最も
絶対値の大きな値を含む行と最も絶対値の小さな
値を含む行 を入れ替える方がよい。
•16
解が求まらない場合
•以下の例のような場合は、連立一次方程式は解く
ことができない。
•17
x0+ 2x1+ 2x2+ 2x3= 19
x0+ 2x1+ 2x2+ 2x3= 17
x0+ x1+ x2+2x2= 14
x1+ x2+ x3= 10
x0+ 2x1+ 2x2+ 2x3= 17
x0+ x1+ x2+2x2= 14
x1+ x2+ x3= 10
