2005/07/05 数値解析学と演習 福井 先生 前回は、LU 分解を学んだ 密行列の場合の解法は、ガウスの消去法 ( = LU 分解 ) で決り !! (しかし..)現実的に解かなければならない連立方程式 編微分方程式の差分解法から出てくる係数行列の連立方程式を解く => 一般に疎行列になる。 [例] 規則正い図形(矩形)の偏微分方程式 1 次元 : 前後のみ 0 でない => 3 重対角行列 # まあ、これは、ガウスでも +--O--*--O--+ 2 次元 : 前後 + 上下 => 3 重対角行列 + 二つの対角 | | | | | +--+--O--+--+ | | | | | +--O--*--O--+ | | | | | +--+--O--+--+ | | | | | 3 次元 : 前後 + 左右 + 上下 => 3 重対角行列 + 四つの対角 O O |/ O-*-O /| O O 疎行列で、ガウスの消去法を使うと、fill-in ( 0 の所が 0 でなくなる.. ) という現象が生じる => 折角、 0 になっている ( = 計算しなくてよい部分 ) が、また 0 でなくなる ( = 計算が必要な部分になる ) のは望ましくない n = 10^8 ( 現状で実際に計算される次元数 ) だと、そもそもメモりにも入らない.. # 安定した現象から作られた係数行列は、対角要素の絶対値が大きい # => 安定した現象では、「自分自身から影響 = 対角要素」が一番大きい !! 疎行列の場合は、疎行列の性質 ( 0 の要素が多い ) を利用して、計算量を減らすことのできる、反復方法が望ましい !! [ヤコビ法] a_11 x_1 + a_12 x_2 + a_13 x_3 = b_1 a_21 x_1 + a_22 x_2 + a_23 x_3 = b_2 a_31 x_1 + a_32 x_2 + a_33 x_3 = b_3 に対して ( 対角要素が優位で、答があれば.. )、 x_1 = \frac{1}{a_11} ( b_1 - a_12 x_2 - a_13 x_3 ) x_2 = \frac{1}{a_22} ( b_2 - a_21 x_1 - a_23 x_3 ) x_3 = \frac{1}{a_33} ( b_2 - a_31 x_1 - a_32 x_2 ) の計算を繰り返し、収束すれば Okey => これがヤコビ法 # 固有値の解法でも同じ名前の解法があるが、それとは別法 !! [収束判定の方法] 1) 解で見る 収束しているかどうかで確認 => cf) \Sum_{i=1}^n |x_{k+1}-x_k| < ε # これは、解が同じ程度の大きさななら Okey # 解の間で、絶対値に差があると、結果が、絶対値の大きいものに引ずられる。 # # 絶対値の大きい解が収束しても、小さい方の解は、まだ、収束していなかったりする 2) Ax=b に代入して成立するかどうか 本来の意味からいえばこちらが正い判定方法 ただ、計算量が増えるので、できれば避けたい。 # 収束判定は状況に応じて、選択する # 絶対値が同じ位だったら、解で見る、そうでなければ方程式に代入.. # # 演習問題だったら、固定回数もあり.. ? # このプログラムを素直に書くと、実は、ヤコビ法にならなない。 # # x_1 = \frac{1}{a_11} ( b_1 - a_12 x_2 - a_13 x_3 ) # x_2 = \frac{1}{a_22} ( b_2 - a_21 x_1 - a_23 x_3 ) # ^^^ ここで新しい値を使ってしまう !! # ヤコビ方法にならない ## 確かに、これは「ヤコビ法」としては..間違いなのだが.. ## しかし、どうせ、繰り返すのであれば、新しい値の方が良いのでは ? ## この「誤り」を肯定して取り入れると... ガウス=ザイデル法 ### どうしてこんなことを考えたのか.. ? きっと最初は間違えたんだろうなぁ.. [ガウス=ザイデル法] 普通は、ヤコビ法より、ガウス=ザイデル法の方が収束が速い # せいぜい 2 倍位 ## しかも、プログラムも簡単 ## => ヤコビ法では、古い値を保存しておく必要がある # 一般に収束までの繰り返し回数は、係数行列の性質によって異なる # 最初から対角行列なら、一度で御仕舞い # 全ての要素が同じような大きさだと、なかなか収束しない。 # ヤコビ法の方がガウス=ザイデル法が速い行列も作れる # 昔だったら、これで話が御仕舞いだったが.. 最近は、並列処理がされるようになった ヤコビは、並列化が簡単 => 各々の列での計算が独立 ガウス=ザイデル法は、最新の要素を利用する => 次の行で前の行の値を使う => 並列化が困難 並列処理をするならば、ヤコビ法の方が有利 # 行列の性質を調べて、都合が良い方法を選ぶ [収束条件] |a_ii| > \sum_{j=1..n, i \ne j} a_ij 一般に、偏微分方程式であれば、 2 \pm h ; h は刻み幅 => 一般は成立する # 処が.. 刻み幅を小くすると.. この条件が成立しにくくなくなる # メッシュを小くすれば、Δx も小さくなるので、方程式の係数が似たような値になり、方程式間の関係が強くなるので、数値計算的に不安定になる [プログラミング上の注意] 1 次元であれば、前後で 2 つ 2 次元であれば、前後左右で 4 つ 3 次元であれば、前後左右上下で 6 つ の要素だけを計算すればよい n が大きくなると、これは大変な、高速化になっていることが解る。 # 19 世紀(計算機ができる前)にできた方法だが、未だに現役で頑張っている # 新しいのは、QR 法位 # => 他にも色々新しい方法が生き残っていない.. == # 少し、時間に余裕があるので、反復法の話 # cf. 戸川先生の「マトリックス数値計算法」 連立方程式の一次方程式の反復法 ヤコビ法 ガウス=ザイデル法 ( ここまでは、本日、学んだ ) SOR 法 ( 来週やる : ガウス=ザイデル法の変形 ? ) CG 法 ( 30 年前位に考えられた ) # 一般に、反復方法は、繰り返す回数が事前には解らない。 # 速いものは、直に、遲いものは、なかなか収束しないのが普通 => 原理的に n 回で収束することが解っている 行列の直交性を利用することによって、n 回の反復で計算 # その意味で、厳密な意味では、反復法とは言えないが、反復法的な利用が可能 暫く、もてはやされたが、直に誰も利用しなくなった.. why ? => 収束しない !! # 確かに、理論的には、無限桁で計算すれば n 回で、できる 誤差が入ると、直交性が崩れてしまう !! # アスペクト比 ( 縦横方向の比率 ) が、大きいと誤差が大きくなりすぎる # => 一般には、そうなることが多い 最近 ( 15 年位 ) 復権 プレコンディショニング # 「アスペクト比」を改善することにより、誤差の発生を防ぐ => 問題によって、処理を変える => 色々な論文がある # ということは、ようするに、「決定打」がないということ !! # オランダのファンデア=フォルストの成果が一番使える CG の収束性 # 他の反復法は、回数によって、残差が減るのが普通 n 回目で突然、残差がググッと減るようなパターンも在りうる !! フォルストの方法 CG 法 + 反復の度に、残差が減るような工夫がされている CG 法の対象 対称行列を考える 非対称行列の場合は BCG 法と呼ぶ # フォルストの方法は BCG の改良 == [まとめ] 反復法のメリット 0 の係数が沢山ある ( 疎行列の ) 性質を上手く利用できる手法 => 演算量と、記憶量を削減することができる # 次元の多い問題では、反復法しか選択できない 東大の吉村先生 ( アドベンチャープロジェクト ) 一億次元 # 元々、原子炉の研究者 # 付属物もついた原子炉をまるまるシミューレション # => 一億元になる.. ジャンボジェット機の最初の設計の時の解析 21 万元 ( 2 x 10^5 ) 程度 == [演習] # 去年の演習で利用したテキストの問題 A = ( 5 -4 -2 1 ) b = ( 5 ) 2 3 1 2 21 -2 2 8 2 16 4 1 -3 5 18 # 対角が優位にしてある点がポイント # 一部、成立していない行もあるが、この程度なら大丈夫 .. == [参考] プログラムをどう作るか ? 結局は、 i=1,n x_i = \frac{1}{a_ii} ( b_i - \sum_{j=1}^{i-1} a_ij x_j - \sum_{j=i+1}^n a_ij x_j ) = \frac{1}{a_ii} ( b_i - \sum_{j=1,i \ne j}^{n} a_ij x_j ) の計算がポイントなのだが.. # 対応方法として、 # Loop を分ける => プログラムが複雑 # if 文を入れる => 遅くなるので嬉しくない # 次の変形は、「旨い方法」に見えるが.. . = \frac{1}{a_ii} ( b_i - \sum_{j=1}^{n} a_ij x_j ) - x_i ととすると、確かに 数学的 ( = Coding 的 ? ) には、正いが.. 最後の引き算の所で、対角要素 ( = 絶対値の大きな数 ) で引くので、誤差が大きくなる。 # 数値計算的には、誤差が大きくなるのが大変気になる。