Solve Linear System and Invert Matrix

From Directorforum Collaboration Wiki

Jump to: navigation, search

Taken from a post on Direct-L by Danny Kodicek

I was just looking over some old code written by someone that did some work for us a while back (naming no names, but I was surprised he didn't know better  ;) ) and happened to notice that they were using a really inefficient piece of code for solving a set of linear simultaneous equations. I've rewritten the code but I thought as this is the kind of thing that lots of people need and few people know, some of you might find it useful. It's a fairly standard algorithm and it's very very fast compared to inverting the equation matrix (actually, technically it's also a fast algorithm for inverting the matrix! See below). I'm sure it could be optimised further if you wanted to use it a lot.

To use it, convert equations like this:

2x + 3y = 8
3x - 4y = -5

into a list like this:

[[2,3,8], [3,-4,-5]]

Each element represents one equation, with the n coefficients first, and the constant at the end.

The function returns the solution as a list, with each unknown in order. So in this case it would return [1,2], since the solution is x=1, y=2. 

on solveLinearSystem tMatrix
 -- pass this a list of n rows, each of which has n+1 elements.
 -- NB: this doesn't do any error checking, so make sure you pass it the right kind of data
 
 -- because it needs to modify the matrix, it makes a duplicate first
 tM = []
 tN = tMatrix.count
 tMap = []
 repeat with i = 1 to tN
   tM.add(tMatrix[i].duplicate())
 end repeat
 repeat with i = 1 to tN
   tRow = tM[i]
   -- find the first non-zero element in the row
   repeat with j = 1 to tN
     tEl = tRow[j]
     if tEl<>0 then
       tMap[i] = j
       -- divide by the element to make it 1
       tRow = tRow/tEl
       tM[i] = tRow
       -- subtract multiples of the row from all other rows
       repeat with k = 1 to tN
         if k<>i then
           tE = tM[k][j]
           if tE<>0 then
             tM[k] = tM[k] - tE*tRow
           end if
         end if
       end repeat
       exit repeat
     end if
   end repeat
   if j>tN then return #cannotDivideByZero -- indeterminate solution
 end repeat
 
 -- read out the solution from the last column of the matrix
 tRes = []
 repeat with i = 1 to tN
   tRes.add(tM[tMap[i]][tN+1])
 end repeat
 return tRes
end

The same algorithm can be used to invert a matrix. This time, we apply it simultaneously to the rows of an identity matrix: 

on invertMatrix tMatrix
 -- pass this a list of n rows, each of which has n elements.
 -- NB: this doesn't do any error checking, so make sure you pass it the right kind of data
 
 -- because it needs to modify the matrix, it makes a duplicate first
 tM = []
 tN = tMatrix.count
 tMap = []
 repeat with i = 1 to tN
   tRow = tMatrix[i].duplicate()
   repeat with j = 1 to tN
     if j = i then tRow.add(1)
     else tRow.add(0)
   end repeat
   tM.add(Row)
 end repeat
 repeat with i = 1 to tN
   tRow = tM[i]
   -- find the first non-zero element in the row
   repeat with j = 1 to tN
     tEl = tRow[j]
     if tEl<>0 then
       tMap[i] = j
       -- divide by the element to make it 1
       tRow = tRow/tEl
       tM[i] = tRow
       -- subtract multiples of the row from all other rows
       repeat with k = 1 to tN
         if k<>i then
           tE = tM[k][j]
           if tE<>0 then
             tM[k] = tM[k] - tE*tRow
           end if
         end if
       end repeat
       exit repeat
     end if
   end repeat
   if j>tN then return #cannotDivideByZero -- matrix has no inverse
 end repeat
 
 -- read out the solution from the second half of the matrix
 tRes = []
 repeat with i = 1 to tN
   tRow = tM[tMap[i]]
   tR = []
   repeat with j = tN+1 to tN*2
     tR.add(tRow[j])
   end repeat
   tRes.add(tR)
 end repeat
 return tRes
end
Retrieved from "http://collab.directorforum.com/Solve_Linear_System_and_Invert_Matrix"
Personal tools