Re: Eigenvalues &/or eigenvectors, anyone?
Re: Eigenvalues &/or eigenvectors, anyone?
- Subject: Re: Eigenvalues &/or eigenvectors, anyone?
- From: bryan <email@hidden>
- Date: Thu, 08 May 2003 23:42:20 -0400
- Organization: apex
Mark Butcher wrote:
For the mathematicians on this list:
Does anyone have a routine for calculating the eigenvalues &/or
eigenvectors of a matrix?
MarkB
Here are a few helpful elementary matrix functions that might provide a
start.
Several functions are also available in mathematica which may be
accessible depending on your use.
I have reprinted here an old document followed by the scripts that were
sownloaded (ages ago)
I have been unable to find the location that I got them from, but the
website referred to in the ReadMe
still exists and the author, Deivy Petrescu, may be reachable.
I have not tested thes on newer systems (post sys 8.5) (Yeah I remember
when that was new, and Sys 6 too.)
Some of the characters probably will get trashed by the ListServe, so if
anybody wants the original files on this
please contact me directly.
Bryan Kaufman
ReadMe:
Matematica Matrices
Matematica Matrices does the most commom matrix computations. It has
handlers that adds, subtracts, multiply, finds the inverse, computes the
determinant, put a matrix in the upper or lower diagonal form and find
its transpose.
It also allows one to solve a linear system of the form Ax=b, where A is
a square, invertible matrix.
The matrix
/a b c\
\e f g/
is represented as list in the following form:
{{a,e},{b,f},{c,g}}.
So the first element of the list is the first column of the matrix. You
could write the matrix as
{{a,b,c},{e,f,g}},
but you would have to transpose it to send it to the handlers of
Matematica Matrices.
Disclaimer and Copyright Information
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation - version 2.
This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
Public License in file COPYING for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
675 Mass Ave, Cambridge, MA 02139, USA.
__ Note: I'd appreciate an e-mail or postcard if you use or incorporate
any of the modules in your application or work.__
Support
There is no support for this module. However, if you find any bug or you
are experiencing any problems , please write to email@hidden that
I'll try to fix and work on your problem.
Acknowledgements
The idea of this modules came during discussions with Bill Briggs, who
helped me out with the initial project with many corrections, notations
problems and bug fixes.
Faults and omissions are mine alone.
Deivy Petrescu
http://www.dicas.com/
email@hidden
(*
Matrix manipulation using lists.
A list is considered a matrix where each column is an entry in the list.
For instance a 3x2
matrix in this notation is given by {{1,0,0},{1,0,1}} The first column
of this matrix is the first item {1,0,0}. The
second column is (1,0,1). The first line of this matrix is (1,1), the
second is (0,0) and the third is (0,1).
The product handler (Mprod) computes the (matrix) product of two matrix
(lists) THAT CAN BE MULTIPLIED.
To be able to multiply two matrices A and B, the matrices have to be of
the type LxJ and JxK
respectively. This is checked by the handler. The sum handler (Msum)
sums or subtracts two matrices of the same type.
The mbyscalar handler multiplies a matrix A by a scalar s. The transpose
handler transposes a matrix.
It can also be used if the matrix is written in the notation where the
line is an element of the list. In this case, transpose
the given matrix to have any of the handlers below deal with it.
Upperdiag and lowerdiag leaves the given SQUARE matrix in the
equivalent triangular form. Determinant return the determinant of the
SQUARE matrix sent to the handler. The inverse handler
inverts a SQUARE matrix and solvesystem return the solution of the
system Ax=B where A is a square matrix and B is a vector (or column matrix).
Ismatrix checks that the given list is a numerical matrix . It rejects
the empty set and any matrix containing the empty set and matrices.
It returns the type of the matrix (for instance, 3x1).
The Cronecker(x) handler generates the identity matrix in x dimensions
*)
on ismatrix(Matriz)
set lit to 1
set len to length of Matriz
try
if Matriz = {} then error
if class of item 1 of Matriz is list then set lit to length of item 1 of
Matriz
repeat with itm in Matriz
if lit > 1 then
if length of itm ? lit then error
repeat with anitm in itm
if class of anitm = list then error
set anitm to anitm as number
end repeat
end if
if lit = 1 then set itm to itm as number
if lit = 0 then error
end repeat
on error
display dialog "You either sent an empty list, or a list that can not be
coerced into a matrix." buttons {"OK"} default button 1
error number -128
end try
return {lit, len}
end ismatrix
on Msum(M1, M2, y) -- if y is a negative number then this handler
subtracts M2 from M1. In any other case, it adds them.
set mn to ismatrix(M1)
set nm to ismatrix(M2)
if mn ? nm then
display dialog "Can not add matrices of different type." buttons
{"Cancel"} default button 1
error number -128
end if
set Matriz to {}
repeat with m from 1 to length of M1
set coluna to {}
if (item 1 of mn) = 1 then
if y < 0 then set (item m of M2) to -(item m of M2)
set end of coluna to (item m of M1) + (item m of M2)
else
repeat with k from 1 to (item 1 of mn)
if (y as text) contains "-" then set (item k of item m of M2) to -(item
k of item m of M2) -- error here, have changed to the appropriate line
set end of coluna to ((item k of item m of M1) + (item k of item m of M2))
end repeat
end if
set end of Matriz to coluna
end repeat
return Matriz
end Msum
on Mprod(M1, M2)
set mn to ismatrix(M1)
set nm to ismatrix(M2)
if (item 2 of mn) ? (item 1 of nm) then
display dialog "You can not multiply these two matrices." buttons {"OK"}
default button 1
error number -128
end if
set Matriz to {}
repeat with m from 1 to item 2 of nm
set coluna to {}
repeat with k from 1 to item 1 of mn
set ckm to 0
repeat with i from 1 to item 2 of mn
--set aki to item k of item i of a
--set bim to item i of item m of b
if item 1 of mn = 1 and item 1 of nm ? 1 then set ckm to ckm + ((item i
of M1) * (item i of item m of M2))
if item 1 of nm = 1 and item 1 of mn ? 1 then set ckm to ckm + ((item i
of item k of M1) * (item i of M2))
if item 1 of nm = 1 and item 1 of mn = 1 then set ckm to ckm + ((item i
of M1) * (item i of M2))
if item 1 of nm ? 1 and item 1 of mn ? 1 then set ckm to ckm + ((item k
of item i of M1) * (item i of item m of M2))
end repeat
set end of coluna to ckm
end repeat
set end of Matriz to coluna
end repeat
return Matriz
end Mprod
on cronecker(k) -- this handler writes the identity matrix in k
dimensions in the list form. Items of the list are columns.
if (k ? 0) then return {}
try
set k to k as integer
on error
display dialog "k must be an integer." buttons {"OK"} default button 1
error number -128
end try
set Matriz to {}
repeat with l from 1 to k
set coluna to {}
repeat with i from 1 to k
if i = l then
set end of coluna to 1
else
set end of coluna to 0
end if
end repeat
set end of Matriz to coluna
end repeat
return Matriz
end cronecker
on mbyscalar(Matriz, Escalar)
set mn to ismatrix(Matriz)
repeat with m from 1 to length of Matriz
if (item 1 of mn) = 1 then
set (item m of Matriz) to Escalar * (item m of Matriz)
else
repeat with k from 1 to (item 1 of mn)
set (item k of item m of Matriz) to Escalar * (item k of item m of Matriz)
end repeat
end if
end repeat
return Matriz
end mbyscalar
on upperdiag(Matriz) -- this handler finds the upperdiagonal matrix U
similar to the square matrix Matriz.
set mn to ismatrix(Matriz)
if item 1 of mn ? item 2 of mn then
display dialog "This is not a square matrix." buttons {"OK"} default
button 1
error number -128
end if
if item 1 of mn = 1 then return Matriz
set Nmatriz to transpose(Matriz) -- novamatriz has lines as items and
matiz has columns as items.
repeat with j from 1 to (length of Matriz) - 1 -- both matrices have the
same length. they are square matrices.
set fator to 0 -- know when to leave
repeat with k from j + 1 to length of Matriz
if (item j of item j of Nmatriz) ? 0 then
set fator to ((item j of item k of Nmatriz) / (item j of item j of Nmatriz))
else
repeat with m from j + 1 to length of Matriz
if (item j of item m of Nmatriz) ? 0 then
log "1 " & (item j of item m of Nmatriz)
set {item j of Nmatriz, item m of Nmatriz} to {item m of Nmatriz, item j
of Nmatriz}
log "2 " & (item j of item m of Nmatriz)
set fator to ((item j of item k of Nmatriz) / (item j of item j of Nmatriz))
exit repeat
end if
end repeat
if item j of item j of Nmatriz = 0 then return Nmatriz
end if
repeat with l from 1 to length of Matriz
set (item l of item k of Nmatriz) to (item l of item k of Nmatriz) -
fator * (item l of item j of Nmatriz)
end repeat
end repeat
end repeat
return transpose(Nmatriz)
end upperdiag
on determinant(Matriz)
set Nmatriz to upperdiag(Matriz)
set det to 1
repeat with j from 1 to length of Nmatriz
set det to det * (item j of item j of Nmatriz)
end repeat
return det
end determinant
on inverse(Matriz)
set mn to ismatrix(Matriz)
if item 1 of mn ? item 2 of mn then
display dialog "This is not a square matrix." buttons {"OK"} default
button 1
error number -128
end if
if length of Matriz = 1 and item 1 of Matriz ? 0 then return {1 / (item
1 of Matriz)}
if length of Matriz = 1 and item 1 of Matriz = 0 then
display dialog "This matrix is not invertible." buttons {"OK"} default
button 1
error number -128
end if
set Nmatriz to transpose(Matriz) -- novamatriz has lines as items and
matiz has columns as items.
set IDT to cronecker(length of Matriz)
repeat with j from 1 to (length of Matriz) - 1 -- both matrices have the
same length. they are square matrices.
set fator to 0 -- know when to leave -- initialize
repeat with k from j + 1 to length of Matriz
if (item j of item j of Nmatriz) ? 0 then
set fator to ((item j of item k of Nmatriz) / (item j of item j of Nmatriz))
else
repeat with m from j + 1 to length of Matriz
if (item j of item m of Nmatriz) ? 0 then
set {item j of Nmatriz, item m of Nmatriz} to {item m of Nmatriz, item j
of Nmatriz}
set {item j of IDT, item m of IDT} to {item m of IDT, item j of IDT}
set fator to ((item j of item k of Nmatriz) / (item j of item j of Nmatriz))
exit repeat
end if
end repeat
if item j of item j of Nmatriz = 0 then
display dialog "This matrix is not invertible." buttons {"OK"} default
button 1
error number -128
end if
end if
repeat with l from 1 to length of Matriz
set (item l of item k of Nmatriz) to (item l of item k of Nmatriz) -
fator * (item l of item j of Nmatriz)
set (item l of item k of IDT) to (item l of item k of IDT) - fator *
(item l of item j of IDT)
end repeat
end repeat
end repeat
if item (item 1 of mn) of item (item 1 of mn) of Nmatriz = 0 then
display dialog "This matrix is not invertible." buttons {"OK"} default
button 1
error number -128
end if
repeat with j from (length of Matriz) to 1 by -1 -- both matrices have
the same length. they are square matrices.
set fator to 0 -- know when to leave
repeat with n from (length of Matriz) to 1 by -1
set (item n of item j of IDT) to ((item n of item j of IDT) / (item j of
item j of Nmatriz))
end repeat
set (item j of item j of Nmatriz) to 1
if j ? 1 then
repeat with k from j - 1 to 1 by -1
set fator to (item j of item k of Nmatriz)
repeat with l from length of Matriz to 1 by -1
set (item l of item k of Nmatriz) to (item l of item k of Nmatriz) -
fator * (item l of item j of Nmatriz)
set (item l of item k of IDT) to (item l of item k of IDT) - fator *
(item l of item j of IDT)
end repeat
end repeat
end if
end repeat
return transpose(IDT)
end inverse
on transpose(Matriz)
set mn to ismatrix(Matriz)
set Nmatriz to {}
repeat with k from 1 to item 1 of mn
set coluna to {}
repeat with i from 1 to item 2 of mn
if item 1 of mn = 1 then
set coluna to coluna & (item i of Matriz) --aik
else
set coluna to coluna & (item k of item i of Matriz)
end if
end repeat
set end of Nmatriz to coluna
end repeat
return Nmatriz
end transpose
on solvesystem(Matriz, coluna) -- solves the linear system Matriz*x=coluna
set z to inverse(Matriz)
return Mprod(z, coluna)
end solvesystem
on lowerdiag(Matriz) -- this handler finds the upperdiagonal matrix U
similar to the square matrix Matriz.
set mn to ismatrix(Matriz)
if item 1 of mn ? item 2 of mn then
display dialog "This is not a square matrix." buttons {"OK"} default
button 1
error number -128
end if
if item 1 of mn = 1 then return Matriz
set Nmatriz to transpose(Matriz) -- novamatriz has lines as items and
matiz has columns as items.
repeat with j from (length of Matriz) to 2 by -1 -- both matrices have
the same length. they are square matrices.
set fator to 0 -- know when to leave
repeat with k from j - 1 to 1 by -1
if (item j of item j of Nmatriz) ? 0 then
set fator to ((item j of item k of Nmatriz) / (item j of item j of Nmatriz))
else
repeat with m from j - 1 to 1
if (item j of item m of Nmatriz) ? 0 then
log "1 " & (item j of item m of Nmatriz)
set {item j of Nmatriz, item m of Nmatriz} to {item m of Nmatriz, item j
of Nmatriz}
log "2 " & (item j of item m of Nmatriz)
set fator to ((item j of item k of Nmatriz) / (item j of item j of Nmatriz))
exit repeat
end if
end repeat
if item j of item j of Nmatriz = 0 then return Nmatriz
end if
repeat with l from 1 to length of Matriz
set (item l of item k of Nmatriz) to (item l of item k of Nmatriz) -
fator * (item l of item j of Nmatriz)
end repeat
end repeat
end repeat
return transpose(Nmatriz)
end lowerdiag
_______________________________________________
applescript-users mailing list | email@hidden
Help/Unsubscribe/Archives:
http://www.lists.apple.com/mailman/listinfo/applescript-users
Do not post admin requests to the list. They will be ignored.