第九章-偏微分方程差分方法_高等教育-微积分.pdf

《第九章-偏微分方程差分方法_高等教育-微积分.pdf》由会员分享，可在线阅读，更多相关《第九章-偏微分方程差分方法_高等教育-微积分.pdf（14页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精品资料 欢迎下载 第 9 章 偏微分方程的差分方法 含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。9.1 椭圆型方程边值问题的差分方法 9.1.1 差分方程的建立 最典型的椭圆型方程是 Poisson（泊松）方程 Gyxyxfyuxuu),(),()(2222 （9.1）G 是 x,y 平

2、面上的有界区域，其边界为分段光滑的闭曲线。当 f(x,y)0 时，方程（9.1）称为 Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件 第一边值条件 ),(yxu （9.2）第二边值条件 ),(yxnu （9.3）第三边值条件 ),()(yxkunu （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 的一些离散节点（xi，yi）上的近似值 u

3、i,j(xi，yi）。差分方法的基本思想是，对求解区域 G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。设 G=0 xa,0 y0,B(x,y)Bmin 0,E(x,y)0。引进半节点,12121hxxii,22121hyyii利用一阶中心差商公式，在节点（i,j）处可有)(2),1(),1(),()(),1(),(),(),1(1)(),21)(),21)(1),)(211211,211,211211hOhjiujiujixuhOhjiujiuAhjiujiuAhhOji

4、xuAjixuAhjixuAxjiji精品资料 欢迎下载 对yuyuBy),(类似处理，就可推得求解方程（9.9）的差分方程 hjijijijijijijijijijiGjijifuauauauaua),(),(,1,1,1,1,1,1,1,1 （9.10）其中 jijijijijijijijijijijijijijijijijijiEBBhAAhaDhBhaDhBhaChAhaChAha,21,21,22,21,2121,221,221,221,221,1,2121,1,1,2121,1)()()2()2()2()2(（9.11）显然，当系数函数A(x,y)=B(x,y)=1,C(x,y)=

5、D(x,y)=E(x,y)=0 时，椭圆型方程（9.9）就成为 Poisson方程（9.1），而差分方程（9.10）就成为差分方程（9.6）。容易看出，差分方程（9.10）的截断误差为)(2221hhO阶。9.1.2 一般区域的边界条件处理 前面已假设 G 为矩形区域，现在考虑 G 为一般区域情形，这里主要涉及边界条件的处理。考虑 Poisson 方程第一边值问题),(),(),(),(yxyxuGyxyxfu （9.12）其中 G 可为平面上一般区域，例如为曲边区域。仍然用两组平行直线：x=x0+ih1,y=y0+jh2,i,j=0,1,对区域 G 进行矩形网格剖分，见图 9-3。如果一个内

6、节点（i,j）的四个相邻节点（i+1,j），（i-1,j），（i,j+1）和（i,j-1）属于 GG，则称其为正则内点，见图 9-3 中打“。”号者；如果一个节点（i,j）属于G且不为正则内点，则称其为非正则内点，见图 9-3 中打“.”号者。记正则内点集合为hG，非正则内点集合为h。显然，当 G 为矩形区域时，精品资料 欢迎下载 hhhhGG,成立。在正则内点（i,j）处，完全同矩形区域情形，可建立五点差分格式 hjijijijijijijiGjifuuuhuuuh),(,2121,1,1,22,1,121 （9.13）在方程（9.13）中，当（i,j）点临近边界时，将出现非正则内点上的未知

7、量，因此必须补充非正则内点处的方程。若非正则内点恰好是边界点，如图 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)进行线性插值)()()()(21hOCuxxxxAuxxxxBuACABACBC 则得到点 B 处的方程

8、 ABCBxxuhAhhu,)(111 线性插值法精度较高，为二阶近似。对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般精品资料 欢迎下载 区域上边值问题（9.12）的差分近似解。对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。第二、三边值条件的处理较为复杂，这里不再讨论。9.2 抛物型方程的差分方法 本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。9.2.1 一维问题 作为模型，考虑一维热传导的初边值问题 Ttlxtxfxuatu0,0),(22 （9.14

9、）lxxxu0),()0,(（9.15）Tttgtlutgtu0),(),(),(),0(21 （9.16）其中 a 是正常数，)()(),(),(21tgtgxtxf和都是已知的连续的函数。现在讨论求解问题（9.14）-(9.18)的差分方法。首先对求解区域 G=0 xl,0tT进行网格剖分。取空间步长 h=l/N,时间步长=T/M,其中 N,M 是正整数，作两族平行直线 MkkttNjjhxxkj,1,0,1,0,将区域 G 剖分成矩形网格，见图 9-5，网格交点（xj，tk）称为节点。用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解 u(x,t)在每个节点（xj，tk）

10、处的近似值),(kjkjtxuu。为简化记号，简记节点（xj，tk）=u(j,k)。利用一元函数的 Taylor展开公式，可推出下列差商表达式 精品资料 欢迎下载 )(),()1,(),(Okjukjukjtu （9.17）)()1,(),(),(Okjukjukjtu （9.18）)(2)1,()1,(),(2Okjukjukjtu （9.19）)(),1(),(2),1(),(2222hOhkjukjukjukjxu （9.20）1.古典显格式 在区域 G 的内节点(j,k)处，利用公式（9.17）和（9.20），可将偏微分方程（9.14）离散为)(),1(),(2),1(),()1,(2

11、2hOfhkjukjukjuakjukjukj 其中),(kikjtxff。舍去高阶小项)(2hO，就得到节点近似值（差分解）kju所满足的差分方程 kjkjkjkjkjkjfhuuuauu21112 （9.21）显然，在节点(j,k)处，差分方程（9.21）逼近偏微分方程（9.14）的误差为)(2hO，这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。现将（9.21）式改写为便于计算的形式，并利用初边值条件（9.15）与（9.16）补充上初始值和边界点方程，则得到 MktgutguNjxuMkNjfruurruukkNkkjjkjkjkjkjkj,1,0),(),(1,2,1),(

12、1,1,0,1,2,1)21(2100111 （9.22）其中2har称为网比。与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程（9.22），当第 k 层节点值kju已知时，可直接计算出第 k+1 层节点值1kju。这样，从第 0 层已知值)(0ijxu开始，就可逐层求出各时间层的节点值。差分方程（9.22）的求解计算是显式的，无须求解方程组，故称为古典显格式。此外，在式（9.22）中，每个内节点处方程仅涉及 k 和 k+1 两层节点值，称这样的差分格式为精品资料 欢迎下载 双层格式。差分方程（9.22）可表示为矩阵形式 011,1,0,uMkFAuukkk （9.23）其

13、中 rrrrrrrrA212121 TNTkkNkNkkkkTkNkkxxtrgffftrgfFuuu)(,),()(,),(),(1121221111 2.古典隐格式 在区域 G 的内节点(j,k)处，利用公式（9.18）和（9.20），可将偏微分方程（9.14）离散为)(),1(),(2),1()1,(),(22hOfhkjukjukjuakjukjukj 舍去高阶小项)(2hO，则得到如下差分方程 kjkjkjkjkjkjfhuuuauu21112 （9.24）它的截断误差为)(2hO，逼近精度与古典显格式相同。改写（9.24）式为便于计算的形式，并补充上初始值与边界点方程，则得到 Mk

14、tgutguNjxuMkNjfuruurrukkNkkjjkjkjkjkjkj,1,0),(),(1,2,1),(,1,0,1,2,1)21(2100111 （9.25）精品资料 欢迎下载 与古典显格式不同，在差分方程（9.25）的求解中，当第 k-1 层值1kju已知时，必须通过求解一个线性方程组才能求出第 k层值kju，所以称（9.25）式为古典隐格式，它也是双层格式。差分方程（9.25）的矩阵形式为 01,2,1,uMkFuBukkk （9.26）其中 rrrrrrrrB212121 向量,kkFu同（9.23）式中定义。从（9.26）式看到，古典隐格式在每一层计算时，都需求解一个三对角

15、形线性方程组，可采用追赶法求解。3.Crank-Nicolson格式（六点对称格式）利用一元函数 Taylor 展开公式可得到如下等式 )()1,(),(21)21,()(),()1,()21,(22222222OkjxukjxukjxuOkjukjukjtu 使用这两个公式，在)21,(kj点离散偏微分方程（9.14），然后利用（9.20）式进一步离散二阶偏导数，则可导出差分方程 2121121111112221kjkjkjkjkjkjkjkjkjfhuuuhuuuauu （9.27）其截断误差为)(22hO，在时间方向的逼近阶较显格式和隐格式高出一阶。这个差分格式称为 Crank-Nico

16、lson格式，有时也称为六点对称格式，它显然是双层隐式格式。改写（9.27）式，并补充初始值和边界点方程得到 精品资料 欢迎下载 MktgutguNjxuMkNjfruurruruurrukkNkkjjkjkjkjkjkjkjkj,1,0),(),(1,2,1),(1,1,0,1,2,12)1(2)1(22100211111111 （9.28）它的矩阵形式为 1,1,0,2,1,)()(0211MkuMkFuAIuBIkkk （9.29）在每层计算时，仍需求解一个三对角形方程组。4.Richardson格式 利用公式（9.19）和（9.20），可导出另一个截断误差为)(22hO阶的差分方程 k

17、jkjkjkjkjkjfhuuuauu2111122 称之为 Richardson 格式。可改写为 kjkjkjkjkjkjfuuuruu2)2(21111 (9.32)这是一个三层显式差分格式。在逐层计算时，需用到1kju和kju两层值才能得到 k+1层值1kju。这样，从第 0 层已知值)(0jjxu开始，还须补充上第一层值1ju，才能逐层计算下去。可采用前述的双层格式计算1ju。除上述四种差分格式外，还可构造出许多逼近偏微分方程（9.14）的差分格式，但并不是每个差分格式都是可用的。一个有实用价值的差分格式应具有如下性质：（1）收敛性。对任意固定的节点（xj，tk）,当剖分步长0,h时，

18、差分解kju应收敛到精确解 u（xj，tk）。（2）稳定性。当某一时间层计算产生误差时，在以后各层的计算中，这些误差的传播积累是可控制而不是无限增长的。理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。因此，这里主要研究差分格式的稳定性。精品资料 欢迎下载 作为例子，先考查 Richardson格式的稳定性。设kju是当计算过程中带有误差时，按 Richardson 格式（9.30）得到的实际计算值。kju是理论值，误差kjkjkjuue。假定右端项kjf的计算是精确的，网比21r，则kje满足 )2(1111kjkjkjkjkjeeeee （9.31）设前 k-1 层计算时精确的，误

19、差只是在第 k层0j点发生，即 01,0,00jjeeekjkjkj。则利用（9.31）式可得到误差的传播情况，见表 9-1。表 9-1 r=1/2时 Richardson 格式的误差传播 j k j0-4 j0-3 j0-2 j0-1 j0 j0+1 j0+2 j0+3 j0+4 k 0 0 0 0 0 0 0 0 k+1 0 0 0 -2 0 0 0 k+2 0 0 -4 7 -4 0 0 k+3 0 -6 17 -24 17 -6 0 k+4 -8 31 -68 89 -68 31 -8 k+5-10 49 -144 273 -388 273 -144 49 -10 k+6 71 -26

20、0 641 -1096 1311 -1096 641 -260 71 从表中看出，误差是逐层无限增长的。表中的计算虽然是就网比21r进行的，实际上对任何 r0 都会产生类似现象，所以 Richardson 格式是不稳定的。利用误差传播图表方法考查差分格式的稳定性虽然直观明了，但只能就具体取定的 r 值进行，并且也不适用于隐式差分格式。9.2.2 差分格式的稳定性 前节构造的几种双层差分格式都可以表示为如下的矩阵方程形式 01uFHuukkk （9.32）其中 H 称为传播矩阵。对于显格式 H=A,隐格式 H=B-1,六点对称格式 H=(I+B)-1(I+A)。一般的三层格式也可以转化为双层格式

21、。为了讨论方便，设在初始层产生误差0，且假定右端项kF的计算是精确的。用ku精品资料 欢迎下载 表示当初始层存在误差0时，由差分格式（9.32）得到的计算解，则ku满足方程 001uFuHukkk （9.33）记误差向量kkkuu，则k满足方程 为初始误差01,2,1,kHkk （9.34）定义 9.1 称差分格式（9.32）是稳定的，如果对任意初始误差0，误差向量k在某种范数下满足 000,0,kCk （9.35）其中 C 为与 h,无关的常数。这个定义表明，当差分格式稳定时，它的误差传播是可控制的。从（9.34）式递推得到 TkHkk0,0 因此，差分格式稳定的充分必要条件是 TkCHk0

22、,（9.36）定理 9.3 （稳定性必要条件）差分格式稳定的必要条件是存在与 h,无关的常数 M,使谱半径 MH 1)(（9.37）定理 9.4 （稳定性充分条件）设 H 为正规矩阵，即HHHH*,则（9.37）式也是差分格式稳定的充分条件。下面讨论几种差分格式的稳定性。为便于讨论，引进 N-1 阶矩阵 0110110110S 这个特殊矩阵的特征值为 精品资料 欢迎下载 1,2,1,cos2NjNjSj （9.38）例 9-1 古典显格式 此时 H=A=(1-2r)I+rS。利用（9.38）式和三角函数公式，可求得 H 的特征值为 1,2,1),2(sin412NjNjrj 为使稳定性条件式（

23、9.39）成立，必须且只须21r。由于 H=A 为实对称矩阵，所以古典显格式稳定的充分必要条件是网比21r 例 9-2 古典隐格式 此时 H=B-1,B=(1+2r)I-rS。利用（9.38）式可求得 H 的特征值为 1,2,1,)2(sin41 12NjNjrj 显然，对任意 r0，条件（9.37）成立。注意，H=B-1仍为实对称矩阵，所以古典隐格式对任何网比 r0 都是稳定的，称为绝对稳定。例 9-3 六点对称格式 此时 H=(I+B)-1(I+A)，利用矩阵 A 和 B 的特征值可得到矩阵 H 的特征值为 1,2,1,)2(sin42)2(sin4222NjNjrNjrj 则对任意 r0

24、,条件（9.37）成立。由于 A 和 B 均为实对称矩阵，且 AB=BA，则可验证 H 也是实对称矩阵。所以六点对称格式是绝对稳定的。习题九 9-1 试用五点差分格式求解 Poisson 方程的边值问题 ),(,0|),(,162222yxuGyxyuxu 其中 G=-1x,y1。取步长 h=0.5 求解。9-2 试写出求解 Laplace 方程边值问题 精品资料 欢迎下载 40,0)3,(,4sin)0,(30,0),4(),3(),0(30,40,162222xxuxxuyyuyyyuyxyuxu 的五点差分格式，取步长 h=1,并写出差分方程的矩阵形式。9-3 试用古典显格式求解热传导方程定解问题 TttutuxxxxuTtxxutu0,0),1(),0(10),1(4)0,(0,10,22 只计算 k=1,2 两层上的差分解，取网比 r=1/6,h=0.2。9-4 用古典隐格式求解热传导方程定解问题 3.00,0),1(),0(10,sin)0,(3.00,10,22ttutuxxxutxxutu 取2.0,1.0h精确解为 xetxutsin),(2。9-5 导出求解热传导方程的差分格式 )2(1)1(11211kjkjkjkjkjkjkjuuuhuuuu 的截断误差，并选取使其达到二阶。