Ｃプログラミング ー連立１次方程式の解法一 　　　　早稲田大学

本日の目標 ●連立一次方程式の解法 　　●ガウスの消去法 　　　　●前進消去 　　　　●後退代入 　　●ＮａＮ・ＩＮＦＩＮＩＴＹ 　　・部分ピホット選択付きカウスの消去法

連立１次方程式の解法 ●逆行列を用いる方法 ●クラメール（Ｃｒｅｍｅｒ）の公式 ●直接解法（消去法とも呼はれる） 　　●カウス〔Ｇａｕｓｓ）の消去法 　　．ＬＵ分解 ●反復解法 　　●ヤコビ（」ａｃｏｂ；）法 密行列の計算に向く． 大規模な疎行列の解法に向く， ・ガウス・ザイデル（Ｇａｕｓｓ－Ｓｅｌｄｅｌ）法

ガウスの消去法 例題１ 次の３元１次連立方程式をガウスの消去法を用いて解け． 　　　　　　　　　｛ｉ妊ｉ縷≡ｉｉ

ガウスの消去法 ｛ 　２Ｔ　十４ｔｌ 　３砿　十Ｓｖ 　ｓｔ　十７む 十６ｚ 十７ｓ 十２１ｚ ＝　　６ ＝　　１５ ＝　　２４ ‘１） ② （３） 妙　〔２）一（１）Ｘ量　【〔２）からｘを消去する】 ｛ 　２ｘ　十４∬ 　Ｏｘ　十２ｙ 　ｓｘ　十７Ψ 十６ｘ －２ｚ 十２１ｘ ＝　６０００（１） ＝　　６　ｔｉ＆ｂ（４） ＝　２４　ｔｒＥ・〔３） 年　（３）一（１）ｘ巷　【〔３）からｚを消去する】 ｛ 　２￠　十４ｙ 　Ｏｚ　十２ｔＪ 　Ｏａ　－３ｔＪ 十６ｚ －２ｚ 十６ｚ ）〕） １４５ （〔〔 いＯ　曾 　ｂ Ｏ　Ｏ　ワ ６６９ ＝；＝ ↓　〔５）一（４）ｘ：ヂ【（５＞からｙを消去する】 ｛ 　２α　＋４馴 　０＝　÷２ｙ 　Ｏ＝　十〇ｙ ＋６ｚ －２ｚ 十３ｚ ＝　６・・（１＞ ＝　　６　叩。（４＞ ＝　１８　。月口（６＞

ガウスの消去法 （ｉｌ茎）（の一㈲ 妙　１行目から０行目קを引く ） 　　４ ６６ 　　２ （ 　＝ ー ￠びＺ （ ー ６辺２１ ４ワ一７ ２０５ ー 年　２行目から０行目ｘ巷を引く ー ６６９ ー 　＝ ） τΨＺ ー ー 　２ ６一６ ４２弓 ２００ ー ↓　２行目から１行目ｘＴ３を引く ） 　　８ ６６ 　　１ （ ） ￠ＶＺ ー ー 　２ ６一３ ４２０ ２００ （ 前頁の連立方程式の変形を行列形式で書き直し た　以後，この書き方で説明する

前進消去 　　　　　　　　　　　り ●連立１次方程式Ａｍ－＝ｂを変形して，行列Ａを上三角行列の形に 　する． ●０列目について１行目以降の要素を０に，１列目について２行目以降 　の要素をＯに，という操作を繰リ返す．今，Ｏ，，ｋ－１列目の処理 　は済んでいるとして，ｋ列目のｋ＋１行目以降をＯにすることを考え 　る．（次の例はＮ＝５，ｋ＝１） 　　　　ｌｅ　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　Ｌ 催ｉｉ：離ｉｉｉ　ｉｉｉｉ〕ｉｉ〕一〔ｉｉ｝・聡１１ｉｉｉ　ｉｉｉｉ　ｌｉｉｉ〕〔ｉｉ〕・ 勿論、変形の前後でＡとｂの値は変化する（同じ記号を使っているので注意）。 ｉｉ〕

前進消去 ●更にｋ列目の中のｋ＋１，，ｚ－１行目の処理は済んでいるとして．　ｚ 　行目の処理を考える．（次の例はＮ＝５，　ｋ－＝⊥，ｔ＝３） 　　　　　ｋ 、ｆｌ：ｉ；：ｌｉ；ｉミｉ；ｉ　ｌｉ；ｉ ０１２３ 四昭Ｚ昭 ０１２３４ ｂＯＯＯｂ ⇒ 　　　Ｌ 鮒ｉ：ｌｌｉ；ｉ Ａｏ．３ Ａｔ　３ Ａ２．Ａ Ａ３．３ Ａｄ．３ Ａｏ．己１ 滝１，己Ｉ Ａ２．４ Ａ３コ己１ 鴻４コ珊 ０１２３４ 躍冨冨冨躍 ０１２３４ しｈいｈし

前進消去 ●Ａ３，１を０にするには．３行目から１行目の会１：｝倍を引けばよい． 、．α ↓ｋ⇒ ｌｉ：渥ｉｉ藝ｉｉｉ　ｌｉｉｉ〕〔ｉｉ〕一〔１〕 ↑ｓ⇒ ｒ；会ｉ，｝として．７ｚ１，，４についてＡ３“←Ａ３．２－Ａｌ．ｌｒ．　ｂ３←ｂ３－ｂ１７とすれ ばよい． これを一般化すると，ｒ；会慧として」・・　ｋ，，Ｎ－１について Ａ。．ｓ←Ａ，．」－ＡＬ．ｓｒ，　ｂ、←ｂ、－ｂｔ．ｒとすればよい．

前進消去 プログラムの概形は次のような３重ループになる： ｆｏｒ　（ｋ　＝　０，　ｋ　＜　Ｎ　－　１，　ｋ＋＋）　｛　　　／＊　ｋ＝０，Ｑ。。，Ｎ－２　＊／ 　ｆｏｒ　（１　＝　ｋ　＋　１，　⊥　く　Ｎ，　１＋＋）　｛　　／＊　１＝ｋ＋１，０　０，Ｎ－１　＊／ 　　ｒ昌Ａ［１］［ｋ］／Ａ［ｋ］［ｋ］， 　　ｆｏｒ（〕＝ｋ，　ｊくＮ，　Ｊ＋＋）｛　／＊Ｊ＝ｋ，・。。，Ｎ－１＊／ 　　　Ａ［Ｌ］［コ］一＝Ａ［ｋ］［ｊ］＊ｒ， 　　｝ 　　ｂ［エ］一＝ｂ［ｋ］＊ｒ， 　｝ ｝

後退代入 ●連立１次方程式且５＝６を更に変形して，Ａを上三角行列から単位行列の 　形にしたい　そうすれば，至＝ｂとなり，ｂが答になる ●必要なのはｂの値だけなので，実際のアルゴリズムではＡの値の更新は省 　略ざれる， ●Ｎ－１，，ｚ＋⊥行目の処理は済んでいるとして，　ｔ行目の処理を考える（次 　の例はｚ＝２） ＬＣ轟：ｉｌ…鱗……〕〔ｉｉ〕一〔１〕 　　　　　　　　　　　↑フ⇒ ｂ２←（ｂ２－Ａ２、３ｂ３－Ａ２，４ｂ４）１Ａ２，２とすればよい． ・れを一般化すると妊（ｂ，一Σ鑑‡ｌＡｚ“ｂｓ）／Ａ・，ｔとすればよ…

後退代入 プログラムは次のような２重ループになる． ｆｏｒ　（１　＝　Ｎ　－　１，　ユ　〉＝　０，　ユー一）　｛　　／＊　１＝Ｎ－１，　。。，０　＊／ 　ｆｏｒ　（コ　＝　１　＋　１，　コ　く　Ｎ，　」＋＋）　｛　　／＊　」＝ユ＋１，０。。，Ｎ－１　＊／ 　　ｂ［１］一＝Ａ［１］［Ｊ］＊ｂ［Ｊ］， 　｝ 　ｂ［１］　／昌Ａ［１］［１］， ｝

例題１ ●最後にｂに解が入ることになる

例題１のヒント （概形） ＃⊥ｎＧｌｕｄｅ　＜３ｔｄ⊥ｏ　ｈ＞ ＃⊥ｎＧｌｕｄｅ　く３ｔｄｌ±ｂゆ ｖＯｌｄ　ＰｒｌｎｔＡｂ（ｄｏｕｂｌ臼　＊Ａ，　ｄｃｕｂｌｅ　ｓｂ，１Ｅｔ　Ｅ）｛ 　　／１行列・ベクトルを表示する関数＊ｔ ｝ ｌｎｔ　皿ａ１皿（ＶＯＩｄ）　｛ 　　ｍ七　１馳　」Ｌ　ｋ、　Ｎ． 　　ｄｏｕｂ１８　ｒｐ 　　ｄｏｕｂｌｅ　掌ａ五，　零誼印 　　ｄｏｕｂｌｅ　掌Ａ，　＃ｂ， ｐｒｌｎｔｆρ１可ｉ剛、． ｓｃ皿ｆ（ｌｌ旧ｌｌ，ｔｈＮ）， Ａ　＝　（ｄｏｕｂｌｅ　掌）　国ａ上１００（ＮＩＶＳｓ±乙ｅｅ±（ａｏｕｂ⊥ｅ））， ±ｉ（Ａ＝＝ＮＵＬＬ）｛ 　　ｐｒ・鵬ｆ（１’Ｃ皿’ｔ　ａ１１・ｃ＆ｔ８皿肥。ｒｙ≒ｎ唱’）， 　　８×１ｔ（１）， ｝ ｘ　＝　〔ｄｏｕｂｌｅ　零）　ｍａ⊥１００（Ｎ＊ｓ⊥ｚｅｅ：（ｄｏｕｂｌｅ））， ±ｉ〔ｘ＝＝肌Ｌ）｛ 　　ｐｒｚｎｔｆ〔．’ＣａｎＰ七　ａ⊥１００邑ｔｅ　ｍｅｍｅｒｙ瓶”）， 　　旧Ｘｌｔ（１）， ｝ ｝ ｅ⊥　＝八． ｆＤｒ（　ｌ　ｉ　Ｏ　．　１　〈　Ｎ　．　ユ中ウ　〉｛ 　　　ｆＤｒ（ＪｉＯ．」く閥．」＋ウ）｛ 　　　　　ｐｒｌｎｔｆ⊂１印貞［祖］鵬司　ｉ　Ｉ／ｔ　１　．　」　）． 　　　　　５ｃ皿（’．Ｘｌｆ”，　ａ±＋」　）， 　　　｝ 　　　ａ±＋＝Ｎ， ｝ ｔｅ＝（　±　＝　Ｏ　，　±　く　畑　，　⊥十十　）｛ 　　　ｐｒｉｎ±±ｒ．．ｂ［盟ｄ］　＝　１．．　±　］」 　　　５Ｇ㎝￡（［Ｘｌｆ．．，　ｂ十ｔ）， ｝ ＰｒｍｔＡｂ〔Ａ．　ｂ．　Ｎ）， １＊ｔｈ進消去＊１ 〆ホ後退代入１！ ｒＥｔｕｍ　Ｏ．

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

例題２

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

ＮａＮ・ＩＮＦＩＮＩＴＹ ●ｄ。ｕｂｌｅ型のような浮動小数点を扱う際多くの計算機ではｌＥＥＥ７５４ 　と呼ばれる２進浮動小数点の規格に従って処理が行われる． ●］ＥＥＥ７５４では，計算過程において発生する可能性がある数の中で， 　通常の浮動小数点数では表現できない特別な値が定義されている， 　例ＮａＮ，＋ＩＮＦＩＮＩＴＹ，－ＩＮＦＩＮＩＴＹ

ＮａＮ・ＩＮＦＩＮＩＴＹ 非数ＮａＮ〔Ｎｏｔ　ａ　Ｎｕｍｂｅｒ） 　●０１０，Ａなとの無効演算の結果として生成される特殊な数で、　ｎａｎ 　　と表示ざれる 　●ｎａｎを含む計算の結果はｎａｎとなる．例⊥＋ｎａｍ＝ｎａｎ

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

部分ピボット選択付きガウスの消去法 ●何故正しい解が計算ざれないのか？ ●前進消去の第ｋ列消去において，Ａｋ，ｋ（例題ではＡｌ，１）＝０の場合， ｒ一舞の計算で・での割算が発生し・計算結果か一・ｎ・となる・ 　　　Ｌ “９・°譜　：１・：：：・１：？・ｌ ｉｌｉ’備ｉｌｉ：ｉｌｉｉｉ ０１２３４ 匹四四四罵 ⇒ ０ ０００００ ＾ ↓ｋ ＡＯ脚１　　ＡＯ脚２　　ＡＯβ　　　ＡＯｒ４ 　ユユ　　　ｔ，２　　　ｔ　　　Ａｔ，ｒｉ ｌ　盆・：震：１盆：１ ０　　　Ａ【脚２　　Ａ【β　　Ａ【脚４ ０１２３４ 匹躍匹四罵 Ｄ１２３４ しｂｈｈし ●方程式の順番を入れ替えても解は変わらないので，行を入れ替えて 　Ａｌ，１≠０となるようにすればよい． ●絶対値が０に近い値で割ると，一般に計算精度が悪くなるため，常 　に１行目とＡ１１，Ａ２１，Ａ３、１，Ａ４，１の中で絶対値が最大なものの行との 　入れ替えを行うことにすれば更によい

部分ピボット選択付きガウスの消去法 部分ピボット選択 前進消去の第ｋ列消去において，（Ａｋｋ≠Ｏであっても常に）ｋ行目と 塩，κ，，ＡＮ＿１，ｋの中で絶対値か最大なものの行との入れ替えを行う． 即ち，ｋ行目とｍ≡ｍａｘＦ鳶，，Ｎ＿Ｌｌ．４剥行目を入れ替える． ●例えは，Ａｌ．１，Ａ２、⊥，Ａ３１，Ａ４１がそれぞれＯ，３，－４，２だとしたら，１行 　目と３行目を入れ替えてから，第ｋ列消去を行う． 』：鰐ｉｉ灘ｉｉｉ　ｌｉｉｉ〕〔１〕一〔ｉｉ） 　　　　ｋ行目とｍ行目の入れ替えを行う 勿論Ａの行たけでなく，ｂの対応する要素も入れ替える 列は入れ替えないので，ｆの要素は入れ替えてはならない ●たまたまｋ＝ｍの場合は，入れ替え作業は不要，

配列内の最大要素の探索 要素数Ｎのｄｏｕｂｌｅ型配列の何番目に最大値が入っているかを求めるプロ グラム ●次のように探索する方法が一般的．例文内のコメントをよく理解す 　ること． 　＃ｌｎｃｌｕｄｅ　＜Ｓｔｄｌｏ　ｈ＞ 　＃ｄｅｔｌｎｅ　Ｎ　５ 　１ｎｔ　ｍａｌｎ（ＶＯｌｄ）　｛ 　　ユｎｔ　ｍ７　１ｔ 　　ｄｏｕｂｌｅ　Ａ脚］　＝　｛２，　１，　５，　３，　４｝， 　　皿＝０，　　　　　　　　　　〆＊仮に０番目が最大だとする掌／ 　　ｆｏｒ（ユ＝１，１《Ｎｔ　１＋＋）｛　〆＊１番目から最後まで調べる＊／ 　　ｌｉ（Ａ［１］＞Ａｉｎ］）｛　　　〆＊もし１番目の方か大きいなら＊！ 　　　ｍ旨１／　　　　　　　　〆＊エ番を記憶＊／ 　　｝ 　　｝ 　　ｐｒｌｎｔｆ（１１Ａ［Ｘｄ］　＝　Ｘｆ　Ｉｓ　ｍａｘ、ｎｌｌ」　皿：　Ａ回）． 　　ｒ巳七ｕｒｎ　Ｏ， 　｝

部分ピボット選択付きガウスの消去法 ●プログラムの概形 　　／ホ前進消去ホ／ 　　ｆｏｒ（ｋ＝０，　ｋくＮ－１，　ｋ＋＋）｛ 　　　／＊部分ピボノト選択＊／ 　　　ｍ＝ｋ，　ん仮の最大番号＊／ 　　　ｆｏｒ（ユｉｋ＋ｉ，１く眺１＋＋）｛ 　　　　　省略 　　　｝ 　　　ｐｒｌｎｔｆ（「ｌｃｏ１皿’Ａｄ　ｒｏｗ　Ｚｄ　ｌｓ　皿ａｘ、ｎ１㌧　ｋ．　ｍ）， 　　　ｌｉ（ｍ　ｂ＝ｋ）｛　／＊ｖａ≠ｋの時入れ替え＊！ 　　　　　省略（入れ替えの関数　ｓｖａｐ） 　　　｝ 　　　〆＊第ｋ列｝自去＊ノ 　　　ｆｏｒ　（１　＝　ｋ　＋　１，　１　く　Ｎ，　ユ＋＋）　｛ 　　　　ｆ。ｒ（ｊ＝ｋ，」＜Ｎ，）＋＋）｛ 　　　　　　省略 　　　　｝ 　　　｝ 　　　ＰｒｌｎｔＡｂ（Ａ｝　ｂ）， 　　｝

アルゴリズムの実装手順

本日のまとめ ●連立一次方程式の解法 　　●ガウスの消去法 　　　　●前進消去 　　　　●後退代入 　　●ＮａＮ・ＩＮＦＩＮＩＴＹ 　　・部分ピホット選択付きカウスの消去法