有限差分法

1有限差分法求解SchrÖdinger方程摘要:针对量子力学中大多数量子体系的哈密顿算符都比较复杂,薛定谔方程均不能得到严格解或分析解的问题,提出了用数学中的有限差分法来解决计算量子力学中薛定额方程的本征问题.对普通的径向薛定谔方程和含时的薛定谔方程进行了有限差分法的分析,给出了两种薛定谔方程的有限差分法的离散方程,并以线性谐振子为例,进行了计算机编程推算.结果表明,该方法在研究量子力学问题中具有广泛的应用前景.关键词:有限差分法;本征值;波函数2引言在量子力学中,对于一些简单的量子体系,如一维无限深势阱、线性谐振子、氢原子等体系的薛定谔(SchrÖdinger)方程可以严格求解,得到描述体系的精确的状态波函数和能量.但大多数量子体系的哈密顿算符(Hamiltonian)都比较复杂,薛定谔(SchrÖdinger)方程均不能得到严格解或解析解,有时能级可以给出解析表达式,但却无法得到波函数的解析表达式.因此,研究和发展薛定谔(SchrÖdinger)方程的计算方法就具有重要的意义.由于有限差分法可以处理几乎所有形式的势能函数,并且编制的计算机主程序不依赖于势能函数的具体形式.因此可以进行相对精确的计算,对于确定量子体系的能级和相应的波函数有着重要的意义.3第一章有限差分法1一维薛定谔方程的有限差分解法根据有限差分法中的二阶微分中心差分算符(其中忽略(3x)及更高阶项)22212dfxfxxfxfxxdxx(1.1)可以将一维薛定谔方程2222dxVxxExdx(1.2)化为2222xxxxxxVxEx(1.3)以点0,1,2,3,,mxamxmM,baxM(其中,ab为左右边界点,计算时将此点波函数设为零)将坐标分为M个相等的间隔,当M充分大时,x就足够小.将第m个分点的波函数简记为mmx.这样,(1.3)式可以化简为21122mmmmmxE(1.4)式中2222mmxVx令22,22mmRRVxx(1.3)可改写为11mmmmmRRE(1.5)2边界条件与计算区间当取0,1,2,3,,.mM并且注意到满足条件00M,则由(5)式得到一系列线性方程式11211223223343322122111MMMMMMMMMRERRERRERRERE这样将本征值方程离散化为矩阵方程SE(1.6)其中411223322110000000000000,000000000MMMMRRRRRSRRR(1.7)我们将相对复杂的方程就转化为矩阵S的对角化问题,利用Matlab可以容易将矩阵S的本征值和本征函数同时求出.薛定谔方程一维束缚态解的特征是当x时,波函数0.实际情况是在有限远处波函数已经为零,由其对低能级的状态而言.当能量逐渐增大,波函数不为零的区域也在逐渐扩大.针对某一能级计算时区间的选取要满足大于波函数不为零的区域,也就是说,在未达到计算区间处波函数应该已经为零,这才能保证满足束缚态要求的边界条件.为保证这点,当计算得到的波函数仅仅是在计算区间的起始和终点为零而不是已经为零,就要进一步扩大计算区间再进行计算,一直到在未到计算的始终点位置波函数已经明显为零.并且这时候增大或减少计算区间(保证计算区间包含了波函数不为零的区域)的长度不再影响本征值和波函数不为零区域的形态.条件00M相当于在计算起点和终点处存在一无限高势垒,只有当它足够远时才不会影响所计算粒子的所处的势能及粒子的运动状态.这里的讨论不包括计算无限深势阱中运动的粒子,因为这种情况计算的起点和终点处实际上就存在无限高势垒.3波函数的归一化由有限差分法计算所得的波函数并不是归一化的波函数,为了求得归一化波函数,可以利用数值积分求得归一化系数,然后除以波函数.当然数值积分不可能像解析解一样积分区间取为,,但因为在计算区间之外波函数全为零,所以积分区域可以选择为0,Mxx.这样经过归一化后的波函数的曲线与解析解完全一致.5第二章简谐振子的有限差分解2.1算法我们以线性谐振子为例,线性谐振子的能量本征值方程为22222022dExdx(2.1)为方便起见,引入没有量纲的变量代替x,它们的关系是,xx(2.2)并令2E(2.3)薛定谔方程可改写为下面无量纲的形式222dd(2.4)令,,1mmmmammM,考虑到(1.3)式和边界条件,0xx(2.5)(2.4)式可写为211222121mmmmm(2.6)方程(2.6)式,实际上可以将(1.5)式两边同乘2得到.因为2222222222211,2222222222mmmmRxxRVxxE边界条件(2.5)可写为00,0M(2.7)考虑到(2.7),式(2.6)可写为S其中62122222221223222322212222122210000012100001210000,12100001200000MmMm(2.8)2.2计算区间的选择如果我们取计算区域为2,2,这相当于粒子在图2.1所示的势场内运动.图1程序:clearx1=-2:0.01:2;y1=x1.^2;x2=-2;y2=x2.^2:0.001:5;x3=2;y3=x3.^2:0.001:5;plot(x1,y1,'-k',x2,y2,'-k',x3,y3,'-k','LineWidth',2)%第n个本征函数axis([-4405])set(gca,'XTick',-5:1:5)set(gca,'YTick',0:0.5:5)xlabel('\fontsize{14}\rm\xi','Color','k')ylabel('\fontsize{14}\rmV(\xi)','Color','k')而在此势场中运动的粒子与波函数可由如下程序计算得到:cleara=-2;b=2;N=1000;M=N+1;Dx=(b-a)/N;form=1:Mx(m)=a+(m-1)*Dx;end图2.1计算区间取为2,2时,粒子实际所处的势.7F=zeros(M);form=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endform=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);forn=1:1E(n)=d(n,n)/2%第n个本征值figureplot(x,v(:,n),'-k','LineWidth',2)%第n个本征函数axis([-2200.06])set(gca,'XTick',-2:0.5:2)set(gca,'YTick',0:0.01:0.06)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];gridonxlabel('\fontsize{14}\rm\xi','Color','k')ylabel(t,'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')end这时得到的本征值为0.5369,而波函数如图2所示,我们看到在计算的边界点上波函数才为零,而在边界点内波函数并不为零.显然本征值与波函数的形态与解析解相差较远,问题的关键就是所取的计算区间太小造成的.如果计算的区间选为4,4,再重新计算其程序修改为:cleara=-4;b=4;N=2000;M=N+1;Dx=(b-a)/N;form=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);图2计算区间取为2,2时,粒子波函数的形状.8form=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endform=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);forn=1:1E(n)=d(n,n)/2%第n个本征值figureplot(x,v(:,n),'-k','LineWidth',2)%第n个本征函数axis([-4400.05])set(gca,'XTick',-4:1:4)set(gca,'YTick',0:0.01:0.05)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];gridonxlabel('\fontsize{14}\rm\xi','Color','k')ylabel(t,'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')end得到的本征值为0.5000,本征函数的曲线如图3所示,这时我们看到波函数在计算的边界点内已经为零,这说明我们在计算边界点的无限高势垒并未影响此能级粒子的运动,所以得到了与解析解一样的本征值,但本征函数与解析解还差一个常数因子,这是因为有限差分法给出的波函数并没有归一化.实际上现在我们适当改变计算区间并不改变其本征值与本征函数的形态.下面我们再选计算区间为5,5重新计算,其计算程序为:cleara=-5;b=5;N=2500;M=N+1;Dx=(b-a)/N;form=1:Mx(m)=a+(m-1)*Dx;图3计算区间取为4,4时,粒子波函数的形状.图4计算区间取为5,5时,粒子波函数的形状.9endF=zeros(M);form=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endform=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);forn=1:1E(n)=d(n,n)/2%第n个本征值figureplot(x,v(:,n),'-k','LineWidth',2)%第n个本征函数axis([-5500.05])set(gca,'XTick',-5:1:5)set(gca,'YTick',0:0.01:0.05)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];gridonxlabel('\fontsize{14}\rm\xi','Color','k')ylabel(t,'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')end得到的本征值为0.5000,本征函数的曲线如图4所示,这时我们看到波函数在计算的边界点内已经为零的区域和图3基本一致.如果把两张图绘制在一起,发