2014年7月 2日（水） - 青柳研究室 - 九州大学

2014年7月2日（水）
数値解析講義
前期水曜３時限
情報理学コース
青柳睦
aoyagi@cc.九州大学のドメイン
研究室
箱崎キャンパス情報基盤研究開発センター５階５０２，５００
伊都キャンパスウエスト2号館 1006室右ドア側
1
2．連立一次方程式の解法
:
2.5
2.6
2.7
2.8
2.9
講義内容
｛特異値分解法｝
共役勾配法 →演習
反復法 2.7.1 ヤコビ法 →演習
ガウス・ザイデル法
｛ＳＯＲ法（概要のみ解説）｝
3．固有値と固有ベクトルの数値計算
3.1 固有値問題の序論
3.2 ヤコビの方法（対角行列への帰着） →演習
3.3 ギブンスの方法（３重対角行列への帰着）→演習
3.4 DKA法（３重対角行列の固有値の計算 )
3.5 ｛ハウスホルダー法（３重対角行列への帰着）｝
今日の講義
3.6 ｛バイセクション法（３重対角行列の固有値の計算）｝
4. 数値計算の応用
4.1 差分法による熱伝導方程式の数値計算
4.2 ・・・
2
A1 ≡ A
対角形に収束す
るまで繰り返す

A2 = P −1 A1 P 
 相似変換の繰り返し
−1
A3 = P A2 P 


一旦、３重対角
形に変換する
3.2 ヤコビ法
*



*
0




*







*
0


* 

対角行列
* *


0 
* * *


* * 


   

 0
 * *


*
*


3.5 ハウスホルダー法
3重対角行列
3.4 DKA法
対角要素→固有値
3.3 ギブンス法
3.6 バイセクション法
固有値は３重対角行列に適した
別の反復法により求める
3
3．5 ハウスホルダー法
（参考程度）
ギブンス法にならぶ三重対角化の代表的な方法で，演算回数が
行列の次数（ｎ）と明確に関係付けられる．
A : n × n 対称行列
変換
* *


0 
* * *


* * 







 0
 * *


*
*


3重対角行列
三重対角化のプロセスは固有値問題以外にも有効に使われ，
「ハウスホルダー変換」とも呼ばれる．
まず，最初のステップ・・
相似変換 B = P −1 AP を考える
*

*
B = 0


0

* 0  0

* *  *
* *  *

   
* *  * 
行または列を一括で0に
なるようにする
c.f. ギブンス法は要素毎
変換行列に次の形を仮定する
P= I − 2ut u
u = t (u1 , u2 , , un ) : n 次元ベクトル
u=
uu
=
2
2
t
n
2
=
u
∑ i 1 : 正規化
i =1
P は対称な直交行列であることを示す。
t
t
t
t
( I 2ut u) =
P =−
I − 2t (ut u) =
I − 2ut u =
P
( I − 2ut u)( I − 2ut u)
PP =
PP =
I − 2ut u − 2ut u + 4ut uu t u
=
I − 2ut u − 2ut u + 4u(t uu) t u
=
I − 2ut u − 2ut u + 4u t u =
I
=
∴ P = t P = P −1
P = P −1 は対称だから
t
1
=
B t ( P −=
AP)
t
−1
−1
PA t=
( P −1 ) PAP
=
P=
AP B
よって B = P −1 AP も対称
u として
 0  第1成分が0
 
*

u=

 
*
の形のものを選ぶことができたと仮定すると，
このとき P の形は，
P= I − 2ut u
1

0
−1

=
P P=


0
0  0

*  *
  

*  *
ここで，任意の行列
GP
G に対して
 g11 g12  g1n   1 0  0 



g
g
g

0
*
*

21
22
2n  
=

 
        

 

 g n1 g n 2  g nn   0 *  * 
 g11 *  * 


*
*
g

 21

    


 g n1 *  * 
GP の1列目はG の1列目と同じとなることがわかる．
相似変換 B = P −1 AP では
B = ( P −1 A) P の1列目は P −1 Aの1列目と同じであるから，
a : Aの第1列目のベクトル
b : Bの第1列目のベクトル
とすると
P −1a =
Pa =
( I − 2ut u)a =
b
*

*
B = 0


0

* 0 
* * 
* * 
  
* * 
0

*
*  となるための
 b の形は

* 
*
 
*
=
b =
0
 

0
 
従って，
勝手な a = t ( a11 , a21 , , an1 ) と上の形の b に対して
0
 
*
t

( I − 2u u)a =
b かつ u =

 
*
となる u があればよい
 b11 
 
 b21 
 0
 
  
 0
 
x 2 = y 2 であるような
任意の n 次元ベクトル x, y に対して単位ベクトル
x−y
u=
x−y 2
u 2 =1
は次を満たす
( I − 2ut u)x =
y
（注）行列式と区別するためベクトルノルムを、 x
2
と表記した。
【証明】
=
xx
t
t
=
xy
x=
2
2
t
y=
yy と
2
2
n
t
x
y
=
∑ i i yx を用いると
i =1
t
x
y
2(
)
(x − y )
−
t
( I − 2u u)x =
(I − t
)x
(x − y )(x − y )
2(x − y )(t xx − t yx)
= x− t
xx −t yx −t xy + t yy
2(t xx −t yx)
x− t
=
(x − y)
t
2( xx − yx)
=y
b = t ( b11 , b21 , 0,  , 0 ) において
2
、
=
b11 a=
b
11
21
a
2
2
n
=∑a
j =1
2
j1
n
2
a
∑ j1 とすれば
j =2
=b
2
11
+b
2
21
= b
2
2
が満たされるので
a−b
u=
a−b 2
とすれば
0
 
*

t
u=
となる。
かつ
−
u
u)a =
b
(
I
2

 
*
a−b
具体的な P (= I − 2u u) , u =
の形は、以下のようになる。
a−b 2
t
まず、 b11 , b21
を前ページのように選ぶと、
 a11   b11 
   
 a21   b21 

w ≡ a −=
b  a31  −  0=
   
     
a   0 
 n1   
 0 


a
s
+
 21 
 a31 


  
 a 
 n1 
2
,=
s 2 b=
21
n
2
a
∑ j1
j =2
s の符号は a21と同符合になるようにとる
(桁落ちを避けるため)
w
2
2
= a−b
2
2
n
= a + 2a21s + s + ∑ a j12 = 2a21s + 2 s 2
2
21
2
j =3
w を用いると P は次のように表せる
P=
I − 2u u =
I −αw w
t
t
1
α= 2
s + a21s
B = P −1 AP は次のようにして効率的に計算できる
B = P −1 AP
(
) (
I − α wt w A I − α wt w
=
)
A − α w t wA − α Aw t w + α 2 w t wAw t w
=
(
A w t q + qt w
=−
)
ここでp 、qはそれぞれ
p = α Aw
q= p −
α
2
w t pw
*
 
*
以上から B = P −1 AP の1列目は b =  0  の形になり
 

0
 
B は対称だから次の形になる
*

*
B = 0


0





* * 
*
*
*

0
*
*

0

*
*


* 
A2 ≡ B と書くことにしてさらに変換することを考える
A3 = P2 −1 A2 P2
A2 の要素を aij(2) とすると
 a11(2)
 (2)
 a21
 0
A2 = 
 0
 

 0

a12(2)
(2)
a22
0
0
(2)
a23
(2)
a24
(2)
a32
(2)
a42

(2)
a33
(2)
a43

(2)
a34

an(2)2
an(2)3

 0 
(2) 
 a2 n 

 a3(2)
n

 
  

(2) 
 a11 
この部分はその
ままにして
ここを0にしたい
そのための変換行列 P2= I − 2u 2t u 2の形は
1

0
−1
0
=
P2 P=
2


0





0 * 
0
1
0

0
0
*

0

0
*


* 
したがって u 2 の形は
0
  第1、第2成分が0
0
u2 =  * 
 

*
 
(3)
A2 , A3 の第2列目のベクトルを a(2)
,
a
2
2 とすると
(2)
(3)
a
−
a
(3)
2
2
=
a(2)
a
かつ
u
2
2 2
2
(3)
2
a(2)
−
a
2
2
が
2
u 2 =t ( 0, 0,*, ,*)
の形になるようにしたい
そのためには
(
t
(3)
(3)
(3)
a(3)
=
a
,
a
,
a
12
22
32 , 0,  , 0
2
)
において
n
(3)
12
a
( )
=
a ,a =
a ,a =
− s2 , s2 =
± ∑ a
(2)
12
とすればよい
(3)
22
(2)
22
(3)
32
j =3
(
(2)
j2
2
(2)
a32
と同符合になるようにする
)
そのとき
w2 ≡ a
w2
2
2
(2)
2
−a
=
a
(3)
2
(2)
2
 0 


 0 
(2)
 a32
+ s2 
=


(2)
 a42 





(2)
 an 2 
−a
(3) 2
3 2
等となる
( )
s =∑ a
j =3
(2)
j2
2
(2)
=
2 s22 + 2 a32
s2
P2 =
I − 2u 2 u 2 =
I − α 2w 2 w 2
t
n
2
2
t
1
α2 = 2
(2)
s2 + a32
s2
変換後の A3 は下のような形になる
A3
*

*
0
−1
P=

2 A2 P2
0


0
* 0 0 
* * 0 
* * * 
0 * * 
   
0 * * 
0

0
*

*


*
【ハウスホルダー変換の図形的な意味】
ハウスホルダー変換は「鏡映」変換
Jacobi,Givensは「回転」
以下同様の変換を繰り返すと
Ak +1 = Pk −1 Ak Pk
Pk =
I − 2u k u k =
I − αk wk wk
t
t
0









0
 (k )

(k )
( k +1)
w k ≡ a k − a k =  ak +1, k + sk 
 ak( k+)2, k 





 a(k ) 
nk


1
αk = 2
sk + ak( k+)1, k sk
s =
2
k
∑ (a )
n
j= k +1
n − 2 回の変換で3重対角行列に変換される
(k )
jk
2
3.6 バイセクション法
（参考程度）
 の固有値を求める
実対称な3重対角行列 A
 α1 β1

0
 β1 α 2 β 2

β 2 α3 β3
1
−

A P=
AP 
=
  


 α n −1
0

β n −1







β n −1 

αn 
固有値は固有方程式
λ I − A
λ − α1 − β1
− β1 λ − α 2 − β 2
− β 2 λ − α3
0
− β3
= 0



− β n − 2 λ − α n −1 − β n −1
0
− β n −1 λ − α n
の解。
ここで行列λ I − A の一連の主小行列式を考える
⇒DKA法の資料(6/25) を参照のこと．
主小行列式
p1 (λ )= λ − α1
λ − α1 − β1
p2 (λ ) =
=(λ − α1 )(λ − α 2 ) − β12
− β1 λ − α 2

pk (λ ) =
λ − α1 − β1
− β1 λ − α 2 − β 2
− β 2 λ − α3 − β3

0



0

λ − α k −1 − β k −1
− β k −1 λ − α k

pn (λ ) =
λ − α1 − β1
− β1 λ − α 2 − β 2
− β 2 λ − α3

0
0
− β3

− β n−2

λ − α n −1 − β n −1
− β n −1 λ − α n
= λI − A
固有方程式は pn (λ ) = 0 となる
pn (λ ), , p1 (λ ) は以下で示すように漸化式
pk (λ ) =
(λ − α k ) pk −1 (λ ) − β k2−1 pk − 2 (λ ), k =
2, , n
を満たす(ただし p0 (λ ) = 1 )。
「スツルム列とスツルムの定理」
スツルム列
実係数を持つ多項式 f ( x)と区間[a, b]が与えられたとする。
このとき次の4条件を満足する実係数多項式の列
f ( x), f1 ( x), f 2 ( x), , f l ( x)
は区間[a, b]においてスツルム列をなすという。
ただし f 0 ( x) = f ( x)
条件 :
(1) 任意の点x において f k ( x)と f k −1 ( x) は同時に0にならない
(2) ある点
=
x0 において f k ( x0 ) 0 ならば f k −1 ( x0 ) f k +1 ( x0 ) < 0
(3) fl ( x) は定符号
(4) ある点
=
x0 において f ( x0 ) 0 ならば f ′( x0 ) f1 ( x0 ) > 0
スツルムの定理
f ( x), f1 ( x), f 2 ( x), , f l ( x)が区間[a, b]で
スツルム列をなし、f (a ) f (b) ≠ 0とする
x を固定して { f ( x), f1 ( x), f 2 ( x), , f l ( x)} を左から右に
見ていったときの符号の変化の回数を N ( x) とすると
N ( x) は x より大きな零点の数となる
したがって区間[a, b]に存在するf ( x) の零点の個数n0は
=
n0 N (a ) − N (b)
pn (λ ), , p1 (λ ), p0 (λ ) は実軸上の閉区間において
スツルム列の条件を満たす
(前の f ( x), f1 ( x), , fl ( x) とは添え字の順序が逆)
区間[a, b]が与えられて pn (a ) ≠ 0 、pn (b) ≠ 0 のとき
λを固定して { pn (λ ), pn −1 (λ ),  , p1 (λ ), p0 (λ )} を左から右に
見ていったときの符号の変化の回数をN (λ )とすると
N (λ ) は λ より大きな固有値の数となる
（参考）
補足
pn (λ ), , p1 (λ ), p0 (λ ) がスツルム列の4つの条件を
満たすことを以下に示す
漸化式を
pk (λ ) =
(λ − α k ) pk −1 (λ ) − β k2−1 pk − 2 (λ )
として
(a )
(1) 隣り合う pk (λ ), pk −1 (λ ) は同時に0にならない
(a ) において=
pk (λ ) p=
0 であるとすると
k −1 (λ )
β k ≠ 0から pk − 2 (λ ) が0となり、以下同様にして
すべての p j (λ ), j ≤ k が0になる。
しかしこれは p0 (λ ) = 1に矛盾する
（参考）
補足
(2) ある点
=
λ0 において pk (λ0 ) 0 ならば pk −1 (λ0 ) pk +1 (λ0 ) < 0
(a ) のk を一つ増やした式
(λ − α k +1 ) pk (λ ) − β k2 pk −1 (λ )
pk +1 (λ ) =
においてpk (λ0 ) = 0 とすると
pk +1 (λ0 ) = − β k2 pk −1 (λ0 )
よってβ k ≠ 0から pk −1 (λ0 ) pk +1 (λ0 ) < 0
(3) p0 (λ ) が定符号であることは p0 (λ ) = 1から明らか
=
λ0 において pn (λ0 ) 0 ならば pn′ (λ0 ) pn −1 (λ0 ) > 0.
(4) ある点
pk (λ ) =
(λ − α k ) pk −1 (λ ) − β k2−1 pk − 2 (λ ) (a)
をλで微分すると
pk′ (=
λ ) pk −1 (λ ) + (λ − α k ) pk′ −1 (λ ) − β k2−1 pk′ − 2 (λ ) (b)
（参考）
補足
(4) (b) × pk −1 − (a) × pk′ −1 を計算すると
pk′ (λ ) pk −1 (λ ) − pk (λ ) pk′ −1 (λ )
= β k2−1 { pk′ −1 (λ ) pk − 2 (λ ) − pk −1 (λ ) pk′ − 2 (λ )} + pk2−1 (λ ).
ここで
qk ≡ pk′ (λ ) pk −1 (λ ) − pk (λ ) pk′ −1 (λ )
とおくと上式は
qk (λ ) =
pk2−1 (λ ) + β k2−1qk −1 (λ ), k =
2,3, , n (c).
ところが
q1 (λ ) =
p1′(λ ) p0 (λ ) − p1 (λ ) p0′ (λ ) =
p1′(λ ) =
1> 0
であるから(c)より
qk (λ ) > 0, k =
2,3, , n
補足
特に k = n のとき
qn (λ ) = pn′ (λ ) pn −1 (λ ) − pn (λ ) pn′ −1 (λ ) > 0
であるが,仮定よりpn (λ0 ) = 0であるから
=
qn (λ0 ) pn′ (λ0 ) pn −1 (λ0 ) > 0
(4)
（参考）
大きいほうからk 番目の固有値λk を
求めるバイセクション法の手順
出発の区間[ a0 , b0 ]を適切な方法によって選ぶ
(1) c = a j −1 + b j −1
j
2
とする
(2) N (c j ) ≥ k ならば a j =
c j , bj =
b j −1
(λk は [c j , b j −1 ]の中にある)
N (c j ) < k ならば a j =
a j −1 , b j =
cj
(λk は [a j −1 , c j ]の中にある)
(3) b j − a j が十分小さくなったら終了。
λk の近似値としてc j を採用する
k 番目の固有値λｋについて
j = 2,3, , に
ついて繰り返す
例バイセクション法途中経過 n = 15、λ7 を求める
[ a j , b j ] j = 3, , 7 の図
出発の区間[ a0 , b0 ]の選び方
固有値の存在範囲について次が成り立つ
x t ( x1 , , xn ) 、n × n 行列 Aに対して
=
n 次元ベクトル
x 1 、 A 1 をそれぞれ次のように定義する
n
Ax 1
x 1 = ∑ xi
A 1 = max
x≠0
x1
i =1
次が成り立つ
n
(1)
A 1 = max ∑ aik
(2)
A 1 ≥ max λi
1≤ k ≤ n
i
i =1
(1)は次のように示される
Ax 1
=
n
n
∑ ∑a x
=i 1 =j 1
ij
n
j
n
≤ ∑∑ aij x j
=i 1 =j 1
n
n

 n


= ∑  x j ∑ aij  ≤ ∑ x j  max ∑ aik 
1≤ k ≤ n
i 1
i 1
=j 1 =
=

 j 1 =

n
n
= x 1 max ∑ aik
1≤ k ≤ n
i =1
したがって
Ax 1
x1
ここで
n
≤ max ∑ aik
1≤ k ≤ n
i =1
右辺がk = mで最大になったとすれば
m番目
x = (0, , 0, 1 , 0,  , 0)とすると
t
=
Ax 1
より
n
a , x
∑=
i =1
im
1
1
n
n
Ax 1
= ∑
=
aim max ∑ aik
1≤ k ≤ n
x1
=
i 1 =i 1
となって等号が成立し、したがって
n
Ax 1
=
= max ∑ aik
A 1 max
1≤ k ≤ n
x≠0
x1
i =1
(2)は次のように示される
λi を A の固有値、xi をその固有値に属する固有ベクトルとすると
Axi = λi xi
から
=
Axi 1
λ=
i xi 1
n
λx
∑=
n
λ ∑
=
x ji λi xi
i ji
i
=j 1 =i 1
したがって
Axi 1
= λi
xi 1
よって
A 1 =max
x≠0
Ax 1
x1
≥ max
i
Axi
xi
1
=max λi
i
1
1
A 1 ≥ max λi であるから
i
出発の区間[ a0 , b0 ]に対しては a0 = − A 1 , b0 = A 1 とすればよい
3重対角行列
 α1 β1

0
 β1 α 2 β 2

β 2 α 3 β3
A=
  


 α n −1
0

β n −1







β n −1 

α n 
ゲルシュゴリン(Gershgorin)の定理
に対しては
n
A 1 =max ∑ aik =max { β k −1 + α k + β k } , β 0 =β n =0
1≤ k ≤ n
と求められる
i =1
1≤ k ≤ n

Download Report