2003/07/08 前回、ヤコビ法とガウズ＝ザイデル法の二つがある。 ここらへんを少しまとめたい.. 連立方程式の解法 直接法 => ガウズの消去法 ( 汎用版 ) が判れば十分 反復法 => ヤコビ法, ガウズ＝ザイデル法, etc.. # 色々あるのは、「決定的な方法がない」ってこと ガウズの消去法 n^3/2 で終了 安心できる方法 スパース行列でも時間がかかってしまう 反復法 収束するかどうかは解らない 上手く行けば得するが、失敗すると悲惨 対角優位でないと収束しない ( 正値 ) スパース行列の利点が利用できる 0 の要素 ( は、そもそも計算が不要なので.. ) そのまま もし、スパース行列が使われないのであれば、ガウズの消去法で、決りだが、 アプリケーション ( 偏微分方程式.. ) では、スパースになる。 # 1000x1000 なら、どっちでもにたようなもの ( Memory に入る ) # 1000000x1000000 なら、反復法しか利用できない( Memory に入らない ) == 残差 : Ax - b のグラフ ( 横軸を反復回数、縦軸を残差の対数を取る ) 性質がよけれあば、 グラフは、だいたい、右さがりの直線的なグラフ 運が悪いと収束しない ( 行列の性質による ) 右さがりの直線的なグラフ => 誤差が 1/2, 1/2 となっている。 x_{n+1} = x_{n} + \delta x_{n} \delta x_{n} = 1/2 \delta x_{n-1} もし、\delta x_{n} の振舞いがいつもこうなら.. 加速できないか ? x_{n+1} = x_{n} + w \delta x_{n} で加速することが考えられる ( SOR 法 )。 # ちなみに、w = 1 の時が、ヤコビ、ガウス＝ザイデル法 A が正則なら 0 < w < 2 で収束することが証明されている。 オストロフスキーの定理 サイエンス社、バーガー著、渋谷 政昭 訳 「大型行列の反復解法」 w をどのような値にすれば、よいかは A による 経験的に、偏微分方程式の場合は w = 1.80 - 1.99 が良い SOR 法 : Successive Over-Relaxation method ( 逐次過剰緩和法 ) ガウス＝ザイデル (G-S) 法では、3 次元の解法では、 a11 x1 + a12 x2 + a13 x3 = b1 a21 x1 + a22 x2 + a23 x3 = b2 a31 x1 + a32 x2 + a33 x3 = b3 で、 x1' = 1/a11 ( b1 - a12 x2 - a13 x3 ) x2' = 1/a22 ( b1 - a21 x1 - a23 x3 ) x3' = 1/a33 ( b1 - a31 x1 - a32 x2 ) とし、 x_1' = x_1 + \delta x_1 x_2' = x_2 + \delta x_2 x_3' = x_3 + \delta x_3 なので、 x_1" = x_1 + w \delta x_1 x_2" = x_2 + w \delta x_2 x_3" = x_3 + w \delta x_3 で、計算する。 S-R 法は、戸川先生の本(絶版)がよい # 本の数字は、活字のものは、信用できない # => プリンタの出力したものそのものの方が安心できる # # 活字は、人間が拾うので、間違うことがある 流体問題は、S-R 法が用いられる。 # ここまでが、シラバス的な内容、以下は、トピックス(試験にでない) 数値計算法の研究の流れ 直接法 ほとんどすんでいる ( イアン・ダラ [ イギリス、スコットランド ] ) 反復法 ヤコビ法 G-S法 SOR 法 共役勾配法 ( CG 法 ) # 新しい方法は、発見当時、流行るが、問題が見付かって # 廃れることがある。 CG 法も、一旦廃れた (収束しないことがあることがわかった) が、最近、 「CG 法 + 前処理」が見直されてきている。 反復法の原理 与えられた方程式 Ax = b を変形して反復方法の式としている ようするに、残差 r = Ax - b を小くするように次の x を决める # r = 0 ならば、答が出たことになる。 r を直接 0 に近付けるのは難しい r^2 を 0 に近付けることを考える。 => CG 法 F(x1,..,xn) = \Sum_{i=1}^n ( b_i - ai1 x1 - .. - ain xn )^2 F(x1,..,xn) >= 0 F(x1,..,xn) = 0 の時、x1, .. xn が解 # これは、問題の焼き直しを行っている # 直線の 0 点 <=> 放物線の最小点 # このように 問題の焼き直しを行うと、色々と便利なことが.. A : 正定値、対称 の時.. r^t A r = (x-x')^tA(x-x') を考える ( r は残差ベクトル ) x^{(1)} : 初期値 r^{(1)} = b - A x^{(1)} P^{(1)} = r^{(1)} # P^{(1)}は修正ベクトル \alpha^{(1)} = \frac{(P^{(1)},r^{(1)})}{(P^{(1)},A P^{(1)})} x^{(k+1)} = x^{(k)} + \alpha^{(1)}P^{(k)} r^{(k+1)} = r^{(k+1)} - \alpha^{(1)}AP^{(k)} \beta^{(k)} = \frac{ (r^{(k+1)},r^{(k+1)}) }{ (r^{(k)},r^{(k)}) } P^{(k+1)} = r^{(k+1)} + \beta^{(k)}P^{(k)} # 行列の計算の回数を減らすための工夫がなされている。 # Ap を一度計算すればよい # ( 二度同じものがでるから、覚えている ) CG 法は、n 回の反復で答がでる !! # もし、無限精度で計算ができたらね.. ( <= 落し穴 ) なぜ、n 回でできる ? x は n 次元上の点 A の直交性を利用して、x の座標を一つずつ、確定する n 回で確定 A の性質によって誤差が累積しなければ、良い方法 P の取りかたは、その n 番目の値を確定するためにきめられる => 解の方向を向いているという保障はない # n - 1 回目までは、あまり残差が減らない可能性がある # => n 回目にいきなり収束するようなことがあっ # ても不思議はない 解の方向を向くように P を定める ( 最急下法 ) ただし、無駄が生じる可能性がある # CG 法では、A の直交性を性を利用する # ので無駄がない CG 法の変形 ( 前処理 ) 一般に.. Ax = b の A が E なら、直解ける A が E に近いほど、解きやすい そこで、 BA x = Bb を考える。 もし、B = A^{-1} なら.. BA = E なので、解けている もし、B = E なら.. BA = A なので、なにもしていない B として、E と A^{-1} の間のものを考えてあげる。 # もちろん、B が逆行列なら解けているわけだが... # B として、なにを考えると良いか ? スパース行列の場合 逆行列は、スパース行列にならない => スパース行列だけど逆行列に近いものを考える # cf. 普通に消去法を適用して、逆行列を解くが.. # もとの行列の要素が 0 なら、逆行列の方も # その要素を強制的に 0 として、計算しない # => B : スパースなので計算がらく # もちろん AB もスパース B の取りかたは色々ある。 # 村田先生 => 対角要素の列を一つ余分に計算するとか.. [まとめ] CG 法は、n 回で答が得られる but.. A が 対称正定値 無限精度 今、Hot な研究者 ファン・デ・アフォルスト BGG - stab 法 == [演習] なし。 [来週(07/15)] まとめ 09:00 - 09:30 試験 09:30 - 10:30 [成績] 7 月 T1 1-2 月 T2 成績 = ((T1+T2)+ΣR)/2 R は 15 回