2024年2月8日发(作者：)

WORD格式可编辑

第9章 偏微分方程的差分方法

含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。

9.1椭圆型方程边值问题的差分方法

9.1.1 差分方程的建立

最典型的椭圆型方程是Poisson（泊松）方程

2u2u

u(22)f(x,y),(x,y)G （9.1）

xyG是x,y平面上的有界区域，其边界Γ为分段光滑的闭曲线。当f(x,y)≡0时，方程（9.1）称为Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件

第一边值条件

u(x,y) （9.2）

第二边值条件

u(x,y) （9.3）

nuku)(x,y) （9.4）

n 第三边值条件

(这里，n表示Γ上单位外法向，α(x,y),β(x,y),γ(x,y)和k(x,y)都是已知的函数，k(x,y)≥0。满足方程（9.1）和上述三种边值条件之一的光滑函数u(x,y)称为椭圆型方程边值问题的解。

用差分方法求解偏微分方程，就是要求出精确解u(x,y)在区域G的一些离散节点（xi，yi）上的近似值ui,j≈(xi，yi）。差分方法的基本思想是，对求解区域G做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。

设G={0

x=ih1, i=0,1,„,N1, h1=a/N1

y=jh2, j=0,1,„,N2, h2=b/N2

将G剖分为网格区域，见图9-1。h1，h2分别称为x方向和y方向的剖分步长，网格交点（xi，yi）称为剖分节点（区域内节点集合记为Gh={（xi，yi）; （xi，yi）∈G}）,网格线与边界Γ的交点称为边界点，边界点集合记为Γh。

专业知识整理分享

WORD格式可编辑

现在将微分方程（9.1）在每一个内节点（xi，yi）上进行离散。在节点（xi，yi）处，方程（9.1）为

2u2u[2(xi,yi)2(xi,yi)]f(xi,yi),(xi,yi)Gh （9.5）

xy需进一步离散（9.5）中的二阶偏导数。为简化记号，简记节点（xi，yi）=(i,j)，节点函数值u（xi，yi）=u(i,j)。利用一元函数的Taylor展开公式，推得二阶偏导数的差商表达式

2u1(i,j)[u(i1,j)2u(i,j)u(i1,j)]0(h12)22xh1u12(i,j)[u(i,j1)2u(i,j)u(i,j1)]0(h2)22yh2代入（9.5）式中，得到方程（9.1）在节点(i,j)处的离散形式

2

11[u(i1,j)2u(i,j)u(i1,j)][u(i,j1)2u(i,j)u(i,j1)]22h2

h1

2fi,j0(h12h2),(i,j)Gh2其中fi,jf(xi,yi)。舍去高阶小项0(h12h2)，就导出了u(i,j)的近似值ui,j所满足的差分方程

11

[u2uu][ui,j12ui,jui,j1]fi,j,(i,j)Gh （9.6）i1,ji,ji1,j22h1h222在节点(i,j)处方程（9.6）逼近偏微分方程（9.1）的误差为O(h1h2)，它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。

在差分方程（9.6）中，每一个节点(i,j)处的方程仅涉及五个节点未知量ui,j，ui+1,j，ui-1,j，ui,j+1，ui,j-1，因此通常称（9.6）式为五点差分格式，当h1= h2=h时，它简化为

专业知识整理分享

WORD格式可编辑

1[ui1,jui1,jui,j1ui,j14ui,j]fi,j,(i,j)Gh

2h差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值ui,j ,(i,j)∈Gh外，还包括边界点值。例如，点(1,j)处方程就含有边界点未知量u0,j。因此，还要利用给定的边值条件补充上边界点未知量的方程。

对于第一边值条件式（9.2），可直接取ui,j=α（xi，yi）， (i,j)∈Γh （9.7）

对于第三（k=0时为第二）边值条件式（9.4），

以左边界点(1,j)为例，见图9-2,

利用一阶差商公式

uu(0,j)u(1,j)(0,j)O(h1)

nh1则得到边界点(0,j)处的差分方程

u0,ju1,jh1k0,ju0,jr0,j （9.8）

联立差分方程（9.6）与（9.7）或（9.8）就形成了求解Poisson方程边值问题的差分方程组，它实质上是一个关于未知量{ui,j}的线性代数方程组，可采用第2，3章介绍的方法进行求解。这个方程组的解就称为偏微分方程的差分近似解，简称差分解。

考虑更一般形式的二阶椭圆型方程

[uuuu(A)(B)CDEu]f(x,y),(x,y)G （9.9）

xxyyxy其中A(x,y)≥Amin>0, B(x,y) ≥Bmin >0, E(x,y) ≥0。引进半节点xi121xih1,

2yi121yih2,利用一阶中心差商公式，在节点（i,j）处可有

2u1u1u1(A)(i,j)[(A)(i,j)(A)(i,j)]O(h12)xxh1x2x21u(i1,j)u(i,j)u(i,j)u(i1,j)[A1A1]O(h12)i,jh1i2,jh1h12uu(i1,j)u(i1,j)(i,j)O(h12)x2h1 专业知识整理分享

WORD格式可编辑

对uu类似处理，就可推得求解方程（9.9）的差分方程

(B),yyy[ai1,jui1,jai1,jui1,jai,j1ui,j1ai,j1ui,j1ai,jui,j]f(i,j),(i,j)Gh其中

（9.10）

h12ah(ACi,j)11i1,ji,j22h12ah(ACi,j)11i1,ji,j22h2 （9.11）

(B12Di,j)ai,j1h2i,j22h22ah(BDi,j)i,j121i,j2222ai,jh1(A1A1)h2(B1B1)Ei,ji,ji,ji,ji,j2222显然，当系数函数A(x,y)=B(x,y)=1, C(x,y)=D(x,y)=E(x,y)=0时，椭圆型方程（9.9）就成为Poisson方程（9.1），而差分方程（9.10）就成为差分方程（9.6）。容易看2出，差分方程（9.10）的截断误差为O(h12h2)阶。

9.1.2 一般区域的边界条件处理

前面已假设G为矩形区域，现在考虑G为一般区域情形，这里主要涉及边界条件的处理。

考虑Poisson方程第一边值问题

uf(x,y),(x,y)G （9.12）

u(x,y),(x,y)其中G可为平面上一般区域，例如为曲边区域。仍然用两组平行直线：x=x0+ih1,y=y0+jh2,i,j=0,±1,„,对区域G进行矩形网格剖分，见图9-3。

如果一个内节点（i,j）的四个相邻节点（i+1,j），（i-1,j），（i,j+1）和（i,j-1）属于GG，则称其为正则内点，见图9-3中打“。”号者；如果一个节点（i,j）属于G且不为正则内点，则称其为非正则内点，见图9-3中打“.”号者。记正则，非正则内点集合为h。显然，当G为矩形区域时，内点集合为Gh 专业知识整理分享

WORD格式可编辑

Gh,hh成立。

Gh

在正则内点（i,j）处，完全同矩形区域情形，可建立五点差分格式

11 （9.13）[u2uu][ui,j12ui,jui,j1]fi,j,(i,j)Gh

i1,ji,ji1,j22h1h2在方程（9.13）中，当（i,j）点临近边界时，将出现非正则内点上的未知量，因此必须补充非正则内点处的方程。

若非正则内点恰好是边界点，如图9-4中

D点，则利用边界条件可取uD=α(D)

对于不是边界点的非正则内点，如图9-4中B点，一般可采用如下两种处理方法。

a.直接转移法.取与点B距离最近的边界点（如图9-4中E点）上的u的值作为u(B)的近似值uB，即uB=u(E)=α(E)

直接转移法的优点是简单易行，但精度较低，只为一阶近似。

b.线性插值法.取B点的两个相邻点（如图9-4中边界点A和正则内点C作为插值节点对u(B)进行线性插值

u(B)xCxBxxAu(A)Bu(C)O(h12)

xCxAxCxA则得到点B处的方程

uBh1(A)uC,xBxA

h1h1线性插值法精度较高，为二阶近似。

对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般 专业知识整理分享

WORD格式可编辑

区域上边值问题（9.12）的差分近似解。

对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。

第二、三边值条件的处理较为复杂，这里不再讨论。

9.2 抛物型方程的差分方法

本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。

9.2.1 一维问题

作为模型，考虑一维热传导的初边值问题

u2ua2f(x,t),0xl,0tT （9.14）

txu(x,0)(x),0xl （9.15）

u(0,t)g1(t),u(l,t)g2(t),其中a是正常数，f(x,t),0tT （9.16）

(x),g1(t)和g2(t)都是已知的连续的函数。

现在讨论求解问题（9.14）-(9.18)的差分方法。首先对求解区域G={0≤x≤l,

0≤t≤T}进行网格剖分。取空间步长h=l/N,时间步长τ=T/M,其中N,M是正整数，作两族平行直线

xxjjh,j0,1,,Nttkk,k0,1,M

将区域G剖分成矩形网格，见图9-5，网格交点（xj，tk）称为节点。

用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解u(x,t)在每个节点（xj，tk）处的近似值uju(xj,tk)。为简化记号，简记节点（xj，tk）=u(j,k)。

利用一元函数的Taylor展开公式，可推出下列差商表达式

专业知识整理分享

k

WORD格式可编辑

uu(j,k1)u(j,k)(j,k)O() （9.17）

tuu(j,k)u(j,k1)(j,k)O() （9.18）

tuu(j,k1)u(j,k1)(j,k)O(2) （9.19）

t2

2uu(j1,k)2u(j,k)u(j1,k)2

2(j,k)O(h) （9.20）

2xh1.古典显格式

在区域G的内节点(j,k)处，利用公式（9.17）和（9.20），可将偏微分方程（9.14）离散为

u(j,k1)u(j,k)kau(j1,k)2u(j,k)u(j1,k)k2fO(h)

j2hk其中fjf(xi,tk)。舍去高阶小项O(h2)，就得到节点近似值（差分解）uj所满足的差分方程

1ukukjjakkukj12ujuj1h2fjk （9.21）

显然，在节点(j,k)处，差分方程（9.21）逼近偏微分方程（9.14）的误差为O(h2)，这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。现将（9.21）式改写为便于计算的形式，并利用初边值条件（9.15）与（9.16）补充上初始值和边界点方程，则得到

1kkkukrukjj1(12r)ujruj1fjj1,2,,N1,k0,1,,M1 （9.22）

0uj(xj),j1,2,,N1ukg(t),ukg(t),k0,1,M1kN2k0其中ra2称为网比。

h与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程（9.22），当第k层节点值{uj}已知时，可直接计算出第k+1层节点值{uj}。这样，从第0层已知值uj(xi)开始，就可逐层求出各时间层的节点值。差分方程（9.22）的求解计算是显式的，无须求解方程组，故称为古典显格式。此外，在式（9.22）中，每个内节点处方程仅涉及k和k+1两层节点值，称这样的差分格式为 专业知识整理分享

0kk1

WORD格式可编辑

双层格式。

差分方程（9.22）可表示为矩阵形式

uk1AukFk,k0,1,,M1

0 （9.23）

u

其中

r12rr12rr



Arr12rkTuk(u1k,,uN1)

F(f1rg1(tk),f2,,fN2,fN1rg2(tk))

kkkkkT((x1),,(xN1))T

2. 古典隐格式

在区域G的内节点(j,k)处，利用公式（9.18）和（9.20），可将偏微分方程（9.14）离散为

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

2h舍去高阶小项O(h)，则得到如下差分方程

k1ukjujkkukj12ujuj1

ah22fjk （9.24）

它的截断误差为O(h)，逼近精度与古典显格式相同。改写（9.24）式为便于计算的形式，并补充上初始值与边界点方程，则得到

kkk1rukfjkj1(12r)ujruj1ujj1,2,,N1,k0,1,,M

0 （9.25）

uj(xj),j1,2,,N1ukg(t),ukg(t),k0,1,M1kN2k0

专业知识整理分享

WORD格式可编辑

与古典显格式不同，在差分方程（9.25）的求解中，当第k-1层值{uj}已知时，必须通过求解一个线性方程组才能求出第k层值{uj}，所以称（9.25）式为古典隐格式，它也是双层格式。

差分方程（9.25）的矩阵形式为

kk1Bukuk1Fk,k1,2,M （9.26）

0u其中

r12rr12r

B

kkr

rr12r向量u,F,同（9.23）式中定义。从（9.26）式看到，古典隐格式在每一层计算时，都需求解一个三对角形线性方程组，可采用追赶法求解。

-Nicolson格式（六点对称格式）

利用一元函数Taylor展开公式可得到如下等式

u1u(j,k1)u(j,k)(j,k)O(2)t2

2

22u11uu(j,k)[(j,k)(j,k1)]O(2)22222xxx使用这两个公式，在(j,k)点离散偏微分方程（9.14），然后利用（9.20）式进一步离散二阶偏导数，则可导出差分方程

1ukukjjk1k1k1kkk1k1uj12ujuj1uj12ujuj1a[]fj2 （9.27）

222hh12其截断误差为O(h)，在时间方向的逼近阶较显格式和隐格式高出一阶。这个差分格式称为Crank-Nicolson格式，有时也称为六点对称格式，它显然是双层隐式格式。改写（9.27）式，并补充初始值和边界点方程得到

22 专业知识整理分享

WORD格式可编辑

1k11rukrukj12(1r)ujj11kkkk2ru2(1r)uruj1jj12fj （9.28）

j1,2,,N1,k0,1,,M10uj(xj),j1,2,,N1ukg(t),ukg(t),k0,1,M1kN2k0它的矩阵形式为

1kk1k(IB)u(IA)uF2,k1,2,M （9.29）

0u,k0,1,,M1在每层计算时，仍需求解一个三对角形方程组。

4. Richardson格式

利用公式（9.19）和（9.20），可导出另一个截断误差为O(2h2)阶的差分方程

11ukukjj2akkukj12ujuj1h2fjk

称之为Richardson格式。可改写为

11kkukuk2r(ukfjk (9.32)

jjj12ujuj1)2这是一个三层显式差分格式。在逐层计算时，需用到uj和uj两层值才能得到k+1层值uj。这样，从第0层已知值uj(xj)开始，还须补充上第一层值uj，才能逐层计算下去。可采用前述的双层格式计算uj。

除上述四种差分格式外，还可构造出许多逼近偏微分方程（9.14）的差分格式，但并不是每个差分格式都是可用的。一个有实用价值的差分格式应具有如下性质：

（1）收敛性。对任意固定的节点（xj，tk）,当剖分步长,h0时，差分解ujkk1kk1011应收敛到精确解u（xj，tk）。

（2）稳定性。当某一时间层计算产生误差时，在以后各层的计算中，这些误差的传播积累是可控制而不是无限增长的。

理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。因此，这里主要研究差分格式的稳定性。

专业知识整理分享

WORD格式可编辑

作为例子，先考查Richardson格式的稳定性。设uj是当计算过程中带有误差kk时，按Richardson格式（9.30）得到的实际计算值。误差ekuj是理论值，jujuj。kkk假定右端项fj的计算是精确的，网比r1，则ekj满足

211kk

ekek(ekjjj12ejej1) （9.31）

设前k-1层计算时精确的，误差只是在第k层j0点发生，即

1kek0,ekjj0,ej0,jj0。

则利用（9.31）式可得到误差的传播情况，见表9-1。

表9-1 r=1/2时Richardson格式的误差传播

j

k

k

k+1

k+2

k+3

k+4

k+5

k+6

j0-4

0

0

0

0

ε

-10ε

71ε

j0-3

0

0

0

ε

-8ε

49ε

-260ε

j0-2

0

0

ε

-6ε

31ε

-144ε

j0-1

0

ε

-4ε

17ε

-68ε

273ε

j0

ε

-2ε

7ε

-24ε

89ε

-388ε

j0+1

0

ε

-4ε

17ε

-68ε

273ε

j0+2

0

0

ε

-6ε

31ε

-144ε

j0+3

0

0

0

ε

-8ε

49ε

-260ε

j0+4

0

0

0

0

ε

-10ε

71ε 641ε -1096ε 1311ε -1096ε 641ε

从表中看出，误差是逐层无限增长的。表中的计算虽然是就网比r1进行的，2实际上对任何r>0都会产生类似现象，所以Richardson格式是不稳定的。

利用误差传播图表方法考查差分格式的稳定性虽然直观明了，但只能就具体取定的r值进行，并且也不适用于隐式差分格式。

9.2.2 差分格式的稳定性

前节构造的几种双层差分格式都可以表示为如下的矩阵方程形式

ukHuk1Fk （9.32）

0u其中H称为传播矩阵。对于显格式H=A, 隐格式H=B,六点对称格式H=(I+B)

(I+A)。一般的三层格式也可以转化为双层格式。

为了讨论方便，设在初始层产生误差，且假定右端项F的计算是精确的。用u 专业知识整理分享

0k-1 -1k

WORD格式可编辑

表示当初始层存在误差0时，由差分格式（9.32）得到的计算解，则uk满足方程

ukHuk1Fk （9.33）

00u记误差向量uu，则k满足方程

kkkkHk1,k1,2, （9.34）

0为初始误差定义9.1 称差分格式（9.32）是稳定的，如果对任意初始误差0，误差向量k在某种范数下满足

k0

C,k0,00 （9.35）

其中C为与h,τ无关的常数。

这个定义表明，当差分格式稳定时，它的误差传播是可控制的。

从（9.34）式递推得到

kHk0,0k因此，差分格式稳定的充分必要条件是

T

HkC,0kT （9.36）

定理9.3 （稳定性必要条件）差分格式稳定的必要条件是存在与h,τ无关的常数M,使谱半径

(H)1M （9.37）

**定理9.4 （稳定性充分条件）设H为正规矩阵，即HHHH,则（9.37）式也是差分格式稳定的充分条件。

下面讨论几种差分格式的稳定性。为便于讨论，引进N-1阶矩阵

01101

S10110这个特殊矩阵的特征值为

专业知识整理分享

WORD格式可编辑

j2cosSj,Nj1,2,,N1 （9.38）

例9-1古典显格式 此时H=A=(1-2r)I+rS。 利用（9.38）式和三角函数公式，可求得H的特征值为

2

j14rsin(j),2Nj1,2,,N1

1。由于H=A为实对称矩阵，所2为使稳定性条件式（9.39）成立，必须且只须r以古典显格式稳定的充分必要条件是网比r1

2例9-2 古典隐格式 此时H=B-1,B=(1+2r)I-rS。利用（9.38）式可求得H的特征值为

2

j[14rsin(j1)],2Nj1,2,,N1

-1显然，对任意r>0，条件（9.37）成立。注意，H=B仍为实对称矩阵，所以古典隐格式对任何网比r>0都是稳定的，称为绝对稳定。

例9-3 六点对称格式 此时H=(I+B)-1(I+A)，利用矩阵A和B的特征值可得到矩阵H的特征值为

j)2N

j,j24rsin2()2N24rsin2(j1,2,,N1

则对任意r>0,条件（9.37）成立。由于A和B均为实对称矩阵，且AB=BA，则可验证H也是实对称矩阵。所以六点对称格式是绝对稳定的。

习题九

9-1 试用五点差分格式求解Poisson 方程的边值问题

2u2u2216,(x,y)G

xy(x,y)u|0,其中G={-1

9-2 试写出求解Laplace方程边值问题

专业知识整理分享

WORD格式可编辑

2u2u2216,0x4,0y3yxu(0,y)y(3y),u(4,y)0,0y3

u(x,0)sinx,u(x,3)0,0x44的五点差分格式，取步长h=1,并写出差分方程的矩阵形式。

9-3 试用古典显格式求解热传导方程定解问题

u2u,0x1,0tTtx2

u(x,0)4x(1x),0x1u(0,t)u(1,t)0,0tT只计算k=1,2两层上的差分解，取网比r=1/6, h=0.2 。

9-4 用古典隐格式求解热传导方程定解问题

u2u,0x1,0t0.3tx2

u(x,0)sinx,0x1u(0,t)u(1,t)0,0t0.3取0.1,h0.2精确解为

u(x,t)e2tsinx。

9-5 导出求解热传导方程的差分格式

(1)1ukukjjk1ukjuj1kkk(u2uuj1jj1)

2h的截断误差，并选取θ使其达到二阶。

专业知识整理分享