2022年第九章-偏微分方程差分方法.docx

精选学习资料 - - - - - - - - - 精品资料 欢迎下载第 9 章 偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程；由于变量的增多和区域的复杂性,求 偏微分方程的精确解一般是不行能的,常常采纳数值方法求方程的近似解；偏微分 方程的数值方法种类较多,最常用的方法是 差分方法 ；差分方法具有格式简洁,程 序易于实现,运算量小等优点,特殊适合于规章区域上偏微分方程的近似求解；本 章将以一些典型的偏微分方程为例,介绍差分方法的基本原理和详细实现方法；9.1 椭圆型方程边值问题的差分方法 9.1.1 差分方程的建立 最典型的椭圆型方程是 Poisson （泊松）方程u x 2 u 2 y 2 u 2 f x , y , x , y G（9.1 ）G 是 x, y 平面上的有界区域,其边界 为分段光滑的闭曲线；当 fx, y 0 时,方程（9.1 ）称为 Laplace 拉普拉斯）方程；椭圆型方程的定解条件主要有如下三种边 界条件第一边值条件uux ,yy x ,y （9.2 ）其次边值条件（9.3 ）x ,nkuu第三边值条件（9.4 ）n这里, n 表示 上单位外法向, x,y, x,y, x,y和 k x,y 都是已知的函数,k x,y 0；满意方程（ 9.1 ）和上述三种边值条件之一的光滑函数 u x,y 称为椭圆型方程边值问题的解；用差分方法求解偏微分方程,就是要求出精确解ux,y 在区域 G 的一些离散节点（ xi,yi）上的近似值ui, j xi,y i）；差分方法的基本思想是,对求解区域G 做网格剖分,将偏微分方程在网格节点上离散化,导出精确解在网格节点上近似值所满足的差分方程,最终通过求解差分方程,通常为一个线性方程组,得到精确解在离散节点上的近似值；设 G=0<x<a, 0< y<b 为矩形区域,在x, y 平面上用两组平行直线x=ih1, i=0,1, , N1, h1=a/ N1y=jh 2, j=0,1, , N2, h2=b/ N2将 G 剖分为网格区域,见图9-1 ；h1,h2分别称为 x 方向和 y 方向的剖分步长,网格名师归纳总结 交点（xi,yi）称为剖分节点 （区域内节点集合记为Gh=（xi,yi）; （xi,yi）G ）,第 1 页,共 14 页网格线与边界 的交点称为边界点,边界点集合记为h；- - - - - - -精选学习资料 - - - - - - - - - 现在将微分方程（精品资料欢迎下载xi,9.1 ）在每一个内节点（xi,yi）上进行离散；在节点（yi）处,方程（ 9.1 ）为2 2x u2 x i , y i y u2 x i , y i f x i , y i , x i , y i G h（9.5 ）需进一步离散（9.5 ）中的二阶偏导数；为简化记号,简记节点（xi,yi）= i, j ,节点函数值 u（xi,yi）=u i, j ；利用一元函数的 Taylor 绽开公式,推得二阶偏导数的差商表达式2ui,j1ui,j,1j2 ui,jui,j,1j02 h 1x22 h 12ui,j1ui12 ui,jui1 02 h 2y2h2 2代入（ 9.5 ）式中,得到方程（9.1 ）在节点 i, j 处的离散形式1ui,1j2 h 22 ui,jui,1j1ui,j12 ui,jui,j1 2 h 12 h 2fi,jj02 h 1i,i,jGh02 h 1h2,就导出了u i, j 的近似值 ui,j所满其中ffxy i；舍去高阶小项i,2足的差分方程1ui,1j2 ui,jui,1j1ui,j12ui,jui,j1fi,j,Oi,jhGh（9.6 ）2 h 12 h 2在节点 i, j 处方程（ 9.6 ）靠近偏微分方程（9.1 ）的误差为2 h 12 2,它关于剖分步长是二阶的；这个误差称为差分方程靠近偏微分方程的 影响近似解的精度；截断误差 ,它的大小将在差分方程（ 9.6 ）中,每一个节点 i, j 处的方程仅涉及五个节点未知量ui, j,ui+1,j, ui-1 ,j,ui,j+1,ui,j-1,因此通常称（它简化为9.6 ）式为 五点差分格式 ,当 h1= h2=h 时,名师归纳总结 - - - - - - -第 2 页,共 14 页精选学习资料 - - - - - - - - - 1u i1,ju i,1ju i,j精品资料j14 u欢迎下载fi,j, i,jGh1u i,i,jh2差分方程（ 9.6 ）中,方程个数等于内节点总数,但未知量除内节点值 ui,j , i, j Gh外,仍包括边界点值；例如,点 1, j 处方程就含有边界点未知量 u0, j；因此,仍要利用给定的边值条件补充上边界点未知量的方程；对于第一边值条件式（9.2 ）,可直接取ui, j= （xi, yi）, i, j h（9.7 ）对于第三（ k=0 时为其次）边值条件式（9.4 ）,以左边界点 1, j 为例,见图9-2, 利用一阶差商公式u0 ,ju,0jh 1u ,1jOh 1n就得到边界点 0, j 处的差分方程u 0 , j u 1 , jk ,0 j u 0 , j r 0 , j（9.8 ）h 1联立差分方程（9.6 ）与（ 9.7 ）或（ 9.8 ）就形成了求解 Poisson 方程边值问题的差分方程组, 它实质上是一个关于未知量 ui, j 的线性代数方程组,可采纳第 2,3章介绍的方法进行求解；这个方程组的解就称为偏微分方程的 分解 ；考虑更一般形式的二阶椭圆型方程差分近似解 ,简称 差名师归纳总结 xAuyBuCuDuEufx ,y,x,yx iGj（9.9 ）第 3 页,共 14 页xyxy其中 A x, y Amin>0, B x, y Bmin >0, Ex, y 0；引进半节点1x i1 h ,12y i1y i1 h 22,利用一阶中心差商公式,在节点（i, j）处可有22xAui,j1Aui1,jAu i1,jO 2 h 1xh 1x2x2,1O 2 h 11A i1,jui,1ju i,jA i1,ju i,ju ih 1h 1h 122u i,ju i,1ju i,1jO 2 h 1x2h 1- - - - - - -精选学习资料 - - - - - - - - - 对yBu,u精品资料欢迎下载类似处理,就可推得求解方程（9.9 ）的差分方程yya ii1,jui,1ji,jai,1ju i,1ja i,j1 u i,j1a i,j1ui,j1a i,jui,j（9.10 ）f,j,G h其中a i,1,jh 12A i1,jh 1Ci,j2Bi,j1Bi,j1E i,j（9.11 ）22a i1,jh 12A i1,jh 1Ci,j22a ij1h 22Bi,j1h2Di,j22a i,j1h 22Bi,j1h2Di,j22a i,jh 12A i1 2,jA i1 2,jh 222明显,当系数函数 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 h 1 2h 2 2 阶；9.1.2 一般区域的边界条件处理前面已假设 G 为矩形区域, 现在考虑 G 为一般区域情形,这里主要涉及边界条件的处理；考虑 Poisson 方程第一边值问题名师归纳总结 uuxfx,y ,x ,y G（9.12 ）第 4 页,共 14 页,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 ）处,完全同矩形区域情形,可建立五点差分格式1ui,1j2ui,jui,1j1ui,j12ui,ju i,j1fi,j,i,jGh（9.13 ）2 h 12 h 2在方程（ 9.13 ）中,当（ i, j）点接近边界时,将 显现非正就内点上的未知量,因此必需补充非正就内点处的方程；如非正就内点恰好是边界点,如图 9-4 中 D点,就利用边界条件可取 uD= D 对于不是边界点的非正就内点,如图 9-4 中 B 点,一般可采纳如下两种处理方法；a. 直接转移法 . 取与点 B 距离最近的边界点 （如图 9-4 中 E 点）上的 u 的值作为 u B 的近似值 uB,即 uB=uE= E 直接转移法的优点是简洁易行,但精度较低,只为一阶近似；b. 线性插值法 . 取 B 点的两个相邻点（如图9-4 中边界点 A 和正就内点C 作为插值节点对uB 进行线性插值2 1hu Bx CxBuA xBxAuCOAAx Cxx Cx就得到点 B 处的方程uBh 1h 1Ah 1uC,xBxA线性插值法精度较高,为二阶近似；对每一个非正就内点进行上述处理,将所得到的方程与（9.13 ）式联立,就组成了方程个数与未知量个数相一样的线性代数方程组；求解此方程组就可得到一般名师归纳总结 - - - - - - -第 5 页,共 14 页精选学习资料 - - - - - - - - - 区域上边值问题（精品资料欢迎下载9.12 ）的差分近似解；对于一般区域上二阶椭圆型方程（9.9 ）的第一边值问题,可完全类似处理；其次、三边值条件的处理较为复杂,这里不再争论；9.2 抛物型方程的差分方法 本节介绍抛物型方程的差分方法,重点争论差分格式的构造和稳固性分析；9.2.1 一维问题 作为模型,考虑一维热传导的初边值问题ua2ufx ,t,0xl,0tT（9.14 ）x2tux,0x,0xl（9.15 ）2t,0tT（9.16 ）u0,tg 1t,ul,tg其中 a 是正常数,fx,t,x ,g 1t和g2t都是已知的连续的函数；现在争论求解问题（9.14 ）-9.18的差分方法；第一对求解区域G=0 x l, 0tT 进行网格剖分；取空间步长 作两族平行直线h=l/ N, 时间步长 =T/M, 其中 N,M 是正整数,xxjkjh,kj,1,0,Nttk,0 ,1,M将区域 G 剖分成矩形网格,见图9-5 ,网格交点（ xj,tk）称为节点；名师归纳总结 用差分方法求解初边值问题（9.14 ）-（ 9.16 ）就是要求出精确解u x, t 在每第 6 页,共 14 页个节点 （xj,t k）处的近似值ukuxj,tk；为简化记号, 简记节点 （xj,tk）=u j, k ；j利用一元函数的Taylor绽开公式,可推出以下差商表达式- - - - - - -精选学习资料 - - - - - - - - - u tj,kuj,k精品资料欢迎下载2kOh2（9.17 ）1u j,kO u tj,kuj,kuj,k1O （9.18 ）u tj,kuj,k1 u j,k1 O（9.19 ）22uj,kuj,1k2uj,kuj,1（9.20 ）x2h21. 古典显格式在区域 G 的内节点 j, k 处,利用公式（9.17 ）和（9.20 ）,可将偏微分方程 （9.14 ）离散为其中uj,kf1,tuj,kauj,1k2 u j,kuj,1kfjkO h 2h2fkx ik；舍去高阶小项Oh2,就得到节点近似值（差分解）u 所 kj满意的差分方程uk1ukauk12ukuk1fk（9.21 ）jjjjjh2j明显,在节点 j, k 处,差分方程（9.21 ）靠近偏微分方程 （9.14 ）的误差为 这个误差称为 截断误差 ,它反映了差分方程靠近偏微分方程的精度；现将（Oh2,9.21 ）式改写为便于运算的形式,并利用初边值条件（边界点方程,就得到9.15 ）与（ 9.16 ）补充上初始值和k 1 k k k ku j ru j 1 1 2 r u j ru j 1 f j0 j ,1 2 , , N ,1 k ,1,0 , M 1（9.22 ）u j x j , j ,1 ,2 , N 1k ku 0 g 1 t k , u N g 2 t k , k ,1,0 M其中 r a 2 称为 网比 ；h与时间相关问题差分方程的求解通常是按时间方向逐层进行的；对于差分方程k k 1（9.22 ）,当第 k 层节点值 ju 已知时,可直接运算出第 k+1 层节点值 ju ；这0样,从第 0 层已知值 u j x i 开头,就可逐层求出各时间层的节点值；差分方程（9.22 ）的求解运算是显式的,无须求解方程组,故称为古典显格式 ；此外,在式名师归纳总结 （9.22 ）中,每个内节点处方程仅涉及k 和 k+1 两层节点值,称这样的差分格式为第 7 页,共 14 页- - - - - - -精选学习资料 - - - - - - - - - 精品资料 欢迎下载双层格式 ；差分方程（ 9.22 ）可表示为矩阵形式uk1AukFk,k,1,0,M1（ 9.23 ）u0其中1r2 r1rr2 rA rukk u 1,uk1Tfk,rf12,2 r1rg2tkTNFkf1krg1tk,kfk N2Nx1,xN1T2. 古典隐格式 在区域 G 的内节点 j, k 处,利用公式（9.18 ）和（9.20 ）,可将偏微分方程 （9.14 ）离散为uj,ku j,k1auj,1k2 u j,kuj,1kfkOh2舍去高阶小项Oh22 hj,就得到如下差分方程ukuk1auk12ukuk1fjkjjjjj（9.24 ）h2它的截断误差为Oh2,靠近精度与古典显格式相同；改写（9.24 ）式为便于计算的形式,并补充上初始值与边界点方程,就得到名师归纳总结 uruk j121,2ruk jruk j1uk j1,1,fk jM（9.25 ）第 8 页,共 14 页Mj,1,N,1k0,1,0xj,j,12 ,Nk10jkkug 1 tk,ug2 tk,0N- - - - - - -精选学习资料 - - - - - - - - - 与古典显格式不同,在差分方程（精品资料欢迎下载k-1 层值uk1已知时,9.25 ）的求解中,当第j必需通过求解一个线性方程组才能求出第k 层值uk,所以称（ 9.25 ）式为 古典隐j格式 ,它也是双层格式；差分方程（ 9.25 ）的矩阵形式为k Buuk1Fk,k2,1,M（9.26 ）u0其中12 r1rrr2 rBr向量uk Fk,同（ 9.23 ）式中定义；从（r12 r9.26 ）式看到,古典隐格式在每一层运算时,都需求解一个三对角形线性方程组,可采纳追逐法求解；3.Crank-Nicolson格式（六点对称格式）O2利用一元函数Taylor绽开公式可得到如下等式uj,k1uj,k1 uj,kO2t22uj,k112uj,k2uj,k1 x222x2x2使用这两个公式,在j,k1点离散偏微分方程（9.14 ）,然后利用（ 9.20 ）式进2一步离散二阶偏导数,就可导出差分方程u k j 1 u k j a 12 u k j 1 1 2 uh k j2 1 u k j 1 1 u k j 1 2h u2 k j u k j 1 f j k 12（9.27 ）其截断误差为 O 2 h 2 ,在时间方向的靠近阶较显格式和隐格式高出一阶；这个差分格式称为 Crank-Nicolson 格式 ,有时也称为 六点对称格式 ,它明显是双层隐式格式；改写（ 9.27 ）式,并补充初始值和边界点方程得到名师归纳总结 - - - - - - -第 9 页,共 14 页精选学习资料 - - - - - - - - - ruk122 1r uk j1精品资料fk1欢迎下载（9.28 ）ruk1j1j1ruk12 1rukruk1221jjjjj,1,N,1k0 ,1,Mu0 jgxj,jk N,12 ,2,kNk10 ,1,Muk 01tk,ug t,它的矩阵形式为IB uk1IAukFk1,k,12,M（9.29 ）2u0,k01, ,M1在每层运算时,仍需求解一个三对角形方程组；4. Richardson格式O2h2阶的差分利用公式（ 9.19 ）和（ 9.20 ）,可导出另一个截断误差为方程u k j 1 u k j 1 u k j 1 2 u k j u k j 1 k a 2 f j 2 h称之为 Richardson 格式 ；可改写为u k j 1 u k j 1 2 r u k j 1 2 u k j u k j 1 2 f j k 9.32 k 1 k 这是一个三层显式差分格式；在逐层运算时,需用到 ju 和 u j 两层值才能得到 k+1k 1 0 1 层值 u j；这样,从第 0 层已知值 u j x j 开头,仍须补充上第一层值 u ,才能逐层运算下去；可采纳前述的双层格式运算 u ；1除上述四种差分格式外,仍可构造出很多靠近偏微分方程（9.14 ）的差分格式,但并不是每个差分格式都是可用的；一个有有用价值的差分格式应具有如下性质：（ 1）收敛性 ；对任意固定的节点（xj,tk）, 当剖分步长,h0时,差分解ukj应收敛到精确解u（xj,tk）；（ 2）稳固性 ；当某一时间层运算产生误差时,在以后各层的运算中,这些误差的传播积存是可掌握而不是无限增长的；理论上可以证明,在肯定条件下,稳固的差分格式必定是收敛的；因此,这里主要争论差分格式的稳固性；名师归纳总结 - - - - - - -第 10 页,共 14 页精选学习资料 - - - - - - - - - 作为例子,先考查精品资料欢迎下载k ju 是当运算过程中带有误差Richardson格式的稳固性；设k k k k时,按 Richardson 格式（9.30 ）得到的实际运算值；ju 是理论值, 误差 e j u j u j；假定右端项 jf k的运算是精确的,网比 r 1,就 je 满意 k2e kj 1e kj 1 e kj 1 2 e kj e kj 1 （9.31 ）设前 k-1 层运算时精确的,误差只是在第 k 层 0j 点发生,即k 1 k ke j ,0 e j 0 , e j 0 , j j 0；就利用（ 9.31 ）式可得到误差 的传播情形,见表 9-1 ；表 9-1 r=1/2 时 Richardson 格式的误差传播jj0-4 j0-3 j0-2 j0-1 j0 j0+1 j0+2 j0+3 j0+4 kk 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-260 641-1096 1311-1096 641-260 71从表中看出,误差是逐层无限增长的；表中的运算虽然是就网比 r 1进行的,2实际上对任何 r>0 都会产生类似现象,所以 Richardson 格式是不稳固的；利用误差传播图表方法考查差分格式的稳固性虽然直观明白,但只能就详细取定的 r 值进行,并且也不适用于隐式差分格式；9.2.2 差分格式的稳固性前节构造的几种双层差分格式都可以表示为如下的矩阵方程形式ukHuk1FkH=A, 隐格式 H=B-1 , 六点对称格式（9.32 ） -1u0其中 H 称为 传播矩阵 ；对于显格式H= I +BI+ A ；一般的三层格式也可以转化为双层格式；名师归纳总结 为了争论便利,设在初始层产生误差0 ,且假定右端项Fk的运算是精确的；用uk第 11 页,共 14 页- - - - - - -精选学习资料 - - - - - - - - - 表示起初始层存在误差精品资料欢迎下载uk满意方程0 时,由差分格式（9.32 ）得到的运算解,就ukHuk1Fk（9.33 ）（ 9.34 ）u00记误差向量kukuk,就k 满意方程kHk1,k,1 ,20为初始误差定义 9.1称差分格式 （9.32 ）是稳固的, 假如对任意初始误差0 ,误差向量k在某种范数下满意kC0,k0 ,00（9.35 ）其中 C 为与 h, 无关的常数；这个定义说明,当差分格式稳固时,它的误差传播是可掌握的；从（ 9.34 ）式递推得到kHk0,0kT因此,差分格式稳固的充分必要条件是定理 9.3 HkC,0kT（9.36 ）h, 无关的（稳固性必要条件）差分格式稳固的必要条件是存在与常数 M, 使谱半径定理 9.4 H1MHH*H*H（9.37 ）（稳固性充分条件）设 H 为正法规阵,即, 就（9.37 ）式也是差分格式稳固的充分条件；下面争论几种差分格式的稳固性；为便于争论,引进 N-1 阶矩阵0 11 0 1S10110这个特殊矩阵的特点值为名师归纳总结 - - - - - - -第 12 页,共 14 页精选学习资料 - - - - - - - - - S2cosj精品资料欢迎下载,j,1 ,2,N1（ 9.38 ）jN例 9-1 古典显格式此时 H=A=1-2rI+rS； 利用（9.38 ）式和三角函数公式,可求得 H 的特点值为j 1 4 r sin 2 j , j ,1 2 , , N 12 N为使稳固性条件式（9.39 ）成立,必需且只须 r 1；由于 H =A 为实对称矩阵,所2以古典显格式稳固的充分必要条件是网比 r 12例 9-2 古典隐格式 此时 H=B-1,B=1+2rI-rS ；利用（ 9.38 ）式可求得 H 的特征值为j 1 4 r sin 2 j 1 , j 1 2, , , N 12 N明显,对任意 r>0 ,条件（ 9.37 ）成立；留意,H=B-1 仍为实对称矩阵,所以古典隐格式对任何网比 r>0 都是稳固的,称为 肯定稳固 ；例 9-3 六点对称格式 此时 H=I+B-1I+A ,利用矩阵 A 和 B 的特点值可得到矩阵 H 的特点值为j24rsin2jN,j,12,N1AB=BA,就可2 j24r