数値計算法

Transcription

1 12.1 電気回路網に関するキルヒホッフの法則による解法 1 工学的諸問題を多元連立 1 次方程式で表現することができる. 例えば, 荷物を最短の時間と最低のコストで輸送するためにはどのようなルートで物流を行うか という問題, 工場の部品の在庫の状況からいかに最小のコストで製品をつくるか という問題, 機械要素の運動の問題, 電気回路の解析の問題など, いくつか挙げられる. つまり, 計算機で多元連立方程式を解くことができれば, 幅広く工学の諸問題の解決に役立つ. ここではガウスの消去法で多元連立 1 次方程式を解く方法を学ぶ. 例題としてブリッジ T 型とホーイストンブリッジの電気回路を扱う ガウスの消去法による多元連立 1 次方程式の解法ガウスの消去法は 2 つの手順で進められる. 初めに前進消去 (forward elimination) を行い, 次に後退挿入 (back substitution) を行う. まず計算手順を具体的な例によって把握する. 次にそれらの処理の内容を式で示す 前進消去と後退挿入の手順多元連立 1 次方程式は四則演算 ( 加減乗除 ) と行や列の入替えを行って解ける. 前進消去では行列の左下の三角の部分を 0 にする演算を行う. 後退挿入では行列の下から解を求める演算を行う. 例題 1 次の 3 元連立 1 次方程式をガウスの消去法で解く. x 1 + 2x 2 + 3x 3 = 14 { 4x 1 + 5x 2 + 6x 3 = 32 7x 1 + 8x 2 = 前進消去 回数 番号 x1 x2 x3 b 演算内容 後退挿入 演 算 9x 3 = 27 x 3 = 27 9 = 3 3x 2 6(3) = x 2 = = 2 3 x 1 + 2(2) + 3(3) = 14 x 1 = = 1 答えは x 1=1,x 2=2,x 3=3 である.

2 しかしながら式の係数によっては, 前進消去ができない場合がある. 例題 2 ( 前進消去ができない場合 ) 次の 3 元連立 1 次方程式をガウスの消去法で解く. x 1 + 2x 2 + 3x 3 = 14 { 4x 1 + 8x 2 + 6x 3 = 38 7x 1 + 8x 2 = 23 前進消去を行っていくと,5のところで x 2 の係数が 0 となり, 計算を続けることができなくなる. そのような場合では,5と6 行を入れ換える ( つまり8と9). この操作は方程式の上下の位置を変える操作であるので, 結果に影響しない 前進消去 回数 番号 x1 x2 x3 b 操作 との交換 5 との交換 後退挿入 計 算 6x 3 = 18 x 3 = 18 6 = 3 6x 2 21(3) = x 2 = = 2 6 x 1 + 2(2) + 3(3) = 14 x 1 = = 1 答えは x 1=1,x 2=2,x 3=3 である. 5のように, 対角要素 a i i (i=1, 2,...) が 0 のとき, それ以降の前進消去はできなくなる. また, 対角要素が 0 に近い非常に小さい値のときは丸め誤差が大きくなる. それらを避けるためには対角要素の選択, つまりピボットティング (pivoting) が行われる. 具体的には a i j (i = j, j + 1,..., n ) の最大行を探して, その行と i 行を入れ換える. この様にすれば前進消去を進めることができる. 更に精度を上げるためには列方向にも同様のことを行う. 行および列の両方向についてのピボッティングを行うことを完全ピボッティングという.

3 ガウスの消去法の数式による表現 n 元連立 1 次方程式を次式で示す.(x n 以外は定数である.) a 11 a 12 a x 1n 1 b 1 [ a21 x a2n 2 b ] [ ] = [ 2 ] a n1 a n2 a nn x n b n (1) 要素 a 11, a 12,, a nn から構成される n n 行列に, 要素 b 1, b 2, b n から構成される列を付加し た次の行列について考える. a 11 a 12 a 1n b 1 [ a n1 a n2 a nn b n ] (2) 手順 1: 前進消去 各要素に対し, 次の演算を行う. a ij = a ij a im a mj a mm (3) ただし,m = 1, 2,, n - 1, i = m + 1, m + 2,, n, j = m, m + 1, m + 2,, n + 1. 一通りの計算 が終わると, 次式を得る. ただし 2 行以降の要素は値が変わるので, ダッシュ ( ) がつく. [ a 11 a 12 a 13 a 14 a 15 a 16 a 1n b 1 0 a 22 a 23 a 24 a 25 a 26 a 2n b a 33 a 34 a 35 a 36 a 3n b a 44 a 4n b 4 a 45 a a nn b n ] (4) 式 (4) は下の式 (5) の意味を持っている. [ a 11 a 12 a 13 a 14 a 15 a 16 a 1n 0 a 22 a 23 a 24 a 25 a 26 a 2n 0 0 a 33 a 34 a 35 a 36 a 3n a 44 a 4n a 45 a 46 x 1 x 2 x n] a nn ] [ = [ b 1 b 2 b n ] (5) 手順 2: 後退挿入 式 (5) の n 行に注目する.a nn x n = b n の関係からx n が得られる. 次に n - 1 行では, a n 1 n 1 x n 1 + a nn x n = b n 1 の関係からx n 1 が得られる. 同様に下から順に繰り返すことによっ て,x 1 から x n が得られる. これらの演算は次のように表現できる. x n = a n n+1 a n n x i = a i n+1 n i j=1 a i i+j x i+j a ij (6) (7) ただし,i = n - 1,, 2, 1,a n n + 1 = b n である.

4 プログラムの処理手順の確認とその実行プログラム ( 付録参照 ) の処理内容を確認する. プログラムには例題 1 についてのデータがセットされている. 次にプログラムのコメント行,make procedure に従ってコンパイルを行い, プログラムを実行する. ガウスの消去法の演算過程と解はテキスト画面に表示される. 出力された結果と例題 1 の結果とを比較して計算が正しいことを確認する. 出力について以下に説明する. 0 回目 x 1 x 2 x 3 ( ピボッティング : 行方向 :1 行と 3 行の入れ替えをする ) 回目 x 1 x 2 x 3 ( ピボッティング : 列方向 :1 列と 2 列の入れ替えをする ) 回目 x 2 x 1 x 3 (2,3 行について前進消去をする ) 回目 x 2 x 1 x 3 ( ピボッティング : 列方向 :2 列と 3 列の入れ替えをする ) 回目 x 2 x 3 x 1 (3 行について前進消去をする ) 回目 x 2 x 3 x 1 ( 前進消去終わり )

5 実行結果 連立方程式を解く際には解の種類が一意に決まるか否かあらかじめ調査するとよい. 前述の例題 1 と 2 の係数の rank はそれぞれ 3 であり,x の数 n と同じであった. つまり 1 種類の解が存在することになる. 連立方程式の係数によっては rank<n となることもある. そのような場合では無数の解が存在する. その様な問題についてこのプログラムを実行すると, そのうちの一つだけの解が求まることになる. 例えば次の問題を考える. x 1 + 2x 2 + 3x 3 = 14 { 4x 1 + 5x 2 + 6x 3 = 32 7x 1 + 8x 2 + 9x 3 = 23 係数部分をまとめると次のようになる. ( 4 5 6) この rank は 2 であり, 解は無数にある. 例えば次のようなものがある. x 1=1,x 2=2,x 3=3 x 1=2,x 2=0,x 3=4 x 1=10,x 2=-16,x 3=12 従って,n 元連立 1 次方程式の問題を解く際には rank を求めておくとよい. 次に例題 2 についても計算してみよ.( プログラムの編集については付録を参照のこと.) このプログラムが完全ピボッティングに対応していることを確認せよ. レポート C n 元連立 1 次方程式で表現できる問題を見つけ, プログラムを使って解を求めよ. 問題は自 身で行っている研究の内容や工学的内容が望ましい. レポートの書式はレポート A に従うこと. 参考文献 1) 佐藤次男ほか,C による理工学問題の解法, 日刊工業新聞社,ISBN

6 付録ガウスの消去法のプログラム 黄色い部分にセットされているデータは例題 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; double epsilon=1e-18; double x[nunber_of_cf]; // 未知変数 (x1,x2,x3) の個数 // ピボット 対角要素の下限値をセットする // 連立方程式の解 ( 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 // // x[0] =x1, x[1] =x2,..., x[n-1] =xn

7 // 下の 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) の計算

8 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 は, 列の入れ替えのピボッティングのとき, 列の入れ替えの位置を記憶するために使われる. 列の入替では解の位置が変わってしまう. 具体的例を以下に示す. 例 : x 1 x 2 x 3 x n (x 2 と x 3 の列を入れ替えると ) x 1 x 3 x 2 x n この状態では後退消去で連立方程式の解 (x 1, x 2, ) を求めることはできない. つまり列の入れ替えのピボッティングが終わった後に, また列の要素の位置を交換し, その後, 後退消去を計算することになる. そのためピボッティング以前の解の位置を変数 LN に保存する.