8-高阶紧致格式

1§10.高阶紧致差分格式先考虑导数的差分近似。若某一差分近似的精度是p阶的，则近似的误差就是()ph。要想进一步提高精度，通常有两种途径：减小h（h-version）或是提高p（p-version）。但由于计算机资源的限制，h不可能无限地减小，因此在需要高精度流场计算的情形（如，粘性边界层、湍流等），就要考虑采用高阶格式。通常情形，构造高阶格式需要更多的点。例如：两点差分近似()()()fxhfxfxh+-¢»只有一阶精度。而使用三个点，就可以构造出二阶近似()()()()2432fxhfxhfxfxh-+++-¢»精度越高，需要的点就更多。对于中心差分近似也有类似的结果。但是这种高阶近似用在差分格式中，除了计算公式更加复杂，计算量增加之外，还会造成其他困难。例1：以一个简单的常微分方程初值问题为例。设0a。0duaudx+=（01x?），()0u=取M个网格，空间步长1hM=，网格点记作jxjh=（0,1,2,,jM=L），网格点上的近似解记作()jjuux»。2因0a，导数采用向后差分近似，就有10jjjuuauh--+=（1,2,3,,jM=L）实际的计算方案为0u=，111jjuuha-=+（1,2,3,,jM=L）上述格式用到两个点，但只有一阶精度。如果采用二阶差分近似，则成为12340jjjjuuuauh---++=（2,3,,jM=L）这个格式具有二阶精度。可是由于涉及三个点，所以只能从2j=开始计算。而初始条件只提供了0u=。因此1u的计算就需要补充另外的等式。对于更为复杂的流动控制方程以及更复杂、精度更高的数值格式，这种问题就更加严重。现在我们从另外一个角度来考察上述问题。将导数的近似值记作jjduudx¢»，则差分格式就可写成0jjuau¢+=我们刚才所做的不过是用不同的差分来代替ju¢。因此，我们遇到的困难就是：用高阶差分代替ju¢，就会涉及更多的点。而我们的问题也就是：有没有不涉及更多点的高阶差分？3我们借助算符演算来讨论这一问题。例2：由1IE-?-可推出1EI-=-?，于是有()1234234111lnlnln11112341111234DEEIhhhhh-==-=--?骣÷ç÷=--????ç÷ç÷ç桫骣÷ç÷=????ç÷ç÷ç桫LL上式右端取第一项，就得到一阶差分近似()111DIEhh-谎=-，即()()()fxfxhfxh--¢»如果取前两项，就得到二阶近似()()()()22111121211122222234DIEIEhhIEIEEIEE-------骣轾÷ç÷谎+?-+-犏ç÷ç÷ç犏桫臌=-+-+=-+即()()()()2432fxhfxhfxfxh-+++-¢»这些就是前面用到的向后差分近似。但如果继续演算，有2341111234Dh骣÷ç÷=????ç÷ç÷ç桫L4234232342434111111122342111111234211124611111121202IhIhIhI骣骣鼢珑鼢=-蜒+???珑鼢珑鼢珑桫桫-?é骣÷çê÷=????ç÷êç÷ç桫ë-?ù骣÷çú÷-???ç÷úç÷ç桫û骣÷ç÷=?+??ç÷ç÷ç桫-ÑÑ×LLLL上式中2Ñ的系数为零，因此取第一项相当于取了前两项，也能得到二阶精度的近似。即112DhIÑ»-?注意到此式中只出现了Ñ的一次方，因此只涉及两个点。上面导出了一个新的差分近似，是用差分算子的有理分式表示的，因此称为微分算子的有理函数近似（Pade逼近）。而通常的差分近似都是用多项式表示的。例3：由()()()()1122233233331314124IIEEIEIE--轾犏=--=--+犏臌轾犏=--++犏臌5()3133511244EIE-轾犏=---+=+犏臌()()()()()()112245545555212551212552455791161464161221464161111216216IIEEIEEIEEEIEEEIEE------轾犏=--=--+犏臌轾犏=--++++犏臌éê=---+êëùú--+-+úû骣÷ç÷=++=++ç÷ç÷ç桫于是有()()()35135579351113sh6401113116440216111630Dhhhh-轾犏==-+-犏臌轾骣骣鼢珑犏鼢=-++++-珑鼢犏珑鼢珑桫桫臌骣÷ç÷=-+-ç÷ç÷ç桫LLL上式右端取第一项，就得到二阶精度的中心差分近似()()1112DEEhh-?-，即()()()2fxhfxhfxh+--¢»而取前两项，就能得到四阶精度的中心差分近似6()321211188612DEEEEhh--骣÷ç÷?=--+-ç÷ç÷ç桫即()()()()()288212fxhfxhfxhfxhfxh+-++---¢?但又有352352352352111630111111663061111163061161603611DhIhIhIhI骣÷ç÷=-+-ç÷ç÷ç桫骣骣鼢珑鼢=+-+-珑鼢珑鼢珑桫桫+é骣÷çê÷=-+-ç÷êç÷ç桫ë+ù骣÷çú÷+-+ç÷úç÷çèû=-+×øLLLL531180骣÷ç÷+-ç÷ç÷ç桫L和前一个例子一样，上式中只取第一项，就能得到四阶精度的中心差分近似2116DhI»+而且该差分近似只涉及三个点。7以上的讨论表明，有理函数近似可以达到我们原来的目的，即：有理函数近似具有更高的精度，又不涉及更多的点。下面考虑微分算子有理函数近似在数值格式中的应用。这种有理函数的表达式只是一种算符操作，在实际应用中就需要将有理分式化为整式，过程如下。例4：由112DhIÑ»-?有112IDh骣÷ç÷-鸦?ç÷ç÷ç桫作用在函数()fx上，()()112IDfxfxh骣÷ç÷-鸦?ç÷ç÷ç桫即()()112Ifxfxh骣÷ç¢÷-鸦?ç÷ç÷ç桫将算子展开，就是()()()()2fxfxhfxfxhhⅱ+---»对中心差分近似也有类似地的结果。8例5：由2116DhI»+有2116IDh骣÷ç÷+?ç÷ç÷ç桫作用在函数()fx上，()()2116IDfxfxh骣÷ç÷+?ç÷ç÷ç桫即()()2116Ifxfxh骣÷ç¢÷+?ç÷ç÷ç桫将算子展开，就是()()()()()462fxhfxfxhfxhfxhhⅱ?+++-+--»以上两个例子表明，有理函数给出的差分近似，会同时有多个点处的导数值出现，需联立求解。而通常的差分近似，只出现一个点处的导数值，可逐点计算。这两者之间的区别，类似于隐式格式与显式格式的区别。正因为如此，微分算子的有理函数近似也称为隐式差分近似。同时，由于涉及较少的点，通常又称为紧致差分近似。例6：将例4中的紧致差分近似应用于例1中给出的初值问题，91102jjjjjjuauuuuuh--ì¢ï+=ïïïïíⅱ+-ïï=ïïïî（1,2,3,,jM=L）整理后，得到未知解的近似ju及其导数值近似ju¢的联立方程组11022jjjjjjuauhuuhuu--ì¢ï+=ïïïíïⅱï-+=+ïïî解得11112222jjjjjjuhuuhaauhauuha----ì¢ï+ïï=ï+ïïíï¢+ïï¢=-ïï+ïî（1,2,3,,jM=L）对于0j=，利用原方程可给出初值0u=，00uaua¢=-=-由此可见，在紧致差分格式的求解过程中，未知解的近似及其导数值的近似都是未知量，是需要联立在一起求解的。上面的例子是一个两点紧致格式，最终得到了一个递推关系式，逐点计算。对于涉及三个甚至更多点的高阶紧致格式，就需要将未知解的近似1u、2u、L、Mu及其导数值的近似1u¢、2u¢、L、Mu¢（如果原方程还包括二阶导数，则还有二阶导数值的近似1uⅱ、2uⅱ、L、Muⅱ）全部放在一起联立求解。因此，高阶紧致格式中需要求解的未知量比较多，这是它的一个弱点。10下面列出一阶导数和二阶导数高阶紧致差分近似的一些结果。1.Pade逼近（三点四阶）11111111246210212jjjjjjjjjjjuuuuuhuuuuuuh+-+-+-+-ⅱ?++-=ⅱⅱⅱ++-+=2.对称紧致格式（五点六阶）11221111141339294jjjjjjjuuuuuuuhh+-+-+---ⅱ?++=+3.对称紧致格式（五点八阶）2112112214413699364025272544jjjjjjjjjuuuuuuuuuhh++--+-+-ⅱⅱ?++++--=+4.迎风紧致格式（三点三阶）1114521336jjjjjuuuuuh+--+-ⅱ+=5.迎风紧致格式（五点五阶）211211236443325560jjjjjjjuuuuuuuh++----++--ⅱ+=6.广义紧致格式（对称三点六阶）上面给出的紧致差分近似，计算一阶导数的紧致差分里不会出现11二阶导数的近似值，计算二阶导数的紧致差分里也不会出现一阶导数的近似值。如果突破这个限制，就成为广义紧致差分近似。例如()()()()()()2111112111111488242413948904871673901515jjjjjjjjjjjjjjjhhuuuuuuuuhhuuuuuuu++-+-+-+-+-ⅱⅱⅱ?--+-+-+=ⅱⅱⅱ?--+++-=7.广义紧致格式（迎风两点三阶）1112111201220122jjjjjjjjjjjjuuuuuuhhuuuuuuhh+++---ⅱⅱⅱ-+--+=ⅱⅱⅱ-+--+=12最后给出一个实例。例7：考虑Burgers方程（对流扩散方程）两点边值问题2210,01,xxuuuaxtxxuuuUx==ìï抖?ï+=ïï抖¶ïïïíï镦?¶ï÷ç÷?==ç÷ïç÷ç¶镨?ïî将空间区域0,1轾犏臌均匀划分成M个网格，则空间网格的尺寸为1hM=，网格点坐标为jxjh=（0,1,2,,jM=L）。在ntnt=D时刻（1,2,3,n=L，tD是时间步长），将未知函数及其空间导数在网格点上的近似值分别记作(),njjuuxt»，njttjxxuFx==¶»¶，22njttjxxuSx==¶»¶现假设上一时刻()11ntnt-=-D的近似解已经求出，记成()1,njjbuxt-»，在计算过程中视为已知。于是，在空间区域0,1轾犏臌内的第j个网格点（1,2,,1jM=-L）处，有原方程的差分近似0jjjjubaFSt-+-=D和空间导数的Pade逼近131111111124026210012jjjjjjjjjjjuuFFFhuuuSSSh+-+-+-+--++-=-+++-=在左边界0j=处，有边界条件00uF+=原方程的差分近似00000ubaFSt-+-=D以及空间导数的广义紧致格式10101020212uuFFSShh-+--+=在右边界jM=处，有边界条件MuU=在此处，原方程成为0MMaFS-=还有空间导数的广义紧致格式11120212MMMMMMuuFFSShh----+--+=将所有这些集成在一起，就得到线性方程组140000002200011111112221111211021221240333520126120212jjjjjjjjjjjjjjjMMMMMuFuatFtSbthhhhuFSuFSuatFtSbthhhuFFuFhhhuSuSuSuUatFtShhuF--++--++--+=+D-D=D---+-+=+D-D=D---+-=---+-==D-D=---MM210212MMMMhhSuFS-ìïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïíïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïï+-+=ïïî