Microsoft Word - COMP-MATH-2016-FULLTEXT.doc

#13 2017/1/17 16 連立一次方程式を解くその 2 1) 対角線上に 0 が現れる場合を回避する 0x + 3y = 9 2x + 4y = 16 これを kadai-26 で解いてみなさい 以下のようにエラーになる 0.00000000 3.00000000 9.00000000 2.00000000 4.00000000 16.00000000-1.#IND0000-1.#IND0000-1.#IND0000-1.#IND0000-1.#IND0000-1.#IND0000 でも 式 (0) と (1) をひっくり返したら エラーにならない 2x + 4y = 16 0x + 3y = 9 2.00000000 4.00000000 16.00000000 0.00000000 3.00000000 9.00000000 1.00000000 0.00000000 2.00000000 0.00000000 1.00000000 3.00000000 だから 対角線上に 0 が現れたら 式の順序を変えて 対角線上に 0 が来ないようにすればよい 上の例 では 式 (0) と (1) を交換する (Swap) ことによって 対角線上の 0 を回避している ( 式の順番を変えても 解 には影響がないから ) [i] 番目の式と [k] 番目の式を交換 (Swap) するとは a[i][0] a[i][1] a[i][j] a[k][0] a[k][1] a[k][j] ところで 二つの数を交換 (Swap) することをプログラムではどう書くか? a = 1.0; b=2.0; // a と b を交換するぞ a = b; b = a; これではダメ元の a の内容が残っていない そこで c = a; a = b; b = c; 待避場所として 一個の変数 ( ここでは C) を用いて 3 ステップの手順が必要になる 1/12/2017 P-111

これを係数の配列で行う (Ex24.c Kadai-27.c をコピーして修正する ) #define N 2 A[N][N+1] = 0,3,9,2,4,16 ; pra(),swapa(); // 交換する関数の宣言を追加する P,Q; i,j,k; for(i = 0; i < N; ++i) if(a[i][i] == 0.0) // 対角線項が 0 なら次の行と交換する P = A[i][i]; swapa(i,i+1); A[i][j] /= P; // 確認のために印字 // 改めて 対角線上の係数を P に置く // 以下 Kadai-27.c と同じ for(k = 0; k < N; ++k) if(k!= i) Q = A[k][i]; A[k][j] -= Q*A[i][j]; swapa(i,k) // i 行と k 行を交換する関数を追加する i,k; j; c; // 待避場所として c を使う for(j = 0; j < ; ++j) c = ; // c に何かを退避する A[i][j] = A[k][j]; = c; // c を何かに戻す pra() 以下 同じなので省略 上の ex24.c を完成させ 以下を解けることを確認しなさい 1/12/2017 P-112

0x + 3y = 9 2x + 4y = 16 ところで a[i][i] が 0 の場合 その下の係数 つまり a[i+1][i] も 0 かもしれない 1 0 a[0][i] a[0][n-1] a[0][n]= 右辺 [0] 0 1 0 0 a[i][i] a[i+1][i] 0 0 a[n-1][i] そこで a[i][i] から下に見ていって 絶対値が一番大きな係数を持つ行を見つけて その行と入れ替 えるという方針にする もし 上の網掛け部分が全部 0 なら それ以上掃き出しはできないのだから 解が求められない として終了して良い ( 実際には 解がない 場合と 解は決まらない 場合があるが これを区別するのは この 課題の範囲を超える ) まず findmax(i) という関数を用意する これは A[i][i] から下に係数を見て 最大の絶対値を持 つ係数 a[k][i] の k を求める もし 全部 EPSIRON(1.0e-10) 以下なら 全部ゼロと見なして処理を 中止する ------------ 練習 : 最大のもの を見つける 配列 A にデータの数値が入っている その中の最大の数値と 配列の番号 ( 添字 ) を求めなさい ex25.c #define N 10 A[10] = 23,45,21,44,68,96,13,77,62,35; i,imax,valmax; valmax = A[0]; // 仮の最大値 imax = 0; for( i = 0; i < N; ++i) if(a[i] > valmax ) valmax = A[i]; imax = i; prf(" 最大は A[%d] = %d\n",imax,valmax); 結果は最大は A[5] = 96 要するに 仮の最大値を適当に決めて それより大きいものが現れたら 置き換えていく これを踏まえて 1/12/2017 P-113

ex24. に以下の太字部分を 変更または加える #define N 3 A[N][N+1] = 0,1,0,1,0,0,1,2,1,0,0,3 ; pra(),swapa(); i,j,k; P,Q; for(i = 0; i < N; ++i) k = findmax(i); // i 式以降で絶対値が最大になる式の番号 k を求める if( k!= i) swapa(i,k); // それが [i] でなければ [i] と [k] 行を交換する // 確認の印字 P = A[i][i]; // 以下同様 A[i][j] /= P; for(k = 0; k < N; ++k) if(k!= i) Q = A[k][i]; A[k][j] -= Q*A[i][j]; #define EPSIRON (1.0e-10) // これ以下の値はゼロと見なす #include <math.h> // fabs() 絶対値を使うために必要 findmax(i) i; // A[i][i] から下を調べる k,kmax; amax = EPSIRON; for(k = i; k < N; ++k) if(fabs(a[k][i]) > amax) if(amax <= EPSIRON) return(kmax); 以下 ex29.c と同じなので省略 swapa(i,k) の本体 pra() の本体 // 今までの最大よりも kmax = k; // 大きい場合 amax = fabs(a[k][i]); // それを記録する prf("******** n 解が求められません n"); exit(); // ここでプログラムを終了させる 1/12/2017 P-114

ex24 を実行すると こうなる 0.00000000 1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 1.00000000 2.00000000 1.00000000 0.00000000 0.00000000 3.00000000 1.00000000 0.00000000 0.00000000 3.00000000 //0 行と2 行が交換された 0.00000000 0.00000000 1.00000000 2.00000000 0.00000000 1.00000000 0.00000000 1.00000000 1.00000000 0.00000000 0.00000000 3.00000000 0.00000000 1.00000000 0.00000000 1.00000000 // 1 行と2 行が交換された 0.00000000 0.00000000 1.00000000 2.00000000 1.00000000 0.00000000 0.00000000 3.00000000 0.00000000 1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 1.00000000 2.00000000 // 最後なので交換はしない 意図的に 解が決まらないようにしてみると #define N 3 A[N][N+1] = 0,1,0,1,0,0,1,2,0,0,0,3 ; 0.00000000 1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 1.00000000 2.00000000 0.00000000 0.00000000 0.00000000 3.00000000 ******** 解が求められません // 0 列目は全部 0 なので 処理できない give up! ここまでで 一応掃き出し法の完成とする 仕上げに 係数をデータとして読むようにしよう 1/12/2017 P-115

以下を最終形 Kadai28.c とする 太字以外は ex24.c と同じ #define NMAX 30 // 連立の最大は 30 までとする A[NMAX][NMAX+1]; // データは不要 N; // N は入力させる pra(),swapa(); P,Q; i,j,k; prf(" 式の個数 N は?"); scanf("%d", ); // 整数 N を入力 for(i=0; i < ; ++i) prf(" 第 %d 式の係数をスペースで区切って入力しなさい : ",i); for(j = 0; j < ; ++j) scanf("%lf", & ); //(%lf はエル エフ long float) // 以下 ex29.c と同じ for(i = 0; i < N; ++i) k = findmax(i); if( k!= i) swapa(i,k); P = A[i][i]; A[i][j] /= P; for(k = 0; k < N; ++k) if(k!= i) Q = A[k][i]; A[k][j] -= Q*A[i][j]; #define EPSIRON (1.0e-10) #include <math.h> findmax(i) i; k,kmax; amax = EPSIRON; for(k = i; k < N; ++k) if(fabs(a[k][i]) > amax) kmax = k; 1/12/2017 P-116

amax = fabs(a[k][i]); if(amax <= EPSIRON) prf("******** n 解が求められません n"); exit(); return(kmax); swapa(i,k) i,k; j; c; for(j = 0; j < N+1; ++j) c = A[i][j]; A[i][j] = A[k][j]; A[k][j] = c; pra() i,j; prf(" n"); for( i = 0; i<n;++i) for( j=0; j < N+1; ++j) prf("%12.8f ",A[i][j]); prf(" n"); 実行例 式の個数 N は?3 第 0 式の係数をスペースで区切って入力しなさい : 27 9 3 14 第 1 式の係数をスペースで区切って入力しなさい : 8 4 2 5 第 2 式の係数をスペースで区切って入力しなさい : 1 1 1 1 27.00000000 9.00000000 3.00000000 14.00000000 8.00000000 4.00000000 2.00000000 5.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 0.00000000 0.00000000 0.33333333 0.00000000 1.00000000 0.00000000 0.50000000 0.00000000 0.00000000 1.00000000 0.16666667 h 秀丸でデータをファイル ( 例えば test.dat) を作っておく 3 27 9 3 14 8 4 2 5 1 1 1 1 実行は Kadai28 < test.dat とすれば 毎回データを入力しなくてもよい Kadai28.c を完成させ 下を解きなさい 1/12/2017 P-117

a + b + c + d = 10 2a + b + 3c + d = 17 3a + 5b + c + 2d = 24 a b + c d = 2 挑戦 : 検算もしてくれるとうれしい 3. 検算の仕方入力された係数行列 A のコピーを B として保持する A の N 列に解が求まったあと B の係数と解から方程式の右辺を計算する Kadai29.c(Kadai28.c の修正 ) #define NMAX 30 N; A[NMAX][NMAX+1]; B[NMAX][NMAX+1]; // 係数をコピーしておくための配列 pra(),swapa(),copyab(),kenzan(); // 関数を 2 こ追加 i,j,k; P,Q; prf(" 式の個数 N は?"); scanf("%d",&n); for(i=0; i<n; ++i) prf(" 第 %d 式の係数をスペースで区切って入力しなさい : ",i+1); for(j = 0; j < N+1; ++j) scanf("%lf",&a[i][j]); copyab(); // この時点のAをBにコピーする for(i = 0; i < N; ++i) k = findmax(i); if( k!= i) swapa(i,k); P = A[i][i]; A[i][j] /= P; for(k = 0; k < N; ++k) if(k!= i) Q = A[k][i]; A[k][j] -= Q*A[i][j]; 1/12/2017 P-118

kenzan(); // 検算の実行 copyab() // AをBにコピーする関数 i,j; for(i =0; i < N; ++i) for(j=0; j < N+1; ++j) B[i][j] = A[i][j]; kenzan() // 検算する関数 i,j; S; prf(" 検算 --------\n"); for(i = 0; i < N; ++i) for(j = 0, S = 0; j < N; ++j) prf("%10.6f * %10.6f ",B[ ][ ],A[ ][ ]); if(j!= N-1) prf("+ "); S += B[ ][ ]*A[ ][ ]; prf("= %f\n", ); このあとは Kadai28.c とおなじなので 省略 実行例 式の個数 N は?3 第 1 式の係数をスペースで区切って入力しなさい : 27 9 3 14 第 2 式の係数をスペースで区切って入力しなさい : 8 4 2 5 第 3 式の係数をスペースで区切って入力しなさい : 1 1 1 1 27.00000000 9.00000000 3.00000000 14.00000000 8.00000000 4.00000000 2.00000000 5.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 0.00000000 0.00000000 0.33333333 0.00000000 1.00000000 0.00000000 0.50000000 0.00000000 0.00000000 1.00000000 0.16666667 検算 -------- 27.000000*0.333333 + 9.000000*0.500000 + 3.000000*0.166667 = 14.000000 8.000000*0.333333 + 4.000000*0.500000 + 2.000000*0.166667 = 5.000000 1.000000*0.333333 + 1.000000*0.500000 + 1.000000*0.166667 = 1.000000 1/12/2017 P-119