（计算数学专业论文）padé型小模板紧致格式.pdf

摘要本文主要介绍了各类高精度紧致格式及其性质。紧致格式最大的优点在于模板比较小从而计算起来比较简单高效，全局计算的性质又使之可以避免振荡对于边界条件的处理，紧致格式也具有一定的优越性，体现在不依赖于边界条件的类型的约束上，以及从整格点计算到半格点的高精度格式，避免了平均算子的耗散性紧致格式包括很多内容，其中比较常见的有三种第一种是为p a d 6形式的有限差分紧致格式，这里我们特别使用了一点小的技巧使高精度与小模板这样一对数值计算的矛盾得到部允的化解，数值实验证明结果令人满意；利用它来处理高阶非线性发展方程的办法相对其他差分格式也比较简单。第二种是有限体积的紧致格式，这种方法由于融入了有限体积格式的守恒性，间断解等优点，在对浅水问题的模拟上表现出比较好的性质，可惜在当前为了使两种方法得到融合，一定的代价被付出了，就是原来有限体积法对网格的适应性被破坏了，在这里只能对正交网格进行处理；第三种是本质无振荡的紧致格式，在处理间断问题的时候，如s o d 激波管问题，具有比普通e n o 、w e n o 更高的效率也保证了解在间断处的锐利，是一种具有高分辨率的特性的数值模拟方法；我们还简要介绍了发展较新的群速度控制紧致格式，这种方法针对从物理学角度看间断的产生原因，提出了群速度控制抑制振荡的解决办法，结合紧致格式的优点，与本质无振荡格式具有同样的优点。关键词：紧致格式，有限差分方法，有限体积方法，本质无振荡格式群速度控制格式，孤立子，溃坝问题，激波，p a d 6 逼近1 1 1a b s t r a c ti nt h i st h e s i sw em a i n l yf o c u so ni n t r o d u c i n gs e v e r a lk i n d so fc o m p a c ts c h e m e sa n dt h e i rc h a r a c t e r s a sw ek n o w ，t h eb e s ta d v a n t a g eo fc o m p a c ts c h e m e si ss m a l ls t e n c i la n dt h ee a s ym e t h o dt oa c t u a l i z e a n db e c a u s et h es t e n c i l i sb a s e do nt h et o t a ln o d e s t h eo s c i l l a t i o nc a nb ea v o i d e d w h e nw ed e a lw i t ht h eb o u n d a r yc o n d i t i o n ，w ec a na l s of i n dt h ea d v a n t a g eo fc o m p a c ts c h e m et h a tt h es c h e m ei si n d e p e n d e n to ft y p e so fb o u n d a r yc o n d i t i o na n dd i s s i p a t i o na v o i d a n c e t h e r ea r es om a n yc o m p a c ts c h e m e s ；h o w e v e rw ew i l lj u s ti n t r o d u c et h r e ek i n d so ft h o s e t h ef i r s to n ei sc o m p a c tf i n i t ed i f i e r e n c es c h e m ei np a d ef o r m h e r ew eu s es o m es p e c i a ls k i l l st os o l v et h ec o n f l i c tb e t w e e ns t e n c i la n dp r e c i s i o na n dt h en u m e r i c a lr e s u l t ss h o wt h a tt h em e t h o di ss a t i s f i e d a sac o n s e q u e n c e ，t h el a s tn u m e r i c a ls c h e m ei sv e r ye a s y t h es e c o n do n ei sc o m p a c tf i n i t ev o l u m es c h e m ew h i c ha l s oh a sa d v a n t a g e so ff i n i t ev o l u m en l e t h o d ，s u c ha sc o n s e r v a t i o na n dd i s c o n t i n u o u ss o l u t i o ns ot h a tw h e nw eu s et h i ss c h e m et os o l v es h a l l o ww a t e rp r o b l e m ，w ec a nf i n dt h a tt h es i m u l a t i o ni sv e r yg o o d h o w e v e rw eh a v et 0l o s et h ea d v a n t a g eo fg e n e r a lf v mt h a ti tc a ns o l v en o n u n i f o r mg r i d t h et h i r do n ei sc o m p a c te s s e n t i a ln o n o s c i l l a t i o ns c h e m e w h e nw es o l v ed i s c o n t i n u o u sp r o b l e m s u c ha ss h o c kw a v ep r o b l e m t h i ss c h e m eh a sh i g h e re f f i c i e n c yt h a ne n oa n dw e n oa n ds t i l lh a sh i g hr e s o l u t i o ni nt h ed i s c o n t i n u o u sp o i n t sa se n oa n dw e n o w ea l s oi n t r o d u c eac o m p a c ts c h e m ew i t hg r o u pv e l o c i t yc o n t r o lw h i c hh a st h es a m ec h a r a c t e rw i t hc o m p a c te n o i ti sb r o u g h tf o r w a r df r o mt h ep h y s i c a lp o i n to fv i e w b e c a u s et h en u m e r i c a lo s c i l l a t i o ni sp r o d u c e db yt h ed i f f e r e n tg r o u pv e l o c i t i e s lh e r ew ec a nc o n t r o lt h ea l lg r o u pv e l o c i t ya n dc o m b i n et h ec o m p a c ts c h e m et oa v o i dt h eo s c i l l a t i o n k e y w o r d s ：c o m p a c ts c h e m e s ，f i n i t ed i f f e r e n c em e t h o d ，f i n i t ev o l u m e nm e t h o d ，e s s e n t i a l l yn o n - o s c i l l a t i o n ，g r o u pv e l o c i t yc o n t r o ls c h e m e s ，s o l i t o n ，d a r n b r e a k i n gp r o b l e m ，s h o c kw a v e ，p a d 6a p p r o a c h 致谢本文是在刘儒勋教授的悉心指导下完成的，在此向他致以最诚挚的谢意。研究生三年，刘老师不仅在学习上不厌其烦地指导、帮助和督促，在生活上刘老师与师母也给予我们这群学生无微不至的关心。刘老师的研究生课程与讨论班都使我受益匪浅，他对学术的严谨和执着更是我应该效仿的。在数学系的这七年学习是我一生最宝贵的财富，感谢老师们精彩的讲课，使我一步步走近数学。特别是冯玉瑜老师，陈祖墀老师，舒其望老师，陈发来老师的关于计算数学的课给我打好了基础。这里也要感谢班主任黄稚新老师多年的爱护我还要感谢师姐融华、徐岩、师兄张瑞、张磊，师弟夏银华和实验室的各位同学在我学习和助教工作上的帮助，感谢p b 9 7 0 1 和s a o i o i 所有同窗，感谢他们给了我在科大最美好的回忆最后感谢我的父母一直以来对我的默默支持、呵护以及为我的学业所付出的艰辛第一章引论紧致格式是在1 9 世纪3 0 年代提出来的，经过几十年的发展，融合了一些其他的数值方法，它已经发展成为一个方法类本章简单地介绍了紧致方法的发展简史以及它的一些特性。1 1c o m p a c t 方法简史最早的紧致方法在刚提出的时候是为了解决湍流和声学模拟问题，它保证了短波的类谱性质，并且不受网格拓扑，边界条件类型的约束，具有一些优点，但是也有很多缺点。比如，和所有的中心格式一样，它容易引起振荡，稳定性较差，而且振幅并不随着网格尺度变小而减少，随后发展成g i b b s 现象随后有一些工作对紧致方法的各方面性质，稳定性，一致性等等进行了分析和数值模拟，9 0 年代初l e l e 3 2 1总结了前人的工作，在j c p 上发表了一篇关于紧致格式的文章，从这个时候开始，紧致格式变得非常受欢迎并陆续被引入不同的数值计算领域，融合了很多其他的一些数值方法的优点。众所周知，模板大不仅会引起数值振荡，系统矩阵处理起来也很麻烦，尤其在处理有多阶导数的发展方程的时候。p a x i 6 形式的紧致格式是一个可以考虑的解决办法，但是在处理很多不同阶的导数时候，它的优点也被抵消了很多。所以我们在这篇文章中采用了一点小的技巧，巧妙的选择了紧致格式的系数，将它化为一种分母相同的有理逼近的紧致格式。不仅保持了原来紧致格式的小模板、易于处理，不产生振荡的优点，而且仍然能达到空间四阶精度，对于一般的问题，显然只要控制时间的精度这就已经足够了。1 8 4 4 年j s r u s s e l l 把这种始终保持在水面上，向前平移的孤立水峰，叫做孤立波( s o l i t a r yw a t e rwa _ v e s ) 1 8 9 5 年k o r t e w e g 和他的学生d ev r i e s 提出了流体中单向波传播的数学模型：毗+ 6 u u z + 札z z = 0即通常所说的k o r t e w e g d ev r i e s 方程，简称k d v 方程本文利用对孤立子的模拟证明了紧致格式的一些计算上的优点有限体积法是上个世纪八十年代以来发展起来的一种新型的微分方程的离散方法长期以来，微分方程数值求解的离散方法中，使用得最多的是有限差分法和有限元方法。有限差分法着眼于求解区域剖2 0 0 4 年中国科学技术大学硕士毕业论文第2 页至三兰! ! 兰! ! ；! 呈! 竺! 兰! 童兰重兰分的节点上的函数值，方法简便、灵活，离散的格式丰富多样，在收敛性、稳定性等理论研究方面也比较完善。但由于计算中对求值节点的分布要求比较规则，不能适应复杂的几何求解域另外，它不能求许多实际物理问题中出现的弱解有限元方法基于微分方程的弱解形式和广义变分原理，利用能适应复杂的几何形状的求解区域网格剖分，在剖分单元上用形函数插值逼近来求解有限元方法能适应复杂的几何求解域，但有时需要求解大型的线性方程组，计算上没有差分方法那样灵活方便，并且处理大变形间断问题较困难。f v m 在一定程度上吸收了f d m 与f e m 的长处，同时又克服了它们的缺点f v m 从控制体的积分形式出发，对求解区域的剖分同f e m 一样具有单元特征，能适应复杂的求解区域，离散方法具有差分方法的灵活性，具有间断解的适应性。f v m 与紧致格式的结合虽然有一些优点，但是在处理不同阶导数的时候不可避免要解两个以上的大型矩阵。m a r c e oh k o b a y a s h i 2 7 1 在1 9 9 9 年提出了耦合方法很好的把多阶导数一次处理，方法的精度，稳定性，谱性质都利用f o u r i e r 分析证明。因为气动力学方程的解随时间推移会产生间断，即使初始条件是光滑的，所以对激波的捕捉成为计算流体力学的一个重要研究内容。如何改善激波附近的数值解成为首要的问题，有t v d 格式，e n o 格式，n n d 格式等。其中特另h 是1 9 8 6 年以来发展起来的e n o 方法，它利用巧妙的选择模板的算法使得捕捉的激波很锐利，但是大量的选择比较占用了较多的计算的时间，1 9 9 6 年a d a m s i i 较早的提出将紧致格式与e n o 格式耦合起来解决激波模拟问题，在精度上不逊于普通e n o 格式，面紧致格式也不会产生其他的数值振荡，同时时间上大大节省，提高了效率。正如我们所知，w e n o 格式在同样的模板上可以得到比e n o 格式更高精度的插值，这个思想与紧致格式利用小模板得到更高精度的插值有相似之处，所以这两种方法的结合也是很顺理成章的。近几年来有一些这方面的工作，结果都很不错虽然这些激波捕捉方法都有较高的精度和激波分辨率，但是仍然是传统型的差分格式，数值解中总是存在非物理振荡傅德薰和马延文在1 9 9 6 讨论了关于数值解的非物理振荡的产生原因，归结为非均一的波群的群速度 1 1 ) ，随后傅德薰和马延文提出了四阶群速度直接控制的迎风紧致格式 1 0 】，并证明该格式具有良好的激波捕捉能力他们又利用这类群速度控制的紧致格式格式计算了双曲守恒律方程，不可压扰动方程，n a v i e r - s t o k e s 方程，并作了相当详细的误差分析2 0 0 3 年他们又在六阶精度的紧致格式和五阶精度的迎风紧致格式的基础上，引入了一个控2 0 0 4 芷第一章引论中国科学技术大学硕士毕业论文第3 页1 2c o m p a c t 方法特性制数值解的群速度的函数，推导出六阶的群速度控制紧致格式【14 朱庆勇等2 0 0 4 年有一篇文章 4 1 关于迎风群速度控制的紧致格式对三维可压流的模拟，结果非常令人满意，这也标志紧致格式在越来越多，越来越复杂的流体力学问题中得到应用1 2c o m p a c t 方法特性紧致方法的特性主要有以下几点：( 1 ) p a d s 方法是联系全域总体而不是局部的从实质上说，所有求解p d e 的数值方法都是由已知的离散节点上的值或者区间的积分平均值来重构函数各阶导数值和半格点值如有限差分方法就是利用差商来求差分值，有限元方法是利用单元上的播值逼近，它们本质上走的是统一的路子。p a d s 方法却不是这样，它是利用全域上的节点值构造出全域上的导数值和半格点值。这样做的好处是避免了差商与多项式逼近，正如我们所知，高阶多项式必然引进振荡机制，比如g i b b s 现象，这是精度与振荡的不可调和的矛盾。( 2 ) p a d s 方法建立了整格点与半格点的逼近关系，一次生成全域的半格点值相比较。一般的算数平均，引入了很强的数值耗散，有光滑效果，不能达到高分辨率，p a d s 方法则不存在上述矛盾。( 3 ) 可以利用偏心的边界处理得到带状矩阵，带状矩阵的计算已经拥有较为成熟的理论系统，避免了大稀疏矩阵的难于处理，整个系统矩阵相对易于计算，计算效率很高( 4 ) 通过一定的技巧，可以统一处理各阶导数，真正达到小模板的优势，这部分我们做了一些有意义的工作。对于具有不同阶导数的发展方程，我们的小模板p a d s 格式更有优势。第二章p a d 6 型有限差分方法和高阶发展方程高阶非线性发展方程具有广泛而深刻的物理背景和现实意义。发展方程数值解的困难在于含有高阶导数而且人们很看重孤立子守恒数值研究，历史上有很多有名的数值数学家从事各种高阶、守恒、保形格式的研究尽管有很多成功的算例【3 6 】，【1 9 ，但是不尽满意的工作也是有的。问题在于高精度与数值振荡的矛盾使得精度难以提高，也在于模板过宽和实施的难度这一章在前人有理逼近和紧致格式研究和应用的基础上提出了一种简单有效的p a d 6 格式对于方程中不同的各阶导数采用同样分母的p a d 6 逼近形式，使节点模板达到尽可能的压缩，而且又引入了型格式求三对角阵的解法。使精度既保持了四阶的要求，而且可以实现隐式求解效果，对稳定性、代数方程组的收敛保证都有很好的效果2 1 紧致格式的引人在偏微分方程的数值计算中，有限差分方法一直以来都是最基本的重要方法在编程上它具有相对较为容易实现的优点，但是在模板与精度的关系中，它却不是很具有优越性。可以证明：为了线性逼近一个函数的m 阶导数，( ，并且要求达到阶精度，所需要的该函数的模板的节点数，至少应当等于最大导数阶与要求的精度阶之和，即m + 。但是如果采用有理逼近或者p a d 6 逼近，在同样的要求下，所要求的节点模板则可以大为节省。例如在线性差分逼近情况下，一般地( 厶) = n ，f 一2 + 占 一1 + c ， + 蛳+ 1 + e + 2 + o ( a x 4 )其中，”t = “( ) 是未知函数在空间网格z = 龟上的函数近似值，差分嗣格空间步长h = z 利用t a y l o r 展开，比较系数，我们可以得到：a + b + c + d + e = 0 一。a a + a + z e = 忐4 n + b + d + 4 e = 0 ，一8 a b + d + 8 e = 0 42 0 0 4 年中国科学技术大学硕士毕业论文第5 页苎三圭! ! 竺兰童! ! 兰垒童童皇童竺垒曼奎兰! 兰! 竺兰兰垒竺! 坠从上面的式子可以看出为了达到四阶精度，我们至少需要一个5 节点模板如此等等，不必赘述。当方程中具有更高阶的导数的时候，需要的模板更大，在边界处采用延拓点的处理使编程赋值很不方便，如果采用内插值的方法，想要达到很高的精度也有困难。在这里，缩小模板成为必须的选择l e l e 3 2 ，k o b a y a s h i 2 7 等给出了一般情况下的c o m p a c t 格式形式，p ( 一一2 + 矗2 ) + n ( 丘l + 十1 ) + 鼻= e 互十6 学+ 。学，( 2 1 1 )它的系数可以由t a y l o r 展开确定这里以一阶导数的系数为例说明：下口lo + 2 2 6 + 3 2 c = 2 鲁( n + 2 2 卢) ，( f o u r t h o r d e r )。+ 2 4 6 + 3 4 c = 2 象( 。+ 2 4 卢) ，( s i x t h - - o r d e r )n + 2 6 6 + 3 6 c = 2 磊陋+ 2 6 既( 扰h t 一。r d e r )n + 2 8 b + 3 8 c2 等( a + 2 8 既( 劬t 7 i - - o f 撕)( 2 1 2 )( 2 1 ，3 )f 2 i 4 1( 2 1 j 5 )( 2 1 ，6 )在四阶情况下，系数为卢=，。=；(n+2)，=；(40r-0buc = 0( 2 1 7 )卢= ，o = ；( q + 2 ) ，= i1 ) ，c =( 2 1 - 7 j在六阶情况下，系数为n = ；肛咖= 警，。= ；，c = o( 2 埘)当取方程左端为五节点模板，即卢0 的时候，四阶精度系数如。= ；( 4 + 2 。一1 6 卢+ 5 c ) ，6 = ；( 一1 + 4 + 2 2 3 8 c )( 2 1 9 )六阶精度系数如下：n = ；( 9 + a - 2 蚴，6 = 去( 一9 + 3 2 c r + 6 2 , 3 ) ，c = 击( 1 3 。+ 1 2 卢) ( 2 1 l 。)2 0 0 4 年中国科学技术大学硕士毕业论文第6 页董三兰兰竺! 竺童兰兰垒奎童塞童竺垒苎奎量墼；! 篁兰竺皇竺! ! 垒八阶精度系数如下：卢= k o ( 一3 + 8 。) ，。= ；( 1 2 7 。) ，6 = 1 击( 5 6 s 。一1 8 3 ) ，c = 丽1 ( 9 q 一4 1f 2 1 1 1 1十阶精度系数如下：a = ；，卢= 嘉，。= 西1 7 ，6 = 而1 0 1 ，c = 而1( 2 1 r 1 2 )下面给出系数与截断误差的表格式左模板右模板截断误差( 2 1 7 )35备( 3 a 1 ) 4 ，( 5 )( 2 1 - 8 )35可4 6 ，( 7 )( 2 1 9 )57帚( 一1 + 3 a 一1 2 # + l o c ) h 4 ，( 5 )( 2 1 1 0 )57鲁( 3 8 0 + 2 0 # ) h 6 ，( 7 )( 2 11 1 )57等( 2 0 1 ) 8 ，( 9 ( 2 1 1 2 )57揣 1 0 ，( 1 1 )l e l e 关于二阶导数的关系式如下：( 一：。十一j 。) + a ( 一：。+ 一j 。) + 7 = c 盎芏二竽+ 6 血挚+ 。血学，( 2 。)导数的关系a 十b4 - c = 1 + 2 q + 2 卢，( s e c o n d o r d e r )n + 2 2 6 + 3 2 c = 备( 。+ 2 2 鼽( 细m 一”圳n + 2 4 6 + 3 4 c = 筹( q 十2 4 鼽( s 诚 - o r d e r )。十2 6 6 十3 6 c = 等( n + 2 6 既( e 溉忱一。r 训。+ 2 8 6 + 3 8 c = 等( 。+ 2 8 既( t e n t h - o r d e r )在左端为三节点模板的情况下，四阶精度系数为( 2 1 1 4 )( 2 1 1 5 )( 2 1 1 6 )( 2 1 ，1 7 )( 2 1 1 8 )( 2 1 1 9 )卢= o ，n = ；( 1 一n ) ，6 = ；( 1 0 a - 1 ) ，c = o( 2 1 舶)2 0 0 4 年中国科学技术大学硕士毕业论文第7 页第二章p a d 6 型有限差分方法和高阶发展方程21 紧致格式的引入六阶精度系数为a = ；l , f l = o ，n = 五1 2 ，6 = 暑，c = o( 2 1 2 1 )当取方程左端为五节点模板，即卢0 的时候，四阶精度的系数如下：n = ；( 4 - 4 n 一4 0 f l + 5 c ) ，6 = ；( 一1 + 1 0 。+ 4 6 f l 一8 c )( 2 1 2 2 )六阶精度系数如下n = 掣，s ：八阶精度系数如下3 + 2 4 a 一6 臼5。：坠訾( 2 1 删。2 1 r 驯3 8 a 一96 9 6 1 1 9 1 a ，2 4 5 4 a 一2 9 41 1 7 9 d 一3 4 42 狐一，o2 1 r ，。2 1 矿，。可矿( 2 1 2 4 )十阶精度系数如下3 3 4 4 31 0 6 5 ，。2 丽，2 丽，o2 丽，o下面给出系数与截断误差的表一1 0 3 8 ，c ：旦( 2 1 2 5 ) 8 9 91 7 9 8，。“。，格式左模板右模板截断误差( 2 12 0 )35寻( 1 l a 一2 ) h 4 ，( 6 ( 2 1 2 1 )35 高h 6 ，( 8 ( 2 1 2 2 )57寻( 一2 + l l a 一1 2 4 f l + 2 0 c ) h 4 ，( 6 ( 2 1 2 3 )57鲁( 9 3 8 a + 2 1 4 f 1 ) h 6 ，( 8 )( 21 2 4 )57 笔3 氟浮h 8 o )( 2 1 - z 5 )57 西丽6 。1 9 。h l o ，( 1 2 )l e l e 关于半格点插值的逼近关系式如下；口( ， 一2 + 五+ 2 ) + a ( ，f l + + 1 ) + ， = ；( 十2 + f i 一 )+ ；( 五+ + 一 ) + ；( + + 一 )( 2 1 2 6 )当取方程左端为五节点模板，即卢0 的时候，四阶精度系数如下：。= ；( 9 + l o 。一1 4 f l + 1 6 c ) ，6 = ；( 一1 + 6 。+ 3 0 f f - 2 4 c ) ( 2 1 2 7 )2 0 0 4 年中国科学技术大学硕士毕业论文第8 页苎三兰! ! 竺兰童! ! 塞坌童兰皇苎竺兰竺奎兰塑；! 竺兰兰垒竺! ! 全六阶精度系数如下n = 击( 7 5 + 7 0 a - 4 2 既6 = 去( 2 7 0 p + 1 2 6 一2 5 ) ，c = 两1 ( 7 0 卢- l o a + 3 )( 2 1 2 8 )八阶精度系数如下：卢= 1 1 4 ( 1 r - - 5 ，n = 丁7 a + i 0 ，6 = 可1 8 9 a - 5 0 ，c = 等( 2 ，1 )十阶精度系数如下n = 努，口= 面5 ，a = ；，6 = 面5 ，c = 面1( 2 1 3 0 )下面给出系数与截断误差的表格式左模板右模板截断误差( 2 1 2 7 )56去( 3 一l o a + 7 0 z 一1 2 8 c ) h 4 ，( 4 j( 2 1 2 8 )56i 去( 5 1 4 a + 4 2 芦) h 6 ，( 6 j( 2 1 2 9 )56 吕 8 ，( 8 )( 2 1 3 0 )562 5 1 0 4 8 h l o ，( 1 0 )考虑到较小的模板也能达到我们运算需要的四阶精度，可以令原格式中的一些系数为零，可得各阶导数的小模板紧致格式的常用形式，这里我们将原格式的n 取的一致，这么做的优点将在后面体现：o ( 一o ( 一：，十矗。) + 一卸粤+ 。学，十一j 。) + 一= 5 五学点! 二! 盘墨二! 2n ( 一：。+ 群。) + 冉6 血型掣我们看如下的例子毗十“22 瓦“2 2第一步：解五对角阵“? + ( “。) ?第二步：解五对角阵“? + ( u x x ) ?+ a 业( 2 1 3 1 )，( 2 13 2 )2 五十l + 2 一l 一， 一22 h 3( 2 1 3 3 )2 0 0 4 年中国科学技术大学硕士毕业论文第9 页第二章p a d 6 型有限差分方岳和高阶发展方程2 2p a d 6 有理逼近格式第三步；( u d 产u n + m lr u n = 壶( u 。) ? 一n ( “。) ? ，得到预估值u n - t - 1 第四步；用辞“求( u 。) 。n - t - 1 ，( “。：) r 1第五步：c - n 格式：u n + 面l - - 一u 9 + n 垫立掣：，( 垫凿= 掣)以u r l 为中间层，这里相当于一个时间中心差，具有2 阶精度这里虽然得到了较小的模板，但是要解五对角方程，极大的影响了数值计算的的难度和效率，所以我们希望通过引进p a 挺有理逼近格式这种结构紧凑且精度高的有理分式逼近，在不减小模板的情况下，避免解五对角方程2 2p a d 6 有理逼近格式在这里我们采用一种紧凑的分式p a 挺逼近形式代替首先引入以下算子符号：t =，一1= 地“p = ( e ；十e j ，声= ( e + e 一1 ) ：6 = e 1 一一e ，6 u = ( t e q + 一ie - 1 ) u i 它们分别是移位算子，1 ee平均算子和中心差分算子。利用上述算子符号，上式又可以转化为p a d 6 逼近形式。1肺。儿一无丁巧i 利用p a d d 逼近，对于守恒方程“十，( u ) z = 0 ，( 2 21 )我们就有的半离散p a d 6 型格式毗+ ；# 舞 = o ，( 2 删毗+ 万南ko ，( 2 2 2 这时可以采用r u n g e - k u t t a 方法直接进行积分求解或者全离散p a d d型格式盟a 坐t + 三h 岛1 叫( 2 2 + ；6 2 “”p ”7可以利用差分方法【3 0 j 和三对角方程组的追赶法计算问题在于当方程中存在多项不同阶的数值导数时对于诸如b u r g -e 强k d v ，r l w 和k d v b u r g e r s 等非线性方程，如果不能够推导出相同商式的p a 雒逼近形式，则节点模板并不能节省，甚至可能增大。2 0 0 4 年中国科学技术大学硕士毕业论文第1 0 页第二章p a d 6 型有限差分方法和高阶发展方程2 3 小模板简化p a d 6 格式现以下例b u r g e r s 型方程加以说明：1“+ ，( “) z2 瓦“一( 2 2 4 )这是n a v i e r s t o k e s 方程数值研究中一个很好的模型方程如果采用一阶和二阶不同商算子的p a d 6 形式，对方程的空间导数进行逼近，则差分方程为+ ；南， = 上h 2 r e 南u ，( 2 。5 )啦。+ 元研拍2 一巧弘m 5 )和( 1 + 1 6 固( 1 + 去d 2 ) 蛳+ 元1 p 斯= 志乳，( 2 2 6 )显然，方程左边模板比较大，需要5 个节点，并没有达到前面所述的减小节点模板的效果，实际上仍然是五节点的模板，需要求解五对角系数阵的常微分方程组。或者需要采用分数步方法( 1 + 5 2 ) v i + ；p 哦= 志执。，( 1 + 去铲) u 。= 同样是大为增加了求解的困难。对于更高阶的发展方程，模板的结点数还要进一步的扩充分析格式，我们可以发现，模板增大的原因是p a d 6 逼近的分式中分母算子不统一，为此我们试图寻求一种统一的形式以便得到压缩的模板。2 3 小模板简化p a d 6 格式一般有理逼近或紧致格式的解法，是全隐式的方法譬如对于守恒方程 + c 厶= 0 ，( 2 3 1 )的求解过程，一般说来，可以分为下面的两个步骤：( 1 ) 在t = t “时间层上，根据前面的有理逼近方程，它可以写成矩阵形式a f = ；b f2 0 0 4 年中国科学技术大学硕士毕业论文第1 1 页第二幸p a d d 型有限差分方法和高阶发展方程2 3 小模板简化p a d 6 格式其中a ，b 是n n 矩阵，为t = t 。时间层上的未知节点 z = x d个数，f ，f 分别是相应的未知函数导数，及函数值向量。对于上述简单线性方程情况那是不必预先求解的；然而在通常，特别是非线性问题时，就必须预先解出f ( 2 ) 如果是上述的线性方程情况，可以直接写出有理逼近或者紧致( c o m p a c t ) 格式a j d f r ：；b f即使再进行时间离散，它也是完全隐式型差分格式。如果采用常微分方程的r u n g e - k u t t a 方法求解，它也是不方便的在非线性的情况下，则根据解出的f ，或者一m ) 代入方程( 8 ) 的差分格式，u 。n “= 工( u ? ，( n ? ) ) ，根据在第一节的讨论和分析，我们可以提出以下的小模板简化的p a d d 四阶逼近格式根据论文 3 2 】， 2 7 ，要求减小节点模板可以选取卢= 0 ，c = 0 。而且为了保证对于一阶和二阶微商，“得到同样分母算子的p a 拍逼近形式，下面取左端相同的系数n ，右端的差商系数则分别取n ，b 和i ，5a 一。一+ a 矗。= 。学+ n 学。一二。+ 一+ d 爿i ，= 5 盟挚+ i 纽学2o 1 万一+ o 万一( 23 ，2 )( 2 33 )又如果仅仅须要达到4 阶逼近精度，就只要满足下述系数关系n = ；( a + 2 ) ，6 = ；( 4 a - 1 ) ，a = ；( 1 一n ) ，5 = ；( 1 0 c r - 1 ) 可见，o l 的不同选择能够导致各种不同格式，我们考虑其中两种比较特殊的格式。( 1 ) 选b = 0 。则o = ，n = ，b = 0 ，a = 1 ，5 = ，就有p a d 6逼近形式一；南= 去鬻江。啕毗2 元研滞地。2 萨诱蚴心。4 )( 2 ) 选5 = 0 则o t = 南，o = ；，b = 一 ，a = 2 ，5 = 0 ，又有另一组p a d 6 逼近形式一；鬻t l i ，? - t x x 萨1 南协。2 0 0 4 年中国科学技术大学硕士毕业论文第1 2 页第二章p a d 6 型有限差分方岳和高阶发展方程23 小模板简化p a d 6 格式由于对于一阶和二阶微商f ，“”有了相同分母算子的p a d 6 逼近形式，所以可以对b u r g e r s 方程采用上面得到的两组p a d 6 逼近进行离散考虑到通量函数的复杂性，取第一组，得( 1 + ；固+ ；卢= 志( ；庐+ j 2 固“。( 2 3 6 )这是三对角的方程组也可以采用方程组的r k 迭代法求解现在，我们介绍一种”方案的全离散解法。记v i = “r 1 一“?”方法即进行如下两步过程( 1 + 兰6 6 2 ) q = 百a t p j f + 3 + 丢鼍( ；庐+ ；6 2 ) 恢，( 2 a 7 )“y 1 = “? + v i ，( 2 3 8 )其中+ j = ( ，( u ? ) + ，( “? + 1 ) ) ，u ? + = ( 婶+ u r l ) 。开始的，( “n + - )和“? “可以用，( “? ) 和n ? 代替，作为予估值；而后进行迭代校正格式是时间二阶，空间四阶的。考虑到方程有二阶微商项，一般要求a t = o ( h 2 ) ，所以可以说是空间四阶的。对于三阶的发展方程，我们依然采用 3 2 一文中的形式，选择达到四阶精度的卢= 0 三对角形式同前面的处理一样，对于三阶微商也取类似形式a ( 一+ 一+ - ) + 一= 。学+ n 学n ( 一! 。+ 一j 。) + 一= 5 墨生二+ in ( 一：l + 群。+ f ”k + i 一2 + 一15 血型喘攀+ a 血h 2_ 2 五+ 1 十2 ，i 一1 一 一22 h s( 2 3 9 )显然可以看出，如果。取值相同，在计算中可以缩小模板利用t a y l o r展开比较两边系数，为达到四阶精度，则需要满足如下关系：。= ；( 2 + n ) ，6 = ；( 一1 十4 。) ，a = ；( 1 一n ) ，6 = ；( 一1 + 1 0 n ) ，i ：2 ： n 一12 0 0 4 年中国科学技术大学硕士毕业论文铭1 3 页量三兰! 竺! 兰查! ! 圭坌童童皇童竺墨竺童兰! ! 兰苎堡望篓竺矍为得到较小的模板，如果我们选择6 = 0 ，代入上式解出相应的各系数得。= 互1 ，n = ；，a = ；，五= ；，5 = ；，a = 。，5 = o ( 2 3 a 一0 )。= 互，8 = 百，。= i ：岔= 吾1 。= ，8 = z d = u j从而可以得到一组p a d 逼近形式，写成算子形式就是，7 = ；南c 扣+ 和j2 瓦i 雨( 百p 6 + 石p 6 ) ，i= 去南c 跏厂7 = 嘉南扩 ，2 4高维问题的推广( 2 3 1 1 )( 2 3 1 2 )( 2 3 1 3 )其实，这种紧致格式向高维问题的的推广是很顺理成章的，我们可以很自然的想到利用分裂步法来处理这里先考虑简单的二维b u r g e r s方程：札f + “z + “v = 扩“2 z 十u v ( 2 4 1 )采用算子分裂方法，将上式写为li u + “u 。= z 。1三u 十“1 咕= 可v( 2 4 2 )然而我们必须注意到，继续采用上文的格式和四阶p a 挺格式我们可以达到时间二阶和空间四阶，对于此方程式本来并没有没有掉阶的问题。但是如果采用常用的分裂步法，对x 、y 方向各走 t ，这时候时间只有一阶精度，等于浪费了原格式处理方法的优点。所以此处我们采用s t r m l gs p l i t 的方法，x 方向先走 t ，然后y 方向走 础，然后x 方向再走 t ，依此类推直至最后一步x 方向走 t 这个方法的优点在于算子是可交换的，可以保持时间的二阶精度，避免掉阶问题的发生。每一步的半隐式格式如下构造：掣州h 瓦1 研# 5u 印击帮u ? ( 2 a 3 )2 0 0 4 年中国科学技术大学硕士毕业论文第1 4 页量三兰! 型! 竺三! ! 兰垒奎兰童童竺兰竺奎兰! 兰! 童竺望兰竺苎整理后可以化为如下等式左边为三对角矩阵，等式右边为五对角矩阵的形式：招n + l + ( j 2 + 石a t 学增+ 1 - f i l u 譬1 =西1 瓦a x “孙。+ ( 石1 + j 2 瓦a t ) u 耳。一；念u ? + ( 百1 + 。2 。a 。t ) u 。，+ 击瓮u 暑。( 2 4 4 )写成算子矩阵形式如下；a u n + 1 = b “?( 2 4 5 )考虑到边界条件的插值计算，00oooooo000000 0 ；十盛堡学 0 000ooo0ooooo0o000。00 吉蕊a t + ；爰一j 3 墨a t + ；爰南基0000000接下来我们考虑既是非线性又互相耦合的不可压n a v i e r s t o k e s 方程。由于密度p 是常数，连续性方程就变得简单了：赛十赛+ 丝o z = 。( 2 a 6 )a z a ，u 。、。此时粘性应力张量也大大简化，动量方程简化为：罾再( = 一j 1 v p 十舻v( 2 4 7 )2 0 0 4 年中国科学技术大学硕士毕业论文第1 5 页第二章p a d 6 型有限差分方法和高阶发展方程24 高维问题的推广假设动力粘度p 为常数，而v = 肛p 为运动粘性系数这里我们考虑二维平面流动方程简化为：丝o x + 骞= 。( 2 4 固锄“尝+ 。祟+ 。祟：一三关+ ，v 。瓦+ “瓦+ ”瓦2 一石磊+ 州祟+ 。宴+ 。丝：一三塑+ 。v 2 。o瓦+ u 磊+ 一y p o y + p v ”可以看出它有三个未知函数，三个方程为了便于处理常的涡一流函数方法，引进一个满足：0 世0 曲“2 面”2 一瓦的流函数砂( z ，y ) ，从而得到涡传输方程：裳+ 筹+ 筹= 舻e抚如a 或筹+ “塞+ “髻= v v 2 e其中a ua u、o yo x将( 2 ，4 1 1 ) 带入( 2 4 1 4 ) 可以得到流函数的p o i s s o n 方程：v z 世；雾+ 雾= e( 2 4 9 )( 2 a 1 0 )我们采用通( 2 4 1 1 )( 2 - 4 1 2 )( 2 4 1 3 )( 2 4 1 4 )( 2 4 1 5 )通过以上步骤，我们减少了一个方程和一个未知函数，现在仅仅需要解一个对流扩散方程和一个p o i s s o n 方程即可。对流扩散方程直接利用前面的格式即可接下来我们简要介绍一下多重网格法，详细请参阅【24 。在矩阵计算中，迭代法对摆动误差分量很有效，但是由于光滑分量的影响，收敛速度往往会放慢多重网格法协助迭代法消除光滑分量的影响，它的基本要点在于( 1 ) 误差可以划分为高频摆动分量和低频光滑分量；( 2 ) 设计特殊的地带方法消除摆动分量；( 3 ) 利用粗网格修正消除光滑分量。2 0 0 4 年中国科学技术大学硕士毕业论文第二章p a d 型有限差分方法和高阶发展方程第1 6 页25 算例所谓光滑分量与摆动分量是相对的，与网格尺度有关系在细网格上视为光滑的分量在粗网格上可能是摆动分量，从而可以用迭代法在粗网格上消除掉一旦迭代在细网格上收敛放慢，就表示误差已经光滑了，则转道较粗的一层网格上去，消除在粗网格上为摆动分量的误差这样一层一层直道消除各种误差分量，再一层一层返回到细网格上去实际的计算中往往采用几重网格，这里我们以二层网格的v循环为例看一下对方程的处理步骤。( 1 ) 在细网格上用松弛迭代法求解，l h u “= f “迭代”l 步后，计算所得近似值扩的残余：“= f “一l h v “( 2 ) 在粗网格上求解误差方程：l 2 h u 2 = 理“矿( 3 ) 进行粗网格修正： 6 + 喙口2 “一矿( 4 ) 以新的”为初值，在细网格上再迭代”2 步，然后回到( 1 ) ，开始下一个v 循环，直到达到一个收敛标准为止。其中工 代表细网格上差分算子，工2 h 代表粗网格上差分算子；毋是把细网格上的残余限制到粗网格上的算子；喙是把粗网格上的结果插值到细网格上的算子多重网格法在加速收敛的过程中并没有给计算带来很大的存储量一维问题时候，总存储量最多是单网格的两倍，；二维问题时候，总存储量不超过单网格的 倍。我们对p o i s s o n 方程作经典的中心差分并利用多重网格的解法。2 5 算例l 。b u r g e r s 方程的初边值问题：它可以作为流体动力学n a v i e r - s t o k e s 方程的简单模型方程，又可以代表浅水波问题的洪水数学模型，而且也是当代交通流动力学的模型方程通常，我们称为非线性输运扩散模型方程。如果在输运项再引入2 0 0 4 年中国科学技术大学硕士毕业论文第1 7 页篓三兰兰型! 墨童堡兰竺童兰皇苎竺兰竺童兰矍兰兰竺系数，那么通过调节方程的两个系数，可以进行输运和扩散的效应试验。如果采用下面的初边值条件，可以作为一种形成激波、边界层的难度很大的数值实验模型首先考虑如下特殊初值的b u r g e r s 方程初边值问题口21毗+ ( 号) z2 壶“一( 2 5 1 )u ( o ，z ) = ；s i n ( t r x ) + s i n ( 2 7 r x ) ，( 2 5 2 )u ( o ，t ) = “( 1 ，t ) = 0 ( 2 5 3 )这歌方程既有非线性对流项( 萼) z ，又含有扩散项击“。，r e = 1 0 4 。由于在上述的初值条件下解将迅速地形成激波，而且随着时间的进行发展成边界层，所以构造上述b u r g e r s 定解的差分格式要求很高，技巧性很强我们利用上面的p a 出逼近并采用”方法来进行为了保证移项后格式左边是三对角的对角优矩阵，我们选择前面提到的格式( 1 4 5 ) ，即选择。= 南，保证5 = 0 。令“。= “n + 酉l n = 矗，方程可以写成两种不同形式( 1 ) 紧致( c o m p a c t ) 形式杀+ 砜+ 去7 2 。+ 1