1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
| #! /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
@a = Marshal.load(Marshal.dump(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
|