ガウスの消去法
(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も入
れかえなければならない)
• 式の定数倍:行列(と右辺のベクトル)上のある行を定数
倍する
• 2つの式の和:2つの式の和をとり、一方の式をその和と
入れかえる
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
この例はガウスの消去法の2つの基本的段階を示している.
第1段階は前進消去段階で,系統的に変数を処理していき,
対角要素より下の要素がすべてゼロになるよう変換する.
この処理は,三角化とよばれることもある.
第2段階は後退代入段階で,第1段階でえられた三角行列を
使って変数の値を計算する.
5
• 方法の概略
一般に,N個の未知数をもつN個の方程式からなる解く系を解く
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
これらの方程式は行列を使うと1つの式
 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
前進消去は次のようにまとめられる.
はじめに,第1式に適当な数を掛けたものを第1式以外の
各式に対し加えて,第1式以外から第1変数を消去する.
次に,第2式に適当な数を掛けたものを第3式から第N式の
それぞれに対して加えて,第1,2以外の各式から第2変数
を消去する.
次にはじめの3式以外のすべての式から第3変数を消去す
るというようにする.
第i変数を第j式 (i 1  j  N) から消去する時は,第i式を
a ji 倍したものを第j式から引く.
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]はゼロ,またはゼロによる割算が発生した場合は...
その外側のループでi行を適当な(i+1とNの間の)行と入
れかえて,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は,1から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/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の代わりに4つの配列を保持することによっ
てNに関して線形にできる.非ゼロの3つの対角線そ
れぞれに対して1つ,第(N+1)列に対して1つの計4つ
である.このプログラムは最大要素で掃き出していな
いので,ゼロによる割算や計算誤差の累積に対する
保証はない.
3重対角行列のあるものに対しては,その心配がない
ことが証明できる.
15
ダウンロード

ガウスの消去法