第7回 連立方程式の解法(その2)

• 今週はオンデマンドでの実施とします。
• Oh-o! Meijiのクラスウェブに、解説動画へのリンクを掲載してあります。この動画に沿って作業を進めてください。
• プログラムを作成する箇所では、一旦動画を止め、プログラムが完成して動作を確認してから、また動画を再生してください。

7.1 ピボット選択

アルゴリズム

1. a₁₁〜a_n1の中で絶対値が最大となる係数がある行を探す。ここでは s行目の a_s1 の絶対値が最大であるものとする。
• すなわち、先ず何らかの方法で絶対値が最大となる行を探す必要がある。

2. 1行目とs行目を交換する。
   
3. a_s1をピボットにして掃き出す。
   
4. 次の行に移る。以降繰り返し。

7.2 例題

7.3 ガウス消去法

7.3.1 前進消去

7.3.2 後退代入

7.4 ガウス消去法のコーディング

7.4.1 前進消去のアルゴリズム

1. ピボットを1行1列からn-1行n-1列に移しながら以下の操作を繰り返す。このとき、7.1 のピボット選択を行いながら操作を進める。
2. a_kk をピボットにし、i 行について a_ik/a_kk を求め、
   i行 - k行×a_ik / a_kk を行う。

7.4.2 後退代入のアルゴリズム

1. b'_n から b'₁ に向かって以下の操作を繰り返す。
2. b'_i を初期値として、b'_i+1 から b'_n まで a'_ij×b'_j の値を引く。
3. 2.で求めた値を a'_ii で割る。

課題 (提出〆切：金曜22時)

1. 行列 A とベクトル bのデータが A4_mat.csv, b4_vec.csv に記述されている。 ガウス・ジョルダン法およびガウス消去法を用いて Ax = b を x について解き、方程式の解を求めるプログラムを作成せよ。 以下の指示に従うこと。
   • プログラムのファイル名は7_1.cppとすること。
   • 解の個数は既知であるものとしてプログラムを作成して構わない。
   • ファイルはコマンドライン引数から読み込むようにすること。
   • ガウス・ジョルダン法およびガウス消去法の両方にピボット選択のアルゴリズムを適用すること。
   • 計算過程が把握できるように、1行の処理が終わるごとに係数行列を画面に出力すること。 また、行を入れ替えた際は、何行目と何行目を入れ替えたのかを画面に表示させること。
   • ガウス消去法のプログラムは、以下のような関数構成で作成すること。

     void forward_elimination(Matrix A){
       // 前進消去のコードだけを書くこと
     }
     
     void backward_substitution(Matrix A){
       // 後退代入のコードだけを書くこと
     }
     
     // ガウス消去法を使って連立方程式を解く
     void gauss_elimination(Matrix A){
       forward_elimination(A); // 前進消去
       backward_substitution(A); // 後退代入
     }
     

   • 以下は実行例である。実際の計算結果とは異なるため注意すること。

     Z:\compmech\ex7>7_1.exe A4_mat.csv b4_vec.csv
        2.000   1.000   3.000   3.000  |  35.000
        3.000   1.000   0.000   2.000  |   9.000
        1.000  -2.000   3.000   5.000  |  32.000
        1.000   1.000  -3.000   1.000  | -23.000
     
     ----- Solution by the Gauss Jordan method -----
     Changed row 1 and row 2.
       3.000   1.000   0.000   2.000  |   9.000
       2.000   1.000   3.000   3.000  |  35.000
       1.000  -2.000   3.000   5.000  |  32.000
       1.000   1.000  -3.000   1.000  | -23.000
     ...
     
     ----- Answers -----
     x[1] =  1.000
     x[2] =  2.000
     x[3] =  3.000
     x[4] =  4.000
     
     ----- Solution by the Gauss Elimination method -----
     < Forward elimination >
     Changed row 1 and row 2.
       3.000   1.000   0.000   2.000  |   9.000
       2.000   1.000   3.000   3.000  |  35.000
       1.000  -2.000   3.000   5.000  |  32.000
       1.000   1.000  -3.000   1.000  | -23.000
     ...
     
     ----- Answers -----
     x[1] =  1.000
     x[2] =  2.000
     x[3] =  3.000
     x[4] =  4.000
     
     

2. ガウス・ジョルダン法とガウス消去法のプログラムの実行時間を計測するプログラムを作成し、実行時間を比較せよ。 以下の指示に従うこと。
   • プログラムのファイル名は7_2.cppとすること。
   • 解の個数は既知としてプログラムを作成して構わない。
   • 行列 A とベクトル bのデータは A100_mat.csv, b100_vec.csv を用いること。これらのファイルはコマンドライン引数から読み込むようにすること。
   • 同じ連立方程式を1000回繰り返して解くことにより、実行時間を評価すること。
     これは、連立方程式を1回解くだけでは実行時間が短すぎて、正しい比較が出来ないためである。
   • 関数timeGetTimeは、Windowsが起動してから経過した時間をミリ秒単位で返す関数である。
     例えば、以下のようなプログラムにより、足し算を1億回繰り返したときの実行時間を計測することが出来る。
     このサンプルプログラムを参考にして、ガウス・ジョルダン法とガウス消去法のプログラムの実行時間を計測すること。
     なお、MacOS などWindows 環境以外でプログラムを作成する人は、clock 関数など使用するとよい。

     // time_test.cpp
     #include <stdio.h>
     #include <windows.h>
     // ↓ clock 関数を用いて時間計測する人向け
     // #include <time.h>
     
     int main(void) {
       DWORD time1, time2;
       // ↓ clock 関数を用いて時間計測する人向け
       // clock_t time1, time2;
     
       float a = 0;
             
       time1 = timeGetTime();
       // ↓ clock 関数を用いて時間計測する人向け
       // time1 = clock();
     
       // 足し算を1億回繰り返す
       for (int i=0;i<100000000;++i) {
         a += 0.001;
       }
       time2 = timeGetTime();
       // ↓ clock 関数を用いて時間計測する人向け
       // time2 = clock();
             
       printf("Calculation time : %ld [ms]\n", time2-time1);
       // ↓ clock 関数を用いて時間計測する人向け
       // printf("Calculation time : %f [s]\n", (double)(time2-time1)/CLOCKS_PER_SEC);
             
       return 0;
     }

   • プログラムは7_1.cpp をもとに作成するとよい。ただし、画面には、計算のプロセスおよび計算結果を一切表示させず、1000回繰り返し計算を行ったときの最終的な実行時間のみを表示させること。 画面に文字を出力する処理 (printf, print_mat())は、連立方程式を解く処理と比較して非常に重い。このため、計算の経過を画面に出力すると正しい比較が出来ない。
   • 以下は実行例である。OS はwindows10、CPUは12th Gen Intel(R)Core(TM)i7-12700、メモリは16GB の環境において計測した結果である。 各自のプログラムやコンパイルする環境に応じて実行時間の計測結果が異なることが予想されるため、参考程度に確認するとよい。
     上記の計測環境においては、関数timeGetTime を用いた場合でも、Clock 関数を用いた場合でも、ほぼ同じ計測結果を得ることができた。

     Z:\compmech\ex7>7_2.exe A100_mat.csv b100_vec.csv
     Number of Repetitions : 1000
     
     ----- Gauss Jordan method -----
     Calculation time : 636 [ms]
     
     ----- Gauss Elimination method -----
     Calculation time : 407 [ms]
     
     

   • 解の個数を n とすると、ガウス・ジョルダン法の計算量 (乗算の数)は約n^3/2、ガウス消去法の計算量は約n^3/3 であることが知られている。
     多自由度のシステム (自動車など) を1秒間に1000 回計算するようなシミュレーションを実施するにあたって、計算速度は重要な要因となることから、連立方程式の解を求めるためにガウス・ジョルダン法よりもガウス消去法が用いられている。

課題サンプルプログラム

• 7_1.cpp
• 7_2.cpp