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

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

(13.1)即

(13.2)

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

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

(13.3)

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

(13.4)因此有

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

(13.6)

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

(13.7)在一維情況下有

(13.8)

通常又寫(xiě)為算子形式

(13.9)

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

對(duì)于一維方勢(shì)阱問(wèn)題，勢(shì)函數(shù)u(x)=

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

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

即

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

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

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

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

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

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

根據(jù)自然邊界條件，x→±∞時(shí)，波函數(shù)ψ(x)→0，因此對(duì)于勢(shì)阱左側(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

可以分別通過(guò)

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

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

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

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

(2)計(jì)算

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

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

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

【例13.1】

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

【例13.1】

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

計(jì)算程序:

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)!對(duì)應(yīng)式(13.12(a))

phi1=r2*cos(W*r1)-r1*sin(W*r1)!對(duì)應(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ù)以上計(jì)算程序可得表13.1，該表給出了節(jié)點(diǎn)數(shù)隨能量的變化情況。從表中可知，為基態(tài)能級(jí)，可以再取Emax=-15.2，Emin=-15.6，M=51，可得到更為精確的結(jié)果E0=-15.412，…。同理可得第一激發(fā)態(tài)能級(jí)為

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

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

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

(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))可以表示為以下矩陣方程

即

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

定理13.1

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

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

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

為了使非對(duì)角元素為0，即

只要選擇角θ，滿(mǎn)足即可。因此可以選取。

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

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

經(jīng)過(guò)旋轉(zhuǎn)變換，很容易將建立起AK中的元素與AK-1中的元素的對(duì)應(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可以看做是一個(gè)對(duì)角陣(非主對(duì)角元素接近于零)，根據(jù)式(13.24)可得

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

V0=I

(13.26)

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

(13.27)記Vk=［］，則

(13.28)

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

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

【例13.3】

計(jì)算對(duì)稱(chēng)矩陣

的本征值及本征向量。

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

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

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

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

C=0.707107，tanθ=0.517638

cosθ=0.888074，sinθ=0.459701

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

④可見(jiàn)矩陣A的近似本征值為

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

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

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

參考計(jì)算程序:

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!找出矩陣中非對(duì)角線模值最大的元素,以及其所處的行數(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