下载文档原格式
有限差分法finite difference method用差分代替微分，是有限差分法的基本出发点。
是一种微分方程和积分微分方程数值解的方法。
把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。
然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。
此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。
对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。
另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。
此外，还有一个重要的概念必须考虑，即差分格式的稳定性。
因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。
前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。
只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。
最常用的方法是数值微分法，比如用差商代替微商等。
另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。
此外还可以用待定系数法构造一些精度较高的差分格式。
龙格库塔龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
有限差分法有限差分法是数学领域的一项最新成果，它在某些特定情况下能得到非常好的结果。
所谓有限差分方程就是利用积分和求差公式将差分方程化成为多个等价的偏微分方程组的组合形式，然后再应用最优化方法求解这种方程组，从而得出未知数的近似值。
当已知方程组的每个参数及其变量代入数据计算后的误差时，只要对其进行必要的调整或者修改后，就可获得满意的精度与效率的估计值。
此外，还可以通过有限差分方程的求解来了解其物理背景。
比如说在物体碰撞问题中，两个质点之间距离的测量往往涉及到很复杂的三维几何关系。
即使是一个小的距离误差也会引起很大的误差。
因此，对于碰撞问题中两个质点之间的相互位置误差测量，必须考虑它们之间的三维几何关系，并根据具体问题建立相应的坐标系统。
有限差分方程可以用来描述许多不同类型的实际问题，例如质量、压力、速度、温度、流动、热传导、声音和电磁场等。
但是由于数学模型本身的复杂性，使得有限差分方程在求解上遇到了困难。
因此，人们开始寻找一种更加直观的方法来解决问题。
有限差分法正是基于此原理提出的。
利用有限差分方程求解偏微分方程，我们首先要给出所求解的偏微分方程的数学表达式，这样才能够在有限差分方程的数学模型中寻找解析解。
有限差分方程的解析解，需要借助解析函数的理论来确定。
但是在自然科学和工程技术领域里，对于一般的实际问题，很少会存在着某种数学模型完全适合于所有的具体问题，那么对于任意一个偏微分方程，总是存在着一个解析解。
当把偏微分方程的解析解用适当的坐标表示出来后，有限差分方程的求解就转化为如何寻找与这个解相对应的函数值的问题。
通常，解析函数的形式是比较复杂的，因此需要运用数值方法进行拟合，从而得到符合实际的数学表达式。
然后通过对这个数学表达式的求解来确定所求偏微分方程的解析解。
这种数值求解方法称为数值积分法。
在研究有限元法和边界元法时都可以采用一些简单易行而且计算机可能很容易处理的函数作为边界条件，而这些函数本身又是很容易计算的。
有限差分法推导【最新版】目录1.有限差分法的基本概念2.有限差分法的推导方法3.有限差分法的应用实例4.有限差分法的优缺点正文一、有限差分法的基本概念有限差分法是一种数值计算方法，主要应用于求解偏微分方程的初值问题。
它是通过将连续的函数值用有限个离散点上的函数值来代替，从而将偏微分方程转化为关于这些离散点上的代数方程组。
这种方法可以有效地降低问题的复杂度，使得求解过程更加简便。
二、有限差分法的推导方法有限差分法的推导过程主要包括以下几个步骤：1.对边界条件进行离散处理，将边界上的函数值用有限个离散点上的函数值来代替。
2.对偏微分方程进行离散处理，将偏微分方程转化为关于这些离散点上的代数方程组。
3.求解代数方程组，得到离散点上的函数值。
4.通过插值方法，将离散点上的函数值还原为连续函数。
三、有限差分法的应用实例有限差分法广泛应用于各种物理、工程和数学问题中，例如求解热传导方程、波动方程和亥姆霍兹方程等。
下面以求解一维热传导方程为例，展示有限差分法的应用过程。
假设我们要求解如下的热传导方程：u/t = k * ^2u/x^2x = [0, 1]t = [0, T]边界条件：u(0, t) = f(t), u(1, t) = 0初始条件：u(x, 0) = 0我们可以通过以下步骤应用有限差分法：1.对边界条件进行离散处理，将边界上的函数值用有限个离散点上的函数值来代替。
2.对偏微分方程进行离散处理，将偏微分方程转化为关于这些离散点上的代数方程组。
3.求解代数方程组，得到离散点上的函数值。
4.通过插值方法，将离散点上的函数值还原为连续函数。
四、有限差分法的优缺点有限差分法具有以下优点：1.适用范围广泛，可以应用于各种偏微分方程的初值问题。
2.推导过程相对简单，容易理解和实现。
3.计算精度较高，可以通过增加离散点数来提高精度。
然而，有限差分法也存在以下缺点：1.计算量较大，需要处理大量的代数方程组。
2.对于某些问题，可能需要进行特殊的处理，例如处理不稳定的代数方程组。
有限差分法有限差分法（Finite Differential Method, FDM ）什么是有限差分法 有限差分法是指用泰勒技术展开式将变量的导数写成变量，在不同时间或空间点值的差分形式的方法。
按时间步长和空间步长将时间和空间区域剖分成若干网格，用未知函数在网格结(节)点上的值所构成的差分近似代替所用偏微分方程中出现的各阶导数，从而把表示变量连续变化关系的偏微分方程离散为有限个代数方程，然后解此线性代数方程组，以求出溶质在各网格结(节)点上不同时刻的浓度。
有限差分法的基本步骤（1）剖分渗流区，确定离散点。
将所研究的水动力弥散区域按某种几何形状(如矩形、任意多边形等)剖分成网络系统。
（2）建立水动力弥散问题的差分方程组。
（3）求解差分方程组。
采用各种迭代法，如点逐次超松驰方法(SOR)、线逐次超松驰方法(LSOR)、迭代的交替方向隐式方法(IADI)及强隐式方法(SID)等。
(1) 现在分别对时间（从0时刻到到期日）和股票价格（S max )为可达到的足够高的股票价格）进行分割，即\triangle S=S_{max}/M,\triangle T/N,这样就分别有N+1个时间段和M+1个股票价格，建立如图(所示的坐标方格，将定解区域网格化，坐标方格上的点（i,j ）对应时刻和股票价格，用变量f i ,j 表示（i,j ）点的期权价格。
2.建立差分格式（1）内含的有限差分方法其步骤可分为以下几步：(1)求前向差分近似:(2) 后向差分格式：(3)将(2),(3)式平均可更加对称地求出的近似，即(4)(2)求用前向差分近似：(5)(3)求(6)(4)将(4),(5),(6)式代入(1)式可得到内含有限差分公式：+ b j f i,j−c j f i,j + 1 = f i + 1,j(7)aj f i,j− 1其中：i=0,1,…，N-1。
j=0,1…，M-1针对看跌期权和看涨期权可分别求出方程的边界条件：看跌期权：看涨期权：(5)利用边界条件和(7)式可以给出M-1个联立方程组：+ b j f N− 1,j + c j f N− 1,j + 1j=1,2…，M-1aj f N− 1,j− 1求解这M-1个联立方程组即可以求出期权价格，但对美式看跌期权时我们必须考虑其提前执行的情况。
1. 引言有限差分法（Finite Difference Method，FDM）是一种求解微分方程数值解的近似方法，其主要原理是对微分方程中的微分项进行直接差分近似，从而将微分方程转化为代数方程组求解。
有限差分法的原理简单，粗暴有效，最早由远古数学大神欧拉（L. Euler 1707-1783）提出，他在1768年给出了一维问题的差分格式。
1908年，龙格（C. Runge 1856-1927）将差分法扩展到了二维问题【对，就是龙格-库塔法中的那个龙格】。
但是在那个年代，将微分方程的求解转化为大量代数方程组的求解无疑是将一个难题转化为另一个难题，因此并未得到大量的应用。
随着计算机技术的发展，快速准确地求解庞大的代数方程组成为可能，因此逐渐得到大量的应用。
发展至今，有限差分法已成为一个重要的数值求解方法，在工程领域有着广泛的应用背景。
本文将从有限差分法的原理、基本差分公式、误差估计等方面进行概述，给出其基本的应用方法，对于一些深入的问题不做讨论。
2. 有限差分方法概述首先，有限差分法是一种求解微分方程的数值方法，其面对的对象是微分方程，包括常微分方程和偏微分方程。
此外，有限差分法需要对微分进行近似，这里的近似采取的是离散近似，使用某一点周围点的函数值近似表示该点的微分。
下面将对该方法进行概述。
2.1. 有限差分法的基本原理这里我们使用一个简单的例子来简述有限差分法的基本原理，考虑如下常微分方程\begin{cases} u'(x)+c(x)u(x)=f(x), \quad x \in [a, b]; \\u(x=a) = d \end{cases} \tag{1}微分方程与代数方程最大的不同就是其包含微分项，这也是求解微分方程最难处理的地方。
有限差分法的基本原理即使用近似方法处理微分方程中的微分项。
为了得到微分的近似，我们最容易想到的即导数定义u'(x)=\lim_{\Delta x\rightarrow 0} \frac{u(x+\Delta x)-u(x)}{\Delta x}\approx \frac{u(x+\Delta x)-u(x)}{\Delta x} \tag{2}上式后面的近似表示使用割线斜率近似替代切线斜率，\Delta x 即为步长，如图 1(a)所示。
第四章有限差分方法4.1引言有限差分法：数值求解常微分方程或偏微分方程的方法。
物理学和其他学科领域的许多问题在被分析研究之后, 往往可以归结为常微分方程或偏微分方程的求解问题。
一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。
有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。
在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值。
但是从原则上说，这种方法仍然可以达到任意满意的计算精度。
因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值插值计算来近似得到。
这种方法是随着计算机的诞生和应用而发展起来的。
其计算格式和程序的设计都比较直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分。
有限差分法的具体操作分为两个部分：(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式； (2)求解差分方程组。
在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。
通常采用的是规则的分割方式。
这样可以便于计算机自动实现和减少计算的复杂性。
网络线划分的交点称为节点。
若与某个节点P 相邻的节点都是定义在场域内的节点，则P 点称为正则节点；反之，若节点P 有处在定义域外的相邻节点，则P 点称为非正则节点。
在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。
有限差分法的差分格式：一个函数在x 点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。
如对一个单变量函数f(x)，x 为定义在区间[a,b]的连续变量。
以步长h=Δx 将[a,b]区间离散化，我们得到一系列节点x = a , x = x + h , x = x + h = a + 212132Δx , ..., x = x + h = b , 然后求出 f(x)在这些点上的近似值。
有限差分法求解SchrÖdinger方程摘要: 针对量子力学中大多数量子体系的哈密顿算符都比较复杂, 薛定谔方程均不能得到严格解或分析解的问题,提出了用数学中的有限差分法来解决计算量子力学中薛定额方程的本征问题.对普通的径向薛定谔方程和含时的薛定谔方程进行了有限差分法的分析, 给出了两种薛定谔方程的有限差分法的离散方程, 并以线性谐振子为例, 进行了计算机编程推算. 结果表明, 该方法在研究量子力学问题中具有广泛的应用前景.关键词: 有限差分法; 本征值; 波函数引言在量子力学中,对于一些简单的量子体系,如一维无限深势阱、线性谐振子、氢原子等体系的薛定谔(SchrÖdinger)方程可以严格求解, 得到描述体系的精确的状态波函数和能量. 但大多数量子体系的哈密顿算符(Hamiltonian)都比较复杂, 薛定谔(SchrÖdinger)方程均不能得到严格解或解析解, 有时能级可以给出解析表达式, 但却无法得到波函数的解析表达式. 因此,研究和发展薛定谔(SchrÖdinger)方程的计算方法就具有重要的意义.由于有限差分法可以处理几乎所有形式的势能函数, 并且编制的计算机主程序不依赖于势能函数的具体形式. 因此可以进行相对精确的计算, 对于确定量子体系的能级和相应的波函数有着重要的意义.第一章 有限差分法1 一维薛定谔方程的有限差分解法根据有限差分法中的二阶微分中心差分算符(其中忽略(3x ∆)及更高阶项)()()()()()22212d f x f x x f x f x x dx x ∆∆∆=+-+-⎡⎤⎣⎦ (1.1) 可以将一维薛定谔方程()()()()2222d x V x x E x dx ψψψμ-+= (1.2) 化为()()()()()()2222x x x x x x V x E x μψ∆ψψ∆∆ψ+-+-=-⎡⎤⎣⎦ (1.3)以点()0,1,2,3,,m x a m x m M ∆=+=, b ax M∆-=(其中,a b 为左右边界点,计算时将此点波函数设为零)将坐标分为M 个相等的间隔, 当M 充分大时, x ∆就足够小. 将第m 个分点的波函数简记为()m m x ψψ=.这样, (1.3)式可以化简为()21122m m m m m x E μψβψψ∆ψ+--+-=(1.4)式中()()2222m m x V x μβ∆=+令()()22,22m m R R V x x αμ∆==+(1.3)可改写为11m m m m m R R E ψαψψψ-+-+-= (1.5)2 边界条件与计算区间当取0,1,2,3,,.m M = 并且注意到满足条件00M ψψ==, 则由(5)式得到一系列线性方程式 11211223223343322122111M M M M M M M M M R E R R E R R E R R E R E αψψψψαψψψψαψψψψαψψψψαψψ----------=-+-=-+-==-+-=-+=这样将本征值方程离散化为矩阵方程S E ψψ= (1.6)其中1122332211000000000000,000000000M M M M R R R R R S R R Rαψαψαψψαψαψ-----⎛⎫⎛⎫⎪ ⎪-- ⎪ ⎪ ⎪ ⎪--==⎪ ⎪ ⎪ ⎪ ⎪ ⎪-- ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭ (1.7) 我们将相对复杂的方程就转化为矩阵S 的对角化问题, 利用Matlab 可以容易将矩阵S 的本征值和本征函数同时求出.薛定谔方程一维束缚态解的特征是当x →±∞时, 波函数0ψ→. 实际情况是在有限远处波函数已经为零, 由其对低能级的状态而言. 当能量逐渐增大, 波函数不为零的区域也在逐渐扩大. 针对某一能级计算时区间的选取要满足大于波函数不为零的区域, 也就是说, 在未达到计算区间处波函数应该已经为零, 这才能保证满足束缚态要求的边界条件. 为保证这点, 当计算得到的波函数仅仅是在计算区间的起始和终点为零而不是已经为零, 就要进一步扩大计算区间再进行计算, 一直到在未到计算的始终点位置波函数已经明显为零. 并且这时候增大或减少计算区间(保证计算区间包含了波函数不为零的区域)的长度不再影响本征值和波函数不为零区域的形态. 条件00M ψψ==相当于在计算起点和终点处存在一无限高势垒, 只有当它足够远时才不会影响所计算粒子的所处的势能及粒子的运动状态. 这里的讨论不包括计算无限深势阱中运动的粒子, 因为这种情况计算的起点和终点处实际上就存在无限高势垒.3 波函数的归一化由有限差分法计算所得的波函数并不是归一化的波函数, 为了求得归一化波函数, 可以利用数值积分求得归一化系数, 然后除以波函数. 当然数值积分不可能像解析解一样积分区间取为(),-∞∞, 但因为在计算区间之外波函数全为零, 所以积分区域可以选择为()0,M x x . 这样经过归一化后的波函数的曲线与解析解完全一致.第二章 简谐振子的有限差分解2.1 算法我们以线性谐振子为例, 线性谐振子的能量本征值方程为22222022d E x dx ψμωψμ⎛⎫+-= ⎪⎝⎭(2.1)为方便起见,引入没有量纲的变量ξ代替x ,它们的关系是,x x μωξαα≡≡=(2.2)并令2Eλω=(2.3)薛定谔方程可改写为下面无量纲的形式222d d ψξψλψξ-+= (2.4) 令(),,1m m m m a m m M ψψξξ∆ξ==+≤≤, 考虑到(1.3)式和边界条件(),0x x ψ→∞→ (2.5)(2.4)式可写为()()()211222121m m m m m ψξψψλψ∆ξ∆ξ∆ξ-+⎛⎫ ⎪-++-= ⎪⎝⎭(2.6) 方程(2.6)式, 实际上可以将(1.5)式两边同乘2ω得到. 因为()()()()()()2222222222211,2222222222mm mmR x x RV x x Eμωωωμ∆∆ξ∆μωαξωωωω∆ξ∆ξλω====+=+=+= 边界条件(2.5)可写为00,0M ψψ== (2.7)考虑到(2.7), 式(2.6)可写为S ψλψ=其中()()()()()()()()()()()()()212222222122322232221222212221000001210000121000,12100001200M m M m ξ∆ξ∆ξξψ∆ξ∆ξ∆ξψξψ∆ξ∆ξ∆ξψψξψ∆ξ∆ξ∆ξξ∆ξ∆ξ----⎛⎫+-⎪ ⎪ ⎪ ⎪-+-⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪-+-⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-+-⎝⎭ ⎪⎪⎪-+ ⎪ ⎪⎝⎭(2.8)2.2 计算区间的选择如果我们取计算区域为[]2,2-, 这相当于粒子在图2.1所示的势场内运动. 图1程序: clear x1=-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([-4 4 0 5]) 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') 而在此势场中运动的粒子与波函数可由如下程序计算得到: clear a=-2; b=2; N=1000; M=N+1; Dx=(b-a)/N; for m=1:Mx(m)=a+(m-1)*Dx;图 2.1 计算区间取为[]2,2-时, 粒子实际所处的势.end F=zeros(M); for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2; endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end[v,d]=eig(F); for n=1:1E(n)=d(n,n)/2 %第n 个本征值 figureplot(x,v(:,n),'-k', 'LineWidth',2) %第n 个本征函数 axis([-2 2 0 0.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)'];grid onxlabel ('\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-, 再重新计算其程序修改为:clear a=-4; b=4; N=2000; M=N+1; Dx=(b-a)/N; for m=1:Mx(m)=a+(m-1)*Dx; end图2 计算区间取为[]2,2-时, 粒子波函数的形状.F=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E(n)=d(n,n)/2 %第n个本征值figureplot(x,v(:,n),'-k', 'LineWidth',2) %第n个本征函数axis([-4 4 0 0.05])set(gca,'XTick',-4:1:4)set(gca,'YTick',0:0.01:0.05)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];grid onxlabel ('\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;for m=1:M 图3 计算区间取为[]4,4-时, 粒子波函数的形状.x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2; endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end[v,d]=eig(F); for n=1:1E(n)=d(n,n)/2 %第n 个本征值figureplot(x,v(:,n),'-k', 'LineWidth',2) %第n 个本征函数axis([-5 5 0 0.05]) set(gca,'XTick',-5:1:5) set(gca,'YTick',0:0.01:0.05)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];grid onxlabel ('\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基本一致. 如果把两张图绘制在一起, 发现它们在区间[]4,4-内是完全重合的(条件是计算的精度相同, 即每单位长度的计算点数相同).选取计算精度分别为1000、2000、3000、4000、5000、6000计算本征值与绘制本征函数, 程序如下: clear a=-4; b=4; N=1000; M=N+1; Dx=(b-a)/N;图5计算区间分别取为[]4,4-和[]5,5-时,粒子波函数在[]4,4-区间由点线绘制,在[]5,5-区间由实线绘制,其中在区间[]4,4-内是完全重合.for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E1(n)=d(n,n)/2figureplot(x,v(:,n),'-k', 'LineWidth',2)axis([-4 4 0 0.07])set(gca,'XTick',-4:1:4)set(gca,'YTick',0:0.01:0.07)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)']; grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k')ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')endcleara=-4;b=4;N=2000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E2(n)=d(n,n)/2hold onplot(x,v(:,n),'-k', 'LineWidth',2) endcleara=-4;b=4;N=3000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E3(n)=d(n,n)/2hold onplot(x,v(:,n),'-k', 'LineWidth',2) endcleara=-4;b=4;N=4000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E4(n)=d(n,n)/2hold onplot(x,v(:,n),'-k', 'LineWidth',2) endcleara=-4;b=4;N=5000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E5(n)=d(n,n)/2hold onplot(x,v(:,n),'-k', 'LineWidth',2)endcleara=-4;b=4;N=6000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E6(n)=d(n,n)/2hold onplot(x,v(:,n),'-k', 'LineWidth',2)end图见图6+其中从上到下分别为精度为1000、2000、3000、4000、5000、6000对应的本征函数. 对于相同计算区间[]4,4-, 采用不同计算精度分别为1000、2000、3000、4000、5000、6000点得到的本征值为:E1 =0.49999846150440E2 =0.49999997592365E3 =0.50000025861460 E4 =图6计算区间取为[]4,4-, 计算点数分别为1000、2000、3000、4000、5000和6000.对应的波函数曲线在图中从上到下分布.0.50000035829738E5 =0.50000040479853E6 =0.50000043022122由图6分析可知, 提高其计算精度, 达到一定的程度后对本征值影响不大, 但波函数的变化仍然明显, 不过逐渐趋于相同, 差别在缩小.2.3波函数的归一化我们选择计算区间为[]-, 计算点分别取1000、2000、3000、4000、5000、6000进行计算, 波函数4,4利用语句A=sqrt(trapz(x,v(:,n).*v(:,n)))来进行归一, 计算程序为:cleara=-4;b=4;N=1000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E1(n)=d(n,n)/2figureA=sqrt(trapz(x,v(:,n).*v(:,n)));plot(x,v(:,n)/A,'-k', 'LineWidth',2)axis([-4 4 0 0.8])set(gca,'XTick',-4:1:4)set(gca,'YTick',0:0.1:0.8)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k') ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-') set(gca,'xgrid','on','GridLineStyle','-') endcleara=-4;b=4;N=2000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E2(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n)));hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2)endcleara=-4;b=4;N=3000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E3(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) endcleara=-4;b=4;N=4000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E4(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) endcleara=-4;b=4;N=5000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E5(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) endcleara=-4;b=4;N=6000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E6(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n)));hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2)end计算得到的本征值为:E1 =0.5000, E2 =0.5000, E3 =0.5000, E4 =0.5000, E5 =0.5000, E6 =0.5000计算所得的波函数如图7所示, 所有波函数的曲线均重合为一条曲线.下面我们设计的程序计算区间分别取[][][][][]-, 计算点分别取9,9-----和[]4,4,5,5,6,6,7,7,8,82000、2500、3000、3500、4000和4500, 这样可保证计算精度相同. 计算程序:cleara=-4;b=4;N=2000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E1(n)=d(n,n)/2figureA=sqrt(trapz(x,v(:,n).*v(:,n)));plot(x,v(:,n)/A,'-k', 'LineWidth',2)axis([-9 9 0 0.8])set(gca,'XTick',-9:1:9)set(gca,'YTick',0:0.1:0.8)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)']; grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k')ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')endcleara=-5;b=5;N=2500;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E2(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n)));hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2)endcleara=-6;b=6;N=3000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E3(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) endcleara=-7;b=7;N=3500;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E4(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) endcleara=-8;b=8;N=4000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E5(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) endcleara=-9;b=9;N=4500;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2; endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end[v,d]=eig(F); for n=1:1 E6(n)=d(n,n)/2A=sqrt(trapz(x,v(:,n).*v(:,n))); hold onplot(x,v(:,n)/A,'-k', 'LineWidth',2) end计算结果本征值为: E1 =0.5000, E2 =0.5000, E3 =0.5000, E4 =0.5000, E5 =0.5000, E6 =0.5000 计算所得的波函数如图8所示, 所有波函数在相同的区间内是重合的. 在区间[][][]4,4,5,5,6,6,---[][]7,7,8,8--和[]9,9-分别有6、5、4、3、2条波函数的曲线均重合为一条曲线.数值计算所得的波函数是不归一的, 所以和解析解相差一归一化因子. 这样数值计算所得的波函数图形与解析解的波函数图形不一致. 为了由数值计算所得的波函数和解析解一样, 我们可以将数值计算所得的波函数归一化, 方法是对数值计算所得的波函数的积进行数值积分, 然后波函数除此积分结果的开平方值即可. 积分区域为计算所进行的区域即可, 因为在此区域外波函数为零, 这样得到的归一化因子除以波函数即可得到与解析解一致的波函数. 求归一化因子的程序是A=sqrt(trapz(x,v(:,n).*v(:,n))), 这里的n 是指求出的归一化因子是属于第n 个本征值的本征波函数的归一化因子. 当原来数值计算得到的波函数除以相应的归一化因子后, 所得到的波函数与解析解完全一致.2.4 数值解与解析解对照我们知道简谐振子的波函数的解析解为:()()()222221122221122112!2!nn n n n n n n n d d e e e e e d d n n ξξξξξααψξξξππ---⎛⎫⎛⎫ ⎪ ⎪=-=- ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭. 由此可得解析解绘图程序: clear n=1;图7计算区间取为[][][]4,4,5,5,6,6,---[][][]7,7,8,8,9,9---,计算点数分别为2000、2500、3000、3500、4000和4500.所有的波函数曲线在相同的区间重合为一条曲线.syms xy=sym('exp(-x^2)');y1=sym('exp(x^2/2)');dydxn=diff(y,'x',n);Hn=(-1)^n*dydxn;Psix=(1/(2^n*pi^0.5*factorial(n))^(1/2))*y1*Hn;X=-5:0.01:5;Y1=subs(Psix,x,X);hold onplot(X,Y1,'-r','LineWidth', 2)%axis([-5 5 -1 1])%set(gca,'XTick',-5:1:5)%set(gca,'YTick',-1:0.2:1)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)']; grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k')ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')数值计算程序:clearn=1a=-5;b=5;N=1000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);E=d(n,n)/2A1=sqrt(trapz(x,v(:,n).*v(:,n)))figureplot(x,v(:,n)/A1,'-k', 'LineWidth',2)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)']; hold onxlabel ('\fontsize{14}\rm\xi', 'Color','k')ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')将两张图画在一起发现它们完全重合如图8所示.%绘制第n个本征函数并给出第n个本征值程序:clearn=1a=-5;b=5;N=2000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:4d(n,n) %第n个本征值A=sqrt(trapz(x,v(:,n).*v(:,n))); %归一化系数figureplot(x,v(:,n)/A,'-k', 'LineWidth',2) %第n个本征函数axis([-5 5 -1 1])set(gca,'XTick',-5:1:5)set(gca,'YTick',-1:0.2:1)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k')ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')end将波函数归一化后, 发现在达到一定的计算精度后, 不论计算区域变化或计算精度变化, 在计算区域重合的区域, 波函数图像是完全重合的. 图7中的实线、点线和虚线分别是计算区间取[]-计算精度5,5为2000和3000及计算区间取[]-, 计算精度取1600. 所得的本征值为4,4E1 =0.49999921883387E2 =0.49999965283684E3 =0.49999969105011而波函数如图8所示, 三者完全重合.计算程序:cleara=-5;b=5;N=2000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E1(n)=d(n,n)/2 %第n个本征值A=sqrt(trapz(x,v(:,n).*v(:,n))); %归一化系数figureplot(x,v(:,n)/A,'-k', 'LineWidth',2) %第n个本征函数axis([-5 5 0 0.8])set(gca,'XTick',-5:1:5)set(gca,'YTick',0:0.1:0.8)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)'];grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k')ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-')set(gca,'xgrid','on','GridLineStyle','-')endcleara=-5;b=5;N=3000;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E2(n)=d(n,n)/2 %第n个本征值A=sqrt(trapz(x,v(:,n).*v(:,n))); %归一化系数hold onplot(x,v(:,n)/A,':k', 'LineWidth',2) %第n个本征函数endcleara=-4;b=4;N=1600;M=N+1;Dx=(b-a)/N;for m=1:Mx(m)=a+(m-1)*Dx;endF=zeros(M);for m=1:M;F(m,m)=2/((Dx)^2)+x(m)^2;endfor m=1:M-1;F(m+1,m)=-1/((Dx)^2);F(m,m+1)=-1/((Dx)^2);end[v,d]=eig(F);for n=1:1E3(n)=d(n,n)/2 %第n个本征值A=sqrt(trapz(x,v(:,n).*v(:,n))); %归一化系数hold onplot(x,v(:,n)/A,':k', 'LineWidth',2) %第n个本征函数end简谐振子的解析解为:()()()222221122221122112!2!nn n n n n n n n d d e e e e e d d n n ξξξξξααψξξξππ---⎛⎫⎛⎫ ⎪ ⎪=-=- ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭. 由此可得解析解绘图程序: clear n=1; syms xy=sym('exp(-x^2)'); y1=sym('exp(x^2/2)'); dydxn=diff(y,'x',n); Hn=(-1)^n*dydxn;Psix=(1/(2^n*pi^0.5*factorial(n))^(1/2))*y1*Hn; X=-5:0.01:5;Y1=subs(Psix,x,X); hold onplot(X,Y1,'-r','LineWidth', 2) %axis([-5 5 -1 1])%set(gca,'XTick',-5:1:5) %set(gca,'YTick',-1:0.2:1)t=['\fontsize{14}\rm\Psi','_{',int2str(n-1),'}','(\xi)']; grid onxlabel ('\fontsize{14}\rm\xi', 'Color','k') ylabel (t, 'Color','k')set(gca,'ygrid','on','GridLineStyle','-') set(gca,'xgrid','on','GridLineStyle','-')数值计算所得的波函数是不归一的, 所以和解析解相差一归一化因子. 这样数值计算所得的波函数图形与解析解的波函数图形不一致. 为了由数值计算所得的波函数和解析解一样, 我们可以将数值计算所得的波函数归一化, 方法是对数值计算所得的波函数的积进行数值积分, 然后波函数除此积分结果的开平方值即可. 积分区域为计算所进行的区域即可, 因为在此区域外波函数为零, 这样得到的归一化因子除以波函数即可得到与解析解一致的波函数. 求归一化因子的程序是A=sqrt(trapz(x,v(:,n).*v(:,n))), 这里的n 是指求出的归一化因子是属于第n 个本征值的本征波函数的归一化因子. 当原来数值计算得到的波函数除以相应的归一化因子后, 所得到的波函数与解析解完全一致.参考文献:[1] 刘建军，翟利学. 有限差分法解能量本征方程[J]. 北京工业大学学报, 2008, 34(3): 325-331.[2] 陈皓，高明，汪青杰. 用有限差分法解薛定谬方程[J]. 沈阳航空工业学院学报, 2005, 22(l):87-88.[3]RlVA C, ESCORCI R A, GOVOROV A O, et al.Charged donors in quantum dots: Finite differenceand fraetional dimensions results[J]. Phys. Rev.B. 2004, 69(24): 245306·245313.[4]AN xing-tao, LIU Jian-jun. Hydrogenic impurities in parabolic quantutn-well wires in amagnetic field[J]. Appl. phys. 2006, 99(12): 123713-123717.[5]ZHAI Li-xue, LIU Jian-jun. Negatively charged donors in paralic quantum-wellwires undermagnetic fields[J]. J.Appl. Phys., 2007, 102(5): 053706-053710.[6] RlVAC, SCHWEIGERT V A., PEETERS F M. Off-center D-centers in a quantum well in the presenceof a perpendieular magnetic field: Angular-momentum transitions and magnetic evaporation.Phs, Rev.B, 1998, 57(24):15392-15399.[7] JOLY Y. Finite-difference method for the caleulation of low-energy, positron diffraction[J].Phys. Rev. B, 1996.53(19): 13029.13037.[8] 曾谨言. 量子力学[M]. 北京: 科学出版社, 2000.[9] 孙志忠. 偏徽分方程数值解法[M]. 北京: 科学出版社, 2005:1-2.（注：可编辑下载，若有不当之处，请指正，谢谢!）。
