Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active June 1, 2019 05:44
Show Gist options
  • Save komasaru/ae84e81f9498c3f2e120f117f5e45303 to your computer and use it in GitHub Desktop.
Save komasaru/ae84e81f9498c3f2e120f117f5e45303 to your computer and use it in GitHub Desktop.
Ruby script to solve simultaneous equations by LU-decomposition(outer-product form).
#! /usr/local/bin/ruby
#*********************************************
# 連立1次方程式の解法 ( LU 分解(外積形式ガウス法) )
#*********************************************
#
class SleLu
#A = [
# [1, 2],
# [4 ,5]
#]
#B = [3, 6]
A = [
[3.0, 1.0, 1.0, 1.0],
[1.0, 2.0, 1.0, -1.0],
[1.0, 1.0, -4.0, 1.0],
[1.0, -1.0, 1.0, 1.0]
]
B = [0.0, 4.0, -4.0, 2.0]
def initialize
@a = Marshal.load(Marshal.dump(A)) # 行列 A
@n = @a.size # 次元の数
@x = Array.new(@n, 0.0) # 解
end
# 主処理
def exec
exit if @n < 1
puts "A ="
display_mtx(@a)
puts "b ="
display_vec(B)
lu_decompose
#puts "(LU) ="
#display_mtx(@a)
solve
puts "x ="
display_vec(@x)
rescue => e
$stderr.puts "[#{e.class}] #{e.message}"
e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
end
private
# 行列 A の LU 分解
def lu_decompose
0.upto(@n - 1) do |k|
if @a[k][k] == 0
puts "Can't decompose A to LU ..."
exit
end
tmp = 1.0 / @a[k][k]
(k + 1).upto(@n - 1) { |i| @a[i][k] *= tmp }
(k + 1).upto(@n - 1) do |j|
tmp = @a[k][j]
(k + 1).upto(@n - 1) { |i| @a[i][j] -= @a[i][k] * tmp }
end
end
rescue => e
raise
end
# 連立1次方程式の解法
# * @a は LU 分解済み
# * L の対角成分が 1
def solve
# 前進代入
# * Ly = b から y を計算
y = Array.new(@n)
0.upto(@n - 1) do |i|
sum = B[i]
0.upto(i - 1) { |j| sum -= @a[i][j] * y[j] }
y[i] = sum
end
# 後退代入
# * Ux = y から x を計算
(@n - 1).downto(0) do |i|
sum = y[i]
(i + 1).upto(@n - 1) { |j| sum -= @a[i][j] * @x[j] }
@x[i] = sum / @a[i][i]
end
rescue => e
raise
end
def display_mtx(mtx)
mtx.each do |row|
puts row.map { |a| "%8.2f" % a }.join
end
rescue => e
raise
end
def display_vec(vec)
puts vec.map { |a| "%8.2f" % a }.join
rescue => e
raise
end
end
SleLu.new.exec if __FILE__ == $0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment