はじめMath! Javaでコンピュータ数学

第44回行列の数学 逆行列と連立一次方程式 [後編]

前回学習したガウス-ジョルダン法のアルゴリズムをJava言語で実装しましょう。少々長いサンプルコードですが、⁠コード・リーディング」の訓練だと思ってじっくり読み込んでから演習に取り組んでみてください。

問題 ガウス-ジョルダン法で連立一次方程式を解く次のプログラムを完成させましょう。

ガウス-ジョルダン法によるプログラムを次に示します。肝心のピボット処理、掃き出し処理のメソッドが未完成です。不足を補って完成させましょう。ソースコードにはスモークテスト[1]が含まれています。未完成版をそのままコンパイルしても、一応動作し次のような出力を得ることが出来ます。

c:\temp>java Sample_GaussJordanElimination
Gauss-Jordan Elimination Test Start
  関数が望んだtrue を返しませんでしたabxROH
  関数が望んだtrue を返しませんでしたabxNROH
  関数が望んだfalse を返しませんでしたa2b2xROH
  関数は真を返しましたが、正しい解ではありませんでした。
a2b2xROH
Test Done

正しくコードを完成させると、次のような出力を得ることが出来るでしょう。

C:\temp>java Sample_GaussJordanElimination
Gauss-Jordan Elimination Test Start
  関数は真を返しましたが、正しい解ではありませんでした。
abxNROH
  関数が望んだtrue を返しませんでしたa2b2xROH
Test Done

先ずは未完成版のソースをじっくり読んで、大まかなコードの流れをつかみましょう。それが出来たら、ピボット処理、掃き出し処理のメソッドを順に完成させてください。コツは「常に動作するようにコードを加えていくこと」です。コードを加えてはコンパイルし、実行するのです。コードが動かないというのは、非常にモチベーションを奪います。しかし、コードが動いているうちは何となく気持ちが前向きです。投げ出さずに続けるための、自分のためのサプリメントだと思ってください。では、健闘を祈ります!

ソースコード:Sample Gauss Jordan Elimination.java未完成版
//サンプルコード
//Gauss-Jordan Elimination(ガウス-ジョルダン法)で
//方程式の解を求める。配列版。
//filename : Sample_GaussJordanElimination.java

class Sample_GaussJordanElimination {

  public static void main(String[] args) {

    //---------------------------------
    //テスト用のデータを準備
    //---------------------------------
    //一意な解が得られる場合
    double a[][] = //係数行列A
      { { 0, 1, 2 },
        { 1, 2, 1 },
        { 1,-1, 1 } };
    double b[][] = //定数ベクトルB
      { { 3 },
        { 6 },
        { 3 } };
    double x[][] = //解ベクトルX
      { { 0 },
        { 0 },
        { 0 } };
    double ROH[][] = //望まれる解ベクトル(Result Of Hope)
      { { 3 },
        { 1 },
        { 1 } };
    double NROH[][]= //間違った解ベクトル(Not Result Of Hope)
      { { 3 },
        { 1 },
        { 2 } };
    //一意な解が得られない場合
    double a2[][] = //係数行列A2
      { { 1, 0, 3 },
        { 2, 3, 4 },
        { 1, 3, 1 } };
    double b2[][] = //定数ベクトルB2
      { { 1 },
        { 3 },
        { 2 } };
    //上記の行列サンプルデータは、
    // 『Java2 による数値計算』
    // 赤間世紀 著
    // 技報堂出版 刊 P.91
    // 『新版代数と幾何』
    // 古屋茂 監修 
    // 飯田・上田・田河・田中 共著
    // 大日本図書 P.115
    //を利用させていただきました。 


    //---------------------------------
    //ガウスジョルダン法を実行します
    //---------------------------------
    System.out.println("Gauss-Jordan Elimination Test Start");

    //Test 1
    //計算結果が望んだ解と等しい。isDoRight_ROH もtrue を指定。
    //正しく終了し、特別なメッセージを出力しないはず。
    test_isGJE_DoRight(a,"a",b,"b",x,"x",ROH,"ROH",true);

    //Test 2
    //計算結果が望んだ解と異なる。isDoRight_ROH にtrue を指定。
    //関数はfalse を返すので、矛盾する条件指定である。
    //異常ととらえ、メッセージを出力するはず。
    test_isGJE_DoRight(a,"a",b,"b",x,"x",NROH,"NROH",true);
    //Test 3
    //一意な解が得られない。isDoRight_ROH がfalse。
    //予想通りの結果であるため、特別なメッセージを出力しないはず。
    test_isGJE_DoRight(a2,"a2",b2,"b2",x,"x",ROH,"ROH",false);

    //Test 4
    //一意な解が得られない。isDoRight_ROH がtrue。
    //関数はfalse を返し、矛盾する条件指定のため、
    //特別なメッセージを出力するはず。
    test_isGJE_DoRight(a2,"a2",b2,"b2",x,"x",ROH,"ROH",true);

    System.out.println(" Test Done");

  }// end of main


/*
 * test_isGJE_DoRight
 * 目的: isGJE_DoRight メソッドのテスト
 * 引数: a 係数行列
 *      b 定数ベクトル
 *      x 解ベクトル
 *      ROH 解ベクトルの正解
 *      isDoRight_ROH 正しく解が得られたか。
 *      valName* 引数に渡した変数名
 * 戻り値: boolean
 *        evokeGaussJordanElimination の結果と
 *        isDoRight の結果が等しいときに真を返す。
 *        evokeGaussJordanElimination が真、
 *        isDoRight の値が真であっても、
 *        望まれる解が得られなかった場合は偽
 */
  static boolean test_isGJE_DoRight(double a[][],String valNameA,
                                    double b[][],String valNameB,
                                    double x[][],String valNameX,
                                    double ROH[][],String valNameROH,
                                    boolean isDoRight_ROH){
    boolean GJEresult = evokeGaussJordanElimination(a,b,x);
    if (GJEresult == isDoRight_ROH) {
      if (GJEresult == false){
        //System.out.println(" 望んだ結果であるが、"+
        // "一意な解は得られませんでした。"
        // + valNameA + valNameB + valNameX);
        return true;
      }
      if (isMatrixEqualsDouble(x,ROH) == true) {
        //System.out.println(" 望んだとおりの結果であり、" +
        // "関数の返した答えは正しい値でした。"
        // + valNameA + valNameB + valNameX);
        return true;
      } else {
        System.out.println(" 関数は真を返しましたが、" +
             "正しい解ではありませんでした。"
             + valNameA + valNameB + valNameX + valNameROH);
        return false;
      }
    } else {
      System.out.println(" 関数が望んだ"+isDoRight_ROH+
             "を返しませんでした"
             + valNameA + valNameB + valNameX + valNameROH);
      return false;
    }
  }// end of test_isGJE_DoRight


/*
 * evokeGaussJordanElimination
 * 目的:引数で与えた行列について
 *      Gauss-Jordan Elimination を実施する。
 * 引数: a[][] 係数行列。
 *    : b[][] 定数ベクトル。
 *    : x[][] 解ベクトル。
 * 戻り値: boolean 一意な解が得られれば真
 *         一意な解が得られなければ偽を返す。
 *         引数に与えた変数に問題があった場合も偽を返す。
 */
  static boolean evokeGaussJordanElimination(
         double a[][], double b[][], double x[][]){
    //解ベクトルをゼロクリア
    for(int row = 0; row < a.length; ++row) x[row][0] = 0;
    //係数行列、定数・解ベクトルの次元数は適切か
    if ( isSquareMatrixDouble(a) &&
         (a.length == b.length) &&
         (a.length == x.length) ){
      //引数に与えた変数に問題なし。
      for (int k=0; k < a.length; ++k){
        //k は現在注目しているピボット位置
        regulatePivot(a,b,k); //ピボットの調整
        if (a[k][k] == 0){
          //ピボットがゼロのとき、正則ではない。異常終了
          return false;
        } else {
          //ピボットの値が正常なので掃き出し実行。
          doSweep(a,b,k);
        }
      }// of for k
      //ここまで来たと言うことは、正しく計算が終了しているということ。
      for(int row = 0; row < a.length; ++row) x[row][0] = b[row][0];
    } else {
      //引数に与えた変数に問題がある。
      return false;
    }
    //ここまで到達したということは成功
    return true;
  }// end of evokeGaussJordanElimination


/*
 * isMatrixEqualsDouble
 * 目的:引数で与えた二つの行列が等しいか
 * 引数: a[][]
 *    : b[][]
 * 戻り値: boolean 等しい行列であるなら真を返す。
 */
  static boolean isMatrixEqualsDouble(double x[][], double ROH[][]){
    if ( (x.length == ROH.length) ){
      for (int row = 0; row < x.length; ++row){
        if (x[row].length != ROH[row].length) return false;
        for (int col = 0; col < x[row].length; ++col){
          if (x[row][col] != ROH[row][col]){
            //Value not equal
            return false;
          }
        }
      }
      return true;
    } else {
      //Matrix is not equal
      return false;
    }
  }// end of isMatrixEqualsDouble


/*
 * isSquareMatrixDouble
 * 目的:引数に与えた行列を表す配列が、正方行列かどうか判定する。
 *      行数と一行に入っている要素数が等しいか確認する。
 * 引数: m[][] 判定する行列
 * 戻り値: boolean 正方行列なら真
 */
  static boolean isSquareMatrixDouble(double m[][]){
    int row = m.length; //行数取得
    int cul = m[0].length; //一行目の要素数(列数)の取得
    boolean result = true;
    if (row == cul) {
      for (int i = 0 ; i < row ; ++i ){
        if (cul == m[i].length){
          //列数は妥当である。
        } else {
          //列数が妥当ではない。
          result = false;
          break;
        }
      }
    } else {
      result = false;
    }
    return result;
  }// end of isSquareMatrixDouble


/*
 * regulatePivot
 * 目的: 係数行列のピボット値が適正になるように調整する。
 *      具体的には、現在注目しているピボットのある列について、
 *      最も大きな数値のある行と、現在のピボット行を交換する。
 * 引数: a[][] 係数行列。
 *    : b[][] 定数配列。
 *    : k 現在注目しているピボット位置
 */
  static void regulatePivot(double a[][],double b[][] ,int k){

  // ここから 目的のコードを加えましょう。
  //
  // ここまで

  }// end of regulatePivot


/*
 * doSweep
 * 目的:掃き出し処理。
 * 引数: a[][] 係数行列。
 *    : b[][] 定数配列。
 *    : k ピボット位置
 */
  static void doSweep(double a[][],double b[][],int k){

  // ここから 目的のコードを加えましょう。
  //
  // ここまで

  }// end of doSweep

}// end of class Sample_GaussJordanElimination

解説

コードの不足していた二つのメソッドのみを示します。完成後はサンプルコードで示したテストデータ以外の様々な値で計算させてみると良いでしょう。

ソースコード: Gauss Jordan Elimination.java完成版
/*
 * regulatePivot
 * 目的: 係数行列のピボット値が適正になるように調整する。
 *      具体的には、現在注目しているピボットのある列について、
 *      最も大きな数値のある行と、現在のピボット行を交換する。
 * 引数: a[][] 係数行列。
 *    : b[][] 定数配列。
 *    : k 現在注目しているピボット位置
 */
  static void regulatePivot(double a[][],double b[][] ,int k){
    double max = Math.abs( a[k][k] );
    double temp = 0;
    int rowOfMaxPivot = k;//最大のピボット候補値の格納された行の値
    if(k != (a.length - 1)){//最終行でなければ続きを実行
      for (int row = k+1; row < a.length; ++row){//ピボット行の次の行から
        if( Math.abs(a[row][k]) > max ){//より大きな値を見つけて
          max = Math.abs(a[row][k]);//最大値として保管し
          rowOfMaxPivot = row;//その行が何行目か保管しておく。
        }
      }
    }
    if(rowOfMaxPivot != k){//最大値のある行が、現在注目しているピボットのある行でなければ
      //最大値のある行と、現在のピボット行とで値の交換
      //現在のピボットのある列を含めて右側の値について交換を実施。
      for (int col = k; col < a.length; ++col){
        temp = a[k][col];
        a[k][col] = a[rowOfMaxPivot][col];
        a[rowOfMaxPivot][col] = temp;
      }
      temp = b[k][0];
      b[k][0] = b[rowOfMaxPivot][0];
      b[rowOfMaxPivot][0] = temp;
    }
  }// end of regulatePivot


/*
 * doSweep
 * 目的:掃き出し処理。
 * 引数: a[][] 係数行列。
 *    : b[][] 定数配列。
 *    : k ピボット位置
 */
  static void doSweep(double a[][],double b[][],int k){
    double pivotVal = a[k][k]; //掃き出し前のピボット値を保持
    double pivotColVal = 0; //これから掃き出しを行う行のピ
                            //ボット列の値を保持しておく。
    //ピボット値のある行の各要素を、ピボット値で割る。
    for( int col=k; col < a.length; ++col)
      a[k][col] = a[k][col] / pivotVal;
    b[k][0] = b[k][0] / pivotVal;
    //全ての行について、ピボット位置以外のピボット列の要素をゼロにする。
    //その副作用として、ピボット列以外の要素同士も減算を行う。
    for( int row = 0; row < a.length; ++row){
      pivotColVal = a[row][k];
      if( row != k ){
        for( int col = k; col < a.length; ++col)
          a[row][col] = a[row][col] - pivotColVal * a[k][col];
        b[row][0] = b[row][0] - pivotColVal * b[k][0];
      }
    }
  }// end of doSweep

ガウス-ジョルダン法の注意点

ガウス-ジョルダン法に限らないことですが、浮動小数点数を用いた計算の結果は近似値です。必ず誤差を含んでいます。今回の問題のために用意した係数行列・定数行列は、解として整数値が得られるような行儀の良いものを準備しました。ですからテストがおとなしく終了したのです。ところが、そのほかの値を用いたとき、期待する結果が1なのに、得られた計算結果が0.9999になるかもしれません。人間の目から見れば、ほぼ同じ値、まあ、計算誤差だろう、と推察できるのですが、コンピュータにとっては異なる値にしか見えません。テストは失敗だ、と判定してしまうことでしょう。妥当なテストを行うためには、コンピュータにどの程度の誤差を許容させるかという条件付けを指示してやる必要があります。実数型を用いた計算では、常にこの誤差の存在を忘れないようにしましょう。プログラムが正しく動作しているのに、バグを疑う無駄をしてしまいかねません。

今回はここまで

今回の課題はいかがだったでしょうか。ほんのわずかなコード追加なのですが、正しく動作させるまでにはずいぶん努力を要したことでしょう。解説に示したコードはまだまだ改善の余地が多々残されています。オーバーフロー、アンダーフローといった計算に関する問題を回避すること、メモリの節約に関すること、考えれば数多くあることでしょう。皆さんがガウス-ジョルダン法を実際に利用したい場合には、この解説で示したコードで満足することなく、さらなる正確性、パフォーマンスアップを追い求めてください。

おすすめ記事

記事・ニュース一覧