Numpyの行列計算eig()とeigh()の出力の違い
numpyの行列計算でよくわからない仕様があったのでメモ。
どうしてこの問題に行き着いたかはこちらの記事参照
http://buhin-blog.blogspot.jp/2018/03/blog-post_18.html
行列Aの固有値計算には一般行列用のeig(A)とエルミート行列用のeigh(A)がある。
エルミート行列Aをeigh(A)で計算した場合とeig(A)で計算した場合、固有値リスト内の順番が違うようだ。
おそらく仕様なのだけれど、よくわからない（つらい）。
順番をSortすると同じになるので、出てくる固有値は一致しているみたい。
In [16]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
#1.1 束縛状態を持つ1次元シュレーディンガー方程式
#1.1.1 ポテンシャルを一つ持つ場合の束縛状態
V0=-5.0
def V(x): #ポテンシャル
V=V0*np.exp(-x**2)
return V
N=100
xmax=10.0
xmin=-10.0
xvec=np.linspace(xmin, xmax, N) #横軸を示すリスト
#ポテンシャルのプロット
plt.plot(xvec, V(xvec), label="Potential")
Out[16]:
In [17]:
def make_H1d(N, xmax, xmin, V): # ポテンシャルVのときのNxN行列化
mat_H=np.zeros([N,N]) #空の行列を作る
#V(x)=-1*V0*np.exp(-1*x*x/(ξ*ξ))
a = (xmax-xmin)/(N-1)
for i in range(1,N+1):
x= (i-1)*a + xmin
for dx in range(-1, 2):
j=i+dx
v=0.0
if dx==0:
v=(2/a**2 + V(x))
elif dx==1:
v=-1/a**2
elif dx== -1 :
v=-1/a**2
if 1 <= j <= N:
mat_H[i-1,j-1]=v #ハミルトニアンの[i-1,j-1]成分
return mat_H
In [18]:
#N=100の場合は途中から振動してうわ～という感じ
N=100
mat_H = make_H1d(N,xmax,xmin,V)
energy,mat_v = np.linalg.eigh(mat_H) #エルミート行列の固有値と固有関数、ここが問題
energy2,mat_v2 = np.linalg.eig(mat_H)
xvec2=np.linspace(0, 100, N) #横軸を示すリスト
plt.plot(xvec2,energy, c="blue", label="eigh")
plt.plot(xvec2,energy2, c="red", label="eig")
plt.legend()
Out[18]:
In [19]:
#N=10の場合でもズレが有る。
N=10
mat_H = make_H1d(N,xmax,xmin,V)
energy,mat_v = np.linalg.eigh(mat_H) #エルミート行列の固有値と固有関数、ここが問題
energy2,mat_v2 = np.linalg.eig(mat_H)
xvec2=np.linspace(0, 100, N) #横軸を示すリスト
plt.plot(xvec2,energy, c="blue", label="eigh")
plt.plot(xvec2,energy2, c="red", label="eig")
plt.legend()
Out[19]:
In [20]:
#N=1000の場合はもはや固有値の順番が逆順になってる
N=1000
mat_H = make_H1d(N,xmax,xmin,V)
energy,mat_v = np.linalg.eigh(mat_H) #エルミート行列の固有値と固有関数、ここが問題
energy2,mat_v2 = np.linalg.eig(mat_H)
xvec2=np.linspace(0, 100, N) #横軸を示すリスト
plt.plot(xvec2,energy, c="blue", label="eigh")
plt.plot(xvec2,energy2, c="red", label="eig")
plt.legend()
Out[20]:
In [21]:
#N=1000の場合でsort()した場合は一致する(きれい)
N=1000
mat_H = make_H1d(N,xmax,xmin,V)
energy,mat_v = np.linalg.eigh(mat_H) #エルミート行列の固有値と固有関数、ここが問題
energy2,mat_v2 = np.linalg.eig(mat_H)
xvec2=np.linspace(0, 100, N) #横軸を示すリスト
energy.sort()
energy2.sort()
plt.plot(xvec2,energy, c="blue", label="eigh")
plt.plot(xvec2,energy2, c="red", label="eig")
plt.legend()
Out[21]:
コメント
コメントを投稿