部品の在庫の状況からいかに最小のコストで製品をつくるか」という問題，機械要素の運動の 問題，電気回路の解析の問題など，いくつか挙げられる．つまり，計算機で多元連立方程式を 解くことができれば，幅広く工学の諸問題の解決に役立つ． ここではガウスの消去法で多元連立 1 次方程式を解く方法を学ぶ．例題としてブリッジ T 型 とホーイストンブリッジの電気回路を扱う．
12.1.1 ガウスの消去法による多元連立 1 次方程式の解法
ガウスの消去法は 2 つの手順で進められる．初めに前進消去（forward elimination）を行い， 次に後退挿入（back substitution）を行う．まず計算手順を具体的な例によって把握する．次に それらの処理の内容を式で示す． 12.1.1.1 前進消去と後退挿入の手順 多元連立 1 次方程式は四則演算（加減乗除）と行や列の入替えを行って解ける．前進消去で は行列の左下の三角の部分を 0 にする演算を行う．後退挿入では行列の下から解を求める演算 を行う． 例題 1 次の 3 元連立 1 次方程式をガウスの消去法で解く． { 𝑥4𝑥11+ 2𝑥+ 5𝑥22+ 3𝑥+ 6𝑥33= 14= 32 7𝑥1+ 8𝑥2 = 23 前進消去 後退挿入 答えは x1=1，x2=2，x3=3 である． 回数 番号 x1 x2 x3 b 演算内容 0 ① ② ③ 1 2 3 4 5 6 7 8 0 14 32 23 1 ④ ⑤ ⑥ 1 2 3 0 -3 -6 0 -6 -21 14 -24 -75 ②－①×4 ③－①×7 2 ⑦ ⑧ ⑨ 1 2 3 0 -3 -6 0 0 -9 14 -24 -27 ⑥－⑤×2 演 算 −9𝑥3= −27 → 𝑥3= −27 −9 = 3 −3𝑥2− 6(3) = −24 → 𝑥2= −24 + 18 −3 = 2 𝑥1+ 2(2) + 3(3) = 14 → 𝑥1= 14 − 4 − 9 = 1しかしながら式の係数によっては，前進消去ができない場合がある． 例題 2 （前進消去ができない場合） 次の 3 元連立 1 次方程式をガウスの消去法で解く． { 𝑥4𝑥11+ 2𝑥+ 8𝑥22+ 3𝑥+ 6𝑥33= 14= 38 7𝑥1+ 8𝑥2 = 23 前進消去を行っていくと，⑤のところで x2の係数が 0 となり，計算を続けることができなくな る．そのような場合では，⑤と⑥行を入れ換える（つまり⑧と⑨）．この操作は方程式の上下 の位置を変える操作であるので，結果に影響しない． 前進消去 後退挿入 回数 番号 x1 x2 x3 b 操 作 0 ① ② ③ 1 2 3 4 8 6 7 8 0 14 38 23 1 ④ ⑤ ⑥ 1 2 3 0 0 -6 0 -6 -21 14 -18 -75 ②－①×4 ③－①×7 2 ⑦ ⑧ ⑨ 1 2 3 0 -6 -21 0 0 -6 14 -75 -18 ⑥との交換 ⑤との交換 答えは x1=1，x2=2，x3=3 である． ⑤のように，対角要素 ai i (i=1, 2, ...) が 0 のとき，それ以降の前進消去はできなくなる．また， 対角要素が 0 に近い非常に小さい値のときは丸め誤差が大きくなる．それらを避けるためには 対角要素の選択，つまりピボットティング（pivoting）が行われる．具体的には ai j (i = j, j + 1, ..., n ) の最大行を探して，その行と i 行を入れ換える．この様にすれば前進消去を進めることがで きる．更に精度を上げるためには列方向にも同様のことを行う．行および列の両方向について のピボッティングを行うことを完全ピボッティングという． 計 算 −6𝑥3= −18 → 𝑥3= −18 −6 = 3 −6𝑥2− 21(3) = −75 → 𝑥2= −75 + 63 −6 = 2 𝑥1+ 2(2) + 3(3) = 14 → 𝑥1= 14 − 4 − 9 = 1
12.1.1.2 ガウスの消去法の数式による表現 n 元連立 1 次方程式を次式で示す．（xn以外は定数である．） [ 𝑎11 𝑎12 ⋯ 𝑎1𝑛 𝑎21 ⋮ ⋱ 𝑎2𝑛⋮ 𝑎𝑛1 𝑎𝑛2 ⋯ 𝑎𝑛𝑛 ] [ 𝑥1 𝑥2 ⋮ 𝑥𝑛 ] = [ 𝑏1 𝑏2 ⋮ 𝑏𝑛 ] …(1) 要素 a11, a12, …, annから構成される n × n 行列に，要素 b1, b2 …, bnから構成される列を付加し た次の行列について考える． [ 𝑎11 𝑎12 ⋯ 𝑎1𝑛 𝑏1 ⋮ ⋱ ⋮ 𝑎𝑛1 𝑎𝑛2 ⋯ 𝑎𝑛𝑛 𝑏𝑛 ] …(2) 手順 1： 前進消去 各要素に対し，次の演算を行う． 𝑎𝑖𝑗 = 𝑎𝑖𝑗− 𝑎𝑖𝑚 𝑎𝑚𝑗 𝑎𝑚𝑚 …(3) ただし，m = 1, 2, …, n - 1, i = m + 1, m + 2, …, n, j = m, m + 1, m + 2, …, n + 1．一通りの計算 が終わると，次式を得る．ただし 2 行以降の要素は値が変わるので，ダッシュ（’）がつく． [ 𝑎11 𝑎12 𝑎13 𝑎14 𝑎15 𝑎16 ⋯ 𝑎1𝑛 𝑏1 0 𝑎22′ 𝑎23 ′ 𝑎24 ′ 𝑎25 ′ 𝑎26 ′ ⋯ 𝑎2𝑛′ 𝑏2′ 0 0 𝑎33 ′ 𝑎34 ′ 𝑎35 ′ 𝑎36 ′ ⋯ 𝑎3𝑛′ 𝑏3′ 0 0 0 𝑎44 ′ 𝑎45 ′ 𝑎46 ′ ⋯ 𝑎4𝑛′ 𝑏4′ ⋮ 0 0 0 0 0 … 0 𝑎𝑛𝑛′ 𝑏𝑛′ ] …(4) 式(4)は下の式(5)の意味を持っている． [ 𝑎11 𝑎12 𝑎13 𝑎14 𝑎15 𝑎16 ⋯ 𝑎1𝑛 0 𝑎22′ 𝑎23 ′ 𝑎24 ′ 𝑎25 ′ 𝑎26 ′ ⋯ 𝑎2𝑛′ 0 0 𝑎33 ′ 𝑎34 ′ 𝑎35 ′ 𝑎36 ′ ⋯ 𝑎3𝑛′ 0 0 0 𝑎44 ′ 𝑎45 ′ 𝑎46 ′ ⋯ 𝑎4𝑛′ ⋮ 0 0 0 0 0 … 0 𝑎𝑛𝑛′ ][ 𝑥1 𝑥2 ⋮ 𝑥𝑛 ] = [ 𝑏1 𝑏2′ ⋮ 𝑏𝑛′ ] …(5) 手順 2： 後退挿入 式(5)の n 行に注目する．𝑎𝑛𝑛′ 𝑥𝑛= 𝑏𝑛′の関係から𝑥𝑛が得られる．次に n - 1 行では， 𝑎𝑛−1 𝑛−1′ 𝑥𝑛−1+ 𝑎𝑛𝑛′ 𝑥𝑛= 𝑏𝑛−1′ の関係から𝑥𝑛−1が得られる．同様に下から順に繰り返すことによっ て，x1 から xnが得られる．これらの演算は次のように表現できる． 𝑥𝑛=𝑎𝑛 𝑛+1𝑎 𝑛 𝑛 …(6) 𝑥𝑖= 𝑎𝑖 𝑛+1−∑𝑛−𝑖𝑗=1𝑎𝑖 𝑖+𝑗 𝑥𝑖+𝑗 𝑎𝑖𝑗 …(7) ただし，i = n - 1, …, 2, 1，ann + 1 = b nである．
12.1.2
プログラムの処理手順の確認とその実行 プログラム（付録参照）の処理内容を確認する．プログラムには例題 1 についてのデータが セットされている．次にプログラムのコメント行，make procedure に従ってコンパイルを行い， プログラムを実行する．ガウスの消去法の演算過程と解はテキスト画面に表示される．出力さ れた結果と例題 1 の結果とを比較して計算が正しいことを確認する．出力について以下に説明 する． 0 回目 x1 x2 x3 (ピボッティング：行方向：1 行と 3 行の入れ替えをする) 1 2 3 14 4 5 6 32 7 8 0 23 1 回目 x1 x2 x3 (ピボッティング：列方向：1 列と 2 列の入れ替えをする) 7 8 0 23 4 5 6 32 1 2 3 14 2 回目 x2 x1 x3 (2，3 行について前進消去をする) 8 7 0 23 5 4 6 32 2 1 3 14 3 回目 x2 x1 x3 (ピボッティング：列方向：2 列と 3 列の入れ替えをする) 8 7 0 23 0 -0.375 6 17.625 0 -0.75 3 8.25 4 回目 x2 x3 x1 (3 行について前進消去をする) 8 0 7 23 0 6 -0.375 17.625 0 3 -0.75 8.25 5 回目 x2 x3 x1 (前進消去終わり) 8 0 7 23 0 6 -0.375 17.625 0 0 -0.563 -0.563実行結果 連立方程式を解く際には解の種類が一意に決まるか否かあらかじめ調査するとよい．前述の例 題 1 と 2 の係数の rank はそれぞれ 3 であり，x の数 n と同じであった．つまり 1 種類の解が存 在することになる．連立方程式の係数によっては rank<n となることもある．そのような場合で は無数の解が存在する．その様な問題についてこのプログラムを実行すると，そのうちの一つ だけの解が求まることになる．例えば次の問題を考える． { 𝑥1+ 2𝑥2+ 3𝑥3= 14 4𝑥1+ 5𝑥2+ 6𝑥3= 32 7𝑥1+ 8𝑥2+ 9𝑥3 = 23 係数部分をまとめると次のようになる． (1 2 34 5 6 7 8 9 ) この rank は 2 であり，解は無数にある．例えば次のようなものがある． x1=1，x2=2，x3=3 x1=2，x2=0，x3=4 x1=10，x2=-16，x3=12 従って，n 元連立 1 次方程式の問題を解く際には rank を求めておくとよい． 次に例題 2 についても計算してみよ．（プログラムの編集については付録を参照のこと．） このプログラムが完全ピボッティングに対応していることを確認せよ． レポート C n 元連立 1 次方程式で表現できる問題を見つけ，プログラムを使って解を求めよ．問題は自 身で行っている研究の内容や工学的内容が望ましい．レポートの書式はレポート A に従うこと． 参考文献 1) 佐藤次男ほか，C による理工学問題の解法，日刊工業新聞社，ISBN4-526-0363203.
付録 ガウスの消去法のプログラム
黄色い部分にセットされているデータは例題 1 のものである．他の問題を計算するときには 黄色い部分を変更する．
/* Gauss elimination method */
// make procedure: gcc 12_GEM.c -o 12_GEM -lm #include <stdio.h>
#include <math.h>
#define NUMBER_OF_DATA 3 //連立方程式の解（x1,x2,x3,…）の個数
#define NUNBER_OF_CF 11 // 扱う解の数の上限の数
static double data[NUNBER_OF_CF][NUNBER_OF_CF+1] = {
//ここでは例題 1 のデータをセットする {1, 2, 3, 14}, /* a11, a12, ..., a1n, b1 */
{4, 5, 6, 32}, /* a21, a22, ..., a2n, b2 */ {7, 8, 0, 23} /* an1, an2, ..., ann, bn */ }; int GEM(); //ガウスの消去法を行う関数のプロトタイプ宣言 void show_procedure(); //配列変数の内容を画面へ表示するための関数のプロトタイプ宣言 int n = (int)NUMBER_OF_DATA; //未知変数（x1,x2,x3）の個数 double epsilon=1e-18; //ピボット、対角要素の下限値をセットする double x[NUNBER_OF_CF]; //連立方程式の解（ x1, x2, ..., xn）がセットされる配列 int main() { int i; int n=(int)NUMBER_OF_DATA; // 連立方程式の解、つまり x の個数をセットする int pivot; //ピボットの状態を示すフラグ（0：対角要素が 0 に近くて計算不可，1：計算可） pivot = GEM(); //ガウスの消去法を実行し、ピボットの状態を戻り値として得る if(pivot == 1) //1：計算結果（解 x1、x2,…）を表示、1!=：計算不可の表示
for(i=0 ;i<n; i++) { printf("x %d = %lf \n", i+1, x[i]); }
else { printf("This process was stopped, because a pivot is to small.\n"); } //対角要素が 0 に近接した場合には計算終了となる
i=getchar(); //計算終了後に直ぐ画面が消えないようにキー入力をする
return 0; //プログラムを正常終了する
}
// Gauss Elimination Method // 連立方程式の要素と変数の関係 // AX=B
//
// A = | a11, a12, ..., a1n |, X = | x1 |, B = | b1 | // | a21, a22, ..., a2n | | x2 | | b2 | // | ... | | ... | | ... | // | an1, an2, ..., ann | | xn | | bn | //
// data = | a11, a12, ..., a1n, b1 |, data[0][0] =a11 data[行][列] // | a21, a22, ..., a2n, b2 |, data[1][0] =a21
// | ... |
// | an1, an2, ..., ann, bn |, data[n-1][n-1] =ann //
//下の3 つの行は原理で述べられた数式の添え字とプログラム中の変数の添え字との対応をとるためのマクロ #define A(x, y) data[x-1][y-1]
#define X(z) x[z-1] #define LN(z) ln[z-1] int GEM() //ガウスの消去法を実行する副関数 { double max; //ピボッティングのとき、要素の最大値を記憶するための変数 double tmp, tmp1; //テンポラリ変数 int ln[NUMBER_OF_DATA+1]; //列の入れ替えの情報を保存する配列変数
int column, i, itmp, ipj, j, m, nm1, l, row; //row:行
show_procedure(); //配列変数の内容を画面へ表示する
for(i=1; i<=n; i++) { LN(i)=i; } //入れ替え前の解 x の添え字（i=1,2,…）をセット
for(m=1; m<=n-1; m++) { //前進消去の開始
max = 0.0; //要素の最大値を記憶するための変数をリセット
for(i=m; i<=n; i++) {
for(j=m; j<=n; j++) { //要素が最大のものを探索する：ピボッティング
if(max >= fabs(A(i, j))) { continue; } //con.：}}の前へ
max = fabs(A(i, j)); //最大値を記録する
row = i; //最大の要素の行と列の位置を記録する
column = j; // row：行、column：列
}}
if(max <= epsilon) { return 0; } //もし最大のものが 0 に近ければ演算停止
if(m != row) { //m==rowは行の最後を意味する
for(l=m; l<=n+1; l++) {
tmp = A(row,l); //行方向のピボッティング.行の入れ替え
A(row, l) = A(m, l); A(m, l) = tmp; }
printf("Coefficient and constant were changed in row-direction: %d
<-> %d lines.\n", m, row); //入れ替えの情報を画面へ表示する
show_procedure(); //配列変数の内容を画面へ表示する
} /*もし m==column：列の最後*/
if(m != column) { //m==columnは列の最後を意味する
for(l=1; l<=n; l++) { //列方向のピボッティング.列の入れ替え
tmp = A(l, column); A(l, column) = A(l, m); A(l, m) = tmp;
}
printf("Coefficient and constant were changed in column-direction: %d
<-> %d lines.\n", m, column); //入れ替えの情報を画面へ表示する show_procedure(); //配列変数の内容を画面へ表示する //後の後退挿入で解(x1,x2,…)の存在する列の情報が必要になるため tmp1 = LN(m); //列の入れ替えの位置を LN に記録する LN(m) = LN(column); LN(column) = tmp1; }
for(i=m+1; i<=n; i++) { //前進消去（ガウスの消去法）による計算
for(j=n+1; j>=m; j--) {
A(i, j) -= A(i, m) * A(m, j) / A(m, m); //式(3)の計算
printf("Coefficient and constant after the forward elimination, %d line below.\n", m+1);
show_procedure(); //配列変数の内容を画面へ表示する
}
if(fabs(A(n, n)) <= epsilon) { return 0; } //もし Ann が 0 に近ければ演算打ち切り
//後退挿入の開始
X(n) = A(n, n+1) / A(n, n); //後退挿入の計算,式(6), xn is set.
itmp = LN(n);
A(n, itmp) = X(n); // a[n-1][itemp-1]=xn (=LN(n) )
//もし左端からガウスの消去法を行うと、左端の係数が0となり、計算できない //つまり右端から計算するために、i--とする
for(i=n-1; i>=1; i--) {
nm1 = n - i; // nm1 means n minus 1:(n-1)
X(i) = 0.0;
//式(7)の計算
for(j=1; j<=nm1; j++) { // sigma(aii+j*xi+j)
ipj = i + j; // ipj means i plus j:(i+j)
X(i) += A(i, ipj) * X(ipj); // X(i) = A(i, ipj) * X(ipj) + X(i);
}
X(i) = (A(i, n+1) - X(i)) / A(i, i); /* x[i] <-> xi */ itmp = LN(i);
A(n, itmp) = X(i); }
for(i=1; i<=n; ++i) { X(i) = A(n, i); } //入れ替えを元に戻す
return 1; }
void show_procedure() //配列変数の内容を画面へ表示する関数
{
int i, j;
for(i=0; i<n; i++) {
for(j=0; j<n+1; j++) { printf("%8.3lf", data[i][j]); } printf("\n"); } printf("\n"); } メモ： プログラムで変数 LN は，列の入れ替えのピボッティングのとき，列の入れ替えの位置を記 憶するために使われる．列の入替では解の位置が変わってしまう．具体的例を以下に示す． 例： x1 x2 x3 … xn (x2と x3の列を入れ替えると) x1 x3 x2 … xn この状態では後退消去で連立方程式の解（x1, x2, …）を求めることはできない．つまり列の入 れ替えのピボッティングが終わった後に，また列の要素の位置を交換し，その後，後退消去を 計算することになる．そのためピボッティング以前の解の位置を変数 LN に保存する．
