（计算数学专业论文）对流扩散方程的紧差分格式.pdf

，；、 厂谲i 丽1 1 1 1 1 1 111, 本人声明所呈交的学位论文是本人在导师指导下进行的研究工 作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地 方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含 为获得电子科技大学或其它教育机构的学位或证书而使用过的材料。 与我一同工作的同志对本研究所做的任何贡献均己在论文中作了明 确的说明并表示谢意。 签名：圣 曰期：乃p 年6 月p 曰 论文使用授权 本学位论文作者完全了解电子科技大学有关保留、使用学位论文 的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁 盘，允许论文被查阅和借阅。本人授权电子科技大学可以将学位论文 的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或 扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后应遵守此规定) 签名：墨丝导师签名： 日期：纠 t 摘要 摘要 对流扩散方程是偏微分方程中非常重要的一类基本运动方程，是用来描述黏 性流体的非线性方程的线性化模型方程。本文主要针对一维和二维的对流扩散方 程，根据已有的某些有限差分理论中的观点、理论和方法，分别推导了它们的高 阶精度的紧差分格式。 对于一维对流扩散方程来说，本文主要通过研究差分算子与微分算子的关系， 推导出一维对流扩散方程的四阶和六阶紧差分格式，并对紧差分格式做了误差分 析，数值实验结果验证了格式的有效性；然后，分析了三对角托普里兹矩阵的谱 和谱条件数；最后，本文利用r i c h a r d s o n 外推公式及陈一林推广公式对方程的数 值解做数值外推和误差修正，使其达到六阶精度，数值实验结果表明该种方法具 有很好的效果。 对于二维对流扩散方程来说，本文主要从一维对流扩散方程的二阶和四阶紧 差分格式入手，通过研究差分算子和微分算子的关系，推导出二阶精度的五点差 分格式和四阶精度的九点紧差分格式；然后，根据差分格式的计算模板，采用正 弦变换追赶法、b i c g 、b i c g s t a b 等算法进行数值求解，数值实验验证了格式的 有效性。 关键词：对流扩散方程，紧差分格式，差分算子，微分算子，外推技术 r _ _ - 。_ _ - _ _ _ _ _ - - 1 1 。一 i a b s t r a c t a b s t r a o t c o n v e c t i o n - d i f f u s i o ne q u a t i o n sa r ef u n d a m e n t a le q u a t i o n sf o rd y n a m i c sw h i r o c c u p ya ni m p o r t a n tp o s i t i o ni np a r t i a ld i f f e r e n t i a le q u a t i o n s ，a n da r et h e1 i n e a r i z e d m o d e lo ft h en o n l i n e a re q u a t i o n sc o m i n gf r o mt h ev i s c o u sf l u i dd y n a m i c s i nt h i s p a p e r , b a s e do ns o m eo ft h ef i n i t ed i f f e r e n c et h e o r ) rp o i n to fv i e w , t h e o r i e sa n d m e t h o d s ，t h eh i 曲a c c u r a c yc o m p a c td i f f e r e n c es c h e m e sf o ro n ea n dt w od i m e n s i o n a l c o n v e c t i o n d i f f u s i o ne q u a t i o n sa r ep r o p o s e dr e s p e c t i v e l y f o ro n e d i m e n s i o n a lc o n v e c t i o n d i f f u s i o n e q u a t i o n ，af o u r t ha n ds i x t ho r d e r c o m p a c td i f f e r e n c es c h e m ei sd e r i v e db ys m d 妒n gt h er e l a t i o nb e t w e e nt h ed i f f e r e n c e o p e r a t o ra n dt h ed i f f e r e n t i a lo p e r a t o r a n dt h ee r r o ra n a l y s i sf o rt h ec o m p a c td i f f e r e n c e s c h e m ei sa l s og i v e n t h en u m e r i c a le x p e r i m e n t s p r o v et h o s es c h e m e s e f f e c t i v e t h e n t h es p e c t r u ma n dc o n d i t i o nn u m b e ro ft r i d i a g o n a l t o e p l i t zm a t r i xi sa n a l y z e d s e p a r a t e l y f i n a l l y ，t h er i c h a r d s o na n dc h e n l i ne x p a n df o r m u l ai su s e dt oc o m p u t ea s i x t ho r d e rn u m e r i c a ls o l u t i o no ft h ee q u a t i o n t h en u m e r i c a le x p e r i m e n t ss h o wt h i s m e t h o dh a sa g o o de f f o r t f o rt w o - d i m e n s i o n a lc o n v e c t i o nd i f f u s i o ne q u a t i o n ，a c c o r d i n gt ot h es e c o n da n d f o u r t hc o m p a c td i f f e r e n c es c h e m eo fo n e - d i m e n s i o n a le q u a t i o n ，t h e f i v e p o i n ta n d n i n e - p o i n td i f f e r e n c es c h e m ei sd e r i v e db ys m d y i n gt h er e l a t i o nb e t w e e nt h ed i f f e r e n c e o p e r a t o ra n dt h ed i f f e r e n t i a lo p e r a t o r t h e nt h ea l g o r i t h m s ，s u c ha s ，c h a s i n gm e t h o do f s i n et r a n s f o r m ，b i c ga n db i c g s t a b ，a r eu s e dt ot h en u m e r i c a ls o l u t i o no fe q u a t i o n a c c o r d i n gt ot h ec a l c u l a t i o nt e m p l a t e so ft h ec o m p a c td i f f e r e n c es c h e m e s t h e n u m e r i c a le x p e r i m e n t sp r o v et h o s es c h e m e s e f f e c t i v e k e y w o r d ：c o n v e c t i o n d i f f u s i o ne q u a t i o n ，c o m p a c td i f f e r e n c es c h e m e ，d i f f e r e n c e o p e r a t o r , d i f f e r e n t i a lo p e r a t o r , e x t r a p o l a t i o nt e c h n i q u e n 2 4 四阶紧差分格式的误差分析1 3 2 5 常系数线性微分算子的六阶紧差分逼近1 4 2 6 六阶紧差分格式的直接解算子1 7 2 7 六阶紧差分格式的误差分析1 7 2 8 六阶数值外推1 7 2 9 数值实验1 9 第三章二维对流扩散方程的紧差分格式2 5 3 1 对流扩散方程的二阶差分格式2 6 3 2 对流扩散方程的四阶紧差分格式2 9 3 3 九点紧差分格式的矩阵形式3 2 3 4 九点紧差分格式的求解算法3 3 3 5 数值实验3 4 小结4 1 致谢4 2 参考文献4 3 i i i 第一章绪论 第一章绪论弟一早珀t 匕 程是偏微分方程中非常重要的一类基本运动方程，是用来描述黏 方程的线性化模型方程，它可以用来描述河流、大气、核废物等 的分布，以及流体的流动和流体中热的传导等众多物理现象。在 技术领域中的关于流体力学乃至生物体中的些物理现象的数值 都可以归结为对流扩散方程。因此，对对流扩散方程数值解的研 要的理论和实际应用价值的。 程的数值解的求解方法有多种，像有限体积法、有限差分法、有 种方法都有其独特的优势。目前，关于对流扩散方程的定解问题 有有限差分法和有限元法【l 。3 两大类。有限元方法的优势在于处理 灵活，求解的精度高，但其也有不足之处，主要表现在：往往需 要求解大型的带状稀疏矩阵，其计算量和存储量的需求也很大，实现隐格式比较 困难，在编制计算机程序时的工作量比较大。因此，这在一定程度上阻碍了有限 元方法的进一步发展和应用。 有限差分方法作为求解对流扩散方程数值解的一种非常传统的数值计算方 法，在经过了几十年的发展之后，已经取得了很大的成功，尤其是最近十几年来 发展迅速，研究成果颇多。有限差分方法的主要思想就是在微分方程中用离散的 函数值的线性组合来代替导数，从而得到相应的差分方程，通过解差分方程得到 微分方程解的近似值。而我们的主要目的就是通过不断的改进差分格式来降低误 差，使精度得到提高。 本文主要研究的是一维和二维常系数对流扩散方程的紧差分格式及其数值解 法。关于对流扩散方程的差分格式的建立是有多种的，其中隐格式稳定性好，但 是计算量大，需要好的求解算法。而显格式虽然稳定性较差，但比较容易得到， 也易于计算，所以在实际问题中也有很多采用该种格式。现在，对于对流扩散方 程的差分方法的数值计算研究主要集中在构造高精度的差分格式上。与以往文献 中构造差分格式的方法有所不同，本文主要从差分算子与微分算子的关系的角度 出发，最终分别建立三点、五点和九点的紧差分格式。这种方法的运用，避免了 在建立差分格式的过程中，对方程右端项的重复利用，使得推导过程更加简洁， 减少了很多的工作量。而在对流扩散方程的数值求解方面，我们主要通过差分格 电子科技大学硕士学位论文 式的特点，提取出计算模板，然后根据模板对一些经典的迭代算法进行改进，使 得求解过程的工作量有所减小，最终求出方程的数值解。论文的主要结构如下： 第一章介绍了对流扩散方程的研究背景及其研究方法。 第二章从算子的角度出发推导了一维对流扩散方程的四阶和六阶紧差分格 式，并对紧差分格式做了误差分析；然后，分析了三对角t o e p l i t z 矩阵的谱和谱 条件数；最后，对六阶数值外推技术做了一些讨论，并根据r i c h a r d s o n 外推公式 推导出了一个可以用来作为误差修正的公式。 第三章从算子的角度出发推导了二维对流扩散方程的二阶和四阶紧差分格 式，并分别给出了这两种差分格式的计算模板。最后，根据计算模板对方程进行 数值求解。 2 。在流体力学 乃至生物体中的一些物理现象也可以归结为扩散问题。因此，研究对流扩散方程 的数值解法有着非常重要的理论和实际应用价值。 关于对流扩散方程的差分格式的建立是有多种的。经典的有限差分方法是用 导数的有限差分逼近公式去取代微分方程中的导数，从而将微分方程化为只含一 组未知数的代数方程，然后解方程组计算出网格点上的函数值的数值逼近值。而 本章主要考虑的是从算子的角度出发研究一维常系数对流扩散方程的边值问题： r 却u “( a 嚣o lu 八( b 2 芦6 沼， 1) = ，) = 。 其中，f ( x ) 为已知函数，且有p 0 ，q 0 。 为了方便起见，本文中对求解区域进行网格剖分时，采用的是等距网格剖分。 把求解区域剖分成矩形网格点后，这些做剖分的直线称作网格线，网格线的交点 就称为网格点。 取正整数尼，令= 2 ，则记 = b - f a 为步长，可得到区间 口，6 上的网格点为： = a + 乃，( j = 0 , 1 ，加， 函数”( 石) 在网格点处的值记为u ( x ，) ，其数值逼近值记为“，。 记一阶中心差分算子和二阶中心差分算子分别为覆和正2 ，使得 8 ，u ( x j ) = u ( x j + 1 ) 一u ( x j 一】) 正2 u ( x a = “( _ + 1 ) - 2 u ( x 2 ) + u ( x ，一1 ) 记微分算子符号d ：；。定义位移算子e 的意义为： e u ( x j ) = “( x ，+ 1 ) ，e “( x ，) = “( 工，- 1 ) 电子科技大学硕士学位论文 故一阶和二阶中心差分算子司用位移算于e 表不为： 五= e e ，瓯2 = e - 2 i + e 。 其中，i 表示单位算子( 在不引起符号混乱的情况下，通常省略为1 ) 。利用t a y l o r 级数，可得： e = p _ | l d = 喜去c 五。，。，e = e 一柚= 喜去c h d ) 。 所以，将差分算子做泰勒级数展开，并取其中的有限项，可得到： 8 ，u ( x j ) = ( 2 h d + 庇3 d 3 ) “( 乃) + d ( 5 ) ( 2 2 ) 0 。 瓯2 “( _ ) = ( 砌2 + 壶砌4 ) “( _ ) + o ( h 6 ) ( 2 - 3 ) 最吣户( 2 h d + 2 6 h 3 d 3 + 去砌5 ) 吣抄o ( h 7 ) ( 2 - 4 ) 正2 “( x ) = ( 五2 d 2 + 壶办4 d 4 + 去五6 d 6 ) “( x ) + o ( h 8 ) ( 2 5 ) 2 1 常系数线性微分算子的四阶紧差分逼近 关于一维对流扩散方程四阶精度的有限差分格式的研究比较多。文献 1 0 1 运 用泰勒展开公式和数项级数的收敛性给出了一维线性对流扩散问题的一类高精度 三点紧致差分格式，具体给出了四阶精度的差分格式，如不考虑计算的复杂性， 利用该思想理论上可继续推导出更高精度的格式。文献 1 1 1 贝1 将一阶和二阶导数 的四阶差分逼近公式去替代微分方程中的导数，从而得到一种具有四阶精度的差 分格式。但该种格式涉及到相邻的五个点，并不紧凑，故相较于用相邻的三个点 建立起来的紧致差分格式来说，其效果要差很多。文献 1 2 1 从中心差分公式出发， 利用二阶微分的四阶差分公式，并结合方程本身，得到了数值求解两点边值问题 的一种紧致差分格式，且具有四阶精度。该格式仅涉及到中心点及相邻两个网格 点，相较以往的格式不但具有更高的精度，而且计算简便。 下面考虑将对流扩散方程( 2 1 ) 中的一阶和二阶导数项看作一个整体微分算 子，并从微分算子与差分算子的关系方面入手，建立微分方程( 2 1 ) 的四阶精度 的紧差分格式。 首先记微分算子符号： 4 第二章一维对流扩散方程的紧差分格式 l = - d 2 + p d ( 2 6 ) 胜用公瓦( 2 - 2 ) 和( 2 - 3 ) ，结苗瓦( 2 - 6 ) ，口j 得剑一个微分算于与爱分算于的夫 系式： ( 一古文2 + 轰戈m _ ) = ( 三一击j l z 2 d 4 + 詈五2 d 3 m _ ) + o ( h 4 ) ( 2 - 7 ) 为得到四阶精度的紧差分格式，需要将式( 2 7 ) 中右端项的高阶导数用微分算子 符号l 表示出来。由( 2 - 3 ) 式，可得： 去p 2 跏( 一) = 击嘞2 此( _ ) + o ( h 4 ) ( 2 - 8 ) 将式( 2 7 ) 和式( 2 8 ) 两端分别相加，可得： ( 一吉4 2 + 轰最一击p 2 暝2 ) “( 。) = 陋一击职一。2 + 印) 2 “( 勺) + o ( h 4 ) 整理上式，可得： - ( 古+ 苦) 正2 + 磊戈 “( _ ) = 犯一1 1 2 五2 r ) 甜( _ ) + d ( n ( 2 - 9 ) 去掉截断误差项o ( h 4 ) ，可得到四阶的紧差分算子逼近： 【_ ( 嘉+ 苦) 疋2 + 轰屯- m = l “，一西1j l z 2 r “， ( 2 - 1 0 ) 由微分方程( 2 1 ) 可得： l u j = f j q uj e u j = 巧j q hj = 巧j 一时j 七q 3 uj 代入式( 2 1 0 ) 整理，然后 记差分算子： 弘一( 矿1 + 昏疋2 + 轰最怕+ 1 z h ：q 2 ) 电子科技大学硕士学位论文 记微分算子的多项式为： 讹) = ( 1 + 击哟) 一击耽 故可得到微分方程( 2 1 ) 的四阶紧差分格式： l h u j = p ( c ) f j 将式( 2 1 1 ) 展开整理，得到四阶紧差分格式对应的三对角方程组： 其中： ( 2 1 1 ) 雄- 1 + 及跣j + 肛+ 1 = 奶+ 6 奶 ( 2 1 2 ) 州吾+ 争怕+ _ 1 h 2 q2 ， 由式( 2 1 2 ) 可以看出，四阶紧差分格式中仅涉及吩1 ，“，“j + 1 三个结点处 的未知函数值，故称为三点紧差分格式。 2 2 四阶紧差分格式的直接解算子 结合微分方程( 2 1 ) 中的边界条件，由三对角方程组( 2 1 2 ) 可得到四阶紧 差分格式( 2 1 1 ) 对应的线性方程组： 彳u = f 6 ( 2 1 3 ) p 一拍 p 一拍 + 一 一 矿一心 矿一挖 一屹 矿一挖 。一矿 。一旷 一 扣 一 一 口 = = 8 y 第二章一维对流扩散方程的紧差分格式 其中，a = t r i d i a g y ，口，f l 是- e n 阶三对角t o e p l i t z 矩阵，u - u 1 “2 “。 7 是一个，z 维列向量，f 是一维数组按行自然排序且结合边界条件所构成的列向 一 里。 对于实三对角t o e p l i t z 矩阵a ，可以通过对角相似变换转化为对称或t 反对 称”的形式，以得出线性方程组( 2 1 3 ) 的直接解算子。 首先取咒阶三对角矩阵 4 】： d = d i a g ( 1 ，万，万2 ，万”。1 ) ( 万 0 ) 对矩阵彳做对角相似变换，则有： d a d = 伐 8 6 万一1 7 口 。 p 6 万一1 ， 口 记彳= d a d 。下面分两种情况进行讨论： 1 ) 若励 0 ，n 令p d - d 1 7 ，故有万= 圻劢，且可得到： 2 - 3 1 仰= 乃谢昭p 劬( 万，s 劬( 历万) 显然，j 是一个对称的三对角t o e p l i t z 矩阵，故可得它的刀个特征值为4 】： 以= 舢历c o s 熹，2 ，z ( 2 设对应特征值五的特征向量为工。= ( 而，j c 2 n ，) ，则： 牡辱s 协鲁 川幺，咒 构造矩阵q = _ ，_ 劲，_ ” ，j = l ，2 ，z 。显然，q 是一个对称的正交矩阵， 即： 7 电子科技大学硕士学位论文 q = q r = q _ 2 ) 若励 o ，则令胪= 一万一1 y ，故有万= 丽，且可得到： 彳：d 一以d ：t r i d i a g s i g n ( y ) 、：2 - f l y ，口，s 切( ) 厢) 可得彳的n 个特征值为【4 】： 矗= 口捌厨c o s 熹， 设对应特征值友的特征向量为x 。- - ( x l ( k ) 恐，。) r ，则： 牡矿1 后s m 等 一幺，z 构造矩阵q = _ ，_ 2 1 ，乃”】，j = 1 ，2 ，咒。 则不论上述哪种情况，都有： = 。1 a mqada d d i a g ( i q 一，乃) o = “ ，五) e d 一= d i a g ( ；q ，五。) ，则由上式可得： a = d o _ 一1 d , 4 q d 一1 式带入线性方程组( 2 1 3 ) ，整理可得： d q d 4 q d u = f 阶紧差分格式的直接解算子可表示为： u = d q 一1 d a - q d 一1 f 第二章一维对流扩散方程的紧差分格式 2 3 紧差分格式的谱及谱条件数分析 在文献 4 中，作者讨论了三对角t o e p l i t z 矩阵的特征值和特征向量。本节将 在上述基础上，对三对角t o e p l i t z 矩阵的谱和谱条件数作进一步的分析。 2 3 1 三对角t o e p l i t z 矩阵的谱分析 由上一节对三对角t o e p l i t z 矩阵彳的分析，可知： 1 当，7 同号时，矩阵彳是一个对称的三对角t o e p l i t z 矩阵。 设矩阵彳为非奇异矩阵，则： 1 ) 当口0 时，由( 2 1 4 ) 式可得到矩阵彳的谱半径： p ( 面。峄陋) i = 口+ 2 历c o s 熹 同理，由( 2 - 1 4 ) 式可知，存在一个正整数毛k ，使得l 友l 取得最小值，不妨设 为r 1 。即： ”嘤n 件+ 2 历c o s 矧 则可得到矩阵彳一t 的谱半径： p ( 神= , m i n l 以( 彳) i = 1 , 7 。 2 ) 当口 0 时，由( 2 1 4 ) 式可得到矩阵彳的谱半径： p ( 两= 峄) i = - ( 口+ 2 历c o s i n + x 1 ) 同理，由( 2 1 4 ) 式可知，存在一个正整数尼：k ，使得l 以l 取得最小值，不妨设 为r 2 。即： 电子科技大学硕士学位论文 铲刊删= 卜历c o s 矧 则可得到矩阵彳。1 的谱半径： p ( 神= 1 哑n i 以( 彳) l = 1 刀： 2 当，7 异号时，矩阵彳是一个非对称的三对角t o e p l i t z 矩阵。 由式( 2 1 5 ) 可知，矩阵j 显然是一个非奇异矩阵。 故可得到矩阵彳和彳- 1 的谱半径分别为： 彪) = 峄l 以( 彳) l = 口2 _ 4 办耐熹 夕( 排州硎= 1 、 a , 2 _ 4 f l y c o s = ( n + 1 ) 2 由上节的分析可知，矩阵么相似于矩阵彳，故矩阵么与彳具有相同的谱半径。 即： ( 2 1 6 ) t o e p l i t z 矩阵时，其 第二章一维对流扩散方程的紧差分格式 当彳为非对称三对角t o e p l i t z 矩阵时，其谱条件数为： h a ) = a 2 - 4 f l y e o s 2 熹 4 励耐k 胛暑攀 证明 当矩阵彳为对称三对角t o e p l i t z 矩阵时，则有： 0 彳0 ：= p ( 彳) ，0 彳_ 19 ：= p ( 彳- 1 ) 故有：砥j ) = p ( 彳) p ( 3 1 ) ，结论得证。 当互为非对称三对角t o e p l i t z 矩阵时，不妨设： ( 2 1 7 ) 其中，为，2 阶单位矩阵，s = t r i d i a g 一1 ，0 ，1 ) 为n 阶反对称三对角矩阵。故有： 且s 的n 个特征值 5 1 为： 则有： s 胃= 一s 五，= 一2 i c o s 兰( ：1 ，玎) 。 ，z + 1 ”7 彳日彳= ( 甜+ 二万弦片) ( 甜+ 廊) = 口2 ，+ 口厮( s + s ) 一筘s = 饯2 i 一口冷hs = 铲i + 8 7 s 2 = i j h 故根据式( 2 1 5 ) 可得，矩阵彳彳的r 个特征值为： 7 = 及, 2 - 4 办c o s 2 鬲j x ( = 1 ，z ) 则可求出矩阵彳彳和矩阵( 彳月彳) 一1 的谱半径分别为： 电子科技大学硕士学位论文 p ( 彳胃彳) = 乎l 乃( 彳胃昝饼2 4 办曲s 2 三n + l = p ( 研 肌彳面_ l 】= 1 眄n l 乃( j 片彳) i = ( a 2 - 4 , 3 脚2 等掣班湫一) 2 所以矩阵彳和矩阵彳一1 的谱范数分别为： 故有： 0 j 0 ：= p ( 彳) ，l i 彳- 10 ：= 夕( 彳- 1 ) 砥彳) = ：p i | 2 = 户( 彳) p ( 彳一) = a 2 - 4 f l y c o s 2 熹 扛4 办c o s 2k 行暑斗 ( 2 1 8 ) 结论得证。 由于矩阵的谱范数是酉不变范数，故它是相容的矩阵范数，结合式( 2 - 1 7 ) 和式( 2 1 8 ) 可得： i i a l l ：= 0 砌。1 l l ：l i d i i ：p l l 2 = p ( d ) 户( 彳) p ( d 。) ( 2 - 1 9 ) 同理可得： p l l 2 = i l 旃d 1 0 ：- = l l m l l ：0 么一1l i ： p ( d ) 2 p ( d 一1 ) 2p ( j ) ，p ( j 1 ) = p ( d ) 2 p ( d 一1 ) 2 砥彳) ( 2 2 1 ) 当万1 时，即i 州i 剧，有： 胭川= 肝h 协2 2 ， 当0 万 o ，使得：l i a - f l ：c ：a 所以 0 e l l ：q c 2 h 4 由范数等价性可知，存在常数c 。 0 ，使得： 故有： 其中，c = c l c 2 c 3 ，定理得证。 乎l “( 乃) - u j c 忍4 2 5 常系数线性微分算子的六阶紧差分逼近 关于一维对流扩散方程六阶精度的有限差分格式的研究并不是很多。在文献 1 3 和 1 4 中，作者通过埃尔米特插值多项式以及泰勒公式，推导出了具有六阶精 度的三点紧致差分格式。文献 1 5 1 贝u 在四阶紧差分格式的基础上，通过外推技术 得到了一种具有六阶精度的差分格式。在2 1 节中推导的三点紧差分格式虽然在 精度上有所提高，但是其误差精度仅是四阶的( o ( h 4 ) ) 。本节将根据建立四阶紧 差分格式的思路推导六阶精度的紧差分格式，不但较四阶的紧差分格式在精度上 有了进一步的提高，而且也具有良好的稳定性。 首先由微分算子( 2 6 ) l = 一d 2 + 加 可以得到： l 2 = ( 一d 2 + p d ) 2 = d 4 2 p d 3 + p 2 d 2 l 3 = ( - d 2 + p d ) 3 = 一d 6 + 3 p d 5 3 p 2 d 4 + p 3 d 3 由( 2 4 ) 式和( 2 5 ) 式可得到一个微分算子与差分算子的关系式： 1 4 ( 2 2 4 ) ( 2 2 5 ) 子符号三表示出来。由( 2 3 ) 式可知： 瓯2 u ( x ，) = h ：d 2 u ( x ，) + d ( j i z 4 ) 疋2 “( _ ) = ( h ：d ：h l i 4 。4 ) “( x ，) + 。( 五6 ) 为得到具有o ( h 6 ) 的格式，可将上面两式化简为： 譬咖户苦珊户等嘶抄 h 4 p 4 d 2 u ( x ，) = h 2 p 4 疋2 u ( x ) + d ( 6 ) 将上面两式代入( 2 2 7 ) 式化简整理，可得： 嗡h 2 p 4 + 轰珊( 乃) - 州等h 协2 ) l 2 + 嘉卸( + d ( 哟 ( 2 2 8 ) 去掉截断误差项o ( h 6 ) ，可得到六阶的紧差分算子逼近： ( 面h 2 p 4 音 ) 瓯2 + 轰珏，= + ( 面h 4 p 2 2 u j + 嘉5 0 l 3 “， 电子科技大学石贞士学位论文 由微分方程( 2 1 ) 可知： l uj = j q u 3 i , 2 uj = l j q l uj = u j q j j 七t u j u i = e j q l jj 七t l uj = t j q l j3 七t j q 3 u j 带入上式整理，然后 记差分算子： 厶= ( 等一古一鲁) 正2 + 轰瓦+ q _ q 2 ( 面h 4 p ：一西h 2 ) + g 。丽h 4 记微分算子的多项式为： m ) = 1 - g ( 面h 4p 2 一刮h 2 十9 2 爿h 4 + ( 面h 4 p 2 h i 2 ) 一g 翥工+ 羔r 故可得到微分方程( 2 1 ) 的六阶紧差分格式： l h “，= e ( l ) f j ( 2 - 2 9 ) 将式( 2 2 9 ) 展开整理，得到六阶紧差分格式对应的三对角方程组： 其中： 弘j 、+ 僦j + 8 u j n = j 礼u j + c ef j 0 2 - 3 0 ) 口= 一导一吾一譬) + q _ q 2 ( 面h 4 p 2 一西h 2 ，十9 3 丽h 4 尸= 一吉+ 轰 + 等 1 pp 2 h 2 p 4 y = h 2 2 h 一1 2 + 面 口= 1 一g ( 面h 4 p 2 一h 1 2 2 q 2 j h 。4 石 1 6 第二章一维对流扩散方程的紧差分格式 6 ：一生+ 盟一口旦 1 2 7 2 0 1 3 6 0 h 4 3 6 0 2 6 六阶紧差分格式的直接解算子 结合方程( 2 1 ) 式中的边界条件，由三对角方程组( 2 3 0 ) 可得到六阶紧差 分格式( 2 2 9 ) 对应的线性方程组： 么u = f ( 2 3 1 ) 其中，a = t r i d i a g ) i ，口，历是一个，2 阶三对角t o e p l i t z 矩阵，u = “lu 2 ”。】r 是一个玎维列向量，f 是一维数组f 按行自然排序且结合边界条件所构成的列向 量。 同四阶紧差分格式的直接解算子推导过程一样，我们可以得出方程( 2 3 1 ) 的直接解算子为： u = d q 一1 d a 一1 q d 一1 f 2 7 六阶紧差分格式的误差分析 定理2 3 设线性方程组( 2 3 1 ) 对应的系数矩阵么为非奇异矩阵，则对差分 格式( 2 - 2 9 ) 的解“，存在常数c 0 ，有如下的误差估计： 平l “( _ ) - u j c h 6 证明过程与定理2 2 的证明过程类似，略。 2 8 六阶数值外推 外推法实际上是一种古老的加速收敛方法，其本质是依赖于近似误差对单个 电子科技大学硕士学位论文 离散参数的渐进展开。它已成为数值分析各个领域中具有普遍意义的加速收敛技 术。文献 1 5 ，1 7 ，1 8 ，1 9 就是在四阶差分格式的基t i t _ l ，通过外推技术得t i 六阶精度 的有限差分格式。本节将从数值方面来讨论外推技术。 已知r i c h a r d s o n 外推公式【6 】可表示为： 鹂= 竿 其中，p 是外推之前的精度，外推之后的精度应该增加到p + 2 。“，2 6 和“易分别是 外推之前粗网格q ：。和细网格q 。上偶数点的计算结果。磁是外推之后在细网格 q 。上的偶数点的计算结果。 1 禾i i ii i i , 4 专淼一林群推广公式作误差修正 由r i c h a r d s o n 外推公式可得： 踢叫，= 等竿呓= 错 将上式化简可得： 磅= “乞+ 锊( 2 - 3 2 ) 式( 2 3 2 ) 表明：外推之后细网格q 。上偶数点的计算结果即为外推之前该点 的值加上一个误差修正值，该误差修正值即为外推前细网格q 。上的偶数点的值减 去相应的粗网格q ：。上的点的值再乘以一个倍数上2 - 1 。 由于在做数值外推时，细网格q 。上的奇数点在粗网格q ：。上没有相应的对应 点，则根据上述思想，细网格上奇数点的误差修正值可i y , 其相邻两偶数点的误差 修正值的平均值，即： u 2 h _ l - - u 乞一一+ 互_ 南 ( “：h j 一：一“2 一h t ) + ( “：h ，一“j 2 6 ) 】( 2 - 3 3 ) 第二章一维对流扩散方程的紧差分格式 我们把式( 2 3 3 ) 称为陈传淼一林群推广公式。 当外推之前的精度p = 2 时，便可得到陈一林公式7 1 ： y ( 弓) = y ；( 弓) + 吉 y 尝( x ，) 一y 一( 一) + y ：( z ，卅) 一y 一( z ，卅) 2 基于四阶紧差分格式的六阶数值外推算法 1 ) 取定一个值，利用四阶紧差分格式分别求出粗网格q ：。和细网格q 。上 的方程的近似解u 。和。 2 ) 对u 和u 2 ，利用r i c h a r d s o n 外推公式，计算出细网格q 。偶数点上的外 推结果。 3 ) 对于细网格q 上奇数点上的外推结果，可利用陈一林推广公式( 2 - 3 3 ) 求得结果，此时取p = 4 ，即： 磅一，= h + 爿1 ( “乞- 2 - “川2 h ) + ( “务一“2 6 ) 2 9 数值实验 前面给出了一维对流扩散方程的四阶和六阶两种紧差分格式，并从理论上分 析了它们的误差。下面将通过数值实验来进一步验证上面结论的正确性。考虑对 流扩散型的微分方程 一u 目+ b u 。+ “= b c o s x + 2 s i n x 0 x 刀 具有边界条件：u ( o ) = 0 ，u ( j r ) = 0 。这一问题的解析解为：u ( x ) = s i n x 。 分别取b = l 和b = 1 0 0 时计算，利用四阶和六阶紧差分格式求数值解，并比较 结果，差分解与准确解的最大绝对误差见表2 1 。 从表2 1 中可以看出两种格式都具有良好的稳定性和较高的误差精度，是求 解对流扩散方程的一种行之有效的方法。 1 9 电子科技大学硕士学位论文 表2 1三点紧差分格式实验数据表 b = lb = 1 0 0 步长h误差误差 四阶精度六阶精度 四阶精度六阶精度 1 0 4o 0 0 1 03 3 4 17 e 0 0 50 0 5 2 50 0 5 3 6 86 2 0 4 0 e 0 0 55 3 1 6 1 e 0 0 7o 0 1 1 7o 0 1 3 5 1 63 8 5 2 8 e 0 0 68 4 5 71e 0 0 90 0 0 1 60 0 0 3 7 3 22 4 0 4 0 e 0 0 71 3 2 0 8 e 0 101 2 2 6 4 e 一0 0 42 3 5 3 4 e 一0 0 4 矧6 41 5 0 2 5 e 0 0 82 0 8 6 6 e 0127 9 0 31e 0 0 63 3 6 1 1 e 0 0 6 1 2 89 3 9 5 3 e 0 104 9 7 3 8 e 0 1 44 9 5 3 6 e 0 0 7 4 9 7 6 8 e 0 0 8 2 5 65 8 1 6 9 e 0 1 18 6 3 3 1 e 0 1 33 0 9 6 8 e 0 0 87 7 7 6 6 e 0 1 0 5 1 2 1 9 0 4 1 e 0 1 2 3 9 9 6 8 e 0 13 1 9 35 4 e 0 0 91 2 1 0 6 e 0 1 1 1 0 2 4 4 7 5 7 3 e 012 7 6 2 0 5 e 012 1 2 0 5 4 e 0 106 10 6 2 e 0 1 4 由实验数据可知，误差与h = 有关。记误差数据为e ( h ，) ，( j = 1 , 2 ，z ) 。 j 设e ( h ) zc h p 。为了利用实验数据来求得数值阶数，将等式两端取对数，得： l n e ( h ) = l n c + p l n h 其中p 是差分逼近阶数，由实验数据拟合确定的p 就是数值阶数。在对数尺度下， 尸是线性函数的斜率。步长的对数尺度数据是单调增数列： 一l n h l ，- l n h 2 ，一l n h 。 该数列的差为：l n h ，- l n h j 一= l n ( h ，h i - 1 ) 。对应的误差对数尺度数列为： 一1 h e ( h , ) ，一l n e ( h 2 ) ，一l n e ( h 。) 用两组数据做线性拟合，一次项系数即为差分逼近的数值阶数尸，它是理论阶数 的近似值。 图2 1误差对数尺度线性拟合图图2 - 2 误差对数尺度线性拟合图 图2 1 和图2 2 分别是当b = l 和b = 1 0 0 时针对表2 1 中四阶紧差分格式的误 差数据所做的误差对数尺度的线性拟合，直线斜率分别为：p = 4 0 0 4 1 ，p = 3 9 4 6 7 ， 分别说明的是四阶紧差分格式解的数值阶数。 图2 3误差对数尺度线性拟合图图2 - 4 误差对数尺度线性拟合图 图2 3 和图2 4 则说明的是六阶紧差分格式的误差对数尺度的线性拟合，其 当b = l 和b = 1 0 0 时六阶紧差分格式的数值阶数分别为p = 5 8 9 9 9 ，p = 5 9 9 8 5 。 实验结果表明本文中的四阶和六阶紧差分格式的数值解的数值阶数与差分理 论解的阶数是一致的，从而验证了格式的有效性。 2 l 电子科技大学硕士学位论文 图2 5 数值解与精确解的曲线图图2 - 6 数值解与精确解的曲线图 图2 5 和图2 - 6 所示的分别是b = l 时四阶紧差分格式取1 6 个点和1 2 8 个点时 的数值解和精确解的曲线图。 图2 7 数值解与精确解的曲线图 图2 - 8 数值解与精确解的曲线图 图2 7 和图2 8 所示的则分别是b = l 时六阶紧差分格式取1 6 个点和1 2 8 个点 时的数值解和精确解的曲线图。 图2 5 到图2 8 中圆圈表示数值计算时所取的点，实线表示精确解。从图中 可以看出，随着网格点的增加，数值解与精确解的曲线图也越来越吻合，说明了 数值解与精确解之间的误差越来越小，精确度也越来越高。 第二章维对流扩散方程的紧差分格式 表2 2 矩阵彳和么的谱条件数实验 b = l b - - 1 0 0 步长h k ( 么) 硬彳) k ( 么) h a ) 霭43 2 7 1 03 0 0 9 35 8 2 8 05 8 2 8 0 81 2 4 8 7 8 1 1 6 1 5 52 5 2 7 6 02 5 2 7 3 3 1 6 4 9 4 2 0 3 4 6 1 8 9 89 8 8 2 8 49 5 6 9 9 7 3 21 9 7 1 7 3 8 18 4 5 2 4 5 3 2 3 6 2 01 2 6 6 6 1 0 6 47 8 8 1 9 4 77 3 7 8 7 2 6 7 1 7 8 4 9 1 1 1 4 4 霭1 2 83 1 5 2 3 e + 0 0 32 9 513 e + 0 0 3l8 9 2 0 4 13 2 9 0 5 积2 5 61 2 6 0 9 e + 0 0 41 1 8 0 5 e + 0 0 4 5 8 5 0 8 8 51 1 2 8 5 5 5 1 25 0 4 3 4 e + 0 0 44 7 2 1 9 e + 0 0 42 14 9 8 e + 0 0 34 3 1 3 0 8 刀1 0 2 42 0 17 4 e + 0 0 51 8 8 8 8 e + 0 0 58 4 0 4 1e + 0 0 31 7 0 5 1 9 4