数値解析 第 1 回 15 年 7 月 8 日 水 ) 理学部物理学科情報理学コース 1
講義内容 1. 非線形方程式の数値解法 1.1 はじめに 1. 分法 1.3 補間法 1.4 ニュートン法 1.4.1 多変数問題への応用 1.4. ニュートン法の収束性. 連立 1 次方程式の解法.1 序論と行列計算の基礎. ガウスの消去法.3 3 重対角行列の場合の解法.4 LU 分解法.5 特異値分解法.6 共役勾配法.7 反復法.7.1 ヤコビ法.7. ガウス ザイデル法.7.3 SOR 法
3. 固有値と固有ベクトルの数値計算 3.1 固有値問題の序論 3. ヤコビの方法 対角行列への帰着 ) 3.3 ギブンスの方法 3 重対角行列への帰着 ) 3.4 ハウスホルダー法 3 重対角行列への帰着 ) 3.5 DKA 法 3 重対角行列の固有値の計算 ) 3.6 バイセクション法 3 重対角行列の固有値の計算 ) 4. その他の数値計算法 3
成績評価, その他の連絡事項 出席点 5 割 + 期末試験 5 割 講義資料は毎回配布予定です. 欠席した場合, 各自でダウンロードしてください. http://yebisu.cc.kyushu-u.ac.jp/~watanabe PDF 版をその週の木曜朝までには公開予定 ) 講義開始後約 3 分で出欠を取ります. 研究室 : 9-64-95 E-mail: watanabe@cc. 九州大学のドメイン ) Subject に, 必ず 数値解析 と記入してください. 4
今後の予定 7 月 8 日 水 ) 講義 7 月 15 日 水 ) 休講 7 月 日 水 ) 講義 7 月 9 日 水 ) 講義 8 月 5 日 水 ) 3 限目 13: 14:3) 定期試験 8 月末日成績締切 5
対角形に収束するまで繰り返す 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
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: 339-34 7
最初のステップ 相似変換 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
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 1 1 1 よって B = 1 P AP も対称 9
u として u = 第 1 成分が の形のものを選ぶことができたと仮定すると, このとき P の形は, P = P 1 = P = I uu t 1 1
ここで, 任意の行列 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
B = となるための b の形は b = = b 11 b 1 ) と上の形の したがって, a = t a 11,a 1,,a n1 となる t I uua ) = b u があればよい. かつ u = b に対して 1
ベクトル u の取り方 x = y 任意の n 次元ベクトル は u = であるような x y x y t I uu) x= y u xy, = 1 に対し, ベクトル をみたす. 注 ) 行列式と区別するためベクトルノルムをと表記. x 13
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
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
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
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
以上から B = 1 P AP の1 列目は b = の形になり B は対称だから次の形になる B = 18
次のステップ 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
そのための変換行列 P = P 1 = P = I u t u 1 1 の形は したがって u = u の形は 第 1, 第 成分が
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
そのとき 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 等となる
変換後の A 3 は下のような形になる A 3 = P 1 A P = Householder 変換の図形的な意味 Householder 変換は 鏡映 変換 Jacobi,Givens は 回転 3
以下同様の変換を繰り返すと 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
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
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
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