二维抛物方程的有限差分法

《二维抛物方程的有限差分法》由会员分享，可在线阅读，更多相关《二维抛物方程的有限差分法（34页珍藏版）》请在装配图网上搜索。

1、华北电力大学本科毕业设计（论文）二维抛物方程的有限差分法摘要二维抛物方程是一类有广泛应用的偏微分方程，由于大部分抛物方程都难以求得解析解，故考虑采用数值方法求解。有限差分法是最简单又极为重要的解微分方程的数值方法。本文介绍了二维抛物方程的有限差分法。首先，简单介绍了抛物方程的应用背景，解抛物方程的常见数值方法，有限差分法的产生背景和发展应用。讨论了抛物方程的有限差分法建立的基础，并介绍了有限差分方法的收敛性和稳定性。其次，介绍了几种常用的差分格式，有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等，重点介绍了古典显式

2、格式和交替方向隐式格式。进行了格式的推导，分析了格式的收敛性、稳定性。并以热传导方程为数值算例，运用差分方法求解。通过数值算例，得出古典显式格式计算起来较简单，但稳定性条件较苛刻；而交替方向隐式格式无条件稳定。关键词：二维抛物方程；有限差分法；古典显式格式；交替方向隐式格式FINITE DIFFERENCE METHOD FOR TWO-DIMENSIONAL PARABOLIC EQUATION Abstract Two-dimensional parabolic equation is a widely used class of partial differential equation

3、s. Because this kind of equation is so complex, we consider numerical methods instead of obtaining analytical solutions. finite difference method is the most simple and extremely important numerical methods for differential equations. The paper introduces the finite difference method for two-dimensi

4、onal parabolic equation.Firstly, this paper introduces the background and common numerical methods for Parabolic Equation, Background and development of applications. Discusses the basement for the establishment of the finite difference method for parabolic equation And describes the convergence and

5、 stability for finite difference method.Secondly, Introduces some of the more common simple differential format,for example, the classical explicit scheme, the classical implicit scheme, Crank-Nicolson implicit scheme, Douglas difference scheme, weighted six implicit scheme and the alternating direc

6、tion implicit format. The paper focuses on the classical explicit scheme and the alternating direction implicit format. The paper takes discusses the derivation convergence,and stability of the format . The paper takes And the heat conduction equation for the numerical example, using the differentia

7、l method to solve. Through numerical examples, the classical explicit scheme is relatively simple for calculation, with more stringent stability conditions; and alternating direction implicit scheme is unconditionally stable.Keywords: Two-dimensional Parabolic Equation; Finite-Difference Method; Ecl

8、assical Explicit Scheme; Alternating Direction Implicit SchemeII华北电力大学本科毕业设计（论文）目 录摘要IAbstractII1绪论11.1课题背景11.2发展概况21.2.1抛物型方程的常见数值解法21.2.2有限差分方法的发展21.3差分格式建立的基础31.3.1区域剖分31.3.2差商代替微商41.3.3差商代替微商格式的误差分析51.4本文主要研究内容52显式差分格式72.1常系数热传导方程的古典显式格式72.1.1古典显式格式格式的推导72.1.3古典显式格式的算法步骤83隐式差分格式103.1古典隐式格式103.2

9、Crank-Nicolson隐式格式113.3 Douglas差分格式133.4加权六点隐式格式133.5交替方向隐式格式143.5.1 Peaceman-Rachford格式143.5.2 Rachford-Mitchell格式153.5.3 Mitchell-Fairweather格式153.5.4交替方向隐式格式的算法步骤154实例分析与结果分析164.1算例164.1.1已知有精确解的热传导问题164.1.2未知精确解的热传导问题184.2结果分析195稳定性探究与分析205.1稳定性问题的提出205.2 几种分析稳定性的方法205.3 r变化对稳定性的探究215.3.1 古典显式格式

10、的稳定性215.3.2 P-R格式格式的稳定性22结语24参考文献25附录 P-R格式的C+实现代码26致谢28 华北电力大学本科毕业设计（论文）1绪论1.1课题背景抛物方程是一类特殊的偏微分方程，二维抛物方程的一般形式为 (1-1)其中。1渗流、扩散和热传导、静电场等很多领域经常会遇到求解二维抛物型方程。微分方程愈复杂，找出解的解析表达式愈困难。对大部分的抛物方程而言，找出解的解析表达式及其困难。因此求出抛物方程的近似解是很有意义的1。本文研究的近似方法是数值方法，目标在于给出解在一些离散点上的近似解值。常见的数值求解方法有有限差分法，有限元法，有限体积法等。有限差分法(Finite Dif

11、ference Methed)是在1966年Yee提出的，用于解决电机工程中电磁场的初值和边值问题2。有限差分法得到迅速发展后，不仅广泛应用于自然科学，在社会科学的各领域也在越来越多地被应用着2。某些常用的数值解法，欧拉方法，隐式欧拉法，一般单步法，Crank-Nicolson隐式格式，Runge-Kutta方法等，已成为微分方程数值解领域的经典方法.，在工程计算中的应用随处可见1；基于这些方法，人们也在做更深的研究。采用有限差分法解决二维抛物方程，一些研究者采用一维抛物型方程的古典显格式，古典隐格式，Crank-Nicolson隐式格式，加权六点隐格式等的直接推广；还有有一些研究者采用交叉方

12、向的隐式差分格式(Alternating Direction Implicit Scheme)等基于二维的方法。1.2发展概况1.2.1抛物型方程的常见数值解法抛物型方程的常见数值解法有有限差分法，有限元法，有限体积法，有限单元法，无网格方法等。有限差分法(Finite Difference Methed)的主要思想是将问题离散化，将问题离散化，用差商代替微商，将微分方程和定界问题都用代数方程来代替。然后解这些代数问题构成的方程组。该方法包括区域剖分和差商代替导数两个过程。具体地，首先将求解区域划分为差分网格，用有限个网格节点代替连续的求解区域。其次，利用Taylor级数展开等方法将偏微分方程

13、中的导数项在网格节点上用函数值的差商代替来进行离散，从而建立以网格节点上的值为未知量的代数方程组1。有限元方法(Finite Element Methods)的基础是变分原理和分片多项式插值。该方法的构造过程包括以下三个步骤。首先，利用变分原理得到偏微分方程的弱形式(利用泛函分析的知识将求解空间扩大)。其次，将计算区域划分为有限个互不重叠的单元(三角形、四边形、四面体、六面体等)。再次，在每个单元内选择合适的节点作为求解函数的插值点，将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式，得到微分方程的离散形式1。有限体积法(Finite Volume Me

14、thods)又称为控制体积法。其基本思路是：将计算区域划分为一系列互不重叠的控制体，并使每个网格点周围有一个控制体；将待求解的微分方程对每一个控制体积积分，便得出一组离散方程。该方法的未知量为网格点上的函数值。为了求出控制体积的积分，须假定函数值在网格点控制体边界上的变化规律。从积分区域的选取方法看来，有限体积法属于有限元方法中检验函数取分片常数插值的子区域法；从未知量的近似方法看来，有限体积法属于采用局部近似的离散方法1。有限单元法的基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，借助于变分原理或加权余量法，将微分方程离散求解1。无网格

15、方法的近似函数建立在一些离散的节点上，不需要借助网格，在涉及到网格畸变，网格移动等问题中显示出了明显的优势1。1.2.2有限差分方法的发展随着人们对显示格式的进一步研研究，提出了很多新的高精度，绝对稳定的差分格式。文献2讨论了变系数热传导方程的两层绝对稳定差分格式。研究了二维变系数非齐次热传导方程的两层绝对稳定的差分格式问题。首先运用Pade逼近导出了差分格式,给出了差分格式的截断误差;讨论了差分格式的绝对稳定性和收敛性,且收敛阶为;最后给出了数值例子,数值结果和理论结果是吻合的。文献3第一部分用待定系数法对维抛物型方程构造出了高精度(截断误差达)能显式计算，稳定性较好()的三层(特殊情况下可

16、以是两层)显式差分格式，在空间变量更多的情况下，指出了构造高精度差分格式的一般方法。文献3第二部分用算子方法对二维和三维抛物型方程构造出了高精度(截断误差达)的绝对稳定的交替方向隐式差分格式，在空间变量更多的情况下，也指出了构造高精度交替方向隐式差分格式的一般方法。并附有数值例子，它表明理论分析的正确性和所建立的差分格式的有效性。文献4对二阶抛物型方程构造了一族含参数高精度三层差分格式。当参数满足一定的条件时,差分格式绝对稳定，其局部截断误差阶数最高可达。适当地调节参数，可以得到一个七点显式差分格式和一个两层六点隐格式。数值例子表明,对稳定性所作的分析是正确的。文献5多维抛物型方程和双曲方程的

17、差分解法对一般m维热传导方程，波动方程的初、边值问题，提出若干个交替方向的差分格式。1.3差分格式建立的基础6 1.3.1区域剖分14构造一维线性抛物方程 (1-2) 的有限差分逼近，首先将求解区域用两组平行于轴和轴的直线构成网格覆盖，网格边长在方向为，方向为，网格节点为，为整数6。网格如图1-1。图1-1 一维抛物方程的网格剖分二维抛物方程的区域剖分将剖分区域扩展到三维空间。首先将求解区域用三组平行于轴，轴和轴的直线构成网格覆盖，网格边长在方向为，方向为，方向为，网格节点为，为整数3。网格如图1-26。图1-2 二维抛物方程的网格剖分1.3.2差商代替微商差分方法的基本思想就是以差商代替微商

18、。考虑如下两个Taylor公式： (1-3) (1-4)从式(1-3)得到： 从式(1-4)得到：从式(1-3) -式(1-4)得到：从式(1-3) +式(1-4)得到： 下表1-1面介绍常用的几种差商6。表1-1：常用的几种差商形式有限差分近似公式误差阶（向前差分） （向后差分） （中心差分）由式(1-1)可得： (1-5)1.3.3差商代替微商格式的误差分析为了分析分析数值方法的精确度，考察收敛性。常常在成立的假定下，估计误差这种误差称为“局部截断误差”局部截断误差是以点的精确解为出发值，用数值方法推进到下一个点而产生的误差8。若如下： (1-6) (1-7)局部截断误差是关于的极小函数，

19、则收敛13。为了分析分析数值方法的精确度，考察收敛性。整体截断误差是以点的初始值为出发值，用数值方法推进i+1步到点，所得的近似值与精确值 的偏差：称为整体截断误差8。若如下： 14整体截断误差不随无限扩大，则收敛。1.4本文主要研究内容本文主要研究二维抛物方程的有限差分方法，研究工作分为以下四个方面：1)差分格式的构造方法有限差分法法解二维抛物方程包括区域剖分和差商代替导数两个过程。差分格式就是在网格结点上求出微分方程的近似解的一种方法，因此又称为网格法。2)差分格式的分析方法数值方法是近似计算方法，需从收敛性，相容性，稳定性方面考察。收敛性是考察步长足够小时，数值解是否无限逼近 解析解。稳

20、定性主要最初产生的小误差在以后的计算中虽然会传递下去，但不会无限制的扩大。3)显式差分格式显式差分格式是差分方法中可逐层逐点分别求解的格式，是一种有限差分近似方法。显式差分格式的优点在于计算简单，但它并不总是稳定的。4)隐式差分格式与显式差分格式不同，隐式差分格式中包括了（n+1）时间层上二个或二个以上结点处的未知值（例如 ），使用隐式差分格式和使用显式差分格式求解完全不同。相对而言，使用隐式差分格式求解，每时间层包含有较多的计算工作量。从后面对差分格式的稳定性分析可知，隐式格式的优点在于，其稳定性要求对步长比的限制大为放宽，而这正是我们所期望的。 2显式差分格式现在，对二维抛物方程式(1-1

21、) ，从方程式(1-4)出发，构造有限差分的显式格式。由于一维热传导方程 (2-1)二维热传导方程 (2-2)在热传导，磁扩散等许多领域有重要的应用。实际导热问题必然涉及边值条件，在有限差分法中它们也必须差分化。因此，我们需要研究的不仅是差分方程本身，而且是包括全部内部区域和所有边界上的差分方程所组成的代数方程组。又称为差分格式3。 2.1常系数热传导方程的古典显式格式首先考虑热传导方程的边值问题离散化152.1.1古典显式格式格式的推导现在对热传导方程式(2-1) 推导其最简单的显隐式差分逼近古典显隐式格式。由故式中右边如果仅保留二阶导数项，且以替代，替代，则得差分格式 或者 (2-3)这是

22、一个显式格式(四点格式)，如图(2-1)所示。图2-1：古典显式格式格式(2-3)应用在一维热传导问题中，得到古典显式格式为： (2-3)每一层各个节点上的值是通过一个方程组求解得到的。显一隐格式区域分解方法就是以显格式计算出相邻子区域相交内边界的近似值的一种方法。显一隐格式区域分解方法综合了二者的优点，借助前一层数值解的信息，用显格式给出在这层的子问题的未知内边界条件，把一个整体区域上的问题化为若干个子区域上的子问题，在每个子区域上用隐式方法求解，从而实现了并行。这可以从下面的计算过程看出来8。2.1.3古典显式格式的算法步骤古典显式格式，以m=1为例，列成矩阵有如下形式： 系数矩阵为为了得

23、到二维问题的误差估计，我们做一些改动。因此，这个算法从到的计算步骤为：第一步，根据给出的初边值条件，得出t=0时；第二步，根据古典显式格式的公式，由第得出6。下文中将通过一个具体的数值例子说明了计算的方法，体现了这种格式的实用性和优越性。3隐式差分格式与显式差分格式不同，隐式差分格式中包括了(n+1)时间层上二个或二个以上结点处的未知值(例如 )，使用隐式差分格式和使用显式差分格式求解完全不同。相对而言，使用隐式差分格式求解，每时间层包含有较多的计算工作量。从后面对差分格式的稳定性分析可知，隐式格式的优点在于，其稳定性要求对步长比的限制大为放宽，而这正是我们所期望的6。对二维抛物方程式(1-1

24、) ，从方程式(1-4)出发，构造有限差分的显式格式。由于热传导方程式(2-1)构造差分格式。我们需要研究的不仅是差分方程本身，而且是包括全部内部区域和所有边界上的差分方程所组成的代数方程组。又称为差分格式。3.1古典隐式格式现在对热传导方程推导其最简单的隐式差分逼近古典隐式格式。由故式中右边如果仅保留二阶导数项，且以替代，且以替代，替代，则得差分格式或者 (3-1)将格式(2-3)应用于一维热传导方程，古典显式格式为： (3-2)显式与隐式的比较如下6：(1)同一阶数下，隐式的局部截断误差的系数的绝对值比显式的要小；(2)显式的计算工作量比隐式的小；(3)隐式的稳定范围比显式的大。局部截断误

25、差的阶最高是3，式(3-2) 是允许函数任意变化情况下截断误差最小的二阶方法。要再提高阶就必须增加计算函数值的次数。上述式(3-2)又称为古典隐式格式。故格式用图3-1表示，其截断误差阶为,与古典差分格式相同。这是一个四点差分格式，如图3-1所示。图3-1：古典隐式格式为了求得第(n+1)时间层上的的值，必须通过解线性代数方程组。这是一个隐式差分格式，必须联合其初边值条件求解。格式(3-2)通常称为古典隐式格式。我们也可以通过直接用差分算子代替 的方即：代入微分方程，得到格式(3-1)。古典隐式格式的方程组和矩阵形式如下：当知道第n层上的时，要确定第n+1层上各点值必须通过求解一个线性代数方程

26、组。以m=1为例：其系数矩阵如下：3.2 Crank-Nicolson隐式格式Crank-Nicolson隐式差分格式是解热传导方程(2-1)的常用的差分格式，近年来，关于抛物方程的区域分解算法作为并行计算的一种有效工具，吸引了很多学者的注意。这类算法的主要困难在于；如何定义内边界点的值和在子区域上选取合理的计算解去近似。为了推导它，由式（1-4），有由得 (3-3)两边仅保留前二项，用代替，替代，则得差分格式这是一个隐式差分格式，称为Crank-Nicolson差分格式，截断误差阶为，也可写为 (3-4)由于格式（3-4）中包括六个结点，故也可称为六点格式（如图3-2所示）。图3-2： Cr

27、ank-Nicolson隐式格式也可将 代入微分方程(2-1)，得到Crank-Nicolson格式。解微分方程(2-1)，根据Crank-Nicolson格式得到的方程组：其矩阵表达式为：3.3 Douglas差分格式6基于如同Crank-Nicolson格式一样的六个网格结点可获得另一精度较高的差分格式，如在前式(3-3)中仅保留直到，的项，即有：可令则可得：代入上式，则有如下差分格式 (3-5)它称为Douglas差分格式，具有截断误差阶。精度较高的差分格式(Douglas格式),并通过一个具体的数值例子说明了计算的方法,体现了这种格式的实用性和优越性。3.4加权六点隐式格式9前面，我们

28、已经推导了热传导方程(2-1)的古典显式显示格式，古典隐式格式及Crank-Nicolson格式等。实际上，它们都可以作为本节推导的加权六点隐式格式的特殊情形。由得即两边去掉高于二阶导数的项，且用代替，替代，则得差分格式或者 (3-5)这是一个六点差分格式（如图3-3所示），称为加权六点差分格式。图3-3：加权六点差分格式3.5交替方向隐式格式交替方向隐式格式是关于x方向的显格式，关于y方向的隐格式。在n+1/2层上用古典显式格式计算出过度值，再在n+1层上用古典隐格式校正预测值。得到的稳定性和收敛性是相近的。但是，只讨论了关于常系数热传导方程的在z和可方向都是相同的时间步长的情形，所以我们在

29、这里试图做一些扩展。我们首先考虑在内边界点用一个方向为显格式，另一个方向为隐格式的情形，再考察在z和y方向都用显格式的情形，而且我们讨论在z和可方向用不同的时间步长的情形6。3.5.1 Peaceman-Rachford格式考虑二维抛物方程(2-1)的差分解法，以上对显式格式的稳定性分析发现，的要求比一维情形更苛刻。分析表明，维数越高时，要求时间步长越小，计算工作量越大3。(P-R)ADI格式是由Peaceman-Rachford在1995年发表的，从第n层到n+1层，P-R格式分两步进行，每一步只需解具有三对角系数的线性方程组6。P-R格式为 (3-6) (3-7)，是中心差算子，P-R格式

30、对任意无条件稳定。但预测P-R格式对三个变量空间问题，无条件稳定性不再成立3。3.5.2 Rachford-Mitchell格式Douglas和Rachford在1956年提出另一个隐式差分格式，即Douglas-Rachford格式6。D-R格式是第一个能被推广到三维情形的交替方向隐式格式，二维D-R格式是 (3-7) (3-8) 3.5.3 Mitchell-Fairweather格式Mitchell和Fairweather在1964年推导了一个高精度ADI差分格式，称为M-F格式6。M-F格式截断误差达较之P-R格式和D-R格式更好，且无条件收敛3。3.5.4 交替方向隐式格式的算法步骤

31、 ADI格式的算法步骤与古典显式格式的算法步骤相似，只是ADI格式是隐式格式，需用追赶法求解。追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。第一步，根据给出的初边值条件，得出t=0时;第二步，根据各个交替方向隐式格式中的第一个公式，由第求得出如式(3-6)，采用追赶法，确定追赶因子a=1+r/2,b= r/2,c=1+r/2,d=,=b-a*,=(d-*a)/w,然后根据追赶因子确定。=,= d -*;第三步，再由追赶法由第求得出6。由该算法，计算的近似解的过程。下文中将通过一个具体的数值例子说明了计算的方法,体现了这种格式的实用性和优越性。4实例分析与结果分析4.1

32、算例4.1.1已知有精确解的热传导问题例1是二维Dirichlet边值问题的热传导方程方程的精确解为采用古典显示格式，以求解得，时，关于，的函数如图4-1所式。图4-1：古显格式解例1，关于，的图像（，）分别以古典显示格式，P-R(ADI)格式，以相同的步长(，)求解同一处()的与解析解间的差值，列在误差表中如表4-1，表4-1所示。表4-1 古显格式解例1的误差（，）误差y=0y=0.1y=0.2y=0.3y=0.4y=0.5x=0000000x=0.10-0.00013-0.00025-0.00035-0.00041-0.00043x=0.20-0.00025-0.00048-0.0006

33、6-0.00077-0.00081x=0.30-0.00035-0.00066-0.0009-0.00106-0.00112x=0.40-0.00041-0.00077-0.00106-0.00125-0.00131x=0.50-0.00043-0.00081-0.00112-0.00131-0.00138x=0.60-0.00041-0.00077-0.00106-0.00125-0.00131x=0.70-0.00035-0.00066-0.0009-0.00106-0.00112x=0.80-0.00025-0.00048-0.00066-0.00077-0.00081x=0.90-0.

34、00013-0.00025-0.00035-0.00041-0.00043x=1000000续表4-1 古显格式解例1的误差（，）误差y=0.6y=0.7y=0.8y=0.9y=1x=000000x=0.1-0.00041-0.00035-0.00025-0.000130x=0.2-0.00077-0.00066-0.00048-0.000250x=0.3-0.00106-0.0009-0.00066-0.000350x=0.4-0.00125-0.00106-0.00077-0.000410x=0.5-0.00131-0.00112-0.00081-0.000430x=0.6-0.00125

35、-0.00106-0.00077-0.000410x=0.7-0.00106-0.0009-0.00066-0.000350x=0.8-0.00077-0.00066-0.00048-0.000250x=0.9-0.00041-0.00035-0.00025-0.000130x=100000表4-2 P-R(ADI)格式解例1的误差（，）误差y=0y=0.1y=0.2y=0.3y=0.4y=0.5x=0000000x=0.108.86E-071.19E-061.47E-066.44E-071.18E-06x=0.201.19E-062.05E-062.19E-069.66E-079.93E-0

36、7x=0.301.47E-062.19E-061.48E-069.21E-073.37E-09x=0.406.44E-079.66E-079.21E-07-8.6E-07-1.4E-06x=0.501.18E-069.93E-073.37E-09-1.4E-06-1.4E-06x=0.606.44E-079.66E-079.21E-07-8.6E-07-1.4E-06x=0.701.47E-062.19E-061.48E-069.21E-073.37E-09x=0.801.19E-062.05E-062.19E-069.66E-079.93E-07x=0.908.86E-071.19E-061

37、.47E-066.44E-071.18E-06x=1000000续表4-2 P-R(ADI)格式解例1的误差（，）误差y=0.6y=0.7y=0.8y=0.9y=1x=000000x=0.16.44E-071.47E-061.19E-068.86E-070x=0.29.66E-072.19E-062.05E-061.19E-060x=0.39.21E-071.48E-062.19E-061.47E-060x=0.4-8.6E-079.21E-079.66E-076.44E-070x=0.5-1.4E-063.37E-099.93E-071.18E-060x=0.6-8.6E-079.21E-0

38、79.66E-076.44E-070x=0.79.21E-071.48E-062.19E-061.47E-060x=0.89.66E-072.19E-062.05E-061.19E-060x=0.96.44E-071.47E-061.19E-068.86E-070x=1000004.1.2未知精确解的热传导问题例2是二维Dirichlet边值问题的热传导方程。采用古典显式格式，以求解得，时，关于，的函数如图4-2所式。图4-2：古显格式解例2，关于的图像（，）采用P-R格式，由0变化到1时，不同的t, 关于，的图像（，）对比图如图4-3。图4-3：P-R格式解例2，变化 关于，的图像（，）4.

39、2结果分析由例1，可看出，古显格式和P-R格式精度都较高，且P-R格式较古显格式更接近解析解。由例2，可由结果图4-2看出的大致图形，由图4-3可看出随时间的推移，温度渐渐靠近边界的温度。可见，数值方法也能一定程度上反应解的情况。5稳定性探究与分析5.1稳定性问题的提出考虑数值算例例1，例2，用古典显示格式和P-R格式解时，不同的r可能导致影响格式的稳定性。以例2分析，采用古典显式格式，以求解得，变化时，关于，的函数如图5-1所示。图5-1：古显格式解例2，关于，的图像（，）下面我们先研究准确解和微分方程的解之差，时随着时间层数n的增大而无限增大还是有所控制。如果这种差别是无限增加的，则差分格

40、式不稳定，不稳定的格式是不可用的。5.2 几种分析稳定性的方法随着人们对差分格式的深入研究，发展了以下几种稳定性研究方法。-图研究法是在假定在固定的某节点上引入一个误差，即把改成了+，而在这一层上其他节点的值保持不变，且假定计算时没有引入误差，我们探究此时误差是否会无限增大的方法。稳定性分析的矩阵方法，采用矩阵的2范数，来衡量稳定性。Fourior级数法引进Fourior级数，通过考察增长因子来考察稳定性。Fourior级数法又称为Von-Neumann方法。在判断有限差分近似的稳定性方法中，以Fourior方法使用较为广泛，它仅适用于线性常系数的有限差分近似。以一维热传导方程的显式格式为例，

41、过程如下：第一步：首先，要研究的差分方程可写为：其中，一维热传导方程的古典显式格式则表示为，第二步：其次，对进行变量分离：第三步：将替代为代入所考察的有限差分方程。当即对所有的，时，格式稳定。由于，故时一维热传导方程的古典显式格式稳定。运用Fourior级数法可推导出，一维热传导方程的古典隐式格式、Crank-Nicolson隐式格式无条件稳定。即对任意的r,该格式都稳定。加权六点隐式格式的稳定性条件与有关，如下：如果，则稳定性条件为；如果，则稳定性条件为；如果，则无条件稳定。二维热传导方程的几种差分格式的稳定性也可采用Fourior级数法求得，得出古典显式格式的稳定性条件为，古典隐式格式、C

42、rank-Nicolson隐式、P-R格式、D-R格式、M-F格式无条件稳定。5.3 r变化对稳定性的探究5.3.1 古典显式格式的稳定性根据5.1中对古典显式格式稳定性的分析可知，古典显式格式的稳定性条件是。采用古典显式格式解例2，以求解得，变化时，关于，的函数如上图5-1所示。采用古典显式格式解例1，以求解得，变化时，数值解与解析解间误差如下表5-1所示。表5-1 古显格式解例1的误差（，）误差r=1/8r=2/8r=3/8r=4/8x=0,y=00000x=0.1,y=0.1-0.00013-0.00026-2.4E+26-5.9E+29x=0.2,y=0.2-0.00048-0.000

43、95-4.3E+26-6.1E+29x=0.3,y=0.3-0.0009-0.0018-5.5E+26-5E+29x=0.4,y=0.4-0.00125-0.00248-6.1E+26-5.8E+29x=0.5,y=0.5-0.00138-0.00273-6.3E+26-7.7E+29x=0.6,y=0.6-0.00125-0.002481.11E+27-5E+29x=0.7,y=0.7-0.0009-0.0018-1E+27-5.4E+29x=0.8,y=0.8-0.00048-0.00095-7.8E+26-1E+29x=0.9,y=0.9-0.00013-0.00026-4.3E+26-

44、7.8E+29x=1.0,y=1.00000误差r=5/8r=6/8r=7/8r=1x=0,y=00000x=0.1,y=0.19E+283.42E+324.17E+294.48E+29x=0.2,y=0.26.1E+294.72E+328.93E+293.31E+28x=0.3,y=0.39E+293.86E+329E+295.95E+28x=0.4,y=0.44.4E+292.39E+324.48E+299.57E+29x=0.5,y=0.51.3E+291.74E+324.84E+298.93E+29x=0.6,y=0.64.4E+292.39E+324.84E+297.3E+29x=0

45、.7,y=0.71.7E+293.25E+328.16E+291.54E+29x=0.8,y=0.87.4E+292.58E+322.06E+293.31E+28x=0.9,y=0.99.3E+291.53E+325.41E+298.16E+29x=1.0,y=1.00000根据图5-1，表5-1可看出，当时稳定。否则，格式可能不稳定。且r越大，越容易发生震荡。5.3.2 P-R格式格式的稳定性采用P-R格式解例2，由0变化到1时，不同的时, 关于，的图像（，）对比图如图5-2。采用P-R格式解例1，由0变化到1时，不同的时, 数值解与解析解间误差（，变化）对比如表5-2。图5-2：P-R格式

46、解例2，变化 关于，的图像（，）表5-2 P-R格式解例1的误差（，）误差r=1/8r=2/8r=3/8r=4/8x=0,y=00000x=0.1,y=0.18.86E-071.19E-06-0.000464-0.000489x=0.2,y=0.22.05E-061.19E-06-0.000883-0.000928x=0.3,y=0.31.48E-069.21E-07-0.001214-0.001275x=0.4,y=0.4-8.6E-079.66E-07-0.001424-0.001497x=0.5,y=0.5-1.4E-061.47E-06-0.001497-0.001574x=0.6,y

47、=0.6-8.6E-071.18E-06-0.001424-0.001497x=0.7,y=0.71.48E-061.47E-06-0.001214-0.001275x=0.8,y=0.82.05E-062.05E-06-0.000883-0.000928x=0.9,y=0.98.86E-071.47E-06-0.000464-0.000489x=1.0,y=1.00000续表5-2 P-R格式解例1的误差（，）误差r=5/8r=6/8r=7/8r=1x=0,y=00000x=0.1,y=0.1-0.000464-0.000396-0.000288-0.000152x=0.2,y=0.2-0.

48、000883-0.000753-0.000548-0.000288x=0.3,y=0.3-0.001214-0.001034-0.000753-0.000396x=0.4,y=0.4-0.001424-0.001214-0.000883-0.000464x=0.5,y=0.5-0.001497-0.001275-0.000928-0.000489x=0.6,y=0.6-0.001424-0.001214-0.000694-0.000464x=0.7,y=0.7-0.001214-0.001034-0.000753-0.000396x=0.8,y=0.8-0.000883-0.000753-0.

49、000548-0.000288x=0.9,y=0.9-0.000464-0.000396-0.000288-0.000152x=1.0,y=1.00000由图5-2，表5-2可验证，对P-R格式而言，r取01中的任何数，该格式都是稳定的。结语本文主要讨论抛物方程的有限差分法，系统阐述了抛物方程的的几种数值解法，介绍了有限差分法的基本思想，步骤，和广泛的应用。着重讨论了抛物方程的有限差分解法，分析了这些方法的收敛性和稳定性。以热传导方程为例，编程实现了使用古典显式格式和P-R格式解二维抛物方程。通过数值算例，验证了这两种格式的稳定性条件，当步长比，古典显式格式稳定，否则可能不稳定；对任意步长比r

50、，P-R格式都是稳定的。进一步确定了自己的研究方向：实现差分格式步长自适应调整，可以提高计算结果的准确性，在不知道准确解的情况下，得到较准确的结果。探究新的稳定性好，计算简单的差分格式。并进一步学习有限元法，有限体积法，有限单元法，无网格方法等算法。参考文献1 李荣华偏微分方程数值方法M 北京：高等教育出版社，20052 李国生有限差分法求解静电场问题J电气电子教学学报 2005.05(05):34-383 马小霞抛物型方程的高精度差分格式的构造和理论研究J数学季刊，2002,17(3):57-614 林鹏成多维抛物型方程和双曲方程的差分解法J福州大学学报，1964(2):56-595 张星，

51、单双荣解二阶抛物型方程的一族高精度恒稳格式J华侨大学学报（自然科学版），2009,30(2):69-726 戴嘉尊，邱建贤微分方程数值方法M南京：东南大学出版社 20027 王元明，工程数学 数学物理方程与特殊函数M北京：高等教育出版社.20048 黎杨抛物型方程的有限差分解法及其在电磁场中的应用D 武汉：武汉理工大学理学院，2010.9 单恺婷Continuous Block -Methed And Their AplicationsD济南：山东大学数学学院，200910 赵光明无网格方法及其在冲击动力学中的研究D成都：西南交通大学材料工程学院，200811 李晓红常微分方程数值解法及其应用

52、D南京：东南大学数学学院，200812 马明书解抛物型方程的一个高精度两层显格式J河南师范大学学报(自然科学版)，1996.02:81-9013 孙志忠偏微分方程数值方法J高校计算数学学报，1996,24(2):190-19314 闰秋峰抛物型方程的区域分解直接方法及其应用J北京理工大学学报，2000,20(10): 539-54315 王晓峰，王守印变系数热传导方程的两层绝对稳定差分格式J新乡学院学报(自然科学版). 2009,26(1):69-7616 袁光伟，杭旭登热传导方程基于界面修正的迭代并行计算方法J 数值计算与计算机应用，2006,24(1):67-80附录 P-R格式的C+实现

53、代码#include#includevoid main() / 变量解释： u12121为第n层,u22121为第n+1/2层,u32121为第n+1层 / 变量解释： d20三对角方程右端项,g20,w20追赶法中的变量float u12121=0,u22121=0,u32121=0,d20,g20,w20; float h=3.14/20,K=3.14*3.14/(5.0*20*20),r=1.0/5,a=0,b=0,c=0;int i,j,k;a=r/2,b=1+r,c=r/2; / t=0赋初值for(j=0;j21;j+) for(i=0;i21;i+)u1ij=sin(i*h)+s

54、in(j*h);for(k=0;k1/K;k+) /x=0,20赋初值for(i=0;i21;i+) u20i=u220i=u2i0=u2i20=u30i=u320i=u3i0=u3i20=0;/开始第n层到第n+1/2层的追赶计算（x方向）for(j=1;j20;j+) /追for(i=1;i20;i+)di=(1-r)*u1ij+r/2*(u1ij+1+u1ij-1);for(w1=c/b,g1=d1/b,i=2;i0;i-)u2ij=gi+wi*u2i+1j;for(j=1;j20;j+) /开始第n+1/2层到第n+1层的追赶计算（y方向） /追for(i=1;i20;i+)di=(1-r)*u2ij+r/2*(u2i+1j+u2i-1j);for(w1=c/b,g1=d1/b,i=2;i0;i-)u3ji=gi+wi*u3ji+1; /将第n+1层赋给第n层层