有限差分法模拟一维谐振子

TOC\o"1-5"\h\z\o"CurrentDocument"第1章概述 4\o"CurrentDocument"第2章有限差分方法 5\o"CurrentDocument"2.1有限差分法基本思想 5\o"CurrentDocument"2.2差分方程组的求解 3\o"CurrentDocument"2.2.1高斯-赛德尔迭代法 32.2.2逐次超松弛法 6\o"CurrentDocument"第3章求解谐振子的微分方程 7\o"CurrentDocument"3.1一维谐振子 7\o"CurrentDocument"3.2二维各向同性谐振子 9\o"CurrentDocument"第4章总结 10\o"CurrentDocument"参考文献 11\o"CurrentDocument"附录 12\o"CurrentDocument"附1一维线性谐振子的程序设计 12\o"CurrentDocument"附1.1基态一维线性谐振子 12\o"CurrentDocument"附1.2第一激发态一维线性谐振子 13\o"CurrentDocument"附1.3第二激发态一维线性谐振子 13\o"CurrentDocument"附2二维线性谐振子的程序设计 14第1章概述微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似； 把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。有限差分法可广泛用来求解偏微分方程的近似解，在电磁场中求解点位函数的拉普拉斯方程时，可采用有限差分法的基本思想是：用网格将场域进行分割，再把拉普拉斯方程用以各网格点处的点位作为未知数的差分方程式来进行代换，将求解拉普拉斯方程解得问题变为求联立差分方程组的解得问题［1］，在差分网格非常多和情况下，利用并行计算方法对其进行区域分解，每个进程负责运算一部分区域，区域边界之间进行必要地通信可有效提高计算速度，解决更大规模的问题。往往只讨论它在静态场中的应用，即泊松方程或拉普拉斯方程的有限差分形式，很少涉及到它在时谐场(即亥姆霍兹方程)中的应用。本文重点讨论亥姆霍兹方程的有限差分形式以及它在时谐场中的应用。同时，有限差分法(finitedifferencemethod)是基于差分原理的一种数值计算方法，在求解微分方程定解问题中广泛应用。有限差分法是以差分原理为基础的一种数值计算法。它用离散的函数值构成的差商来近似逼近相应的偏导数，而所谓的差商则是基于差分的应用的数值微分表达式。用离散的只含有有限个未知量的差分方程组去近似代替连续变量的微分方程和定解条件，并把差分方程组的姐作为威风方程定解问题的近似解.有限差分法可以处理几乎所有形式的势函数，且主程序不依赖于势函数的具体形式，对于多数两字体都可以进行相对准确的计算。因此，将有限差分法应用于量子力学本征值问题的计算，有助于相对准确地进行量子体系和形象直观地教学研究［23］。量子力学教程中队一维无限深势阱、线性谐振子、氢原子等量子体系的薛定谔方程进行了严格的求解，得到了描述体系状态的波函数和能量的精确解。多数量子体系的哈密顿算数比较复杂，薛定谔方程不能严格求解，因此，研究和发展薛定谔方程的数值计算方法具有重要意义⑷。

第2章有限差分方法2.1有限差分法基本思想有限差分法是解偏微分方程的主要数值方法之一，其基本思想是把连续的问题离散化，即首先对求解区域作网格剖分，用有限个网格节点代替连续区域；其次将微分算子离散化，从而把微分方程的定解问题化为代数方程组的求解问题，解方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。参照文献［5］，给出有限差分法数值计算的基本思想：（1） 区域的离散或子区域的划分。（2） 插值函数的选择。（3） 方程组的建立。（4） 方程组的求解。2.2差分方程组的求解利用有限差分要面临求解的问题。在实际应用中，可以采用逐次迭代的迭代方法求解。这里介绍两种常用的迭代方法：高斯赛-德尔迭代法和逐次超松弛迭代法。2.2.1高斯-赛德尔迭代法中（k+1）=1中（k+1）+中（k中（k+1）=1中（k+1）+中（k+1）+中（k）+中（k）+h2二i,j4 -(2-1)%I其中i,j=1,2,(2-1)Xi-1,j 'i,jT 'i+1,j 'i,j+1 8J反复迭代（k=0,1......），一直进行到对所有节点满足下列条件为止。(2-2)"以+1）一平以）＜W(2-2)式中，W是预定的最大允许误差。在高斯-赛德尔迭代中，网格节点一般按“自然顺序”排列，即先“从左到右”，再“从上到下”排列，如图1-1所示；

图1-1图1-1网格节点排列2.2.2逐次超松弛法则可得将式(2-1)在同一点上相邻两次迭代的差值计为r(n)(i,j)则可得、1R(n)(l、1R(n)(l,j)=中(k+1)—中(k)=•i，j l,j 4—(P(k+1)+中(k+1)+中(k)+中(k)+4"l—1,j i,j-1 l+1,j l,j+1P(k+1)=P(k)+ap(n)a4(2-4)P(k+1)=P(k)+ap(n)a4(2-4)= P(k+1)+P(k+1)+P(k)+P(k)+h2—7 4p(k)4"i—1，ji，j—1i+1，ji，j+1 8i，j|其中a称为“加速收敛因子”，是一个供选择的参数，其值在1<a<2之间，该方法的快慢与a有着明显的关系。实践表明，如果a选得好，可以较快的加速迭代的速度。经验表明正方形场域由正方形网格划分，每边节点数为(P+1)时，最佳的收敛因子为a0(2-5)

第3章求解谐振子的微分方程3.1一维谐振子对一维线性谐振子，其能量本证方程为h2d2甲1 广(3-1)+口①2x2Q=Eq(3-1)2曰dx22边界条件(3-2)去长度单位为力s'2(3-2)去长度单位为力s'2引入无量纲参量&=oxa=以①.:方，则式则式(3-1)化为无量纲形式令V=qG)，ii式令V=qG)，ii式(3-2)可写为二qG)=&2dg2q^=Xq(E)式(2-3)可写为db+g2h-，Ei+i=w(3-3)(3-4)V0=0,V0=0,v=0(3-5)考虑式(3-5)，式(3-4)可写为SVi以Vi，其中Vm=2VJ、m-1/,2,2、+g2+101wx2（A&*也2211w2 +g（Ag）2 m-21（A&）2S为三对角矩阵用Matlab编程计算可同时得到矩阵S

步长值时一维谐振子能量S为三对角矩阵用Matlab编程计算可同时得到矩阵S

步长值时一维谐振子能量E的数值计m安N=0N=1N=2500.20.99752.98744.96731000.10.99942.99694.99192000.050.99982.99924.9980表1一维谐振子能量（以"2单）计算结果比较精确，且计算精度随求解区域内格点数目增加而增高，主量子1 3 5数为0、1、2时。能量精确解分别为-方①、-方①、-方①，简并度f为1.表1」 」 」中N为量子数，取m=200，△户0.05时波函数图形如图3-1所示.表1中N表1中N为量子数，取m=200，^e=0.05时波函数图形如图3-1所示，Matlab程序见附1.1，附1.2，附1.3.图3-1一维谐振子波函数3.2二维各向同性谐振子对于二维各向同性谐振子，薛定谔方程为边界条件一宗言+等)+2岫2边界条件一宗言+等)+2岫2+y"网XT3,中G,y)=0；yT3,中(x,y)=0(3-6)(3-7)取长度单位呻，能量单位为力W；2，引入无量纲参量&=ax,匚=ay,a=寸呻：力，入=E；2力①,则式(3-6)可写为(3-8)+岂邱,・)+£*2从,・)=抑&,°)(3-8)令^=中(，°)&=i-—△&,°=j-—A°,1<i<m-1;1<j<n-1,式(3-8)可写1,j ijiI2JjI2J为：一r八巾-r：\© +12La：\+h八1+(&2+°2)0-© -f人八巾 =叫(A°力 i,j-1 (A&力 i-1,j |_ 1也&力(A°>fe J J i,j 也&力 i+1,j 也°力 i,j+1 i,j(3-9)式(3-7)可写为80,=0,8.0=0,0,=0,8.=0 (3-10)则在是(2-10)下，式(2-11)可写为D©.1+C©.+D©.1=人0. (3-11)102\r©102\r©)0,1,j©2,j:，©=:1j©(A°)2/©—-Lj"m-1,jJC=2^^-+二］+(2C=2^^-+二］+(2+C@)1j1一@2)0)_山:晶通P01'2J+0+C2/1j输 再+@1一也1)(尸+^2+匚2

1j/则式（3-11）可进步写成矩阵式=人e，其中jjfCD0 0)DC00S二：:::00.CDk00.DC)ks为三对角块矩阵，是大型稀疏矩阵，用Matlab编程计算可同时得到稀疏矩阵s的本征值和本征矢.表2给出了取不同格点、步长值时二维谐振子能量E的数值计算结果与精确解的比较，从表2可以得到与表1相似的结论，主量子数为0、1、2时，能量精确解分别为力①、2力①、3力①，简并度分别为1、2、3.mn氐N=0N=1N=220200.50.51.96953.90465.777020200.50.5—3.91765.805420200.50.5——5.881040400.250.251.99593.98005.952540400.250.25—4.01455.987040400.250.25——6.116880800.1250.1252.00404.00005.993380800.1250.125—4.05076.043780800.1250.125——6.2166表2二维谐振子的能量（以"2单位）取m=n=20，Ae=A^0.5二维谐振子的波函数图形如图3-3所示，其中横向的2个坐标分别为无量纲参量旨。纵向坐标为波函数平。,。），Matlab程序见附2。-0；25图3-3二维谐振子波函数第4章总结计算物理学中，有限差分法（简称差分法），它以概念清晰，方法简单，直观的特点应用于数值分析领域，因而应用广泛。无论是常微分方程还是偏微分方程，各种类型的二阶线性方程，以全高阶或非线性方程，均可利用差分法转换为代数方程组，然后利用计算机求其数值解。有限差分法是以差分原理为基础的一种数值解法，通过学习计算物理学的有限差分法，我们对很多问题都可以得到足够高的计算精度。有限差分法是我们较容易掌握的数值解法，我们在物理学中，可以利用有限差分法解任何偏微分方程。利用有限差分法解边值问题时，首先将求解区域分为很多个网格和节点，并用差商代替微商，然后，使区域中的偏微分方程转化为以节点的数值为未知量的差分方程组，最后，解该方程组便可得到各离散点待求的数值解。该数值解是近似解，但逼近区域的真实解。在我们的实验设计中，如果离散化的点选择的足够密的话，我们得到的数值解与真实解的误差就能减小到可接受的程度。参考文献[1]吴连坳，井孝功，丁慧明等.径向薛定谔方程的有限差分解法[J],吉林大学自然科学报，1994，(03):67-70.⑵张志涌.精通MATLAB6.5北京：北京航空航天大学出版社，2003.刘建军，翟利学.有限差分法解能量本征方程，北京工业大学学报，2008,(34)：325-328.曾谨言.量子力学[M].北京：科学出版社,，2000.附录附1一维线性谐振子的程序设计附1.1基态一维线性谐振子clcclearM=200;H=20;h=H/(2*M);R=1./h人2;fori=1:2*M-1forj=1:2*M-1A(i,j)=0;endendforp=1:2*M-1A(p,p)=2*R+((p-M)*h)人2;endfori=1:2*M-2A(i,i+1)=-R;endfori=2:2*M-1A(i,i-1)=-R;end[v,d]=eig(A);d=eig(A);y=v(:,1);z=linspace(-H/2,H/2,2*M-1);y11=y.A2;M2=trapz(z,y11);y12=1/M2*y11;%plot(z,y12);y=(1/M2)A0.5*y;plot(z,y);holdon附1.2第一激发态一维线性谐振子clcclearM=200;H=20;h=H/(2*M);R=1./hA2;fori=1:2*M-1forj=1:2*M-1A(i,j)=0;endendforp=1:2*M-1A(p,p)=2*R+((p-M)*h)A2;endfori=1:2*M-2A(i,i+1)=-R;endfori=2:2*M-1A(i,i-1)=-R;end[v,d]=eig(A);d=eig(A);y=v(:,2);z=linspace(-H/2,H/2,2*M-1);y11=y.A2;M2=trapz(z,y11);y12=1/M2*y11;%plot(z,y12);y=(1/M2)A0.5*y;plot(z,y);holdon附1.3第二激发态一维线性谐振子clcclearM=200;H=20;h=H/(2*M);R=1./hA2;fori=1:2*M-1forj=1:2*M-1A(i,j)=0;endendforp=1:2*M-1A(p,p)=2*R+((p-M)*h)A2;endfori=1:2*M-2A(i,i+1)=-R;endfori=2:2*M-1A(i,i-1)=-R;end[v,d]=eig(A);d=eig(A);y=v(:,3);z=linspace(-H/2,H/2,2*M-1);y11=y.A2;M2=trapz(z,y11);y12=1