2005/10/18 数値解析学と演習 福井 先生 ヤコビ法は、対称行列を想定している。 累乗法は、(対称行列にも適用可能だが..) 非対称行列を対象とした( 最も簡単な.. ) 方法 [累乗法] 与えられた行列 A に対して、適当な x_0 ( \ne 0 : 「適当な」というのはようするに、何が最適かが解らないということ.. ) を取り、 x_{k+1} = A x_{k} を繰り返し計算する。 [累乗法の数学的な原理] n x n 行列 A の (相異なる..) n 個の固有値を \lamda_i ( i = 0 .. n-1 ) とし、 i < j なら \lambda_i > \lamnda_j とする。 また、 \lambda_i に対応した固有ベクトルを v_i とする。 固有ベクトルの性質により、 x_0 = c_0 v_0 + c_1 v_1 + .. + c_{n-1} v_{n-1} となるので、 x_1 = A x_0 = A (c_0 v_0 + c_1 v_1 + .. + c_{n-1} v_{n-1}) = c_0 A v_0 + c_1 A v_1 + .. + c_{n-1} A v_{n-1} = c_0 \lamnda_0 v_0 + c_1 \lambda_1 v_1 + .. + c_{n-1} \lambda_{n-1} v_{n-1} となる。 一般に、 x_k = A^k (c_0 v_0 + c_1 v_1 + .. + c_{n-1} v_{n-1}) = c_0 A^k v_0 + c_1 A^k v_1 + .. + c_{n-1} A^k v_{n-1} = c_0 \lambda_0^k v_0 + c_1 \lambda_0^k v_1 + .. + c_{n-1} \lambda_0^k v_{n-1} = \lambda_0^k ( c_0 v_0 + c_1 (\frac{\lambda_1}{\lambda_0})^k v_1 c_1 (\frac{\lambda_1}{\lambda_0})^k v_1 + .. + c_{k} (\frac{\lambda_{k}}{\lambda_0})^k v_{k} c_{k} (\frac{\lambda_{k}}{\lambda_0})^k v_{k} ) となる。 \lambda_0 > \lamnda_i より、 \frac{\lambda_i}{\lamnda_0} < 1 よって、 \limit{k\rightarrow\infty} \frac{\lambda_i}{\lamnda_0}^k \rightarrow 0 つまり、 x_{k+1} = (\lambda_0)^k x_k となるので、 \lambda_0 = \frac{|x_{k+1}|}{|x_{k}}| で、固有値 \lambda_0 ( 最大の固有値 ) が得られることになる。 # こんなに、簡単に行くの ? .. 「巧く行くときには、巧く行く」 \lambda_0 と \lambda_{k} の比が 1 に近いとだめ # 倍違えば、だいたい大丈夫 !! 現実の問題 ( 例 : 建物の固有振動数. ) 等は、( 対称性を取り除いておけば.. ) 近い値にはならない。 弦の振動の場合 \lambda_1 は、\lambda_0 の半分になることが解っている 利点 行列とベクトルの掛け算だけで Okey => スパースの場合などは、更に、計算がサボれる 欠点 最大の固有値しか計算できない => 「減次 ( text p.231 )」を行う。 [減次] \lambda_0 が求まれば、 v_0 が得られるので、行列 A から、\lambda_0 成分を取り除いた A' を作る ( text p.338 ) A v_0 = \lambda_0 v_0 A^t w_0 = \lambda_0 w_0 s = \frac{(w,v)} として、 A' = A - svw^tA x_0 = x - \frac{(w,x_0)}{(w,v)}v とすればよい。 A' の固有値は、 \lambda_1, \lambda_2, .., \lambda_{n-1}, 0 となる。実際には、誤差のために \lambda_1, \lambda_2, .., \lambda_{n-1}, \eps なので、 \eps が、0 でないのでこれが成長するとだめ。 # 誤差が生れないように、固有ベクトルの直交性を満すように注意深く計算を行うと、確かに巧く行くのだが、このように直交性を満す方法を取ると結果的に、計算量が、膨大になり、ヤコビ法の方が、速くなってしまう !! ## 直交性を満すようにすると、元々の行列がスパースでも、途中の計算で、デンスな行列の計算が必要になり、羃乗法のメリットが失われてしまう !!! # 無限桁 (誤差なし..)ならば、これで Okey # 現実の世界では、誤差があり、これが、成長することになるので、最初の数個以外は、実用にならない。 ## 最初の数個でも意味はある ## cf. 地震の例 : 固有値が大きいということは、周波数が低いもの。建物が壊れるかどうかどうかは、周波数が低いものが影響するので、これで十分 ### 先日の地震は、関東の地下の構造が複雑なので、面白い現象がおきている ### 遠い方が近い所より、震度がおおきい ### => 地下の岩盤の形状が複雑なため ### ロサンジェルスの地震 ### 新しい建物 ( 高い建物 ) だけが壊れた ### 低い周波数が共振した ### 関東地震の数値シミュレーション結果 ### 高い振動数の波は、減数が激しい ### 低い周波数が遠くまで伝達する !! ### 高い建物だけが壊れる可能性が高い 逆に、絶対値の小さい固有値 ( 高周波成分 ) が知りたい !! どうするか ? (数学的には..) A^{-1} の固有値は、\frac{1}{\lamda_0}, \frac{1}{\lamda_1} .. となる A^{-1} で、反復法を行えば、\frac{1}{\lamda_{n-1}} が得られる => 逆反復法 実際には、A^{-1} が、数値的に不安定なので、 x_{k+1} = A^{-1}x_k から、 A x_{k+1} = x_k を導き、この連立方程式を ( LU 分解等を用いて.. ) 解く # 実は、この為に「LU 分解」の話をした この方法は、収束が速い.. ( たった、三度で収束してしまう !! ) # 有名な先生の論文のプログラムに loop がないので変だと思ったが、確かめてみたら、確かに、三回で収束し、loop が不要だった !! [原点移動] (text p.238) A' = (A-\sigma I)^{-1} とすれば、固有値は、 \frac{1}{\lambda_0-\sigam}, \frac{1}{\lambda_1-\sigam},.., \frac{1}{\lambda_{n-1}-\sigam} となる。 この A' で、反復法を行えば、\sigam に近い \lambda_i に対して、\frac{1}{\lambda_i-\sigam} が最大になるので、これが得られる。 \sigam に近い \lambda_i が欲しい時に得られる # どのような固有値が欲しいかによって、手法を選択する !! # 全部必要 => ヤコビ # 最大/最小が必要 => 羃乗法 # ある値に近い固有値 => 原点移動 [羃乗法を使う時の注意] 同じ固有値があるとだめ !! 問題そのものに対称性があると、自然に重複固有値が出る => これを「取り除く」という作業をしないと駄目 [例] 正方形の場合は、対称性が 3 軸くらいある 軸の所で、折り曲げて、対称性を消去した問題を考える ## ヤコビ法の場合は、重複していても問題ない == 前回の演習に関する注意 text の p.276 のプログラムに誤りがある 1) kmanx の値が初期化されていない kmax = 2 * n * n 位が妥当 ( もちろん、問題による ) 2) 下から 7 行目に意味のない部分がある if ( fabs(apg=agp)