今回から、しばらく線型方程式の解法 ( 前期一杯 ) の話になる。 解く対象 A x = b A \in M[n,n], x, b \in V[n] 一般 ( 理論的 ) には、公式 (クラーメル) があり x_i = \frac{det(A_i)}{det(A)} where A_i は A の i 列目を b に書き換えた者 で求めることができる。 でも、実際の計算では利用しない Why det(A) の計算に n^3 時間かかるから.. ( 演算量が大きすぎる ) # 理論的には正しい # 解の存在 ( det(A) = 0 ? ) # 解の性質 # を考える場合には、利用する。 じゃあ、どう解くか ? 直接法 ガウスの消去法 # 計算機が出てくる前から利用されていた.. 反復法 SOR etc.. # まだ、発展途上.. ガウスの消去法 直接法としては、これで決り ( 特殊な状況がなければ ) # 他の方法は、ガウスの消去法の変形に過ぎない == 問題 ( n = 4 の場合 ) A x = b a_11 a_12 a_13 a_14 x_1 b_1 a_21 a_22 a_23 a_24 x_2 b_2 a_31 a_32 a_33 a_34 x_3 b_3 a_41 a_42 a_43 a_44 x_4 b_4 が与えられたらどうするか ? # つまり、「解ける」というのはどうゆうことか ? a_11' x_1 b_1' a_22' x_2 b_2' a_33' x_3 b_3' a_44' x_4 b_4' なら、解けたと考えてよい。 要するに、 a_11' x_1 b_1' a_22' x_2 b_2' a_33' x_3 b_3' a_44' x_4 b_4' なので、直に、計算できる。 じゃあ、どうやって、この形にするか ? 簡単だが、間違った方法 単に、a_ij ( i != j ) を、無条件に 0 にしてしまう 形式上では解けたと同じ でも、誤っている.. Why ? 方程式の意味を変ている 方程式の意味を変えずに、a_ij を 0 にできれば良い 意味を変更しない、操作をして、都合の良い形 ( a_ij = 0 ) にしたい 基本変形を利用する # 行列の変換 / 行列の操作 基本変形 本の行列 A で、在る行のα倍を他の行に加えた結果えられる行列 A' を考えると Ax = b の根と A'x = b' の根は同じ 例: 一行目をα倍して二行目に加える (α a_11 + a_21)x_1 + (α a_12 + a_22)x_2 + (α a_13 + a_23)x_3 + (α a_14 + a_24)x_4 = α b_1 + b_2 a_21' を 0 にしたければ.. a_21' = (α a_11 + a_21) なので、 α a_11 + a_21 = 0 つまり、 α = - \frac{a_21}{a_11} とすれば良い。 他の a_ij' も、同様にして 0 にできる 順番を気を付ける必要がある # 左上から下の方向に消せば Okey この操作で、一旦 0 にした要素が、また、0 以外になったら困る ガウスの消去法に関しては ( 間違わない限り.. ) そうはならない # この性質は、「連立方程式を解く」場合のとっても嬉しい性質 # たの場合 ( 例えば、後期に行う固有値の場合は、成立しない .. ) # もし、一旦 0 になった項が再び、計算の途中で 0 じゃなくなる場合は、有限回数で # 終るという保証ができなくなる可能性がある ( 無限 Loop になるかも !! )。 == 基本的な操作は、「基本操作 : ある行を定数倍して、別の行に加える」だけ !! # これを繰り返すだけで、希の形にできる ## 計算量を更に減らすために、幾つか工夫がある。 このような形で、a_ij を 0 にして行く方法 消去法 ( ガウス = ジョルダン 法 ) == 消去法 ( ガウス = ジョルダン 法 ) の計算量 一つの a_ij を消すには.. 最初は、n 回のかけ算と足し算 次は、n - 1 回 結局、大体 \frac{2}{3}n^3 となる。 # クラーメルは n^3 だったことに注意 !! # もっと、計算量を減らすには ? とりあえず、下半分を 0 にしてから.. \Sum{i=1}{n-1} n^2 x_n から順に代入法で解く \frac{1}{2} n^2 で解ける 下半分を 0 にする (前進消去) 計算量 Do K = 2, N do I = K+1, N do J = K, N A(I,J) = A(I,J) - (A(I,K)/A(K,K))*A(K,J) end do end do end do # これは、三角垂の体積と同じオーダの量 # \frac{1}{3} N^3 となる x_n から順に代入法で解く(後退代入)の計算量 三角形の面積と同じ \frac{1}{2}N^2 ガウス法 ( = 先進消去 + 後退代入 ) 両方を加えるので、全体として \frac{1}{3} N^3 数値的にも安定 # そもそも計算量が少いので、誤差が入りにくい # 代入法は、内積の計算になるので、安定した方法を利用している # 消去法 ( ガウス = ジョルダン 法 ) は、前進消去 + 後退消去 # 計算量は \frac{2}{3} N^3 なので、約 2 倍遅くなる ## 普段、演習で解く問題であれば、2 倍は大した差ではないが.. ## 本当の問題は、7000 x 7000 の組行列 ( 境界要素法 ) なら 5 時間掛る ## 2 倍だと 10 時間となるので、大変 !! == 世の中には、様々な解法がある # でも、騙されてはいけない !! 基本は、ガウスの消去法で、その変形でしかない !! 例: 消去法の三重ループは、色々と変形可能 例えば、Loop の K, I, J の順番 ( 6 通り ) を変えた物が別の名前に.. # 数学的には全く同じ # 計算機のアーキテクチャによっては、速度に影響がある場合もある ## メモリの構造によって、Cache がヒットするか/しないかという差が.. # ガウスの時代 ( = ベートベンの時代 ) 200 年位前の人 == これまでは、密(デンス)行列 ( 全ての要素が基本的に 0 でない ) 物を考えてきた # 本当の問題では、「境界要素法」の場合位しかない。 一般的な問題では、組(スパース)行列が殆ど 帯行列 ( 連立微分方程式 ) 対角要素とその前後だけが非零 一次元の微分方程式を差分化して方程式 対角要素とその前後だけ 三重対角行列になる 二次元の微分方程式を差分化して方程式 帯が 5 本になる # 連続体に関係する微分方程式を解く場合 ランダムスパース行列 非零要素が疎にある # ネットワークに関係する場合 通信ネットワーク 電気回路 # ○○網 対称か非対称かによって、解法が異る # 対称 def a_ij = a_ji ## これまでは、基本的に、非対称を考えている ## 物理的な問題をモデル化すると、通常は、対称になっていることが多い ## 物理現象には、「保存則」があり、それが、「対称」と関係がある ## # 右から見ても左からみても同じ ## # 吸い込み口とか、掃き出し口があると、保存則が成立しない ## ## 対称性しょうじない 対称なデータを扱う方法 コレスキー法 # 原理的には、(ガウスの消去法に比べ..) 半分で済む ## ただ、メモリをサボる ( 1/2 ) のために、プロ ## グラムは少し、複雑になる == 線型方程式の胆 「基本操作の繰り返しでできる」ということ == 演習 ( 2004/06/15 ) ガウスの消去法で、連立方程式を解く データ text p.38 を利用する [教科書への注意] 教科書のプログラム matrix.h を別ファイルで作成し include しないといけない ピボット選択 消去法を行う場合、分母に来る要素 ( ピボット ) の値が 0 だと計算できない # 0 でない要素を探す必要がある # ピボット選択 ## ピボットが選択できない場合 ## 実は、根がない場合 ( 行列式が 0 ) 例: 10 x_1 + x_3 = 1 10 x_2 + 3 x_3 = 10 3 x_1 + 10 x_3 = 20 は普通に解けるが ( 一行目と二行目を交換した.. ) 10 x_2 + 3 x_3 = 10 10 x_1 + x_3 = 1 3 x_1 + 10 x_3 = 20 は、素直には解けない ( 途中で、0 割がおきる ) # 数学的に同じ方程式なので、解けたり解けなかったりするのは変 # 都合が悪い場合は、行を入れ替えてよい # ピボット選択 !! ## ピボットとしては、要素の中で絶対値の一番大きい物にする ## 数値的に安定な結果になる ## 小さい物で割り算をすると数値的に不安定になる ## # 物理的な問題を定式化した場合は、普通安定になる ## # 勝手に作ると、不安定な結果になる可能性がある ピボット選択では、数値の正規化も行う必要がある 10 x_2 + 3 x_3 = 10 10 x_1 + x_3 = 1 3 x_1 + 10 x_3 = 20 なら、二行めの 10 だ (正しい) が、( 三行目を 1000 倍した ) 10 x_2 + 3 x_3 = 10 10 x_1 + x_3 = 1 3000 x_1 + 10000 x_3 = 20000 だと、単純に、絶対値の大きいものを選択すると、 (誤った 3000 ) を選んでしまう # スケーリングを行う必要がある # 各行の係数の最大値で行の係数を割り、 # スケーリング ( 正規化 ) する # ピボット選択をするならスケーリングする !! # ピボット選択しないならスケーリングには、意味がない