5.1.1. 位置表示の行列表現

(Fourier grid Hamiltonian Method というらしい)

位置表示の基底 \(|q\rangle\) は完全系

\[\int dq |q\rangle \!\langle q | = \mbox{1}\hspace{-0.25em}\mbox{l}\]

をなすので, ハミルトニアン \(\hat{H} = T(\hat{p}) + V(\hat{q})\) のq-表示の行列表現を求めると

\[\langle q_i | \hat{H} | q_j \rangle = \langle q_i | T(\hat{p}) | q_j\rangle + \langle q_i | V(\hat{q}) | q_j \rangle\]

であり,第二項は

\[\langle q_i | V(\hat{q}) | q_j \rangle = V(q_j)\delta_{i,j}\]

である.第一項目は

\[\begin{split}\langle q_i | T(\hat{p}) | q_j\rangle & = \sum_{p_k} \langle q_i | T(\hat{p}) | p_k\rangle\!\langle p_k| q_j\rangle & = \frac{1}{2\pi\hbar}\sum_{p_k} e^{\sqrt{-1}p_k(q_i-q_j)/\hbar} T(p_k)\end{split}\]

でとしてFourier変換を通じて与えられる. こうして、ハミルトニアンのq-表示の行列表現が与えられるので,数値的に対格化すれば良い.

かなりいいかげんに作ったが,sample コードは以下の通りである.

# fourier_grid.py

import numpy as np

twopi = 2.0*np.pi
class Harmonic(object):
    def V(self, x):
        return 0.5*x*x
    
    def T(self, x):
        return 0.5*x*x

class FGHM(object):
    def __init__(self, dim, domain, H):
        self.dim = dim
        self.domain = domain
        self.__set_domain(domain)
        self.H = H
        
    def __set_domain(self,domain):
        qmin,qmax = domain[0][0],domain[0][1]
        pmin,pmax = domain[1][0],domain[1][1]
        if (qmin > qmax) or (pmin > pmax):
            raise ValueError("Domain Error p in [%f,%f], p in [%f %f]\n" % (qmin,qmax, pmin,pmax))
        if pmax != -pmin:
            raise ValueError("p in (min,max) and max=-min\n")
        q = np.linspace(qmin, qmax,self.dim,endpoint=False)
        p = np.linspace(pmin, pmax,self.dim,endpoint=False)
        self.x = [q,p]
        self.h = (qmax - qmin)*(pmax - pmin)/self.dim

    def V(self):
        eye = np.eye(self.dim, dtype=np.complex128)
        return self.H.V(self.x[0])*eye

    def _T(self):
        mat = np.zeros([self.dim,self.dim],dtype=np.complex128)
        q,p = self.x[0],self.x[1]
        hbar = self.h/twopi
        for i in range(self.dim): # loop of q_i
            for j in range(self.dim): # loop of q_j
                sum = np.sum(np.exp(1.j*p*(q[i] -q[j])/hbar)*self.H.T(p))/self.dim #loop of p
                mat[i][j] = sum
        return mat
    
    def T(self):
        eye = np.eye(self.dim)
        mat = np.zeros([self.dim,self.dim],dtype=np.complex128)
        q,p = self.x[0],self.x[1]
        for i,e in enumerate(eye):
            pbase = np.fft.fft(e)
            qbase = np.fft.ifft(pbase*self.H.T(p))
            mat[i] = qbase
        return mat.transpose()
    
    def getMatrix(self, fft=True):
        if fft:
            Ham = self.T() + self.V()
        else:
            Ham = self._T() + self.V()
        return Ham 

    def eigen(self,fft=True):
        evals, evecs = np.linalg.eig(self.getMatrix(fft))
        sort_index = [i[0] for i in sorted(enumerate(evals.real), key=lambda x:x[1])]
        evecs = evecs.transpose()[sort_index]
        return evals[sort_index], [self.normalize(vec) for vec in evecs]
        
    def normalize(self, x):
        norm = np.abs(np.sum(x*np.conj(x)))
        return x/np.sqrt(norm)
        
    

上記ファイルを作って,別のファイルから読み込んで実行してみよう.例えば、調和振動子

\[\hat{H} = \frac{1}{2}\hat{p}^2 + \frac{1}{2}\hat{q}^2\]

の固有エネルギーは

\[E_n = \hbar(n + \frac{1}{2})\]

で与えられる.ので \(E_n/\hbar\) vs \(n\) の次元依存性を見てみよう.

# eiven_vals.py

from fourier_grid import *
import matplotlib.pyplot as plt
import numpy as np

twopi = 2.0*np.pi
dim_list = [10, 50, 100,200]
qmin, qmax = -5.0,5.0
pmin, pmax = -5.0,5.0
domain = [[qmin, qmax],[pmin,pmax]]
for dim in dim_list:
    Ham = FGHM(dim, domain, Harmonic())
    hbar = Ham.h/twopi
    evals, evecs = Ham.eigen(fft=True)
    x = np.linspace(0.0,dim, dim)
    plt.plot(x,x,'--')
    plt.plot(np.arange(dim),evals.real/hbar, '-o',label="dim=%d" % dim)
    print("dim:%d" %dim, "E_n/hbar",evals.real[0:5]/hbar)
plt.xlabel(r"$n$ (quantum number)")
plt.ylabel(r"$E_n/\hbar$")
plt.legend(loc=2,numpoints=1)
plt.savefig("energy_vs_n.png")
plt.show()

実行結果は

dim:10 E_n/hbar [ 0.4999988   1.5000347   2.499525    3.50382845  4.47835611]
dim:50 E_n/hbar [ 0.5  1.5  2.5  3.5  4.5]
dim:100 E_n/hbar [ 0.5  1.5  2.5  3.5  4.5]
dim:200 E_n/hbar [ 0.5  1.5  2.5  3.5  4.5]
_images/energy_vs_n.png

また,固有関数は

# eiven_vecs.py

from fourier_grid import *
import matplotlib.pyplot as plt
import numpy as np

twopi = 2.0*np.pi
dim = 100
qmin, qmax = -5.0,5.0
pmin, pmax = -5.0,5.0
domain = [[qmin, qmax],[pmin,pmax]]
Ham = FGHM(dim, domain, Harmonic())
hbar = Ham.h/twopi
evals, evecs = Ham.eigen(fft=True)
x = np.linspace(0.0,dim, dim)
for i,vec in enumerate(evecs):
    plt.plot(Ham.x[0],Ham.H.V(Ham.x[0]), '-k', linewidth=2)
    plt.plot(Ham.x[0], np.abs(vec*np.conj(vec)) + evals[i].real, '-', linewidth=3)
plt.ylabel(r"$V(q)$,$|\psi_n(q)|^2$",fontsize=20)
plt.xlabel(r"$q$",fontsize=20)
plt.ylim(0,2)
plt.xlim(-3,3)
plt.savefig("eigen_state_with_n.png")
plt.show()
_images/eigen_state_with_n.png

課題

nが3/4程度を超えると厳密解からずれ始める. 固有関数は位相空間上で円形で与えられる関わらず, 位相空間を長方形でgridを切って近似している為であると思われるので考察せよ.