2005/06/28 数値解析学と演習 福井 先生 先週は、連立方程式の解法 : ガウスの消去法 連立方程式 Ax = b が与えられた時に、A を対角行列 A' の形に変形し、 A'x = b' となれば、「解けた」とみなす ( ガウス = ジョルダン 法 )。 # 前進消去 + 後退消去 ポイントは、 Ax = b の根と A'x = b' の根が等しいように A から A' を作る手続き => 基本変形 [注意] 基本変形を旨い順番で適用すれば、一度 0 にした係数は、0 のままにできる # 後期に行う、固有値問題では、こうは上手く行かない.. ガウスの消去法 前進消去 + 後退代入 [計算量の話] ガウス = ジョルダン 法 前進消去 ( \frac{1}{3}n^3 ) + 後退消去 ( \frac{1}{3}n^3 ) ガウスの消去法 前進消去 ( \frac{1}{3}n^3 ) + 後退代入 ( \frac{1}{2}n^2 ) # n が大きくなれば、ガウスの消去法は O(\frac{1}{3}n^3), == [本日の話] LU 分解 # 前回、「連立方程式を解くために、逆行列を求める」のは間違い # 計算量は多い / 誤差も多い # # 本当に、逆行列の係数が欲しいくない限り、良いところはない !! ## しかし、Ax=b_i のように b_i が沢山あった場合はどうする.. ? ## => LU 分解を使う !! Ax = b_i で b_i が沢山あったら... もし、予め b_i が解っている場合.. A( x_1, x_2 .., x_n ) = ( b_1, b_2, .., b_n ) と並行して行えばよい.. しかし、 Ax=b_0 => x = b_1 Ax=b_1 => x = b_2 .. のような形で、一旦、方程式を解いた結果を次の問題で使うような場合もある => 並列にできない !! # すると、やっぱり、逆行列を使いたくなるが.. => 不要 LU 分解を使う [LU 分解] 行列 A を、二つの行列 L と U の積の形にする => \frac{1}{3}n^3 だけかかる L は、下三角行列 U は、上三角行列 LU 分解が済んでいれば.. => b_i が与えられたときに、解は n^2 で求めることができる。 # どうやって解くか ? A = LU と分解できたとして.. Ly = b は、( A が四次とすれば.. ) l_11 y_1 = b_1 l_21 y_1 + l_22 y_2 = b_2 l_31 y_1 + l_32 y_2 + l_33 y_3 = b_3 l_41 y_1 + l_42 y_2 + l_43 y_3 + l_44 y_4 = b_4 の形なので、(前進代入を使って..) y を求めることができる => \frac{1}{2}n^2 更に、 Ux = y は、 u_11 x_1 + u_12 x_2 + u_13 x_3 + u_14 x_4 = y_1 u_21 x_1 + u_22 x_2 + u_23 x_3 = y_2 u_31 x_1 + u_32 x_2 = y_3 u_11 x_1 = y_4 の形なので、やはり x を求めることができる => \frac{1}{2}n^2 # 計算的にも安定性が良い。 # A はどうやって、LU の分解をするか ? LU 分解はどうする ? ( Tex には、クラウト法を使えばよい書いてあるが.. ) => 本質的には、ガウスの消去法と同じ !! => 手順を覚えておくだけ !! まず、Ax=b の形をガウスの消去法で解くことを考え、前進消去を行う。 a_11 x_1 + a_12 x_2 + a_13 x_3 + a_14 x_4 = b_1 a_21 x_1 + a_22 x_2 + a_23 x_3 + a_24 x_4 = b_2 a_31 x_1 + a_32 x_2 + a_33 x_3 + a_34 x_4 = b_3 a_41 x_1 + a_42 x_2 + a_43 x_3 + a_44 x_4 = b_4 a_21 を消すには、 (a_21 + α_21 a_11)x_1 = 0 となるようにするとよいので、 α_21 = - \frac{a_11}{a_21} となる。 ... 一般に、前進消去で、項目 a_ij を消すために、利用される係数 α_ij をとっておけば、 それが l_ij になり、 また、 消去した残りの部分が u_ij になっている。 ただし、l_ii = 1 となる。 # L の係数の個数は、\frac{n(n+1)}{2}, U の係数の個数は、 # \frac{n(n+1)}{2} なので、全部で n(n+1) 個あるが、A # の要素は nn 個。しかし、l_ii = 1 なので、これは記憶 # しなくてよいと思うと、元の A の下に L の要素 ( ただ # し、対角要素は保存しない ) 、A の上に U の要素を保存 # すれば、元の A の中に、L と U を保存することができる。 # クラウト法の中身を良く読むと、実は、ガウスの消去法の係数を保存しているだけ.. ## LU 分解のプログラムを沢山作っている間に、LU 分解とガウスの消去法の関係に気が付いた.. ガウスの消去方法 一列目は、一回の計算の結果 二列目は、二回の計算の結果 .. => n 列目は n 回の計算を累積 (凝縮) した結果を作ってうる => 計算量が減るタネ # 結局、ガウスの消去法と、クラウト法は同じような物 # 係数を処理する順番や、L と U の何方の対角を 1 にするかのバリエーションがある.. [枢軸選択の問題(ピボット選択)] 計算の途中で、a_ii が 0 になってしまったらどうするか ? # 見かけ上、これは、解けなくなってしまうように見える => どうする ? 方程式の順番を変更する => ピボット選択 : 対角が零にならないようにする # ピボット選択は、誤差を小くするためにも利用されるが # => できるだけ絶対値の大きな値を a_ii に使うのが望ましい。 # 誤差の事を考えると、更に、係数の中で、最も絶対値の大きいものにすると、最も望ましい。 # 式の順番の変更は、全く問題ないが、係数を変更すると、解の順序が変るので注意する # ピボット選択はコストがかかるので、問題を作る段階で、上手く係数を選択し、そもそもピボット選択が不要になるように方程式をたてる工夫をすることが多い。 # ここらへんの話を始めると、(得意な分野なので..) 1 年位続いちゃうので、ここらへんんで止める.. # 上は、数学的な話、以下は、物理的な話 川の流れは、上の方から水が流れてくるので、上流の方が原因になっている 上流の方を上の式にして、連立方程式を解くようにする 物理的に、情報が伝わる方向に方程式を解く => 答えが正しく出やすい 紐を引っ張る問題 力を加えるのは、紐の下の端 下の方 ( 力の係る方 : 原因のある方 ) から計算すると解ける => 逆向きだと、上手く解けない ( 不安定に成りやすい ) # 原因から、結果に向けて計算するとよい !! # 後期で行う偏微分方程式でも、同様に、流れの上 ( 風上 ) から下へ計算した方がよい 計算を行う順番は、「数学的」には、どの順番でも同じはずなのだが.. 「物理的」に正い順番に計算した方が、途中の計算結果を眺めていると、安定していることが解る。 [メッシュの例] --------- => ( 川の流れ ) --------- |/|/|/|/|/| |/|/|/|/|/| の切り方 : 上下の対称性がないので変な結果に.. |\|\|\|\|\| |/|/|/|/|/| の切り方 : 上下の対称性があるので、良い結果に.. # 元の問題に対称な性質があればメッシュも.. # 近似は、Δx を 0 に収束させれば、どちらも同じ ( 数学的には.. ) # 実際は、0 でないので、メッシュの影響が生じる # => 「自然な形(物理学的な意味)」が望ましい。 みんなは、もう、ガウスの消去法や、LU 分解を作ることはない 連立方程式をたてることが仕事 => 連立方程式をたてる時は、物理的な性質を着目する !! # ときたい問題が特殊な構造をもっていて、その性質を利用すると早くなる場合 # 自分で作ることになるかも.. # クラウト法を変形するより、ガウスの消去法を変形した方が楽 == [演習] ( Text ではクラウト法 ) # 本当は、先週のガウス法のプログラムを変形させて、LU にして欲しいが... [問] Text のクラウト法のプログラムを作成して、連立方程式を解く [係数] 前回の問題と同じ