第九章 偏微分方程差分方法汇总(完整资料).doc

【最新整理，下载后即可编辑】

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

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

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

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

G y x y x f y

u

x u u ∈=??+??-≡?-),(),,()(2222

（9.1）

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

第二边值条件 ),(y x n

u

β=??Γ

（9.3）

第三边值条件 ),()(

y x ku n

u

γ=+??Γ

（9.4）

这里，n 表示Γ上单位外法向，α(x,y ),β(x,y ),γ(x,y )和k (x,y )

都是已知的函数，k (x,y )≥0。满足方程（9.1）和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。

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

x =ih 1, i =0,1,…,N 1, h 1=a /N 1 y =jh 2, j =0,1,…,N 2, h 2=b /N 2

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

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

h

i i i i i i i i G y x y x f y x y

u

y x x u ∈=??+??-),(),,()],(),([2222

（9.5）

需进一步离散（9.5）中的二阶偏导数。为简化记号，简记节点（x i ，y i ）=(i ,j )，节点函数值u （x i ，y i ）=u (i ,j )。利用一

元函数的Taylor 展开公式，推得二阶偏导数的差商表达式

)(0)]1,(),(2)1,([1),()(0)],1(),(2),1([1

),(2

222

22

212122h j i u j i u j i u h j i y u h j i u j i u j i u h j i x u +-+-++=??+-+-++=??

代入（9.5）式中，得到方程（9.1）在节点(i ,j )处的离散形式

h

j i G j i h h f j i u j i u j i u h j i u j i u j i u h ∈++=-+-+--+-+-

),(),(0)]1,(),(2)1,([1)],1(),(2),1([12221,2221

其中),(,i i j i y x f f =。舍去高阶小项)(02221h h +，就导出了u (i ,j )的近似值u i ,j 所满足的差分方程

h j i j i j i j i j i j i j i G j i f u u u h u u u h ∈=+--+--

-+-+),(,]2[1]2[1,1,,1,22

,1,,121

（9.6）

在节点(i ,j )处方程（9.6）逼近偏微分方程（9.1）的误差为

)(2

221h h O +，它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。

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

h j i j i j i j i j i j i G j i f u u u u u h

∈=-+++-

-+-+),(,]4[1

,,1,1,,1,12 差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值u i ,j ,(i ,j )∈G h 外，还包括边界点值。例如，点(1,j )处方程就含有边界点未知量u 0,j 。因此，还要利用给定的边值条件补充上边界点未知量的方程。 对于第一边值条件式（9.2），可直接取u i ,j =α（x i ，y i ）， (i ,j )∈Γh （9.7）

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

以左边界点(1,j )为例，见图9-2, 利用一阶差商公式

)(),1(),0(),0(11

h O h j u j u j n u +-=??

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

j

j j j

j r u k h u u ,0,0,01

,1,0=+-

（9.8）

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

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

G y x y x f Eu y

u D x u C y u B y x u A x ∈=+??+??+????+????-),(),,(])()([

（9.9）

其中A (x ,y )≥A m in >0, B (x ,y ) ≥B m in >0, E(x ,y ) ≥0。引进半节点,12

1

2

1

h x x i i ±=±

,22

12

1

h y y

i i ±=±利用一阶中心差商公式，在节点（i ,j ）处可有

)(2),1(),1(),()

(])

,1(),(),(),1([1)()],2

1

)((),21)([(1),)((211211,2

11,211211h O h j i u j i u j i x u h O h j i u j i u A h j i u j i u A h h O j i x u A j i x u A h j i x u A x j i j i +--+=??+----+=

+-??-+??=????-+对y

u

y u B y ??????),(类似处理，就可推得求解方程（9.9）的差分

方程

h

j i j i j i j i j i j i j i j i j i j i G j i j i f u a u a u a u a u a ∈=-+++---++---+),(),,(][,,1,1,1,1,,1,1,1,1

（9.10） 其中

???

????

?????

???++++=-=+=-=+=-+--+----+-+---+-+j

i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i E B B h A A h a D h B h a D h B h a C h A h a C h A h a ,21,21,2

2,21,2121,,221

,221,,22

1,2

2

1,,1,21

21,1,1,212

1,1)()()2()2()2()2(

（9.11）

显然，当系数函数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）。容易看出，差分方程（9.10）的截断误差为)(2221h h O +阶。

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

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

考虑Poisson 方程第一边值问题

??

?Γ

∈=∈=?-),(),,(),(),,(y x y x u G y x y x f u α

（9.12）

其中G 可为平面上一般区域，例如为曲边区域。仍然用两组平行直线：x =x 0+ih 1,y =y 0+jh 2,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 '，非正则内点集合为h Γ'。显然，当G 为矩形区域时，h h h h

G G Γ=Γ'=',成立。

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

h j i j i j i j i j i j i j i G j i f u u u h u u u h '∈=+--+--

-+-+),(,]2[1

]2[1,1,,1,2

2

,1,,121

（9.13）

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

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

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

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

a.直接转移法.取与点B 距离最近的边界点（如图9-4中E 点）上的u 的值作为u (B )的近似值u B ，即u B =u (E)=α(E) 直接转移法的优点是简单易行，但精度较低，只为一阶近似。

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

)()()()(21h O C u x x x x A u x x x x B u A

C A

B A

C B C +--+--=

则得到点B 处的方程

A B C B x x u h A h h u -=+++=

δδ

δ

αδ,)(111 线性插值法精度较高，为二阶近似。

对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般区域上边值问题（9.12）的差分近似解。

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

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

9.2 抛物型方程的差分方法

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

9.2.1 一维问题

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

T t l x t x f x

u

a t u ≤<<<+??=??0,0),,(22

（9.14）

l

x x x u <<=0),()0,(?

（9.15）

T t t g t l u t g t u ≤≤==0),(),(),(),0(21 （9.16）

其中a 是正常数，)()(),(),,(21t g t g x t x f 和?都是已知的连续的函数。

现在讨论求解问题（9.14）-(9.18)的差分方法。首先对求解区域G ={0≤x ≤l, 0≤t ≤T }进行网格剖分。取空间步长h =l/N ,时间步长τ=T /M,其中N ,M 是正整数，作两族平行直线

M

k k t t N

j jh x x k j ,1,0,,,1,0,

======τ

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

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

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

O k j u k j u k j t u +-+=??

（9.17）

)()1,(),(),(ττ

O k j u k j u k j t u +--=?? （9.18） )(2)1,()1,(),(2ττ

O k j u k j u k j t u +--+=??

（9.19）

)(),1(),(2),1(),(2

222h O h

k j u k j u k j u k j x u +-+-+=??

（9.20）

1.古典显格式

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

)()

,1(),(2),1()

,()1,(22

h O f h

k j u k j u k j u a

k j u k j u k j +++-+-+=-+ττ 其中

),(k i k j t x f f =。舍去高阶小项)(2h O +τ，就得到节点近似

值（差分解）k j u 所满足的差分方程

k j k

j k j k j k j

k j f h

u u u a

u u ++-=--++2

1

112τ

（9.21）

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

?????

??===-==-=-=++-+=-++M

k t g u t g u N j x u M k N j f ru u r ru u k k N k k j j k j k j k j k j k j ,1,0),(),(1

,,2,1),(1,,1,0,1,,2,1)21(210

111?τ

（9.22） 其中2

h a

r τ

=称为网比。

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

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

?

??=-=+=+φ0

11

,,1,0,u M k F Au u k k k

（9.23） 其中

??

?

??

??

?

????????---=r r r r r r r r A 212121

T

N T k k N k N k k k k T

k N k k x x t rg f f f t rg f F u u u ))(,),(())

(,,,),(()

,,(1121221111----=++==???ττττ

2. 古典隐格式

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

)(),1(),(2),1()

1,(),(2

2

h O f h

k j u k j u k j u a

k j u k j u k j +++-+-+=--ττ

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

k j k

j k j k j k j

k j f h

u u u a

u u ++-=--+-2

1

11

2τ

（9.24）

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

?????

??===-===-=+=-++---+M

k t g u t g u N j x u M k N j f u ru u r ru k k N k k j j k j k j

k j k j k j ,1,0),(),(1

,,2,1),(,,1,0,1,,2,1)21(210

111?τ

（9.25）

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

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

???==+=-φ

1,2,1,u M k F u Bu k k k

（9.26） 其中

??

?

??

??

?????????+---+--+=r r r r r r r r B 212121

向量?,,k k F u 同（9.23）式中定义。从（9.26）式看到，古典隐格式在每一层计算时，都需求解一个三对角形线性方程

组，可采用追赶法求解。

3.Crank-Nicolson 格式（六点对称格式）

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

)()]1,(),([21)21,()

(),()1,()21,(2

2

2

22222τττ

O k j x

u k j x u k j x u O k j u k j u k j t u ++??+??=+??+-+=+?? 使用这两个公式，在)2

1,(+k j 点离散偏微分方程（9.14），然

后利用（9.20）式进一步离散二阶偏导数，则可导出差分方程

2

1

2

1

12111111]22[21+-++-++++++-++-=-k j k j k j k j k j k j k j k j

k j f h

u u u h u u u a u u τ

（9.27）

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

?????

????===-==-=-=++-+=-++-+-++-+++M

k t g u t g u N j x u M k N j f ru u r ru ru u r ru k k N k k j j k j k j k j k j k j k j k j ,1,0),(),(1

,,2,1),(1,,1,0,1,,2,12)1(2)1(2210

21

1111

111?τ

（9.28）

它的矩阵形式为

?????-===++=+++1

,,1,0,,2,1,)()(0211M k u M

k F u A I u B I k k k φ

（9.29）

在每层计算时，仍需求解一个三对角形方程组。

4. Richardson 格式

利用公式（9.19）和（9.20），可导出另一个截断误差为)(22h O +τ阶的差分方程

k j k

j k j k j k j

k j f h u u u a

u u ++-=--+-+2

1

11122τ

称之为Richardson 格式。可改写为

k j k

j k j k j k j k j f u u u r u u τ2

)2(21111++-+=-+-+

(9.32)

这是一个三层显式差分格式。在逐层计算时，需用到1-k j u 和

k j u 两层值才能得到k +1层值1

+k j u 。这样，从第0层已知值

)(0j j x u ?=开始，

还须补充上第一层值1j u ，才能逐层计算下去。可采用前述的双层格式计算1j u 。

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

（1）收敛性。对任意固定的节点（x j ，t k ）,当剖分步长

0,→h τ时，差分解k j u 应收敛到精确解u （x j ，t k ）

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

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

作为例子，先考查Richardson 格式的稳定性。设k j u 是当计算过程中带有误差时，按Richardson 格式（9.30）得到的实际计算值。k j u 是理论值，误差k j k j k j u u e -=。假定右端项k j f 的计算是精确的，网比2

1=r ，则k j e 满足 )2(1111k j k j k j k j k j e e e e e -+-++-+=

（9.31）

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

01,0,,00

j j e e e k j k j k j ≠===-ε。

则利用（9.31）式可得到误差ε的传播情况，见表9-1。 表9-1 r=1/2时Richardson 格式的误差传播

然是就网比2

1=r 进行的，实际上对任何r>0都会产生类似

现象，所以Richardson 格式是不稳定的。

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

9.2.2 差分格式的稳定性

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

???=+=-φ

1u F Hu u k k k

（9.32）

其中H 称为传播矩阵。对于显格式H =A , 隐格式H =B -1,六

点对称格式H =(I +B ) -1 (I+A )。一般的三层格式也可以转化为双层格式。 为了讨论方便，设在初始层产生误差0ε，且假定右端项k F 的计算是精确的。用k u 表示当初始层存在误差0ε时，由差分格式（9.32）得到的计算解，则k u 满足方程

???+=+=-001ε

φu F u H u k k k

（9.33） 记误差向量k k k

u u -=ε，则k ε满足方程

???==-为初始误差

1,2,1,εεε k H k k

（9.34）

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

000,0,ττεε<<≥≤k C k （9.35）

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

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

从（9.34）式递推得到

τ

εεT

k H k k ≤

≤=0,0

因此，差分格式稳定的充分必要条件是

τ

T

k C H k ≤

≤≤0,

（9.36）

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

τ

ρM H +≤1)(

（9.37）

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

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

??

?

????

?????????=0110110110 S 这个特殊矩阵的特征值为 1,,2,1,cos

2-==N j N

j S j

π

λ

（9.38）

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

sin 412-=-=N j N

j r j

πλ

为使稳定性条件式（9.39）成立，必须且只须2

1≤r 。由于

H =A 为实对称矩阵，所以古典显格式稳定的充分必要条件

是网比2

1≤r

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

2-=+=-N j N

j r j

πλ

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

例9-3 六点对称格式 此时H=(I+B)-1(I+A)，利用矩

阵A 和B 的特征值可得到矩阵H 的特征值为

1,,2,1,)

2(sin 42)2(

sin 4222-=+-=

N j N

j r N j r j π

πλ

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

习题九

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

?????Γ∈=∈=??+??Γ),(,

0|),(,162222y x u G

y x y

u

x u 其中G ={-1

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

???

????≤≤==≤≤=-=<<<<=??+??4

0,0)3,(,4sin )0,(30,0),4(),3(),0(3

0,40,162222x x u x x u y y u y y y u y x y u

x

u π

的五点差分格式，取步长h =1,并写出差分方程的矩阵形式。

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

???????≤≤==≤≤-=≤<<

t t u t u x x x x u T t x x u

t u 0,0),1(),0(10),1(4)0,(0,10,2

2

只计算k =1,2两层上的差分解，取网比r=1/6, h =0.2 。

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

3

.00,0),1(),0(10,sin )0,(3.00,10,2

2≤≤==≤≤=≤<<

-=。

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

)2(1)

1(1121

1k

j k j k j k j

k j k j

k j u u u h

u u u u -+-++-=

---+τ

θ

τ

θ

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