LU分解による連立一次方程式の解法とC言語での実装。
【学習内容の構成】
1. LU分解の原理：行列Aを下三角行列Lと上三角行列U
に分解する手法
2. 連立一次方程式の解法：Lc=bとUx=cの2段階で解を
算出
3. プログラム実装：ファイルからデータを読み込みLU
分解を実行するCプログラム
4. ピボッティング：対角要素が0または小さい値の場
合に行を入れ替えて計算を安定化
• 前提：行列の基礎知識，C言語によるファイル入出力
• 意義：Gaussの消去法より高速な数値計算手法の習得
2
LU分解
• 連立一次方程式を解く手段
• 通常、連立方程式を解く場合は、変数を減らして
いくという方針で解いていくが、LU分解もその１
つ
• Gaussの消去法とLU分解の比較
• 計算機で実装した場合、Gaussの消去法はLU分解に比
べて計算速度がかなり「遅い」
• 逆行列を求める方法も「遅い」
•3
連立一次方程式の行列表現
• 一般に連立一次方程式は、
• これを行列表現すると、
•4
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
の二つに分解する。
•5
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
•6
LU分解で連立一次方程式を解く(2)
• Ux=c とおいて、まず、
Lc = b
を解く
• c を求めたら、
Ux = c
を解くことで、ｘが求まる
•7
• LU分解する。
•8
A
0,0
A
0,1
A
0,2
A
1,0
A
1,1
A
1,2
A
2,0
A
2,1
A
2,2
=
１
0 0
ｌ
1,0
１
0
ｌ
2,0
ｌ
2,1
１
u
0,0
u
0,1
u
0,2
0 u
1,1
u
1,2
0 0 u
2,2
A
00
=u
00
A
01
=u
01
A
02
=u
02
A
10
=l
10
*u
00
A
11
=l
10
*u
01
+u
11
A
12
=l
10
*u
02
+u
12
A
20
=l
20
*u
00
A
21
=l
20
*u
01
+l
21
*u
11
A
22
=l
20
*u
02
+l
21
*u
12
+u
22
u
00
=A
00
u
01
=A
01
u
02
=A
02
l
10
=A
10
/u
00
u
11
=A
11
-l
10
*u
01
u
12
=A
12
-l
10
*u
02
l
20
=A
20
/u
00
l
21
=(A
21
-l
20
*u
01
)/u
11
u
22
=A
22
-l
20
*u
02
-l
21
*u
12
u
00
=A
00
u
01
=A
01
u
02
=A
02
l
10
=A
10
/u
00
l
20
=A
20
/u
00
u
11
=A
11
-l
10
*u
01
u
12
=A
12
-l
10
*u
02
l
21
=(A
21
-l
20
*u
01
)/u
11
u
22
=A
22
-l
20
*u
02
-l
21
*u
12
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
u
0,k
= A
0,k
where(0≦k≦n-1)
l
i,0
= A
i,0
/u
0,0
where(1≦i≦n-1)
u
1,k
= A
1,k
-u
0,k
・l
1,0
where(1≦k≦n-1)
0
0
1
1
2
2
3
3
一般化すると、
•9
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・ci where(0≦k≦n-1)
i=0
k-1
(6) の Lc = b
•10
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 について解くと、
•11
プログラム例
入力行列は既にファイル（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の中身
3x
0
+ 2x
1
+ 6x
2
+ x
3
= 17
2x
0
+ 4x
1
+ x
2
+ 6x
3
= 23
5x
0
+ 4x
1
+ x
2
+3x
2
= 23
3x
0
+ 2x
1
+ 5x
2
+ 6x
3
= 26
これは、以下の方程式を表しているものとする。
•12
#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;
}
•13
実行結果
•14
LU分解できない場合
• Ai,iやuk,k等の対角上の要素が0であれば、0徐算が
発生しLU分解できない。
•15
0 3 4 1 20
3 2 2 2 19
5 6 2 3 27
4 8 5 2 14
例
ピボッティング(1)
• これを避けるために、0である要素を含む行と0で
ない要素を含む行を入れ替える（ピボッティン
グ）。
•16
3 2 2 2 19
5 6 2 3 27
4 8 5 2 14
0 3 4 1 20
0 3 4 1 20
3 2 2 2 19
5 6 2 3 27
4 8 5 2 14
例
ピボッティング(2)
• 0でなくても、0に近い小さな値で割ったときは丸
め誤差が大きくなる。
• この誤差を小さくするために、0でなくても、最も
絶対値の大きな値を含む行と最も絶対値の小さな
値を含む行 を入れ替える方がよい。
•17
解が求まらない場合
• 以下の例のような場合は、連立一次方程式は解く
ことができない。
•18
x
0
+ 2x
1
+ 2x
2
+ 2x
3
= 19
x
0
+ 2x
1
+ 2x
2
+ 2x
3
= 17
x
0
+ x
1
+ x
2
+2x
3
= 14
x
1
+ x
2
+ x
3
= 10
x
0
+ 2x
1
+ 2x
2
+ 2x
3
= 17
x
0
+ x
1
+ x
2
+2x
3
= 14
x
1
+ x
2
+ x
3
= 10
