12.pptx

Transcription

1 数値解析 第 1 回 15 年 7 月 8 日 水 ) 理学部物理学科情報理学コース 1

2 講義内容 1. 非線形方程式の数値解法 1.1 はじめに 1. 分法 1.3 補間法 1.4 ニュートン法 多変数問題への応用 1.4. ニュートン法の収束性. 連立 1 次方程式の解法.1 序論と行列計算の基礎. ガウスの消去法.3 3 重対角行列の場合の解法.4 LU 分解法.5 特異値分解法.6 共役勾配法.7 反復法.7.1 ヤコビ法.7. ガウス ザイデル法.7.3 SOR 法

3 3. 固有値と固有ベクトルの数値計算 3.1 固有値問題の序論 3. ヤコビの方法 対角行列への帰着 ) 3.3 ギブンスの方法 3 重対角行列への帰着 ) 3.4 ハウスホルダー法 3 重対角行列への帰着 ) 3.5 DKA 法 3 重対角行列の固有値の計算 ) 3.6 バイセクション法 3 重対角行列の固有値の計算 ) 4. その他の数値計算法 3

4 成績評価, その他の連絡事項 出席点 5 割 + 期末試験 5 割 講義資料は毎回配布予定です. 欠席した場合, 各自でダウンロードしてください. PDF 版をその週の木曜朝までには公開予定 ) 講義開始後約 3 分で出欠を取ります. 研究室 : watanabe@cc. 九州大学のドメイン ) Subject に, 必ず 数値解析 と記入してください. 4

5 今後の予定 7 月 8 日 水 ) 講義 7 月 15 日 水 ) 休講 7 月 日 水 ) 講義 7 月 9 日 水 ) 講義 8 月 5 日 水 ) 3 限目 13: 14:3) 定期試験 8 月末日成績締切 5

6 対角形に収束するまで繰り返す A 1 A A = P 1 1 A 1 P 1 A 3 = P 1 A P 相似変換の繰り返し 固有値問題の数値解法 3. ヤコビ法 対角行列 対角要素 固有値 3 重対角行列 一旦,3 重対角形に変換する 3.3 ギブンス法 固有値は 3 重対角行列に適した 別の反復法により求める 3.4 ハウスホルダー法 3.5 DKA 法 3.6 バイセクション法 6

7 3.4 Householder 法 Givens 法にならぶ3 重対角化の代表的な方法で, 演算回数が行列の次数 n) と明確に関係付けられる. A: n n対称行列 変換 3 重対角行列 3 重対角化のプロセスは固有値問題以外にも有効に使われ, Householder 変換 とも呼ばれる. 論文 Alston S. Householder, Unitary triangularization of a nonsymmetric matrix, Journal of the ACM JACM), Volume 5, Issue 4 Oct. 1958) Pages:

8 最初のステップ 相似変換 B = B 1 = P APを考える 行または列を一括で になるようにする 参考 Givens 法は要素毎 変換行列に次の形を仮定する P u = I uu t u = t u 1,u,,u n ) n t uu u i i= 1 = = = : n次元ベクトル 1 : 正規化 8

9 P は対称な直交行列 証明 P= I uu) = I uu) = I uu= P t t t t t t t PP = PP = I uu) I uu) t t t t t t t = I uu uu+ 4uuu u t t t t = I uu uu+ 4 u uu) u t t t = I uu uu+ 4u u= I à P= t P = P 1 P= P 1 および A は対称だから B= P AP) = PA P ) = PAP = P AP= B t t 1 t t よって B = 1 P AP も対称 9

10 u として u = 第 1 成分が の形のものを選ぶことができたと仮定すると, このとき P の形は, P = P 1 = P = I uu t 1 1

11 ここで, 任意の行列 G に対して g 11 g 1 g 1n g GP = 1 g g 1 n g n1 g n g nn = GP の 1 列目は G の 1 列目と同じとなることがわかる. g 11 g 1 g n1 相似変換 B= B 1 ) P A P 1 = P APでは 1 P A の 1 列目はの 1 列目と同じであるから, a : Aの第 1 列目のベクトル b : Bの第 1 列目のベクトルとすると 1 P P I a= a= uu t ) a= b 11

12 B = となるための b の形は b = = b 11 b 1 ) と上の形の したがって, a = t a 11,a 1,,a n1 となる t I uua ) = b u があればよい. かつ u = b に対して 1

13 ベクトル u の取り方 x = y 任意の n 次元ベクトル は u = であるような x y x y t I uu) x= y u xy, = 1 に対し, ベクトル をみたす. 注 ) 行列式と区別するためベクトルノルムをと表記. x 13

14 t I uu) x= y の証明 t t xx = x = y = yy と n t t xy = xy i i = yx i= 1 を用いると I u t u)x = I x y) t x y) t x y)x y) )x = x x y)t xx t yx) t xx t yx t xy + t yy = x t xx t yx) t xx t yx) x y ) = y 14

15 b = t b 11,b 1,,,) b 11 = a 11 a u = とすれば u = a b a b において n b 1 = a j1 とすれば j= n = a j1 = b11 + b1 = b j= 1 が満たされるので かつ t I uua ) = b となる. 15

16 t P = I uu), u = a b a b 具体的なの形は, 以下のようになる. まず w a b = b11, b1 を前ページのように選ぶと, a 11 a 1 a 31 a n1 n = = 1 j1 j=, s b a n = = a 1 + a1s+ s + aj1 = a1s+ s j= 3 w a b b 11 b 1 = a 1 + s a 31 a n1 s の符号は a 1 と同符合になるようにとる 桁落ちを避けるため ) 16

17 wを用いるとp は次のように表せる P= I uu t = I αww t α = s 1 + a s 1 B 1 = P APは次のようにして効率的に計算できる 1 B = P AP t ) t I αww A I αww) = = A α A αa + α A t t t t ww ww ww ww t t wq qw) = A + ここで,p, q はそれぞれ p= α Aw α t q= p w pw 17

18 以上から B = 1 P AP の1 列目は b = の形になり B は対称だから次の形になる B = 18

19 次のステップ A B と書くことにしてさらに変換することを考える A A = P A P A = 1 3 ) の要素を a ij ) a 11 ) a 1 ) a 1 ) a ) a 3 ) a 3 ) a 33 ) a 4 ) a 34 ) a n ) a 3n ) ) a 4 a 43 ) a n とすると ) a n3 ) a 11 この部分はそのままにして ここを にしたい 19

20 そのための変換行列 P = P 1 = P = I u t u 1 1 の形は したがって u = u の形は 第 1, 第 成分が

21 A A の第 列目のベクトルを a が, 3 u = t,,,, の形になるようにしたいそのためには 3) = t a 3) 1,a 3),a 3),,, 3 a ) ), a ) 3) ) 3) ) 3) a a = かつ = ) 3) a a a a u において とすると 3) ) 3) ) 3) n ) 1 = 1, =, 3 =, =± j= 3 a a a a a s s a とすればよい ) j ) a ) 3 と同符合になるようにする 1

22 そのとき w a ) a 3) = a ) 3 + s ) a 4 ) a n s n ) = j= 3 a ) j ) 3) ) = 3 = s + a3 s w a a 1 P = I u t t u = I αw w α = ) s + a s 3 等となる

23 変換後の A 3 は下のような形になる A 3 = P 1 A P = Householder 変換の図形的な意味 Householder 変換は 鏡映 変換 Jacobi,Givens は 回転 3

24 以下同様の変換を繰り返すと A P A P = 1 k+ 1 k k k P I I α t t k = uk uk = kwk wk k w k a ) k a k+1) k = k ) a k+1,k k ) a k+,k k ) a nk + s k n ) k = j= k+ 1 n 回の変換で3 重対角行列に変換される α k s = 1 s + a s k ) k k+ 1, k k k a ) jk 4

25 Householder 変換の Fortran コード例 1/3) =========================================================================== THIS SUBROUTINE TRIDIAGONARIZES THE GIVEN SYMMETRIC MATRIX BY HOUSEHOLDERS METHODS A... MATRIX TO BE TRIDIAGONARIZED USED LOWER LEFT TRIANGULAR PART) N... SYSTEM DIMENSION AL.. DIAGONAL ARRAY OF TRIDIAGONARIZED MATRIX BE.. SUBDIAGONAL ARRAY OF TRIDIAGONARIZED MATRIX CO.. NORMALIZATION FACTOR =========================================================================== SUBROUTINE TRIDIAGONARIZATIONA,N,AL,BE,CO) IMPLICIT NONE INTEGER,INTENTIN) :: N REALKIND1D)),DIMENSIONN,N),INTENTINOUT) :: A REALKIND1D)),DIMENSIONN),INTENTOUT) :: AL,CO REALKIND1D)),DIMENSIONN-1),INTENTOUT) :: BE REALKIND1D)),DIMENSION:),ALLOCATABLE :: W,P,Q INTEGER :: I,J,K REALKIND1D)) :: S,T,EPS INTRINSIC SQRT IFN<=) RETURN ALLOCATEWN),PN),QN)) EPS=1.D-5 PARAMETER CHECK ALLOCATE WORK ARRAYS COLUMN THRESHOLD 5

26 DO K=1,N- K-TH HOUSEHOLDER TRANSFOMATION T= DO I=K+1,N T=T+AI,K) END DO S=SQRTT) IFAK+1,K)<.D) S=-S TO AVOID ROUNDING ERROR ALK)=AK,K) BEK)=-S IFT<EPS) CYCLE COK)=1.D/T+AK+1,K)S) WK+1)=AK+1,K)+S DO I=K+,N WI)=AI,K) END DO AK+1,K)=WK+1) DO I=K+1,N T=.D DO J=K+1,I T=T+AI,J)WJ) END DO DO J=I+1,N T=T+AJ,I)WJ) END DO PI)=COK)T END DO Householder 変換の Fortran コード例 /3) 6

27 Householder 変換の Fortran コード例 3/3) T=.D DO I=K+1,N T=T+PI)WI) END DO S=.5DCOK)T DO I=K+1,N QI)=PI)-SWI) END DO DO J=K+1,N DO I=J,N AI,J)=AI,J)-WI)QJ)-QI)WJ) END DO END DO END DO ALN-1)=AN-1,N-1) ALN)=AN,N) BEN-1)=AN,N-1) 前ページのループの終了 DEALLOCATEW,P,Q) END SUBROUTINE TRIDIAGONARIZATION 7