ガウスの消去法を使って、解を求めるプログラム。不定解・解なしの場合も対応。

gauss.rb

require 'rational'


class AugmentedMatrix

	def initialize a, b
		@matrix = []

		@i_size = a.size
		@j_size = a[0].size

		@i_size.times do |i|	
			vector = []
			
			if @j_size == a[i].size then
				@j_size.times do |j|
					vector.push a[i][j]			
				end
				vector.push b[i]
			
				@matrix.push vector
			
			else
				raise "Argument Error: size of matrix is irregular."

			end

		end
		
		@gaussed = false
	end


	def inverse a
		#return 1.0 / a
		return Rational(1, 1) / a
	end
	
	def zero
		#return 0.0 
		return Rational(0, 1)
	end

	def one 
		#return 1.0
		return Rational(1, 1)
	end


	def print_matrix
		puts "### matrix: "
			
		@matrix.size.times do |i|
			if i == 0
				print "[ [ " 
			else
				print "  [ "
			end
			@matrix[i].size.times do |j|
				a_ij = @matrix[i][j]
				print "#{a_ij}, "	
			end
			if i == (@matrix.size-1)
				puts "] ]"
			else
				puts "]"
			end
		end

	end


	def print_solutions x
		puts "### solutions: "
		
		if x.size > 0 then
			x[0].size.times do |j|
				print "x#{j+1} = "
				
				x.size.times do |i|
					c = " + c#{i}*"
					if i==0 
						c = ""
					end
					print "#{c}(#{x[i][j]})"
				end
	
				puts ""
			end
		else
			puts "# no solutions."
			
		end
	
	end




	# 行基本変形 1: 
	# (i,j) 行を入れ替える
	def swap_line i, j
		temp = @matrix[i]
		@matrix[i] = @matrix[j]
		@matrix[j] = temp
	
	end


	# 行基本変形 2:
	# i 行をスカラー倍する
	def scale(i, scalar)
		vector = @matrix[i]
		vector.size.times do |j|
			vector[j] = vector[j] * scalar
		end
	end
	
	
	# 行基本変形 3: 
	# i 行を j 行に scalar 倍して足し合わせる
	def add(i, j, scalar)
		@matrix[j].size.times do |k|
			@matrix[j][k] = @matrix[j][k] + @matrix[i][k] * scalar
		end
		
	end




	# ガウスの消去法を実行する
	# 以降, 行を i, 列を j とします
	def gauss
		current_i = 0
		(@matrix[0].size - 1).times do |current_j|
		
			# ピボットを決定（すなわち、先頭行が非ゼロな行ベクトルをみつける
			i_pivot = current_i
				
			while i_pivot < @matrix.size && @matrix[i_pivot][current_j] == zero
				i_pivot += 1
			end

		
			# ピボット行が存在する
			if i_pivot < @matrix.size then
							
				# ピボットを一番上に移動
				self.swap_line( i_pivot, current_i )
						
				# ピボット行の逆元
				inverse_pivot = inverse( @matrix[current_i][current_j] )
						
				# 逆元をかける
				self.scale( current_i, inverse_pivot )
										
				# スカラー倍して足し合わせる
				@matrix.size.times do |i|
					if i != current_i 
						self.add( current_i, i, - @matrix[i][current_j])
					end
				end

				current_i += 1
			end
			
			
			puts "->"
			self.print_matrix
		
		end

		@gaussed = true
	end



	# 後退代入
	def solve
		j_max = @matrix[0].size - 1
		x = [Array.new( j_max, nil )]
		
	
		@matrix.size.times do |k|
			i = @matrix.size - 1 - k
			
					
			# 階段行列の左端を見つける ("1" を探索)
			j = 0
			while j < j_max && @matrix[i][j] != one
				j += 1
			end
			
			if j < j_max
				# 1 をみつけた, 代入開始
				current_j = j
				x[0][current_j] = @matrix[i][j_max]			
				
				j += 1
				while j < j_max
					# nil であれば不定解なのでベクトルを追加
					if x[0][j] == nil
						x[0][j] = zero
						
						x.push( Array.new( j_max, zero ) )
						x[x.size-1][j] = one

					end
				
					x.size.times do |l|
						x[l][current_j] -= @matrix[i][j] * x[l][j]
					end
					
					j += 1
				end
				
			else 
				if @matrix[i][j_max] != zero then
					# 解なし
					return []
				end
			
			end
	
		end
	
		x
	end

end




if __FILE__ == $0 then

def test_a
	
	a = [ [1, 3, 1], [1, 1, -1], [3, 11, 5] ]
	b = [9, 1, 35]
	
	matrix = AugmentedMatrix.new a,b
	matrix.print_matrix
	matrix.gauss

	puts ""

	x = matrix.solve
	matrix.print_solutions x
	
	puts ""
	puts ""

end
test_a




def test_b
	a = [ [2, 4, 2, 2], [4, 10, 3, 3], [2, 6, 1, 1], [3, 7, 1, 4] ]
	b = [8, 17, 9, 11]

	matrix = AugmentedMatrix.new a,b
	
	matrix.print_matrix
	matrix.gauss

	puts ""

	x = matrix.solve
	matrix.print_solutions x
	
	puts ""
	puts ""

end
test_b


def test_c
	# null space 
	a = [ [2, 4, 2, 2, 1, 1], [4, 10, 3, 3, 1, 1], [2, 6, 1, 1, 1, 1], [3, 7, 1, 4, 1, 1] ]
	b = [0, 0, 0, 0]

	matrix = AugmentedMatrix.new a,b
	
	matrix.print_matrix
	matrix.gauss

	puts ""

	x = matrix.solve
	matrix.print_solutions x
	
	puts ""
	puts ""

end
test_c



def test_d
	# no solutions
	
	a = [ [1, 0, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 0] ]
	b = [0, 0, 1]

	matrix = AugmentedMatrix.new a,b
	
	matrix.print_matrix
	matrix.gauss

	puts ""

	x = matrix.solve
	matrix.print_solutions x
	
	puts ""
	puts ""

end
test_d

end