一维谐振子Schrdinger定态波动方程的数值求解

一维谐振子Schrdinger定态波动方程的数值求解..o作者：丁长林，指导教师：彭芳麟单位：北京师范大学物理系关键词：一维谐振子数值解摘要：本文运用科学数学软件MATLAB6.1研究一维谐振子定态波动方程的数值求解。一维谐振子定态解是属于有边界条件的常微分方程问题，MATLAB软件中的bvp4c指令正好可以用来求解物理上的这一类问题，本文就用它求出了一维谐振子的本征值和本征函数的数值解，与我们已经知道的解析解非常接近。全文分四部分：一、一维谐振子定态波动方程的处理；二、有关的MATLAB指令；三、数值求解程序；四、用户界面设计。AbstractInthispaper,Iinvestigatethenumericalsolutionoftheone-dimensionalharmonicoscillatorbyusingthemathematicsoftwareMATLAB6.1.Theproblemaboutthestationarystateoftheone-dimensionalharmonicoscillatorbelongstoBVPs(boundaryvalueproblems)ofODE(ordinarydifferentialequation),whichcanbesolvedbyusingthefunctionBVP4CofMATLABonthebeam.Igetthenumericalsolutionoftheone-dimensionalharmonicoscillatorfromgroundstatetotheeleventhexcitedstate,andcontrastitwiththeanalyticsolutionknowntousineigenvalueandeigenfunction,theyaremostadjacent.Thepaperconsistsoffourparts.Thefirstpartisthesimplificationofstationarystateequationoftheone-dimensionalharmonicoscillator;ThesecondpartisabouttheusageofsomeMATLABfunctionsneededinprogramming;Thethirdpartisthe――――――――――――――――――――毕业论文――――――――――――――――――――――programme;Inthefourthpart,Idesignedthegraphicaluserinterface.一、一维谐振子定态波动方程的简化处理一维谐振子定态波动方程为：[+−2222)2/(dxdmhπV(x)]Ψ(x)=EΨ(x)-----------------(1)势函数V(x)=21m2ω2x边界条件是：|x|∞时，Ψ(x)0.(1)式可化为：22dxdΨ(x)＋22/2π）（hm[E－21m]Ψ(x)＝0-----------------(2)2ω2x为便于数值求解，取m＝1，2ω22/2π）（hm＝50.则hω/2π=22/50（能量单位）.（2）式化为可数值求解的常微分方程：22dxdΨ(x)＋50（E－2212x）=0--------------------------------------------------------(3)在数值计算中，不可能利用∞处边界条件来求解，因此应作一定的处理。如图：从物理上看，要求的是粒子在势阱内的束缚态，预期本征函数在经典物理学容许的区域（EV(x)）内，波函数是振荡的，在经典物理学禁戒的区域（EV(x)）内，波函数是呈指数衰减的。若所求的能量本征值0EV(±1)=0.5，在|x|=1时，Ψ(x)0.因此，可将边界条件简化为：Ψ(－1)＝0，Ψ(＋1)＝0-------------------------------------------------------------------------(4)这样，常微分方程（3）和边界条件（4）构成的常微分方程边值问题就可利用数学软件MATLAB中的bvp4c指令来进行数值求解。二、有关的MATLAB指令1。bvp4c指令bvp4c可以用来解决具有两个边界条件的常微分方程问题，求解的区间为[a,b],边2――――――――――――――――――――毕业论文――――――――――――――――――――――界是a和b。指令的语句格式为：sol=bvp4c(@odefun,@bcfun,solinit)sol=bvp4c(@odefun,@bcfun,solinit,options)sol=bvp4c(@odefun,@bcfun,solinit,options,p1,p2,…)其含义为：odefun----求解的常微分方程的函数名，方程的形式为y’=f(x,y)。例如常微分方程y”=f(x,y,lambda)(lambda为未知参数，即所需求的本征值，以下同)，MATLAB要求先将方程写成一阶的常微分方程组，令y(1)=y,y(2)=dy/dx,得odefun就应写成{)2(/)1(),,(/)2(ydxdylambdayxfdxdy==functionddydx=[y(2)ydx=odefun(x,y,lambda)f(x,y,lambda)]bcfun----关于边界条件的函数bc(ya,yb),若有未知参数lambda，则应为bc(ya,yb,lambda)。由于bvp4c的核心是在区间[a,b]上求形式为y’=f(x,y)的常微分方程的积分，故边界条件函数中还要给出函数在起始点a的导数值，以便于积分。Options----用指令bvpset建立的优化选项，通常用默认值则不必输入，详细细节请查看bvpset和bvpget的帮助信息,这里不再赘述。p1,p2…----传递已知参数p1,p2,…给函数odefun,bcfun以及在options中定义的函数，这时用[]表示默认的优化选项。solinit----用指令bvpinit建立的结构数组。sol----是一个关于常微分方程数值解的结构数组。用bvp4c指令可得到在区间[a,b]上连续的解，并且其一阶导数连续，这正好符合物理上对波函数的要求。sol.x为bvp4c指令选择的把区间[a,b]分化的数据网格；sol.y为在这些网格点上的近似函数值y(x);sol.yp为在这些网格点上的函数的近似导数值y’(x);sol.parameters为经过数值计算求得的本征值；sol.solver为求解常微分方程的指令bvp4c.2。bvpinit指令bvpinit用来为bvp4c建立初始的猜测解solinit。该指令的语句格式为：solinit=bvpinit(x,yinit)solinit=bvpinit(x,yinit,lambda)solinit=bvpinit(x,yinit,lambda,p1,p2…)其含义为：x----对区间[a,b]初始网格化所得的有序节点，规定x(1)=a,x(end)=b,x(1)x(2)…x(end)。指令bvp4c在求解过程中还有可能改变这个数据网格，通常取x=linspace(a,b,10)就足够了，但网格节点应在解变化迅速的地方多取一些。若选择不当，就有可能得不到正确的解或出现如下的Warning:Unabletorefinethemeshanyfurther--theJacobianofthecollocationequationsissingular.yinit----对y(1),y(2)提供初始的猜测，y(1)的猜测尤为关键，y(2)随便取一个通常的3――――――――――――――――――――毕业论文――――――――――――――――――――――函数即可，因任意一个函数均可泰勒展开，故通常猜测y(1)为：y(1)=ax+ax+…+a1x+ann1−n1−n0这样的函数形式。y(1),y(2)也可猜测为常数，此时yinit为一个常数矩阵。lambda------边界值问题（BVP）中的未知参数，必须为其提供一个猜测值。p1,p2…---传递给猜测函数yinit(x,p1,p2,…)的已知参数，当无未知参数lambda时，调用的语句格式为solinit=bvpinit(x,yinit,[],p1,p2…)。参数p1,p2…只有在yinit为函数时才使用。3。deval指令deval用来求常微分方程的数值解在指定自变量数值处的值。该指令的语句格式为：sxint=deval(sol,xint)其含义为：sol----调用指令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb或bvp4c来求解常微分方程时所返回的结构数组，sol.solver就是所调用的指令名，对应有序的行矢量sol.x中的元素sol.x(i)，sol.y(:,i)为常微分方程的数值解在此处的值（通常sol.y(1,i)为函数值和sol.y(2,i)为函数的导数值）。这个结构数组还包含着在(sol.x(i),sol.x(I+1))里进行插值所需要的数据。这种数据的形式依赖于产生sol所调用的指令。xint----指定的自变量数值。sxint-------sxint(1,i)为自变量在xint(i)处的函数值，sxint(2,i)为自变量在xint(i)处的函数的导数值。三、数值求解程序下面是数值求解程序。初始的猜测解为Ψ(x)＝a,a,b,c,d为函数oscillator的输入变量，即不同的a,b,c,d给出不同的初始猜测解。p在子函数mat4bc中定义为全局变量，因此在MATLAB的指令窗（CommandWindow）中也要将p定义为全局变量并赋值后，才可调用函数oscillator。不同的p,a,b,c,d就有可能得到不同的本征值和本征函数。dxbxx+++23functionoscillator(a,b,c,d)lambda=0;％初始猜测值，求得的lambda的本征值将在0附近solinit=bvpinit(linspace(-1,1,10),@mat4init,lambda,a,b,c,d);％为bvp4c建立初始的猜测解solinitsol=bvp4c(@mat4ode,@mat4bc,solinit);％求解一维谐振子定态波动方程，sol.y就是Ψ(x)，sol.yp就是Ψ’(x)，sol.parameters为经过数值计算求得的lambda的本征值fprintf('Theeigenvalueisapproximately%7.5f.\n',(sol.parameters)*25*2^(1/2));％输出以hω/2π为单位的能量本征值xint=linspace(-1,1,250);％把区间[-1,1]再细化成251个等间隔的节点，以便于画图Sxint=deval(sol,xint);％插值法求在这251个节点处的解，Sxint(1,i)就是Ψ(xint(i)),Sxint(2,i)就是Ψ’(xint(i))4――――――――――――――――――――毕业论文――――――――――――――――――――――figure;％打开一个新的图形窗口plot(xint,Sxint(1,:));％在打开的图形窗口中作本征函数Ψ(x)的图象title('Eigenfunctionofharmonicoscillatorofonedimension');％在图的上端加标题Eigenfunctionofharmonicoscillatorofonedimensionxlabel('X');％加上x轴的标注ylabel('Ψ(X)');％加上y轴的标注functiondydx=mat4ode(x,y,lambda)dydx=[y(2)-50^2*(lambda-x^2/2)*y(1)];%一维谐振子定态波动方程functionyinit=mat4init(x