数値計算法I

数値計算法
第6回
反復解法
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
ガウスの消去法の問題点
• ガウスの消去法
– 特異行列で無い限り必ず解ける
– 方程式の数が大きくなると計算回数が増加する。
3
n
に比例 nは未知数の数（元数）
3
• 反復解法
– ある程度の戦略をもてば、あてずっぽうにやっ
 ていった方が早く計算を終える時もある。
2
ヤコビの反復法
a11x1 a12x2 a13x3 a14x4  b1

a21x1 a22x2 a23x3 a24x4  b2

a31x1 a32x2 a33x3 a34x4  b3

a41x1 a42x2 a43x3 a44x4  b4
行列の対角項aii≠0の時、解を以下の繰り返しで求める。


(k )
(k )
(k )
b

a
x

a
x

a
x
13 3
14 4
x1(k1)  1 12 2
a11

(k )
(k )
(k )

b

a
x

a
x

a
x
23 3
24 4
x2(k1)  2 21 1

a22

 (k1) b3  a31x1(k )  a32x2(k )  a34x4(k )

x3
a33

(k )
(k )
(k )
 (k1) b4  a41x1  a42x2  a43x3
x4


a44

3
ヤコビの反復法

(k )
(k )
(k )
b

a
x

a
x

a
x
13 3
14 4
x1(k1)  1 12 2
a11

(k )
(k )
(k )

b

a
x

a
x

a
x
23 3
24 4
x2(k1)  2 21 1

a22

 (k1) b3  a31x1(k )  a32x2(k )  a34x4(k )

x3
a33

(k )
(k )
(k )
 (k1) b4  a41x1  a42x2  a43x3
x4


a44

収束判定は


(1)
n
x(k1)  x(k )  
i
i
i1
は十分に小さな数（自分で決める）
4
ヤコビの反復法の収束条件
Ax  b
(2)
係数行列Aを下三角行列、上三角行列、対角行列の和で表す
A L DU
(3)
a11
a12


a21 a22

A


an1 an2

a11 0

0 a22

D


0
 0
a1n 

a2n 


ann
 0
0

a21 0

L


an1 an2
0 a12
0 


0  U  0 0
 



0 0
ann
0

0


0
a1n 

a2n 


0 
5
ヤコビの反復法の収束条件
(1)式の繰り返しは、以下の式で表せる。
(k1)
(k)
x
 Hx  z
(4)
ここで、
(5)
H  D1(LU)

1
zD b
(6)
反復回数Kで真値に達したとすると、

(7)
x(K)  Hx(K)  z

(4)-(7)より、
(k1)
(K)
(k )
(K)
x
 x  H(x  x )
k回目の誤差
 k+1回目の誤差
6
ヤコビの反復法の収束条件
x
(k1)
x
(K)
 H(x
k+1回目の誤差
(k )
x
(K)
)
k回目の誤差
よって、誤差は反復回数が一回増える毎に
 H倍になる。Hが全体的に1より小さければよい。
そのためには、

aij
1
ji aii
(i 1,2, ,n)
行列の対角要素aiiが他の要素より遥かに
大きければよい。対角優位

7
ヤコビの反復法の収束条件
より正確には、ヤコビの反復法で収束するためには、
行列Hの全ての固有値の絶対値が1より小さくないと
いけない。
収束の早さは、絶対値が最大の固有値λmaxの値で決ま
る。最大の固有値が小さい程、収束が早い。
max 1

8
ヤコビの反復法の例
 4x1  x2  x3  4

2x1  6x2  x3  8
x  2x 5x  7
2
3
 1
4

A  2

1

1
6
2
1

1
5

0 1 1
4 0 0




D  0 6 0 LU  2 0 1




1
2
0
0
0
5




9
ヤコビの反復法の例

0

1
4 0 0 0 1 1 

 
 1
1
H J  D (LU)  0 6 0 2 0 1
3




0 0 5 1 2 0 1

5
1 1

4 4
1
0
6
2 
0
5 
HJの固有値は、0.5281079,-0.26405±0.19156i
なので収束する。
10
ヤコビの反復法の例
回数誤差
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
2.51667
0.80000
0.26167
0.12264
0.06761
0.03543
0.01863
0.00990
0.00520
0.00276
0.00145
0.00077
0.00041
0.00021
0.00011
0.00006
x1
x2
x3
1.00000
1.25000
1.01667
1.05000
0.99292
1.00931
0.99650
1.00216
0.99891
1.00058
0.99969
1.00016
0.99991
1.00004
0.99998
1.00001
-1.00000
-1.66667
-2.01667
-1.97500
-2.01722
-1.99431
-2.00449
-1.99814
-2.00113
-1.99944
-2.00030
-1.99984
-2.00008
-1.99996
-2.00002
-1.99999
0.00000
1.60000
1.81667
2.00333
1.98000
2.00831
1.99586
2.00249
1.99883
2.00067
1.99966
2.00018
1.99990
2.00005
1.99997
2.00001
11
Jacobiの反復法の例
implicit real*8 (a-h,o-z)
parameter(nmax=3)
dimension a(nmax,nmax)
real*8 lu
dimension lu(nmax,nmax),dinv(nmax,nmax)
dimension h(nmax,nmax)
dimension x(nmax),b(nmax),z(nmax),xnew(nmax)
do i=1,n
do j=1,n
if(i.ne.j)then
lu(i,j)=a(i,j)
dinv(i,j)=0
else
lu(i,j)=0
dinv(i,j)=1/a(i,j)
end if
end do
end do
do i=1,n
z(i)=b(i)*dinv(i,i)
end do
c
a(1,1)=4
a(1,2)=1
a(1,3)=1
a(2,1)=2
a(2,2)=6
a(2,3)=1
a(3,1)=1
a(3,2)=2
a(3,3)=5
b(1)=4
b(2)=-8
b(3)=7
x(1)=1
x(2)=-1
x(3)=0
c
n=3
itermax=100
errlim=1.d-4
c
do i=1,n
do j=1,n
sum=0
do ii=1,n
sum=sum+dinv(i,ii)*lu(ii,j)
end do
h(i,j)=-sum
end do
end do
c
c
iter=0
10 continue
iter=iter+1
do i=1,n
xnew(i)=z(i)
do j=1,n
xnew(i)=xnew(i)+h(i,j)*x(j)
end do
end do
c
error=0
do i=1,n
error=error+abs(xnew(i)-x(i))
end do
c
write(6,200)iter,error,(x(i),i=1,n)
200 format(i4,4f12.5)
if(error.lt.errlim)goto 100
if(iter.lt.itermax)then
do i=1,n
x(i)=xnew(i)
end do
goto 10
end if
100 continue
stop
end
12
Gauss Seidel法
a11x1 a12x2 a13x3 a14x4  b1

a21x1 a22x2 a23x3 a24x4  b2

a31x1 a32x2 a33x3 a34x4  b3

a41x1 a42x2 a43x3 a44x4  b4
行列の対角項aii≠0の時、解を以下の繰り返しで求める。


(k )
(k )
(k )
b

a
x

a
x

a
x
13 3
14 4
 x1(k1)  1 12 2
a11

(k1)
(k )
(k )

b

a
x

a
x

a
x
23 3
24 4
 x2(k1)  2 21 1

a22

 (k1) b3  a31x1(k1)  a32x2(k1)  a34x4(k )

 x3
a33

(k1)
 a42x2(k1)  a43x3(k1)
 (k1) b4  a41x1
x4


a44

13
Gauss-Seidel法
• Jacobiの方法では、k+1番目の解(x1,x2,x3,--xn)を全
て、k番目の解に基づいて決める。
• Gauss-Seidel法では、k+1番目の解(x1,x2,x3,--xn)を
求める際、k+1番目が既に求まっているものは、
それを使う。
– より最新の解を使って計算する
– 収束が速い
• Jacobiの反復法で収束するなら、Gauss-Seidel法
でも収束する
• 通常は、Jacobi法を使うくらいなら、GaussSeidel法を使う。
14
Gauss-Seidel法の例
 4x1  x2  x3  4

2x1  6x2  x3  8
x  2x 5x  7
2
3
 1
回数誤差
1
2
3
4
5
6
7
2.85000
0.63333
0.05681
0.00877
0.00106
0.00011
0.00001
x1
x2
1.00000 -1.00000
1.25000 -1.75000
0.97500 -1.96667
0.99375 -1.99653
0.99917 -1.99970
0.99991 -1.99998
0.99999 -2.00000
x3
0.00000
1.85000
1.99167
1.99986
2.00005
2.00001
2.00000
15
Gauss-Seidel法の例
implicit real*8 (a-h,o-z)
parameter(nmax=3)
dimension a(nmax,nmax)
dimension x(nmax),b(nmax),xold(nmax)
c
a(1,1)=4
a(1,2)=1
a(1,3)=1
a(2,1)=2
a(2,2)=6
a(2,3)=1
a(3,1)=1
a(3,2)=2
a(3,3)=5
b(1)=4
b(2)=-8
b(3)=7
x(1)=1
x(2)=-1
x(3)=0
c
n=3
itermax=100
errlim=1.d-4
c
do i=1,n
xold(i)=x(i)
end do
iter=0
c
10 continue
iter=iter+1
c
do i=1,n
rhs=b(i)
do j=1,n
if(i.ne.j)then
rhs=rhs-a(i,j)*x(j)
end if
end do
x(i)=rhs/a(i,i)
end do
c
error=0
do i=1,n
error=error+abs(x(i)-xold(i))
end do
c
write(6,200)iter,error,(xold(i),i=1,n)
200 format(i4,4f12.5)
if(error.lt.errlim)goto 100
if(iter.lt.itermax)then
do i=1,n
xold(i)=x(i)
end do
goto 10
end if
100 continue
stop
end
16
Jacobi法とGauss-Seidel法の収束の比較
101
Gauss-Seidel
Jacobi
100
10-1
-2
10
誤差
10-3
10-4
10-5
10-6
0
5
10
反復回数
15
20
17
加速緩和法
• Gauss-Seidel法よりも収束を更に早める方法
• Gauss-Seidel法でxik+1を求める。
• xikとxik+1の差を使って、次式で補正したyiを使っ
て、次のステップに進む。
x
x
x
x
y1  x1(k)  

y2  x2(k)  

(k)
y

x
 3
3 

(k)
y

x
4 
 4




(k 1)
1
 x1(k)
(k 1)
2
 x2(k)
(k 1)
3
 x3(k)
(k1)
4
 x4(k)
αは加速係数
α=1でGauss-Seidel法と一致する。
0<α<2の範囲で最適な加速係数を決める
18
加速緩和法の例
 4x1  x2  x3  4

2x1  6x2  x3  8
x  2x 5x  7
2
3
 1
α=0.9の時、
回数誤差
 1 2.56500
2
3
4
5
6
7
0.78150
0.09553
0.01632
0.00252
0.00033
0.00004
x1
1.00000
1.22500
1.00000
0.99438
0.99869
0.99979
0.99997
x2
-1.00000
-1.67500
-1.93750
-1.99062
-1.99879
-1.99986
-1.99999
x3
0.00000
1.66500
1.95900
1.99577
1.99962
1.99997
2.00000
19
加速緩和法の例
α=1.5の時、
回数誤差
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
4.27500
2.01250
1.09146
0.53257
0.26777
0.13375
0.06688
0.03344
0.01672
0.00836
0.00418
0.00209
0.00105
0.00052
0.00026
0.00013
0.00007
x1
x2
x3
1.00000 -1.00000 0.00000
1.37500 -2.12500 2.77500
0.77500 -1.88750 1.60000
1.10313 -2.05104 2.19979
0.94719 -1.97403 1.90017
1.02628 -2.01295 2.04993
0.98685 -1.99352 1.97504
1.00657 -2.00324 2.01248
0.99671 -1.99838 1.99376
1.00164 -2.00081 2.00312
0.99918 -1.99960 1.99844
1.00041 -2.00020 2.00078
0.99979 -1.99990 1.99961
1.00010 -2.00005 2.00020
0.99995 -1.99997 1.99990
1.00003 -2.00001 2.00005
0.99999 -1.99999 1.99998
20
加速係数
 4x1  x2  x3  4

2x1  6x2  x3  8
x  2x 5x  7
2
3
 1
誤差が10-4より小さくなるのに要する反復回数
20

Gauss-Seidel
15
反復回数
10
5
0
0.4
0.6
0.8
1
1.2
加速係数、α
1.4
1.6
最適なαの値は、解くべき式によって異なる。
21
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
• 今までに習ったニュートン法と行列計算で、実際の研究
上の問題の多くに対応できる。
• 例：反応計算
• 以下の反応があったとする
– 2H2+O2<->2H2O
– この反応は反応係数が温度に依存する。ある温度での
H2,O2,H2Oの割合は以下の3個の連立式で成り立つ
ax12 x2  bx32  0

2x2  x3  c
2x  2x  d
3
 1
22
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
ax12 x2  bx32  0

2x2  x3  c
2x  2x  d
3
 1
x1:H2分子の数
x2:O2分子の数
x3:H2O分子の数
a:2H2+O2->2H2Oの反応係数
b:2H2O->2H2+O2の反応係数
c:O2原子の数
d:H2原子の数
23
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
ax12 x2  bx32  0

2x2  x3  c
2x  2x  d
3
 1
• 上の連立方程式は、そのままでは行列の形に表せない。
• Newton法を用いるために、
 f (x1, x2 , x3 )  ax12 x2  bx32

g(x1, x2 , x3 )  2x2  x3  c
h(x , x , x )  2x  2x  d
1
3
 1 2 3
とおいて、3個の関数、f,g,hが全て0になる
x1,x2,x3の組み合わせを探す。
24
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
• 各関数をx1,x2,x3の周りにTaylor展開して、

f
f
f
 f (x1   x1, x2   x2 , x3   x3 )  f (x1, x2 , x3 )  x  x1  x  x2  x  x3 
1
2
3


g
g
g
 x1 
 x2   x3 
g(x1   x1, x2   x2 , x3   x3 )  g(x1, x2 , x3 ) 
x1
x2
x3


h
h
h
 x1 
 x2   x3 
h(x1   x1, x2   x2 , x3   x3 )  h(x1, x2 , x3 ) 
x1
x2
x3

2次以上の項を無視し、左辺を0とすると、
25
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
 f
f
f
 x  x1  x  x2  x  x3   f (x1, x2 , x3 )
2
3
 1
 g
g
g
 x2   x3  g(x1, x2 , x3 )
  x1 
x2
x3
 x1
 h
h
h
 x2   x3  h(x1, x2 , x3 )
  x1 
x2
x3
 x1
• kステップ目からk+1ステップ目に動かす解の移動
量δx1,δx2,δx3を求める連立一次方程式になる。
• 但し、係数はその都度変化する。
26
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
 f (x1, x2 , x3 )  ax12 x2  bx32

g(x1, x2 , x3 )  2x2  x3  c
h(x , x , x )  2x  2x  d
1
3
 1 2 3
において、
 f
f
2 f
 x  2ax1x2 , x  ax1 , x  2bx3
2
3
 1
 g
g
g
 0,
 2,
1

x2
x3
 x1
 h
h
h
 2,
 0,
2

x2
x3
 x1
27
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
δx1,δx2,δx3は以下の式により求まる。

2ax1x2 x1  ax12 x2  2bx3 x3   ax12 x2  bx32

2 x2   x3   2x2  x3  c

2 x1  2 x3   2x1  2x3  d 

行列形式に直すと、
2ax1k x2k ax1k2 2bx3k   x1 
 ax1k2 x2k  bx3k2 

 
 k

k
2
1   x2     2x2  x3  c 
 0
k
k
 2




0
2   x3 
2x1  2x3  d 
28
高次連立方程式の解法
（ニュートン法と行列計算の組み合わせ例）
2ax1k x2k ax1k2 2bx3k   x1 
 ax1k2 x2k  bx3k2 

 
 k

k
2
1   x2     2x2  x3  c 
 0
 2
2x1k  2x3k  d 
0
2   x3 
• 行列の係数や、右辺はk番目のx1,x2,x3とその他の定数
a,b,c,dで決まる。上の行列から、 δx1,δx2,δx3を求め
て、k+1番目の値を以下のように決める。
x1k1  x1k  x1
 k1
k
x2  x2  x2
k
xk1

x
3  x 3
 3
k+1番目の値とk番目の値が変わらなくなるまで計算を行う。
29

Download Report