10-高阶紧致格式

§10. 高阶紧致差分格式 10.1 高阶差分

先考虑导数的差分近似。若某一差分近似的精度是 p 阶的，则近似的误差就是 ()

p h O 。要想进一步提高精度，通常有两种途径：减小 h （h -version ）或是提高 p （p -version ）。但由于计算机资源的限制，h 不可能无限地减小，因此在需要高精度流场计算的情形（如，粘性边界层、湍流等），就要考虑采用高阶格式。

构造高阶格式需要用到导数的高阶差分近似。通常情形，这需要更多的点。例如：两点差分近似

()()()

f x h f x f x h

+-￠?

只有一阶精度。而使用三个点，就可以构造出二阶近似

()()()()

2432f x h f x h f x f x h

-+++-￠?

精度越高，需要的点就更多。对于导数的中心差分近似，也有类似的结果。

但是这种高阶近似用在差分格式中，除了计算公式更加复杂，计算量增加之外，还会造成其他困难。

例1：以一个简单的常微分方程初值问题为例。设 0a > 。

0du

au dx

+= （01x < ） ， ()0u =α

取 M 个网格，空间步长 1

h M

=

，网格点记作 j x jh =（0,1,2,,j M =L ），网格点上的近似解记作 ()

j j u u x ? 。

因 0a > ，导数采用向后差分近似，就有

1

0j j j u u au h

--+= （1,2,3,,j M =L ）

实际的计算方案为

0u =α ， 11

1j j u u ha

-=

+ （1,2,3,,j M =L ）

上述格式用到两个点，但只有一阶精度。如果采用二阶差分近似，则成为

12

340j j j j u u u au h

---++= （2,3,,j M =L ）

这个格式具有二阶精度。可是由于涉及三个点，所以只能从 2j = 开始计算。而初始条件只提供了 0u =α 。因此 1u 的计算就需要补充另外的等式。对于更为复杂的流动控制方程以及更复杂、精度更高的数值格式，这种问题就更加严重。

现在我们从另外一个角度来考察上述问题。将导数的近似值记作

j j

du

u dx ￠?

，则差分格式就可写成 0j j u au ￠+=

我们刚才所做的不过是用不同的差分来代替 j u ￠ 。而我们遇到的困难是：用高阶差分代替 j u ￠ ，就会涉及更多的点。因此，我们的问题就是：有没有不涉及更多点的高阶差分？

我们借助算符演算来讨论这一问题。 例2：由 1I E -?

- 可推出 1E I -=- ，于是有

()

1

2

3

4

2

3

4

111ln ln ln 11112

3

4

11

1

12

34

D E E I h h h h h -==-=--

骣÷?÷=--?

????÷÷

?桫骣÷?÷=?????÷?÷

桫

L L 上式右端取第一项，就得到一阶差分近似

()

1

1

1D I E h

h -谎

=- ，即 ()()()f x f x h f x h

--￠? 如果取前两项，就得到二阶近似

()()

()()

2

2

1111212

11

1

22

222234D I E I E h h

I E I E E I E E -------骣轾

÷?÷谎+?-+-犏?÷?÷

犏桫臌

=-+-+=-+

即

()()()()

2432f x h f x h f x f x h

-+++-￠?

这些就是前面用到的向后差分近似。

但如果继续演算，有

2

3

4

2

3

4

234

232

4

11112

3

4

11

111112234

2111

1112

342

111 24

6

1

12

01D h I h I h

I h

I 骣÷?÷=?????÷?÷

桫骣骣鼢珑鼢=-蜒+???珑鼢珑鼢桫桫-

é骣÷?ê÷=?????÷ê?÷桫?

-

ù骣÷?ú÷-????÷ú?÷桫?

=?-

籽L L

L L 3

4

1

112

12

骣÷?÷+?

??÷?÷

桫

L

上式中 2? 的系数为零，因此取第一项相当于取了前两项，也能得到二阶精度的近似。即

1

12

D h

I ?

?

-

注意到此式中只出现了 ? 的一次方，因此只涉及两个点。

上面导出了一个新的差分近似，是用差分算子的有理分式表示的，因此称为微分算子的有理函数近似（Pade 逼近）。而通常的差分近似都是用多项式表示的。

再举一个中心差分的例子。

例3：由

()

(

)

()

()

1

1

22

2

3

3

2

33

3

3

13

141

24I I E E I E I E --轾犏=--=--+犏臌

轾犏=--++犏臌μδμδμμδμδμδμδμδ

()

3

133511244E I E -轾犏=---+=+犏臌

μδμδμδμδ

()

()

()

()

()

(

)

112

2

4

5

54555

5

2125

5

1

21255

2455791161

464161

221

46416

1111216216I I E E I E E I E E E I E E E I E E ------轾犏=--=--+犏臌

轾犏=--++++犏臌

éê=---+ê?

ùú--+-+ú?

骣÷?÷=++=++?÷÷?桫

μδμδμμδμδμδμδμδμδμδμδμδδδμδμδμδ

于是有

()()()35135579351113sh 640

1113116440216111630

D h h h h -轾犏==-+-犏臌轾骣骣鼢珑犏鼢=-++++-珑鼢犏珑鼢桫桫臌骣÷?÷=-+-?÷?÷桫L L L μδμδμδμδμδμδμδμδμδμδμδμδμδ

上式右端取第一项，就得到二阶精度的中心差分近似

()()

1

1

12D E E h

h -?-μδ ，即 ()()()2f x h f x h f x h

+--￠? 而取前两项，就能得到四阶精度的中心差分近似

()

3212

11188612D E E E E h h

--骣

÷?÷?=--+-?÷÷?桫

μδ

μδ 即

()()()()()

288212f x h f x h f x h f x h f x h

+-++---￠?

随着精度的提高，差分近似涉及到了更多的点。

但是根据算符验算，又有

352352352

352111

6301

1111166306111116306

11 616

0361

1D h I h I h

I h

I 骣÷?÷=-+-?÷?÷桫骣骣鼢珑鼢=+-+-珑鼢珑鼢桫桫

+é骣÷?ê÷=-+-?÷ê?÷桫?

+ù骣÷?ú÷+-+?÷ú?÷è?=-+×?L L L L μδμδμδδμδμδμδδ

μδμδμδδ

μδμδμδδ

531180

骣÷?÷+-?÷?÷桫L μμδδ

和前一个例子一样，上式中只取第一项，就能得到四阶精度的中心差分近似

2

116

D h

I ?

+μδ

δ

而且该差分近似只涉及三个点。

以上的讨论表明，有理函数近似可以达到我们原来的目的，即：

有理函数近似具有更高的精度，又不涉及更多的点。

下面考虑微分算子有理函数近似在数值格式中的应用。这种有理函数的表达式只是一种算符操作，在实际应用中就需要将有理分式化为整式，具体的方法如下。

例4：由

1

12

D h

I ?

?

-

有

11

2I D h

骣÷?÷-鸦 ?÷?÷

桫 作用在函数 ()

f x 上，

()()11

2I Df x f x h

骣÷?÷-鸦 ?÷?÷

桫 即

()()11

2I f x f x h

骣÷?￠÷-鸦 ?÷?÷

桫 将算子展开，就是

()()

()()

2

f x f x h f x f x h h

ⅱ+---?

对中心差分近似也有类似地的结果。 例5：由

21

16

D h

I ?

+

μδ

δ

有

2116

I D h 骣÷?÷+ ?÷÷?桫δμδ 作用在函数 ()

f x 上，

()()2116

I Df x f x h 骣÷?÷+ ?÷?÷桫δμδ 即

()()2116I f x f x h 骣÷?￠÷+ ?÷?÷

桫δμδ 将算子展开，就是

()()()

()()

46

2f x h f x f x h f x h f x h h

ⅱ +++-+--?

以上两个例子表明，有理函数给出的差分近似，会同时有多个点处的导数值出现，需联立求解。而通常的差分近似，只出现一个点处的导数值，可逐点计算。这两者之间的区别，类似于隐式格式与显式格式的区别。正因为如此，微分算子的有理函数近似也称为隐式差分近似。同时，由于涉及较少的点，通常又称为紧致差分近似。

例6：将例4中的紧致差分近似应用于例1中给出的初值问题，

1102j j j j j j u au u u u u h --í￠?+=????ìⅱ+-??=????

（1,2,3,,j M =L ） 整理后，得到未知解的近似 j u 及其导数值近似 j u ￠ 的联立方程组

11022j j j j j j u au hu u hu u --í￠?+=??ì

?ⅱ?-+=+???

解得

11112 222j j j j j j u hu u ha au hau u ha ----í￠?+??=?+??ì

?￠+??￠=-??+??

（1,2,3,,j M =L ） 对于 0j = ，利用原方程可给出初值

0u =α ， 0

0u au a ￠=-=-α 由此可见，在紧致差分格式的求解过程中，未知解的近似及其导数值的近似都是未知量，是需要联立在一起求解的。

上面的例子是一个两点紧致格式，最终得到了一个递推关系式，逐点计算。对于涉及三个甚至更多点的高阶紧致格式，就需要将未知

解的近似 1u 、2u 、L 、M u 及其导数值的近似 1

u ￠ 、2u ￠ 、L 、M u ￠ （如果原方程还包括二阶导数，则还有二阶导数值的近似 1u ⅱ 、2u ⅱ 、L 、M u ⅱ）全部放在一起联立求解。因此，高阶紧致格式中需

要求解的未知量比较多，这是它的一个弱点。

下面列出一阶导数和二阶导数高阶紧致差分近似的一些结果。 1. Pade 逼近（三点四阶）

11

11

11

11

2

46

210212

j j j j j j j j j j j u u u u u h

u u u u u u h +-+-+-+-ⅱ ++-=ⅱⅱⅱ++-+=

2. 对称紧致格式（五点六阶）

1122

11

11141339294j j j j j j j u u u u u u u h h

+-+-+---ⅱ ++=+ 3. 对称紧致格式（五点八阶）

211

211221441

3699364025272544j j j j j j j j j u u u u u u u u u h h

++--+-+-ⅱⅱ ++++--=+

4. 迎风紧致格式（三点三阶）

1114521

336j j j j j u u u u u h

+--+-ⅱ+=

5. 迎风紧致格式（五点五阶）

21121123644332

5560j j j j j j j u u u u u u u h

++----++--ⅱ+=

6. 广义紧致格式（对称三点六阶）

上面给出的紧致差分近似，计算一阶导数的紧致差分里不会出现二阶导数的近似值，计算二阶导数的紧致差分里也不会出现一阶导数的近似值。如果突破这个限制，就成为广义紧致差分近似。例如

()()()()()()21111121111114882424

1 3948904871673901515

j j j j j j j j j j j j j j j h h u u u u u u u u h h

u u u u u u u ++-+-+-+-+-ⅱⅱⅱ --+-+-+=ⅱⅱⅱ --+++-= 7. 广义紧致格式（迎风两点三阶）

11

12

1

1

1

2

01220

12

2j j j j j j

j j j j j j u u u u u u h h u u u u u u h

h

+++---ⅱⅱⅱ-+--+=ⅱⅱⅱ-+--

+

=

10.2 紧致格式

例7：考虑Burgers 方程（对流扩散方程）两点边值问题

2210 , 01 , x x u u u a x t x x u u u U x ==í?抖 ?+=<

?镦 ??÷?÷?==?÷??÷?镨 ??

νβα 将空间区域 0,1轾犏

臌

均匀划分成M 个网格，则空间网格的尺寸为 1

h M

=

，网格点坐标为 j x jh =（0,1,2,,j M =L ）。 在 n t n t =D 时刻（1,2,3,n =L ，t D 是时间步长），将未

知函数及其空间导数在网格点上的近似值分别记作

()

,n j j u u x t ? ， n

j

t t j x x u

F x

==??

? ， 2

2

n

j

t t j x x u S x ==??

?

现假设上一时刻 ()

11n t n t -=-D 的近似解已经求出，记成

()

1,n j j b u x t -? ，在计算过程中视为已知。

于是，在空间区域 0,1轾犏

臌

内的第j 个网格点（1,2,,1j M =-L ）处，有原方程的差分近似

0j j j j u b aF S t

-+-=D ν

和空间导数的Pade 逼近

11

11

11

11

2

40

26

2100

12

j j j j j j j j j j j u u F F F h

u u u S S S h

+-+-+-+--++-

=-+++-=

在左边界 0j = 处，有边界条件

00u F +=βα

原方程的差分近似

00

000u b aF S t

-+-=D ν 以及空间导数的广义紧致格式

10

101020212

u u F F S S h h -+--+= 在右边界 j M = 处，有边界条件

M u U =

在此处，原方程成为

0M M aF S -=ν

还有空间导数的广义紧致格式

1

11

20212

M M M M M M u u F F S S h h ----+--+=

将所有这些集成在一起，就得到线性方程组

0000002200011111112221111211 021221240333520126120212j j j j j j j j j j j j j j j M M M M M u F u a tF tS b t h h h h u F S u F S u a tF tS b t h h h u F F u F h h h u S u S u S u U a tF tS h h u F --++--++--+=+D -D =D ---+-+=+D -D =D ---+-=---+-==D -D =---M M βαννν210

212M M M M h h S u F S -í??????????????????????????????????????ì????????????????????????????????????+-+=???