浙江大学 实用数值计算方法 1
第七章 偏微分方程
7.1 一般介绍
7.2 一阶双曲型方程的差分求解法
7.3 一阶双曲型方程的特征线求解法
7.4 一阶双曲型方程的线上求解法
7.5 二阶椭圆型方程的差分求解法
7.6 二阶椭圆型方程的有限元求解法
7.7 二阶椭圆型方程的加权残差求解法
7.8 二阶抛物型方程的差分求解法
7.9 二阶抛物型方程的线上求解法
7.10 二阶双曲型方程的特征线求解法浙江大学 实用数值计算方法 2
7.1 偏微分方程的一般介绍
Partial Differential Equations(PDEs)
自变量数 至少 2个
,,,,
0,,,,,,,,
,,
2
2
2
yx
u
u
x
u
u
y
u
u
x
u
u
uuuuuuyxF
yxuu
xyxxyx
yyxyxxyx
阶数 方程中导数的最高阶数
0
0
0
3
3
yyyx
yxx
yx
uu
uu
ubu
三阶二阶一阶
性态以一阶方程为例
0 cubua yx
浙江大学 实用数值计算方法 3
7.1
yyxyxxyx
yx
yxyyxyxx
yx
yx
yx
yx
uuuuuuyx
uuuyx
yx
CBA
FuEuDuCuBuA
uu
xuuu
buu
uuuyx
uyx
yx
yxccyxbbyxaa
cba
,,,,,,,:
,,,,:
,:
:,:,:
0:::
0
0
0
,,,,
,,
,
,,,,,
,,
2
2
：
对二阶方程例为常数非线性拟线性线性
N onlinear
rQ uasilinea
L inear
非线性拟线性线性非线性拟线性线性浙江大学 实用数值计算方法 4
7.1
类型一阶栓区型方程
xt
yx
uvu
cubua
0流动方程
Advection Equation(AE)
二阶线性方程
0 GCuBuAu yyxyxx
e llip t icACB 椭圆型042
方程，传热方程方程
Lap la c euu
P o is s o nyxfuu
yyxx
yyxx
0
,
p a r a b o licACB 抛物型042
方程扩散方程
B u r g e ruKuuu
uu
yyyx
yyx
h y p e r b o licACB 双曲型042
波动方程yyxx uu?
浙江大学 实用数值计算方法 5
7.1
求解方法有限差分法 Method of Finite Differences (MFD)
特征线法 Method of Characteristics (MOC)
线上求解法 Method of Lines (MOL)
有限元素法 Method of Finite Elements (MFE)
加权残差法 Method of Weighled Residuals (MWR)
问题收敛性 Convergence
当采取的步骤趋于无限时，数值结果是否趋于理论值？
稳定性 Stability
在某一步引入的误差，经多步数值计算后，会扩大或抑制？
dx
dy
x
y?
k
k
ey
y
ey
y?
1
1
浙江大学 实用数值计算方法 6
7.2 一阶双曲型方程的差分求解法
0 xuvtu
或称流动方程 Advective
Advection Equation (AE)
v为流速因子该方程的介折解
vtxfu
dw
du
x
w
dw
du
x
u
dw
du
v
t
w
dw
du
t
u
vtxwwfu
,因为求具体解时需要提供 2个辅助条件
vtxftxu
vtxftxu
00
00
,
,
浙江大学 实用数值计算方法 7
7.2
assuming the forcing function is a Rump
The solution of is shown below.
sW
Ws
s
W
W
Wf
,0.0
0,1
0,0.1
0 xt vuu
wf
W0s?
0.1
0.0
u
t
x
0tt?
0xx?
jx
nt
ntt?
jxx? tvx
txu j,
txu,0
0,txu
ntxu,
图 7.1 Propagation of the Wave Front
浙江大学 实用数值计算方法 8
7.2.1 最简单的差分化格式构想
t
x0x 1?jx jx 1?jxx?
1?nt
nt
0t
t?
tbxtu
xaxtu
x
uv
t
u
0
0
,
,xjxx j 0
tnttn 0
211
,
1
,
2
xO
x
uu
x
u
tO
t
uu
t
u
n
j
n
j
nj
n
j
n
j
nj
x
uuv
t
uu njnjnjnj
2
111
njnjnjnj uuxtvuu 111 2 1nn
图 7.2
浙江大学 实用数值计算方法 9
7.2.1
以上方法称为 时间镶嵌空间中心 的差分表达
Forward Time Centered Space
FTCS represetation
实际上这个方法不能用,不稳定的方法
Unstable Method
考虑数据误差 r
由于原方程为线性，故误差的传播关系是与原方程完全相同的差分方程差分方程独立解的一般形式
Independent Solutions of Difference Equations
njnjnjnjnjnjnjnj ruruxtvruru 111111 2
njnjnjnj rrxtvrr 111 2
xjkir
xjkir
nn
j
nn
j
e x p
e x p
11?
F a c t o rio nA m p lif ic a t称为放大因子?
是否上升序列决定,,,,,,1110 njnjnjjj rrrrr
浙江大学 实用数值计算方法 10
7.2.1
线性微分方程的解
xFQydxdyPdx yd22
xFQyDyPyD2
应为补充解和特殊解之和 补充解系由下式求出
0
0
0
21
2
2
yPDPD
yQDPD
QyDyPyD
e q u a t io na u x ilia r y
QPxxPP 的两个根为方程式 0,221
00 21 yPDyPD
补充解系由两个 独立解组成
xPAyxPAy 222111 e x pe x p
为复数时当 iPQP
yPxPPADy
4
e x p
2
111111
pyyyy 21
决定。由解，特殊解解组成的补充为微分方程的两个独立
xFy
yy
p
21?
浙江大学 实用数值计算方法 11
7.2.1
差分方程的解可用算符运算方法 Operator Calculus 导出差分算符 Difference Operator
jjjjj jjj xfxfyyy
xxx
11
1
jjjjj Eyyyyy 11
,1E
2
2
21
1
jj
jj
jj
yyE
yEy
yEy
它和微分算符一样，是一种线性算符用于线性二阶差分方程
xFQyPyy nnn 12
和微分方程类似，它的 补充解 可由下式得到
0
0
0
21
2
12
n
n
nnn
yPEPE
yQPEE
QyPyy
浙江大学 实用数值计算方法 12
7.2.1
故补充系由两个独立解组成
（ Independent Solutions）
0
0
2
1
n
n
yPE
yPE
两个独立解为
nnnn PAyPAy 2211 21
nnnnn yPEyyPPAy,11因为：
nnn PAPAy 2211
ixiPQP e x p42时，当
xniAy nn e x p? 差分方程的一个独立解 （ Eigenmode）
nn
n
n
n
nn
n
PyEy
in xAixpy
xniAyEy
aa
ai
aa
aiaa
iaia
n
aa
aa
e x pe x p
1e x p
!5!3!4!2
1
!4!3!2
1e x p
!!2
1e x p
1
1
5342
432
2
浙江大学 实用数值计算方法 13
7.2.1
差分方程独立解的一般形式
in xAy nn e x p
用于本题的情况
xik jr
xik jr
nn
j
nn
j
e x p
e x p
11?
乘以随时间坐标的推进不断njr
F a c t o rio nA m p lif ic a t称为“放大因子”?
时，为不稳定。当 1
将独立解代入差分表达式
njnjnjnj rrx tvrr 111 2
得到
xki
x
tv
xikxik
x
tv
xik j
xjikxjik
x
tv
xjikxjik
x
tv
xik j
nn
n
s i n2
2
e x pe x p
2
e x p
1e x p1e x p
2
1
1e x p1e x p
2
e x p
1
浙江大学 实用数值计算方法 14
7.2.2 差分格式的改进
的条件是使相应的放大因子为：
格式为：
得到的差分而代之以不用
1
s inc os
22
1
2
1
2
1111
1
11
i
Xk
x
tv
iXk
uu
x
tv
uuu
uu
u
M e t hodL ax
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
1xtv
Courant Condition
t
x0x 1?jx 1?jxjx
nt
0t
1?nt
I
R?
tk
1xtv?
图 7.3
图 7.4
浙江大学 实用数值计算方法 15
7.2.2
Courant 条件的物理意义波形传递系沿 x=vt线
t节点的选取
当节点取在线上：
当节点取在线外：
当节点取在线内：
Lax差分格式也写成以下形式可以看成为以下偏微分方程的 FTCS差分式
1,2 xtvvxt
不稳定1,3 xtvvxt
稳定1,1 xtvvxt
t
uuu
x
uuv
t
uu njnjnjnjnjnjnj 11111 2
2
1
2
2
22
x
u
t
x
x
uv
t
u
dissipative term 耗散项
Numerical Viscosity 数值黏度
t
x1?jx jx 1?jx
x?
nt
1t?
3t?
2t?
图 7.5
浙江大学 实用数值计算方法 16
7.3 一阶双曲型方程的特征线求解法
Method of Characteristics (MOC)
dy
y
u
dx
x
u
du
yxu
A
y
u
BC
x
u
yxCBA
的全微分线性的函数，可以是
,
,,,
dy
y
udx
A
y
uBC
du
CyuBxuA
0 Ad uC d xyuB d xAd y
这是原方程的转换方程，它们的解相同。
为原方程的特征线方程0 B d xAd y
yxFABdxdy,
在特征线上，满足 的为解。0 A duC dx
浙江大学 实用数值计算方法 17
7.3.1 Method of Characteristics
(MOC)
yxFdxdy,?
yxGdxdU,?
,,,,,,,
,,,,,,,
,,,,,,,
,
,,,,,,,
,
,,,,
10
11110
00100
01000
10
0
njnn
j
j
n
n
yxuyxuyxu
yxuyxuyxu
yxuyxuyxu
yxG
dx
du
yxuyxuyxu
yxF
dx
dy
yyyy
xx
yx
积分，可以得到再对条件轨迹。根据给定的初始线的积分可以得到一族特征方程作为初值，对常微分别由时，分从平面上在
y
x
ny
1y
0y
y?
x?0,0u
yxu,
0x jx1x
图 7.6
浙江大学 实用数值计算方法 18
7.3.1 Method of Characteristics (MOC)
yxFdxdy,?
yxGdxdU,?
y
x
ky0
01y
00y
y?
x?0,0u
yxu,?
0x jx1x
ky1
11y
10y
jky
0jy
1jy
kL
0L
1L
,2,1,0,
,
,,,,001000
kL
yxF
dx
dy
yyyxxyx
k
k
到一族特征线的轨迹积分可以得程作为初值，对常微分方时，分别由平面上，从在
,,,,,,,,,:
,,,,,,,,,:
,,,,,,,,,:
221100
12121110101
02021010000
jkjkkkk
jj
jj
yxyxyxyxL
yxyxyxyxL
yxyxyxyxL
图 7.7
浙江大学 实用数值计算方法 19
7.3.1
故积分上式，可以得到
，满足常微分方程因为在每一条特征线上可以计算出根据给定的初始条件，
yxG
dx
du
yxuyxuyxu
yxu
k
,
,,,,,,,
,,
00010000
0
,,,,,,,
,,,,,,,
,,,,,,,
1100
1111010
0101000
jkjkk
jj
jj
yxuyxuyxu
yxuyxuyxu
yxuyxuyxu
进行插值计算。方向上对函数值则还应进行上的解，方向上也规则的离散点此，若需要知道上。因方向不规则的离散点是定义在方法求得的离散解由此可知，用
uy
y
yxy
yxuM O C
jkj
jkj
,
,
浙江大学 实用数值计算方法 20
7.4 一阶双曲型方程的线上求解法
Method of Lines (MOL)
有限差分法：偏微分方程完全离散成为一组差分方程用线性代数方程组求解线上求解法：偏微分方程部分离散成为一组常微分方程用常微分方程积分方法求解
CxuBtuA
ni
x
u
A
B
A
C
dt
du
nixx
tgtxu
xftxu
x
u
A
B
A
C
t
u
tx
i
i
i
,,1,0
,,1,0,
,
,
,
0
0
得到常微分方程组：
轴离散化，将浙江大学 实用数值计算方法 21
线上求解法 Method of Lines (MOL)
nitutxu i,,1,0,,视为把部分离散化方法
t
x0x ix nxx?
1t
0t
t?
线间距积分步长
。常微分方程组的数值解分方法计算然后用任何一种数值积各种差分，样条，的近似值。
导方法得到：可以用任何一种数值求为已知。或初始条件
x
u
tutxu
i
ii
00
,
nixuABACtu ii,,1,0
7.4
图 7.8
浙江大学 实用数值计算方法 22
7.5 二阶椭圆型方程的差分求解法
02
2
2
2
yuxu
称为稳态热传导方程，通式为
Dirichlet 问题
Neumann 问题算符称为方程称为
L a p la c e
x
L a p la c eu
k k
2
2
2
2 0
02
2
2
2
yuxu
02
2
2
2
yuxu
yfyxu
yfyxu
yfyxu
yfyxu
n
m
41
301
2
10
,
,
,
,
ygxu
ygxu
ygxu
ygxu
n
m
y
y
x
x
4
3
2
1
0
0
浙江大学 实用数值计算方法 23
y
x
ny
0y
jy
0x ix mx
xfyxu n 4,?
xfyxu 30,?
02222 yuxu
yxQ, u(x m
,y)=
f 2(y
)
u(x
0,y
)=
f 1(y
)
xgyu
y
3
0
xgyu
ny
4
xgyu
mx
2
xg
yu x 10
Laplace 方程 的 Dirichlet 边界条件 和 Neumann 边界条件 和 Poisson 方程方程Po isson
yxyuxu,2222
边界条件也需 4个，有 3类给定方法
Dirichlet 边 界 条 件
Neumann 边 界 条 件
混合 边 界 条 件
7.5
图 7.9
浙江大学 实用数值计算方法 24
7.5.1 Laplace算符的差分表达
022222 yuxuu
简化表示为级数导出。精度可由二阶导数的差商表示及
2
4
''
2
4
4
3
'''
2
''
'
4
4
3
'''
2
''
'
12
2
2462
2462
h
f
xf
h
hxfxfhxf
h
f
h
xf
h
xf
hxfxfhxf
h
f
h
xf
h
xf
hxfxfhxf
T ay lor
i
iii
iii
iii
iii
iii
22 11'' 2 hOh ffff iiii
用于 Laplace算符
2
11
2
112
,,2,
,,2,
,
y
yxuyxuyxu
x
yxuyxuyxu
yxu
jijiji
jijiji
ii
并用简化表示当,hyx
jijijijijiij uuuuuhu,1,1,,1,122 41
浙江大学 实用数值计算方法 25
1?ix ix 1?ix
1?jy
jy
1?jy
1,?jiu
jiu,jiu,1?jiu,1?
1,?jiu
0
11
161
11
1
0
1
141
1
1
,
2
2
2
2
2
2
2
2
2
2
ijk
ij
u
h
z
u
y
u
x
u
u
u
h
yxu
L a pl ac e
在三维空间内方程的差分表达
7.5.1
图 7.10
浙江大学 实用数值计算方法 26
例,Laplace 方程的 Dirichlet 边界问题
co0
co0co0co0
co0 co0 co0
co100
y
x
cm10
cm20
5?yh
5?xh
解得或出一个线性方程对每个内部节点均可写
1 0 0
0
0
410
141
014
04001 0 0
0400
04000
3
2
1
32
231
12
T
T
T
TT
TTT
TT
cT
cT
cT
o
o
o
786.26
143.7
786.1
3
2
1
7.5.1
图 7.11
浙江大学 实用数值计算方法 27
为了提高精度需要加密网络
1 8 15 5 12
2 9 16 6 13
3 10 17 7 14
4 11 18 1 8 15
5 12 2 9
6 13 3 10
7 14 4 11
987654321 151413121110
9
8
7
6
5
4
3
2
1
100000141000001
100000141000001
10000004100000
1000001410000
100000141000
10000014100
1000001410
100000141
10000014
5,15121 非零元素数带宽ha
7.5.1
图 7.12
浙江大学 实用数值计算方法 28
Laplace 方程 Dirichlet边界问题的差分求解
消去法
直接迭代 Liebmann 方法
相继松弛 S.O.R,方法
交替方向 A.D.I.方法
04,,1,1,,1,12 jijijijijiij uuuuuu
411,1,,1,1, kjijijijikji uuuuu
1,1,,
1,
11,1,,1,1
1,4
kjikjikji
kji
kjijijiji
kji
uRuu
uuuuuuR
1
,1,1,
,,1,1,
1
,
1
,1,1,
,,1,1
1
,,
2
4
2
4
2
4
2
4
k
jijiji
k
jijiji
k
ji
k
ji
k
jijiji
k
jijiji
k
ji
k
ji
uuu
P
uuu
P
uu
uuu
P
uuu
P
uu
7.5.1
浙江大学 实用数值计算方法 29
7.6 二阶椭圆型方程的有限元素法求
Method of Finite Elements (MFE)
以 Laplace 方程的 Dirichlet 问题为例根据变分原则
Variational
Principles
等价性定理以上方程的解将使以下泛函为最小。
yxyxgyxu
Dyx
y
u
x
u
,,,,
,,02
2
2
2
d x d y
y
u
x
uuJ
D
22
2
1
S
D
me?
D
图 7.13
浙江大学 实用数值计算方法 30
7.6
将 D进行剖分，常用的是三角剖分法对任何一个元素用二原线性函数近似在三个顶点上可得到
SDeD
m
m,
yaaayxW e 321,
kkk
jjj
iii
Wyaxaa
Wyaxaa
Wyaxaa
321
321
321
kk
jj
ii
kk
jj
ii
kkk
jjj
iii
wx
wx
wx
e
a
yw
yw
yw
e
a
yxW
yxW
yxW
e
a
1
1
1
2
1
1
1
1
2
1
2
1
3
21
其中
kk
jj
ii
yx
yx
yx
e
1
1
1
2
浙江大学 实用数值计算方法 31
7.6
X
Y
u
ix
jx
kx iy
ky
jy
yxu,
e
yxW,
i
j
k
e?
Ui=Wi
Uk=Wk
Uj=Wj
图 7.14
浙江大学 实用数值计算方法 32
7.6
所以其中既然顶点坐标均为规定，所以并有
kkkk
jjjj
iiii
e
wydxcb
wydxcb
wydxcb
e
yxW
2
1
,
ijkjikijjik
kijikjkijkj
jkikjijkkji
xxdyycyxyxb
xxdyycyxyxb
xxdyycyxyxb
,,
,,
,,
kjie WWWyxW,,,
kkjjii
e
e
y
kkjjii
e
e
x
wdwdwd
ey
W
yxW
wcwcwc
ex
W
yxW
2
1
,
2
1
,
浙江大学 实用数值计算方法 33
7.6
使泛函最小的问题，即对近似为对求极值，或因此得到：
可解得
e e
e
y
e
x
yx
dx dyww
dx dyuuuJ
22
22
2
1
2
1
dx dy
e
wdwdwd
e
wcwcwc
wJ
kkjjii
e e
kkjjii
2
2
2
22
1
nmwJW
m
,,2,1,0
rWA
nmuw mm,,2,1,
n为内部节点数边界上的 W为给定浙江大学 实用数值计算方法 34
7.6
2
2
121
1
,
,
,,c os,c os,
,,,
,
,
,,,,
yx
yxh
yxuyxh
y
u
yxq
x
u
yxp
yxyxgyxu
Dyx
yxf
yxuyxr
y
u
yxq
yx
u
yxp
x
对于更为一般性的情况需要极小化的泛函将是也可剖分为有限个元素后求解
2
2
12
2
22
2
1
,
,,,
2
1
dsuhuhdx d yuyxf
uyxr
y
u
yxq
x
u
yxpuJ
y
x
D
1?
2?
sx
syS
normal
genttan
1?2
图 7.15
浙江大学 实用数值计算方法 35
7.8 二阶抛物型方程的差分求解法动态扩散方程对于一维空间用差商代替微商，可以有各种选择，例如所以有需要另有更方便的方法
t
C
DC?
12
xgtxC
xgtxC
xgtxC
x
C
D
t
C
n 3
20
10
2
2
,
,
,
,
iix
x
jijii
tx
jjt
jiji
tx
xxh
h
CCC
t
C
tth
ht
CC
t
C
ji
ji
12
,1,1
2
2
1
1,1,
,
2
,
2
jijiji
x
tjiji CCChhDCC,1,,1
21,1,22
时的数值和 ji tt 1?
浙江大学 实用数值计算方法 36
7.8
显式方法得到或者：
则有：
22,1,,1
,
2
2
,1,
,
2
x
x
jijiji
tx
t
t
jiji
tx
hO
h
ccc
x
c
hO
h
cc
t
c
ji
ji
jijiji
x
tjiji ccc
Dh
hcc
,1,,12,1,2
jijijiji crccrc,,1,11,21
2
122
2 rDhhh
Dhr
tx
x
t 时，，当其中
jijiji ccc,1,11,21
0x ix 1?ix1?ix nx x
t
1?jt
jt
0t
图 7.16
浙江大学 实用数值计算方法 37
7.8
示例：
取得到的数值解与以下解析解比较
2.0,20
0,0
04.00,
119.0
2
2
tc
tc
xc
x
c
t
c
s e c2.67
1 1 9.0
42
s e c1 1 9.04
2
2
t
x
h
cmDcmh
1
2
1
2
20
12
s i n1200294.0e x p
12
10
s i n01175.0e x p
20
2
,
n
n
xn
tn
xn
tn
x
txC
饱和蒸汽
C2H5OH
32.0,20 cmmgtc304.00,cmmgxc?
s e c11904.0 2cmD?
cm2020 0 x
30.0,0 cmmgtc
空气
2cmA
图 7.17
浙江大学 实用数值计算方法 38
0
1
2
3
4
5
6
840 24201612
%c
cmx 4?
cmx 12?
Number of time steps
Analytical Solutions
Numerical Solutions
s t e pt im ep e rt
cmxs e t
cmD
x
t
Dr
s e c6.334
119.0
25.0
0.4
s e c119.0,25.0
2
2
2
Analytical versus Numerical Solutions
Diffusion Dynamics
r?0.25
7.8
图 7.18
浙江大学 实用数值计算方法 39
0
1
2
3
4
5
6
420 121086
%c
cmx 4?
cmx 12?
Number of time steps
Analytical Solutions
Numerical Solutions
s t e pt im ep e rt
cmxs e t
cmD
x
t
Dr
s e c2.674
119.0
5.0
0.4
s e c119.0,5.0
2
2
2
Analytical versus Numerical Solutions
Diffusion Dynamics
r?0.5
7.8
图 7.19
浙江大学 实用数值计算方法 40
7.8
显式法的稳定性分析所以
jijiji
ji
Wce
tx
,,,
,
时存在误差若在
时，有在令所以代入前式，并根据级数展开
2
1
0,
,
,
,
2
,
2
21
21
2
,1,
2
2
22
,
,,1
2
1
22
,
,,1
1,,,,1,1
,,1,11,
r
MDEe
D
t
x
r
t
xc
tww
x
tcx
x
c
xww
x
tcx
x
c
xww
T ayl or
wwrwwr
ereere
jj
i
jiji
ji
jiji
ji
jiji
jjjijiji
jijijiji
2
2
,,1,11,
,,
21
x
tcD
t
xct
ereere
ji
jijijiji
浙江大学 实用数值计算方法 41
7.8
因此有
tME
tMErrEE
j
jjj
2121
0100
11
1
2
ttMEtMjE
tMEtMEE
j
jjj
则若当所以时，在
0,0,
2
1
,00
11
00
txr
MtE
Et
jj
0
,,
,
2
2
,
2
2
jiji
ji
x
c
D
t
c
x
tc
D
t
xc
M
稳定条件为故算法为稳定即 01jE
2
1
2
x
tDr
第七章 偏微分方程
7.1 一般介绍
7.2 一阶双曲型方程的差分求解法
7.3 一阶双曲型方程的特征线求解法
7.4 一阶双曲型方程的线上求解法
7.5 二阶椭圆型方程的差分求解法
7.6 二阶椭圆型方程的有限元求解法
7.7 二阶椭圆型方程的加权残差求解法
7.8 二阶抛物型方程的差分求解法
7.9 二阶抛物型方程的线上求解法
7.10 二阶双曲型方程的特征线求解法浙江大学 实用数值计算方法 2
7.1 偏微分方程的一般介绍
Partial Differential Equations(PDEs)
自变量数 至少 2个
,,,,
0,,,,,,,,
,,
2
2
2
yx
u
u
x
u
u
y
u
u
x
u
u
uuuuuuyxF
yxuu
xyxxyx
yyxyxxyx
阶数 方程中导数的最高阶数
0
0
0
3
3
yyyx
yxx
yx
uu
uu
ubu
三阶二阶一阶
性态以一阶方程为例
0 cubua yx
浙江大学 实用数值计算方法 3
7.1
yyxyxxyx
yx
yxyyxyxx
yx
yx
yx
yx
uuuuuuyx
uuuyx
yx
CBA
FuEuDuCuBuA
uu
xuuu
buu
uuuyx
uyx
yx
yxccyxbbyxaa
cba
,,,,,,,:
,,,,:
,:
:,:,:
0:::
0
0
0
,,,,
,,
,
,,,,,
,,
2
2
：
对二阶方程例为常数非线性拟线性线性
N onlinear
rQ uasilinea
L inear
非线性拟线性线性非线性拟线性线性浙江大学 实用数值计算方法 4
7.1
类型一阶栓区型方程
xt
yx
uvu
cubua
0流动方程
Advection Equation(AE)
二阶线性方程
0 GCuBuAu yyxyxx
e llip t icACB 椭圆型042
方程，传热方程方程
Lap la c euu
P o is s o nyxfuu
yyxx
yyxx
0
,
p a r a b o licACB 抛物型042
方程扩散方程
B u r g e ruKuuu
uu
yyyx
yyx
h y p e r b o licACB 双曲型042
波动方程yyxx uu?
浙江大学 实用数值计算方法 5
7.1
求解方法有限差分法 Method of Finite Differences (MFD)
特征线法 Method of Characteristics (MOC)
线上求解法 Method of Lines (MOL)
有限元素法 Method of Finite Elements (MFE)
加权残差法 Method of Weighled Residuals (MWR)
问题收敛性 Convergence
当采取的步骤趋于无限时，数值结果是否趋于理论值？
稳定性 Stability
在某一步引入的误差，经多步数值计算后，会扩大或抑制？
dx
dy
x
y?
k
k
ey
y
ey
y?
1
1
浙江大学 实用数值计算方法 6
7.2 一阶双曲型方程的差分求解法
0 xuvtu
或称流动方程 Advective
Advection Equation (AE)
v为流速因子该方程的介折解
vtxfu
dw
du
x
w
dw
du
x
u
dw
du
v
t
w
dw
du
t
u
vtxwwfu
,因为求具体解时需要提供 2个辅助条件
vtxftxu
vtxftxu
00
00
,
,
浙江大学 实用数值计算方法 7
7.2
assuming the forcing function is a Rump
The solution of is shown below.
sW
Ws
s
W
W
Wf
,0.0
0,1
0,0.1
0 xt vuu
wf
W0s?
0.1
0.0
u
t
x
0tt?
0xx?
jx
nt
ntt?
jxx? tvx
txu j,
txu,0
0,txu
ntxu,
图 7.1 Propagation of the Wave Front
浙江大学 实用数值计算方法 8
7.2.1 最简单的差分化格式构想
t
x0x 1?jx jx 1?jxx?
1?nt
nt
0t
t?
tbxtu
xaxtu
x
uv
t
u
0
0
,
,xjxx j 0
tnttn 0
211
,
1
,
2
xO
x
uu
x
u
tO
t
uu
t
u
n
j
n
j
nj
n
j
n
j
nj
x
uuv
t
uu njnjnjnj
2
111
njnjnjnj uuxtvuu 111 2 1nn
图 7.2
浙江大学 实用数值计算方法 9
7.2.1
以上方法称为 时间镶嵌空间中心 的差分表达
Forward Time Centered Space
FTCS represetation
实际上这个方法不能用,不稳定的方法
Unstable Method
考虑数据误差 r
由于原方程为线性，故误差的传播关系是与原方程完全相同的差分方程差分方程独立解的一般形式
Independent Solutions of Difference Equations
njnjnjnjnjnjnjnj ruruxtvruru 111111 2
njnjnjnj rrxtvrr 111 2
xjkir
xjkir
nn
j
nn
j
e x p
e x p
11?
F a c t o rio nA m p lif ic a t称为放大因子?
是否上升序列决定,,,,,,1110 njnjnjjj rrrrr
浙江大学 实用数值计算方法 10
7.2.1
线性微分方程的解
xFQydxdyPdx yd22
xFQyDyPyD2
应为补充解和特殊解之和 补充解系由下式求出
0
0
0
21
2
2
yPDPD
yQDPD
QyDyPyD
e q u a t io na u x ilia r y
QPxxPP 的两个根为方程式 0,221
00 21 yPDyPD
补充解系由两个 独立解组成
xPAyxPAy 222111 e x pe x p
为复数时当 iPQP
yPxPPADy
4
e x p
2
111111
pyyyy 21
决定。由解，特殊解解组成的补充为微分方程的两个独立
xFy
yy
p
21?
浙江大学 实用数值计算方法 11
7.2.1
差分方程的解可用算符运算方法 Operator Calculus 导出差分算符 Difference Operator
jjjjj jjj xfxfyyy
xxx
11
1
jjjjj Eyyyyy 11
,1E
2
2
21
1
jj
jj
jj
yyE
yEy
yEy
它和微分算符一样，是一种线性算符用于线性二阶差分方程
xFQyPyy nnn 12
和微分方程类似，它的 补充解 可由下式得到
0
0
0
21
2
12
n
n
nnn
yPEPE
yQPEE
QyPyy
浙江大学 实用数值计算方法 12
7.2.1
故补充系由两个独立解组成
（ Independent Solutions）
0
0
2
1
n
n
yPE
yPE
两个独立解为
nnnn PAyPAy 2211 21
nnnnn yPEyyPPAy,11因为：
nnn PAPAy 2211
ixiPQP e x p42时，当
xniAy nn e x p? 差分方程的一个独立解 （ Eigenmode）
nn
n
n
n
nn
n
PyEy
in xAixpy
xniAyEy
aa
ai
aa
aiaa
iaia
n
aa
aa
e x pe x p
1e x p
!5!3!4!2
1
!4!3!2
1e x p
!!2
1e x p
1
1
5342
432
2
浙江大学 实用数值计算方法 13
7.2.1
差分方程独立解的一般形式
in xAy nn e x p
用于本题的情况
xik jr
xik jr
nn
j
nn
j
e x p
e x p
11?
乘以随时间坐标的推进不断njr
F a c t o rio nA m p lif ic a t称为“放大因子”?
时，为不稳定。当 1
将独立解代入差分表达式
njnjnjnj rrx tvrr 111 2
得到
xki
x
tv
xikxik
x
tv
xik j
xjikxjik
x
tv
xjikxjik
x
tv
xik j
nn
n
s i n2
2
e x pe x p
2
e x p
1e x p1e x p
2
1
1e x p1e x p
2
e x p
1
浙江大学 实用数值计算方法 14
7.2.2 差分格式的改进
的条件是使相应的放大因子为：
格式为：
得到的差分而代之以不用
1
s inc os
22
1
2
1
2
1111
1
11
i
Xk
x
tv
iXk
uu
x
tv
uuu
uu
u
M e t hodL ax
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
1xtv
Courant Condition
t
x0x 1?jx 1?jxjx
nt
0t
1?nt
I
R?
tk
1xtv?
图 7.3
图 7.4
浙江大学 实用数值计算方法 15
7.2.2
Courant 条件的物理意义波形传递系沿 x=vt线
t节点的选取
当节点取在线上：
当节点取在线外：
当节点取在线内：
Lax差分格式也写成以下形式可以看成为以下偏微分方程的 FTCS差分式
1,2 xtvvxt
不稳定1,3 xtvvxt
稳定1,1 xtvvxt
t
uuu
x
uuv
t
uu njnjnjnjnjnjnj 11111 2
2
1
2
2
22
x
u
t
x
x
uv
t
u
dissipative term 耗散项
Numerical Viscosity 数值黏度
t
x1?jx jx 1?jx
x?
nt
1t?
3t?
2t?
图 7.5
浙江大学 实用数值计算方法 16
7.3 一阶双曲型方程的特征线求解法
Method of Characteristics (MOC)
dy
y
u
dx
x
u
du
yxu
A
y
u
BC
x
u
yxCBA
的全微分线性的函数，可以是
,
,,,
dy
y
udx
A
y
uBC
du
CyuBxuA
0 Ad uC d xyuB d xAd y
这是原方程的转换方程，它们的解相同。
为原方程的特征线方程0 B d xAd y
yxFABdxdy,
在特征线上，满足 的为解。0 A duC dx
浙江大学 实用数值计算方法 17
7.3.1 Method of Characteristics
(MOC)
yxFdxdy,?
yxGdxdU,?
,,,,,,,
,,,,,,,
,,,,,,,
,
,,,,,,,
,
,,,,
10
11110
00100
01000
10
0
njnn
j
j
n
n
yxuyxuyxu
yxuyxuyxu
yxuyxuyxu
yxG
dx
du
yxuyxuyxu
yxF
dx
dy
yyyy
xx
yx
积分，可以得到再对条件轨迹。根据给定的初始线的积分可以得到一族特征方程作为初值，对常微分别由时，分从平面上在
y
x
ny
1y
0y
y?
x?0,0u
yxu,
0x jx1x
图 7.6
浙江大学 实用数值计算方法 18
7.3.1 Method of Characteristics (MOC)
yxFdxdy,?
yxGdxdU,?
y
x
ky0
01y
00y
y?
x?0,0u
yxu,?
0x jx1x
ky1
11y
10y
jky
0jy
1jy
kL
0L
1L
,2,1,0,
,
,,,,001000
kL
yxF
dx
dy
yyyxxyx
k
k
到一族特征线的轨迹积分可以得程作为初值，对常微分方时，分别由平面上，从在
,,,,,,,,,:
,,,,,,,,,:
,,,,,,,,,:
221100
12121110101
02021010000
jkjkkkk
jj
jj
yxyxyxyxL
yxyxyxyxL
yxyxyxyxL
图 7.7
浙江大学 实用数值计算方法 19
7.3.1
故积分上式，可以得到
，满足常微分方程因为在每一条特征线上可以计算出根据给定的初始条件，
yxG
dx
du
yxuyxuyxu
yxu
k
,
,,,,,,,
,,
00010000
0
,,,,,,,
,,,,,,,
,,,,,,,
1100
1111010
0101000
jkjkk
jj
jj
yxuyxuyxu
yxuyxuyxu
yxuyxuyxu
进行插值计算。方向上对函数值则还应进行上的解，方向上也规则的离散点此，若需要知道上。因方向不规则的离散点是定义在方法求得的离散解由此可知，用
uy
y
yxy
yxuM O C
jkj
jkj
,
,
浙江大学 实用数值计算方法 20
7.4 一阶双曲型方程的线上求解法
Method of Lines (MOL)
有限差分法：偏微分方程完全离散成为一组差分方程用线性代数方程组求解线上求解法：偏微分方程部分离散成为一组常微分方程用常微分方程积分方法求解
CxuBtuA
ni
x
u
A
B
A
C
dt
du
nixx
tgtxu
xftxu
x
u
A
B
A
C
t
u
tx
i
i
i
,,1,0
,,1,0,
,
,
,
0
0
得到常微分方程组：
轴离散化，将浙江大学 实用数值计算方法 21
线上求解法 Method of Lines (MOL)
nitutxu i,,1,0,,视为把部分离散化方法
t
x0x ix nxx?
1t
0t
t?
线间距积分步长
。常微分方程组的数值解分方法计算然后用任何一种数值积各种差分，样条，的近似值。
导方法得到：可以用任何一种数值求为已知。或初始条件
x
u
tutxu
i
ii
00
,
nixuABACtu ii,,1,0
7.4
图 7.8
浙江大学 实用数值计算方法 22
7.5 二阶椭圆型方程的差分求解法
02
2
2
2
yuxu
称为稳态热传导方程，通式为
Dirichlet 问题
Neumann 问题算符称为方程称为
L a p la c e
x
L a p la c eu
k k
2
2
2
2 0
02
2
2
2
yuxu
02
2
2
2
yuxu
yfyxu
yfyxu
yfyxu
yfyxu
n
m
41
301
2
10
,
,
,
,
ygxu
ygxu
ygxu
ygxu
n
m
y
y
x
x
4
3
2
1
0
0
浙江大学 实用数值计算方法 23
y
x
ny
0y
jy
0x ix mx
xfyxu n 4,?
xfyxu 30,?
02222 yuxu
yxQ, u(x m
,y)=
f 2(y
)
u(x
0,y
)=
f 1(y
)
xgyu
y
3
0
xgyu
ny
4
xgyu
mx
2
xg
yu x 10
Laplace 方程 的 Dirichlet 边界条件 和 Neumann 边界条件 和 Poisson 方程方程Po isson
yxyuxu,2222
边界条件也需 4个，有 3类给定方法
Dirichlet 边 界 条 件
Neumann 边 界 条 件
混合 边 界 条 件
7.5
图 7.9
浙江大学 实用数值计算方法 24
7.5.1 Laplace算符的差分表达
022222 yuxuu
简化表示为级数导出。精度可由二阶导数的差商表示及
2
4
''
2
4
4
3
'''
2
''
'
4
4
3
'''
2
''
'
12
2
2462
2462
h
f
xf
h
hxfxfhxf
h
f
h
xf
h
xf
hxfxfhxf
h
f
h
xf
h
xf
hxfxfhxf
T ay lor
i
iii
iii
iii
iii
iii
22 11'' 2 hOh ffff iiii
用于 Laplace算符
2
11
2
112
,,2,
,,2,
,
y
yxuyxuyxu
x
yxuyxuyxu
yxu
jijiji
jijiji
ii
并用简化表示当,hyx
jijijijijiij uuuuuhu,1,1,,1,122 41
浙江大学 实用数值计算方法 25
1?ix ix 1?ix
1?jy
jy
1?jy
1,?jiu
jiu,jiu,1?jiu,1?
1,?jiu
0
11
161
11
1
0
1
141
1
1
,
2
2
2
2
2
2
2
2
2
2
ijk
ij
u
h
z
u
y
u
x
u
u
u
h
yxu
L a pl ac e
在三维空间内方程的差分表达
7.5.1
图 7.10
浙江大学 实用数值计算方法 26
例,Laplace 方程的 Dirichlet 边界问题
co0
co0co0co0
co0 co0 co0
co100
y
x
cm10
cm20
5?yh
5?xh
解得或出一个线性方程对每个内部节点均可写
1 0 0
0
0
410
141
014
04001 0 0
0400
04000
3
2
1
32
231
12
T
T
T
TT
TTT
TT
cT
cT
cT
o
o
o
786.26
143.7
786.1
3
2
1
7.5.1
图 7.11
浙江大学 实用数值计算方法 27
为了提高精度需要加密网络
1 8 15 5 12
2 9 16 6 13
3 10 17 7 14
4 11 18 1 8 15
5 12 2 9
6 13 3 10
7 14 4 11
987654321 151413121110
9
8
7
6
5
4
3
2
1
100000141000001
100000141000001
10000004100000
1000001410000
100000141000
10000014100
1000001410
100000141
10000014
5,15121 非零元素数带宽ha
7.5.1
图 7.12
浙江大学 实用数值计算方法 28
Laplace 方程 Dirichlet边界问题的差分求解
消去法
直接迭代 Liebmann 方法
相继松弛 S.O.R,方法
交替方向 A.D.I.方法
04,,1,1,,1,12 jijijijijiij uuuuuu
411,1,,1,1, kjijijijikji uuuuu
1,1,,
1,
11,1,,1,1
1,4
kjikjikji
kji
kjijijiji
kji
uRuu
uuuuuuR
1
,1,1,
,,1,1,
1
,
1
,1,1,
,,1,1
1
,,
2
4
2
4
2
4
2
4
k
jijiji
k
jijiji
k
ji
k
ji
k
jijiji
k
jijiji
k
ji
k
ji
uuu
P
uuu
P
uu
uuu
P
uuu
P
uu
7.5.1
浙江大学 实用数值计算方法 29
7.6 二阶椭圆型方程的有限元素法求
Method of Finite Elements (MFE)
以 Laplace 方程的 Dirichlet 问题为例根据变分原则
Variational
Principles
等价性定理以上方程的解将使以下泛函为最小。
yxyxgyxu
Dyx
y
u
x
u
,,,,
,,02
2
2
2
d x d y
y
u
x
uuJ
D
22
2
1
S
D
me?
D
图 7.13
浙江大学 实用数值计算方法 30
7.6
将 D进行剖分，常用的是三角剖分法对任何一个元素用二原线性函数近似在三个顶点上可得到
SDeD
m
m,
yaaayxW e 321,
kkk
jjj
iii
Wyaxaa
Wyaxaa
Wyaxaa
321
321
321
kk
jj
ii
kk
jj
ii
kkk
jjj
iii
wx
wx
wx
e
a
yw
yw
yw
e
a
yxW
yxW
yxW
e
a
1
1
1
2
1
1
1
1
2
1
2
1
3
21
其中
kk
jj
ii
yx
yx
yx
e
1
1
1
2
浙江大学 实用数值计算方法 31
7.6
X
Y
u
ix
jx
kx iy
ky
jy
yxu,
e
yxW,
i
j
k
e?
Ui=Wi
Uk=Wk
Uj=Wj
图 7.14
浙江大学 实用数值计算方法 32
7.6
所以其中既然顶点坐标均为规定，所以并有
kkkk
jjjj
iiii
e
wydxcb
wydxcb
wydxcb
e
yxW
2
1
,
ijkjikijjik
kijikjkijkj
jkikjijkkji
xxdyycyxyxb
xxdyycyxyxb
xxdyycyxyxb
,,
,,
,,
kjie WWWyxW,,,
kkjjii
e
e
y
kkjjii
e
e
x
wdwdwd
ey
W
yxW
wcwcwc
ex
W
yxW
2
1
,
2
1
,
浙江大学 实用数值计算方法 33
7.6
使泛函最小的问题，即对近似为对求极值，或因此得到：
可解得
e e
e
y
e
x
yx
dx dyww
dx dyuuuJ
22
22
2
1
2
1
dx dy
e
wdwdwd
e
wcwcwc
wJ
kkjjii
e e
kkjjii
2
2
2
22
1
nmwJW
m
,,2,1,0
rWA
nmuw mm,,2,1,
n为内部节点数边界上的 W为给定浙江大学 实用数值计算方法 34
7.6
2
2
121
1
,
,
,,c os,c os,
,,,
,
,
,,,,
yx
yxh
yxuyxh
y
u
yxq
x
u
yxp
yxyxgyxu
Dyx
yxf
yxuyxr
y
u
yxq
yx
u
yxp
x
对于更为一般性的情况需要极小化的泛函将是也可剖分为有限个元素后求解
2
2
12
2
22
2
1
,
,,,
2
1
dsuhuhdx d yuyxf
uyxr
y
u
yxq
x
u
yxpuJ
y
x
D
1?
2?
sx
syS
normal
genttan
1?2
图 7.15
浙江大学 实用数值计算方法 35
7.8 二阶抛物型方程的差分求解法动态扩散方程对于一维空间用差商代替微商，可以有各种选择，例如所以有需要另有更方便的方法
t
C
DC?
12
xgtxC
xgtxC
xgtxC
x
C
D
t
C
n 3
20
10
2
2
,
,
,
,
iix
x
jijii
tx
jjt
jiji
tx
xxh
h
CCC
t
C
tth
ht
CC
t
C
ji
ji
12
,1,1
2
2
1
1,1,
,
2
,
2
jijiji
x
tjiji CCChhDCC,1,,1
21,1,22
时的数值和 ji tt 1?
浙江大学 实用数值计算方法 36
7.8
显式方法得到或者：
则有：
22,1,,1
,
2
2
,1,
,
2
x
x
jijiji
tx
t
t
jiji
tx
hO
h
ccc
x
c
hO
h
cc
t
c
ji
ji
jijiji
x
tjiji ccc
Dh
hcc
,1,,12,1,2
jijijiji crccrc,,1,11,21
2
122
2 rDhhh
Dhr
tx
x
t 时，，当其中
jijiji ccc,1,11,21
0x ix 1?ix1?ix nx x
t
1?jt
jt
0t
图 7.16
浙江大学 实用数值计算方法 37
7.8
示例：
取得到的数值解与以下解析解比较
2.0,20
0,0
04.00,
119.0
2
2
tc
tc
xc
x
c
t
c
s e c2.67
1 1 9.0
42
s e c1 1 9.04
2
2
t
x
h
cmDcmh
1
2
1
2
20
12
s i n1200294.0e x p
12
10
s i n01175.0e x p
20
2
,
n
n
xn
tn
xn
tn
x
txC
饱和蒸汽
C2H5OH
32.0,20 cmmgtc304.00,cmmgxc?
s e c11904.0 2cmD?
cm2020 0 x
30.0,0 cmmgtc
空气
2cmA
图 7.17
浙江大学 实用数值计算方法 38
0
1
2
3
4
5
6
840 24201612
%c
cmx 4?
cmx 12?
Number of time steps
Analytical Solutions
Numerical Solutions
s t e pt im ep e rt
cmxs e t
cmD
x
t
Dr
s e c6.334
119.0
25.0
0.4
s e c119.0,25.0
2
2
2
Analytical versus Numerical Solutions
Diffusion Dynamics
r?0.25
7.8
图 7.18
浙江大学 实用数值计算方法 39
0
1
2
3
4
5
6
420 121086
%c
cmx 4?
cmx 12?
Number of time steps
Analytical Solutions
Numerical Solutions
s t e pt im ep e rt
cmxs e t
cmD
x
t
Dr
s e c2.674
119.0
5.0
0.4
s e c119.0,5.0
2
2
2
Analytical versus Numerical Solutions
Diffusion Dynamics
r?0.5
7.8
图 7.19
浙江大学 实用数值计算方法 40
7.8
显式法的稳定性分析所以
jijiji
ji
Wce
tx
,,,
,
时存在误差若在
时，有在令所以代入前式，并根据级数展开
2
1
0,
,
,
,
2
,
2
21
21
2
,1,
2
2
22
,
,,1
2
1
22
,
,,1
1,,,,1,1
,,1,11,
r
MDEe
D
t
x
r
t
xc
tww
x
tcx
x
c
xww
x
tcx
x
c
xww
T ayl or
wwrwwr
ereere
jj
i
jiji
ji
jiji
ji
jiji
jjjijiji
jijijiji
2
2
,,1,11,
,,
21
x
tcD
t
xct
ereere
ji
jijijiji
浙江大学 实用数值计算方法 41
7.8
因此有
tME
tMErrEE
j
jjj
2121
0100
11
1
2
ttMEtMjE
tMEtMEE
j
jjj
则若当所以时，在
0,0,
2
1
,00
11
00
txr
MtE
Et
jj
0
,,
,
2
2
,
2
2
jiji
ji
x
c
D
t
c
x
tc
D
t
xc
M
稳定条件为故算法为稳定即 01jE
2
1
2
x
tDr
