数値・記号処理(4) 方程式を解く(一変数) 慶應義塾大学理工学部 櫻井彰人 線型方程式を解く 線型方程式を解くことは簡単である。 中学の数学に出てくる 小学生高学年だって理解できる Cramer の公式で、一瞬! 本当に簡単か? 例えば、100変数ならどうか。 100 x 100 の行列式? 一般的に解く。解けることは分かる 係数を記号にして、四則のみで、解が記述可能 Mathematica で解いてみよう http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear01.nb クラメル(Cramer)の公式 連立一次方程式の基本的な解法であった 2 2 1 x 5 1 1 3 y 7 5 1 2 z 13 x 5 2 1 2 5 1 2 2 5 7 1 1 7 3 1 1 7 3 13 1 2 2 2 1 y 5 13 2 2 2 1 z 5 2 1 13 2 1 1 1 3 1 1 3 1 1 3 5 2 5 2 5 2 1 1 1 Cramer の公式(一般の場合) Cramer の公式 (Cramer’s rule) を使うと a11 x1 a2 x2 an xn d1 a21 x1 a22 x2 a2 n xn d 2 am1 x1 am 2 x2 amn xn d m a11 a12 a1n a21 a22 a2 n A am1 am 2 am xi A1 d1 d2 a12 a1n a22 a2 n dn an 2 ann A2 a11 d1 a1n a21 d 2 a2 n an1 dn ann Ai A http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear02-Cramer.nb Cramer の公式の問題点 行列式の計算が大変 n行 n列の行列式の計算に (n-1) n! 回の乗算が必要 丸め誤差の影響がありうる n! 個の項の加減 http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear03-det.nb 数値解析で用いる連立1次方程式の特徴 未知数が非常に多い 対称行列になることが多い 帯行列になることが多い * * * 0 0 0 0 0 0 0 * * * * 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 * * * * 0 0 0 0 0 0 0 * * * 解法の分類 直接法 何らかの形で直接的に解く。理論的には有限回演算。丸め誤差 の可能性 Cramer 公式 消去法 直交化法 せいぜい6元 Gauss法 一般に使える 三角化 Crout法 Choleski法 正定値対称の場合は標準的 掃き出し法(Gauss-Jorda法) 分かりやすい方法 安定性良. 複雑 反復法 トライアル・アンド・エラ-で解く。収束すれば丸め誤差の影響な し。計算の打ち切りが必要 Gauss-Seidel法 線形反復法 逐次過剰緩和法 Chebyshev反復 最急降下 勾配法 共役勾配 収束遅い 普通 速いが複雑 ナイーブ故 有限回 Gaussの消去法の例 Step 1 第一行を何倍かし、それを第二行以下に加減算し、 第二行以下の第一列成分を0とする 1 1 0 u1 1 4 4 1 2 0 1 2 u2 1 8 1 0 2 1 2 u3 1 8 0 1 2 1 2 1 u4 1 12 1 0 u1 1 4 4 1 0 7 4 1 4 1 2 u 3 16 2 0 1 4 7 4 1 2 u3 3 16 1 u4 1 12 0 1 2 1 2 Gaussの消去法の例 Step 2 第二行を何倍かし、それを第三行以下に加減算し、 第三行以下の第一、二列成分を0とする 1 0 u1 1 4 4 1 0 7 4 1 4 1 2 u 3 16 2 0 1 4 7 4 1 2 u3 3 16 0 1 2 1 2 1 u4 1 12 1 0 u1 1 4 4 1 0 7 4 1 4 1 2 u 3 16 2 0 0 12 7 4 7 u3 3 14 0 0 4 7 6 7 u4 23 168 Gaussの消去法の例 Step 3 第三行を何倍かし、それを第四行以下に加減算し、第四 行以下の第一、二、三列成分を0とする 1 0 u1 1 4 4 1 0 7 4 1 4 1 2 u 3 16 2 0 0 12 7 4 7 u3 3 14 0 0 4 7 6 7 u4 23 168 0 u1 1 4 4 1 1 0 7 4 1 4 1 2 u 3 16 2 0 0 12 7 4 7 u3 3 14 0 0 0 2 3 u4 5 24 Gaussの消去法の例 Step 4 u4,u3,u2,u1の順に解を求めることができる 0 u1 1 4 4 1 1 0 7 4 1 4 1 2 u 3 16 2 0 0 12 7 4 7 u3 3 14 0 0 0 2 3 u4 5 24 5 3 5 u4 24 2 16 3 14 4 7 u4 3 4 5 7 11 u3 12 7 14 7 16 12 48 3 1 11 1 5 4 11 u2 16 4 48 2 16 7 48 11 11 1 17 1 u1 1 1 48 48 4 96 4 Gaussの消去法のプログラミング 通常のように a11 a12 a a22 21 a31 a32 am1 am 2 a13 a23 a33 am 3 a1m u1 f1 a 2 m u2 f 2 a3 m u 3 f 3 amm um f m 前進消去 左下部分を0に i 1 m 1 a11 a12 a a22 21 a31 a32 am1 am 2 j i 1 m a ji f j f j f i aii k i 1 m a ji a jk a jk aik aii a13 a23 a33 am 3 a1m u1 f1 a 2 m u2 f 2 a3 m u 3 f 3 amm um f m 後退代入 解を、下から、順番に a11 a12 0 a 22 0 0 0 0 fm um amm i m 1 1 ( i が減少する順) j i 1 m f i f i aij f j fi ui aii a13 a1m u1 f1 a23 a2 m u2 f 2 a33 a3m u 3 f 3 0 amm um f m プログラム例 http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear04-GE.nb ガウスの消去法の計算量(1) 演算1 mi=aik/akk i=k+1, k+2, ..., n 演算2 aij = aij - mi×akj i=k+1, k+2, ..., n かつ j=k+1, k+2, ..., n 演算3 bi = bi - mi×bk i=k+1, k+2, ..., n n n n A x = b 0 n 0 k akk n-k ガウスの消去法の計算量(2) 前進部分での実行回数 乗除算 n(n-1) + n(n-1)(2n-1)/6 加減算 n(n-1)/2 + n(n-1)(2n-1)/6 後退部分での実行回数 乗除算 n(n-1)/2 + n 加減算 n(n-1)/2 行列式を求める 先ほど、行列式を求めるのは、計算量が大変といった。 しかし、ガウスの消去法を用いれば簡単。なぜなら a12 a11 a22 A a1n a13 a2 n a23 a3 n a33 ann a22 a33 ann det A det A a11 しかしなぜ? 行列式を求める 答え: 行列 A のある行の何倍かを他の行に加えて得ら れる行列を A´ とすると det(A) = det(A´) それ(「行の何倍か、、、、」)はまた何故? 答え: 行列式の多重線形性 det(a1, a2,.., ai+ai´,.., an) = det(a1, a2,.., ai,.., an) + det(a1, a2,.., ai´,.., an) による ガウス・ジョルダンの消去法 ガウスの消去法との違い: ガウスの消去では 対角成分はもとのまま (対角成分を除く)下三角成分のみを0とする ガウス・ジョルダン消去法では 対角成分を1にし、 対角成分以外を0にする 後退代入は不要となる! http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear05-GJ.nb 逆行列を求める ガウス・ジョルダン消去法を繰り返せばよい (LU分解と呼ばれる方法の方が効率がよい) 実際には、一気に計算する a11 a12 a21 a22 a n1 an 2 a11 a12 a 21 a22 an1 an 2 a1n c11 c12 a2 n c21 c22 ann cn1 cn 2 a1n c11 1 a2 n c21 0 ann cn1 0 c1n 1 0 c2 n 0 1 cnn 0 0 a11 a12 a 21 a22 an1 an 2 0 0 1 0 a1n c1i a2 n c2i 1 ann cni 0 逆行列を計算する 実際には一気に計算する a11 a12 a21 a22 a n1 an 2 1 0 0 0 1 0 a1n c11 c12 a2 n c21 c22 ann cn1 cn 2 0 c11 c12 0 c21 c22 1 cn1 cn 2 c1n 1 0 c2 n 0 1 cnn 0 0 a12 c1n a11 a22 c2 n a21 cnn an1 an 2 0 0 1 a1n a2 n ann http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear05-GJ.nb ピボット選択 ピボット: ガウスの消去法で「対角要素」が 0 で あったらどうしよう? 答え: 別の行を持ってくる(絶対値最大が安全) 悪条件 ill-conditioned の方程式 2x1 – x2 = 3 2.1x1 – x2 = 3
© Copyright 2024 ExpyDoc