1

When I tried to use the function "solve" in gem NMatrix, I find wrong result...

I am trying to solve the linear system a*x = b where

a = [[1,  0, 0], 
     [4, -5, 1],
     [0,  0, 1]]

and 

b = [0, 0, 729]

Which should give the answer

x = [0, 145.8, 729]

I call the solve function doing a.solve(b) (as explained at http://www.rubydoc.info/gems/nmatrix/0.2.1/NMatrix#solve-instance_method) but it returns

x = [0, 1.013500790889141e-30, 1.013500790889141e-30]

Therefore, I have two questions:

  1. Am I doing something wrong?
  2. If not, what would be the other solution (excluding GSL) to solve matrix systems with ruby?

Thank you !

Cary Swoveland
  • 106,649
  • 6
  • 63
  • 100
M. S.
  • 33
  • 1
  • 6
  • 1
    What is that definition for `A` did you actually mean `A = NMatrix.new( [1,0,0],[0,-5,1],[0,0,1], dtype: :float64)`? Please post the actual code you are running because asking if you are doing something wrong cannot be answered based on the existing information – engineersmnky Sep 26 '17 at 15:54
  • The actual code is : `a_mat = NMatrix.eye(n)*(-(c+1)) a_mat[0,0] = 1 a_mat[n-1,n-1] = 1 (1..n-2).each do |i| a_mat[i,i-1]=c a_mat[i, i+1]= 1 end` with c = 4 and `B = NMatrix.new([n,1], Array.new(n-1,0) + [nb_days-1], dtype: :float32)` and n = 3 and nb_days = 730 – M. S. Sep 26 '17 at 15:59
  • It's generally helpful to give an example, but when you do all input objects should be valid Ruby objects. Also, it's helpful to assign a variable to each such object so readers can cut-and-paste and refer to those objects by name in answers and comments (and not having to define them). You have assigned them to constants (names beginning withcapital letters). If you meant them to be variables start each name with a lower-case letter. I edited your answer to make those changes. – Cary Swoveland Sep 26 '17 at 20:59
  • Indeed I had wrotten in mathematical language and not ruby, thank you ! – M. S. Sep 27 '17 at 07:44
  • 1
    @M.S. please update your post with the code not a comment. Also in your comment `#solve` is never called so I have a feeling this isn't your complete code – engineersmnky Sep 27 '17 at 13:30

1 Answers1

0

To answer part 2 of your question, you can use Ruby's built-in Matrix class.

require 'matrix'

m = Matrix[[1, 0, 0], [4, -5, 1], [0, 0, 1]]
  #=> Matrix[[1, 0, 0], [4, -5, 1], [0, 0, 1]]
b = Vector[0, 0, 729]
  #=> Vector[0, 0, 729]

a = m.lup.solve(b).to_a
  #=> [(0/1), (729/5), (729/1)]

If one prefers floats (rather than rational numbers),

a.map(&:to_f)
  #=> [0.0, 145.8, 729.0]

See Matrix#lup and Matrix::LUPDecomposition#solve, which solve the linear system using LU decomposition.

If the system of linear equations has no solution (e.g., [[0,0], [0,0]]*x = b) or an infinite number of solutions (that is, m is singular), Ruby will raise an exception. For example,

Matrix[[0,0],[0,0]].lup.solve [1, 1]
  #=> ExceptionForMatrix::ErrNotRegular: Not Regular Matrix

require matrix makes both the Matrix and Vector classes (and subclasses) available for use.

Cary Swoveland
  • 106,649
  • 6
  • 63
  • 100
  • Thank you for your solution. However, I am not quite confortable with it because it is much more expensive in time (you do two computations instead of one: first compute the inverse, which is already a costly computation, and second compute the multiplication). Furthermore this method is less stable when the matrix m has a determinant close to zero. – M. S. Sep 27 '17 at 08:18
  • I changed my answer to use LU decomposition. – Cary Swoveland Sep 28 '17 at 05:02