Class: Algebra::SquareMatrix

Inherits:
MatrixAlgebra show all
Defined in:
lib/algebra/linear-algebra.rb,
lib/algebra/matrix-algebra.rb,
lib/algebra/gaussian-elimination.rb

Defined Under Namespace

Classes: EigenSystem

Constant Summary collapse

Matrices = 
{}

Class Method Summary collapse

Instance Method Summary collapse

Methods inherited from MatrixAlgebra

#*, *, **, #+, #-, #==, #[], [], #[]=, #coerce, #cofactor, #cofactor_matrix, #collect_column, collect_column, collect_ij, #collect_ij, collect_row, #collect_row, #column, #columns, convert_from, #convert_to, #csize, #diag, #display, #display_by_latex, dsum, #dsum, dsum_c, dsum_r, #dup, #e_deg, #each, #each_column, #each_i, each_ij, each_index, #each_index, #each_j, #each_row, #flatten, #i2o, #initialize, #inspect, #jordan_form, #jordan_form_info, #latex, #map, matrix, #matrix, minor, #minor, #old_mul, #orthogonal_basis, #rank, #replace, #row, #row!, #rows, #rsize, #set_column, #set_row, #simplify, #sizes, #to_ary, #to_quint, #to_s, #to_triplet, #transpose, transpose, #types, #vectors, zero

Methods included from AlgebraCreator

#create, #superior?, #sysvar, #wedge

Methods included from GaussianElimination

#divide_c, #divide_c!, #divide_r, #divide_r!, #kernel_basis, #left_eliminate!, #left_eliminate_euclidian!, #left_inverse, #left_sweep, #mix_c, #mix_c!, #mix_r, #mix_r!, #multiply_c, #multiply_c!, #multiply_r, #multiply_r!, #sswap_r!, #step_matrix?, #swap_c, #swap_c!, #swap_r, #swap_r!

Methods included from ElementaryDivisor

#_coeff, #e_diagonalize, #e_diagonalize!, #e_inverse, #el_sweep!, #el_sweep_old!, #elementary_divisor, factorize, #horizontal_sweep!, #sq_find_nondiv, #vertical_sweep!

Methods included from Enumerable

#all?, #any?, #collecti, #sum

Methods included from AlgebraBase

#*, #**, #+, #+@, #-, #-@, #==, append_features, #coerce, #devide?, #ground, #ground=, #regulate, #unit?, #unity, #unity?, #zero, #zero?

Methods included from Orthogonalization

#normalize, #orthogonalize

Constructor Details

This class inherits a constructor from Algebra::MatrixAlgebra

Class Method Details

.const(x) ⇒ Object 



715
716
717
 
# File 'lib/algebra/matrix-algebra.rb', line 715

def self.const(x)
  matrix{|i, j| i == j ? ground.regulate(x) : ground.zero }
end

.create(ground, n, m = nil) ⇒ Object 



731
732
733
734
735
736
737
738
 
# File 'lib/algebra/matrix-algebra.rb', line 731

def self.create(ground, n, m = nil)
  if m && n != m
    raise "type error to create SquareMatrix"
  end
  klass = super(ground, n, n)
  klass.sysvar(:size, n)
  klass
end

.det ⇒ Object 



759
760
761
762
763
 
# File 'lib/algebra/matrix-algebra.rb', line 759

def self.determinant(aa)
  r = create(Integer, aa.size)
  m = r.new(aa)
  m.determinant
end

.determinant(aa) ⇒ Object 



752
753
754
755
756
 
# File 'lib/algebra/matrix-algebra.rb', line 752

def self.determinant(aa)
  r = create(Integer, aa.size)
  m = r.new(aa)
  m.determinant
end

.matrices ⇒ Object 



723
724
725
 
# File 'lib/algebra/matrix-algebra.rb', line 723

def self.matrices
  Matrices
end

.regulate(x) ⇒ Object 



742
743
744
745
746
747
748
749
750
 
# File 'lib/algebra/matrix-algebra.rb', line 742

def self.regulate(x)
  if x.is_a? superclass #SquareMatrix
    x
  elsif y = ground.regulate(x)
    const(y)
  else
    super
  end
end

.symmetric(*a) ⇒ Object 



49
50
51
52
53
54
55
56
 
# File 'lib/algebra/linear-algebra.rb', line 49

def self.symmetric(*a)
  k = size * (size + 1) / 2
  raise "the size of #{a.inspect} must be <= #{k}" if a.size > k
  matrix do |i, j|
    d = (i - j).abs
    a[(i < j ? i : j) + (2 * size + 1 - d) * d / 2] || ground.zero
  end
end

.unity ⇒ Object 



711
712
713
 
# File 'lib/algebra/matrix-algebra.rb', line 711

def self.unity
  matrix{|i, j| i == j ? ground.unity : ground.zero }
end

Instance Method Details

#/(other) ⇒ Object 



762
763
764
765
766
767
768
769
770
771
772
 
# File 'lib/algebra/matrix-algebra.rb', line 762

def /(other)
  #It might be better to use regulate ...
  case other
#    when ground, Numeric
#      matrix{|i, j| self[i, j] / other}
  when self.class.superclass #SquareMatrix
    self * other.inverse
  else
    super
  end
end

#_char_matrix(mop) ⇒ Object 



839
840
841
842
 
# File 'lib/algebra/matrix-algebra.rb', line 839

def _char_matrix(mop)
  # assume mop.ground.ground == self.ground
  mop.unity * mop.ground.var - self
end

#char_matrix(polyring) ⇒ Object 



844
845
846
847
 
# File 'lib/algebra/matrix-algebra.rb', line 844

def char_matrix(polyring)
  matrix_over_poly = Algebra.SquareMatrix(polyring, size)
  _char_matrix(matrix_over_poly)
end

#char_polynomial(polyring) ⇒ Object 



849
850
851
852
 
# File 'lib/algebra/matrix-algebra.rb', line 849

def char_polynomial(polyring)
  matrix_over_poly = Algebra.SquareMatrix(polyring, size)
  _char_matrix(matrix_over_poly).determinant
end

#char_polynomial0(obj = "x") ⇒ Object 



831
832
833
834
835
836
837
 
# File 'lib/algebra/matrix-algebra.rb', line 831

def char_polynomial0(obj = "x")
  require "algebra/polynomial"
  a = Algebra.Polynomial(ground, obj)
  k = Algebra.SquareMatrix(a, size)
  (k.unity * a.var - self).determinant
  #    (k.unity * a.var - k.matrix{|i, j| self[i,j]}).determinant
end

#const(x) ⇒ Object 



719
720
721
 
# File 'lib/algebra/matrix-algebra.rb', line 719

def const(x)
  self.class.const(x)
end

#det ⇒ Object 



787
788
789
790
791
792
793
794
795
796
797
 
# File 'lib/algebra/matrix-algebra.rb', line 787

def determinant
  sum = ground.zero
  perm do |idx|
    product = ground.unity
    idx.each_with_index do |i, j|
	product *= self[i, j]
    end
    sum += product * sign(idx)
  end
  sum
end

#determinant ⇒ Object 



774
775
776
777
778
779
780
781
782
783
784
 
# File 'lib/algebra/matrix-algebra.rb', line 774

def determinant
  sum = ground.zero
  perm do |idx|
    product = ground.unity
    idx.each_with_index do |i, j|
	product *= self[i, j]
    end
    sum += product * sign(idx)
  end
  sum
end

#determinant! ⇒ Object 



786
787
788
789
790
791
792
793
794
795
796
 
# File 'lib/algebra/matrix-algebra.rb', line 786

def determinant
  sum = ground.zero
  perm do |idx|
    product = ground.unity
    idx.each_with_index do |i, j|
	product *= self[i, j]
    end
    sum += product * sign(idx)
  end
  sum
end

#determinant_by_elimination ⇒ Object

ground ring must be a field



260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
 
# File 'lib/algebra/gaussian-elimination.rb', line 260

def determinant_by_elimination
  m = dup
  det = ground.unity
  each_j do |d|
    if i = (d...size).find{|i| !m[i, d].zero?}
      if i != d
        m.sswap_r!(d, i)
        det = -det
      end
      c = m[d, d]
      det *= c
      (d+1...size).each do |i0|
        r = m.row!(i0)
        q = ground.unity * m[i0, d] / c #this lets the entries be in ground
        (d+1...size).each do |j0|
          r[j0] -= q * m[d, j0]
        end
      end
    else
      return ground.zero
    end
  end
  det
end

#determinant_by_elimination_euclidian ⇒ Object

ground ring must be an euclidian domain



286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
 
# File 'lib/algebra/gaussian-elimination.rb', line 286

def determinant_by_elimination_euclidian#_new
  m = dup
  det = ground.unity
  each_j do |d|
    if i = (d...size).find{|i| !m[i, d].zero?}
      if i != d
        m.sswap_r!(d, i)
        det = -det
      end
      (d+1...size).each do |i0|
        r = m.row!(i0)
        q = m[i0, d] / m[d, d]
        (d...size).each do |j0|; r[j0] -=  q * m[d, j0]; end
        unless m[i0, d].zero?
          m.sswap_r!(d, i0)
          det = -det
          redo
        end
      end
      det *= m[d, d]
    else
      return ground.zero
    end
  end
  det
end

#determinant_by_elimination_euclidian_old ⇒ Object 



249
250
251
252
253
254
255
256
257
 
# File 'lib/algebra/gaussian-elimination.rb', line 249

def determinant_by_elimination_euclidian_old
  m = dup
  inv, k = m.left_eliminate_euclidian!
  s = ground.unity
  (0...size).each do |i|
    s *= m[i, i]
  end
  s / k
end

#determinant_by_elimination_old ⇒ Object

def inverse_euclidian; left_inverse_euclidian; end



239
240
241
242
243
244
245
246
247
 
# File 'lib/algebra/gaussian-elimination.rb', line 239

def determinant_by_elimination_old
  m = dup
  inv, k = m.left_eliminate!
  s = ground.unity
  (0...size).each do |i|
    s *= m[i, i]
  end
  s / k
end

#diagonalize ⇒ Object 



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
106
107
108
109
110
111
112
113
114
115
116
 
# File 'lib/algebra/linear-algebra.rb', line 69

def diagonalize
  chp = char_polynomial(Algebra.Polynomial(ground, 't'))

  facts = chp.factorize
  mdf, mods, fas, roots, proots = chp.decompose(facts)
  nu = Algebra.SquareMatrix(mdf, size)

  espaces = {}
  evalues = []
  evectors = []
  roots0 = []
  k = 0
  facts.each do |f, _n|
    dim = f.deg
    if dim <= 0 # 0 might be occurred
      next
    elsif dim == 1
      evs = [-f[0] / f[1]]
    else
      evs = roots[k, dim] # .reverse
      roots0.push [f, evs]
    end
    k += dim

    evs.each do |evalue|
      ks = nu.const(evalue) - self
      kb = ks.kernel_basis
      espaces[evalue] = kb

      kb.each do |v|
        evalues.push evalue
        evectors.push v
      end
    end
  end

  raise "can't Diagonalize" if evalues.size < size
  ttype = Algebra.SquareMatrix(mdf, size)
  r = ttype.collect_column { |i| evectors[i] }
  #    :extfield, :roots, :tmatrix, :evalues, :addelms, :evectors,
  #                   :espaces, :chpoly, :facts
  e = EigenSystem.new(mdf, roots0, r, evalues, proots, evectors,
                      espaces, chp, facts)
  def e.to_ary
    to_a
  end
  e
end

#initialzie(array, conv = false) ⇒ Object 



727
728
729
 
# File 'lib/algebra/matrix-algebra.rb', line 727

def initialzie(array, conv = false)
  super
end

#inverse ⇒ Object 



789
790
791
792
793
794
795
796
797
798
799
800
801
802
 
# File 'lib/algebra/matrix-algebra.rb', line 789

def inverse
  #gaussian-elimination.rb
  if ground.field? &&
      (!ground.respond_to?(:reducible) || ground.reducible)
    inverse_field
  else
    a = adjoint
    d = determinant
    unless d.unit?
      warn("#{a}/#{d} may not be done.")
    end
    a / d
  end
end

#inverse_field ⇒ Object 



230
231
232
233
234
235
236
 
# File 'lib/algebra/gaussian-elimination.rb', line 230

def inverse_field
  n, k, pi = dup.left_eliminate!
  if pi != size
    raise "Error: invert non invertible SquareMatrix"
  end
  n
end

#perm(n = size, stack = []) ⇒ Object 



817
818
819
820
821
822
823
824
825
826
827
828
829
 
# File 'lib/algebra/matrix-algebra.rb', line 817

def perm(n = size, stack = [])
  (0...n).each do |x|
    unless stack.include? x
	stack.push x
	if stack.size < n
 perm(n, stack) do |y|; yield y; end
	else
 yield stack
	end
	stack.pop
    end
  end
end

#sign(a) ⇒ Object 



804
805
806
807
808
809
810
811
812
813
814
815
 
# File 'lib/algebra/matrix-algebra.rb', line 804

def sign(a)
  n = a.size
  a = a.dup
  b = (0...n).collect{|i| a.index(i)}
  s = 1
  (0...(n-1)).each do |i|
    if (j = a[i]) != i
	a[b[i]], b[j], s = j, b[i], -s
    end
  end
  s
end

#size ⇒ Object 



740
 
# File 'lib/algebra/matrix-algebra.rb', line 740

def size; self.class.size; end

#solve_eigen_value_problem ⇒ Object 



118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
 
# File 'lib/algebra/linear-algebra.rb', line 118

def solve_eigen_value_problem
  puts 'A = '; display; puts

  e = diagonalize
  puts "Charactoristic Poly.: #{e.chpoly} => #{e.facts}"; puts
  puts 'Algebraic Numbers:'
  e.roots.each do |po, rs|
    puts "#{rs.join(', ')} : roots of #{po} == 0"
  end; puts

  puts 'EigenSpaces: '
  e.evalues.uniq.each do |ev|
    puts "W_{#{ev}} = <#{e.espaces[ev].join(', ')}>"
  end; puts

  puts 'Trans. Matrix:'
  puts 'P ='
  e.tmatrix.display; puts

  puts 'P^-1 * A * P = '; (e.tmatrix.inverse * self * e.tmatrix).display; puts
end