《计算物理学》课件第13章

13.1定态薛定谔方程13.2一维方势阱中粒子能级和波函数的计算机求解13.3薛定谔方程的矩阵解法习题十三量子力学与经典力学的差别首先表现在对粒子的状态和力学量的描述及其变化规律上。在量子力学中，微观粒子的运动状态称为量子态，是用波函数ψ(r，t)来描述的，波函数所反映的微观粒子波动性，就是德布罗意波。德布罗意波或波函数ψ(r，t)不代表实际物理量的波动，而是描述粒子在空间的概率分布的概率波。13.1定态薛定谔方程量子力学中描述微观粒子的波函数本身是没有直接物理意义的，具有直接物理意义的是波函数模的平方，它代表了粒子出现的概率。波函数模的平方|ψ(r，t)|2代表时刻t，在r处附近空间单位体积中粒子出现的几率。因此|ψ(r，t)|2也被称为概率密度，即某一时刻出现在某点附近在体积元dV中的粒子的概率为|ψ(x，y，z，t,)|2dxdydz。波函数必须满足标准化条件，即单值、连续、有限；同时波函数还必须满足归一化条件ψ*(r，t)ψ(r，t)dτ=1。当微观粒子处于某一状态时，它的力学量(如坐标、动量、角动量、能量等)具有一系列可能值，每个可能值以一定的几率出现。当粒子所处的状态确定时，力学量具有某一可能值的几率也就完全确定了。量子力学中求解粒子问题常归结为解薛定谔方程或定态薛定谔方程。薛定谔方程广泛地用于原子物理、核物理和固体物理，对于原子、分子、核、固体等一系列问题中求解的结果都与实际符合得很好。薛定谔方程通常表示为

(13.1)即

(13.2)

其中式(13.2)为偏微分方程，这里u(r，t)称为势函数，ψ(r，t)称为波函数，m为粒子的质量。我们的目的是求能级E和波函数ψ。若u=u(r)，不含时间，则有

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

(13.3)

代入含时薛定谔方程(13.2)后，可得

(13.4)因此有

(13.5)根据式(13.3)有

(13.6)

式(13.5)中的第二式称为定态薛定谔方程，即

(13.7)在一维情况下有

(13.8)

通常又写为算子形式

(13.9)

显然该方程为一本征方程，E即为本征值，而ψ为本征函数。

对于一维方势阱问题，势函数u(x)=

如图13.1所示，试根据式(13.9)求解本征方程的本征值E(能级)和本征函数ψ(波函数)。13.2一维方势阱中粒子能级和波函数的计算机求解图13.1一维方势阱示意图

根据式(13.9)，有

即

为简单起见，记(不影响结果的一般性)，因此有

(13.10)阱内波函数满足(0≤x≤W)

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

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

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

这里j=0表示势阱左边的波函数，而j=1表示势阱右边的波函数。待定系数有A1、B1、C0、C1、D0、D1。

根据自然边界条件，x→±∞时，波函数ψ(x)→0，因此对于势阱左侧可得D0=0，而在右侧有C1=0。再利用边界条件，在x=0处波函数和一阶导数连续，有

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

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

C0为常数，一般可取为1。因此波函数可以表示为

(13.11(a))

(13.11(b))

(13.11(c))

同理，在x=W处，波函数和一阶导数连续，因此有

(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。结合方程(13.12(a))~(13.12(d))，可得

(13.14(a))

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

可以分别通过

来确定。要确定本征值E可以通过式(13.14(a))中的C1=0确定，这就要求对E取值有限制，必须取确定值。但这种方法一般不予采用，而通常采用计算节点法来实现。我们知道对势阱而言，一般有Emax=0(否则势阱不起作用，变为自由粒子)，Emin=V0。需要说明的是，与基态能量E0相应的波函数没有节点；第n个激发态相应的波函数ψn有n个节点，而且这些能量值恰好是波函数的节点变化时临界的能量值。我们可以通过计算不同能量E对应解波函数的节点数。确定节点变化时临界的能量值E即为能量本征值。

计算节点时，方势阱内节点数是用阱内的半波长数r1W/π(一般是介于某奇数与相邻偶数间的数)取整数后决定的，即若ψ(0)·ψ(W)>0则取相邻的偶数，若ψ(0)·、ψ(W)<0则取相邻的奇数作为节点数；方势阱外节点数是令ψ(x)=0，由式(13.13)解得的C1和D1异号，并且

，则有一节点，否则就无节点。计算节点和能级的步骤如下:

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

(2)计算

(3)利用以上节点计算方法确定节点数。

(4)由节点数计算结果来定出能级。①若相邻的两个能量分别对应0、1两个节点，则E1必处于Ea、Eb间，取

②若相邻的两个能量分别对应一个和两个节点，则E2必处于这两个能量间，依此类推。

【例13.1】

例如取V0=-20.0，W=1.0，Emin=-20.0，Emax=0.0，M取51。利用计算节点法编程计算出基态和第一激发态能级。

【例13.1】

例如取V0=-20.0，W=1.0，Emin=-20.0，Emax=0.0，M取51。利用计算节点法编程计算出基态和第一激发态能级。

计算程序:

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)!对应式(13.12(a))

phi1=r2*cos(W*r1)-r1*sin(W*r1)!对应式(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为奇数

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

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根据以上计算程序可得表13.1，该表给出了节点数随能量的变化情况。从表中可知，为基态能级，可以再取Emax=-15.2，Emin=-15.6，M=51，可得到更为精确的结果E0=-15.412，…。同理可得第一激发态能级为

表13.1变量E与N的对应值

【例13.2】考虑一维无限深势阱，如图13.2所示，势函数为

13.3薛定谔方程的矩阵解法图13.2一维无限深势阱由量子力学基础知识可知，在阱外，当u为∞时，ψ=0。而在阱内，波函数满足微分方程

(13.15)

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

(13.16)

假设将x轴上［-a，a］区间四等分，等分点为i=0，1，2，3，4。根据波函数连续条件可知ψ0=ψ４＝０。以下仅就i=1，2，3进行讨论。将i=1，2，3分别代入上式，则有

(13.17(a))

(13.17(b))

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

即

这是一个典型的本征方程Aψ=Eψ形式。利用线性代数中的有关知识，通过求解|λI-A|=0可以求得本征值λi。但一般来讲，该方法主要用于低阶阵，当A的阶数较高时，本征方程是一个关于λ的高次方程，可采用数值求解的方法求解本征值。方法主要有幂法、反幂法和雅可比法。其中雅可比法主要针对的是A为实对称阵的本征方程(量子力学中A多为实对称矩阵)。

定理13.1

若A为实对称阵(方阵)，则存在正交矩阵R，使得

(13.18)雅可比法的思想是基于定理13.1，设法用一系列简单的正交矩阵RK，逐步将A对角化，即选择RK，令

(13.19)取A0=A，使当K→∞时，AK→diag(λ1，λ2，…，λn)。设，取平面旋转矩阵，则有

为了使非对角元素为0，即

只要选择角θ，满足即可。因此可以选取。

在此将二阶平面旋转矩阵进行推广。在AK-1中选取非对角元素模为最大的元素

(设p<q)，旋转矩阵可以表示为

经过旋转变换，很容易将建立起AK中的元素与AK-1中的元素的对应关系，即

(13.20)当(2)当，因此有

归纳结果如下:

(13.21)若，可取

设逐次所用的平面旋转矩阵为R1，R2，…，Rk，则有

(13.22)

令

(13.23)则

(13.24)

Ak可以看做是一个对角阵(非主对角元素接近于零)，根据式(13.24)可得

(13.25)从而Vk的第j列向量就是矩阵A的本征值所对应的本征向量，并且得到的本征向量系是标准正交系。记

V0=I

(13.26)

根据式(13.23)得到

(13.27)记Vk=［］，则

(13.28)

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

需要说明的是，AK是中经变换化为零的元素，在AK+1时又成非零元素，不能指望通过有限次旋转变换就把A对角化，但利用范数理论可证:当K→∞时，AK→diag(λ1，λ2，…，λn)，实际计算时，可取非对角元素近似为零时，K迭代即可结束。显然在迭代过程中也无需知道旋转矩阵RK的具体形式。

【例13.3】

计算对称矩阵

的本征值及本征向量。

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

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

。经第一次变换，有

②选主元素=0.707107(p=1，q=3)，经式(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的准确本征值为

对应的本征向量为

参考计算程序:

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

!迭代次数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!找出矩阵中非对角线模值最大的元素,以及其所处的行数和列数

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:求迭代系数

!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