1、目 录第一章 概述 1第二章 有限差分方法 22.1 有限差分法基本思想 .22.2 差分方程组的求解 2第三章 求解谐振子的微分方程 43.1 一维谐振子 .43.2 二维各向同性谐振子 .6第四章 总结 9参考文献 10附录 111第一章 概述微 分 方 程 和 积 分 微 分 方 程 数 值 解 的 方 法 。 基 本 思 想 是 把 连 续 的 定 解 区 域用 有 限 个 离 散 点 构 成 的 网 格 来 代 替 ， 这 些 离 散 点 称 作 网 格 的 节 点 ； 把 连 续 定解 区 域 上 的 连 续 变 量 的 函 数 用 在 网 格 上 定 义 的 离 散 变 量 函
2、数 来 近 似 ； 把 原 方程 和 定 解 条 件 中 的 微 商 用 差 商 来 近 似 ， 积 分 用 积 分 和 来 近 似 ， 于 是 原 微 分方 程 和 定 解 条 件 就 近 似 地 代 之 以 代 数 方 程 组 ， 即 有 限 差 分 方 程 组 ， 解 此 方程 组 就 可 以 得 到 原 问 题 在 离 散 点 上 的 近 似 解 。 然 后 再 利 用 插 值 方 法 便 可 以 从离 散 解 得 到 定 解 问 题 在 整 个 区 域 上 的 近 似 解 。 有限差分法可广泛用来求解偏微分方程的近似解，在电磁场中求解点位函数的拉普拉斯方程时，可采用有限差分法的基本思
3、想是：用网格将场域进行分割，再把拉普拉斯方程用以各网格点处的点位作为未知数的差分方程式来进行代换，将求解拉普拉斯方程解得问题变为求联立差分方程组的解得问题 ，在1差分网格非常多和情况下，利用并行计算方法对其进行区域分解，每个进程负责运算一部分区域，区域边界之间进行必要地通信可有效提高计算速度，解决更大规模的问题。往往只讨论它在静态场中的应用，即泊松方程或拉普拉斯方程的有限差分形式，很少涉及到它在时谐场（即亥姆霍兹方程）中的应用。本文重点讨论亥姆霍兹方程的有限差分形式以及它在时谐场中的应用。同时，有限差分法（finite difference method）是基于差分原理的一种数值计算方法，在求
4、解微分方程定解问题中广泛应用。有限差分法是以差分原理为基础的一种数值计算法。它用离散的函数值构成的差商来近似逼近相应的偏导数，而所谓的差商则是基于差分的应用的数值微分表达式。用离散的只含有有限个未知量的差分方程组去近似代替连续变量的微分方程和定解条件，并把差分方程组的姐作为威风方程定解问题的近似解.有限差分法可以处理几乎所有形式的势函数，且主程序不依赖于势函数的具体形式，对于多数两字体都可以进行相对准确的计算。因此，将有限差分法应用于量子力学本征值问题的计算，有助于相对准确地进行量子体系和形象直观地教学研究 。32量子力学教程中队一维无限深势阱、线性谐振子、氢原子等量子体系的薛定谔方程进行了严
5、格的求解，得到了描述体系状态的波函数和能量的精确解。多数量子体系的哈密顿算数比较复杂，薛定谔方程不能严格求解，因此，研究和发展薛定谔方程的数值计算方法具有重要意义 。42第二章 有限差分方法2.1 有限差分法基本思想有限差分法是解偏微分方程的主要数值方法之一，其基本思想是把连续的问题离散化，即首先对求解区域作网格剖分，用有限个网格节点代替连续区域；其次将微分算子离散化，从而把微分方程的定解问题化为代数方程组的求解问题，解方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。参照文献 ，给出有限差5分法数值计算的基本思想：（1） 区域的离散或
6、子区域的划分。（2） 插值函数的选择。（3） 方程组的建立。（4） 方程组的求解。2.2 差分方程组的求解利用有限差分要面临求解的问题。在实际应用中，可以采用逐次迭代的迭代方法求解。这里介绍两种常用的迭代方法：高斯-赛德尔迭代法和逐次超松弛迭代法。2.2.1 高斯-赛德尔迭代法采用正方形网格分割场域这个方法是，先对节点 选取初值 。其),(jiyx)0(ij中（0）表示 0 次近似值；下角标 表示节点所在位置，即第 行 列的交点。ji, i再按下式其中 (2-1)(1)(1)()()()2, ,1,14ijkkkkijijijijijh ,.21,ji反复迭代 ，一直进行到对所有节点满足下列条
7、件为止。.),0(2-2)Wkjikji)(,)1(,式中， 是预定的最大允许误差。W在高斯-赛德尔迭代中，网格节点一般按“自然顺序”排列，即先“从左到右” ，再“从上到下”排列，如图 2-1 所示；3yx图 2-1 网格节点排列1 2 3 45 6 7 82.2.2 逐次超松弛法将式(2-1)在同一点上相邻两次迭代的差值计为 ，则可得),(jinR(2-3)(,2)(1,)(,1)(,)1(,)(,)1(,)( 4, kjiijkjikjikjikjikjikjin hiR 按式(2-1)迭代理想的收敛情况是所有内点的点位函数的余数为零，但这是不可实现的，余数 时正时负，时大时小，当网格很大
8、时，收敛的速度很慢，),(ji这就需要改小网格，但迭代的次数将随之增加。逐次超松弛法就是在高斯-赛德尔迭代法中引入了加速收敛因子 对其进行校正，即(2-4) )(,2)(1,)(,1)(,)1(,)(,)(,)(, 44 kjijikjikjikjikjinjijiji h其中 称为“加速收敛因子” ，是一个供选择的参数，其值在 之间，该 21方法的快慢与 有着明显的关系。实践表明，如果 选得好，可以较快的加速迭代的速度。经验表明正方形场域由正方形网格划分，每边节点数为 时，)1(p最佳的收敛因子为 ，0(2-5)psin1/204第三章 求解谐振子的微分方程3.1 一维谐振子对一维线性谐振子
9、,其能量本证方程为(3-1)Exh221d边界条件, (3-2)x0去长度单位为 ，能量单位为 ，引入无量纲参量 ，2ax， ，则式(3-1)化为无量纲形式a21E= (3-3)2d2令 = ， = ，1 ，1，式(2-3)可写为iiimi(3-4)iiii 12212式(3-2)可写为 (3-5)0,0m考虑式(3-5)，式(3-4)可写为 ，其中iiS121mi5S= ， 2122222221 100 000 0mm S 为三对角矩阵，是大型稀疏矩阵，用 Matlab 编程计算可同时得到矩阵 S的本征值和本征矢表 3-1 是取不同格点、步长值时一维谐振子能量 E 的数值计算结果与精确解的比
10、较，从表 3-1 可以看出有限差分法m N=0 N=1 N=250 0.2 0.9975 2.9874 4.9673100 0.1 0.9994 2.9969 4.9919200 0.05 0.9998 2.9992 4.9980计算结果较精确，且计算精度随求解区域内格点数目增加而增高，主量子数为 0、1、2 时。能量精确解分别为 、 、 ，简并度 为 1。表21325f3-1 中 N 为量子数，取 m=200， =005 时波函数图形如图 3-1 所示表 3-1 中 N 为量子数，取 m=200， =005 时波函数图形如图 3-1 所示，表 3-1 一维谐振子能量(以 单位)2图 3-1
11、一维谐振子波函数6Matlab 程序见附 1.1，附 1.2，附 1.33.2 二维各向同性谐振子对于二维各向同性谐振子，薛定谔方程为(3-6)Eyxyx)(21)(22边界条件(3-7);0,0,yx取长度单位 ，能量单位为 ，引入无量纲参量2,ay， ,则式(3-6)可写为a21E(3-8) ,- 2令 式(3-8)可写 ,1;1,2, njminjmijijiji为： jijijijijiji ,1,12,22,121, )( (3-9)式(3-7)可写为(3-10)0,0,0 nijmij 则在是(2-10)下，式(2-11)可写为(3-11)jjjjDC11 jmjjD,12,122
12、100 ，7 21222122212221 000 000 jjjjC 则式(3-11)可进一步写成矩阵式 ，其中jjSCD 00 为三对角块矩阵，是大型稀疏矩阵，用Matlab编程计算可同时得到稀疏矩阵 的本征值和本征矢表3-2给出了取不同格点、步长值时二维谐振子能量E的数值计算结果与精确解的比较，从表3-2可以得到与表3-1相似的结论，主量子数为0、1、2时, 能量精确解分别为 简并度分别为l、2、3，、 m n N=0 N=1 N=220 20 0.5 0.5 1.9695 3.9046 5.777020 20 0.5 0.5 3.9176 5.805420 20 0.5 0.5 5.8
13、81040 40 0.25 0.25 1.9959 3.9800 5.952540 40 0.25 0.25 4.0145 5.987040 40 0.25 0.25 6.116880 80 0.125 0.125 2.0040 4.0000 5.993380 80 0.125 0.125 4.0507 6.043780 80 0.125 0.125 6.2166表 3-2 二维谐振子的能量(以 单位)28取 m=n =20， =0.5 二维谐振子的波函数图形如图 3-2 所示，其中横向的 2 个坐标分别为无量纲参量 ，纵向坐标为波函数 ，Matlab 程序,见附 2。图 3-2 二维谐振子波
14、函数9第四章 总 结计算物理学中，有限差分法（简称差分法） ，它以概念清晰，方法简单，直观的特点应用于数值分析领域，因而应用广泛。无论是常微分方程还是偏微分方程，各种类型的二阶线性方程，以至高阶或非线性方程，均可利用差分法转换为代数方程组，然后利用计算机求其数值解。有限差分法是以差分原理为基础的一种数值解法，通过学习计算物理学的有限差分法，我们对很多问题都可以得到足够高的计算精度。有限差分法是我们较容易掌握的数值解法，我们在物理学中，可以利用有限差分法解任何偏微分方程。利用有限差分法解边值问题时，首先将求解区域分为很多个网格和节点，并用差商代替微商，然后，使区域中的偏微分方程转化为以节点的数值
15、为未知量的差分方程组，最后，解该方程组便可得到各离散点待求的数值解。该数值解是近似解，但逼近区域的真实解。在我们的实验设计中，如果离散化的点选择的足够密的话，我们得到的数值解与真实解的误差就能减小到可接受的程度。10参考文献1 吴连坳，井孝功，丁慧明等. 径向薛定谔方程的有限差分解法J. 吉林大学自然科学报. 1994，(03) ：67-70.2 张志涌. 精通 MATLAB6.5M. 北京航空航天大学出版社，2003.3 刘建军，翟利学. 有限差分法解能量本征方程J. 北京工业大学学报. 2008，(34)： 325-328.4 曾谨言. 量子力学M. 北京：科学出版社，2000.11附录附
16、 1 一维线性谐振子的程序设计附 1.1 基态一维线性谐振子clcclearM=200; H=20; h=H/(2*M);R=1./h2; for i=1:2*M-1for j=1:2*M-1A(i,j)=0;endendfor p=1:2*M-1A(p,p)=2*R+(p- M)*h)2;endfor i=1:2*M-2A(i,i+1)=-R;endfor i=2:2*M-1A(i,i-1)=-R;endv,d=eig(A); d=eig(A); y=v(:,1); z=linspace(-H/2,H/2,2*M-1);y11=y.2;M2=trapz(z,y11);y12=1/M2*y11
17、;%plot(z,y12);y=(1/M2)0.5*y;plot(z, y);hold on12附 1.2 第一激发态一维线性谐振子clcclearM=200; H=20; h=H/(2*M);R=1./h2; for i=1:2*M-1for j=1:2*M-1A(i,j)=0;endendfor p=1:2*M-1A(p,p)=2*R+(p-M)*h)2;endfor i=1:2*M-2A(i,i+1)=-R;endfor i=2:2*M-1A(i,i-1)=-R;endv,d=eig(A); d=eig(A); y=v(:,2); z=linspace(-H/2,H/2,2*M-1);y
18、11=y.2;M2=trapz(z,y11);y12=1/M2*y11;%plot(z,y12);y=(1/M2)0.5*y;plot(z, y);hold on附 1.3 第二激发态一维线性谐振子clcclearM=200; H=20; h=H/(2*M);13R=1./h2; for i=1:2*M-1for j=1:2*M-1A(i,j)=0;endendfor p=1:2*M-1A(p,p)=2*R+(p-M)*h)2;endfor i=1:2*M-2A(i,i+1)=-R;endfor i=2:2*M-1A(i,i-1)=-R;endv,d=eig(A); d=eig(A); y=v
19、(:,3); z=linspace(-H/2,H/2,2*M-1);y11=y.2;M2=trapz(z,y11);y12=1/M2*y11;%plot(z,y12);y=(1/M2)0.5*y;plot(z, -y);hold on附 2 二维线性谐振子的程序设计ticclcclearL=20;W=20;N=20;M=20;hx=L/(2*N);hy=W/(2*M);S=zeros(2*M-1)*(2*N-1); for m=1:2*M-114D(m,m)=-1/(hy.2);endfor m=1:2*N-2mx=(2*M-1)*(m-1)+1;my=(2*M-1)*(m-1)+2*M-1;
20、nx=(2*M-1)*(m-1)+(2*M-1)+1;ny=(2*M-1)*(m-1)+2*(2*M-1);S(mx:my,nx:ny)=D; S(nx:ny,mx:my)=D; endfor n=1:2*N-1C=zeros(2*M-1);for m=1:2*M-1C(m,m)=2*(1/hx2+1/hy2)+(n-M)*hx)2+(m-N)*hy)2;endfor m=1:2*M-2C(m,m+1)=-1/hx.2;C(m+1,m)=-1/hx.2;endcx=(2*M-1)*(n-1)+1;cy=(2*M-1)*(n-1)+2*M-1;S(cx:cy,cx:cy)=C; endV,E=eig(S);Y=V(:,2); k=1;for j=1:2*M-1;for i=1:2*N-1;ZZ(i,j)=Y(k);k=k+1;endendy=linspace(-W/2,W/2,2*N-1);x=linspace(-L/2,L/2,2*M-1);x1,y1=meshgrid(x,y);surf(x1,y1,-ZZ);shading interp; zlabel();toc
