-
Cプログラミング ー連立1次方程式の解法一 早稲田大学
-
本日の目標 ●連立一次方程式の解法 ●ガウスの消去法 ●前進消去 ●後退代入 ●NaN・INFINITY ・部分ピホット選択付きカウスの消去法
-
連立1次方程式の解法 ●逆行列を用いる方法 ●クラメール(Cremer)の公式 ●直接解法(消去法とも呼はれる) ●カウス〔Gauss)の消去法 .LU分解 ●反復解法 ●ヤコビ(」acob;)法 密行列の計算に向く. 大規模な疎行列の解法に向く, ・ガウス・ザイデル(Gauss-Seldel)法
-
ガウスの消去法 例題1 次の3元1次連立方程式をガウスの消去法を用いて解け. {i妊i縷≡ii
-
ガウスの消去法 { 2T 十4tl 3砿 十Sv st 十7む 十6z 十7s 十21z = 6 = 15 = 24 ‘1) ② (3) 妙 〔2)一(1)X量 【〔2)からxを消去する】 { 2x 十4∬ Ox 十2y sx 十7Ψ 十6x -2z 十21x = 6000(1) = 6 ti&b(4) = 24 trE・〔3) 年 (3)一(1)x巷 【〔3)からzを消去する】 { 2¢ 十4y Oz 十2tJ Oa -3tJ 十6z -2z 十6z )〕) 145 (〔〔 いO 曾 b O O ワ 669 =;= ↓ 〔5)一(4)x:ヂ【(5>からyを消去する】 { 2α +4馴 0= ÷2y O= 十〇y +6z -2z 十3z = 6・・(1> = 6 叩。(4> = 18 。月口(6>
-
ガウスの消去法 (il茎)(の一㈲ 妙 1行目から0行目×§を引く ) 4 66 2 ( = ー ¢びZ ( ー 6辺21 4ワ一7 205 ー 年 2行目から0行目x巷を引く ー 669 ー = ) τΨZ ー ー 2 6一6 42弓 200 ー ↓ 2行目から1行目xT3を引く ) 8 66 1 ( ) ¢VZ ー ー 2 6一3 420 200 ( 前頁の連立方程式の変形を行列形式で書き直し た 以後,この書き方で説明する
-
前進消去 り ●連立1次方程式Am-=bを変形して,行列Aを上三角行列の形に する. ●0列目について1行目以降の要素を0に,1列目について2行目以降 の要素をOに,という操作を繰リ返す.今,O,,k-1列目の処理 は済んでいるとして,k列目のk+1行目以降をOにすることを考え る.(次の例はN=5,k=1) le L 催ii:離iii iiii〕ii〕一〔ii}・聡11iii iiii liii〕〔ii〕・ 勿論、変形の前後でAとbの値は変化する(同じ記号を使っているので注意)。 ii〕
-
前進消去 ●更にk列目の中のk+1,,z-1行目の処理は済んでいるとして. z 行目の処理を考える.(次の例はN=5, k-=⊥,t=3) k 、fl:i;:li;iミi;i li;i 0123 四昭Z昭 01234 bOOOb ⇒ L 鮒i:lli;i Ao.3 At 3 A2.A A3.3 Ad.3 Ao.己1 滝1,己I A2.4 A3コ己1 鴻4コ珊 01234 躍冨冨冨躍 01234 しhいhし
-
前進消去 ●A3,1を0にするには.3行目から1行目の会1:}倍を引けばよい. 、.α ↓k⇒ li:渥ii藝iii liii〕〔ii〕一〔1〕 ↑s⇒ r;会i,}として.7z1,,4についてA3“←A3.2-Al.lr. b3←b3-b17とすれ ばよい. これを一般化すると,r;会慧として」・・ k,,N-1について A。.s←A,.」-AL.sr, b、←b、-bt.rとすればよい.
-
前進消去 プログラムの概形は次のような3重ループになる: for (k = 0, k < N - 1, k++) { /* k=0,Q。。,N-2 */ for (1 = k + 1, ⊥ く N, 1++) { /* 1=k+1,0 0,N-1 */ r昌A[1][k]/A[k][k], for(〕=k, jくN, J++){ /*J=k,・。。,N-1*/ A[L][コ]一=A[k][j]*r, } b[エ]一=b[k]*r, } }
-
後退代入 ●連立1次方程式且5=6を更に変形して,Aを上三角行列から単位行列の 形にしたい そうすれば,至=bとなり,bが答になる ●必要なのはbの値だけなので,実際のアルゴリズムではAの値の更新は省 略ざれる, ●N-1,,z+⊥行目の処理は済んでいるとして, t行目の処理を考える(次 の例はz=2) LC轟:il…鱗……〕〔ii〕一〔1〕 ↑フ⇒ b2←(b2-A2、3b3-A2,4b4)1A2,2とすればよい. ・れを一般化すると妊(b,一Σ鑑‡lAz“bs)/A・,tとすればよ…
-
後退代入 プログラムは次のような2重ループになる. for (1 = N - 1, ユ 〉= 0, ユー一) { /* 1=N-1, 。。,0 */ for (コ = 1 + 1, コ く N, 」++) { /* 」=ユ+1,0。。,N-1 */ b[1]一=A[1][J]*b[J], } b[1] /昌A[1][1], }
-
例題1 ●最後にbに解が入ることになる
-
例題1のヒント (概形) #⊥nGlude <3td⊥o h> #⊥nGlude く3tdl±bゆ vOld PrlntAb(doubl臼 *A, dcuble sb,1Et E){ /1行列・ベクトルを表示する関数*t } lnt 皿a1皿(VOId) { m七 1馳 」L k、 N. doub18 rp double 掌a五, 零誼印 double 掌A, #b, prlntfρ1可i剛、. sc皿f(ll旧ll,thN), A = (double 掌) 国a上100(NIVSs±乙ee±(aoub⊥e)), ±i(A==NULL){ pr・鵬f(1’C皿’t a11・c&t8皿肥。ry≒n唱’), 8×1t(1), } x = 〔double 零) ma⊥100(N*s⊥zee:(double)), ±i〔x==肌L){ przntf〔.’CanP七 a⊥100邑te memery瓶”), 旧Xlt(1), } } e⊥ =八. fDr( l i O . 1 〈 N . ユ中ウ 〉{ fDr(JiO.」く閥.」+ウ){ prlntf⊂1印貞[祖]鵬司 i I/t 1 . 」 ). 5c皿(’.Xlf”, a±+」 ), } a±+=N, } te=( ± = O , ± く 畑 , ⊥十十 ){ prin±±r..b[盟d] = 1.. ± ]」 5G㎝£([Xlf.., b十t), } PrmtAb〔A. b. N), 1*th進消去*1 〆ホ後退代入1! rEtum O.
-
avge 1 utfieas tt-EludiAtsbdtnhi #t-eludie tsvdltb ta vbtd trOtlb[devb!o --t deuble -bt int n]{ { tncttjt diavma--tit pminttC-Ab-in-} toiCtre 1" i+tH torCl-e jo j++H pmintt{'TGit#' .Cai.j)) } prh-ftltAtf-:-.brttl. atimt ] pminttC'#'} } intuatmCvold){ tmtt,n,U,H, rm rt rm -tit -abt rm -At -ht pmintir-h-t scant(-ht -nt A . (detht- -] =slleca-Hietmer(deth1-})t tt[Amu]{ prdmtT( Can t ableen[e menony h'}, -:ttCl} } I.(doab1-Dulloca-sizmof(dedble)) tttr"HULL){ prtntT(Cantst1-ewhemansrsh'}, {:tt{1), L } at-l, t"[t-o,±c",t"){ ter{j-o j{ff,j+.}{ pricttc'1[bg[abg-pi,j), s:antC:lt al.]} t U-J, } tcr[t-ottc"ttt-){ prilttC'bua - ' i ), .ca-t(xlt ei) ) ltintAbCA b,H), f,bltHt.t at.lt ab-lt tcr (h.O k{Y-i k-){ tor Ci -ktt i{H t++H r-gal+t-H+h)1-Cab-V tw (j.h jcY l-){ .[ti--Hl}.--(suJ]"r, } btil il bix1 t r. ] ennt-bCA b n sk-n } i,vafltelkv ti-akt tcr(t-".T,t)-O,t"]{ ter(1-t-t,J(N,j-}{ li[i]-l(al+j)lb"] } bri7 i- -raf-il enort-bCi b N} al."t } rrmOt
-
avwa1 e RtivaMcStaop at b t[ U4, Ab ' 2 OO 3 OO s oo A b! ' 2 OO o oo o oe A b! ' 2 OO o oo o oe 4 s 7 4 2 -3 4 2 e oe ee oe oe oe oe oe oe ee 6 7 2t 6 -2 6 6 -2 3 oo eo eo oo eo eo oo oo eo ] ] d ] 1 1 ] 1 ] 6 t5 24 6 6 9 6 6 t8 eo oe oe oo oe oe oo oo oe Ab ' 2 OO e oe o oo A b! ' 2 OO e oo o oe A b! ' 2 OO s oo e oe 4 2 o 4 2 o 4 2 o oo oe oe oo oe oe eo oo ee s -2 a 6 -2 3 5 -2 3 oo oo oo oo oo oo oo oo oo 1 1 1 6 6 6 6 9 6 1-33 lg 16 oo oe oe oo oo oe oo oo oe
-
例題2
-
ewwa 2 . RtmeMtStaopat bt[U4, A b- ' 2OO -t oe 4oe soe A b= ' 2OO ooo ooe ooe -b= ' 2oe ooe ooo ooo A b= ' 2oe ooe ooo ooo qoo -2 OO 2OO -4 OO 4 o -B -t4 oo oo co oo 4OO ooo ]an 1ff[ 4OO ooo ]an =an 1eo 2OO -s eo -s oo 1eo 250 -s eo -S EO too 2SO 1ti iti too 2EO 1nf ]an -3 oo 1oee 4oe 1io oo 5oe ]2oo ioe ]eoo -3 oo 1oeo 2sc 1loeo 11 00 /200 ese ]6oo -9 oe ]ooo 2sc 1loeo i]i 1 Snf i]f 1 Saf -s oe 1 om 2SO 1iO OO i]i I Snf nan 1 ]u] A b- ' 2oo ooo ooo ooo A b= ' 2oo ooo ooo ooo - b= ' 2oo ooo ooo ooo A b= ' 2oo ooo ooo ooo 4OO ooo nan nan 4OO ooo ]an nan 4OO ooo ]an ]an 4OO ooo ]an ]an !oo 2SO ttt ]an 1oo 2EO 1nf ]an too 250 ±nf ]a] too 250 ±nf ]an -3 2 -3 2 -s 2 -3 2 oo se tttt eal oo so i]f nan oe so ilf DSI oe se inf Dal 1ooo 1 iO OO 1 tnt 1 nan 1ooo 1loeo 1 1an 1 1an 1Ooo 1 ]an 1 1an 1 1an d nan d ]an 1 1an 1 1an rn ]1 llP rllT t Sfioli7 ffftu LLI ,, xlt.;.-J pt iOh 7
-
NaN・INFINITY ●d。uble型のような浮動小数点を扱う際多くの計算機ではlEEE754 と呼ばれる2進浮動小数点の規格に従って処理が行われる. ●]EEE754では,計算過程において発生する可能性がある数の中で, 通常の浮動小数点数では表現できない特別な値が定義されている, 例NaN,+INFINITY,-INFINITY
-
NaN・INFINITY 非数NaN〔Not a Number) ●010,Aなとの無効演算の結果として生成される特殊な数で、 nan と表示ざれる ●nanを含む計算の結果はnanとなる.例⊥+nam=nan
-
ewpa2 e NtuE LLNfieffstfiSnigLNOh, ? . tttmiNfioza 1 yJt,uxe[tsLsT, enMTcs A,,, = o opee, r - ::/: oporsuTO TO:UfiD#tL. r opblsureMh'-mf kts6. A,b- 2OO -1 OO 4oe soe A, b= 2OO ooo ooo ooo A,b= 2on ooe ooe ooo A b= 2oe ooe ooe ooo 4oo leo -] oo ]eo 2oo -s eo -4 00 -S OO 4 o -6 -1- oo leo oe 2Eo oo -s eo oo -s so 4on ooo IEP ]an 4OO ooo zan =an tno 2EO IUi 1ti teo 250 tEt ]an -3 oo 1oeo 4oo 1lo eo 6oe 12oo toe ]6oo -3 oo 1oeo 2SO 110 OO 11 oo 12eo Ssc16ee -9 0e /ooo 2SO 1iO oo i]f 1 inf i]i 1 ±nf -a oe iooo 2se 1io oo izt 1 int zam 1 ]an A b- ' 2oo ooo ooo ooo A, b= 2oo ooo ooo ooo A,b= 2oo ooo ooo oon A b= ' 2oo ooo ooo ooo aoo ooo nan nan 4OO ooo nan ]an 4on ooo ]EP ]an 4OO ooo nm ]an !oo ]50 tnt ]an 1oo 250 1nf ]a] too 2SO ±nf ]a] too 250 tnt ]a] -3001O 2SO 110 ittt 1 nan 1 -3001O 2SO 110 i]f 1 Dan 1 -son1o 2se d inf 1 um 1 -s oe 1 2se d trt d Dan 1 oo eo tnt nan oo eo 1an la] no ]an -an la] ]an ]an man lal
-
部分ピボット選択付きガウスの消去法 ●何故正しい解が計算ざれないのか? ●前進消去の第k列消去において,Ak,k(例題ではAl,1)=0の場合, r一舞の計算で・での割算が発生し・計算結果か一・n・となる・ L “9・°譜 :1・:::・1:?・l ili’備ili:iliii 01234 匹四四四罵 ⇒ 0 00000 ^ ↓k AO脚1 AO脚2 AOβ AOr4 ユユ t,2 t At,ri l 盆・:震:1盆:1 0 A【脚2 A【β A【脚4 01234 匹躍匹四罵 D1234 しbhhし ●方程式の順番を入れ替えても解は変わらないので,行を入れ替えて Al,1≠0となるようにすればよい. ●絶対値が0に近い値で割ると,一般に計算精度が悪くなるため,常 に1行目とA11,A21,A3、1,A4,1の中で絶対値が最大なものの行との 入れ替えを行うことにすれば更によい
-
部分ピボット選択付きガウスの消去法 部分ピボット選択 前進消去の第k列消去において,(Akk≠Oであっても常に)k行目と 塩,κ,,AN_1,kの中で絶対値か最大なものの行との入れ替えを行う. 即ち,k行目とm≡maxF鳶,,N_Ll.4剥行目を入れ替える. ●例えは,Al.1,A2、⊥,A31,A41がそれぞれO,3,-4,2だとしたら,1行 目と3行目を入れ替えてから,第k列消去を行う. 』:鰐ii灘iii liii〕〔1〕一〔ii) k行目とm行目の入れ替えを行う 勿論Aの行たけでなく,bの対応する要素も入れ替える 列は入れ替えないので,fの要素は入れ替えてはならない ●たまたまk=mの場合は,入れ替え作業は不要,
-
配列内の最大要素の探索 要素数Nのdouble型配列の何番目に最大値が入っているかを求めるプロ グラム ●次のように探索する方法が一般的.例文内のコメントをよく理解す ること. #lnclude <Stdlo h> #detlne N 5 1nt maln(VOld) { ユnt m7 1t double A脚] = {2, 1, 5, 3, 4}, 皿=0, 〆*仮に0番目が最大だとする掌/ for(ユ=1,1《Nt 1++){ 〆*1番目から最後まで調べる*/ li(A[1]>Ain]){ 〆*もし1番目の方か大きいなら*! m旨1/ 〆*エ番を記憶*/ } } prlntf(11A[Xd] = Xf Is max、nll」 皿: A回). r巳七urn O, }
-
部分ピボット選択付きガウスの消去法 ●プログラムの概形 /ホ前進消去ホ/ for(k=0, kくN-1, k++){ /*部分ピボノト選択*/ m=k, ん仮の最大番号*/ for(ユik+i,1く眺1++){ 省略 } prlntf(「lco1皿’Ad row Zd ls 皿ax、n1㌧ k. m), li(m b=k){ /*va≠kの時入れ替え*! 省略(入れ替えの関数 svap) } 〆*第k列}自去*ノ for (1 = k + 1, 1 く N, ユ++) { f。r(j=k,」<N,)++){ 省略 } } PrlntAb(A} b), }
-
アルゴリズムの実装手順
-
本日のまとめ ●連立一次方程式の解法 ●ガウスの消去法 ●前進消去 ●後退代入 ●NaN・INFINITY ・部分ピホット選択付きカウスの消去法