第十章 偏微分方程数值解法

偏微分方程问题，其求解十分困难。除少数特殊情况外，绝

大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。

§1 差分方法的基本概念

1.1 几类偏微分方程的定解问题

椭圆型方程：其最典型、最简单的形式是泊松（Poisson）方程

uuu2f(x,y)2 xy特别地，当

22f(x,y)0时，即为拉普拉斯（Laplace）方程，又称

22为调和方程

uuu220 xyPoisson方程的第一边值问题为

uu22f(x,y)xyu(x,y)(x,y)(x,y) 其中

22(x,y)为以为边界的有界区域，为分段光滑曲线，

f(x,y)(x,y)分别为，上的已知连



称为定解区域，，续函数。

第二类和第三类边界条件可统一表示为

uu(x,y)(x,y)n 其中n为边界的外法线方向。当0时为第二类边界条件， 0时为第三类边界条件。

抛物型方程：其最简单的形式为一维热传导方程

uua0(a0)2 tx方程可以有两种不同类型的定解问题：

初值问题

22uua20xtu(x,0)(x)

初边值问题

t0,xxuutax200tT,0xlu(x,0)(x)0xlu(0,t)g(t),u(l,t)g(t)0tT12 其中

2(x)，g1(t)，g2(t)为已知函数，且满足连接条件

(0)g1(0),边界条件1件。

第二类和第三类边界条件为

(l)g2(0)

u(0,t)g(t),u(l,t)g2(t)称为第一类边界条

ux1(t)ux0g1(t)u

(t)ug(t)22xlx0tT

其中1，2。当

类边界条件，

否则称为第三类边界条件。

双曲型方程:

最简单形式为一阶双曲型方程

(t)0(t)01(t)2(t)0时，为第二

uua0 tx物理中常见的一维振动与波动问题可用二阶波动方程

22222uatux

描述，它是双曲型方程的典型形式。方程的初值问题为

22uu22a2txu(x,0)(x)ut0(x)tt0,xx

边界条件一般也有三类，最简单的初边值问题为

22uu2t2ax200tT,0xlu(x)0xlu(x,0)(x),tt0u(0,t)g(t),u(l,t)g(t)0tT12

1.2 差分方法的基本概念

差分方法又称为有限差分方法或网格法，是求偏微分方程定 解问题的数值解中应用最广泛的方法之一。

它的基本思想是：先对求解区域作网格剖分，将自变量的连 续变化区域用有限离散点（网格点）集代替；将问题中出现的连 续变量的函数用定义在网格点上离散变量的函数代替；通过用网 格点上函数的差商代替导数，将含连续变量的偏微分方程定解问 题化成只含有限个未知数的代数方程组（称为差分格式）。如果 差分格式有解，且当网格无限变小时其解收敛于原微分方程定解 问题的解，则差分格式的解就作为原问题的近似解（数值解）。 因此，用差分方法求偏微分方程定解问题一般需要解决以下问题： （1）选取网格；

（2）对微分方程及定解条件选择差分近似，列出差分格式； （3）求解差分格式；

（4）讨论差分格式解对于微分方程解的收敛性及误差估计。

下面，用一个简单的例子来说明用差分方法求解偏微分方程 问题的一般过程及差分方法的基本概念。

设有一阶双曲型方程初值问题。

uu0axtu(x,0)(x)t0,x

（1） 选取网格：

－2h -h 0 h 2h 3h 首先对定解区域分，最简单

常用一种网格是用两族分别平行于D{(x,t)x,t0}作网格剖x轴与t轴的等距直线 ，xxkkhttjj(k0,1,2,j0,1,2,)将

D分成许

多小矩形

区域。这些直线称为网格线，其交点称为网格点，也称为节点，

h和分别称作

x方向和t方向的步长。这种网格称为矩形网格。

（2） 对微分方程及定解条件选择差分近似，列出差分格式： 如果用向前差商表示一阶偏导数，即

u(xk1,tj)u(xk,tj)huu(xk1h,tj)x2x(xk,tj)h2u(xk,tj1)u(xk,tj)uut2(xk,tj2)t(xk,tj)2其中

01,21。

uua0方程 t x(xk,tj)处可表示为

在节点

u(xk,tj1)u(xk,tj)

au(xk1,tj)u(xk,tj)h

ah(xk1h,tj) ut(xk,tj2)ux2222

R(xk,tj)其中

(k0,1,2,,j0,1,2,)

u(xk,0)(xk)(k0,1,2,)。由于当

足够小时，在式

h,中略去

R(xk,tj)，就得到一个与方程相近似的差分方程

uk,j1uk,j

auk1,juk,jh0

此处，

uk,j可看作是问题的解在节点(xk,tj)处的近似值。同初值

条件

uk,0(xk)(k0,1,2,)

结合，就得到求问题的数值解的差分格式。 式

ahR(xk,tj)ut2(xk,tj2)ux2(xk1h,tj)22

O(h)称为差分方程的截断误差。

qpRO(h)，则称差分方 如果一个差分方程的截断误差为

程对

t是

q阶精度，对x是p阶精度的。显然，截断误差的阶数

越大，差分方程对微分方程的逼近越好。

若网格步长趋于0时，差分方程的截断误差也趋于0，则称 差分方程与相应的微分方程是相容的。这是用差分方法求解偏微 分方程问题的必要条件。

如果当网格步长趋于0时，差分格式的解收敛到相应微分方 程定解问题的解，则称这种差分格式是收敛的。

§2 椭圆型方程第一边值问题的差分解法

本节以Poisson方程为基本模型讨论第一边值问题的差分方法。 2.1 差分格式的建立

考虑Poisson方程的第一边值问题

2u2u22f(x,y)yxu(x,y)(x,y)(x,y)

(x,y)

取线

h,分别为

x方向和y方向的步长，如图所示，以两族平行

xxkkh，yyjj(k,j0,1,2,)将定

点

的

全

体

记

为

解区域剖分成矩形网 格。节

R{(xk,yj)xkkh,yjj,k,j为整数}。定解区

域内部的节点称为内点，记内点集

R为h。边界

与网格

线的交点称为边界点，边界点全体记为

h。与节点

(xk,yj)沿x方

向或

y方向只差一个步长的点

(xk1,yj)和(xk,yj1)称为节点

(xk,yj)的

相邻节点。如果一个内点的四个相邻节点均属于则内点，正内点的全体记为

，称为正



(1)，至少有一个相邻节点不属于

的内点称为非正则内点，非正则内点的全体记为求出第一边值问题在全体内点上的数值解。

为简便，记

(2)。问题是要

(k,j)(xk,yj)，

u(k,j)u(xk,yj)，

fk,jf(xk,yj)。对正则

(1)(k,j)内点，由二阶中心差商公式

ux22(k,j)u(k1,j)u(k,j)u(k,j)u(k1,j)2h(hhuxh1222uy2

(k,j)u(k1,j)2u(k,j)u(k1,j)h(4)ux4(x2h12u(k,j1)2u(k,j)u(k,j1)2(4)uy4(xk,yj212Poisson方程

uuf(x,y)22xy22 在点

(k,j)处可表示为

u(k1,j)2u(k,j)u(k1,j)u(k,j1)2u(k,j)u(k,22h

其中

h2(4)2(4)2R(k,j)ux4(xk1h,yj)ux4(xk,yj2)O(h1212

为其截断误差表示式，略去

R(k,j)，即得与方程相近似的差分方程

uk1,j2uk,juk1,jh2uk,j12uk,juk,j12fk,j

式中方程的个数等于正则内点的个数，而未知数正则内点处解

uk,j则除了包含

u的近似值外，还包含一些非正则内点处u的近

u值作为该点上u似值，因而方程个数少于未知数个数。在非正则内点处Poisson 方程的差分近似不能按上式给出，需要利用边界条件得到。 边界条件的处理可以有各种方案，下面介绍较简单的两种。 （1）直接转移

用最接近非正则内点的边界点上的近似，这就是边界条件的直接转移。

例如，点P(k,j)为非正则内点，其最接近的边界点为

值的

Q点，则有

uk,ju(Q)(Q)求出，其截断误差为O(h个数相等。

（2）线性插值

(k,j)

上式可以看作是用零次插值得到非正则内点处u的近似值，容易

)。将上式代入，方程个数即与未知数

(2)这种方案是通过用同一条网格线上与点内点作线性插值得到非正则内点P(k,P相邻的边界点与

值的近似。由点

j)处uR与

T的线性插值确定

u(P)的近似值uk,j，得

u k,j其中dRPhd(R)(T)

hdhd2O(h)。将其与方程相近似的差分，其截断误差为

方

程联立，得到方程个数与未知数个数相等的方程组，求解此方程 组可得Poisson方程第一边值问题的数值解。

上面所给出的差分格式称为五点菱形格式，

uk1,j2uk,juk1,jh

实际计算时经常取

2uk,j12uk,juk,j12fk,jh，此时五点菱形格式可化为

1(uk1,juk1,juk,j1uk,j14uk,j)fk,j 2h简记为

12hufk,jk,j ◇

其中◇

uk,juk1,juk1,juk,j1uk,j14uk,j。

2u2u(x,y)22f(x,y)yx u(x,y)lg[(1x)2y2]例1 用五点菱形格式求解拉普拉斯（Laplace）方程第一边值问题

1}。取h。 其中{(x,y)0x,y13

（0,0） (1,0) (2,0) (3,0) [解]网格中有四个内点，均为正则内点。由五点菱形格式，得方 程组

1h21h21h212h(u2,1u0,1u1,2u1,04u1,1)0(u3,1u1,1u2,2u2,04u2,1)0(u2,2u0,2u1,3u1,14u1,2)0(u3,2u1,2u2,3u2,14u2,2)0 代入边界条件

16u1,0lg9,ulg10,0,19ulg25,1,3937u3,1lg,9其解为 u1,125u2,0lg913u0,2lg934u2,3lg940u3,2lg9，

0.2756919u2,10.4603488，

u1,20.3467842当

，u2,20.5080467时

，

对

h1(uk1,juk1,juk,j1uk,j14uk,j)fk,j2 h利用点(k,式，称为五点 矩形格式，简记为

j)，(k1,j1)，(k1,j1)构造的差分格

其

12h2□

uk,jfk,j

□

中

uk,juk1,j1uk1,j1uk1,j1uk1,j14uk,j，其

截断误差为

huuu4R(k,j)6O(h)O424412xxyy(k,j)

2O(h)，称它们具有 五点菱形格式与矩形格式的截断误差均为

2444二阶精度。如果用更多的点构造差分格式，其截断误差的阶数可

以提高，如利用菱形格式及矩形格式所涉及的所有节点构造出的 九点格式就是具有四阶精度的差分格式。

§3 抛物型方程的差分解法

以一维热传导方程

uua20tx2(a0)

为基本模型讨论适用于抛物型方程定解问题的几种差分格式。

3.1 差分格式的建立

首先对

xt平面进行网格剖分。分别取h,为x方向与

t方向

的步长，用两族平行直线

xxkkh(k0,1,2,)，

平面剖分成矩形网格，节点

ttjj(j0,1,2)，将xt为(xk,tj)(k0,1,2,,j0,1,2)。为

，

简便，记

(k,j)(xk,tj)，

u(k,j)u(xk,tj)，，

k(xk)g1jg1(tj)g2jg2(tj),1j1(tj),2j2(tj)。

（一）微分方程的差分近似

在网格内点

(k,j)处，对

(k,j)ut分别采用向前、向后及中心差商公式

utututu(k,j1)u(k,j)u(k,j)u(k,j1)O()O()

(k,j)(k,j)u(k,j1)u(k,j1)O(2)2一维热传导方程

u2ua02tx可分别表示为

(a0)

u(k,j1)u(k,j)

u(k1,j)2u(k,j)u(k1,j)aO(2hu(k,j1)u(k,j1)u(k1,j)2u(k,j)u(k1,j)aO(2h

u(k,j1)u(k,j1)u(k1,j)2u(k,j)u(k1,j)22aO(h22h

由此得到一维热传导方程的不同差分近似

uk,j1uk,j

auk1,j2uk,juk1,jh20

uk,juk,j1

auk1,j2uk,juk1,jhh20

uk,j1uk,j12auk1,j2uk,juk1,j20

2O(h)，上述差分方程所用到的节点各不相同。其截断误差分别为

O(h2)和O(2h2)。因此，它们都与一维热传导方程相容。

如果将式

u(k,j1)u(k,j1)u(k1,j)2u(k,j)u(k1,j)2aO(2h2

1(uk,j1uk,j1)代替，则可得到又一种差分近似

中的uk,j用

2

uk,j1uk,j12auk1,juk,j1uk,j1uk1,jh20

差分方程用到四个节点。由Taylor公式容易得出

uk,j1(uk,j1uk,j1)O(2)

2O(22h2)Oh2故其的截断误差为

。因而不是对任意的h,0，

此差分方程都能逼近热传导方程

u2ua02tx仅当

(a0)，

o(h)时，才成立。

综上可知，用不同的差商公式可以得到微分方程的不同的差

分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳 定性。前面三个方程都具有相容性，而此方程则要在一定条件下 才具有相容性。

（二）初、边值条件的处理

为用差分方法求解定解问题初值问题

2uua20xtu(x,0)(x)t0,xx

初边值问题

u2utax200tT,0xlu(x,0)(x)0xlu(0,t)g(t),u(l,t)g(t)0tT12

还需对定解条件进行离散化。

对初始条件及第一类边界条件，可直接得到

uk,0u(xk,0)k或 k0,1,

(k0,1,,n)

u0,ju(0,tj)g1j

un,ju(l,tj)g2j

(j0,1,,m1)

lTn,m其中

h对第二、三类边界条件

ux1(t)uu2(t)ux0tT

x0g1(t)g2(t)

xl需用差分近似。下面介绍两种较简单的处理方法。 （1）

在左边界(x在右边界

0)处用向前差商近似偏导数

uxux，

(xl)处用向后差商近似

uxux，即

(0,j)u(1,j)u(0,j)O(h)hu(n,j)u(n1,j)O(h)h

(n,j)(j0,1,,m)

则得边界条件的差分近似为

u1,ju n,j其截断误差为

u0,jhun1,jh1ju0,jg1j2jun,jg2j

(j0,1,,m)

O(h)。

ux,即

（2）用中心差商近似

uu(1,j)u(1,j)O(h2)x(0,j)2huu(n1,j)u(n1,j)O(h2)x(n,j)2h

(j0,1,,m)

则得边界条件的差分近似为

u1,ju1,j1ju0,jg1j2hun1,jun1,ju2jn,jg2jh

(j0,1,,m)

2O(h)。误差的阶数提高了，但出现定解区域外的 其截断误差为

,节点(1j)和(n1,j)，这就需要将解拓展到定解区域外。可以通

过用内节点上的

u值插值求出

u1,j和un1,j，也可以假定热传导方

程在边界上也成立，将差分方程扩展到边界节点上，由此消去u1,j

和

un1,j。

（三）几种常用的差分格式

以热传导方程的初边值问题

u2utax200tT,0xlu(x,0)(x)0xlu(0,t)g(t),u(l,t)g(t)0tT12

为例给出几种常用的差分格式。

（1）古典显式格式

令

ar2h，则

uk,j1uk,j可

auk1,j2uk,juk1,jh改

20

成

写

uk,j1ruk1,j(12r)uk,jruk1,j

将其与初始条件及第一类边界条件

uk,0u(xk,0)k或 k0,1,

(k0,1,,n)

u0,ju(0,tj)g1j

un,ju(l,tj)g2j

(j0,1,,m1)

结合，我们得到求解此问题的一种差分格式

uk,j1ruk1,j(12r)uk,jruk1,j(k1,2,,n1,j0,1,,m(k0,1,,n)uk,0kug,ug(j0,1,,m)1jn,j2j0,j

由于第0层

(j0)上节点处的u值已知(uk,0k)，由此即可算出

u在

第一层

(j1)上节点处的近似值uk,1。重复使用此式，可以逐层计

算出所有的

uk,j，因此此差分格式称为古典显式格式。又因式中

式

只出现相邻两个时间层的节点，故此式是二层显式格式。 （2）古典隐式格式

将

uk,juk,j1auk1,j2uk,juk1,jh20

整理并与初始条件及第一类边界条件式联立，得差分格式如下

uk,j1uk,jr(uk1,j12uk,j1uk1,j1)(k1,2,,n1,j0,(k0,1,,n)uk,0kug,ug(j0,1,,m)0,j1jn,j2j

ar2其中

h。虽然第0层上的

u值仍为已知，但不能由上式直

接计算以上各层节点上的值

uk,j，必须通过解下列线性方程组

uk,j1r(uk1,j12uk,j1uk1,j1)ruk1,j1(12r)uk,j1ruk1,j

r12rr12r0

r0r12rrr12r0u1,j1u1,jrguu2,j12,jun2,un2,j1un1,j1un1,jrg (j0,1,,m1)

才能由

uk,j计算uk,j1，故此差分格式称为古典隐式格式。此方程

组是三对角方程组，且系数矩阵严格对角占优，故解存在唯一。

（3）Richardson格式

Richardson格式是将式

uk,j1uk,j12

auk1,j2uk,juk1,jh20整理后与初始条件及第一类边界条件式联立。其计算公式为

uk,j1uk,j12r(uk1,j2uk,juk1,j)(k1,2,,n1,(k0,1,,n)uk,0kug,ug(j0,1,,m)0,j1jn,j2j

这种差分格式中所涉及的节点出现在

j1,j,j1三层上，故为三层

显式格式。Richardson格式是一种完全不稳定的差分格式，因此它在实际计算中是不能采用的。

（4）杜福特-弗兰克尔（DoFort-Frankel）格式

DoFort-Frankel格式也是三层显式格式，它是由式

uk,j1uk,j12

auk1,juk,j1uk,j1uk1,jh2与初始条件及第一类边界条件式结合得到的。具体形式如下：

2r12ruk,j112r(uk1,juk1,j)12ruk,j1(k1,2,,n(k0,1,,n)uk,0kug,ug(j0,1,,m)0,j1jn,j2j

用这种格式求解时，除了第0层上的值

uk,0由初值条件得到，必

须先用二层格式求出第1层上的值

uk,1，然后再按上式逐层计

算

uk,j(j2,3,,m)。

（5）六点隐式格式

对二阶中心差商公式

2u2x

(k,j)u(k1,j)2u(k,j)u(k1,j)2hu如果用

x2均值

2在点

(k,j1)与点(k,j)处的二阶中心差商的平

近似

2ux21(k,j)处的值，即 在

2ux2

21(k,j)2uk1,j12uk,j1,uk1,j1uk1,j2uk,ju2h2同时

ut1(k,j)处的值也用中心差商近似，即 在点

21(k,j)2这样又得到热传导方程的一种差分近似

u xuk,j1uk,j

uk,j1uk,j

a2(uk1,j12uk,j1,uk1,j1uk1,j2uk,2h22O(h)，将上式与初始条件及第一类边界条件式 其截断误差为

联立并整理，得差分格式

rruk,j1uk,j2(uk1,j12uk,j1uk1,j1)2(uk1,j2uk,juk(k1,2,,n1,j0,1,,m1)uk,0k(k0,1,,n)(j0,1,,m)u0,jg1j,un,jg2j

此格式涉及到六个节点，它又是隐式格式，故称为六点隐式格式。与古典隐式格式类似，用六点格式由第

j层的值uk,j计算第j1

层的值uk,j1时，需求解三对角方程组

1rr/20r/21rr/200r/21rr/2r/21r

(j0,1,,m1)

此方程组的系数矩阵严格对角占优，故仍可用追赶法求解。

例2 用古典显式格式求初边值问题

u2u00t3,0x3tx2u(x,0)x20x3 u(0,t)0,u(3,t)90t3 的数值解，取

h1,0.5。

[解] 这里

a1,rah20.5，

(x)x2，

u1,j1u2,j1un2,j1un1,j1g1(t)0,g2(t)9。由格式

uk,j1ruk1,j(12r)uk,jruk1,j(k1,2,,n1,(k0,1,,n)uk,0kug,ug(j0,1,,m)1jn,j2j0,j

可得到

uk,j10.5(uk1,juk1,j)22uk,0xkku0,u93,j0,j 将初值

(k1,2,j0,1,,5)(k0,1,2,3)(j0,1,,6)uk,0代入上式，即可算出

u1,10.5(u2,0u0,0)0.5(40)2 u2,10.5(u3,0u1,0)0.5(91)5

将边界条件

u0,10，u3,19及上述结果代入又可求得

u0,20,u1,22.5,u2,25.5,u3,29

如此逐层计算，得全部节点上的数值解为

u0,30,u1,32.75,u2,35.75,u3,39,u0,40,

§4 双曲型方程的差分解法

对二阶波动方程

22uu2a22tx

uv如果令1t组

uv2，x，则方程可化成一阶线性双曲型方程

v12v2taxv2v1 xt记

v(v1,v2)T，则方程组可表成矩阵形式

2v0avvA t xx10矩阵

A有两个不同的特征值a，故存在非奇异矩阵P，使得

a0PAP

0a1作变换

wPv(w1,w2)T，方程组可化为

wwtx

方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨 论一阶双曲型方程的差分解法。

4.1几种简单的差分格式

考虑一阶双曲型方程的初值问题

uua0xtu(x,0)(x) 将网

t0,xxxt平面剖分成矩形网格，取x方向步长为h,t方向步长为，

格

线

为

xxkkh(k0,1,2,)，

ttjj(j0,1,2)。为简便，记

(k,j)(xk,tj)，u(k,j)u(xk,tj),k(xk)。

以不同的差商近似偏导数，可以得到方程的不同的差分近似

uk,j1uk,j

huk,j1uk,juk,juk1,ja0

huk,j1uk,j

auk1,juk,j0

auk1,juk1,j2h0

2O(h)O(h)O(h)。结合离散截断误差分别为，与

化的初始条件，可以得到几种简单的差分格式

uk,j1uk,jar(uk1,juk,j)uk,0k

(k0,1,2,,j0,1,2,)uk,j1uk,jar(uk,juk1,j)uk,0k

(k0,1,2,,j0,1,2,)aruk,j1uk,j(uk1,juk1,j)(k0,1,2,,j0,1,2,)2uk,0k

其中

rh。如果已知第

j层节点上的值

uk,j，按上面三种格式就

可求出第

j1层上的值uk,j1。因此，这三种格式都是显式格式。

ut 如果对

u采用向后差商，x采用向前差商，则方程可化成

u(k,j)u(k,j1)u(k1,j)u(k,j)aO(h)0h

相应的差分格式为

uk,j1uk,jar(uk1,j1uk,j1)uk,0k

(k0,1,2,,j0,1,2,)此差分格式是一种隐式格式，必须通过解方程组才能由第

j层节

uu点上的值k,j求出第j1层节点上的值k,j1。

例对初值问题

uu0txu(x,0)(x)t0,xx 其中

1x01(x)x0 20x0用差分格式求其数值解

uk,j(j1,2,3,4)，取

1rh2 [解] 记

。

xkkh(k0,1,2,)，由初始条件

k1k1,2,1(xk)k0 20k1,2,按差分格式

uk,j1uk,jar(uk1,juk,j)uk,0k

(k0,1,2,,j0,1,2计算公式为

uk,j131uk,juk1,j

22计算结果略。

如果用差分格式

uk,j1uk,jar(uk,juk1,j)uk,0k

求解，计算公式为

(k0,1,2,,j0,1,2uk,j1计算结果见表略。

1(uk,juk1,j)

21xt1u(x,t)xt2与准确解 0x0比较知，按前一个差分格式所求得的数值解不收敛到初值问题

的解，而后一个差分格式的解收敛到原问题的解。