《計算物理學(xué)》課件第13章

13.1定態(tài)薛定諤方程13.2一維方勢阱中粒子能級和波函數(shù)的計算機(jī)求解13.3薛定諤方程的矩陣解法習(xí)題十三量子力學(xué)與經(jīng)典力學(xué)的差別首先表現(xiàn)在對粒子的狀態(tài)和力學(xué)量的描述及其變化規(guī)律上。在量子力學(xué)中，微觀粒子的運(yùn)動狀態(tài)稱為量子態(tài)，是用波函數(shù)ψ(r，t)來描述的，波函數(shù)所反映的微觀粒子波動性，就是德布羅意波。德布羅意波或波函數(shù)ψ(r，t)不代表實(shí)際物理量的波動，而是描述粒子在空間的概率分布的概率波。13.1定態(tài)薛定諤方程量子力學(xué)中描述微觀粒子的波函數(shù)本身是沒有直接物理意義的，具有直接物理意義的是波函數(shù)模的平方，它代表了粒子出現(xiàn)的概率。波函數(shù)模的平方|ψ(r，t)|2代表時刻t，在r處附近空間單位體積中粒子出現(xiàn)的幾率。因此|ψ(r，t)|2也被稱為概率密度，即某一時刻出現(xiàn)在某點(diǎn)附近在體積元dV中的粒子的概率為|ψ(x，y，z，t,)|2dxdydz。波函數(shù)必須滿足標(biāo)準(zhǔn)化條件，即單值、連續(xù)、有限；同時波函數(shù)還必須滿足歸一化條件ψ*(r，t)ψ(r，t)dτ=1。當(dāng)微觀粒子處于某一狀態(tài)時，它的力學(xué)量(如坐標(biāo)、動量、角動量、能量等)具有一系列可能值，每個可能值以一定的幾率出現(xiàn)。當(dāng)粒子所處的狀態(tài)確定時，力學(xué)量具有某一可能值的幾率也就完全確定了。量子力學(xué)中求解粒子問題常歸結(jié)為解薛定諤方程或定態(tài)薛定諤方程。薛定諤方程廣泛地用于原子物理、核物理和固體物理，對于原子、分子、核、固體等一系列問題中求解的結(jié)果都與實(shí)際符合得很好。薛定諤方程通常表示為

(13.1)即

(13.2)

其中式(13.2)為偏微分方程，這里u(r，t)稱為勢函數(shù)，ψ(r，t)稱為波函數(shù)，m為粒子的質(zhì)量。我們的目的是求能級E和波函數(shù)ψ。若u=u(r)，不含時間，則有

ψ(r，t)=ψ(r)f(t)

(13.3)

代入含時薛定諤方程(13.2)后，可得

(13.4)因此有

(13.5)根據(jù)式(13.3)有

(13.6)

式(13.5)中的第二式稱為定態(tài)薛定諤方程，即

(13.7)在一維情況下有

(13.8)

通常又寫為算子形式

(13.9)

顯然該方程為一本征方程，E即為本征值，而ψ為本征函數(shù)。

對于一維方勢阱問題，勢函數(shù)u(x)=

如圖13.1所示，試根據(jù)式(13.9)求解本征方程的本征值E(能級)和本征函數(shù)ψ(波函數(shù))。13.2一維方勢阱中粒子能級和波函數(shù)的計算機(jī)求解圖13.1一維方勢阱示意圖

根據(jù)式(13.9)，有

即

為簡單起見，記(不影響結(jié)果的一般性)，因此有

(13.10)阱內(nèi)波函數(shù)滿足(0≤x≤W)

其中(V0<E<0)，上面微分方程的解即為

在阱外(x≤0，x≥W)，由于u(x)=0，因此有

其中r2=，上式微分方程的解為

這里j=0表示勢阱左邊的波函數(shù)，而j=1表示勢阱右邊的波函數(shù)。待定系數(shù)有A1、B1、C0、C1、D0、D1。

根據(jù)自然邊界條件，x→±∞時，波函數(shù)ψ(x)→0，因此對于勢阱左側(cè)可得D0=0，而在右側(cè)有C1=0。再利用邊界條件，在x=0處波函數(shù)和一階導(dǎo)數(shù)連續(xù)，有

ψ(0-)=C0=ψ(0+)=B1

ψ′(0-)=r2C0=ψ′(0+)=r1A1因此有

C0為常數(shù)，一般可取為1。因此波函數(shù)可以表示為

(13.11(a))

(13.11(b))

(13.11(c))

同理，在x=W處，波函數(shù)和一階導(dǎo)數(shù)連續(xù)，因此有

(13.12(a))

(13.12(b))

(13.12(c))

(13.12(d))由式(13.12(b))和(13.12(d))可得

(13.13)故由C0、ψ(W)、ψ′(W)可確定C1、D1。結(jié)合方程(13.12(a))~(13.12(d))，可得

(13.14(a))

(13.14(b))式(13.14)中的r1、r2

可以分別通過

來確定。要確定本征值E可以通過式(13.14(a))中的C1=0確定，這就要求對E取值有限制，必須取確定值。但這種方法一般不予采用，而通常采用計算節(jié)點(diǎn)法來實(shí)現(xiàn)。我們知道對勢阱而言，一般有Emax=0(否則勢阱不起作用，變?yōu)樽杂闪Ｗ?，Emin=V0。需要說明的是，與基態(tài)能量E0相應(yīng)的波函數(shù)沒有節(jié)點(diǎn)；第n個激發(fā)態(tài)相應(yīng)的波函數(shù)ψn有n個節(jié)點(diǎn)，而且這些能量值恰好是波函數(shù)的節(jié)點(diǎn)變化時臨界的能量值。我們可以通過計算不同能量E對應(yīng)解波函數(shù)的節(jié)點(diǎn)數(shù)。確定節(jié)點(diǎn)變化時臨界的能量值E即為能量本征值。

計算節(jié)點(diǎn)時，方勢阱內(nèi)節(jié)點(diǎn)數(shù)是用阱內(nèi)的半波長數(shù)r1W/π(一般是介于某奇數(shù)與相鄰偶數(shù)間的數(shù))取整數(shù)后決定的，即若ψ(0)·ψ(W)>0則取相鄰的偶數(shù)，若ψ(0)·、ψ(W)<0則取相鄰的奇數(shù)作為節(jié)點(diǎn)數(shù)；方勢阱外節(jié)點(diǎn)數(shù)是令ψ(x)=0，由式(13.13)解得的C1和D1異號，并且

，則有一節(jié)點(diǎn)，否則就無節(jié)點(diǎn)。計算節(jié)點(diǎn)和能級的步驟如下:

(1)輸入V0、W、Emax、Emin和M。

(2)計算

(3)利用以上節(jié)點(diǎn)計算方法確定節(jié)點(diǎn)數(shù)。

(4)由節(jié)點(diǎn)數(shù)計算結(jié)果來定出能級。①若相鄰的兩個能量分別對應(yīng)0、1兩個節(jié)點(diǎn)，則E1必處于Ea、Eb間，取

②若相鄰的兩個能量分別對應(yīng)一個和兩個節(jié)點(diǎn)，則E2必處于這兩個能量間，依此類推。

【例13.1】

例如取V0=-20.0，W=1.0，Emin=-20.0，Emax=0.0，M取51。利用計算節(jié)點(diǎn)法編程計算出基態(tài)和第一激發(fā)態(tài)能級。

【例13.1】

例如取V0=-20.0，W=1.0，Emin=-20.0，Emax=0.0，M取51。利用計算節(jié)點(diǎn)法編程計算出基態(tài)和第一激發(fā)態(tài)能級。

計算程序:

write(*,*)′inputV0,W,Emin,Emax,M′

read(*,*)V0,W,Emin,Emax,M

N=0

de=(Emax-Emin)/(M-1)

E=Emin-de

do100I=1,M

E=E+de

r1=sqrt(abs(E-V0))

r2=sqrt(abs(E))

phi=(r2/r1)*sin(W*r1)+cos(W*r1)!對應(yīng)式(13.12(a))

phi1=r2*cos(W*r1)-r1*sin(W*r1)!對應(yīng)式(13.12(c))

C1=0.5*exp(-W*r2)*(phi+phi1/r2)

D1=0.5*exp(W*r2)*(phi-phi1/r2)

N=int(W*r1/3.1415926)if(N-int(N/2)*2.ne.0.and.phi.gt.0)N=N+1![KG-*4]N為奇數(shù)

if(N-int(N/2)*2.eq.0.and.phi.lt.0)N=N+1![KG-*4]N為偶數(shù)

if(C1*D1.lt.0.)goto10

5

write(*,*)N，E

goto100

10

if(alog(-D1/C1)/(2*r2).gt.W)N=N+1

goto5

100

continue

end根據(jù)以上計算程序可得表13.1，該表給出了節(jié)點(diǎn)數(shù)隨能量的變化情況。從表中可知，為基態(tài)能級，可以再取Emax=-15.2，Emin=-15.6，M=51，可得到更為精確的結(jié)果E0=-15.412，…。同理可得第一激發(fā)態(tài)能級為

表13.1變量E與N的對應(yīng)值

【例13.2】考慮一維無限深勢阱，如圖13.2所示，勢函數(shù)為

13.3薛定諤方程的矩陣解法圖13.2一維無限深勢阱由量子力學(xué)基礎(chǔ)知識可知，在阱外，當(dāng)u為∞時，ψ=0。而在阱內(nèi)，波函數(shù)滿足微分方程

(13.15)

用有限差分法中的中心差分替代式(13.15)中的微分后可得

(13.16)

假設(shè)將x軸上［-a，a］區(qū)間四等分，等分點(diǎn)為i=0，1，2，3，4。根據(jù)波函數(shù)連續(xù)條件可知ψ0=ψ４＝０。以下僅就i=1，2，3進(jìn)行討論。將i=1，2，3分別代入上式，則有

(13.17(a))

(13.17(b))

(13.17(c))令，式(13.17(a))~(13.17(c))可以表示為以下矩陣方程

即

這是一個典型的本征方程Aψ=Eψ形式。利用線性代數(shù)中的有關(guān)知識，通過求解|λI-A|=0可以求得本征值λi。但一般來講，該方法主要用于低階陣，當(dāng)A的階數(shù)較高時，本征方程是一個關(guān)于λ的高次方程，可采用數(shù)值求解的方法求解本征值。方法主要有冪法、反冪法和雅可比法。其中雅可比法主要針對的是A為實(shí)對稱陣的本征方程(量子力學(xué)中A多為實(shí)對稱矩陣)。

定理13.1

若A為實(shí)對稱陣(方陣)，則存在正交矩陣R，使得

(13.18)雅可比法的思想是基于定理13.1，設(shè)法用一系列簡單的正交矩陣RK，逐步將A對角化，即選擇RK，令

(13.19)取A0=A，使當(dāng)K→∞時，AK→diag(λ1，λ2，…，λn)。設(shè)，取平面旋轉(zhuǎn)矩陣，則有

為了使非對角元素為0，即

只要選擇角θ，滿足即可。因此可以選取。

在此將二階平面旋轉(zhuǎn)矩陣進(jìn)行推廣。在AK-1中選取非對角元素模為最大的元素

(設(shè)p<q)，旋轉(zhuǎn)矩陣可以表示為

經(jīng)過旋轉(zhuǎn)變換，很容易將建立起AK中的元素與AK-1中的元素的對應(yīng)關(guān)系，即

(13.20)當(dāng)(2)當(dāng)，因此有

歸納結(jié)果如下:

(13.21)若，可取

設(shè)逐次所用的平面旋轉(zhuǎn)矩陣為R1，R2，…，Rk，則有

(13.22)

令

(13.23)則

(13.24)

Ak可以看做是一個對角陣(非主對角元素接近于零)，根據(jù)式(13.24)可得

(13.25)從而Vk的第j列向量就是矩陣A的本征值所對應(yīng)的本征向量，并且得到的本征向量系是標(biāo)準(zhǔn)正交系。記

V0=I

(13.26)

根據(jù)式(13.23)得到

(13.27)記Vk=［］，則

(13.28)

按照式(13.28)的迭代公式，在獲得本征值的同時，便可得到本征向量Vk。

需要說明的是，AK是中經(jīng)變換化為零的元素，在AK+1時又成非零元素，不能指望通過有限次旋轉(zhuǎn)變換就把A對角化，但利用范數(shù)理論可證:當(dāng)K→∞時，AK→diag(λ1，λ2，…，λn)，實(shí)際計算時，可取非對角元素近似為零時，K迭代即可結(jié)束。顯然在迭代過程中也無需知道旋轉(zhuǎn)矩陣RK的具體形式。

【例13.3】

計算對稱矩陣

的本征值及本征向量。

解①A0=A，選apq=a12=-1(p=1，q=2)。

由于a11=a22，因此有，從而，

。經(jīng)第一次變換，有

②選主元素=0.707107(p=1，q=3)，經(jīng)式(13.21)計算后有

C=0.707107，tanθ=0.517638

cosθ=0.888074，sinθ=0.459701

③選主元素=-0.627963(p=2，q=3)，同樣可以計算

④可見矩陣A的近似本征值為

λ1≈3.414209，λ2≈0.585986，λ3≈1.999800

而A的準(zhǔn)確本征值為

對應(yīng)的本征向量為

參考計算程序:

programmain

implicitdoubleprecision(T)

doubleprecision::pai=3.1415926535898795d0

doubleprecisionA(3,3),V(3,3)

integeri,j,p,q

dataA/2.0,-1.0,0.0,-1.0,2.0,-1.0,0.0,-1.0,2.0/

!迭代誤差

t_eps=1d-10

!迭代次數(shù)It=1

!初始化本征向量組

doi=1,3

V(i,i)=1.0

doj=1,3

if(i.ne.j)then

V(i,j)=0.0

endif

enddo

enddo!找出矩陣中非對角線模值最大的元素,以及其所處的行數(shù)和列數(shù)

10t_max1=0.0

doi=1,3

doj=1,3

if(i.ne.j)then

if(dabs(A(i,j)).gt.t_max1)then

t_max1=dabs(A(i,j)) p=i

q=j

endif

endif

enddo

enddo

!Step_1:求迭代系數(shù)

!A(p,p)與A(q,q)相等

if(A(p,p).eq.A(q,q))then

if(A(p,q).ge.0.0)then

theta=pai*0.25else

theta=-pai*0.25

endif

t_cos=dcos(theta)

t_sin=dsin(theta)

!A(p,p)和A(q,q)不相等

else

t_C=(A(p,p)-A(q,q))/(2.0*A(p,q))

if(t_C.ge.0.0)then

t=1.0/(dabs(t_C)+dsqrt(t_C**2+1.0))

else

t=-1.0/(dabs(t_C)+dsqrt(t_C**2+1.0))

endif

t_cos=1.0/dsqrt(1.0+t**2)

t_sin=t*t_cos

endif

!Step_2:特殊位置元素的迭代

t_pp=A(p,p)

t_qq=A(q,q)

t_pq=A(p,q)

A(p,p)=t_pp*t_cos**2+t_pq*2.0*t_sin*t_cos+t_qq*t_sin**2

A(q,q)=t_pp*t_sin**2-t_pq*2.0*t_sin*t_cos+t_qq*t_cos**2A(p,q)=(t_qq-t_pp)*t_sin*t_cos+t_pq*(t_cos**2-t_sin**2)

A(q,p)=A(p,q)

!Step_3:一般位置元素的迭代

doi=1,3

if((i.ne.p).and.(i.ne.q))then

t_ip=A(i,p)

t_iq=A(i,q)

A(i,p)=t_ip*t_cos+t_iq*t_sin

A(p,i)=A(i,p)

A(i,q)=-t_ip*t_sin+t_iq*t_cos

A(q,i)=A(i,q)

doj=1,3

if((j.ne.p).and.(j.ne.q))then

A(i,j)=A(i,j)

endif

enddo

endif

enddo!求本征向量

doi=1,3

t_fm=V(i,p)

V(i,p)=t_fm*t_cos+V(i,q)*t_sin

V(i,q)=-t_fm*t_sin+V(i,q)*t_cos

enddo

!Step_4:求誤差

t_max2=0.0

doi=1,3

doj=1,3if(i.ne.j)then

if(dabs(A(i,j)).gt.t_max2)then

t_max2=dabs(A(i,j))

end