一维薛定谔方程的数值解法


1. 物理想法

参照《计算物理学》(Steven E.Koonin著,秦克诚译)第三章的相关内容,我写 了这篇笔记。

我们考虑的问题是求一个在一维位势~{V(x)}~中运动的质量为~{m}~的粒子的量子力 学定态解。我们假定~{V(x)}~的形状如图,是一个三角形势阱,

在~{x=x_{min}}~与~{x=x_{max}}~两点处位势变为无穷,这两点之间则有一个势阱。 我们假定定态不含时的薛定谔方程和边界条件是

\displaystyle  \begin{aligned} \frac{\mathrm{d} \phi}{\mathrm{d} x^2} + k^2(x) \phi(x) = 0; \\ \psi(x_{min}) = \psi(x_{max})=0, \end{aligned} \ \ \ \ \ (1)

其中

\displaystyle  k^2(x)= \frac{2m}{\hbar^2} [E - V(x)].

我们要求是非零的能量本征值。对这个本征值~{E}~来说,在经典力学容许的范围 内({E<V(x)})是振荡的,在经典力学禁戒的区域({E>V(x)})呈现指数行为。整体 来看,本征值小于零的解是束缚态的解,本征值大于零的解是连续解。我们这里 感兴趣的是束缚态。

要算出这个方程束缚态的本征值,我们需要利用打靶法计算。由于积分区域既有 振荡部分,也有指数部分,直接积分会引起数值上的不稳定。所以我们不能把波 函数直接从区域的左边积到右边,再利用边界条件判断是否是本征值。可是考虑 到波函数的连续性以及其导数的连续性,我们可以分别从左边积分到右边,再从 右边积分到左边,在振荡区和指数区的交界处,我们来判断波函数及其导数是否 连续。考虑到波函数归一化时我们总可以选择一种归一化系数,使得两种不同的 积分方式得到的波函数大小在交界处的值相等,所以我们只需要保证导数连续即 可。也就是下式要成立:

\displaystyle  \frac{\mathrm{d}\psi_<}{\mathrm{d}x}\Bigg\vert_{x_m}- \frac{\mathrm{d}\psi_>}{\mathrm{d}x}\Bigg\vert_{x_m} = 0. \ \ \ \ \ (2)

这里的~{\psi_<}~和~{\psi_>}~分别是向右和向左积分得到的波函数。用最简单 的有限差分格式代替这个导数式,我们得到下式:

\displaystyle  f \equiv \frac{1}{\psi} [\psi_<(x_m - h) - \psi_> (x_m - h)] = 0, \ \ \ \ \ (3)

由于归一化保证了~{\psi_< (x_m) = \psi_> (x_m)},上式中的~{\psi}~可以用~ {x_m}~点上的~{\psi_<}~来代替。如果不存在交界区(如具有上图所示位势,而~ {E>0}),那么~{x_m}~可以选择在任何区域,如果有多余两个转折点,那么就必 须把每个区域中的解拼在一起。

之所以用这种拼凑的解法,是因为我们只能保证积分在每个小区域中是准确的, 但是多个区域之间利用连续性条件我们可以凑成一个完整的波函数。

2. 算法分析

我们为了获得精确的积分结果,用到了~Numerov~算法。对薛定谔方程来说,我 们利用~Numerov~算法得到递推关系如下:

\displaystyle  \Bigg(1+\frac{\hbar^2}{12} k^2_{n+1} \Bigg) \psi_{n+1} - 2\Bigg( 1 - \frac{5\hbar^2}{12} k^2_n \Bigg) \psi_n + \Bigg( 1 + \frac{\hbar^2}{12} \Bigg) \psi_{n-1} = O(h^6). \ \ \ \ \ (4)

为了计算积分,我们必须利用泰勒级数算出头两项的数值出来。由于波函数的归 一是任意的,所以波函数在边界处的一阶导数可以随我们任意取,不影响最后得 到的本征值和波函数结果。

为了让计算机更好处理,我们把薛定谔方程无量纲化,写成

\displaystyle  \Bigg[ - \frac{1}{\gamma^2} \frac{\mathrm{d}^2}{\mathrm{d} x^2} + v(x) - \epsilon \Bigg] \psi(x) = 0, \ \ \ \ \ (5)

其中

\displaystyle  \gamma^2 = \frac{ 2m a^2 V_0}{\hbar^2}

是系统经典力学本性的一个无量纲测度,~{\epsilon = E/V_0}~是无量纲能量, 坐标由长度~{a}~标度。

我们用打靶法获得本征值。先从能量的最小可能取值开始往上搜索。如果发现~ {f}~变号,那么可以认为这里可能有一个本征态。于是粗略的获得本 征值的范围,然后减小步长,获得更加精确的结果,直到~{f}~的绝对值小于~ {10^{-6}}~为止。由于算法的问题,~{f}~变号不能一定保证就有本征值在附近。 为了排除这种情况,我们要联合考虑看是否随着搜索的进行~{f}~的绝对值在不断 减少。如果不是如此,那么我们就要去掉这个结果。在得到一个能量本征值之后,我 们以这个能量本征值为新的起点,继续前面的过程,可以获得一系列的能量本征值~ {E_i},直到能量本征值大于零为止。对每一个本征值对应的波函数,我们要分 别归一化。归一化时,我们用~Simpson~法求积分,比梯形法具有更高的精确度。

我们考虑一个具体的问题:半导体硅中低能谷电子有效质量为~{m^*=0.916 m_e},势场 宽为~{35 nm}~,外电场从1伏特变到零,电子在外场中的势能从~{-1}ev~变到零。 我们把积分区域分成500个格子,起步能量步长为~{1/500}~电子伏特,能量尺度 设为~0.1ev,长度标度为~{a=3}nm。最后算出能量本征值为~ -9.2426, -8.6758,-8.2117, -7.8016 , -7.4267, -7.0773 , -6.7477, -6.4340, -6.1336 ,-5.8444, -5.5649, -5.2940, -5.0307,-4.7741, -4.5238, -4.2790,-4.0393, -3.8044, -3.5738, -3.3473, -3.1246, -2.9055, -2.6897,-2.4770, -2.2673, -2.0604, -1.8562, -1.6544, -1.4551,-1.2581, -1.0633, -0.8705, -0.6794, -0.4889, -0.2965, -0.0984, 这里能量的单位是~{0.1}ev。可以看出,随着能级的升高, 能量之间的间隔越来越小。

为了确保算出的能量本征值是正确的,我们把两个不同能量本征值对应的波函数相乘 再积分,发现二者是正交的,正好满足量子力学中不通能量本征值的波函数必须 正交的性质。下面给出几个图,看看不同能量本征值下的波函数有什么特点,并 与前面理论分析的结果相互比较,以确保我们的计算无误。

如图,左边是第三个能级的波函数,右边是第八个能级的波函数。可以看出能级 数与振荡部分波峰波谷之和相等。在经典禁戒区,波函数是指数衰减的。与我们 前面分析的结论是一致的。

下面附有源程序

clear; m0=9.108E-31; ml=0.916*m0; h1=1.055e-34; N=500; M=500; tol=1e-6; a=3e-9; V0=1.602e-20; l=35e-9; dx=l/N; L=l/a; Dx=dx/a; X=0:Dx:L; x=X*a; ga=sqrt(2*ml*V0*a^2/h1^2); Vm=1; V=-Vm:Vm/N:0; V=V*1.602e-19; v=V/V0; eold=v(1)+abs(v(1))/M; for k=1:N/2 U(k,2*k-1:2*k+1)=[1,4,1]; end for k=1:1000 e1=eold+abs(v(1))/M; err=2*tol; v1=-(v-eold); s=find(v1tol & abs(err)<10000 & e1<0 & M<10^10 v1=-(v-e1); s=find(v1<=0);

if length(s)==1 e=e(find(e~=0)); E=e*V0; P=Psi(find(e~=0),:); return end K=-ga^2*(v-e1); psi(N+1)=0; psi(N)=Dx*0.5; for n=N+1:-1:s(1)+1 psi(n-2)=(2*(1-K(n-1)*5*Dx^2/12)*psi(n-1)-(1+K(n)*Dx^2/12)… *psi(n))/(1+K(n-2)*Dx^2/12); end psi1(1)=0; psi1(2)=Dx*0.5; for n=1:s(1)-1 psi1(n+2)=(2*(1-K(n+1)*5*Dx^2/12)*psi1(n+1)-(1+K(n)*Dx^2/12)… *psi1(n))/(1+K(n+2)*Dx^2/12); end psi=psi/psi(s(1)); psi1=psi1/psi1(s(1)); err=psi1(s(1)+1)-psi(s(1)+1); eold=e1; e1=e1+abs(v(1))/M; if err*errold<0 eold=e1-3*abs(v(1))/M; M=M*10; e1=eold+abs(v(1))/M; end end if abs(err)<0.0001 e(k)=eold; end Psi(k,1:s(1))=psi1(1:s(1)); Psi(k,s(1)+1:N+1)=psi(s(1)+1:N+1); Psi(k,:)=Psi(k,:)/sqrt(Dx*sum(U*Psi(k,:).^2')/3); psi=0; psi1=0; M=500; eold=eold+2*abs(v(1))/M; end

Advertisements

发表评论

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 更改 )

Twitter picture

You are commenting using your Twitter account. Log Out / 更改 )

Facebook photo

You are commenting using your Facebook account. Log Out / 更改 )

Google+ photo

You are commenting using your Google+ account. Log Out / 更改 )

Connecting to %s