[SciPy-user] generalized eigenvalue problem for sparse matrices

[abhinav sarkar]

On Tue, Apr 15, 2008 at 6:24 PM, Neilen Marais <nmarais at sun.ac.za> wrote:
> Abinav,
>
>  On Tue, 15 Apr 2008 01:50:20 +0530, abhinav sarkar wrote:
>
>  > Hi
>
> >  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
>  ^^^
>  I think you're looking for -10.
i am looking for the largest eigenvalue
>
>
>  >  B = A - sigma*M
>  >  s = dsolve.splu(B)
>  >  e, v = ARPACK_gen_eigs(M.matvec, s.solve, 2*n, sigma, 2, 'LR')
>                                                              ^^
>  You may try futzing around with this parameter.
>
>  Anyway, I did play around a bit with your example and never really got a
>  good answer. But my experience has been that ARPACK doesn't really work
>  well with small matrices. Try a bigger problem and see if things turn out
>  better, I mean, obviously you don't need a sparse solver to solve a 6x6
>  matrix system ;). You may also find that if you force Matlab to do this
>  problem using a sparse solver you'll run into the same issues.
>  Unfortunately I don't know enough about the ARPACK routines to give more
>  insight. They do work quite well on my problems, where the eigenvalues
>  are always positive.
>
>
>  >  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?
>
>  Seems OK, but I'm no expert ;)
>
>  Regards
>  Neilen
>
>
>  >
>  >  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.
>  > _______________________________________________ SciPy-user mailing list
>
>
> > SciPy-user at scipy.org
>  > http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
>  _______________________________________________
>  SciPy-user mailing list
>  SciPy-user at scipy.org
>  http://projects.scipy.org/mailman/listinfo/scipy-user
>

 I have formulated the problem in a slightly different way and now I
am able to obtain the correct eigenvalues. Here is the new code:

class Solver:
	def __init__(self, A, M, sigma):
		self.A = A
		self.M = M
		self.sigma = sigma
		
	def solve(self, b):
		#(M-sigma*A)*x=A*b
		#return x
		return dsolve.spsolve((self.M-self.sigma*self.A), self.A*b)

A = get_A_mat(n, h, a, Pr, Ra)
M = get_M_mat(n, h, a, Pr, Ra)
I = sparse.eye(2*n,2*n)
sigma = 0.0
sigma_solve = Solver(A, M, sigma).solve
e = 1.0/ARPACK_gen_eigs(I.matvec, sigma_solve, I.shape[0], sigma, 5, 'LR')[0]

Here I have converted the A*x = sigma*M*x to inv(A)*M*x = nu*I*x where
nu = 1/sigma. And replaced s.solve by my own Solver class. Now
solutions obtained are correct.
-- 
Abhinav Sarkar
4th year Undergraduate Student
Deptt. of Mechanical Engg.
Indian Institute of Technology, Kharagpur
India

Mobile:+91-99327-41517
Residence:+91-631-2429887
Web: http://claimid.com/abhin4v
Twitter: http://twitter.com/abhin4v
---------
The world is a book, those who do not travel read only one page.