2013-04-09 10 views
12

大きなスパース行列を含む固有値問題を解くために使用したい小さなコードを作成しました。うまくいきました。今私がしたいのは、スパース行列のいくつかの要素をゼロに設定することです(つまり、境界条件を実装することに相当します)。私はちょうどそれを達成するために、列ベクトル(C0、C1、およびC2)を調整することができます。しかし、より直接的な方法があるかどうかは疑問でした。明らかに、NumPyインデックスはSciPyの疎なパッケージでは機能しません。PythonのSciPyでスパース行列の要素を変更するには?

import scipy.sparse as sp 
import scipy.sparse.linalg as la 
import numpy as np 
import matplotlib.pyplot as 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 from 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 for eigenvalues 
EV = la.eigsh(H,return_eigenvectors = False) 

#check structure of H 
plt.figure() 
plt.spy(H) 
plt.show() 

これは、上記のコードで構築されたマトリックスの視覚化です。最初の行の要素をゼロにしたい。 enter image description here

+0

私は回避策を見つけました。私が使用しているフォーマット(dia_matrix)は、私が達成したいと思っているのは良くありません。代わりにcsr_matrixを使用します。私はこの投稿を閉じなければなりませんか? – seb

+0

これはよく書かれた質問であり、将来他の人に役立つかもしれません。答えとして見つけたものを投稿するのはどうですか? – YXD

+0

OK、私はそれを実行します – seb

答えて

13

コメントに示唆されているように、私は自分の質問に答えを投稿します。 SciPyの疎なパッケージにはいくつかの行列クラスがあり、それらはhereと記載されています。 1つのクラスから別のクラスに疎行列を変換することができます。だから私は何をする必要があるかのために、私は単に

H = sp.csr_matrix(H) 

により、クラスcsr_matrixに私のスパース行列を変換することを選択した後、私は定期的にnumpyの表記を使用して0に最初の行の要素を設定することができます。

H[0,0] = 0 
H[0,1] = 0 
H[0,2] = 0 

完全性のために、以下の完全なコードスニペットを投稿します。

#SciPy Sparse linear algebra takes care of sparse matrix computations 
#http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html 
import scipy.sparse as sp 
import scipy.sparse.linalg as la 

import numpy as np 
import matplotlib.pyplot as 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 from 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) 
+1

列より多くの行がある場合、csrは高速ですが、行より多くの列がある場合、cscは高速です。 –

関連する問題