<4D F736F F F696E74202D C89F090CD8D758B D382E B93C782DD8EE682E890EA97705D>

Size: px

Start display at page:

Download "<4D F736F F F696E74202D C89F090CD8D758B D382E B93C782DD8EE682E890EA97705D>"

あおいひろなが
5 years ago
Views:

1 2005 年 6 月 8 日 ( 水 ) 数値解析講義前期水曜 3 限 3 年 (2005) システム情報学研究院情報理学専攻青柳睦 aoyagi@cc.kyushu-u.ac.jp 研究室 : 情報基盤センター 5 階 502

2 講義内容 { 括弧 } の項目は Skip または参考程度 1. 超越方程式 ( 非線形方程式 ) の解法分法演習 1.2 { 補間法 } 1.3 ニュートン法演習 2. 代数方程式の解法 2.1 グラーフェの方法 ( 演習 ) 2.2 { ベルヌーイの方法 } 3. 連立一次方程式の解法 3.1 行列計算演習 3.2 ガウスの消去法演習重対角行列の場合の解法 3.4 LU 分解法演習講義 vs 演習スケジュールに余裕がある日を見つけて, 数値解析を必要とする理学工学社会学? 分野で活用されている計算科学の概論を紹介する ( 予定 ).

3 3.5 特異値分解法 3.6 共役勾配法演習 3.7 反復法演習ヤコビ法ガウスザイデル法 SOR 法 ( 参考 ) 4. 固有値と固有ベクトルの近似計算 {4.1 レバリエールの方法 } 固有値問題の序 4.2 ヤコビの方法 ( 対角行列への帰着 ) 演習 4.3 ハウスホルダー法 (3 重対角行列への帰着 ) 4.4 バイセクション法 ( 3 重対角行列の固有値と固有ベクトルの計算 ) 4.5 ランチョスの方法 ( 3 重対角行列への帰着 ) 5. ラグランジュの未定乗数法 6. 勾配法 6.1 最急降下法 6.2 最急降下加速法 6.3 ステップ幅の決め方

4 成績評価, その他出席点 5 割 + 前期末試験 5 割講義ノートは, 毎回配りますがもし欠席したら各自でDLしてください server-500.cc.kyushu-u.ac.jp ( 講義ノートはその週の金曜日 17 時までには公開 ) 座席指定でお願いします 30 分後には座席表から出欠を取ります研究室の電話はメールは,aoyagi@cc.kyushu-u.ac.jp Subjectには, 数値解析と記入してください

5 4.1 固有値問題の序固有値 Ax= λx ( A λi) x = 0 A: n n行列 x : 0でない n次元ベクトル ( A λi) det 0 なら A x は0のみ λi の逆行列有 ( A λi) det = 0 : 固有方程式固有多項式 ( n次 )

6 ( A λi) det = 0 : 固有方程式 ( n 次方程式 ) 固有値は n次の代数方程式の解一般に ( n 5) 有限回の代数処理では ( 厳密には ) 求まらない反復法による多数の繰り返し演算が必要 ( 実際には要求される精度に達した時点で打ち切る ) 固有値問題は連立一次方程式同様工学理学社会学の分野で非常に重要 ( 固有振動解析量子力学他 )

7 ( 余談 ) Coupled Simulation ( 連成シミュレーション ) の重要性現実の理学, 工学, 社会学の現象はひとつのモデルで表現できるほど簡単ではない. 現象を理解 ( 説明, 予測 ) するために全体のシステムを部分に分割して考えてみる手法がある. 全体が複雑に絡み合っている様に見える現象の中にも, 部分要素間の相互作用 (Couplings) が本質的であるとみなせる場合がある.

8 ( 余談 ) 連成シミュレーション (1) シミュレーション対象の全体を複数の部分系に分割し部分系の計算途中でお互いのデータを交換しながら全体を解く手法. 部分系 A を解く手法部分系 B を解く手法 ( 例 ) 分子軌道計算 / 分子動力学計算構造力学計算 / 分子動力学計算流体計算 / 構造力学計算統計力学計算 / 分子軌道計算 etc データ交換の頻度とデータ量は部分系間の依存度結合の強さにより弱連成強連成と区別することがある.

9 ( 余談 ) 連成シミュレーション (2) 10~100 サイクル溶媒領域アダティブメッシュ形成 Mediator 離散点登録 r vdw 等のデータ溶質分子属性設定等間隔メッシュ使用メッシュ点数大マルチスケール問題に共通に顕在化する変換技術溶媒分布計算 2 体相関関数有効ホテンシャル又は電荷テータ溶質電子構造計算数 GB RISM FMO 数 GB 共有メモリ SMP 向き, 所要メモリ量 : 数 10GB PC クラスタ向き, 所要メモリ量 : 数 10GB MPICH-G

( 余談 ) 連成シミュレーション (3) RISM FMO 溶媒効果, 水和構造

10 ( 余談 ) 連成シミュレーション (3) RISM FMO 溶媒効果, 水和構造 n 依存 Solvent distribution Solute structure 溶媒緩和溶媒 - 溶質相互作用溶質緩和 Mediator Mediator Grid A In-sphere correlation Grid B MPICH-G2, Globus Fig. RISM-FMO composite simulation RISM: Reference Interaction Site Model FMO: Fragment Molecular Orbital method Fig. Disulfide linkage immersed in H 2 O solvent

11 ( 余談 ) 連成シミュレーションの例 (1/2) 空間, 時間的に局所的な相互作用 (Couplings) でモデル化できる場合に有効 ( 例 ) 2 つの部分系が面で接している構造物とまわりの風 ( 流体 ) など演算量と情報交換量の関係も分散計算向き > 全体は近似モデルを適用するが, 重要な部分は厳密に解く ( 複合シミュレーション ) ( 例 ) 合金の変形, 亀裂生成構造力学 + 分子動力学

12 Tacoma( タコマ ) 橋の崩壊 Washington state, On November 7 of 1940, at approximately 11:00 AM, the Tacoma Narrows suspension-bridge collapsed due to wind-induced vibrations. 橋のねじれ ( 固有 ) 振動 ~ 風による振動

13 単一モデルの適応範囲を超えた? When the twisting motion was at the maximum, elevation of the sidewalk at the right was 28 feet (8.5m) higher than the sidewalk at the left ( It had not been expected). 共鳴共振

14 ( 余談 ) 連成シミュレーションの例 (2/2) 空間時間スケールの異なる現象を ( 多くはむりやり ) 結びつける. Multi-Scale Multi-Physics Simulations 保存量や系の対象性を考慮した物理量の意味変換が必要な場合が多い. ( 例 ) 多粒子の平均運動圧力 ( 例 ) 粒子座標メッシュ ( 格子点 ) ( 例 ) 疎視化動力学

Call join ( フロセス管理 ) あおやぎ研究室では Mediator という通信ライブラリを作っているグリッド連成ミドルウェア (Mediator) は高度セマンティック変換をサポートすることによって指定された相関関係にある離散点上の物理量を交換するグリッド連成ミドルウェア library calls: Call register ( 位置関係, 離散手法, 離散点数, 番号,

15 Call join ( フロセス管理 ) あおやぎ研究室では Mediator という通信ライブラリを作っているグリッド連成ミドルウェア (Mediator) は高度セマンティック変換をサポートすることによって指定された相関関係にある離散点上の物理量を交換するグリッド連成ミドルウェア library calls: Call register ( 位置関係, 離散手法, 離散点数, 番号, 座標, 要素情報, 属性 ) Sim. A Mediator A Mediator B Sim. Sim. B B Sim. A Sim. A Sim. A Mediator A Mediator A Mediator A 離散点の位置関係を探索しデータを自動的に通信変換 Middleware system Mediator B Mediator B Sim. B Sim. B Call send(byte, 物理量 a) Call receive(byte, 物理量 b, 変換関数 ) 最近接球内相関第 1 近接矩形内相関

16 本題に戻ります相似変換 B= 1 P AP P : n n 正則 ( 逆行列のある ) 行列固有多項式は変わらない固有値も変わらない det ( ) ( 1 B λi = det P AP λi) ( 1 1 P AP λp P) = det ( λ ) ( 1 P A I P) = det ( A λ ) ( A λi) ( 1 I PP ) = det = det det det ( 公式 ) 行列積の行列式 ( XY) = det ( YX) ( XYZ ) = det X ( YZ ) = det ( YZ ) X = det ( YZX ) ( ) ( )

17 A A A A = P AP = P A P 相似変換の繰り返し収束するまで繰り返し * * 0 * 0 * * 対角行列ヤコビ法等対角要素固有値有限回の変換 * * * * * 0 * * 0 * * * * 3 重対角行列ハウスホルダー法等固有値別の反復法により求める

18 4.2 ヤコビの方法相似変換の行列 P 1 0 cosφ sin φ P = sinφ cos φ 0 1 p q p q t PP P 1 t = = I P 直交行列変換後の行列 B = 1 P AP の特定の成分を 0 にする

19 相似変換 B 1 = P AP において 1) A が対称 t A = A ならば B = ( P AP) = PA ( P ) = P AP = B t t 1 t t 1 1 よって B も対称 2) A の pq, 行目および pq, 列目のみが変更される

20 1 P A = 1 a11 a1 p a1 q a1 n 0 cosφ sinφ ap 1 app apq apn sinφ cosφ aq 1 aqp aqq aqn 0 1 a a a a n1 np nq nn p q p q p q = a a a a a a a a a a a a a a a a 11 1p 1q 1n p1 pp pq pn q1 qp qq qn n1 np nq nn p q p q a = a cosφ a a = a sinφ + a pk pk qk qk pk qk sinφ cosφ

21 B = = 1 P AP a11 a1 p a1 q a1 n 1 0 a p1 a pp a pq a pn cosφ sinφ a q1 a qp a qq a qn sinφ cosφ 0 an 1 anp anq a nn 1 p q p q = a11 a 1p a 1q a1 n a p1 a pp a pq a pn a q1 a qp a qq a qn an 1 a np a nq a nn p q p q a = a cosφ a a = a sinφ+ a kp kp kq kq kp kq sinφ cosφ k p, q 2 2 a pp = app cos φ + aqq sin φ ( apq + aqp )cosφsinφ 2 2 a qq = app sin φ + aqq cos φ+ ( apq + aqp )cosφsinφ 2 2 a pq = ( app aqq )cosφ sinφ+ apq cos φ aqp sin φ 2 2 a qp = ( app aqq )cosφ sinφ+ apq cos φ aqp sin φ

22 a ij = a ij を用いると結局 bij = aij i, j pq, bpk = bkp = apk cosφ aqk sin φ k p, q bqk = bkq = apk sinφ+ aqk cos φ k p, q app + aqq app aqq bpp = + cos 2φ apq sin 2φ 2 2 app + aqq app aqq bqq = cos 2φ+ apq sin 2φ 2 2 app a qq bpq = sin 2φ+ apq cos 2φ 2 b pq = 2apq 0 とするためには tan 2φ = とすればよい. a a pp qq

23 ( 注 ) もし a pp = a qq ならば ϕ = sign( a pq π ) 4 と選ぶ pq, の選び方絶対値が最大の非対角要素が 0 となるように選ぶ一度 0 になった要素が後の変換で 0 でなくなることもありうる ( 次に収束の保証について考察する )

24 ヤコビ法の収束について (a) 相似変換 B a = b 2 2 ij ij i, j i, j 1 = P AP において次が成り立つ ( 全要素の 2 乗の和は変わらない ) P は直交行列 BB = ( P AP)( P AP) = PA ( P ) P AP = P AAP t t 1 1 t t t 対角和は相似変換によって不変だから tr T A A = tr T B B 成分ごとに書くと 2 2 a = b ij i, j i, j 1 t ( P = P) ij だから

25 補足対角和 ( トレース ) tr X n = i= 1 x ii XY の対角和 = YX の対角和 n n n n n tr XY = ( XY ) = xij y ji y jixij tryx ii = = i= 1 i= 1 j= 1 j= 1 i= 1 相似変換によって変わらない ( ) ( ) ( ) Y XY = Y XY = XY Y = X YY = X tr tr tr tr tr

26 (b) b + b = a + 2a + a pp qq pp pq qq の証明. まず b + 2b + b = a + 2a + a pp pq qq pp pq qq を示す app + aqq app aqq bpp = + cos 2φ apq sin 2φ 2 2 app + aqq app aqq bqq = cos 2φ + apq sin 2φ 2 2 app aqq bpq = sin 2φ+ apq cos 2φ 2 app + aqq app aqq において F =, G =, H = a 2 2 bpp = F + Gcos 2φ Hsin 2φ bqq = F Gcos 2φ + Hsin 2φ bpq = Gsin 2φ+ Hcos 2φ pq とおくと

27 ( ) ( ) x + y + x y = 2x + 2y pp pq qq { ( φ φ) } ( φ φ) ( G φ H φ) 2 ( ) 2 2 b + 2b + b = F + Gcos 2 Hsin 2 + F Gcos 2 Hsin sin2 + cos2 ( φ φ) ( φ φ) F 2 Gcos2 Hsin2 2 Gsin2 Hcos2 = F { G φ H φ GH φ φ} { G φ H φ GH φ φ} ( F G H ) = cos 2 + sin 2 2 cos2 sin2 + 2 sin 2 + cos cos2 sin2 = app + aqq a = app + aqq = 2 + a 2 = a + 2a + a pp pq qq 2 pq pp a 2 qq 2 + a 2 pq

28 b + 2b + b = a + 2a + a pp pq qq pp pq qq において P は b = 0 となるように選んだから pq b + b = a + 2a + a pp qq pp pq qq 以上 (a) (b) から B の非対角要素の 2 乗の和は b = b b + b + b ij ij ii pp qq i j i, j i p, q となる = a a + a + 2a + a ij ii pp pq qq i, j i p, q = a 2a i j 2 2 ij pq

29 b = a 2a ij ij pq i j i j において全てのi a a pq, 2 2 ij pq 和をとると i j よって a = max a は 2 2 pq ij i j jについての項数は i j a n( n 1) a 2 2 ij pq 2 n となるように選んだから nだから 2 nn ( 1) i j b = a 2a 1 2 a ij ij pq ij i j i j nn ( 1) i j a 2a 2 2 ij pq

30 b ij ij i j nn ( 1) i j 2 k 1回相似変換を行った後の行列 Ak の要素を a ( k ) ( k a ) 2 ij 1 ( ) 2 ij i j nn ( 1) i j 2 a k n 3 に対しては 0< 1 2 nn ( 1) < 1 だから k のとき ( ( k ) a ) 2 i j ij 0 a ij と書くとよって ( 上の P による ) 相似変換の繰り返しで非対角成分は全体として 0 に収束する

31 固有値と固有ベクトル k 回目の相似変換の行列を P k と書く 1 A = k+ 1 Pk AkPk, k = 1,2, T = lim PP P k 1 2 を作ると以上の結果から 1 Λ = T AT ここで Λは対角行列 λ1 λ 0 2 Λ= 0 λ n k λ i : Aの固有値

32 1 Λ= T AT から AT = TΛ T t i の第 i列からなる列ベクトルをt が T t1 t2 t n At = (,,, ) k = λ t k k Aの固有値 λ i i とするとに属する固有ベクトルとなる

09.pptx

09.pptx 講義内容数値解析第 9 回 5 年 6 月 7 日水理学部物理学科情報理学コース. 非線形方程式の数値解法. はじめに. 分法. 補間法.4 ニュートン法.4. 多変数問題への応用.4. ニュートン法の収束性. 連立次方程式の解法. 序論と行列計算の基礎. ガウスの消去法. 重対角行列の場合の解法項目を変更しました.4 LU 分解法.5 特異値分解法.6 共役勾配法.7 反復法.7. ヤコビ法.7.