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

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

关键词：二维抛物方程；有限差分法；古典显式格式；交替方向隐式格式

FINITE DIFFERENCE METHOD FOR TWO-DIMENSIONAL PARABOLIC EQUATION
Abstract
Two-dimensional parabolic equation is a widely used class of partial differential equations. 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-dimensional 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 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 direction 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 differential 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; Eclassical Explicit Scheme; Alternating Direction Implicit Scheme

II

华北电力大学本科毕业设计（论文）
目 录
摘要 I
Abstract II
1绪论 1
1.1课题背景 1
1.2发展概况 2
1.2.1抛物型方程的常见数值解法 2
1.2.2有限差分方法的发展 2
1.3差分格式建立的基础 3
1.3.1区域剖分 3
1.3.2差商代替微商 4
1.3.3差商代替微商格式的误差分析 5
1.4本文主要研究内容 5

2显式差分格式 7
2.1常系数热传导方程的古典显式格式 7
2.1.1古典显式格式格式的推导 7
2.1.3古典显式格式的算法步骤 8
3隐式差分格式 10
3.1古典隐式格式 10
3.2 Crank-Nicolson隐式格式 11
3.3 Douglas差分格式 13
3.4加权六点隐式格式 13
3.5交替方向隐式格式 14
3.5.1 Peaceman-Rachford格式 14
3.5.2 Rachford-Mitchell格式 15
3.5.3 Mitchell-Fairweather格式 15
3.5.4交替方向隐式格式的算法步骤 15
4实例分析与结果分析 16
4.1算例 16
4.1.1已知有精确解的热传导问题 16
4.1.2未知精确解的热传导问题 18
4.2结果分析 19
5稳定性探究与分析 20
5.1稳定性问题的提出 20
5.2 几种分析稳定性的方法 20
5.3 r变化对稳定性的探究 21
5.3.1 古典显式格式的稳定性 21
5.3.2 P-R格式格式的稳定性 22
结语 24
参考文献 25
附录 P-R格式的C++实现代码 26
致谢 28

华北电力大学本科毕业设计（论文）
1绪论
1.1课题背景
抛物方程是一类特殊的偏微分方程，二维抛物方程的一般形式为
(1-1)
其中

。[1]
渗流、扩散和热传导、静电场等很多领域经常会遇到求解二维抛物型方程。微分方程愈复杂，找出解的解析表达式愈困难。对大部分的抛物方程而言，找出解的解析表达式及其困难。因此求出抛物方程的近似解是很有意义的[1]。
本文研究的近似方法是数值方法，目标在于给出解在一些离散点上的近似解值。常见的数值求解方法有有限差分法，有限元法，有限体积法等。有限差分法(Finite Difference Methed)是在1966年Yee提出的，用于解决电机工程中电磁场的初值和边值问题[2]。有限差分法得到迅速发展后，不仅广泛应用于自然科学，在社会科学的各领域也在越来越多地被应用着[2]。某些常用的数值解法，欧拉方法，隐式欧拉法，一般单步法，Crank-Nicolson隐式格式，Runge-Kutta方法等，已成为微分方程数值解领域的经典方法.，在工程计算中的应用随处可见[1]；基于这些方法，人们也在做更深的研究。
采用有限差分法解决二维抛物方程，一些研究者采用一维抛物型方程的古典显格式，古典隐格式，Crank-Nicolson隐式格式，加权六点隐格式等的直接推广；还有有一些研究者采用交叉方向的隐式差分格式(Alternating Direction Implicit Scheme)等基于二维的方法。
1.2发展概况
1.2.1抛物型方程的常见数值解法
抛物型方程的常见数值解法有有限差分法，有限元法，有限体积法，有限单元法，无网格方法等。
有限差分法(Finite Difference Methed)的主要思想是将问题离散化，将问题离散化，用差商代替微商，将微分方程和定界问题都用代数方程来代替。然后解这些代数问题构成的方程组。该方法包括区域剖分和差商代替导数两个过程。具体地，首先将求解区域划分为差分网格，用有限个网格节点代替连续的求解区域。其次，利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替来进行离散，从而建立以网格节点上的值为未知量的代数方程组[1]。
有限元方法(Finite Element Methods)的基础是变分原理和分片多项式插值。该方法的构造过程包括以下三个步骤。首先，利用变分原理得到偏微分方程的弱形式(利用泛函分析的知识将求解空间扩大)。其次，将计算区域划分为有限个互不重叠的单元(三角形、四边形、四面体、六面体等)。再次，在每个单元内选择合适的节点作为求解函数的插值点，将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式，得到微分方程的离散形式[1]。
有限体积法(Finite Volume Methods)又称为控制体积法。其基本思路是：将计算区域划分为一系列互不重叠的控制体，并使每个网格点周围有一个控制体；将待求解的微分方程对每一个控制体积积分，便得出一组离散方程。该方法的未知量为网格点上的函数值。为了求出控制体积的积分，须假定函数值在网格点控制体边界上的变化规律。从积分区域的选取方法看来，有限体积法属于有限元方法中检验函数取分片常数插值的子区域法；从未知量的近似方法看来，有限体积法属于采用局部近似的离散方法[1]。
有限单元法的基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，借助于变分原理或加权余量法，将微分方程离散求解[1]。无网格方法的近似函数建立在一些离散的节点上，不需要借助网格，在涉及到网格畸变，网格移动等问题中显示出了明显的优势[1]。
1.2.2有限差分方法的发展
随着人们对显示格式的进一步研研究，提出了很多新的高精度，绝对稳定的差分格式。文献2讨论了变系数热传导方程的两层绝对稳定差分格式。研究了二维变系数非齐次热传导方程的两层绝对稳定的差分格式问题。首先运用Pade逼近导出了差分格式,给出了差分格式的截断误差;讨论了差分格式的绝对稳定性和收敛性,且收敛阶为;最后给出了数值例子,数值结果和理论结果是吻合的。
文献3第一部分用待定系数法对维抛物型方程构造出了高精度(截断误差达)能显式计算，稳定性较好()的三层(特殊情况下可以是两层)显式差分格式，在空间变量更多的情况下，指出了构造高精度差分格式的一般方法。
文献3第二部分用算子方法对二维和三维抛物型方程构造出了高精度(截断误差达)的绝对稳定的交替方向隐式差分格式，在空间变量更多的情况下，也指出了构造高精度交替方向隐式差分格式的一般方法。并附有数值例子，它表明理论分析的正确性和所建立的差分格式的有效性。
文献4对二阶抛物型方程构造了一族含参数高精度三层差分格式。当参数满足一定的条件时,差分格式绝对稳定，其局部截断误差阶数最高可达。适当地调节参数，可以得到一个七点显式差分格式和一个两层六点隐格式。数值例子表明,对稳定性所作的分析是正确的。
文献5多维抛物型方程和双曲方程的差分解法对一般m维热传导方程，波动方程的初、边值问题，提出若干个交替方向的差分格式。
1.3差分格式建立的基础[6]
1.3.1区域剖分[14]
构造一维线性抛物方程
(1-2)
的有限差分逼近，首先将求解区域用两组平行于轴和轴的直线构成网格覆盖，网格边长在方向为，方向为，网格节点为，，为整数[6]。网格如图1-1。

图1-1 一维抛物方程的网格剖分
二维抛物方程的区域剖分将剖分区域扩展到三维空间。首先将求解区域用三组平行于轴，轴和轴的直线构成网格覆盖，网格边长在方向为，方向为，方向为，网格节点为，，为整数[3]。网格如图1-2[6]。

图1-2 二维抛物方程的网格剖分

1.3.2差商代替微商
差分方法的基本思想就是以差商代替微商。考虑如下两个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)
局部截断误差是关于的极小函数，则收敛[13]。
为了分析分析数值方法的精确度，考察收敛性。整体截断误差是以点的初始值为出发值，用数值方法推进i+1步到点，所得的近似值与精确值 的偏差：称为整体截断误差[8]。若如下：

[14]
整体截断误差不随无限扩大，则收敛。
1.4本文主要研究内容
本文主要研究二维抛物方程的有限差分方法，研究工作分为以下四个方面：
1)差分格式的构造方法
有限差分法法解二维抛物方程包括区域剖分和差商代替导数两个过程。差分格式就是在网格结点上求出微分方程的近似解的一种方法，因此又称为网格法。
2)差分格式的分析方法
数值方法是近似计算方法，需从收敛性，相容性，稳定性方面考察。收敛性是考察步长足够小时，数值解是否无限逼近 解析解。稳定性主要最初产生的小误差在以后的计算中虽然会传递下去，但不会无限制的扩大。
3)显式差分格式
显式差分格式是差分方法中可逐层逐点分别求解的格式，是一种有限差分近似方法。显式差分格式的优点在于计算简单，但它并不总是稳定的。
4)隐式差分格式
与显式差分格式不同，隐式差分格式中包括了（n+1）时间层上二个或二个以上结点处的未知值（例如 ），使用隐式差分格式和使用显式差分格式求解完全不同。相对而言，使用隐式差分格式求解，每时间层包含有较多的计算工作量。从后面对差分格式的稳定性分析可知，隐式格式的优点在于，其稳定性要求对步长比的限制大为放宽，而这正是我们所期望的。
2显式差分格式
现在，对二维抛物方程式(1-1) ，从方程式(1-4)出发，构造有限差分的显式格式。由于一维热传导方程
(2-1)

二维热传导方程
(2-2)
在热传导，磁扩散等许多领域有重要的应用。实际导热问题必然涉及边值条件，在有限差分法中它们也必须差分化。因此，我们需要研究的不仅是差分方程本身，而且是包括全部内部区域和所有边界上的差分方程所组成的代数方程组。又称为差分格式[3]。

2.1常系数热传导方程的古典显式格式
首先考虑热传导方程的边值问题

离散化

[15]

2.1.1古典显式格式格式的推导
现在对热传导方程式(2-1) 推导其最简单的显隐式差分逼近——古典显隐式格式。
由

故

式中右边如果仅保留二阶导数项，且以替代，替代，则得差分格式

或者
(2-3)

这是一个显式格式(四点格式)，如图(2-1)所示。

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

系数矩阵为

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

3.1古典隐式格式
现在对热传导方程推导其最简单的隐式差分逼近——古典隐式格式。由

故

式中右边如果仅保留二阶导数项，且以替代，且以替代，替代，则得差分格式

或者
(3-1)

将格式(2-3)应用于一维热传导方程，古典显式格式为：
(3-2)

显式与隐式的比较如下[6]：
(1)同一阶数下，隐式的局部截断误差的系数的绝对值比显式的要小；
(2)显式的计算工作量比隐式的小；
(3)隐式的稳定范围比显式的大。
局部截断误差的阶最高是3，式(3-2) 是允许函数任意变化情况下截断误差最小的二阶方法。要再提高阶就必须增加计算函数值的次数。上述式(3-2)又称为古典隐式格式。
故格式用图3-1表示，其截断误差阶为,与古典差分格式相同。这是一个四点差分格式，如图3-1所示。

图3-1：古典隐式格式
为了求得第(n+1)时间层上的的值，必须通过解线性代数方程组。
这是一个隐式差分格式，必须联合其初边值条件求解。格式(3-2)通常称为古典隐式格式。我们也可以通过直接用差分算子代替 的方即：

代入微分方程，得到格式(3-1)。
古典隐式格式的方程组和矩阵形式如下：
当知道第n层上的时，要确定第n+1层上各点值必须通过求解一个线性代数方程组。以m=1为例：

其系数矩阵如下：

3.2 Crank-Nicolson隐式格式
Crank-Nicolson隐式差分格式是解热传导方程(2-1)的常用的差分格式，近年来，关于抛物方程的区域分解算法作为并行计算的一种有效工具，吸引了很多学者的注意。这类算法的主要困难在于；如何定义内边界点的值和在子区域上选取合理的计算解去近似。为了推导它，由式（1-4），有

由

得
(3-3)
两边仅保留前二项，用代替，替代，则得差分格式

这是一个隐式差分格式，称为Crank-Nicolson差分格式，截断误差阶为，也可写为
(3-4)
由于格式（3-4）中包括六个结点，故也可称为六点格式（如图3-2所示）。

图3-2： Crank-Nicolson隐式格式
也可将

代入微分方程(2-1)，得到Crank-Nicolson格式。解微分方程(2-1)，根据Crank-Nicolson格式得到的方程组：

 其矩阵表达式为：

3.3 Douglas差分格式[6]
基于如同Crank-Nicolson格式一样的六个网格结点可获得另一精度较高的差分格式，如在前式(3-3)中仅保留直到，的项，即有：

可令

则可得：

代入上式，则有如下差分格式
(3-5)
它称为Douglas差分格式，具有截断误差阶。精度较高的差分格式(Douglas格式),并通过一个具体的数值例子说明了计算的方法,体现了这种格式的实用性和优越性。
3.4加权六点隐式格式[9]
前面，我们已经推导了热传导方程(2-1)的古典显式显示格式，古典隐式格式及Crank-Nicolson格式等。实际上，它们都可以作为本节推导的加权六点隐式格式的特殊情形。由

得

即

两边去掉高于二阶导数的项，且用代替，替代，则得差分格式

或者

(3-5)
这是一个六点差分格式（如图3-3所示），称为加权六点差分格式。

图3-3：加权六点差分格式
3.5交替方向隐式格式
交替方向隐式格式是关于x方向的显格式，关于y方向的隐格式。在n+1/2层上用古典显式格式计算出过度值，再在n+1层上用古典隐格式校正预测值。得到的稳定性和收敛性是相近的。但是，只讨论了关于常系数热传导方程的在z和可方向都是相同的时间步长的情形，所以我们在这里试图做一些扩展。我们首先考虑在内边界点用一个方向为显格式，另一个方向为隐格式的情形，再考察在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格式对任意无条件稳定。但预测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 交替方向隐式格式的算法步骤
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算例
4.1.1已知有精确解的热传导问题
例1是二维Dirichlet边值问题的热传导方程

方程的精确解为

采用古典显示格式，以求解得，，时，关于，的函数如图4-1所式。

图4-1：古显格式解例1，关于，的图像（，，）
分别以古典显示格式，P-R(ADI)格式，以相同的步长(，)求解同一处()的与解析解间的差值，列在误差表中如表4-1，表4-1所示。
表4-1 古显格式解例1的误差（，，）
误差
y=0
y=0.1
y=0.2
y=0.3
y=0.4
y=0.5
x=0
0
0
0
0
0
0
x=0.1
0
-0.00013
-0.00025
-0.00035
-0.00041
-0.00043
x=0.2
0
-0.00025
-0.00048
-0.00066
-0.00077
-0.00081
x=0.3
0
-0.00035
-0.00066
-0.0009
-0.00106
-0.00112
x=0.4
0
-0.00041
-0.00077
-0.00106
-0.00125
-0.00131
x=0.5
0
-0.00043
-0.00081
-0.00112
-0.00131
-0.00138
x=0.6
0
-0.00041
-0.00077
-0.00106
-0.00125
-0.00131
x=0.7
0
-0.00035
-0.00066
-0.0009
-0.00106
-0.00112
x=0.8
0
-0.00025
-0.00048
-0.00066
-0.00077
-0.00081
x=0.9
0
-0.00013
-0.00025
-0.00035
-0.00041
-0.00043
x=1
0
0
0
0
0
0
续表4-1 古显格式解例1的误差（，，）
误差
y=0.6
y=0.7
y=0.8
y=0.9
y=1
x=0
0
0
0
0
0
x=0.1
-0.00041
-0.00035
-0.00025
-0.00013
0
x=0.2
-0.00077
-0.00066
-0.00048
-0.00025
0
x=0.3
-0.00106
-0.0009
-0.00066
-0.00035
0
x=0.4
-0.00125
-0.00106
-0.00077
-0.00041
0
x=0.5
-0.00131
-0.00112
-0.00081
-0.00043
0
x=0.6
-0.00125
-0.00106
-0.00077
-0.00041
0
x=0.7
-0.00106
-0.0009
-0.00066
-0.00035
0
x=0.8
-0.00077
-0.00066
-0.00048
-0.00025
0
x=0.9
-0.00041
-0.00035
-0.00025
-0.00013
0
x=1
0
0
0
0
0
表4-2 P-R(ADI)格式解例1的误差（，，）
误差
y=0
y=0.1
y=0.2
y=0.3
y=0.4
y=0.5
x=0
0
0
0
0
0
0
x=0.1
0
8.86E-07
1.19E-06
1.47E-06
6.44E-07
1.18E-06
x=0.2
0
1.19E-06
2.05E-06
2.19E-06
9.66E-07
9.93E-07
x=0.3
0
1.47E-06
2.19E-06
1.48E-06
9.21E-07
3.37E-09
x=0.4
0
6.44E-07
9.66E-07
9.21E-07
8.6E-07
-1.4E-06
x=0.5
0
1.18E-06
9.93E-07
3.37E-09
-1.4E-06
-1.4E-06
x=0.6
0
6.44E-07
9.66E-07
9.21E-07
-8.6E-07
-1.4E-06
x=0.7
0
1.47E-06
2.19E-06
1.48E-06
9.21E-07
3.37E-09
x=0.8
0
1.19E-06
2.05E-06
2.19E-06
9.66E-07
9.93E-07
x=0.9
0
8.86E-07
1.19E-06
1.47E-06
6.44E-07
1.18E-06
x=1
0
0
0
0
0
0

续表4-2 P-R(ADI)格式解例1的误差（，，）
误差
y=0.6
y=0.7
y=0.8
y=0.9
y=1
x=0
0
0
0
0
0
x=0.1
6.44E-07
1.47E-06
1.19E-06
8.86E-07
0
x=0.2
9.66E-07
2.19E-06
2.05E-06
1.19E-06
0
x=0.3
9.21E-07
1.48E-06
2.19E-06
1.47E-06
0
x=0.4
-8.6E-07
9.21E-07
9.66E-07
6.44E-07
0
x=0.5
-1.4E-06
3.37E-09
9.93E-07
1.18E-06
0
x=0.6
-8.6E-07
9.21E-07
9.66E-07
6.44E-07
0
x=0.7
9.21E-07
1.48E-06
2.19E-06
1.47E-06
0
x=0.8
9.66E-07
2.19E-06
2.05E-06
1.19E-06
0
x=0.9
6.44E-07
1.47E-06
1.19E-06
8.86E-07
0
x=1
0
0
0
0
0

4.1.2未知精确解的热传导问题
例2是二维Dirichlet边值问题的热传导方程。

采用古典显式格式，以求解得，，时，关于，的函数如图4-2所式。

图4-2：古显格式解例2，关于的图像（，，）
采用P-R格式，由0变化到1时，不同的t, 关于，的图像（，，）对比图如图4-3。

图4-3：P-R格式解例2，变化 关于，的图像（，）
4.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的增大而无限增大还是有所控制。如果这种差别是无限增加的，则差分格式不稳定，不稳定的格式是不可用的。
5.2 几种分析稳定性的方法
随着人们对差分格式的深入研究，发展了以下几种稳定性研究方法。
-图研究法是在假定在固定的某节点上引入一个误差，即把改成了+，而在这一层上其他节点的值保持不变，且假定计算时没有引入误差，我们探究此时误差是否会无限增大的方法。稳定性分析的矩阵方法，采用矩阵的2范数，来衡量稳定性。
Fourior级数法引进Fourior级数，通过考察增长因子来考察稳定性。Fourior级数法又称为Von-Neumann方法。在判断有限差分近似的稳定性方法中，以Fourior方法使用较为
广泛，它仅适用于线性常系数的有限差分近似。以一维热传导方程的显式格式为例，过程如下：
第一步：首先，要研究的差分方程可写为：

其中，
一维热传导方程的古典显式格式则表示为，

第二步：其次，对进行变量分离：

第三步：将替代为代入所考察的有限差分方程。

当

即对所有的，时，格式稳定。由于，故时一维热传导方程的古典显式格式稳定。
运用Fourior级数法可推导出，一维热传导方程的古典隐式格式、Crank-Nicolson隐式格式无条件稳定。即对任意的r,该格式都稳定。加权六点隐式格式的稳定性条件与有关，如下：
如果，则稳定性条件为；
如果，则稳定性条件为；
如果，则无条件稳定。
二维热传导方程的几种差分格式的稳定性也可采用Fourior级数法求得，得出古典显式格式的稳定性条件为，古典隐式格式、Crank-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/8
r=2/8
r=3/8
r=4/8
x=0,y=0
0
0
0
0
x=0.1,y=0.1
-0.00013
-0.00026
-2.4E+26
-5.9E+29
x=0.2,y=0.2
-0.00048
-