はじめに

今回は、ガウス・ジョルダン消去法プログラムを掲載します。
こちらは、以前掲載したガウスの消去法アルゴリズムを実装したプログラムの一部を変更したものとなっております。

ガウスジョルダン消去法とは

簡単にいうと、行列計算AX = Bにおいて、行列Aを単位行列に変換し、解行列Xを求める手法です。

詳しくは、下記記事などにまとめられています。ここでは、割愛します。。

ガウス・ジョルダン法
 山本研究室/ ガウス・ジョルダン法

プログラム

以前のガウスの消去法プログラムとの主な違いは、calculateメソッドのみです。
calculateメソッドにガウスジョルダン消去法アルゴリズムが記述しています。

メインプログラム

計算を行うメインプログラムファイルを準備します。
このプログラムに必要モジュールをインポートします。

calc.rb

require './gauss_jordan_elimination'

if __FILE__ == $0
  begin
    # 計算クラスインスタンス化
    obj = GaussJordanElimination.new
    # 連立方程式を解く
    obj.display_matrixA_and_matrixB
    obj.calculate
    obj.display_matrixX
  rescue => e
    $stderr.puts "[#{e.class}]"
    $stderr.puts "#{e.message}"
    e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
  else
    puts "finish (no error)"
  ensure
    puts "all finish"
  end
end

モジュール

ガウス・ジョルダン消去法クラスを設定し、その中に必要メソッドを定義します。

gauss_jordan_elimination.rb

require 'csv'

class GaussJordanElimination
=begin
  A. インスタンス変数
    @a: 行列A
    @b: 行列B
    @n_a: 行列Aの次元
    @n_b_col: 行列Bの列数

  B. メソッド
    initialize: コンストラクタ
    display_matrixA_and_matrixB: 行列A, Bを出力
    display_matrixX: 行列X(解)を出力
    calculate: ガウス・ジョルダンの消去法で計算を行う
    partial_pivoting: 部分ピボット選択
=end

  def initialize
    # 初期値
    @a = [
      [-1, 2, 2],
      [2, 3, -1],
      [1, 1, -1]
    ]
    @b = [
      [1, 2],
      [-3, 2],
      [-2, 2]
    ]

    if ARGV[0] && ARGV[1]
      @a = CSV.readlines("./#{ARGV[0]}")
      @b = CSV.readlines("./#{ARGV[1]}")
    end

    @n_a = @a.size
    @n_b_col = @b[0].size
  end

  def display_matrixA_and_matrixB
    # 行列Aを出力
    puts "A ="
    for i in 0...@n_a
      for j in 0...@n_a
        printf("%+.5f\t", @a[i][j])
      end
      puts ""
    end

    # 行列Bを出力
    puts "B ="
    for i in 0...@n_a
      for j in 0...@n_b_col
        printf("%+.5f\t", @b[i][j])
      end
      puts ""
    end
  end

  def display_matrixX
    # 行列X(解)を出力
    puts "X (satisfying AX=B) ="
    for i in 0...@n_a
      for j in 0...@n_b_col
        printf("%+.5f\t", @b[i][j])
      end
      puts ""
    end
  end

  def calculate
    for j in 0...(@n_a)
      partial_pivoting(j)
      # ピボット係数
      coef_p = @a[j][j].to_f

      # ピボット行をピボット係数coefで除算
      # 1. 行列Aについて
      for k in j...@n_a
        @a[j][k] /= coef_p.to_f
      end
      # 2. 行列Bについて
      for k in 0...@n_b_col
        @b[j][k] /= coef_p.to_f
      end

      # ピボット列の掃き出し
      for i in 0...@n_a
        next if i == j
        coef_r = @a[i][j].to_f
        # 各行減算
        # 1.行列Aについて
        for k in j...@n_a
          @a[i][k] -= @a[j][k].to_f * coef_r
        end
        # 2.行列Bについて
        for k in 0...@n_b_col
          @b[i][k] -= @b[j][k].to_f * coef_r
        end
      end
    end
  end

  private
  def partial_pivoting(j)
    #最大値を記録
    amax = @a[j][j].to_f.abs
    imax = j

    for i in (j + 1)...@n_a
      r = @a[i][j].to_f.abs
      if r > amax
        # 最大値を変更
        amax = r
        imax = i
      end
    end

    return if j == imax
    # ピボット行と最大値行が異なれば交換
    # 1.行列Aについて
    for k in 0...@n_a
      r = @a[imax][k]
      @a[imax][k] = @a[j][k]
      @a[j][k] = r
    end
    # 2.行列Bについて
    for k in 0...@n_b_col
      r = @b[imax][k]
      @b[imax][k] = @b[j][k]
      @b[j][k] = r
    end
  end
end

実行

$ ruby calc.rb

行列を定義したCSVファイル(例： a.csv, b.csv)を指定することもできます。

$ ruby calc.rb a.csv b.csv

一応、CSVファイルの例としてa.csvとb.csvの中身を下記に載せておきます。

CSVファイル

a.csv

-1, 2, 2
2, 3, -1
1, 1, -1

b.csv

1, 2
-3, 2
-2, 2

結果

$ ruby calc.rb a.csv b.csv
A =
-1.00000    +2.00000    +2.00000
+2.00000    +3.00000    -1.00000
+1.00000    +1.00000    -1.00000
B =
+1.00000    +2.00000
-3.00000    +2.00000
-2.00000    +2.00000
X (satisfying AX=B) =
+1.00000    -6.00000
-1.00000    +3.00000
+2.00000    -5.00000
finish (no error)
all finish

おわりに

ガウスジョルダン消去法アルゴリズムを実装したプログラムをご紹介しましたが、以前の投稿のガウスの消去法アルゴリズムのプログラムとほとんどの処理（ここでは、メソッド）が共通しています。

ですので、次回は、これら二つのプログラムの共通部分を親クラスに抜き出し、これらのプログラムを使って、オブジェクト指向プログラミングの三大要素であるカプセル化、継承、ポリモーフィズムの理解に繋がるプログラム設計をご紹介したいと思います。

【ガウス・ジョルダン消去法(Gauss-Jordan elimination method)】Rubyによるガウス・ジョルダン消去法アルゴリズム