08数値計算法
2008鳥取大学工学部知能情報工学科
田中美栄子
May 21, 2008
ガウスの消去法
• 連立一次方程式を解く
①式 a11 x1 a12 x2 a13 x3 ... a1n xn b1
②式 a21 x1 a22 x2 a23 x3 ... a2 n xn b2
...
an1 x1 an 2 x2 an 3 x3 ... ann xn bn
a11
①式 新①式 x1 a12 ' x 2 a13' x 3 ... a1n ' x n b1'
②式ー①式×a12
②式 新②式 0 a 22' x 2 a 23' x 3 ... a 2n ' x n b2 '
③式 新③式 0 a 32 ' x 2 a 33' x 3 ... a 3n ' x n b3 '
ガウスの消去法
• 次は右下の（n-1）×（n)行列に注目
①式 a11 x1 a12 x2 a13 x3 ... a1n xn b1
②式 a21 x1 a22 x2 a23 x3 ... a2 n xn b2
...
an1 x1 an 2 x2 an 3 x3 ... ann xn bn
x 2 a 23' ' x 3 ... a 2n " x n b1"
②式ー①式×a12
a 23" x 3 ... a 2n " x n b2 "
ガウスの消去法(4)
• 更新後、元と等価な次の連立方程式ができる
• ②式以下に注目すると(n-1)次の連立方程式
• この最初の式②をその第一係数で割り、次の式の第一係数
を掛けて次の式から引く、という操作を繰り返すと元と等価な
右上3角連立方程式ができる
①式 x1 a12 x2 a13 x3 ... a1n xn b1
②式 x2 a23 x3 ... a2 n xn b2
a33 x3 ... a3 n xn b3
...
ｎ式 • これを下から順に解いてゆく
xn bn
ガウスの消去法(5)
最初から整理。係数だけ書くと次の操作となる。
a 11
a 21
a 12
a 22
a 13 . a 1n
a 23 . a 2 n
b1
b2
1 a 12
0 a 22
a 13
a 23
. a 1n
. a 2n
b1
b2
a 31
.
a 32
.
a 33 . a 3n
. . .
b3
.
0 a 32
0 .
a 33
.
. a 3n
. .
b3
.
a n1 a n 2
a n 3 . a nn
bn
0 a n2
a n 3 . a nn
bn
第k回目の操作は第k行目を a kk で割って、j=k,..,n
列目
( k 1)
: a kj
( k 1)
: b k
a kj
bk
(k)
(k)
/ a kk
/ a kk
(k)
(k)
ガウスの消去法(6)
• K+1行目以下は、 i=k+1,…,n行目について、
j=k,k+1,…,n列をそれぞれ a i ,k 倍して引くと
a
( k 1)
ij
: a
( k 1)
: bi
bi
(k)
ij
(k)
a ik a kj
(k)
a ik b k
(k)
( k 1)
( k 1)
ガウスの消去法(7)
更新k-1回目で係数は
( k 1)
( k 1)
( k 1)
( k 1)
( k 1)
a k ,k
a k ,k 1
a k ,k 2
. a k ,n
bk
( k 1)
( k 1)
( k 1)
( k 1)
( k 1)
a k 1,k
a k 1,k 1
a k 1,k 2
. a k 1,n
b k 1
a k 2,k
( k 1)
a k 2,k 1
.
a n ,k
1
0
( k 1)
a k 2,k 2
.
( k 1)
( k 1)
(k)
. a k 2,n
.
a n ,k 1
a k ,k 1
1
( k 1)
a n ,k 2
(k)
.
( k 1)
a k ,k 2
(k)
a k 1,k 2
bk2
.
. a n ,n
.
.
( k 1)
.
( k 1)
(k)
a k ,n
(k)
a k 1,n
0
0
1
. a k 2,n
.
.
.
.
.
0
0
0
.
a n ,n
( k 1)
(k)
bn
( k 1)
(k)
bk
(k)
b k 1
bk2
(k)
.
(k)
bn
(k)
ガウスの消去法(8)
• 後退代入で解を順に求める
x n bn
(n)
/ a n ,n
(n)
• あとはk=n-1,n-2,…,1の順に
x k bk
(n)
n
a
j k 1
(n)
k, j
xj
•
/** Gauss elimination :written by Mieko Tanaka May31, 2006 */
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
#include <stdio.h>
//#define N 4
//#define N 2
#define N 8
void main(){
double a[N][N+1]={{4.0, 2.0, 1.0, 0.0, 4.0,0.0, 2.0,0.0,3.0},
{1.0,6.0,0.0, 3.0, -1.0, 0.0, -3.0, 2.0, 0.0},
{1.0, -1.0, 2.0, -3.0, 4.0, 5.0, 0.0, 0.0, 3.0},
{2.0, 6.0, 0.0, 3.0,-1.0, 4.0, 0.0, 2.0,3.0},
{1.0, -1.0, 4.0, -3.0, 4.0, 5.0, 3.0, 0.0, 3.0},
{3.0, 6.0, 0.0, 3.0,-1.0, 4.0, 0.0, 2.0,3.0},
{-1.0, -1.0, -4.0, -3.0, 1.0, -2.0, 3.0, 0.0, 3.0},
{1.0, -1.0, -2.0, -3.0, 4.0, 1.0, 0.0, 7.0, 5.0},
};
// double a[N][N+1]={{1.0, 2.0, 3.0},{4.0,1.0,6.0}};
/*
double a[N][N+1]={{1.0,1.0,-3.0,-4.0,-1.0},
{2.0,1.0,5.0,1.0,5.0},
{3.0,6.0,-2.0,1.0,8.0},
{2.0,2.0,2.0,-3.0,2.0}};
*/
int i,j,k, count=0;
double x[N];
double p,q;
//pivot: arrange rows of matrix a to make diagonal nonzero
//forward elimination(always consider small matrix beginning a[k][k]
for(k=0;k<N;k++){
//for all k=0,1…,N-1
p=a[k][k];
//p memorizes a[k][k]
a[k][k]=1;
//Diagonal elements are set to 1
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
for(j=k+1;j<N+1;j++) {a[k][j]/=p;count++;} //Divide right elements (j=k+1,…,N)by p
for(i=k+1;i<N;i++) {
//for every i-th row (i=k+1,…
q=a[i][k]; //q memorizes leftmost element(k-th column)
for(j=k+1;j<N+1;j++) {a[i][j]-=q*a[k][j];count++;} //Subtract k-row times q
a[i][k]=0;
//k-th column below a[k][k]=0
}
}
//backward iteration to compute x[i]
for(i=N-1;i>=0;i--){
x[i]=a[i][N];count++;
for(j=N-1;j>i;j--) {
x[i]-=a[i][j]*x[j];count++;
}
printf("x[%d]=%f\n",i,x[i]);
printf("count=%d\n",count);
}
}
ガウスの消去法の問題点
• まず第一式の第一係数 a11 が０だったり、大変
小さいときは困る。もし a11 1020 のとき a11
で割ると、割らない係数が割った係数より小さ
すぎてアンダーフローしてしまう。ずっと先で
aがゼロになってもまた同じことが起こる。
kk
• ピボット選択（ｐｉｖｏｔ）：係数行列の第ｋ列の第
ｋ行目より下の中で絶対値最大のものを探し
それが第m行目にあればｋ行とｍ行を入れ替
える
正方行列aのLU分解をCrout法
を用いて行う
• Lのk列目とUのｋ行目を計算してすぐ、元のa行
列に上書きして行くようにすると、入力も出力
もaだけで済むのでプログラムが短くなる。
•
• #include <stdio.h>
• #include <math.h>
• #define N 3
• void crout(double a[N][N]);
•
• int main(void){
•
double a[N][N]={{2,4,6},{3,8,7},{5,7,21}};
行列の3角分解
• Ax=bを解くとき、A=LUと分解して、
Ax= LUx = Lｙ= b
• つまり
Lｙ= b を解いてから
Ux =ｙ を解く
• どちらも三角行列だから解きやすい
３角分解の方法
• 例
2 4 6 11
3 8 7
21
5 7 21 31
2 11
3 21
5 31
0
22
32
0 1 u12
0 * 0 1
33 0 0
4 11u12
8 21u12 22
7 31u12 32
u13
u23
1
6 11u13
7 21u13 22 u23
21 31u13 32u23
よって
11 2, 21 3, 31 5, u12 2, u13 3
３角分解の方法
• 例
2 4 6 11 0 0 u11 u12 u13
3 8 7
0 * 0 u 22 u 23
21
22
5 7 21 31 32 33 0 0 u 33
2 11u11
4 11u12
6 11u13
3 21u11 8 21u12 22 u 22 7 21u13 22 u23
5 31u11 7 31u12 32 u22 21 31u13 32u23
u11 u22 u33 1
11 2, 21 3, 31 5, u12 2, u13 3
12変数に対して9方程式なので３つは任意に取れる。
とすると
３角分解の方法
• 以下 22 (8 21u12 ) / u22 (8 3 2) / 1 2
32 (7 31u12 ) / u22 (7 5 2) / 1 3
u23 (7 21u13 ) / 22 (7 3 3) / 2 1
u 23 (21 31u13 32 u 23 ) / u 33
(21 5 3 ( 3) ( 1)) / 1 3
2 0 0
L 3 2 0
5 3 3
1 2 3
U 0 1 1
0 0 1
これを一般化すると
• Crout法(i=1,…,n, j=2,…,n）
k=2,3,…,n に対し
k 1
kk a kk km u mk
m 1
k 1
ik a ik imu mk
11 a11
u kj (a kj km u mj ) / kk
m1
なのでこれを使って
u1 j a1 j / 11 が求められた
m 1
k 1
u kk 1
i1 a i1
• 正方行列aのLU分解をCrout法を用いて行
うとき、Lのk列目とUのｋ行目を計算し
てすぐ、元のa行列に上書きして行くよう
にすると、入力も出力もaだけで済むので
プログラムが短くなる。
•
#include <stdio.h>
• #include <math.h>
• #define N 3
• void crout(double a[N][N]);
•
• int main(void){
•
double a[N][N]={{2,4,6},{3,8,7},{5,7,21}};
•
crout(a);
•
return 0;
• }
•
• void crout(double a[N][N]){
•
int i, j, k;
•
for(k=0;k<N;k++){
•
• /* a[i][k]=L[i][k] for (k<=i) */
• for(i=k;i<N;i++){
•
for(j=0;j<k;j++) a[i][k]-=a[i][j] * a[j][k];
•
}
• /* a[k][j]=U[k][j] for (k<j) */
•
•
for(j=k+1;j<N;j++)
•
for(i=0;i<k;i++) a[k][j] -= a[k][i] * a[i][j];
•
for(j=k+1;j<N;j++) a[k][j] /= a[k][k];
•
•
printf("\n");
•
for(i=0; i<N; i++){
•
for(j=0; j<N; j++) printf("%f ", a[i][j]);
•
printf("\n");
•
}
•
}
• }
•
• ＜問題＞ 出力を表示せよ。（k=0の時は
入力したa行列のままであるが、その後は
ｋ列目とｋ行目が順に書き換わり、最後
はＬとＵを出力している事に注意。）
© Copyright 2024 ExpyDoc
