5.1.3. \(q^{4}\) ポテンシャル

[1] のサンプルにある \(V(\hat{q})=\hat{q}^4/2\) を解いてみる.

せっかくなので (i) 位置表示の行列表現 (ii) 調和振動子の基底を用いた行列表現 を比較してみよう.

#quad_potential_evals.py

class Quad(object):
    def T(self, x):
        return 0.5*x*x
    def V(self, x):
        return 0.5*x*x*x*x
        
from fourier_grid import FGHM
from harmonic_base import Operator
import numpy as np
import matplotlib.pylab as plt

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]]

Ham = Quad()
fig = plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(2,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,3)
ax4 = fig.add_subplot(2,2,4)
for i,dim in enumerate(dim_list):
    OP = Operator(dim,domain, Ham)
    FG = FGHM(dim, domain, Ham)
    evals0 = OP.eigen(sort=True, vector=False)
    evals1, evecs1 = FG.eigen(fft=True)
    card = 1.0*i/len(dim_list)
    x = np.arange(dim,dtype=np.float64)
    ax1.plot(x, evals0.real,'-',label='dim=%d'% dim,color=plt.cm.jet(card),linewidth=4)
    ax1.plot(x, evals1.real,'o',color=plt.cm.jet(card))
    
    ax3.plot(x, evals0.real,'-',label='dim=%d'% dim,color=plt.cm.jet(card),linewidth=4)
    ax3.plot(x, evals1.real,'o',color=plt.cm.jet(card))
    
    ax2.plot(x/dim, evals0.real,'-',label='dim=%d'%dim,color=plt.cm.jet(card),linewidth=4)
    ax2.plot(x/dim, evals1.real,'o',color=plt.cm.jet(card))
    
    ax4.plot(x/dim, evals0.real,'-',label='dim=%d'%dim,color=plt.cm.jet(card),linewidth=4)
    ax4.plot(x/dim, evals1.real,'o',color=plt.cm.jet(card))
    
ax1.set_xlabel(r"$n$",fontsize=20)
ax1.set_ylabel(r"$E_n$",fontsize=20)
ax2.set_xlabel(r"$n/\mathrm{dim}$",fontsize=20)
ax2.set_ylabel(r"$E_n$",fontsize=20)
ax3.set_xlabel(r"$n$",fontsize=20)
ax3.set_ylabel(r"$E_n$",fontsize=20)
ax4.set_xlabel(r"$n/\mathrm{dim}$",fontsize=20)
ax4.set_ylabel(r"$E_n$",fontsize=20)

ax3.loglog()
ax4.loglog()
ax2.legend(loc=0,numpoints=1)
ax4.legend(loc=0,numpoints=1)
plt.tight_layout()
plt.savefig("eigen_vals_quad.png")
plt.show()
_images/eigen_vals_quad.png

曲線が(i)の方法で,ドットが(ii)の方法で解いた結果. 色の違いは行列の次元dimの違い. 高エネルギー準位における両者の不一致が見られるが,低エネルギー準位では両者の結果は良く一致している.

同様に固有状態についても確認する.

#quad_potential_evecs.py

class Quad(object):
    def T(self, x):
        return 0.5*x*x
    def V(self, x):
        return 0.5*x*x*x*x
        
from fourier_grid import FGHM
from harmonic_base import Operator
import numpy as np
import matplotlib.pylab as plt

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 = Quad()
OP = Operator(dim,domain, Ham)
FG = FGHM(dim, domain, Ham)
evals0, evecs0= OP.eigen(sort=True, vector=True)
evals1, evecs1 = FG.eigen(fft=True)    
fig = plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(2,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,3)
ax4 = fig.add_subplot(2,2,4)
x = OP.x[0]
plot_num=10
for i in range(dim):
    if i==plot_num:break
    card = 1.0*i/plot_num
    ax1.plot(x,np.abs(evecs0[i]*np.conj(evecs0[i]))+evals0[i].real, '-',color=plt.cm.jet(card),linewidth=3)
    ax2.plot(x,np.abs(evecs1[i]*np.conj(evecs1[i]))+evals1[i].real, '-',color=plt.cm.jet(card),linewidth=3)
    if i in [0,5,9]:
        ax3.plot(x,np.abs(evecs0[i]*np.conj(evecs0[i])), '-',color=plt.cm.jet(card),linewidth=3)
        ax4.plot(x,np.abs(evecs1[i]*np.conj(evecs1[i])), '-',color=plt.cm.jet(card),linewidth=3)    

ymin,ymax = 0, 2
ax1.plot(x,Ham.V(x),'-k')
ax2.plot(x,Ham.V(x),'-k')
ax1.set_xlabel(r"$q$",fontsize=20)
ax1.set_ylabel(r"$|\Psi_n(q)|^2$",fontsize=20)
ax1.set_ylim(ymin,ymax)
ax1.set_xlim(qmin,qmax)

ax2.set_xlabel(r"$q$",fontsize=20)
ax2.set_ylim(ymin,ymax)
ax2.set_xlim(qmin,qmax)

ax3.set_xlabel(r"$q$",fontsize=20)
ax3.set_ylabel(r"$|\Psi_n(q)|^2$",fontsize=20)
ax3.set_xlim(qmin,qmax)
ax3.semilogy()


ax4.set_xlabel(r"$q$",fontsize=20)
ax4.set_xlim(qmin,qmax)
ax4.semilogy()


ax2.legend(loc=0,numpoints=1)
ax4.legend(loc=0,numpoints=1)
#plt.tight_layout()
plt.savefig("eigen_state_quad.png")
plt.show()
_images/eigen_state_quad.png

エネルギーの低い状態から順に9番目までの固有状態.左図が(i)の方法で,右図は(ii)の方法. 上図はnormal plotで,下図はsemilog plotである(n=0,5,9).

課題

波動関数の形状はほぼ一致しているが,ポテンシャルの外に置けるtunnel tailが両者で異なる. (i)には不自然な振動が生じているが,これは \(\infty<x<\infty\) で定義される基底を使って作っているからであろうか? それとも打ち切っているからであろうか? ちなみに,y軸の \(10^{-32}\) 程度がdoubleの精度限界. 固有関数を求めるのであれば,(ii)の方法の方が良いかもしれない.

課題

厳密解が分るモデルか,摂動論による計算結果の確認.

[1]H. J. Korsch, H. J. and M. Glück, Euro. J. Phys, 23 (2002) 413