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

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

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

摘要

二维抛物方程是一类有广泛应用的偏微分方程，由于大部分抛物方程都难以求得解析解，故考虑采用数值方法求解。有限差分法是最简单又极为重要的解微分方程的数值方法。本文介绍了二维抛物方程的有限差分法。

首先，简单介绍了抛物方程的应用背景，解抛物方程的常见数值方法，有限差分法的产生背景和发展应用。讨论了抛物方程的有限差分法建立的基础，并介绍了有限差分方法的收敛性和稳定性。其次，介绍了几种常用的差分格式，有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等，重点介绍了古典显式格式和交替方向隐式格式。进行了格式的推导，分析了格式的收敛性、稳定性。并以热传导方程为数值算例，运用差分方法求解。通过数值算例，得出古典显式格式计算起来较简单，但稳定性条件较苛刻；而交替方向隐式格式无条件稳定。

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

I

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

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发展概况 ................................................................................................................................... 1 1.2.1抛物型方程的常见数值解法 ................................................................................................ 1 1.2.2有限差分方法的发展 ............................................................................................................ 2 1.3差分格式建立的基础 ............................................................................................................... 3 1.3.1区域剖分 ................................................................................................................................ 3 1.3.2差商代替微商 ........................................................................................................................ 3 1.3.3差商代替微商格式的误差分析 ............................................................................................ 4 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隐式格式 ...................................................................................................... 12 3.3 Douglas差分格式 ................................................................................................................... 13 3.4加权六点隐式格式 ................................................................................................................. 14 3.5交替方向隐式格式 ................................................................................................................. 15 3.5.1 Peaceman-Rachford格式 .................................................................................................... 15 3.5.2 Rachford-Mitchell格式 ....................................................................................................... 15 3.5.3 Mitchell-Fairweather格式 ................................................................................................... 15 3.5.4交替方向隐式格式的算法步骤 .......................................................................................... 16 4实例分析与结果分析 ................................................................................................................ 17 4.1算例 ......................................................................................................................................... 17 4.1.1已知有精确解的热传导问题 .............................................................................................. 17 4.1.2未知精确解的热传导问题 .................................................................................................. 19 4.2结果分析 ................................................................................................................................. 20 5稳定性探究与分析 .................................................................................................................... 21 5.1稳定性问题的提出 ................................................................................................................. 21

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

5.2 几种分析稳定性的方法 ........................................................................................................ 21 5.3 r变化对稳定性的探究 ........................................................................................................ 23 5.3.1 古典显式格式的稳定性 ..................................................................................................... 23 5.3.2 P-R格式格式的稳定性 ....................................................................................................... 24 结语 ............................................................................................................................................... 26 参考文献 ....................................................................................................................................... 27 附录 P-R格式的C++实现代码 ................................................................................................. 28 致谢 ............................................................................................................................................... 30

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

1绪论

1.1课题背景

抛物方程是一类特殊的偏微分方程，二维抛物方程的一般形式为

?u?Lu (1-1) ?t其中

L??u?u?u?u?u?u(a1(x,y,t))?(a2(x,y,t))?b1(x,y,t)?b2(x,y,t)?C(x,y,t) ?x?x?y?y?x?ya1?0,a2?0,C?0。[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]。

1

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

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

有限体积法(Finite Volume Methods)又称为控制体积法。其基本思路是：将计算区域划分为一系列互不重叠的控制体，并使每个网格点周围有一个控制体；将待求解的微分方程对每一个控制体积积分，便得出一组离散方程。该方法的未知量为网格点上的函数值。为了求出控制体积的积分，须假定函数值在网格点控制体边界上的变化规律。从积分区域的选取方法看来，有限体积法属于有限元方法中检验函数取分片常数插值的子区域法；从未知量的近似方法看来，有限体积法属于采用局部近似的离散方法[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?h u(t?h ( O h ) (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?1m 124124?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.2 Crank-Nicolson隐式格式

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

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

22由

22 L?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.3 Douglas差分格式[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.1 Peaceman-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.2 Rachford-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.3 Mitchell-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=0 y=0.1 y=0.2 y=0.3 y=0.4 x=0 0 0 0 0 0 x=0.1 0 -0.00013 -0.00025 -0.00035 -0.00041 x=0.2 0 -0.00025 -0.00048 -0.00066 -0.00077 x=0.3 0 -0.00035 -0.00066 -0.0009 -0.00106 x=0.4 0 -0.00041 -0.00077 -0.00106 -0.00125 x=0.5 0 -0.00043 -0.00081 -0.00112 -0.00131 x=0.6 0 -0.00041 -0.00077 -0.00106 -0.00125 x=0.7 0 -0.00035 -0.00066 -0.0009 -0.00106 x=0.8 0 -0.00025 -0.00048 -0.00066 -0.00077 x=0.9 0 -0.00013

-0.00025

-0.00035

-0.00041

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

误差 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的误差（t?0.1，h?1/20，r?h/k2?1/8）误差 y=0 y=0.1 y=0.2 y=0.3 y=0.4 x=0 0 0 0 0 0 x=0.1 0 8.86E-07 1.19E-06 1.47E-06 6.44E-07 x=0.2 0 1.19E-06 2.05E-06 2.19E-06 9.66E-07 x=0.3 0 1.47E-06 2.19E-06 1.48E-06 9.21E-07 x=0.4 0 6.44E-07 9.66E-07 9.21E-07 -8.6E-07 x=0.5 0 1.18E-06 9.93E-07 3.37E-09 -1.4E-06 x=0.6 0 6.44E-07 9.66E-07 9.21E-07 -8.6E-07 x=0.7 0 1.47E-06 2.19E-06 1.48E-06 9.21E-07 x=0.8 0 1.19E-06 2.05E-06 2.19E-06 9.66E-07 x=0.9 0 8.86E-07 1.19E-06 1.47E-06 6.44E-07 x=1 0 0

0 0 0

18

y=0.5 0 -0.00043 -0.00081 -0.00112 -0.00131 -0.00138 -0.00131 -0.00112 -0.00081 -0.00043

0

y=0.5 0 1.18E-06 9.93E-07 3.37E-09 -1.4E-06 -1.4E-06 -1.4E-06 3.37E-09 9.93E-07 1.18E-06 0

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

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

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

y=0.6 0 6.44E-07 9.66E-07 9.21E-07 -8.6E-07 -1.4E-06 -8.6E-07 9.21E-07 9.66E-07 6.44E-07 0 y=0.7 0 1.47E-06 2.19E-06 1.48E-06 9.21E-07 3.37E-09 9.21E-07 1.48E-06 2.19E-06 1.47E-06 0

y=0.8 0 1.19E-06 2.05E-06 2.19E-06 9.66E-07 9.93E-07 9.66E-07 2.19E-06 2.05E-06 1.19E-06 0

y=0.9 0 8.86E-07 1.19E-06 1.47E-06 6.44E-07 1.18E-06 6.44E-07 1.47E-06 1.19E-06 8.86E-07 0 y=1 0 0 0 0 0 0 0 0 0 0 0

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?1 Vn即对所有的?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.3 r变化对稳定性的探究

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/8 r=2/8 r=3/8 0 0 0 -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.00248 1.11E+27 -0.0009 -0.0018 -1E+27 -0.00048 -0.00095 -7.8E+26 -0.00013 -0.00026 -4.3E+26 0 0 0

r=5/8 0 9E+28 6.1E+29 9E+29 4.4E+29 1.3E+29 4.4E+29 1.7E+29 7.4E+29 9.3E+29 0

r=6/8 0 3.42E+32 4.72E+32 3.86E+32 2.39E+32 1.74E+32 2.39E+32 3.25E+32 2.58E+32 1.53E+32 0

r=7/8 0 4.17E+29 8.93E+29 9E+29 4.48E+29 4.84E+29 4.84E+29 8.16E+29 2.06E+29 5.41E+29 0

误差 x=0,y=0

x=0.1,y=0.1 x=0.2,y=0.2 x=0.3,y=0.3 x=0.4,y=0.4 x=0.5,y=0.5 x=0.6,y=0.6 x=0.7,y=0.7 x=0.8,y=0.8 x=0.9,y=0.9 x=1.0,y=1.0 r=4/8 0 -5.9E+29 -6.1E+29 -5E+29 -5.8E+29 -7.7E+29 -5E+29 -5.4E+29 -1E+29 -7.8E+29 0 r=1 0 4.48E+29 3.31E+28 5.95E+28 9.57E+29 8.93E+29 7.3E+29 1.54E+29 3.31E+28 8.16E+29

0

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

根据图5-1，表5-1可看出，当r?h/k2?1/4时稳定。否则，格式可能不稳定。且r越大，越容易发生震荡。

23

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

5.3.2 P-R格式格式的稳定性

采用P-R格式解例2，r由0变化到1时，不同的r时, u关于x，y的图像（t?0.1，

h??/20，r?h/k2?1/8）对比图如图5-2。采用P-R格式解例1，r由0变化到1时，不

同的r时, 数值解与解析解间误差（t?0.1，h?1/20，r?h/k2变化）对比如表5-2。

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

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

r=1/8 r=2/8 r=3/8 0 0 0 8.86E-07 1.19E-06 -0.000464 2.05E-06 1.19E-06 -0.000883 1.48E-06 9.21E-07 -0.001214 -8.6E-07 9.66E-07 -0.001424 -1.4E-06 1.47E-06 -0.001497 -8.6E-07 1.18E-06 -0.001424 1.48E-06 1.47E-06 -0.001214 2.05E-06 2.05E-06 -0.000883 8.86E-07 1.47E-06 -0.000464

0 0 0

误差

x=0,y=0

x=0.1,y=0.1 x=0.2,y=0.2 x=0.3,y=0.3 x=0.4,y=0.4 x=0.5,y=0.5 x=0.6,y=0.6 x=0.7,y=0.7 x=0.8,y=0.8 x=0.9,y=0.9 x=1.0,y=1.0 r=4/8 0 -0.000489 -0.000928 -0.001275 -0.001497 -0.001574 -0.001497 -0.001275 -0.000928 -0.000489

0

24

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

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

续表5-2 P-R格式解例1的误差（t?0.1，h?1/20）

r=5/8 r=6/8 r=7/8 0 0 0 -0.000464 -0.000396 -0.000288 -0.000883 -0.000753 -0.000548 -0.001214 -0.001034 -0.000753 -0.001424 -0.001214 -0.000883 -0.001497 -0.001275 -0.000928 -0.001424 -0.001214 -0.000694 -0.001214 -0.001034 -0.000753 -0.000883 -0.000753 -0.000548 -0.000464 -0.000396 -0.000288

0 0 0

r=1

0 -0.000152 -0.000288 -0.000396 -0.000464 -0.000489 -0.000464 -0.000396 -0.000288 -0.000152

0

由图5-2，表5-2可验证，对P-R格式而言，r取0～1中的任何数，该格式都是稳定的。

25

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

结语

本文主要讨论抛物方程的有限差分法，系统阐述了抛物方程的的几种数值解法，介绍了有限差分法的基本思想，步骤，和广泛的应用。着重讨论了抛物方程的有限差分解法，分析了这些方法的收敛性和稳定性。以热传导方程为例，编程实现了使用古典显式格式和P-R格式解二维抛物方程。通过数值算例，验证了这两种格式的稳定性条件，当步长比

r?h/k2?1/4，古典显式格式稳定，否则可能不稳定；对任意步长比r，P-R格式都是稳

定的。

进一步确定了自己的研究方向：

实现差分格式步长自适应调整，可以提高计算结果的准确性，在不知道准确解的情况下，得到较准确的结果。

探究新的稳定性好，计算简单的差分格式。并进一步学习有限元法，有限体积法，有限单元法，无网格方法等算法。

26

本文来源：https://www.bwwdw.com/article/wz56.html