Class: Digiproc::Strategies::JacobiStrategy

Inherits:
Object
  • Object
show all
Defined in:
lib/strategies/linear_algebra/jacobi_strategy.rb

Overview

This strategy will solve a system of linear equations using the Jacobi iterative strategy Constructor Inputs: (a_arr, b_arr) correspond to A and B in the equation Ax = B (a should be an nxn 2D array, and b should be a 1D array (even though in the equation it is a column vector))

a = [[4,1,-1],,[1,-3,12]] b = [3,19,31] gs = Digiproc::Strategies::JacobiStrategy.new(a,b) x = gs.calculate # => Matrix[, [2.000021547671973], [3.000054218557957]]

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(a_arr, b_arr) ⇒ JacobiStrategy

Initialize args

a_arr

2D array representing your A matrix

b_arr

1D array representing your B matrix

Where B = Ax defines your series of linear equations



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
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 20

def initialize(a_arr,b_arr)
    # TODO: Raise exception if a_arr is not square and b_arr row_count != a_arr row count
    @b = Matrix.column_vector(b_arr)
    @a = Matrix.rows(a_arr, true)
    d_arr, l_arr, u_arr = [],[],[]
    num_cols = @a.column_count
    @a.row_count.times do 
        d_arr << Array.new(num_cols, 0)
        l_arr << Array.new(num_cols, 0)
        u_arr << Array.new(num_cols, 0)
    end

    @a.each_with_index do |el, row, col| 
        if row > col
            l_arr[row][col] = el
        elsif row < col
            u_arr[row][col] = el
        else
            d_arr[row][col] = el
        end
    end
    @d = Matrix.rows(d_arr)
    @l = Matrix.rows(l_arr)
    @u = Matrix.rows(u_arr)
    #TODO: Ensure no zeros on diagonal
    @dinv = @d.map{ |el| el == 0 ? 0 : 1.0 / el }
end

Instance Attribute Details

#a ⇒ Object (readonly)

Returns the value of attribute a.



13
14
15
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 13

def a
  @a
end

#b ⇒ Object (readonly)

Returns the value of attribute b.



13
14
15
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 13

def b
  @b
end

#d ⇒ Object (readonly)

Returns the value of attribute d.



13
14
15
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 13

def d
  @d
end

#dinv ⇒ Object (readonly)

Returns the value of attribute dinv.



13
14
15
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 13

def dinv
  @dinv
end

#l ⇒ Object (readonly)

Returns the value of attribute l.



13
14
15
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 13

def l
  @l
end

#u ⇒ Object (readonly)

Returns the value of attribute u.



13
14
15
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 13

def u
  @u
end

Instance Method Details

#calculate(threshold: 0.001, safety_net: false) ⇒ Object

Iteratively solves the linear system of equations using the jacobi method accepts an optional parameter which is the threshold value x(n+1) - x(n) should achieve before returning. Must be used with a key-value pair ie gs.calculate(threshold: 0.001) default threshold = 0.001 Returns a column vector Matrix => access via matrix_return If run with the option safety_net: true and the equation diverges, performs A_inverse * B to solve ie gs.calculate(safety_net: true)



59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
 
# File 'lib/strategies/linear_algebra/jacobi_strategy.rb', line 59

def calculate(threshold: 0.001, safety_net: false)
    dinv, b, l, u = @dinv, @b, @l, @u
    c = dinv * b
    t = -1 * dinv * (l + u)
    x_n = Matrix.column_vector(Array.new(@a.column_count, 0))
    counter = 0
    loop do 
        x_n_plus_1 = c + t * x_n
        x_difference = (x_n_plus_1 - x_n).map{ |el| el.abs }
        # puts x_n_plus_1
        should_break = !x_difference.find{ |el| el > threshold }
        x_n = x_n_plus_1           
        break if should_break
        counter += 1
        if counter > 10000000 and safety_net 
            return (@a.inv * b).map{ |el| el.to_f}
        end
    end 
    return (safety_net and x_n.find{ |el| el.to_f.nan? }) ? (@a.inv * b).map{ |el| el.to_f} : x_n
end