Numpyの行列計算eig()とeigh()の出力の違い

eig比較

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]:
[<matplotlib.lines.Line2D at 0x2b7036b1908>]
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]:
<matplotlib.legend.Legend at 0x2b70371db00>
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]:
<matplotlib.legend.Legend at 0x2b703784e48>
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]:
<matplotlib.legend.Legend at 0x2b7037f3b70>
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]:
<matplotlib.legend.Legend at 0x2b703865400>

コメント

このブログの人気の投稿

学振採用者はどこへ消えた?

物理系研究関係者、ツイッターやりすぎランキング(ぶひん調べ)

オレ達はあと何本論文を書けば東大教授になれるんだ?