Re: Matrix inversion

[Jim Hurley]

Michael,

Thank you for these routines. I think if I can get them to work I 
can find the matrix inverse. My source for numerical work is : 
Numerical analysis, by Richard Burden. It gives only  algorithms and 
not program routines.

Unfortunately I am having trouble with your handlers. I tried a 
simple set of linear equations:

x + 2y = 1
x + y = 0
With the solution: x = -1 and y = 1.

But the routine below, LUBKSB, yields  x = -2 and y = 1.

I  used

aMatrix = 1,2
                1,1
 with  n = 2

and

rhsMatrix =1
                  0
Any ideas?

Jim

At 2:40 PM +1000 3/31/03, Michael J. Lew wrote:
Here are a couple of matrix handlers that I translated from Fortran 
routines published in the excellent Numerical Recipes book by Press 
et al. Not enough entirely to invert matrices (that wasn't what I 
needed them for), but the small additional routine is fully 
documented in a downloadable pdf from 
<http://lib-www.lanl.gov/numerical/bookfpdf.html>.

For anyone who does numerical work I highly recommend the Numerical 
Recipes books.

#--****** LUBKSB ******
# Function Revolutionised from Numerical Recipes by MJL
# --Input is aMatrix, an n by n square matrix and rhsMatrix, an n by 1 vector
# also uses the global indx which is built by LUDCMP
# Returns an n by 1 vector solution
# For use with LUDCMP to solve linear equations or invert a matrix
# Note that the input and output parameters are not fully equivalent 
to the originals

function LUBKSB aMatrix,n,rhsMatrix
  global indx
  put aMatrix into a
  put rhsMatrix into b
  set the numberformat to "0.#########################"
  put 0 into ii
  repeat with i =1 to n
     put indx[i] into ip
    put b[ip] into mySum
    put b[i] into b[ip]
     if ii<>0 then
      repeat with j = ii to i-1
        put  mySum-a[i,j]*b[j] into mySum
       end repeat --j
    else
      if mySum <> 0 then put i into ii
     end if
    put mySum into b[i]
   end repeat --i
  repeat with i = n down to 1
    put b[i] into mySum
    if i<n then
      repeat with j = i+1 to n
        put mySum-a[i,j]*b[j] into mySum
      end repeat --j
     end if
    put mySum/a[i,i] into b[i]
   end repeat --i
  return b
end LUBKSB


#--***** LUDCMP *****
# Function Revolutionised by MJL from Numerical Recipes subroutine ludcmp
# input is an n by n square matrix, aMatrix
# returns aMatrix modified to be an LU decomposed version of itself
# also returns (a global) indx which is a vector of the row perms 
affected by partial pivoting (!)
# For use in conjunction with LUBKSB to solve linear equations or 
invert a matrix
# Note that the input and output parameters are not fully equivalent 
to the originals

function LUDCMP amatrix,n
  global indx
  set the numberformat to "0.#########################"
  put 10^-22 into tiny
  put aMatrix into a
  put 1 into d
  repeat with i = 1 to n
    put 0 into big
    repeat with j=1 to n
      if abs(a[i,j]) >big then put a[i,j] into big
    end repeat --j
    if big =0 then
      answer "Singular matrix in LUDCMP. Procedure failed"
      exit to top
    end if
    put 1/big into vv[i]
  end repeat --i
  repeat with j=1 to n
    repeat with i=1 to j-1
      put a[i,j] into mySum
      repeat with k = 1 to i-1
        put mySum-a[i,k]*a[k,j] into mySum
      end repeat --k
      put mySum into a[i,j]
    end repeat --i
    put 0 into big
    repeat with i=j to n
      put a[i,j] into mySum
      repeat with k=1 to j-1
        put mySum-a[i,k]*a[k,j] into mySum
      end repeat --k
      put mySum into a[i,j]
      put abs(mySum)*vv[i] into dum
      if dum>big then
        put dum into big
        put i into imax
      end if
    end repeat --i
    if j<>imax then
      repeat with k=1 to n
        put a[imax,k] into dum
        put a[j,k] into a[imax,k]
        put dum into a[j,k]
      end repeat --k
      put -d into d
      put vv[j] into vv[imax]
    end if
    put imax into indx[j]
    if a[j,j] =0 then put tiny into a[j,j]
    if j<>n then
      put 1/(a[j,j]) into dum
      repeat with i=j+1 to n
        put dum*a[i,j] into a[i,j]
      end repeat --i
    end if
  end repeat --j
  return a
end LUDCMP


--
Michael J. Lew
Senior Lecturer
Department of Pharmacology
The University of Melbourne
Parkville 3010
Victoria
Australia
Phone +613 8344 8304

**
New email address: [EMAIL PROTECTED]
**
_______________________________________________
use-revolution mailing list
[EMAIL PROTECTED]
http://lists.runrev.com/mailman/listinfo/use-revolution