[SciPy-user] generalized eigenvalue problem for sparse matrices
abhinav sarkar
abhinav.sarkar@gmail....
Mon Apr 14 15:20:20 CDT 2008
Hi
I am trying to solve a generalized eigenvalue problem for sparse
matrices. The problem is in form
A*x = 𝜎*M*x
where A and M are sparse matrices and x is a vector and sigma is a
scalar whose value is to be found.
For this I am using the Arpack function ARPACK_gen_eigs provided in
the module scipy.sparse.linalg.eigen.arpack.speigs . However the
solution which I am getting does not match with the solution obtained
from the eig function in MATLAB. For example:
A = [1 0 0 -33 16 0; 0 1 0 16 -33 16; 0 0 1 0 16 -33; 1601 -96 256
-1800 0 0; -96 1601 -96 0 -1800 0; 256 -96 1601 0 0 -1800]
M = [0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; -33 16 0 0 0 0; 16 -33 16
0 0 0; 0 16 -33 0 0 0]
The matrices are written in MATLAB format, spaces separating the row
elements and ";" separating the rows.
When I solve this problem in MATLAB using the function eig, I get the
following eigenvalues:
-155.0345, -56.9898, -45.2209, -31.9376, -28.5367, -9.1611
However when I solve it using ARPACK_gen_eigs to find the two
eigenvalues near 10, I get the following solution:
13.83675153-1.71439075j, 13.83675153+1.71439075j
which is certainly not correct.
I am using the following code to solve the problem:
from scipy import sparse
from scipy.sparse.linalg import dsolve
from scipy.sparse.linalg import eigen
from scipy.sparse.linalg.eigen.arpack.speigs import ARPACK_gen_eigs
n = 3
h = 1.0/(n+1)
a = 1
Pr = 7.0
Ra = 1000
A = get_A_mat(n, h, a, Pr, Ra)
M = get_M_mat(n, h, a, Pr, Ra)
sigma = 10.0
B = A - sigma*M
s = dsolve.splu(B)
e, v = ARPACK_gen_eigs(M.matvec, s.solve, 2*n, sigma, 2, 'LR')
print e
which also seems to be correct to me. Please tell if the method I am
using is correct or not and why am I not getting correct solutions. Is
the ARPACK_gen_eigs broken? Or is there a problem in my code?
Regards
--
Abhinav Sarkar
4th year Undergraduate Student
Deptt. of Mechanical Engg.
Indian Institute of Technology, Kharagpur
India
Web: http://claimid.com/abhin4v
Twitter: http://twitter.com/abhin4v
---------
The world is a book, those who do not travel read only one page.
More information about the SciPy-user
mailing list