偏微分方程数值解法

偏微分方程数值解法 Prepared on 22 November 2020
第十章偏微分方程数值解法
偏微分方程问题，其求解十分困难。

除少数特殊情况外，绝
大多数情况均难以求出精确解。

因此，近似解法就显得更为重要。

本章仅介绍求解各类典型偏微分方程定解问题的差分方法。

§1差分方法的基本概念
几类偏微分方程的定解问题
椭圆型方程：其最典型、最简单的形式是泊松（Poisson）方程
特别地，当
)
,
(≡
y
x
f
时，即为拉普拉斯（Laplace）方程，又
称
为调和方程
Poisson方程的第一边值问题为
其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，
ΓΩ
称为定解区域，
)
,(y
x
f
，
)
,
(y
x
ϕ
分别为Ω，Γ上的已知
连续函数。

第二类和第三类边界条件可统一表示为
其中n
为边界Γ的外法线方向。

当0
=
α时为第二类边界条件，0
≠
α时为第三类边界条件。

抛物型方程：其最简单的形式为一维热传导方程方程可以有两种不同类型的定解问题：
初值问题
初边值问题
其中
)
(x
ϕ
，
)(
1
t
g
，
)(
2
t
g
为已知函数，且满足连接条件
边界条件
)(
),(
),
(
),0(
2
1
t
g
t l
u
t
g
t
u=
=
称为第一类边界条
件。

第二类和第三类边界条件为
其中
)(
1
≥
t
λ
，
)(
2
≥
t
λ。

当
)(
)(
2
1
≡
=t
tλ
λ
时，为第二
类边界条件，
否则称为第三类边界条件。

双曲型方程:
最简单形式为一阶双曲型方程
物理中常见的一维振动与波动问题可用二阶波动方程
描述，它是双曲型方程的典型形式。

方程的初值问题为
边界条件一般也有三类，最简单的初边值问题为
差分方法的基本概念
差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。

它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。

如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解
问题的解，则差分格式的解就作为原问题的近似解（数值解）。

因此，用差分方法求偏微分方程定解问题一般需要解决以下问题： （1）选取网格；
（2）对微分方程及定解条件选择差分近似，列出差分格式； （3）求解差分格式；
（4）讨论差分格式解对于微分方程解的收敛性及误差估计。

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

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

（1） 选取网格： －2h-h0h2h3h
首先对定解区域}0,),{(≥+∞<<∞-=t x t x D 作网格剖
分，最简单
常用一种网格是用两族分别平行于
x
轴与
t
轴的等距直线
kh x x k ==，
(0,1,2
,0,1,2,)j t t j k j τ===±±=将D 分成许
多小矩形
区域。

这些直线称为网格线，其交点称为网格点，也称为节点，
h 和τ
分别称作
x 方向和t
方向的步长。

这种网格称为矩形网格。

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

方程0
u u
a t x
∂∂+=∂∂ 在节点
),(j k t x 处可表示为
其中
(,0)()(0,1,2,)k k u x x k ϕ==±±。

由于
当
τ
,h 足够小时，在式
中略去
),(j k t x R ，就得到一个与方程相近似的差分方程
此处，j k u ,可看作是问题的解在节点),(j k t x 处的近似值。

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

式
称为差分方程的截断误差。

如果一个差分方程的截断误差为)(p
q h O R +=τ，则称差分
方 程对
t
是
q 阶精度，对x 是p 阶精度的。

显然，截断误差的阶数
越大，差分方程对微分方程的逼近越好。

若网格步长趋于0时，差分方程的截断误差也趋于0，则称 差分方程与相应的微分方程是相容的。

这是用差分方法求解偏微 分方程问题的必要条件。

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

§2椭圆型方程第一边值问题的差分解法
本节以Poisson 方程为基本模型讨论第一边值问题的差分方法。

差分格式的建立
考虑Poisson 方程的第一边值问题 取τ
,h 分别为
x 方向和y 方向的步长，如图所示，以两族平行
线
kh x x k ==，(,0,1,2,)j y y j k j τ===±±将定
解区域剖分成矩形网 格。

节点的全体记为
{(,),,,}k j k j R x y x kh y j k j τ===为整数。

定解区
域内部的节点称为内点，记内点集Ω R 为τ
h Ω。

边界
Γ与网
格
线的交点称为边界点，边界点全体记为
τ
h Γ。

与节点
),(j k y x 沿x 方
向或
y
方向只差一个步长的点
),(1j k y x ±和),(1±j k y x 称为节点
),(j k y x 的
相邻节点。

如果一个内点的四个相邻节点均属于ΓΩ ，称为正
则内点，正内点的全体记为
)
1(Ω
，至少有一个相邻节点不属于
ΓΩ
的内点称为非正则内点，非正则内点的全体记为)
2(Ω。

问题是要
求出第一边值问题在全体内点上的数值解。

为简便，记
),(),(j k y x j k =，(,)(,)k j u k j u x y =，
),(,j k j k y x f f =。

对正则
内点)
1(),(Ω∈j k ，由二阶中心差商公式
Poisson 方程22
22(,)u u
f x y x y
∂∂+=∂∂在点
)
,(j k 处可表示为 其中
为其截断误差表示式，略去
),(j k R ，即得与方程相近似的差分方程
式中方程的个数等于正则内点的个数，而未知数j k u ,则除了包含
正则内点处解
u 的近似值外，还包含一些非正则内点处u 的近
似值，因而方程个数少于未知数个数。

在非正则内点处Poisson 方程的差分近似不能按上式给出，需要利用边界条件得到。

边界条件的处理可以有各种方案，下面介绍较简单的两种。

（1）直接转移
用最接近非正则内点的边界点上的u 值作为该点上u
值的
近似，这就是边界条件的直接转移。

例如，点),
(j k P 为非正则内点，其最接近的边界点为Q
点，则有
上式可以看作是用零次插值得到非正则内点处u 的近似值，容易
求出，其截断误差为)(τ+h O 。

将上式代入，方程个数即与未知数
个数相等。

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

由点
R 与
T
的线性插值确定
)(P u 的近似值j k u ,，得
其中RP
d =，其截断误差为
)(2
h O 。

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

上面所给出的差分格式称为五点菱形格式， 实际计算时经常取τ=h ，此时五点菱形格式可化为
简记为
2
1
h ◇
j k j k f u ,,=
其中◇
j k j k j k j k j k j k u u u u u u ,1,1,,1,1,4-+++=-+-+。

例1用五点菱形格式求解拉普拉斯（Laplace ）方程第一边值问题
其中}1,0),{(≤≤=Ωy x y x 。

取3
1
==τh 。

（0,0）(1,0)(2,0)(3,0)
[解]网格中有四个内点，均为正则内点。

由五点菱形格式，得方 程组
⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-+++=-+++=-+++=-+++0
)4(10)4(10)4(10)4(12,21,23,22,12,322,11,13,12,02,221,20,22,21,11,321,10,12,11,01,22u u u u u h
u u u u u h u u u u u h u u u u u h 代入边界条件
1,02,00,10,21,32,33,13,21625
lg ,lg
991013
lg ,lg
992534
lg ,lg
993740
lg ,lg
99
u u u u u u u u ⎧
==⎪⎪
⎪==⎪⎨
⎪==⎪⎪⎪==⎪⎩
其解为
2756919
.01,1=u ，
4603488.01,2=u ，
3467842
.02,1=u ，5080467
.02
,2=u
当
τ=h 时，对
1,1,,1,1,,21(4)k j k j k j k j k j k j u u u u u f h
+-+-+++-= 利用点),(j k ，)1,1(-±j k ，)1,1(+±j k 构造的差分
格式，称为五点
矩形格式，简记为
2
21h □
j k j k f u ,,=
其中□
j k j k j k j k j k j k u u u u u u ,1,11,11,11,1,4-+++=--+--+++，其
截断误差为
五点菱形格式与矩形格式的截断误差均为
)(2h O ，称它们具有 二阶精度。

如果用更多的点构造差分格式，其截断误差的阶数可 以提高，如利用菱形格式及矩形格式所涉及的所有节点构造出的 九点格式就是具有四阶精度的差分格式。

§3抛物型方程的差分解法
以一维热传导方程
为基本模型讨论适用于抛物型方程定解问题的几种差分格式。

差分格式的建立
首先对
xt 平面进行网格剖分。

分别取τ,h 为x
方向与
t
方向
的步长，用两族平行直线
(0,1,2,)k x x kh k ===±±，
)2,1,0( ===j j t t j τ，将xt
平面剖分成矩形网格，节
点为(,)
(0,1,2,
,0,1,2)k j x t k j =±±=。

为
简便，记),(),(j k t x j k =，),(),(j k t x u j k u =，
)(k k x ϕϕ=，)(11j j t g g =，
221122(),(),()j j j j j j g g t t t λλλλ===。

（一）微分方程的差分近似
在网格内点
),(j k 处，对
t
u ∂∂分别采用向前、向后及中心差商公式
一维热传导方程 可分别表示为
由此得到一维热传导方程的不同差分近似
上述差分方程所用到的节点各不相同。

其截断误差分别为
)(2h O +τ，)(2h O +τ和)(22h O +τ。

因此，它们都与一
维热传导方程相容。

如果将式
中的j k u ,用
)(2
1
1,1,-++j k j k u u 代替，则可得到又一种差分近似
差分方程用到四个节点。

由Taylor 公式容易得出
故其的截断误差为
⎪⎪⎭
⎫ ⎝
⎛++2
222
)(h O h O τ
τ。

因而不是对任意的0,→τh ，
此差分方程都能逼近热传导方程
)0(022>=∂∂-∂∂a x
u
a t u ，
仅当
()o h τ=时，才成立。

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

构造差分格式的关键在于使其具有相容性、收敛性和稳 定性。

前面三个方程都具有相容性，而此方程则要在一定条件下 才具有相容性。

（二）初、边值条件的处理
为用差分方法求解定解问题初值问题 初边值问题
还需对定解条件进行离散化。

对初始条件及第一类边界条件，可直接得到
其中
τ
T m h l n ==,
对第二、三类边界条件
需用差分近似。

下面介绍两种较简单的处理方法。

（1）
在左边界)0(=x 处用向前差商近似偏导数
x
u
∂∂，
在右边界
)(l x =处用向后差商近似
x
u ∂∂，即
则得边界条件的差分近似为
其截断误差为
)(h O 。

（2）用中心差商近似
x
u ∂∂,即
则得边界条件的差分近似为 其截断误差为)(2
h
O 。

误差的阶数提高了，但出现定解区域外的
节点),1(j -和),1(j n +，这就需要将解拓展到定解区域外。

可以通
过用内节点上的
u
值插值求出
j u ,1-和j n u ,1+，也可以假定热传导方
程在边界上也成立，将差分方程扩展到边界节点上，由此消去j u ,1- 和
j
n u ,1+。

（三）几种常用的差分格式
以热传导方程的初边值问题 为例给出几种常用的差分格式。

（1）古典显式格式
令
2h
a r τ=，则
022
,1,,1,1,=+----++h
u u u a
u u j
k j k j k j
k j k τ
可改写成
j k j k j k j k ru u r ru u ,1,,11,)21(-+++-+=
将其与初始条件及第一类边界条件
结合，我们得到求解此问题的一种差分格式
由于第0层
)0(=j 上节点处的u 值已知)(0,k k u ϕ=，由此即可算出
u 在
第一层
)1(=j 上节点处的近似值1,k u 。

重复使用此式，可以逐层计
算出所有的
j k u ,，因此此差分格式称为古典显式格式。

又因式中
只出现相邻两个时间层的节点，故此式是二层显式格式。

（2）古典隐式格式
将式
22
,1,,11
,,=+----+-h
u u u a
u u j
k j k j k j k j k τ
整理并与初始条件及第一类边界条件式联立，得差分格式如下
其中
2h
a r τ=。

虽然第0层上的
u
值仍为已知，但不能由上式直
接计算以上各层节点上的值
j k u ,，必须通过解下列线性方程组
才能由
j k u ,计算1,+j k u ，故此差分格式称为古典隐式格式。

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

（3）Richardson 格式
Richardson 格式是将式
整理后与初始条件及第一类边界条件式联立。

其计算公式为 这种差分格式中所涉及的节点出现在
1,,1+-j j j 三层上，故为三层
显式格式。

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

（4）杜福特-弗兰克尔（DoFort-Frankel ）格式
DoFort-Frankel 格式也是三层显式格式，它是由式 与初始条件及第一类边界条件式结合得到的。

具体形式如下：
用这种格式求解时，除了第0层上的值
0,k u 由初值条件得到，必
须先用二层格式求出第1层上的值
1,k u ，然后再按上式逐层计
算
),,3,2(,m j u j k =。

（5）六点隐式格式
对二阶中心差商公式
如果用
22x
u
∂∂在点
)1,(+j k 与点),(j k 处的二阶中心差商的平
均值
近似
2
2x u ∂∂在
)2
1
,(+j k 处的值，即 同时
t
u ∂∂在点
)2
1
,(+j k 处的值也用中心差商近似，即 这样又得到热传导方程的一种差分近似 其截断误差为)(22
h O +τ
，将上式与初始条件及第一类边界条件式
联立并整理，得差分格式
此格式涉及到六个节点，它又是隐式格式，故称为六点隐式格式。

与古典隐式
格式类似，用六点格式由第
j 层的值j k u ,计算第1+j
层的值
1,+j k u 时，需求解三对角方程组
此方程组的系数矩阵严格对角占优，故仍可用追赶法求解。

例2用古典显式格式求初边值问题 的数值解，取
5.0,1==τh 。

[解]这里21,0.5a a r h
τ===，2)(x x =ϕ，9)(,0)(21==t g t g 。

由格式
可得到
将初值
0,k u 代入上式，即可算出
将边界条件
01,0=u ，91,3=u 及上述结果代入又可求得
如此逐层计算，得全部节点上的数值解为
§4双曲型方程的差分解法
对二阶波动方程
如果令1u v t
∂=∂，x u
v ∂∂=2，则方程可化成一阶线性双曲型方程
组
记
T
v )
,(21v v =，则方程组可表成矩阵形式
矩阵
A 有两个不同的特征值a ±=λ，故存在非奇异矩阵P ，使得
作变换
T
w w P )
,(21==v w ，方程组可化为
x
t ∂∂Λ=∂∂w
w
方程组由二个独立的一阶双曲型方程联立而成。

因此本节主要讨 论一阶双曲型方程的差分解法。

几种简单的差分格式 考虑一阶双曲型方程的初值问题 将
xt 平面剖分成矩形网格，取x 方向步长为t h ,方向步长为τ，
网格线为
(0,1,2,)k x x kh k ===±±，
)2,1,0( ===j j t t j τ。

为简便，记
),(),(j k t x j k =，
)(),,(),(k k j k x t x u j k u ϕϕ==。

以不同的差商近似偏导数，可以得到方程的不同的差分近似
截断误差分别为
)(h O +τ，)(h O +τ与)(2
h O +τ。

结合离
散化的初始条件，可以得到几种简单的差分格式
其中
r h
τ
=。

如果已知第
j
层节点上的值
j k u ,，按上面三种格式就
可求出第
1+j 层上的值1,+j k u 。

因此，这三种格式都是显式格式。

如果对
t
u ∂∂采用向后差商，x
u
∂∂采用向前差商，则方程可化成
相应的差分格式为
此差分格式是一种隐式格式，必须通过解方程组才能由第
j 层节
点上的值j k u ,求出第1+j 层节点上的值1,+j k u 。

例对初值问题
⎪⎩
⎪⎨⎧+∞
<<∞-=+∞<<∞->=∂∂+
∂∂x x x u x t x u
t u )()0,(,00ϕ其中
⎪⎪⎩⎪
⎪⎨⎧>=<=0
00
2
101)(x x x x ϕ 用差分格式求其数值解
)4,3,2,1(,=j u j k ，取
2
1
=
=h r τ。

[解]记
),2,1,0( ±±==k kh x k ，由初始条件
按差分格式
计算公式为
j
k j k j k u u u ,1,1
,2
1
23++-=
计算结果略。

如果用差分格式
求解，计算公式为计算结果见表略。

与准确解
1
1
(,)
2
00
x t u x t x t
x
<
⎧
⎪⎪
==
⎨
⎪
>
⎪⎩
比较知，按前一个差分格式所求得的数值解不收敛到初值问题的解，而后一个差分格式的解收敛到原问题的解。