math - How to change elements in sparse matrix in Python's SciPy? -
i have built small code want use solving eigenvalue problems involving large sparse matrices. it's working fine, want set elements in sparse matrix zero, i.e. ones in top row (which corresponds implementing boundary conditions). can adjust column vectors (c0, c1, , c2) below achieve that. however, wondered if there more direct way. evidently, numpy indexing not work scipy's sparse package.
import scipy.sparse sp import scipy.sparse.linalg la import numpy np import matplotlib.pyplot plt #discretize x-axis n = 11 x = np.linspace(-5,5,n) print(x) v = x * x / 2 h = len(x)/(n) hi2 = 1./(h**2) #discretize schroedinger equation, i.e. build #banded matrix difference equation c0 = np.ones(n)*30. + v c1 = np.ones(n) * -16. c2 = np.ones(n) * 1. diagonals = np.array([-2,-1,0,1,2]) h = sp.spdiags([c2, c1, c0,c1,c2],[-2,-1,0,1,2], n, n) h *= hi2 * (- 1./12.) * (- 1. / 2.) #solve eigenvalues ev = la.eigsh(h,return_eigenvectors = false) #check structure of h plt.figure() plt.spy(h) plt.show()
this visualisation of matrix build code above. want set elements in first row zero.
as suggested in comments, i'll post answer found own question. there several matrix classes in in scipy's sparse package, listed here. 1 can convert sparse matrices 1 class another. need do, choose convert sparse matrix class csr_matrix, by
h = sp.csr_matrix(h)
then can set elements in first row 0 using regular numpy notation:
h[0,0] = 0 h[0,1] = 0 h[0,2] = 0
for completeness, post full modified code snippet below.
#scipy sparse linear algebra takes care of sparse matrix computations #http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html import scipy.sparse sp import scipy.sparse.linalg la import numpy np import matplotlib.pyplot plt #discretize x-axis n = 1100 x = np.linspace(-100,100,n) v = x * x / 2. h = len(x)/(n) hi2 = 1./(h**2) #discretize schroedinger equation, i.e. build #banded matrix difference equation c0 = np.ones(n)*30. + v c1 = np.ones(n) * -16. c2 = np.ones(n) * 1. h = sp.spdiags([c2, c1, c0, c1, c2],[-2,-1,0,1,2], n, n) h *= hi2 * (- 1./12.) * (- 1. / 2.) h = sp.csr_matrix(h) h[0,0] = 0 h[0,1] = 0 h[0,2] = 0 #check structure of h plt.figure() plt.spy(h) plt.show() ev = la.eigsh(h,return_eigenvectors = false)
Comments
Post a Comment