第九章偏微分方程差分方法汇总Word文档格式.docx
《第九章偏微分方程差分方法汇总Word文档格式.docx》由会员分享，可在线阅读，更多相关《第九章偏微分方程差分方法汇总Word文档格式.docx（19页珍藏版）》请在冰点文库上搜索。
将G剖分为网格区域，见图9-1。
h1，h2分别称为x方向和y方向的剖分步长，网格交点（xi，yi）称为剖分节点（区域内节点集合记为Gh={（xi，yi）;
（xi，yi）∈G}）,网格线与边界Γ的交点称为边界点，边界点集合记为Γh。
现在将微分方程（9.1）在每一个内节点（xi，yi）上进行离散。
在节点（xi，yi）处，方程（9.1）为
[u2（xi,yi）u2（xi,yi）]f（xi,yi）,（xi,yi）Gh（9.5）
需进一步离散（9.5）中的二阶偏导数。
为简化记号，简记节点（xi，yi）=（i,j），节
2u
2（i,j）x
2（i,j）y
点函数值u（xi，yi）=u（i,j）。
利用一元函数的Taylor展开公式，推得二阶偏导数的差商表达式
12[u（i1,j）2u（i,j）u（i1,j）]0（h12）h12
12[u（i,j1）2u（i,j）u（i,j1）]0（h22）h2
代入（9.5）式中，得到方程（9.1）在节点（i,j）处的离散形式
11
12[u（i1,j）2u（i,j）u（i1,j）]12[u（i,j1）2u（i,j）u（i,j1）]h1h2
fi,j0（h12h22）,（i,j）Gh
其中fi,jf（xi,yi）。
舍去高阶小项0（h12h22），就导出了u（i,j）的近似值ui,j所满
足的差分方程
1
h12[ui1,j2ui,ju
i1,j]h122[ui,j12ui,jui,j1]fi,j,（i,j）Gh（9.6）
在节点（i,j）处方程（9.6）逼近偏微分方程（9.1）的误差为O（h12h22），它关于剖
分步长是二阶的。
这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。
在差分方程（9.6）中，每一个节点（i,j）处的方程仅涉及五个节点未知量ui,j，ui+1,j，ui-1,j，ui,j+1，ui,j-1，因此通常称（9.6）式为五点差分格式，当h1=h2=h时，它简化为
h2[ui1,jui1,jui,j1ui,j14ui,j]fi,j,（i,j）Gh
差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值ui,j,（i,j）∈
Gh外，还包括边界点值。
例如，点（1,j）处方程就含有边界点未知量u0,j。
因此，还
要利用给定的边值条件补充上边界点未知量的方程。
对于第一边值条件式（9.2），可直接取对于第三（k=0时为第二）边值条件式（以左边界点（1,j）为例，见图9-2,利用一阶差商公式
uu（0,j）u（1,j）
（0,j）O（h1）
nh1
则得到边界点（0,j）处的差分方程
u0,ju1,j
k0,ju0,jr0,j
h1
ui,j=α（xi，yi），（i,j）∈Γh（9.7）
9.4），
9.8）
联立差分方程（9.6）与（9.7）或（9.8）就形成了求解Poisson方程边值问题的差分方程组，它实质上是一个关于未知量{ui,j}的线性代数方程组，可采用第2，3章介绍的方法进行求解。
这个方程组的解就称为偏微分方程的差分近似解，简称差分解。
考虑更一般形式的二阶椭圆型方程
uuuu
[（A）（B）CDEu]f（x,y）,（x,y）G（9.9）xxyyxy
≥0。
引进半节点x1xi1h1,
i12i21,
其中A（x,y）≥Amin>
0,B（x,y）≥Bmin>
0,E（x,y）
u1u1u12
x（Axu）（i,j）h11[（Aux）（i12,j）（Aux）（i12,j）]O（h12）
u（i1,j）u（i,j）Au（i,j）u（i1,j）]O（h2）h1i2,jh1
O（h12）
ux（i,j）u（i1,j）2hu（i1,j）
x2h1
uu
对（B）,类似处理，就可推得求解方程（9.9）的差分方程yyy
9.10）
[ai1,jui1,jai1,jui1,jai,j1ui,j1ai,j1ui,j1ai,jui,j]f（i,j）,（i,j）Gh
其中
显然，当系数函数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）的截断误差为O（h12h22）阶。
9.1.2一般区域的边界条件处理
前面已假设G为矩形区域，现在考虑G为一般区域情形，这里主要涉及边界条件的处理。
考虑Poisson方程第一边值问题
9.12）
uf（x,y）,（x,y）Gu（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）属于GG，则称其为正则内点，见图9-3中打“。
”号者；
如果一个节点（i,j）
属于G且不为正则内点，则称其为非正则内点，见图9-3中打“.”号者。
记正则内点集合为Gh，非正则内点集合为h。
显然，当G为矩形区域时，
GhGh,hh成立。
在正则内点（i,j）处，完全同矩形区域情形，可建立五点差分格式
u（B）xCxBu（A）xBxAu（C）O（h12）xCxAxCxA
则得到点B处的方程
线性插值法精度较高，为二阶近似。
对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。
求解此方程组就可得到一般
区域上边值问题（9.12）的差分近似解。
对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。
第二、三边值条件的处理较为复杂，这里不再讨论。
9.2抛物型方程的差分方法
本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。
9.2.1一维问题
作为模型，考虑一维热传导的初边值问题
9.15）
9.16）
u（x,0）（x）,0xl
u（0,t）g1（t）,u（l,t）g2（t）,0tT
其中a是正常数，
f（x,t）,（x）,g1（t）和g2（t）都是已知的连续的函数。
现在讨论求解问题（9.14）-（9.18）的差分方法。
首先对求解区域G={0≤x≤l,0≤t≤T}进行网格剖分。
取空间步长h=l/N,时间步长τ=T/M,其中N,M是正整数，作两族平行直线
xxjjh,j0,1,,N
ttkk,k0,1,M
将区域G剖分成矩形网格，见图9-5，网格交点（xj，tk）称为节点。
用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解u（x,t）在每k
个节点（xj，tk）处的近似值ujku（xj,tk）。
为简化记号，简记节点（xj，tk）=u（j,k）。
利用一元函数的Taylor展开公式，可推出下列差商表达式
u（j,k）u（j,k1）u（j,k）O（）
t
u（j,k）u（j,k）u（j,k1）O（）u（j,k）u（j,k1）u（j,k1）O
（2）t2
2u（j,k）u（j1,k）2u（j,k）u（j1,k）O（h2）x2h2
1.
9.17）
9.18）
9.19）
（9.20）
古典显格式
在区域G的内节点（j,k）处，利用公式（9.17）和（9.20），可将偏微分方程（9.14）离散为
u（j,k1）u（j,k）au（j1,k）2u（j,k）u（j1,k）fkO（h2）
h2j
其中fjkf（xi,tk）。
舍去高阶小项O（h2），就得到节点近似值（差分解）ukj所满足的差分方程
显然，在节点（j,k）处，差分方程（9.21）逼近偏微分方程（9.14）的误差为O（h2），
这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。
现将（9.21）
式改写为便于计算的形式，并利用初边值条件（9.15）与（9.16）补充上初始值和边界点方程，则得到
k1kkkk
ujruj1（12r）ujruj1fj
9.22）
j1,2,,N1,k0,1,,M1u0j（xj）,j1,2,,N1
kk
u0kg1（tk）,uNkg2（tk）,k0,1,M
与时间相关问题差分方程的求解通常是按时间方向逐层进行的。
对于差分方程
kk1
9.22），当第k层节点值{ukj}已知时，可直接计算出第k+1层节点值{ukj1}。
这
样，从第0层已知值u0j（xi）开始，就可逐层求出各时间层的节点值。
差分方程
（9.22）的求解计算是显式的，无须求解方程组，故称为古典显格式。
此外，在式
（9.22）中，每个内节点处方程仅涉及k和k+1两层节点值，称这样的差分格式为
双层格式。
差分方程（9.22）可表示为矩阵形式
12rrr12rr
A
r
r12r
kkkT
u（u1,,uN1）
Fk（f1krg1（tk）,f2k,,fNk2,fNk1rg2（tk））T（（x1）,,（xN1））T
2.古典隐格式
在区域G的内节点（j,k）处，利用公式（9.18）和（9.20），可将偏微分方程（9.14）离散为
u（j,k）u（j,k1）u（j1,k）2u（j,k）u（j1,k）k2
a2fjO（h）h2
舍去高阶小项O（h2），则得到如下差分方程
它的截断误差为O（h2），逼近精度与古典显格式相同。
改写（9.24）式为便于计算的形式，并补充上初始值与边界点方程，则得到
kkkk1kruj1（12r）ujruj1ujfj
9.25）
j1,2,,N1,k0,1,,M
u0j（xj）,j1,2,,N1
u0g1（tk）,uNg2（tk）,k0,1,M
与古典显格式不同，在差分方程（9.25）的求解中，当第k-1层值{ukj1}已知时，
必须通过求解一个线性方程组才能求出第k层值{ukj}，所以称（9.25）式为古典隐格式，它也是双层格式。
差分方程（9.25）的矩阵形式为
9.26）
Bukuk1Fk,k1,2,M
B
rr12r
向量uk,Fk,同（9.23）式中定义。
从（9.26）式看到，古典隐格式在每一层计算时，都需求解一个三对角形线性方程组，可采用追赶法求解。
3.Crank-Nicolson格式（六点对称格式）
利用一元函数Taylor展开公式可得到如下等式
ut（j,k12）u（j,k1）u（j,k）
其截断误差为O（2h2），在时间方向的逼近阶较显格式和隐格式高出一阶。
这个差分格式称为Crank-Nicolson格式，有时也称为六点对称格式，它显然是双层隐式格式。
改写（9.27）式，并补充初始值和边界点方程得到
方程
这是一个三层显式差分格式。
在逐层计算时，需用到ukj1和ukj两层值才能得到k+1
层值ukj1。
这样，从第0层已知值u0j（xj）开始，还须补充上第一层值u1j，才能
逐层计算下去。
可采用前述的双层格式计算u1j。
除上述四种差分格式外，还可构造出许多逼近偏微分方程（9.14）的差分格式，但并不是每个差分格式都是可用的。
一个有实用价值的差分格式应具有如下性质：
（1）收敛性。
对任意固定的节点（xj，tk）,当剖分步长,h0时，差分解ukj应收敛到精确解u（xj，tk）。
（2）稳定性。
当某一时间层计算产生误差时，在以后各层的计算中，这些误差的传播积累是可控制而不是无限增长的。
理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。
因此，这里主要研究差分格式的稳定性。
作为例子，先考查Richardson格式的稳定性。
设ujk是当计算过程中带有误差
时，按Richardson格式（9.30）得到的实际计算值。
ukj是理论值，误差ejkujkukj。
假定右端项fjk的计算是精确的，网比r，则ejk满足
j2j
k1k1kkk
ejej（ej12ejej1）（9.31）
设前k-1层计算时精确的，误差只是在第k层j0点发生，即
k1kkej0,ej0,ej0,jj0。
则利用（9.31）式可得到误差的传播情况，见表9-1。
表9-1r=1/2时Richardson格式的误差传播
kj
j0-4
j0-3
j0-2
j0-1
j0
j0+1
j0+2
j0+3
j0+4
k
ε
k+1
-2ε
k+2
-4ε
7ε
k+3
-6ε
17ε
-24ε
-6ε
k+4
-8ε
31ε
-68ε
89ε
k+5
-10ε
49ε
-144ε
273ε
-388ε
k+6
71ε
-260ε
641ε
-1096ε
1311ε
-260ε
从表中看出，误差是逐层无限增长的。
表中的计算虽然是就网比r12进行的，
实际上对任何r>
0都会产生类似现象，所以Richardson格式是不稳定的。
利用误差传播图表方法考查差分格式的稳定性虽然直观明了，但只能就具体取定的r值进行，并且也不适用于隐式差分格式。
9.2.2差分格式的稳定性
前节构造的几种双层差分格式都可以表示为如下的矩阵方程形式
其中H称为传播矩阵。
对于显格式H=A,隐格式H=B-1,六点对称格式H=（I+B）（I+A）。
一般的三层格式也可以转化为双层格式。
为了讨论方便，设在初始层产生误差0，且假定右端项Fk的计算是精确的。
用uk
表示当初始层存在误差0时，由差分格式（9.32）得到的计算解，则
kk1k
uHuF
00
u
记误差向量kukuk，则k满足方程
kHk1,k1,2,
0为初始误差
定义9.1称差分格式（9.32）是稳定的，如果对任意初始误差0在某种范数下满足
kC0,k0,00
其中C为与h,τ无关的常数。
这个定义表明，当差分格式稳定时，它的误差传播是可控制的。
从（9.34）式递推得到
kk0T
H,0k因此，差分格式稳定的充分必要条件是
定理9.3（稳定性必要条件）差分格式稳定的必要条件是存在与常数M,使谱半径
（H）1M
定理9.4（稳定性充分条件）设H为正规矩阵，即HH*H*式也是差分格式稳定的充分条件。
下面讨论几种差分格式的稳定性。
为便于讨论，引进N-1阶矩阵01101
S
101
uk满足方程
（9.33）
（9.34）
，误差向量k
（9.35）
（9.36）
h,τ无关的
（9.37）
H,则（9.37）
这个特殊矩阵的特征值为
9.38）
Sj2cosj,j1,2,,N1jN
例9-1古典显格式此时H=A=（1-2r）I+rS。
利用（9.38）式和三角函数公式，可求得H的特征值为
2j
j14rsin2（2jN）,j1,2,,N1
j
[14rsin2（）]1,j1,2,,N1
2N
r>
0，条件（9.37）成立。
注意，H=B-1仍为实对称矩阵，所以古典隐
例9-2古典隐格式此时H=B-1征值为
显然，对任意格式对任何网比r>
0都是稳定的，称为绝对稳定。
例9-3六点对称格式此时H=（I+B）-1（I+A），利用矩阵A和B的特征值可得到矩阵H的特征值为
24rsin2（j）
j2N,j1,2,,N1j2j
24rsin2（）2N
则对任意r>
0,条件（9.37）成立。
由于A和B均为实对称矩阵，且AB=BA，则可验证H也是实对称矩阵。
所以六点对称格式是绝对稳定的。
习题九
9-1试用五点差分格式求解Poisson方程的边值问题
2216,（x,y）Gxy
u|0,（x,y）
其中G={-1<
x,y<
1}。
取步长h=0.5求解。
9-2试写出求解Laplace方程边值问题
22
2216,0x4,0y3
x2y2
u（0,y）y（3y）,u（4,y）0,0y3
u（x,0）sinx,u（x,3）0,0x4
4
的五点差分格式，取步长h=1,并写出差分方程的矩阵形式。
9-3试用古典显格式求解热传导方程定解问题
x2
0x1,0tT
u（x,0）4x（1x）,0x1u（0,t）u（1,t）0,0tT
只计算k=1,2两层上的差分解，取网比r=1/6,h=0.2
9-4用古典隐格式求解热传导方程定解问题
u（x,0）sinx,0x1u（0,t）u（1,t）0,0t0.3
取0.1,h0.2精确解为u（x,t）etsinx。
9-5导出求解热传导方程的差分格式
k1kkk1
ujujujuj1（k2kk）
（1）2（uj12ujuj1）
h2
的截断误差，并选取θ使其达到二阶。
升级会员 