第3回連立一次方程式(1) 逆行列 ‐ ガウス‐ジョルダン法

第３回連立一次方程式(１)
逆行列 ‐ ガウス‐ジョルダン法
連立一次方程式
a1,1 x1 + a1, 2 x2 + ... + a1, n xn = b1
a2,1 x1 + a2, 2 x2 + ... + a2, n xn = b2

an,1 x1 + an, 2 x2 + ... + an, n xn = bn
既知
行列
⎡ a1,1 a1, 2  a1, n ⎤ ⎡ x1 ⎤ ⎡ b1 ⎤
⎢a a
⎥ ⎢ x ⎥ ⎢b ⎥
⎢ 2,1 2, 2
⎥⎢ 2 ⎥ = ⎢ 2 ⎥

⎢ 
⎥⎢  ⎥ ⎢  ⎥
⎢a
⎥⎢ ⎥ ⎢ ⎥
a
n, n ⎦ ⎣ xn ⎦
⎣bn ⎦
⎣ n,1
Ax = b
既知
ベクトル
未知
ベクトル
解析的に解けるが n が大きいと時間がかかる
2
基本行演算と基本行列
 基本行演算と基本行列
·  ある行をスカラー倍（０倍は除く）する
i)
ai1
…
ain
i) cai1 …
cain
A
·  この演算は次のような行列を掛けることで実現できる
1
1
i)
0
0
c
1
ai1
…
ain
=
cai1 …
cain
1
Ei(c)
単位行列の第 i 行を c
倍した行列
A
3
基本行演算と基本行列
·  行の交換
ai1
…
ain
i) ak1
…
akn
k) ak1
…
akn
k) ai1
…
ain
i)
A
·  この演算は次のような行列を掛けることで実現できる
1
i)
0
1
ai1
…
ain
k)
1
0
ak1
…
akn
Ei,k
単位行列の第 i 行と第 k
行を交換した行列
A
1
=
i) ak1
…
akn
k) ai1
…
ain
4
基本行演算と基本行列
·  ある行に別の行のスカラー倍を加える
ai1
…
ain
i) ai1+cak1 … ain+cakn
k) ak1
…
akn
k)
i)
ak1
…
akn
A
·  この演算は次のような行列を掛けることで実現できる
1
i)
1
c
ai1
…
ain
k)
0
1
ak1
…
akn
1
Ei,k(c)
=
i) ai1+cak1 … ain+cakn
k)
ak1
…
akn
A
単位行列の第 i 行に第
k 行の c 倍を加えた行
列
5
基本行演算と基本行列
 まとめ
·  行列 A に次の基本行演算を行うことは
 ある行をスカラー倍（０倍は除く）する
 行の交換
 ある行に別の行のスカラー倍を加える
行列 A に次の基本行列を左から掛けることで実現できる．
1
1
i)
0
c
1
1
0
i)
0
1
i)
1
c
1
k)
1
0
k)
0
1
Ei(c)
1
Ei,k
1
1
Ei,k(c)
6
２．逆行列の導出
 行列 A の逆行列 A–1 を計算する． A は正則（|A|≠0）とする．
 すでに分かっていること
逆行列は基本行列の積で表現できる
A–1 = EkEk–1…E2E1
 逆行列計算の基本方針
A–1A = EkEk–1…E2E1 A = I
A
I
E1 A
I を目指して A に基本行演算を 行っていく
E1I
E2E1 A
E2E1I
E3E2E1 A
E3E2E1I
…
…
EkEk–1…E2E1 A = I
同じ基本行演算を I に対して 行っていく
EkEk–1…E2E1I = A–1
7
３．逆行列を求めるアルゴリズム
 例
A
I
1
2
3
1
0
0
1
0
–1
–3
2
0
2
3
4
0
1
0
0
1
2
2
–1
0
3
4
7
0
0
1
0
0
1
0.5
–1
0.5
1
2
3
1
0
0
1
0
0
–2.5
1
0.5
0
–1
–2
–2
1
0
0
1
0
1
1
–1
0
–2
–2
–3
0
1
0
0
1
0.5
–1
0.5
1
2
3
1
0
0
0
1
2
2
–1
0
0
–2
–2
–3
0
1
1
0
–1
–3
2
0
0
1
2
2
–1
0
0
0
2
1
–2
1
1
0
1
k)
0
8
３．逆行列を求めるアルゴリズム
対角要素を１にする
A
1
I
0
1
k)
akk
akk で割る
0
1
0
1
k)
1
0
9
３．逆行列を求めるアルゴリズム
非対角要素を，すべて 0 にする．
A
1 0
aik
1
k)
1
I
aik 倍して引く
0
1
0
1
k)
0
1
0
10
３．逆行列を求めるアルゴリズム注意：数学では行列の行や列の
C言語プログラム
番号は「1」から始まるが，C言語
の配列では「0」から始まる
)
m = 2 * n; [ A | I] をひとつの配列で表し，
for (k=0; k<n; k++){ その要素を
for (j=k+1; j<m; j++) AI[k][j] (k =0,...,n–1, j =0,...,m–1)
AI[k][j] /= AI[k][k]; とする
AI[k][k]=1.0; A
I
for (i=0; i<n ; i++){ 1 0
if (i==k) continue; j
for (j=k+1; j<m; j++) 1
AI[i][j] -= k)
AI[i][k] * AI[k][j]; 0
j
AI[i][k] = 0.0; i)
} } 対角要素を１にする
（基本行演算「行のスカラー倍」を用いる）
)
非対角要素を 0 にする．
（基本行演算「ある行のスカラー倍を別の行に加える」を用いる．）
11
３．逆行列を求めるアルゴリズム
 演習
·  次の行列の逆行列を求めよ
⎛ 2 3 4⎞
⎜
⎟
A = ⎜1 2 3⎟
⎜3 4 7⎟
⎝
⎠
⎛ 2 2 8⎞
⎜
⎟
A = ⎜ 1 2 3⎟
⎜ 2 3 5⎟
⎝
⎠
⎛8 2 2⎞
⎜
⎟
A = ⎜3 1 2⎟
⎜5 2 3⎟
⎝
⎠
A −1
⎛ 1.0 − 2.5 0.5 ⎞
⎜
⎟
= ⎜ 1.0 1.0 − 1.0 ⎟
⎜ − 1.0 0.5 0.5 ⎟
⎝
⎠
A −1
⎛ − 0.25 − 3.5 2.5 ⎞
⎜
⎟
= ⎜ − 0.25 1.5 − 0.5 ⎟
⎜ 0.25 0.5 − 0.5 ⎟
⎝
⎠
A −1
⎛ 0.25 0.5 − 0.5 ⎞
⎜
⎟
= ⎜ − 0.25 − 3.5 2.5 ⎟
⎜ − 0.25 1.5 − 0.5 ⎟
⎝
⎠
次の節の「ピボット選択法」で 対応が可能になる
⎛1 2 3⎞
⎜
⎟
A = ⎜ 2 4 4⎟
⎜3 4 7⎟
⎝
⎠
正則でない． 逆行列が存在しない．
⎛1 2 3⎞
⎜
⎟
A = ⎜ 2 4 4⎟
⎜3 4 7⎟
⎝
⎠
12
４．ガウス‐ジョルダン法（掃出し法による連立一次方程式の解法）
 連立一次方程式 Ax = b を解く．
Ax=b
I を目指して A に基本行演算を 行っていく
E1 A x = E1 b
E2E1 A x = E2E1 b
同じ基本行演算を b に対して 行っていく
E3E2E1 A x = E3E2E1 b
…
EkEk–1…E2E1 A x = EkEk–1…E2E1 b
I
x = EkEk–1…E2E1 b
13
４．ガウス‐ジョルダン法（掃出し法による連立一次方程式の解法）
A
 例
b
1
2
3
1
1
0
–1
–1
2
3
4
1
0
1
2
1
3
4
7
1
0
0
1
0
1
2
3
1
1
0
0
–1
0
–1
–2
–1
0
1
0
1
0
–2
–2
–2
0
0
1
0
1
2
3
1
0
1
2
1
0
–2
–2
–2
1
0
–1
–1
0
1
2
1
0
0
2
0
x
14
４．ガウス‐ジョルダン法（掃出し法による連立一次方程式の解法）
C言語プログラム
逆行列の場合とは
ここが違うだけ
[ A | b] をひとつの配列で表し，
その要素を
Ab[k][j] (k =0,...,n–1, j =0,...,m–1)
とする
A
1
0
1
k)
0
b
j
)
j
i)
)
m = n + 1; for (k=0; k<n; k++){ for (j=k+1; j<m; j++) Ab[k][j] /= Ab[k][k]; Ab[k][k] = 1.0; for (i=0; i<n ; i++){ if (i==k) continue; for (j=k+1; j<m; j++) Ab[i][j] -= Ab[i][k] * Ab[k][j]; Ab[i][k] = 0.0; } } 対角要素を１にする
（基本行演算「行のスカラー倍」を用いる）
非対角要素を 0 にする．
（基本行演算「ある行のスカラー倍を別の行に加える」を用いる．）
15
４．ガウス‐ジョルダン法（掃出し法による連立一次方程式の解法）
 注意
·  計算の手間を考慮するだけならば，ガウス消去法等を用いるべき
 連立一次方程式を解くだけなら，逆行列を求める必要はない
 逆行列を求める方法と同じコストのガウス‐ジョルダン法も必要ない
16
４．ガウス‐ジョルダン法（掃出し法による連立一次方程式の解法）
 演習
·  次の連立方程式を解け
⎛ 2 3 3⎞
⎛1⎞
⎜
⎟
⎜ ⎟
⎜ 1 2 3⎟ x = ⎜ 2 ⎟
⎜ 3 2 1⎟
⎜ 2⎟
⎝
⎠
⎝ ⎠
⎛ 2⎞
⎜ ⎟
x = ⎜ − 3⎟
⎜ 2⎟
⎝ ⎠
⎛2 2 1⎞
⎛1⎞
⎜
⎟
⎜ ⎟
⎜ 3 2 3⎟ x = ⎜1⎟
⎜ 2 3 2⎟
⎜ 3⎟
⎝
⎠
⎝ ⎠
⎛ − 1. 2 ⎞
⎜
⎟
x = ⎜ 1. 4 ⎟
⎜ 0. 6 ⎟
⎝
⎠
17
５．ピボット選択法
 前述のアルゴリズムの問題点
·  ak,k = 0 ⇒ゼロ除算が起こる
ak,j := ak,j/ak,k ⇒ak,j の持つ誤差の拡大
·  ak,k が非常に小さい ⇒
⇒ak,j 自身の値も大きくなる
ai,j := ai,j – ai,k ak,j ⇒情報落ち（ak,jが非常に
大きいため）
1
ピボット
0
1
k)
ai,k
ai,j
ak,k
ak,j
0
18
５．ピボット選択法
 ピボット選択法（部分選択法）
·  基本行演算のうち「行の交換」はまだ使っていなかった．
Ekˆ,k : kˆ = arg max | ai ,k |
i = k ,…, n
1
k)
k̂ )
0
1
0
19
５．ピボット選択法
for (k=0; k<n; k++){ for (j=k+1; j<m; j++) Ab[k][j] /= Ab[k][k]; Ab[k][k] = 1.0; for (i=0; i<n ; i++){ if (i==k) continue; for (j=k+1; j<m; j++) Ab[i][j] -= Ab[i][k] * Ab[k][j]; Ab[i][k] = 0.0; } } k̂ を求め，k行とk̂ 行を交換する
プログラムをここに追加する
1
k)
k̂ )
0
1
0
20

Download Report