09.pptx

Transcription

1 講義内容 数値解析 第 9 回 5 年 6 月 7 日 水 理学部物理学科情報理学コース. 非線形方程式の数値解法. はじめに. 分法. 補間法.4 ニュートン法.4. 多変数問題への応用.4. ニュートン法の収束性. 連立 次方程式の解法. 序論と行列計算の基礎. ガウスの消去法. 重対角行列の場合の解法項目を変更しました.4 LU 分解法.5 特異値分解法.6 共役勾配法.7 反復法.7. ヤコビ法.7. ガウス ザイデル法.7. SOR 法. 固有値と固有ベクトルの数値計算. 固有値問題の序論. ヤコビの方法 対角行列への帰着. ギブンスの方法 重対角行列への帰着.4 KA 法 重対角行列の固有値の計算.5 ハウスホルダー法 重対角行列への帰着 4. その他の数値計算法 内容は変更の可能性があります 成績評価 その他の連絡事項 出席点 5 割 期末試験 5 割 講義資料は毎回配布予定です. 欠席した場合 各自でダウンロードしてください. hp://yesu.cc.yushu-u.c.p/~we P 版をその週の木曜朝までには公開予定 講義開始後約 分で出欠を取ります. 研究室 : ml: we@cc. 九州大学のドメイン Suec に 必ず 数値解析 と記入してください. 4

2 5.7 反復法 連立 次方程式に対する反復解法 定常反復法 [ 概要 ] 適当な初期値 B c R から初めて反復 : を繰り返す. 解が収束したと見なした場合に反復を停止する. 行列 B とベクトル c を固定し 同じ 定常な 反復操作を適用. B と c の取り方を 種類紹介します..7. Jco 法 A : 行列 A : 対角行列 A : 下 半 三角行列 : 上 半 三角行列... A 対角行列下 半 三角行列上 半 三角行列行列 A を に分離する 計算では配列は用意しない. 6 成分ごとに書くと が解ならば 左辺 右辺という関係 を反復法に適用する 一般的によく使われる技法 7 A が連立 次方程式の解なら左辺と右辺は等しい を成分ごとに表示

3 Jco 方法では 以下のように ステップ目の近似解から sep ステップ目の解を求める方法としてとして利用する. Jco 法の orr プログラム例 反復計算部分 } } } } } do m 反復の開始 m 回まで q 前回の反復で得られたをqにコピー do 次の反復ベクトルを決めるO 文の開始 s.;. 変数の初期化 do - ssa*q 前半の行列とベクトルの積の計算 ed do do A*q 後半の行列とベクトルの積の計算 ed do --s/a 反復ベクトルの 番目の要素を決定 ed do 次の反復ベクトルを決めるO 文の終了 収束判定 ed do 反復終了 9 例. 5 / 7 / 5/ 7 / A 9 / /44 5 / 4/ 4 5/ 7 / / / 5 7 ゼロベクトルを初期値として実際に反復計算してみる. 厳密解は この例では 4 回の反復で誤差は /4~/44 程度 Jco 法の収束条件どのような行列に適用可能か B とおくと Jco 反復列は これらより B B B ここで r r B r とおくと

4 r B r r r Br Br BB : 対称行列 以下の線形代数の知識を証明なしに使います BBr 半 正定値対称行列の固有値は非負実数 r BB TΛ T を満たす直交行列 T が存在する Λ は固有値 λ > を対角要素にもつ対角行列 TT I ユークリッド内積 内積の結合則より 証明は次頁参照 補足 内積の結合則の簡単な証明 B を行列 r を縦ベクトルとしたとき 内積 Br r Br Br r lrl l l rrl BB} l rl l l BBr r 和の順番を入れ替え r ここで BB } l は BB の l 要素の意味. 4 固有値と固有ベクトルを用いて式変形 r r BBr r Λ Λ T Tr r Tr Tr したがって z z z r に注意すると m λ r 次に s Tr なる変換を行うと Λs s m λ s s s r r より m λ r s r 5 よって 収束条件 m λ < のときに限り r となることが解る. 注意 共役勾配法と異なり 反復回数 の 理論的な 上限については 一般の 次元問題について 明示的な上限は導出できない. 6 4

5 前に使った例と同じ行列で 収束条件を確認する 5 5 A 7 7 B 4 BB Guss-Sedel 法 / / Jcoの方法とここが異なる 成分ごとの表示は以下の通り 固有値の絶対値は以下 収束条件を満たす ρ m λ を 収束半径 と呼ぶことがある 7 Guss-Sedel法のアイデア/ Guss-Sedel法のアイデア/ Jco法の成分ごとの反復列 上 順 実 行 # % } } } を最新の値: に置き換える } この部分 の計算において } は既にひとつ前までの行演算で得られている 9 } } } } この置き換えは 繰り返し計算中 最新の近似解のほうが 真の解に少しでも近いであろう という考え期待に基づいている } 5

6 Guss-Sedel 法の収束条件 ここで改めて B とおくと Jco の方法で述べたような収束条件を調べることができる 方法は同じなので省略. より Guss-Sedel 法の orr プログラム例 反復計算部分 do m 反復の開始 m 回まで q 前回の反復で得られたをqにコピー do 次の反復ベクトルを決めるO 文の開始 s.;. 変数の初期化 do - ssa* 前半の行列とベクトルの積の計算 ed do do A*q 後半の行列とベクトルの積の計算 ed do --s/a 反復ベクトルの 番目の要素を決定 ed do 次の反復ベクトルを決めるO 文の終了 収束判定 ed do 反復終了 Jco 法との違いはここだけ 例. 5 7 A 加速過緩和 SOR * 法 Guss-Sedel 法による反復式で * Succesve Over-Relo の略 / / 4 / 4/ 4 または / 4 77 / 576 ゼロベクトルを初期値として実際に反復計算 5 7 この例では 回の反復で誤差は /4~/576 程度 確かに Jco 方法より収束が速いと言えそう ω とおき ω という重みつき平均値をとり 段階で ω : ステップ目の近似値を求める. 緩和係数または緩和因子と呼ぶ 4 6

7 ω ω ω を消去すると このページはあくまで式導出ですので省略可 近似解が Guss-Sedel 法により改良される値 SOR 法の式 I ω ω ω ω となる. 次ページにて導出 を消去すると ω I ω を更に ω 倍して解を更新している 上図は ω> の時のイメージ } 5 ω を左から作用 ω I ω ωi ω } ω ω ω ω ω ω ω à ω ω ω ω I ω } } が等しいので 証明できるかな 6 加速過緩和 SOR 法 SOR 法の orr プログラム例 反復計算部分 ω ω ω: 加速係数または緩和係数 因子 ω では Guss-Sedel 法 ω> で 上方緩和 ω< で 下方緩和 と呼ぶ. w. 加速係数の設定 do m 反復の開始 m 回まで q 前回の反復で得られたをqにコピー do 次の反復ベクトルを決めるO 文の開始 s.;. 変数の初期化 do - ssa* 前半の行列とベクトルの積の計算 ed do do A*q 後半の行列とベクトルの積の計算 ed do --s/a 反復ベクトルの 番目の要素を決定 -w*qw* 反復の加速 ed do 次の反復ベクトルを決めるO 文の終了 収束判定 ed do 反復終了 7 Guss-Sedel 法との違いはここだけ 7

8 定常反復法の一般論 方法 B c ρ m λ 反復行列 Jco 法 Guss- Sedel 法 SOR 法 B 興味のある方は応用として I ω ω I ω } 縮小写像 ニュートン法のとき少しふれた ヤコビ法の収束半径と GS の収束半径 SOR 法の緩和係数 ω の最適値 誤差限界など調べてみてください ω ω c 9 一般に 厳密な ρ m λ により ρ の範囲を評価. 実験 m を求めることは困難な場合が多いので 次の不等式 BB} } ρ BB} m BB} } BB} 次の 重対角行列 Aについて Jco Guss-Sedel SOR 法 それぞれの手法で収束を調べた. 次元数 収束条件 eps<.-6 作図の関係で甘めに とした. A # # 式の ρ は.75 < ρ <.75 であった. すなわち Jco 法で 収束 するはず Logerr SOR ω. Iero 反復回数 Guss-Sedel Jco 実験 次の様な行列 Ar 行列と呼ばれる について Jco Guss-Sedel SOR 法それぞれの手法で収束を調べた. 次元数 収束条件 eps<.-6 とした. A # 収束せず 発散 # # # Jco Iero 反復回数 式のρは 次元 のとき < ρ < 99. であった. - ρ> なので Jco 法の収束は保証されず Guss-Sedel SOR 法は数値例では収束傾向にある. Guss-Sedel : 9 回 SORω. : 6 回 Logerr SOR ω. Guss-Sedel