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

本文格式为Word版，下载可任意编辑——二维抛物方程的有限差分法华北电力大学本科毕业设计（论文）

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

有限体积法(FiniteVolumeMethods)又称为控制体积法。其基本思路是：将计算区域划分为一系列互不重叠的控制体，并使每个网格点周边有一个控制体；将待求解的微分方程对每一个控制体积积分，便得出一组离散方程。该方法的未知量为网格点上的函数值。为了求出控制体积的积分，须假定函数值在网格点控制体边界上的变化规律。从积分区域的选取方法看来，有限体积法属于有限元方法中检验函数取分片常数插值的子区域法；从未知量的近似方法看来，有限体积法属于采用局部近似的离散方法[1]。

有限单元法的基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些适合的节点作为求解函数的插值点，借助于变分原理或加权余量法，将微分方程离散求解[1]。无网格方法的近似函数建立在一些离散的节点上，不需要借助网格，在涉及到网格畸变，网格移动等问题中显示出了明显的优势[1]。

1.2.2有限差分方法的发展

随着人们对显示格式的进一步研研究，提出了好多新的高精度，绝对稳定的差分格式。文献2探讨了变系数热传导方程的两层绝对稳定差分格式。研究了二维变系数非齐次热传导方程的两层绝对稳定的差分格式问题。首先运用Pade迫近导出了差分格式,给出了差分格式的截断误差;探讨了差分格式的绝对稳定性和收敛性,且收敛阶为O(k2?h4);最终给出了数值例子,数值结果和理论结果是吻合的。

文献3第一部分用待定系数法对P维抛物型方程(P?1,2,3,4)构造出了高精度(截断误差达O(k2?h4))能显式计算，稳定性较好(r?1/2)的三层(特别状况下可以是两层)显式差分格式，在空间变量更多的状况下，指出了构造高精度差分格式的一般方法。

文献3其次部分用算子方法对二维和三维抛物型方程构造出了高精度(截断误差达

O(k2?h4))的绝对稳定的交替方向隐式差分格式，在空间变量更多的状况下，也指出了构造高精度交替方向隐式差分格式的一般方法。并附有数值例子，它说明理论分析的正确性和所建立的差分格式的有效性。

文献4对二阶抛物型方程构造了一族含参数高精度三层差分格式。当参数满足一定的条件时,差分格式绝对稳定，其局部截断误差阶数最高可达O(k2?h4)。适当地调理参数，可以得到一个七点显式差分格式和一个两层六点隐格式。数值例子说明,对稳定性所作的分析是正确的。

文献5多维抛物型方程和双曲方程的差分解法对一般m维热传导方程，波动方程的初、

2

华北电力大学本科毕业设计（论文）

边值问题，提出若干个交替方向的差分格式。

1.3差分格式建立的基础[6]

1.3.1区域剖分[14]

构造一维线性抛物方程

（x,t）??u??u?u?(（ax,t）)?（bx,t）?（cx,t）u(1-2)?t?t?x?x的有限差分迫近，首先将求解区域?用两组平行于x轴和t轴的直线构成网格覆盖，网格边长在x方向为?x?h，t方向为?t?k，网格节点为(xm,tn)，xm?mh,tn?nk，m,n为整数[6]。网格如图1-1。

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

二维抛物方程的区域剖分将剖分区域扩展到三维空间。首先将求解区域V用三组平行于x轴，y轴和t轴的直线构成网格覆盖，网格边长在x方向为?x?h，y方向为

?y?h，t方向为?t?k，网格节点为(xl,ym,tn)，xl?lh,ym?mh,tn?nk，l,m,n为整数[3]。

网格如图1-2[6]。

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

1.3.2差商代替微商

差分方法的基本思想就是以差商代替微商。考虑如下两个Taylor公式：1113n?n1)?u(t)??u(t)?h??u(2)t?h???u()?t?h?n()u()t?hu(t?h(Oh)(1-3)

2!3!n!3

华北电力大学本科毕业设计（论文）

u(t?h)?u(t)?u?(t)h?111u??(t)h2?u???(t)h3???u(n)(t)hn?O(hn?1)(1-4)2!3!n!u?(ti)?u(ti?1)?u(ti)?O(h)h从式(1-3)得到：

从式(1-4)得到：

u?(ti)?

u(ti?1)?u(ti)?O(h)h从式(1-3)-式(1-4)得到：

u?(ti)?u(ti?1)?u(ti?1)?O(h2)2h从式(1-3)+式(1-4)得到：

u(ti?1)?2u(ti)?u(ti?1)2??u(t)??O(h)i2h下表1-1面介绍常用的几种差商[6]。

表1-1：常用的几种差商形式

有限差分近似公式ur?ur?1（向前差分）hur?1?ur（向后差分）hur?1/2?ur?1/2（中心差分）2h

误差阶O(h)O(h)O(h2)由式(1-1)可得：

n?1n(1-5)um?exp(kL)um1.3.3差商代替微商格式的误差分析

为了分析分析数值方法的确切度，考察收敛性。往往在ui?u(ti)成立的假定下，估计误差ei?1?u(ti?1)?ui?1这种误差称为“局部截断误差〞局部截断误差是以点ti的确切解u(ti)为出发值，用数值方法推进到下一个点ti?1而产生的误差[8]。若如下：

h2h2u(ti?1)?u(ti)?hu?(ti)?u??(?)?u(ti)?hf(ti,u(ti))?u??(?)(1-6)

22h2ei?1?u(ti?1)?ui?1?u??(?)?O(h2)(1-7)

24

华北电力大学本科毕业设计（论文）

局部截断误差是关于hn的微小函数，则收敛[13]。

为了分析分析数值方法的确切度，考察收敛性。整体截断误差是以点t0的初始值的偏差：?i?1?u(ti?1)?ui?1称为整体截断误差[8]。若如下：

为出发值，用数值方法推进i+1步到点ti?1，所得的近似值ui?1与确切值u(ti?1)

ui?1?ui?hf(ti,ui)

?i?1??i?h[f(ti,u(ti)?f(ti,ui)]?Dh2

h2ei?1?u(ti?1)?ui?1?u??(?)?O(h2)

2?i?1??i?h[f(ti,u(ti)?f(ti,ui)]?Dh2

?i?1??i?hL?i?Dh2

?i?1?(1?hL)i?1?0?Dh2[1?(1?hL)?(1?hL)2???(1?hL)i]

?(1?hL)i?1Dh2?0?[(1?hL)i?1?1][14]

hL整体截断误差不随?0无限扩大，则收敛。

1.4本文主要研究内容

本文主要研究二维抛物方程的有限差分方法，研究工作分为以下四个方面：1)差分格式的构造方法

有限差分法法解二维抛物方程包括区域剖分和差商代替导数两个过程。差分格式就是在网格结点上求出微分方程的近似解的一种方法，因此又称为网格法。

2)差分格式的分析方法

数值方法是近似计算方法，需从收敛性，相容性，稳定性方面考察。收敛性是考察步长足够小时，数值解是否无限迫近解析解。稳定性主要最初产生的小误差在以后的计算中虽然会传递下去，但不会无限制的扩大。

3)显式差分格式

显式差分格式是差分方法中可逐层逐点分别求解的格式，是一种有限差分近似方法。显式差分格式的优点在于计算简单，但它并不总是稳定的。

4)隐式差分格式

与显式差分格式不同，隐式差分格式中包括了（n+1）时间层上二个或二个以上结点

n?1n?1n?1处的未知值（例如Um，使用隐式差分格式和使用显式差分格式求解完全不同。,U,U?1mm?1）

相对而言，使用隐式差分格式求解，每时间层包含有较多的计算工作量。从后面对差分格式的稳定性分析可知，隐式格式的优点在于，其稳定性要求对步长比的限制大为放宽，而

5

华北电力大学本科毕业设计（论文）

这正是我们所期望的。

6

华北电力大学本科毕业设计（论文）

2显式差分格式

现在，对二维抛物方程式(1-1)，从方程式(1-4)出发，构造有限差分的显式格式。由于一维热传导方程

2?u2?u?a(2-1)?t?x2

二维热传导方程

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

2.1常系数热传导方程的古典显式格式

首先考虑热传导方程的边值问题

2??u?2u2?u??t?a（?x2??y2)???u(0,y,t)?u(L,y,t)?u(x,0,t)?u(x,L,t)?0?u(x,y,0)?f(x,y)???离散化

u0j，m?u(0,m,jk)?uLj，m?u(L,m,jk)?0

jjul，?u(l,0,jk)?u0l，M?u(M,m,jk)?0

ul0,m?u(lh,mh,0)?f(lh,mh)?fl,m[15]

2.1.1古典显式格式格式的推导

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

?122uln,m?exp(kDx?kDy)uln,m

故

7

华北电力大学本科毕业设计（论文）

11?12244uln,m?uln,m(1?kDx?kDy?k2Dx?k2Dy??)

2211222式中右边假使仅保存二阶导数项，且以2?x2替代Dx，2?y替代Dy，则得差分格式

hhkk?1Uln,m?Uln,m(1?2?x2?2?y2)

hh或者

?1Uln,m?(1-4r)Uln,m?r(Uln?1,m?Uln?1,m?Uln,m?1?Uln,m?1)(2-3)

r?h/k2

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

图2-1：古典显式格式

格式(2-3)应用在一维热传导问题中，得到古典显式格式为：

n?1nnnUm?rUm?1?(1?2r)Um?rUm?1(2-3)

每一层各个节点上的值是通过一个方程组求解得到的。显一隐格式区域分解方法就是以显格式计算出相邻子区域相交内边界的近似值的一种方法。显一隐格式区域分解方法综合了二者的优点，借助前一层数值解的信息，用显格式给出在这—层的子问题的未知内边界条件，把一个整体区域上的问题化为若干个子区域上的子问题，在每个子区域上用隐式方法求解，从而实现了并行。这可以从下面的计算过程看出来[8]。

2.1.3古典显式格式的算法步骤

古典显式格式，以m=1为例，列成矩阵有如下形式：

100000?u1,1?(1?4r)u1,1?r(u0,1?u2,1?u1,0?u1,2)?100000?u2,1?(1?4r)u2,1?r(u1,1?u3,1?u2,0?u2,2)?100000?u3,1?(1?4r)u3,1?r(u2,1?u4,1?u3,0?u3,2)???00000?u1L?1,1?(1?4r)uL?1,1?r(uL?2,1?uL,1?uL?1,0?uL?1,2)?系数矩阵为

8

华北电力大学本科毕业设计（论文）

r?1?4r???r1?4rr???????r1?4rr???r1?4r???为了得到二维问题的误差估计，我们做一些改动。因此，这个算法从un到un?1的计算步骤为：

第一步，根据给出的初边值条件，得出t=0时u0；

?1其次步，根据古典显式格式的公式，Uln,m?(1-4r)Uln,m?r(Uln?1,m?Uln?1,m?Uln,m?1?Uln,m?1)，

由第un得出un?1[6]。

下文中将通过一个具体的数值例子说明白计算的方法，表达了这种格式的实用性和优越性。

9

华北电力大学本科毕业设计（论文）

3隐式差分格式

与显式差分格式不同，隐式差分格式中包括了(n+1)时间层上二个或二个以上结点处的

n?1n?1n?1未知值(例如Um,U,U?1mm?1)，使用隐式差分格式和使用显式差分格式求解完全不同。相

对而言，使用隐式差分格式求解，每时间层包含有较多的计算工作量。从后面对差分格式的稳定性分析可知，隐式格式的优点在于，其稳定性要求对步长比的限制大为放宽，而这正是我们所期望的[6]。

对二维抛物方程式(1-1)，从方程式(1-4)出发，构造有限差分的显式格式。由于热传导方程式(2-1)构造差分格式。我们需要研究的不仅是差分方程本身，而且是包括全部内部区域和所有边界上的差分方程所组成的代数方程组。又称为差分格式。

3.1古典隐式格式

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

?122uln,m?exp(kDx?kDy)uln,m

故

22n?1nexp(?kDx?kDy)um?umun?1m124124?u(1?kD?kD?kDx?kDy??)22nm2x2y式中右边假使仅保存二阶导数项，且以则得差分格式

121212222???替代，且以替代，替代，DDDxxyxxy222hhhk2k2?x?2?y)h2hn?1nUm?Um(1?或者

?11n?1n?1n?1n(1?4r)Uln,m?r(Uln??1,m?Ul?1,m?Ul,m?1?Ul,m?1)?Ul,m(3-1)

r?h/k2

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

n?1nnnUm?rUm?1?(1?2r)Um?rUm?1(3-2)

显式与隐式的比较如下[6]：

(1)同一阶数下，隐式的局部截断误差的系数的绝对值比显式的要小；(2)显式的计算工作量比隐式的小；

10

华北电力大学本科毕业设计（论文）

(3)隐式的稳定范围比显式的大。

局部截断误差的阶最高是3，式(3-2)是允许函数任意变化状况下截断误差最小的二阶方法。要再提高阶就必需增加计算函数值的次数。上述式(3-2)又称为古典隐式格式。

故格式用图3-1表示，其截断误差阶为O(k2?h2),与古典差分格式一致。这是一个四点差分格式，如图3-1所示。

图3-1：古典隐式格式

为了求得第(n+1)时间层上的Un?1的值，必需通过解线性代数方程组。

这是一个隐式差分格式，必需联合其初边值条件求解。格式(3-2)寻常称为古典隐式格

22式。我们也可以通过直接用差分算子代替Dx,Dx,Dy,Dy的方即：

n?1n?un?1ul,m?ul,m()l,m??tkn?1n?1n?1?2un?1ul,m?1?2ul,m?ul,m?1(2)l,m??yh2n?1n?1n?1?2un?1ul?1,m?2ul,m?ul?1,m(2)l,m??xh2代入微分方程，得到格式(3-1)。

古典隐式格式的方程组和矩阵形式如下：

当知道第n层上的uij时，要确定第n+1层上各点值uij?1必需通过求解一个线性代数方程组。以m=1为例：

j?1j?1j?1j?1j?1j?(1?4r)u1,1?r(u0,1?u2,1?u1,0?u1,2)?u1,1?j?1j?1j?1j?1j?1j?(1?4r)u2,1?r(u1,1?u3,1?u2,0?u2,2)?u2?j?1j?1j?1j?1j?1j?(1?4r)u3,1?r(u2,1?u4,1?u3,0?u3,2)?u3,1???j?1j?1j?1j?1j?1j?(1?4r)uN?1,1?r(uN?2,1?uN,1?uN?1,0?uN?1,2)?uN?1,1?其系数矩阵如下：

11

华北电力大学本科毕业设计（论文）

?r?1?4r???r1?4r????????r?????r1?4r?r??r1?4r??3.2Crank-Nicolson隐式格式

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

11?1exp(?kL)uln,m?exp(kL)uln,m

22由

22L?Dx?Dy得

121211221122121211221122n?1[1?kDx?kDy?(kDx)?(kDy)??]um?[1?kDx?kDy?(kDx)?(kDy)222222222222??]uln,m(3-3)两边仅保存前二项，用

121222??y替代Dy代替，，则得差分格式Dxx22hh111?1(1?r?x2?r?y2)Uln,m?(1?r?y2)Uln.m

222这是一个隐式差分格式，称为Crank-Nicolson差分格式，截断误差阶为?(k2?h2)，也可写为

11?1?1n?1n?1n?1nnnnn(1?2r)Uln,m?r(Uln,m?U?U?U)?(1?2r)U?r(U?U?U?U?1l,m?1l?1,ml?1,mml,m?1l,m?1l?1,ml?1,m)

44(3-4)

由于格式（3-4）中包括六个结点，故也可称为六点格式（如图3-2所示）。

12

华北电力大学本科毕业设计（论文）

图3-2：Crank-Nicolson隐式格式

也可将

n?1nu?u?un?1l,ml,m()l,m2??tkn?1n?1n?1nnnu?2u?uu?2u?u?2un?11l?1,ml,ml?1,ml?1,ml,ml?1,m(2)l,m2?(?)22?x2hhn?1n?1n?1nnnu?2u?uu?2u?u?2un?11l,m?1l,ml,m?1l,m?1l,ml,m?1(2)m2?(?)22?y2hh代入微分方程(2-1)，得到Crank-Nicolson格式。解微分方程(2-1)，根据Crank-Nicolson格式得到的方程组：11?1?1n?1n?1n?1n(1?2r)Uln,m?r(Uln,m?U?U?U)?(1?2r)U?r(Uln,m?1?Uln,m?1?Uln?1,m?Uln?1,m)?1l,m?1l?1,ml?1,ml,m44其矩阵表达式为：

j?1j????uu1?r?r/21?rr/211??1，??1，???uj?1????uj???r/21?rr/2?11??r/21?r?r/2??2，??2，????????????????j?1????j?????r/21?r?r/2??uN?1，r/21?rr/2??uN?1，?11?j?1??????j????r1?2r?r/21?2r?????1?1??uN，?uN，3.3Douglas差分格式[6]

基于宛如Crank-Nicolson格式一样的六个网格结点可获得另一精度较高的差分格式，如

22在前式(3-3)中仅保存直到Dx，Dy的项，即有：

112n?1112n(1?kDx2?kDy)ul,m?(1?kDx2?kDy)ul,m

2222可令

Dx2uln,m?2nDyul,m1212n?(1??x)ul,m2xh12

1212n?2?y(1??y)ul,mh12则可得：

1212?1?(1??x)xh212

112Dy?2?y2(1??y2)?1h12Dx2?代入上式，则有如下差分格式

13

华北电力大学本科毕业设计（论文）

11111111?1[1?(r?)?x2?(r?)?y2]Uln,m?[1?(r?)?x2?(r?)?y2]Uln,m(3-5)

26262626它称为Douglas差分格式，具有截断误差阶O(k2?h4)。精度较高的差分格式(Douglas格式),并通过一个具体的数值例子说明白计算的方法,表达了这种格式的实用性和优越性。

3.4加权六点隐式格式[9]

前面，我们已经推导了热传导方程(2-1)的古典显式显示格式，古典隐式格式及Crank-Nicolson格式等。实际上，它们都可以作为本节推导的加权六点隐式格式的特别情形。由

?122uln,m?exp(kDx?kDy)uln,m

得

22?122exp(??kDx??kDy)uln,m?exp[(1??)kDx?(1??)kDx]uln,m,0???1

即

111124?124[1??kDx2??kDy??2k2Dx4??2k2Dy??]uln,m?[1?(1??)kDx2?(1??)kDy?(1??)2k2Dx?22224(1??)2k2Dy??]uln,m

两边去掉高于二阶导数的项，且用

121222??代替，替代，则得差分格式DDxyxy22hh2?12(1??r?x2??r?y)Uln,m?[1?(1??)r?x2?(1??)r?y]Uln,m

或者

?1nnn?1n?1n?1n?1(1?2r?)Uln,m?r(1??)(Uln,m?1?Uln,m?1?Um,l?1?Um,l?1)?r?(Ul,m?1?Ul,m?1?Um,l?1?Um,l?1)?[1?

2r(1??)]Uln,m,0???1(3-5)这是一个六点差分格式（如图3-3所示），称为加权六点差分格式。

?1??

图3-3：加权六点差分格式

14

华北电力大学本科毕业设计（论文）

3.5交替方向隐式格式

交替方向隐式格式是关于x方向的显格式，关于y方向的隐格式。在n+1/2层上用古典显式格式计算出过度值ujn?12，再在n+1层上用古典隐格式校正预计值。得到的稳定性

和收敛性是相近的。但是，只探讨了关于常系数热传导方程的在z和可方向都是一致的时间步长的情形，所以我们在这里试图做一些扩展。我们首先考虑在内边界点用一个方向为显格式，另一个方向为隐格式的情形，再考察在z和y方向都用显格式的情形，而且我们探讨在z和可方向用不同的时间步长的情形[6]。

3.5.1Peaceman-Rachford格式

r?h/k2考虑二维抛物方程(2-1)的差分解法，以上对显式格式的稳定性分析发现，

的要求比一维情形更苛刻。分析说明，维数越高时，要求时间步长越小，计算工作量越大[3]。

(P-R)ADI格式是由Peaceman-Rachford在1995年发表的，从第n层到n+1层，P-R格式分两步进行，每一步只需解具有三对角系数的线性方程组[6]。P-R格式为

(1)?1/2Ul*,n?Uln,mmk/2?1?1/2Uln,m?Ul*,nm?12*n?12n(?U??Ul,m)(3-6)xl,my2h12*n?1/22n?1(?xUl,m??yUl,m)(3-7)2h(2)k/2?nnn，?yul,m?ul?1/2,m?ul?1/2,m?x，?y是中心差算子，?xuln,m?uln,m?1/2?unl,m?1/2P-R格式对任意r>0无条件稳定。但预计P-R格式对三个变量空间问题，无条件稳定性不再成立[3]。

3.5.2Rachford-Mitchell格式

Douglas和Rachford在1956年提出另一个隐式差分格式，即Douglas-Rachford格式

[6]

。D-R格式是第一个能被推广到三维情形的交替方向隐式格式，二维D-R格式是

?1nUl*,nm?Ul,mk?1?1Uln,m?Ul*,nm??1*n?1*n?1Ul*?n1,m?2Ul,m?Ul?1,mh21n?1n?1Uln??1,m?2Ul,m?Ul?1,m?Uln,m?1?2Uln,m?Uln,m?1h2Uln,m?1?2Uln,m?Uln,m?1h2(3-7)

k?h2?(3-8)

3.5.3Mitchell-Fairweather格式

Mitchell和Fairweather在1964年推导了一个高精度ADI差分格式，称为M-F格

15

华北电力大学本科毕业设计（论文）

式[6]。M-F格式截断误差达O(k2?h4)较之P-R格式和D-R格式更好，且无条件收敛[3]。

3.5.4交替方向隐式格式的算法步骤

ADI格式的算法步骤与古典显式格式的算法步骤相像，只是ADI格式是隐式格式，需用追赶法求解。追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。

第一步，根据给出的初边值条件，得出t=0时u0;

其次步，根据各个交替方向隐式格式中的第一个公式，由第un求得出un?1/2如式(3-6)，

1采用追赶法，确定追赶因子a=1+r/2,b=r/2,c=1+r/2,d=(1?2r)Uln,m?r(Uln,m?1?Uln,m?1?Uln?1,m,

4nnn0n+1/2n+1/2n+1/2=d-*;?Ul?1,m)wm=b-a*u,gm=(d-u*a)/w,然后根据追赶因子确定u。uM=gM,umgu?1mm第三步，再由追赶法由第un?1/2求得出un?1[6]。

由该算法，计算un的近似解的过程。下文中将通过一个具体的数值例子说明白计算的方法,表达了这种格式的实用性和优越性。

16

华北电力大学本科毕业设计（论文）

华北电力大学本科毕业设计（论文）

4实例分析与结果分析

4.1算例

4.1.1已知有确切解的热传导问题

例1是二维Dirichlet边值问题的热传导方程

??T?2T?2T?2?2(0?x?1;0?y?1;0?t?T)??t?x?y??(0?x?1;0?y?1)?T(x,y,0)?0?T(0,y,t)?T(1,y,t)?T(x,1,t)?T(x,0,t)?1(0?t?T)???方程的确切解为

exp(??2(k2?l2)t)T(x,y,t)?1?2??sin(k?x)sin(l?y)

?k?1,3,?l?1,3,?kl16??采用古典显示格式，以求解得t?0.1，h?1/20，r?h/k2?1/8时，T关于x，y的函数如图4-1所式。

2图4-1：古显格式解例1，T关于x，y的图像（t?0.1，h?1/20，r?h/k?1/8）

分别以古典显示格式，P-R(ADI)格式，以一致的步长(h?1/20，r?h/k2?1/8)求解

17

华北电力大学本科毕业设计（论文）

同一处(t?0.1)的T与解析解间的差值，列在误差表中如表4-1，表4-1所示。

表4-1古显格式解例1的误差（t?0.1，h?1/20，r?h/k2?1/8）

误差y=0y=0.1y=0.2y=0.3y=0.4x=000000x=0.10-0.00013-0.00025-0.00035-0.00041x=0.20-0.00025-0.00048-0.00066-0.00077x=0.30-0.00035-0.00066-0.0009-0.00106x=0.40-0.00041-0.00077-0.00106-0.00125x=0.50-0.00043-0.00081-0.00112-0.00131x=0.60-0.00041-0.00077-0.00106-0.00125x=0.70-0.00035-0.00066-0.0009-0.00106x=0.80-0.00025-0.00048-0.00066-0.00077x=0.90-0.00013

-0.00025

-0.00035

-0.00041

x=100000续表4-1古显格式解例1的误差（t?0.1，h?1/20，r?h/k2?1/8）

误差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-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-2P-R(ADI)格式解例1的误差（t?0.1，h?1/20，r?h/k2?1/8）误差y=0y=0.1y=0.2y=0.3y=0.4x=000000x=0.108.86E-071.19E-061.47E-066.44E-07x=0.201.19E-062.05E-062.19E-069.66E-07x=0.301.47E-062.19E-061.48E-069.21E-07x=0.406.44E-079.66E-079.21E-07-8.6E-07x=0.501.18E-069.93E-073.37E-09-1.4E-06x=0.606.44E-079.66E-079.21E-07-8.6E-07x=0.701.47E-062.19E-061.48E-069.21E-07x=0.801.19E-062.05E-062.19E-069.66E-07x=0.908.86E-071.19E-061.47E-066.44E-07x=100

000

18

y=0.50-0.00043-0.00081-0.00112-0.00131-0.00138-0.00131-0.00112-0.00081-0.00043

0

y=0.501.18E-069.93E-073.37E-09-1.4E-06-1.4E-06-1.4E-063.37E-099.93E-071.18E-060

华北电力大学本科毕业设计（论文）

2续表4-2P-R(ADI)格式解例1的误差（t?0.1，h?1/20，r?h/k?1/8）

误差x=0x=0.1x=0.2x=0.3x=0.4x=0.5x=0.6x=0.7x=0.8x=0.9x=1

y=0.606.44E-079.66E-079.21E-07-8.6E-07-1.4E-06-8.6E-079.21E-079.66E-076.44E-070y=0.701.47E-062.19E-061.48E-069.21E-073.37E-099.21E-071.48E-062.19E-061.47E-060

y=0.801.19E-062.05E-062.19E-069.66E-079.93E-079.66E-072.19E-062.05E-061.19E-060

y=0.908.86E-071.19E-061.47E-066.44E-071.18E-066.44E-071.47E-061.19E-068.86E-070y=100000000000

4.1.2未知确切解的热传导问题

例2是二维Dirichlet边值问题的热传导方程。

??u?2u?2u??2?2?t?x?y???u(x,y,0)?sin(x)?u(0,y,t)?u(?,y,t)?u(x,?,t)?u(x,?,t)?0???0?x,y??;0?t?1采用古典显式格式，以求解得t?1，h??/20，r?h/k2?1/5时，u关于x，y的函数如图4-2所式。

2图4-2：古显格式解例2，u关于x,y的图像（t?1，h??/20，r?h/k?1/5）

19

华北电力大学本科毕业设计（论文）

采用P-R格式，t由0变化到1时，不同的t,u关于x，y的图像（t?1，h??/20，

r?h/k2?1/8）对比图如图4-3。

2图4-3：P-R格式解例2，t变化u关于x，y的图像（h??/20，r?h/k?1/5）

4.2结果分析

由例1，可看出，古显格式和P-R格式精度都较高，且P-R格式较古显格式更接近解析解。由例2，可由结果图4-2看出u的大致图形，由图4-3可看出u随时间的推移，温度渐渐靠近边界的温度。可见，数值方法也能一定程度上反应解的状况。

20

华北电力大学本科毕业设计（论文）

5稳定性探究与分析

5.1稳定性问题的提出

考虑数值算例例1，例2，用古典显示格式和P-R格式解时，不同的r可能导致影响格式的稳定性。以例2分析，采用古典显式格式，以求解得t?0.1，h?1/20，r?h/k2变化时，T关于x，y的函数如图5-1所示。

图5-1：古显格式解例2，u关于x，y的图像（t?0.1，h??/20）

下面我们先研究确凿解和微分方程的解之差，时随着时间层数n的增大而无限增大还是有所控制。假使这种区别是无限增加的，则差分格式不稳定，不稳定的格式是不可用的。

5.2几种分析稳定性的方法

随着人们对差分格式的深入研究，发展了以下几种稳定性研究方法。

00

即把Um改成了Um+?，而在?-图研究法是在假定在固定的某节点上引入一个误差?，

00

这一层上其他节点的值保持不变，且假定计算时没有引入误差，我们探究此时误差是否会无限增大的方法。稳定性分析的矩阵方法，采用矩阵的2范数，来衡量稳定性。

Fourior级数法引进Fourior级数，通过考察增长因子来考察稳定性。Fourior级数法又称为Von-Neumann方法。在判断有限差分近似的稳定性方法中，以Fourior方法使用较为

21

华北电力大学本科毕业设计（论文）

广泛，它仅适用于线性常系数的有限差分近似。以一维热传导方程的显式格式为例，过程如下：

第一步：首先，要研究的差分方程可写为：

m?N1?aumn?1j?m??bmunj?m

m?N0其中N0=｛-1，0，1｝，N1=｛0｝

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

?1nnnun?ru?(1?2r)u?rujj?1jj?1

其次步：其次，对uij进行变量分开：

nunj??Vpexp[ip2?pxj]ln第三步：将unj替代为Vexp[i?xj]代入所考察的有限差分方程。

Vn?1exp[i?xj]?rVnexp[i?xj?1]?(1?2r)Vnexp[i?xj]?rVnexp[i?xj?1]

Vn?1?rVnexp[i?h]?(1?2r)Vn?rVnexp[?i?h]

Vn?12?h?1?4rsin?r(cos?h?isin?h)?(1?2r)?r(cos?h?isin?h)?1?2r(1?cos?h)

2Vn当

Vn?1?1Vn即对所有的?1?1?4ssin2?h2?1，时，格式稳定。由于0?sin2?h2?1，故r?1时一维热传2导方程的古典显式格式稳定。

运用Fourior级数法可推导出，一维热传导方程的古典隐式格式、Crank-Nicolson隐式格式无条件稳定。即对任意的r,该格式都稳定。加权六点隐式格式的稳定性条件与?有关，如下：

假使r?假使??假使??1，则稳定性条件为??0；211，则稳定性条件为r?；22(1?2?)1，则无条件稳定。21，古典隐式格式、Crank-Nicolson隐式、P-R格式、D-R格式、422

二维热传导方程的几种差分格式的稳定性也可采用Fourior级数法求得，得出古典显式格式的稳定性条件为r?华北电力大学本科毕业设计（论文）

M-F格式无条件稳定。

5.3r变化对稳定性的探究

5.3.1古典显式格式的稳定性

根据5.1中对古典显式格式稳定性的分析可知，古典显式格式的稳定性条件是

r?h/k2?1/4。采用古典显式格式解例2，以求解得t?0.1，h?1/20，r?h/k2变化时，T关于x，y的函数如上图5-1所示。采用古典显式格式解例1，以求解得t?0.1，h?1/20，

r?h/k2变化时，数值解与解析解间误差如下表5-1所示。

表5-1古显格式解例1的误差（t?0.1，h?1/20）

r=1/8r=2/8r=3/8000-0.00013-0.00026-2.4E+26-0.00048-0.00095-4.3E+26-0.0009-0.0018-5.5E+26-0.00125-0.00248-6.1E+26-0.00138-0.00273-6.3E+26-0.00125-0.002481.11E+27-0.0009-0.0018-1E+27-0.00048-0.00095-7.8E+26-0.00013-0.00026-4.3E+26000

r=5/809E+286.1E+299E+294.4E+291.3E+294.4E+291.7E+297.4E+299.3E+290

r=6/803.42E+324.72E+323.86E+322.39E+321.74E+322.39E+323.25E+322.58E+321.53E+320

r=7/804.17E+298.93E+299E+294.48E+294.84E+294.84E+298.16E+292.06E+295.41E+290

误差x=0,y=0

x=0.1,y=0.1x=0.2,y=0.2x=0.3,y=0.3x=0.4,y=0.4x=0.5,y=0.5x=0.6,y=0.6x=0.7,y=0.7x=0.8,y=0.8x=0.9,y=0.9x=1.0,y=1.0r=4/80-5.9E+29-6.1E+29-5E+29-5.8E+29-7.7E+29-5E+29-5.4E+29-1E+29-7.8E+290r=104.48E+293.31E+285.95E+289.57E+298.93E+297.3E+291.54E+293.31E+288.16E+29

0

误差x=0,y=0x=0.1,y=0.1x=0.2,y=0.2x=0.3,y=0.3x=0.4,y=0.4x=0.5,y=0.5x=0.6,y=0.6x=0.7,y=0.7x=0.8,y=0.8x=0.9