ガウスの消去法
（Gaussian elimination）
1
•ガウスの消去法とは…
連立方程式系を解く基本的なアルゴリズム
• 簡単な例
x 3y 4z 8
x y 2z 2
x 2 y 5 z 1
これらを方程式を同時に満たす変数の値を求める．
場合によっては，解がないことや，いくつもあることもある．
(例えば x y 1, x y 2 のように2つの方程式が矛盾する時や，
2つの方程式が同じとか，方程式の数より変数の個数が多い時)
方程式の個数と変数の個数は等しいものとし，解があるときは一意な
解を求める．
2
方程式に対して解を変えないで実行出来る
操作
• 方程式の入れかえ：明らかに，方程式を書き下す順番は
解に影響しない行列の表現では，この操作は行列(と右
辺のベクトル)の行の入れかえに対応する．
• 変数名のつけかえ：この操作は，行列の上では列の入れ
かえに対応する(列iとjを入れかえた時は変数 xiとxjも入
れかえなければならない)
• 式の定数倍：行列(と右辺のベクトル)上のある行を定数
倍する
• ２つの式の和：２つの式の和をとり、一方の式をその和と
入れかえる
xji
3
例
x 3y 4z 8
3つ以上の変数の場
合も扱いやすいように
x y 2z 2
x 2 y 5 z 1
添字を使って変数名
をつける
x1 3x 2 x3 8
x1 x 2 2 x3 2
x1 2 x 2 5 x3 1
3 4 x1 8
1
3 4 x1 8
1
2式=1式-2式
2 2 x 2 6
0
1 2 x 2 2
1
1 2 5 x3 1
1 2 5 x3 1
1 3 4 x1 8
3式=2式-2×3式
3式=1式＋3式
0 2 2 x 2 6
0 1 1 x 3 7
1 3 4 x1 8
0 2 2 x 2 6
0 0 4 x3 8
4
x1 3x 2 4 x3 8
2 x 2 2 x3 6
4 x3 8
2式に代入
2 x2 4 6
x2 5
3式より
x3 2
1式に代入
x1 15 8 8
x1 1
この例はガウスの消去法の２つの基本的段階を示している．
第１段階は前進消去段階で，系統的に変数を処理していき，
対角要素より下の要素がすべてゼロになるよう変換する．
この処理は，三角化とよばれることもある．
第２段階は後退代入段階で，第１段階でえられた三角行列を
使って変数の値を計算する．
5
• 方法の概略
一般に，Ｎ個の未知数をもつＮ個の方程式からなる解く系を解く
a11 x1 a12 x 2 a1NxN b1
a 21 x 2 a 22 x 2 a 2 NxN b2
aN 1x1 aN 2 x 2 aNNxN bN
これらの方程式は行列を使うと１つの式
a11 a12 a1N x1 b1
a
21
a
22
a
2
N
x 2 b 2
aN 1 aN 2 aNN xN bN
すなわち， Ax b のように書ける．ただし A は係数行列，x は変数，
b
は方程式の右辺を表す．
の要素は
Aの行と一緒に扱われるので，
b
( N列と
1) みなし
を b の第
A
N (Nの行列を用いて両者を保
1)
持すると便利である．
6
前進消去は次のようにまとめられる．
はじめに，第１式に適当な数を掛けたものを第１式以外の
各式に対し加えて，第１式以外から第１変数を消去する．
次に，第２式に適当な数を掛けたものを第３式から第Ｎ式の
それぞれに対して加えて，第１，２以外の各式から第２変数
を消去する．
次にはじめの３式以外のすべての式から第３変数を消去す
るというようにする．
第ｉ変数を第ｊ式 (i 1 j N) から消去する時は，第ｉ式を
a ji 倍したものを第ｊ式から引く．
a ii
この処理は次のコードのようにもっと簡潔に書ける．
7
for (i=1;i<=N;i++)
3
#計算時間は N に比例
for (j=i+1;j<=N;j++)
for (k=N+1;k>=i;k--)
a[j][k] -=a[i][k]* a[j][i]/a[i][i];
#同じ行の他の要素を計算する前 に
# a[i][j]を壊さないよう逆方向に進む
• a[i][i]はゼロ，またはゼロによる割算が発生した場合は．．．
その外側のループでｉ行を適当な(ｉ＋１とＮの間の)行と入
れかえて，a[i][i]がゼロでないようにすれば解決
(そのような行が見つからないときには，その行列は特異．す
なわち一意な解が存在しない)
第i列の対角要素より下の要素を消去するのに使われる要素
a[i][i]を枢軸とよぶ．
8
実際には，第i列で非ゼロ要素をもつ行を探すだけでなく，
第i列で絶対値が最大の要素をもつ行を使うのがよい．
なぜなら，行の定数倍に使われる枢軸要素が非常に小さい
と，大きな計算誤差が生じることがあるからである．
a[i][i]が非常に小さいと，第i変数を第j式から消去する時に
使われるスケーリング因子a[j][i]/ a[i][i]が非常に大きくなる．
実際，それは係数が小さく見えるほどまでに大きくなり，
a[j][i]の値が丸め差のために歪んでしまうこともある．
例：
0.001 x1 －6 x 2 =－6.001
3 x1 +5 x 2 =2
9
eliminate()
substitute()
{ int i,j,k,max;
{ int j,k;
float t;
float t;
for (i=1;i<=N;i++)
for (j=N;j>=1;j--)
{ max=i;
for(j=i+1;j<=N;j++)
{ t=0.0;
if (abs(a[j][i])>abs(a[max][i]))max=j;
for (k=j+1;k<=N;k++) t
for(k=i;k<=N+1;k++)
+=a[j][k]*x[k];
{t=a[i][k];
x[j]=(a[j][N+1]-t)/a[j][j];}
a[i][k]=a[max][k]=;
}
a[max][k]=t;}
for(j=i+1;j<=N;j++)
for(k=N+1;k>=i;k--)
a[j][k]-=a[i][k]*a[j][i]/a[i][i];}
}
10
• eliminateは，１からNまでの各iに対して，i列を下に向
かって調べ最大の要素を見つける．その要素を含む行と
I行と入れかえた後に，前と同様にi+1からNの式から第i
変数を消去するというプログラム．
(前のコードに，a[i][i]がゼロの時の処理を加えたもの)
• 枢軸要素a[i][i]を使って第i式以外のすべての式から
（(i+1)式からN式までだけでなく）第i変数を消去するアル
ゴリズムもある．この処理は完全掃き出し（ガウスジョル
ダンの消去法）とよばれる．前進消去では，この一部分を
実行するだけなので部分掃き出しとよばれる．
• substituteは後退代入段階のプログラム
11
性質37.1 N元の連立方程式は N 3 /3回の乗算と加算で解く
ことができる
– substituteの実行時間はO( N 2)なので，ほとんどの仕事
はeliminateで実行される．それを調べてみると，各iに
対してkに関するループがN－i+2回繰り返され，jに関
するループがN－i回繰り返されることがわかる．
これは内側のループが 1 i N 1 (N－i+2) (N－i)=
N 3/３+O(N 2)回実行されることを意味する．
a[j][i]/a[i][i]の値はkに関するループの外で計算される
ので，内側のループは1回の乗算と1回の加算からな
る．
12
• 変形と拡張
今までの方法は，ほとんどの要素がゼロでないN×Nの
行列に適している．要素のほとんどがゼロである疎な行列
に向いた特別な技法がある．
・“帯”行列
– すべての非ゼロ要素が対角部分に近いところにある
行列
この場合，消去法の内側のループの繰り返しがほん
3
の少しなので，全計算時間は N ではなくNに比例する．
・帯行列の中で特別な例：
“三重対角”行列
– 対角上と，対角のすぐ上かすぐ下の要素だけがゼロ
でない行列
13
a11
例 N=5の場合の
三重対角行列の一般形 a 21
0
0
0
for (i=1;i<N;i++)
{ a[i+1][N+1] -=a[i][N+1]*a[i+1][i]/a[i][i];
a[i+1][i+1] -=a[i][i+1]*a[i+1][i]/a[i][i];
}
for (j=N;j>=1;j--)
x[j]=(a[j][N+1]-a[j][j+1]*x[j+1])/a[j][j];
a12
0
0
a 22
a32
a 23
a33
0
a34
0
a 43
a 44
0
0
a54
0
0
0
a 45
a55
• このような行列に対しては前進消去と後退代入は単一の
forループになる．前進消去ではa[i][k]=0(k>i+1)なので
j=i+1とk=i+1の場合だけが含まれる
14
性質37.2 3重対角になっている連立方程式系は線形時間
で解ける
– もちろん3重対角行列に対してはサイズN×Nの2次元
配列は用いない．上記のプログラムに必要な記憶量
は配列aの代わりに４つの配列を保持することによっ
てNに関して線形にできる．非ゼロの3つの対角線そ
れぞれに対して1つ，第(N+1)列に対して1つの計4つ
である．このプログラムは最大要素で掃き出していな
いので，ゼロによる割算や計算誤差の累積に対する
保証はない．
3重対角行列のあるものに対しては，その心配がない
ことが証明できる．
15
© Copyright 2024 ExpyDoc
