1
LU分解を用いた連立一次方程式
LU分解による連立一次方程式の効率的な解法、C言語での
実装
【学習内容の構成】
1. LU分解の原理：係数行列Aを下三角行列Lと上三角行列
Uに分解
2. 解法の手順：Ly=bを解いた後、Ux=yを解くことでxを求め
る
3. プログラム実装：L行列・U行列の生成、前進代入・後退代
入によるC言語実装
4. ピボッティング：対角要素が0または小さい値の場合の行
入れ替えによる誤差対策
• 前提：行列演算の基礎、C言語プログラミング
• 意義：Gaussの消去法や逆行列法より高速な数値計算手法
の習得
2
3
LU分解
• 連立一次方程式を解く手段
• 通常，連立方程式を解く場合は変数を減らしていく
という方針で解いていくが，これもその１つ
• Gaussの消去法とLU分解は別物
– 計算機で実装した場合Gaussの消去法はLU分解に比べ
計算速度がかなり「遅い」
– 逆行列を求めるやり方も「遅い」
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
(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 (2)
5
上三角行列，下三角行列
• LU分解では，行列Aを下三角行列Lと上三角行
列Uの二つに分解する
• 三角行列であれば解を求めることが容易
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 (3)
6
LU分解で連立一次方程式を解く(1)
• 一次方程式 Ax = b を解く
• 行列Aを下三角行列Lと上三角行列Uの二
つにLU分解する
A = LU
• 行列AがLU分解されると求める方程式は
以下のようになる
Ax = LUx = b
7
LU分解で連立一次方程式を解く(2)
• y = Uxとおいて,まず
Ly = b
を解く. (三角行列を係数とする方程式は容
易に解くことができる)
• yを求めたら
Ux = y
を解くことで,ｘを求めることができる
8
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)
ここで
LUx = b (5)
Lc = b (6)
となるcを求める．
(4)
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 = ∑lk,i・ci where(0≦k≦n-1) (8)
最後にUx = cをxについて解くと
(7)
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
=
(9)
i=0
k-1
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) (10)
i=1
k-1
11
プログラム例
#include <stdio.h>
#define M 4
main(int argc, char *argv[])
{
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;
fp = fopen(argv[1], "r"); /* ファイルより入力行列を読み込む */
for(i = 0; i < M; i++){
for(j = 0; j < M; j++){
fscanf(fp,"%lf",&a[i][j]);
}
fscanf(fp,"%lf",&b[i]);
}
fclose(fp);
12
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];
}
}
13
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]);
}
14
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]);
}
}
15
実行結果
3 2 6 1 17
2 4 1 6 23
5 4 1 3 23
3 2 5 6 26
入力行列
3.00000 2.00000 6.00000 1.00000 17.00000
2.00000 4.00000 1.00000 6.00000 23.00000
5.00000 4.00000 1.00000 3.00000 23.00000
3.00000 2.00000 5.00000 6.00000 26.00000
L行列
1.00000 0.00000 0.00000 0.00000
0.66667 1.00000 0.00000 0.00000
1.66667 0.25000 1.00000 0.00000
1.00000 0.00000 0.12121 1.00000
U行列
3.00000 2.00000 6.00000 1.00000
0.00000 2.66667 -3.00000 5.33333
0.00000 0.00000 -8.25000 0.00000
0.00000 0.00000 0.00000 5.00000
x0 = 2.00000
x1 = 1.50000
x2 = 1.00000
x3 = 2.00000
入力行列として次の行列を与える
16
ピボッティング
• もしAi,iやuk,k等の対角上の要素が0であれば，前式
より，0徐算となり解が求まらない
• これを避けるために，0である要素を含む行と0で
ない要素を含む行を入れ替える（ピボッティング）
– しかし、0でなくても，0に近い小さな値で割ったときは丸
め誤差が大きくなる．この誤差を小さくするために，0でな
くても，最も絶対値の大きな値を含む行と最も絶対値の
小さな値を含む行 を入れ替えるのがよい
17
課題
• ピボッティングを行なうプログラムを作りなさい
• 実際に，対角要素の一部が０である行列を読み込
ませて，動作を確認しなさい
