【ガウスの消去法(Gaussian elimination method)】Rubyによるガウスの消去法アルゴリズム(構造化プログラミング/ オブジェクト指向プログラミング)

はじめに

先日、こちらの記事にFortranによるガウスの消去法アルゴリズムをまとめました。
今回は、そのプログラムをrubyで書き換えたものを掲載するとします。

注：ガウスの消去法に関しては上記リンク先に簡単にまとめています。

プログラム

Fortranで記述したプログラムのようにトップダウン式にプログラムを組む1. 構造化プログラミング（全体の流れを考えてから、個々の機能を細分化していく方法）とボトムアップ式にプログラムを組む2. オブジェクト指向プログラミング（必要な部品をオブジェクトという概念で捉えて、そのオブジェクトをモジュールとして組み合わせて全体を構成する方法）のそれぞれで組んだプログラムを下記に載せています。

1. (構造化プログラミングによる)プログラム

メインプログラム

メインプログラムには手続き（処理）のみを記述しています。

main.rb

require 'csv'
require_relative 'gauss_elim'
require_relative 'disp_matrixes'

=begin
変数の意味
@a: 行列A
@b: 行列B
@n_a: 行列Aのサイズ
@n_b_col: 行列Bの列数

変数@a, @bの初期値
@a = [
  [-1, 2, 2],
  [2, 3, -1],
  [1, 1, -1]
]
@b = [
  [1, 2],
  [-3, 2],
  [-2, 2]
]
=end

# 1. 変数の設定
if ARGV[0] && ARGV[1]
  @a = CSV.readlines("./#{ARGV[0]}")
  @b = CSV.readlines("./#{ARGV[1]}")
else
  # 初期値
  @a = CSV.readlines("./a.csv")
  @b = CSV.readlines("./b.csv")
end

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

# 2. 行列A, Bを表示
display_matrixA_and_matrixB(@a, @b, @n_a, @n_b_col)

# 3. ガウスの消去法を適用
gauss_elim(@a, @b, @n_a, @n_b_col)

# 4. 行列X(解)を表示
display_matrixX(@b, @n_a, @n_b_col)

サブルーチン1

サブルーチン1に行列A、B、Xを表示する関数を定義しています。

disp_matrixes.rb

def display_matrixA_and_matrixB(a, b, n_a, n_b_col)
  # 行列Aを出力
  puts "A ="
  n_a.times do |i|
    n_a.times {|j| printf("%+5f\t", a[i][j])}
    puts ""
  end

  # 行列Bを出力
  puts "B ="
  n_a.times do |i|
    n_b_col.times {|j| printf("%+5f\t", b[i][j])}
    puts ""
  end
end

def display_matrixX(b, n_a, n_b_col)
  # 行列X(解)を出力
  puts "X (satisfying AX=B) ="
  n_a.times do |i|
    n_b_col.times {|j| printf("%+5f\t", b[i][j])}
    puts ""
  end
end

サブルーチン2

サブルーチン2にガウスの消去法関数を定義しています。

gauss_elim.rb

require_relative 'ppivot'

# ガウスの消去法
def gauss_elim(a, b, n_a, n_b_col)
  # 前進消去
  (n_a - 1).times do |j|
    partial_pivoting(a, b, n_a, n_b_col, j)
    # ピボット列の掃き出し
    (j + 1).upto(n_a - 1) do |i|
      r = a[i][j].to_f / a[j][j].to_f
      # 各行減算
      # 1.行列Aについて
      j.upto(n_a - 1) do |k|
        a[i][k] = a[i][k].to_f - a[j][k].to_f * r
      end
      # 2.行列Bについて
      0.upto(n_b_col - 1) do |k|
        b[i][k] = b[i][k].to_f - b[j][k].to_f * r
      end
    end
  end

  # 後退代入
  (n_b_col).times do |j|
    (n_a - 1).downto(0) do |i|
      r = b[i][j].to_f
      # シグマ演算
      (i + 1).upto(n_a - 1) do |k|
        # ここでのrはx
        r = r - a[i][k].to_f * b[k][j].to_f
      end
      b[i][j] = r / a[i][i].to_f
    end
  end
end

サブルーチン3

サブルーチン3に部分ピボット選択関数を定義しています。

ppivot.rb

# 部分ピボット選択
def partial_pivoting(a, b, n_a, n_b_col, j)
  #最大値を記録
  amax = a[j][j].to_f.abs
  imax = j

  (j + 1).upto(n_a - 1) do |i|
    r = a[i][j].to_f.abs
    if r > amax
    # 最大値を変更
      amax = r
      imax = i
    end
  end

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

実行

$ ruby main.rb

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

$ ruby main.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 main.rb a.csv b.csv
A =
-1.000000   +2.000000   +2.000000
+2.000000   +3.000000   -1.000000
+1.000000   +1.000000   -1.000000
B =
+1.000000   +2.000000
-3.000000   +2.000000
-2.000000   +2.000000
X (satisfying AX=B) =
+1.000000   -6.000000
-1.000000   +3.000000
+2.000000   -5.000000

2. （オブジェクト指向プログラミングによる）プログラム

メインプログラム

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

calc.rb

require './gauss_elimination.rb'

if __FILE__ == $0
  begin
    # 計算クラスインスタンス化
    obj = GaussElimination.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}" }
  end
end

モジュール1

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

注：今回はモジュールは一つのみですが、複数モジュールを準備できるということを示すために、便宜上モジュール1としています。

gauss_elimination.rb

require 'csv'

class GaussElimination
=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 ="
    @n_a.times do |i|
      @n_a.times {|j| printf("%+.5f\t", @a[i][j])}
      puts ""
    end

    # 行列Bを出力
    puts "B ="
    @n_a.times do |i|
      @n_b_col.times {|j| printf("%+.5f\t", @b[i][j])}
      puts ""
    end
  end

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

  def calculate
    # 前進消去
    (@n_a - 1).times do |j|
      partial_pivoting(j)
      # ピボット列の掃き出し
      (j + 1).upto(@n_a - 1) do |i|
        r = @a[i][j].to_f / @a[j][j].to_f
        # 各行減算
        # 1.行列Aについて
        j.upto(@n_a - 1) do |k|
          @a[i][k] = @a[i][k].to_f - @a[j][k].to_f * r
        end
        # 2.行列Bについて
        0.upto(@n_b_col - 1) do |k|
          @b[i][k] = @b[i][k].to_f - @b[j][k].to_f * r
        end
      end
    end

    # 後退代入
    (@n_b_col).times do |j|
      (@n_a - 1).downto(0) do |i|
        r = @b[i][j].to_f
        # シグマ演算
        (i + 1).upto(@n_a - 1) do |k|
          # ここでのrはx
          r = r - @a[i][k].to_f * @b[k][j].to_f
        end
        @b[i][j] = r / @a[i][i].to_f
      end
    end
  end

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

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

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

実行

$ ruby calc.rb

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

$ ruby calc.rb c.csv d.csv

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

CSVファイル

c.csv

2, 4, 2
2, 0, 1
1, 1, 1

d.csv

6, 1, 0, 0
2, 0, 1, 0
1, 0, 0, 1

結果

$ ruby calc.rb c.csv d.csv
A =
+2.000000   +4.000000   +2.000000
+2.000000   +0.000000   +1.000000
+1.000000   +1.000000   +1.000000
B =
+6.000000   +1.000000   +0.000000   +0.000000
+2.000000   +0.000000   +1.000000   +0.000000
+1.000000   +0.000000   +0.000000   +1.000000
X (satisfying AX=B) =
+3.000000   +0.500000   +1.000000   -2.000000
+2.000000   +0.500000   -0.000000   -1.000000
-4.000000   -1.000000   -1.000000   +4.000000

おわりに

今回、構造化プログラミンング、オブジェクト指向プログラミングそれぞれでガウスの消去法プラグラムを組んでみましたが、臨機応変にどちらの手法を採用すべきか考えていきたいものです。
いずれにしても、メソッド（関数）、クラスの中身はブラックボックスとして捉えることができるので、メインプログラムを見たときに、何をしているのか直感的に理解できるように記述することが大切です。

こちらの記事が役にたったという方は、いいね、よろしくお願いします(^^)