LaxWendroff时间离散的迎风紧致格式求解HJ方程.pdf

中山大学 硕士学位论文 lax-wendroff时间离散的迎风紧致格式求解h-j方程 姓名：范友宝 申请学位级别：硕士 专业：信息计算 指导教师：朱庆勇 20100529 l a x w e n d r o f f 时间离散的迎风紧致格式求解h j 方程 专业：信息计算 硕士生：范友宝 指导教师：朱庆勇教授 摘要 在本文中利用t a y l o r 级数展开式导出了一系列半离散的迎风紧致格式，紧致 格式具有网格基点少、精度高的优点。本文用l a x - w e n d r o f f 的时间离散的迎风 紧致方法求解h a m i l t o n - j a c o b i 方程，这是另一种时间离散方法，简称u c d l w , 从而代替传统的具有t v d 性质的r u n g e - k u t t a 时间离散方法。 u c d l w 格式基于迎风紧致格式计算光滑解所具有的高精度和算法简单，消 耗机时较少的优点而构成的。通过具体算例发现这种格式在处理某些具有不连续 导数等问题上不但减少机时的而且还保持基本无振荡，因此与r u n g e k u t t a 方法 比较，验证了u c d l w 格式具有误差更小，减少机时等优点。 关键词：h a m i l t o n - j a e o b i 方程，拉克斯温德罗夫，迎风紧致格式 u p w i n dc o m p a c ts c h e m e s w i t hl a x w e n d r o f f 略p et i m ed i s c r e t i z a t i o n sf o rh a m i l t o n j a c o b ie q u a t i o n s m a j o r ： n a m e ： s u p e r v i s o r ： i n f o r m a t i o na n dc o m p u t i n gs c i e n c e f a n y o u b a o p r o f e s s o rz h uq i n g y o n g a b s t r a c t i nt h i sp a p e r ，t a y l o re x p a n s i o nm e t h o di sa p p l i e dt od e r i v eac l a s so f s e m i - d i s c r e t eu p w i n dc o m p a c ts c h e m e s ，w h i c ht a k e sa d v 柚t a g e so fh i g ho r d e ra n d l e s sm e s hp o i n t s t h i sp a p e ru s e sl a x w e n d r o f ft i m ed i s c r e t i z a t i o nw i t hu p w i n d c o m p a c ts c h e m e st os o l v eh a m i l t o n j a c o b ie q u a t i o n s t h i si sa na l t e r n a t i v em e t h o d f o rt i m ed i s e r e t i z a t i o nt ot h ep o p u l a rt o t a lv a r i a t i o nd i m i n i s h i n g ( t v d ) r u n g e - k u t t a t i m ed i s c r e t i z a t i o n s ，t e r m e du c d - l w u c d - l ws h o w sa d v a n t a g e so ft h eu p w i n dc o m p a c ts c h e m e s ，s u c ha sh i g h a c c u r a c ya n dc o m p a c t n e s s c o m p a r i n gw i t hr u n g e k u t t am e t h o d ，t h ea d v a n t a g e so f u c d l ws c h e m e sa r el e s se r r o ra n dc p ut i m e a sar e s u l t ，w ef i n dt h i ss c h e m ec a n r e d u c et h ec o s ta n dk e e pn o n o s c i l l a t o r y c a p a b i l i t y f o rs o m ep r o b l e m sw i t h d i s c o n t i n u i t i e s k e y w o r d s ：h a m i l t o n - j a c o b ie q u a t i o n ；l a x w e n d r o f ft i m ed i s e r e t i z a t i o n ；u p w i n d c o m p a c ts c h e m e s n 原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独立进行研究 工作所取得的成果。除文中已经注明引用的内容外，本论文不包含任何其他个人 或集体己经发表或撰写过的作品成果。对本文的研究作出重要贡献的个人和集 体，均已在文中以明确方式标明。本人完全意识到本声明的法律结果由本人承担。 学位论文作者签名：范灰立 嗍蚍岁月7 日 学位论文使用授权声明 本人完全了解中山大学有关保留、使用学位论文的规定，即：学校有权保留 学位论文并向国家主管部门或其指定机构送交论文的电子版和纸质版，有权将学 位论文用于非赢利目的的少量复制并允许论文进入学校图书馆、院系资料室被查 阅，有权将学位论文的内容编入有关数据库进行检索，可以采用复印、缩印或其 他方法保存学位论文。 学位论文作者签名：溅导师签名： 日期脚辟岁月彳日 日期捌拜 日 ，职 中山大学硕十学位论文 第一章绪论 在本章中，我们将介绍h a m i l t o n - j a c o b i 方程( 简称为h j 方程) 基本理论及 其应用，并介绍了有关的发展过程和研究结果。第一节介绍了h j 方程的理论及 其应用第二节介绍h j 方程的发展过程、研究方法和一些相关的结果。第三节 介绍本文所做的主要工作 1 1h a m i l t o n j a c o b i 方程的理论与应用 h a m i l t o n - j a c o b i 方程应用于最优控制、计算流体力学、计算机图形图像、微 分几何、晶体成长、网格生成等许多领域，以及常用的界面追踪技术一l e v e ls e t 方法中所用的控制方程就属于h j 方程近些年来，许多学者对它给予了越来越 多的关注和研究 1 1 1h a m i l t o n j a c o b i 方程的基本理论 在1 1 维空间下一般形式的h j 方程为： l 织+ 日( _ x n ，t ，屯，欢) = 0 1 矽( x ，o ) = 丸( x ) 1 舢 方程( 1 1 ) 解通常是连续的，但是即使给定充分光滑的初值，它的解也可能 是不光滑的，在满足熵条件的情况下，方程( 1 1 ) 具有唯一解。h j 方程解的性 质非常复杂一方面，h j 方程的弱解存在但不唯一；另一方面，即使初值丸( z ) 和函数日足够光滑，其解的导数也可能会在某一时刻出现间断( 类似于守恒律 方程中的激波) 这些性质与守恒型方程的性质有许多相似之处实际上，对于一 中山大学硕十学位论文 维的h j 方程： 矽，+ h ( 矽工) = 0 ( 1 2 ) 令铭= 痧j ，那么方程就可以转换成等价的双曲型守恒律方程： u ，+ h ( “) 工= 0 ( 1 3 ) 而多维的h j 方程与弱守恒律方程组等价，对多维情形虽然没有相应的转换关系， 但也可以从形式上认为h j 方程是双曲守恒律的积分形式。 关于方程( 1 1 ) 解的存在性和唯一性问题，c r a n d a l l 等做了专门的描述。通常若 要确定h j 方程的某个弱解，需要对方程( 1 1 ) 增加一些额外的条件。比较常见 的具有物理意义的是粘性极限解( 简称粘性解) ，是通过粘性消去法【4 1 1 得到的，具 体来说考虑如下粘性方程： 等t + h ( x , t , c c , v c c ) - e a q y = o , ( x , t ) r n x ( 0 , 俐 ( 1 - 4 ) 沙( x ，o ) = 唬( x ) ，z r 这里占是一充分小的j 下常数，由于方程( 1 - 4 ) 属于抛物型方程，其解存在、唯一 并具有较好的正则性如果占趋向于0 ，粘性方程( 1 4 ) 的解就一致收敛到 矽( x ，f ) ，即： 矽( x ，t ) = 卿矽占( x ，f ) ( 1 5 ) 占 u 则c ( x ，t ) 为h j 方程( 1 1 ) 的粘性解，当然它也是h j 方程( 1 1 ) 的弱解 1 1 2h a m i l t o n j a c o b i 方程的应用 目前，h j 方程除了在原有的领域有着广泛的应用以外，例如几何光学，晶体 成长，化学反应，最优控制，图像处理、微分对策等方面在微分几何方面尤其是 平均曲率的曲面流问题都得到了广泛的应用。近些年来，人们已经通过水平集方 法又拓展了h j 方程的许多实际的应用领域 所谓水平集( l e v e ls e t ) 方法就是将n 1 维的界面看作为n 维空间的函数的 2 中山大学硕士学位论文 0 水平集，界面的运动就是用h j 方程来表示的。这种思想最早由物理研究人员 j a s n o w 、o t h a 以及k a w a s a k i 2 5 1 在1 9 8 2 年提出的，用于处理随机界面普适标度的 问题；1 9 8 8 年，o s h c r 和s e t h i a n 9 】运用h j 方程的理论将这种思想系统化并开始 实际用于处理流体界面，随后这几十年，o s h c r 、s u s s m a n 、p e n g 、s e t h i a l l 以及 m u l d e r 等对水平集方法的理论进行了广泛深入的研究，并将此方法用于解决许多 领域的实际问题，例如流体界面处理、网格生成和网格加密、燃烧界面、消去数 值噪声以及生产加工 透过这些领域，我们可以得到这样的结论：h j 方程的研究与应用正与r 俱增。 而由于h j 方程是一类非线性方程，在一般情况下我们很难得到它的解析解，这 将导致h j 方程的数值解方法的深入发展 1 2 h a m i l t o n j a c o b i 方程的数值方法 由于h - j 方程自身的发展以及它在许多领域有着重要的作用，尤其是近十年 来，l e v e ls e t 方法的发展以及应用，极大地促进了h j 方程的数值方法的研究。 这方面开创性的工作是由c r a n d a l l 以及l i o n s 2 明在1 9 8 4 年完成的，使用差分法求 解h j 方程： m + h ( 矽工，矽，) = 0 1 矽( x ，y ，o ) = 矽。( x ，y ) 1 6 令时问步，空间步分别为at ，ax ，ay ，并使用如下的差商记号： 6 ：多m = 吨叱一咖泔、) ia x ；6 二咖i ? k = 婶j k 一咖j k 、) ia x 万多矽，j = ( 矽+ l ，七一矽，j | ) a y ； 万歹矽- ，丘= ( 矽，+ l ，七一痧，七) y 方程( 1 2 1 ) 的差分格式的一般形式为： 咖；：= 咖：，k 一th 嫡：咖f ，k ，6 ：币j ? k ，6 ：毋j ? k ，6 j 多j k 、( 1 - 1 ) a 其中h a m i l t o n 函数日表示数值流通量满足相容性、单调性： 3 中山大学硕十学位论文 ( 1 ) 相容性：h ( a ，a ，b ，b ) = h ( a ，b ) ( 2 ) 单调性：h ( 口，b ，c ，d ) 关于a , c 单调下降，关于b , d 单调上升 这里日可取l a x f r i e d r i c h s 、或者具有熵修正的r o e 格式( i 强) 等类型的数值 h a m i l t o n 函数，例如对应于局部l a x f r i e d r i c h s ( l l f ) 型数值h a m i l t o n 函数的差 分格式为： 觎+ 兀咖) = 畔，半) 专一u - ) 一詈矿卅 ( 1 - 8 ) 其中 q = 删m 。a x 旷) j q ( u , v ) l ，口y = 刚m ( ，a 。x 一) 1 日2 ( “，v ) l q ( 甜，v ) = o h f ( u , v ) ，( “，v ) = o h 瓦( u , 一v ) ，地功= 跣毗6 ) 】 彳，曰】和 e d 】分别是“和矿的取值范围，如果将上式中的，( “+ ，u 一) 和 ，( 1 ，+ ，y 一) 分别替换为【彳，b 】和【c ，d 】，则可得到全局的l a x f r i e d r i c h s 格式 ( l f ) 由于上述差分格式为单调格式，故只能为一阶精度的，为了得到高阶的数值 格式，人们利用h j 方程与守恒律方程( 组) 的紧密联系，在此基础上构造出一 些新的格式。例如c o r r i a s 、蔚喜军【3 4 1 等利用h j 方程构造守恒律方程的数值格式， 而o s h e r 、s e t h i a n 以及l i 0 璐【1 8 1 等则通过应用守恒律方程的t v d t l 4 1 思想得到了求 解h j 方程的二阶差分格式，t v d 格式采用固定基点( s t e n c i l ) ，所以不管光滑区 精度多高，在间断处都降为一阶 19 8 7 年h a r t e n 等提出了e n o 2 l ( e s s e n t i a l l yn o n o s c i l l a t o r ys c h e m e s ) 格式，采 用与所求问题相关的自适应模板，根据各阶差商的绝对值极小选择方式，提高插 值方法精度而实现高分辨率和无振荡的效果，放宽了对t v d 的一些限制，克服了 ”巾格式在极点处降阶的不足，其优点是守恒的、保持t v d 的并且具有一致的 高阶精度，在间断区具有分别率高，在光滑区计算精度高等优点。缺点是：解及 其导数在零点附近产生的误差可能会改变模板的选择；而且在进行模板选择时， 为了得到k 阶精度的格式，需要k 个单元，总体覆盖了2 k 1 个单元，但是仅仅只 有一种形式的模板被使用，甚为浪费 4 中山大学硕士学位论文 1 9 9 4 年l i ux - d 等在e n o 基础上提出了本质加权无振荡格式( w e n o 2 4 1 ) 方法，其主要思想不是选择其中一个模板，而是利用模版的凸组合，通过引入变 化的加权因子，与e n o 格式相比在光滑解的截断误差阶数有了进一步的提高，而 在间断附近仍然保持良好的分辨率。近几年j q i u 、c w s h u 在w e n o 的基础上 提出了h e r m i t e w e n o t 8 1 格式和w e n o l w t 3 1 格式。w e n o 格式是目前最为精细的 计算格式，具有较好的稳定性，特别适合计算包含各类间断和含有多种复杂流动 的流场。但是这类格式计算过程复杂，计算量较大，阶数越高越复杂； 1 9 9 6 年，a b g r a l l 3 3 】在三角形网格上构造了一类解h j 方程的有限体积格式，此 格式是一阶精度格式。在比较弱的限制条件下，有限体积格式的数值解是收敛于 h j 方程的粘性解，并可以借助于e n o ( e s s e n t i a l l yn o n - o s c i l l a t o r ys c h e m e s ) 思想得 到了高阶格式。以后人们利用a b g r a l l 提出的数值通量，在每个三角形单元上构造 高阶插值多项式，得到了一个求解h j 方程的高阶精度格式。数值实验结果表明 该格式具有较高的精度和较好的分辨率 解h j 方程的另一种方法是非结构的有限体积法，相对于结构网格而言，非 结构网格具有如下优点：舍去了对网格节点的结构性限制，易于控制网格单元的 大小、形状及网格节点的位置，因此具有更大的灵活性，对复杂计算域的适应能 力很强，易于高精度处理边界条件，但是它的缺点是编程比较复杂，计算成本也 比较高，但是随着计算机技术的飞速发展，这些缺陷正在逐步消失有限体积法 以积分型守恒方程为出发点，通过对流体运动的体积域的离散来构造积分型离散 方程，这种方法更便于求解复杂的计算区域 1 3 本文的主要工作 本文在求解h j 方程中，介绍了迎风紧致格式的优点：是耗散型的格式，即 随着时间的推移解的幅值产生衰减，因此数值振荡被减少，同时迎风紧致格式利 用较少的网格点数，提高了计算效率。而在时间方向上利用l a x w e n d r o f f 时间离 散方法代替传统的r u n g e k u r a 格式，此格式可以相比同精度的r u n g e k u a a 格式 具有更小的误差，通过具体算列验证了格式的有效性 本文内容安排： 中山大学硕十学位论文 在第二章利用t a y l o r 级数展开得到了一系列半离散的差分格式，包括一般的 迎风格式和紧致迎风型格式，给出了差分格式的误差分析，以及差分格式的分辨率 的性质 在第三章提出了一种新的求解h j 方程的格式，即u c d l w 格式，考虑到迎 风紧致和l a x - w e n d r o f f 格式所具有的不同的特点，这是一种简单高效并且在重构 上具有紧致的优点 第四章我们把这种方法应用到具体的数值算例中，验证了此方法的正确性和 可行性 6 巾山大学硕十学位论文 第二章半离散的高阶迎风紧致差分格式 2 1 迎风紧致型差分逼近式 有限差分法是计算流体力学中最重要的数值方法之一，构造差分格式的方法 很多，近年来为提高计算方法的精度和分辨率，发展了多种新的差分格式的构造 方法，而对于微分方程的空间离散一般有两种方法：半离散与全离散。本文中采 用半离散的方法得到了一系列的迎风紧致型差分逼近式紧致型高精度格式受到 人们的关注在同样逼近精度条件下，紧致型格式更高的分辨率，且所需网格点较 少。对于求解多尺度复杂流动，要求数值方法能很好地捕捉到高波分量，并且要求 高波量的耗散误差及色散误差都比较小 2 1 1 差分逼近式的构造 数值方法的目的之一是寻找恰当的离散方法，并且根据原物理问题的初边值 条件，合理的构造离散化近似方程的初边值问题，使所求得的离散解以一定的精 度逼近原微分方程初边值问题的准确解 考虑如下双曲型方程 l 垒+ 望= 0 窃彘 ( 2 1 ) i “伉o ) = 这里为了讨论方便，假设厂( “) = au ，其中a 为常数，则方程( 2 1 ) 的半离散格式 为t 丝+ 曼= 0 一+ 上= 乱瓠 7 ( 2 - 2 ) 巾山大学硕士学位论文 这里c 血是一阶导数钞苏的逼近，而c 可由下面的式子确定： 口，乃+ ，= zb , ( d + ，+ 一乃+ ，) l| 将上式进行t a y l o r 级数展开，其中系数应该满足： = 岛( 这是相容性所要求) ， ( 2 3 ) ( s + 1 ) 口f z 3 = 岛 ，5 + 1 一( j 一1 ) 5 “ ，s = o ，l ，2 k l ( 2 4 ) f = 一刀l - - - - 一k l 可以利用上述格式构造对称型或者非对称型紧致差分逼近式，构造非对称逼近式 的目的是产生数值耗散，以压制数值解中非物理的高频振荡，它的具体形式依赖 于a 的符号，可利用级数展开方法根据精度要求来确定网格基架点数和式( 2 4 ) 中的系数。下面分两种情况给出几个高阶迎风紧致格式，对z 取不同的值，就会 得到不同的差分逼近式 2 1 2 显式迎风紧致型差分格式 在上式中如果取刀= o ，盔= l ，乞= 2 ，则方程系数满足： f = 睫1 + + 岛 一3 尻。一+ 勿= 0 1 7 n l + + 岛= 0 于是得到三阶显式迎风偏心格式： a 0 时， 1 1 f 2 警：q 3 n 专j l 一3 j o a 0 时， 巧= 丢瓯+ ( 一乃+ 。+ 5 乃一2 d 。) ( 2 5 ) ( 2 - 6 ) 中山大学硕七学位论文 类似的可以得到五阶显式迎风偏斜差分格式： a 0 时， 譬= 去酊( _ 3 如+ 2 7 f j + l + 4 砺1f j - i + ) ( 2 7 ) a 0 时， 巧= 未( 1 砜+ 弼一4 3 f j - - - i + 1 一) ( 2 - 8 ) 从以上逼近式中可以看出，差分方程中只含单个未知函数，则称为显式格式。 以上导数的逼近式中更多的利用了上游点上的函数值( 指相对于波的传播方向而 言) ，这类格式称为迎风偏斜格式这一点从物理上考虑是很容易理解的：因为波 是从上游传下来的 2 1 3 隐式迎风紧致型差分格式 如果令n 0 ，则( 2 - 4 ) 中的f 按隐式计算得到的差分格式就是隐式紧致 差分格式 在这里我们给出两个计算网格基架点更少的迎风型紧致差分格式： 三阶迎风紧致逼近式：( u c d 3 ) 詈譬+ 三弘陲山6x ) 乃 协9 ) 五阶迎风紧致逼近式：( u c d s ) 詈可+ 詈啄= 去( 一乃+ ：+ l1 l + , + 4 7 f j + 3 乃一) ( 2 加) 这里所说的迎风特性与式( 2 1 ) 中的波的传播方向相联系。上标“+ 表示 对应于a 0 嗣 - - 阶导数的逼近式；类似的方法可得到a 0 的逼近式，如此构 造的差分方程具有耗散特征一从上式中可以看出：f 可很容易按_ ，的增长方向求 9 中山大学硕士学位论文 搿竿；i 司样f j 口j 很谷易强y 阴藏j 、力l 司承j 睥。 紧致型差分逼近式与其他的一般差分逼近式相比在形式上可利用较少的网格 基架点得到较高的逼近精度，这是紧致一词的含义所指。在紧致逼近式中t 都以 隐式的形式出现，需要求解一个代数方程组，可用追赶法或者消元法求解。在这 里给出一个证明精度的定理 定理l ：紧差分格式丢哆+ + ；乃+ 丢e 一。= 霹乃是四阶逼近精度的 矾设巧是尼阶髋呜= 血 掣l + d ( 矿1 ) ，猁这里不 考虑时间方向，则乞= 够( _ ) + “妒+ 1 ) ，翌么泰勒展开可得： = 缸p ) + 埘”q ) + 三血y 3 ) + 昙缸y 4 q ) + 去血y 匀q ) + “时) = 缸p q ) 一对”q ) + 主缸y 3 q ) 一昙缸y 4 ) + 去崩y 5 q ) + d ( ) 磊= 鹏) + q ) 毛彬乃q ) 毛甜刃钙) 去鲫蝣) 咭蜊匀q ) + 斛) 扇= 鹏) 一缈吒) 毛鲫动q ) 之蜊q ) 毛彬吒) 蜊$ q ) + 斛) 设尺= 丢乃+ 。+ ；乃+ 吉弓一。一霹乃= 丢乃+ 。+ ；弓+ 丢乃一。一圭( 乃+ 。一乃一。) 代入上面的展开式，可得： r = ( 西1 一a x 5 厂( 。) + 。( 缸川) = 。( 血5 ) 因此k = 4 ，所以格式是四阶的 2 2 差分逼近式的误差分析 1 0 中山大学硕十学位论文 2 2 1 截断误差 将半离散的差分格式： 鲁+ 知 协 进行t a y l o r 级数展开，则得到一个对应的修正方程： 豢+ 望o = 哟缸，一- 粤(212)xa t |斟 、j 差分方程所准确逼近的微分方程，称为该差分方程的修正方程，修正方程与原微 分方程之差就是差分方程的截断误差截断误差的偶次导数项是耗散误差项，奇 次导数项则是色散误差项在截断误差的各阶导数项中，若最低阶导数为偶次则 该差分方程是耗散误差占主导，若最低阶导数为奇次则该差分方程是色散误差占 主导 考虑( 2 9 ) 这是一个三阶迎风紧致格式，- f n 证明从该式中得到的墨以三 x 阶精度逼近予婴，将上式经过t a y l o r 级数展开后，则得到： d 剪 ( f 一等篆+ i 山c 2 丽0 2 f 一百& c 3 万0 3 f + ) = 、 3 融6 融21 8 魂) 纠 ( 出生0 x 等等+ t 缸3 万0 3 f 一等等) ， 协 、 3 苏z6 苏j3 6 础珥 d 对上式两边求_ ，整理得到： 咖 娑= 垒，”一a x 2 ，+ 垡3 ，( 4 ) + 蜕”一a x 2f + 垡厂+ 。( 血4 ) 3 61 836苏 。、 。 把所得到的娑代入( 2 1 2 ) 的左端，整理得到： c + 等窘，”+ 州坊= 洫善+ 等等+ 等参+ 1 亿 中山大学硕士学位论文 对式( 2 t 3 ) 两边求爵3 2 ，并且将求得的丽3 2 f 代入( 2 3 ) 的左端墙整理得到： c = ( a x a 出y + 等等+ ) ， ( 2 - 1 5 ) 于是便知，互以三阶精度逼近于篓，将上式代入( 2 1 1 ) 可得： a x苏 0 u o f a x 3a 4f + 一= 一- = - = + + ( 2 1 6 ) t 3 ta x3 6 础4 一 。 2 2 2 分辨率性质 考虑如f 的模型方程及其半离散方程 害= 0 l 刊 ( 2 - 1 7 ) c 譬反 锄f + = o ( 2 - 1 8 ) a 缸 、 下面我们利用f o u r i e r 分析法对差分方法进行理论分析。 考虑厂似) = a u ，其中口为常数，初值条件“( x ，0 ) ：u o ( x ) ，考虑以2 石为周期 的解，初始条件可进行f o u r i e r 级数展开， 数。不难得到方程( 2 1 7 ) 的精确解 ( x ) = qe x p ( i k x ) ，这里k 代表波 七= 一 “ f ) = 靠e 印 玻 一刃) 】= ( x 一讲) k = - o o ( 2 - 1 9 ) 由此可知，对于具有不同波数的波，他们都以相同的速度a 传播 对于半离散方程，我们设x 方向网格步长为：缸= 等，乃= 缸，j = o ，l ， 设其初始条件为： 1 2 i ，山大学硕十：学位论文 n 2n 2 吁= q 州峰) ，c = 口吃州峨) ，k e = 砖+ 魄( 2 - 2 0 ) k = n 2k = - n 2 其中：砖= 玩 c o s 七( m + 1 ) ) c 一c o s 【加缸】) ， 忍= 尾 咖弘( m + 1 ) 嘲一豳 删) 则半离散方程( 2 - 1 8 ) 的解是 “归。笺：靠e x p ( 一口知e x p 愀( _ 一口惫明 ( 2 2 1 ) 为了方便讨论，我们仅取吩= e x p ( i 奴j ) ，那么此时半离散方程的解为： 啪) = 州一口惫f ) e x p 次( x j 一口去明 对比( 2 2 1 ) 和( 2 - 1 9 ) ，可知对于任意的尼，k a x o 因此格式是正耗散的； 在图( 2 1 ) 中给出了几种差分格式对应的岛和k ，随口的变化曲线。 色散误差k i k + d e l t - x 图( 2 1 ) 差分格式的k f 随a 的变化曲线 1 二阶中心，2 四阶中心，3 二阶迎风，4 三阶迎风紧致，5 五阶迎风紧致 1 4 口3s一口，厶s0o 屹 舶o+口sc j一 ，j 一心 所 以 一他量) l l l l l l 山大学硕十学位论文 耗散误差k r k d e l t - x 图( 2 2 ) 差分格式的k ，随倪的变化曲线 1 两点迎风( 一阶精度)2 三点迎风( 二阶精度) 3 四点迎风偏斜( 三阶精度) 4 六点迎风偏斜( 血阶精度) 5 三点迎风紧致( 三阶精度) 6 五点迎风偏斜( 五阶精度) 图2 1 中的曲线y 爿为t = 口= k a x 曲线它对应于微分方程的准确解。从图 中可以看出，对于低波分量，各种格式都可以较好的逼近准确解，与低阶精度格 式相比，高精度差分逼近式对中高波分量有较好的模拟能力。比较格式3 与5 、4 与6 可以看出，对于同阶精度的差分逼近式，与传统的差分逼近式相比较，紧致 差分逼近式对高波分量有更高的模拟能力。对中心型差分式也可以给出类似的t 随倪的变化曲线，得出类似的结论 在图2 2 中给出了几个差分逼近式对应的k ，随口的变化。它反映了格式的耗 散特性。对准确解及中心型格式k ，= 0 。从图2 2 中可以看出，高精度差分格式 有着更宽的低耗散波段。对于同阶精度的差分逼近格式，紧致型格式的低耗散波 段更宽 l i 山大学硕 ：学位论文 第三章l a x - w e n d ro f f 时间离散的迎风紧致格式 求解h - j 方程 3 1 介绍两种时间离散方法 前面的分析都是针对半离散方程而进行的，迎风紧致格式是一种空间的离散 方法，半离散化后得到的是沿时间方向的常微分方程，而在实际计算时必须对时 间导数进行离散，以下介绍了两种主要的时间离散方法 3 1 1l a x - w e n d r o f f 时间离散法 格式，这个差分格式在实际计算中得到了充分的重视，它与前面格式的推导过程 不同的是：除了采用t a y l o r 级数展开之外，还用到微分方程本身下面以时间与 空间上二阶精度给出具体差分格式 首先考虑双曲型方程： f 丝+ 口, 9 u = 0 c 3 t 砒 ( 3 1 ) 【u ( x ，o ) = u o ( x ) 为简单起见，假设为等距网格剖分，设x ，= j a x ，乙= n a t ，并且设坂五力是 微分方程的光滑解，将u ( x ，乞+ 1 ) 在点( x j , t 。) 处做t a y l o r 展开： “( _ ，乙+ 。) = “( ，乙) + f a a - “t , , , , 2 f 2l r a a t 2 “2 j 砌+ 。( f 3 ) c 3 2 ， 1 6 中山大学硕士学位论文 o uo u o to x 伊优0 ，o u 、，铲“ 可= 瓦一夏) 钏丽 将( 3 - 3 ) 、( 3 4 ) 代x ( 3 2 ) 中有： ( 3 3 ) ( 3 - 4 ) “q ，= 螂，伊础陲l a x $ 一+ 2l 髟上1 刃+ 啦3 ) ( 3 5 ) 再对x 方向上采用中心差分逼近上式中的导数项，例如二阶中心差分得到： 阻疗= 如碱一( ) 】+ ( 3 - 6 ) 匿卜爵1 砧一姚+ 峨舢+ 彬) 俘7 ) 代入( 3 5 ) 化简并略去高阶项，整理可得到如下的差分格式： = 哆一互a 石a tl ，一略) + i a 2 孬a t 2 ( 略。一衅+ 啄，) ( 3 8 ) 从差分格式的构造可以看出( 3 8 ) 式在时间和空间方向均具有二阶精度，我们称 ( 3 8 ) 式即为l a x - w e n d r o f f 格。现在我们令“；= g 玎e x p ( i k j a x ) ，代入( 3 - 8 ) 式，得到 g = l 一2 搠2 s i n 2 竽嘞心n 后缸( 3 - 9 ) 这里名：口会上，于是如果满足条件f 口1 名1 ，那么就有ig l 1 ，这也就是 x 。 l a x - w r e n d r o f r 格式的稳定性条件。 从以上可以看出，可通 z - tt a y l o r 级数展开式达到时间上更高精度的要求，相 比下面介绍的r u n g e - k u t t a 离散方法，l a x w r e n d r o f j f 时间离散法能够用较少的有效 节点达到高精度，但需要多次用到微分方程本身。近期c w s h u 1 6 】等人提出了 l a x w b n d r 0 1 j f 时间离散方法结合高精度的w e n o 格式，即w e n o l w 格式。此后 1 7 中山大学硕士学位论文 3 1 2r u n g e k u t t a 时间离散法 詈乩翩) ( 3 1 0 ) 这里u 为矢量，l 空间离散算子，对式( 3 - 1 0 ) 可采用通常的求解常微分方程 的求解方法r - k 方法，最简单的是时间方向采用一阶逼近： u 斛1 = u 刀+ t l h ( “一) 这一时间格式精度太低，而且在很多情况下格式是不稳定的，一个具有二阶精度 的r - k 方法是： f ，“开+ = “一+ l 。x t l ( m 一) 2 7 i t = “行+ 圮翩斛 ) 如果所求方程为单波方程，而对应的一阶导数取具有二阶精度的中心差分，即： 三 ( u ) = ( ，一盟a x ，) ( 3 1 1 ) 则格式的稳定性条件为c f l 1 。 u 1 = u 打+ f 三| | i ( “刀) ，( 2 ) 一 “ 一 刀+ 1 一 “ 一 四阶非t v d 型的r u n g e - k u t t a 法，简写为( 砌“) ： 1 j 】 ) l ( ，l “ 甜 0 ( 三， 纪 力 f + 十 、，二_ ， 2 “ 甜 一4 2 3 + + “ 铭 34一3 中山大学硕十学位论文 “1 = “一+ 1 l _ _ at l 厅( “疗) “2 = “露+ l ， at l ( “1 ) u 3 = u 疗+ a t l | i l ( “2 ) un + l = 了1 ( 一“以+ 掰( 1 ) + 2 “( 2 ) + “( 3 ) + 1 6at l j | i ( “( 3 ) r u n g e k u t t a 方法利用常微分方程求解器( o d es o l v e r ) ，它具有算法简单、易编程 等特点，而且还具有补巾的稳定性，因此大多数的迎风紧致格式、w e n o 空间 离散化得到的常微分方程均采用r u n g e - k u t t a 时间离敖法， 3 2l a x - w e n d r o f f 时间离散的迎风紧致格式 3 2 1 构造u c d l w 格式求解h - j 方程： 首先给出一维哈密顿一雅闲比万程： 办+ 日( 丸) = 0 矽( x ，o ) = c o ( x ) 3 _ 1 1 ( 1 ) 一维情况： 我们定义矽表示矽在时间方向上的r 阶导数，由泰勒展开式： 妣 锄= 谚+ 锄+ 等谚”+ 等谚m + 等矿+ 斗等妒，2 ， 如果我们要获得在时间上的k 阶精度，则需要分别逼近矽，矽，”矽”矽，这 里时间上我们只展开到三阶精度，当然可按t a y l o r 展开达到更高精度。 从3 1 1 我们可以得到： 1 9 中山大学硕十学位论文 矽= - h ( u ) 矽。= 一h ( “) 丸。 矽。= 一h ”( “) ( 旌) 2 一h ( “) 丸。 3 。1 3 矽4 = - h ”( “) ( 丸) 3 3 h ”( “) ( 丸) 织”- h ( “) 丸 ( 这里“= 丸) 第一步：构造一阶时间微分项：矽= 一日( 跣) 通i 救1 徽 f t a 谚=一局(3-14) 童是数值流通量，满足如下几个性质： ( 1 ) h al i p s c h i t z 蝴； ( 2 ) 具有相容性( c o m i s t e m y ) ，即日( “，u ) = h ( u ) ； ( 3 ) 具有单调性丘( 个，山) ，即对第一个变量非减，第二个变量非增； 这里我们用l a x f r i e d r i c h s 通量分裂格式： e = 日( 华) 一詈( 。卅 ( 3 - 1 5 ) uf + ，uj 一分别表示从左侧， 右侧逼近于“( 再，t ) ，“产= 矽五 ， 口= ma xl h ( “) l ，这里用三阶或者五阶迎风紧致格式逼近于“，+ ，“，一。因此 就可以得到痧在西点的逼近值：( 破) ；= ( u 7 + “i ) 2 ，并且谚就可以通过( 3 1 3 ) 第一式得到 第二步：构造二阶时间微分项：矽”= 一h ( 矽x ) 矽二 由于日 ) 已知，只需求出丸。由于第一步中已经求出办，只需在x 方向 上用中心差分逼近就可以求出败，例如当甜嘉、咭、圬用三阶迎风紧致差分时， 由于应的作用，只需用二阶精度的中心差分逼近丸： 中山大学硕士学位论文 ( n 错 ( 3 1 6 ) 当砖、咭、咭用五阶迎风紧致差分时，由于& 的作用，只需用四阶中心差分逼 近丸： ( 妣立垫警l 血 ( 3 - 1 7 ) 这样就得到了在玉点的逼近值： j = - h ( ( j ) ，) ( j ：) ， 第三步：构造三阶时间微分项矽”= 一日( 矽，) ( 娥) 2 一h ( 矽，) 矽j ( 3 1 8 ) 与第二步类似，只需求出( 苁”) ，利用第二步求出的痧”，f l j - t - a t 2 的作用， 在x 方向只需用中心差分逼近就得到矽”在蓐点的逼近值： 辫”= - h ”( 丸) ( 苁) 2 - h ( 丸) ( ) ( 3 1 9 ) 第四步：四阶时间微分4 ) = 胡( 破) ( ) 3 3 h ”( 丸x ) 一h ( 破) ” 与第三步类似只需要求出，这里表示矽在( x ，t 彪) 点的值，在第三步 中已经求出，由于a t 3 的作用，只需在x 方向上用二阶中心差分逼近即可求出， 最终我们得到痧4 在t 点的逼近值： 谚4 = 一日( 丸) ( 丸) 3 3 h ”( 丸) ( 丸) ( 丸”) 一h ( 丸) ( 丸) ( 3 2 0 ) 如果我们要求在时间方向上达到更高精度的要求，可通过上述类似方法构造。下 一时间步的最终逼近式如下： 地 = 办地“等“等办。+ 等) + + 等) ( 3 捌) ( 2 ) 二维情况： 现在我们继续考虑二维的h j 方程，为简单起见，我们仍然假设为一致网格剖 分，并且血= 毛+ 1 一再，少= y i 十l y i ，从( 3 1 ) 在二维的情况下，我们有： 2 l 中山大学硕士学位论文 矽- - - - h 矽爿锹坳 矽叫1 删酬- 瞰- 呦 妒叫。痧吗砝蟛一嘶屹矽 啤触q 黝嗍j 吗钐雄嘞 ( 3 2 2 ) 这里h i 表示日在f 方向的偏微分，h i ，表示在耖方向上的二阶偏微分，表 示日在沙方向上的三阶偏微分，与一维类似，我们要从点集 办= 痧( 西，乃，f 一) 分 别重构苁，丸，吮，吮，纯，盛和歹，夕，矽f t t ，令“= 识，= 丸那么一阶时间微分：别重构苁，丸，吮，吮，纯，盛和矽，矽，矽，矽钔，令“= 识，= 丸那么一阶时间微分： 矽= 一h ( “，) ( 3 2 3 ) 可通过如下格式逼近： 盒：日( 堕笋，堕笋) 一号( 蟛一啄) 一q 一丐) ( 3 - 2 4 ) 这里“；表示从左侧和右侧逼近于点矽j ( x ，y ，f ) ，咭表示从下面和上面逼近于点 丸( 一，y ，f ) ，其中： 口工= s u p i h l ( “，1 ，) l ，口y = s u p l h 2 ( “，v ) l 材；和污可通过类似一维的情况第一步中用三阶或者五阶紧致差分逼近得到，并且 和一维情况类似我们可以得到矽j 和咖在点( 薯，乃) 的逼近值： ( 苁) 扩( 磅+ 啄) 2 ，( 方) 矿( 咭+ 丐) 2 ( 3 2 5 ) 另一方面我们可以通过中心差分逼近得到如下微分：矽，矽， ，c y , 矽，川， 因此就可以通过( 3 2 2 ) 式求出矽”，矽，矽4 ) 的逼近值 ( 3 ) 下面我们把u c d l w 算法推广到三维情况： 为简单起见，我们仍假设为均匀网格剖分即x = 一再，缈= 心。一咒， z = 气一乃均为常数，在三维的情况下，由( 1 1 ) 式我们有： 中山大学硕士学位论文 牵= 一h 矽”= q 东一月拶一马谚( 3 - 2 6 ) t - - q 。硎) 2 一厶乞) 2 一马，洌) 2 2 q 2 一2 q 3 东刃一2 屹砒一q 一上硝一马 这里h i 表示日在f 方向的偏微分，h i _ ，表示日在耖方向上的二阶偏微分，我 们仍然要从 独= 矽( t ，y j , z k , f 以) 分别重构破，方，唆，九，丸，以，识”，方”，破” 和矽，矽，矽，令“= 织，v = 力，w = 统，那么一阶时间微分可通过l f 通量分裂 格式得到。“；、崂、咭可用三阶或者五阶紧致差分逼近得到，我们可以得到识、 矽，、识在点( x i ，y ，z ) 的逼近值： ( 识) 独( “未+ 啄) 2 、( 吮) 驰( 壤+ 嚎) 2 、( 破) 独( w 彘+ 磊) 2 另一方面通过已求出的矽分别在叉，少，z 方向上用中心差分逼近得到 疵，丸，矽：，例如当“；、