对流扩散方程的六阶指数型混合紧致差分格式

Sixth-Order Compact Exponential Combined Difference Schemes for the Convection-Diffusion Equation

摘要

基于泰勒级数展开，采用待定系数法和四阶Padé公式，构造了两种求解一维定常对流扩散方程的六阶指数型混合紧致差分格式。详细讨论了两种差分格式的性质，并比较了两种差分格式与文献中已有格式的色散误差和耗散误差。通过选取典型算例，验证了文中两种六阶指数型混合紧致差分格式的高精度和非振荡特性。新发展的差分格式特别适用于对流占优问题的求解。

Abstract

Based on the Taylor series expansions， two sixth-order compact exponential combined difference schemes are constructed for solving the one-dimensional steady convection-diffusion equation using the method of undetermined coefficients and the fourth-order Padé formula. This paper discusses the properties of the two proposed difference schemes and compares their dispersion errors and dissipation errors with those of existing schemes in the literature. By selecting typical numeral examples， the high precision and non-oscillatory characteristics of the two sixth-order compact exponential combined difference schemes are verified. These newly developed difference schemes are particularly suitable for solving convection-dominated problems.

Graphical abstract

关键词

对流扩散方程 / 紧致差分格式 / 高精度 / 非振荡特性

Key words

convection-diffusion equation / compact difference scheme / high precision / nonoscillation property

对流扩散方程可以被看作Navier-Stokes方程组的简化模型，在计算流体动力学中起重要作用^［1］。因此，研究对流扩散方程精确、稳定和高效的数值计算方法具有重要的理论意义和应用价值。本文考虑如下定常对流扩散方程

（1）

其中：a是扩散系数；c是对流系数；f(x)是x的足够光滑的函数；u(x)是待求未知量。

求解对流扩散方程的数值方法有很多种，如有限元法、有限体积法、有限差分法、特征线法、谱方法和神经网络等^［2-17］。近年来，研究者们已经发展了许多高精度的紧致差分格式^［3-16］。这些格式依据构成差分格式系数函数类型的不同，可分为两类。一类是多项式型差分格式。其构造一般是基于泰勒级数展开，利用差商逼近微分方程中的微商得到的差分格式。例如，文献［5］基于泰勒级数展开和截断误差修正，推导了对流扩散方程的四阶多项式紧致差分格式。文献［6］采用外推法得到了对流扩散方程的六阶多项式型紧致差分格式。文献［14］采用截断误差修正，建立了椭圆型方程两点边值问题的混合型紧致差分格式。然而，对于大梯度问题的求解，多项式型差分格式往往难以达到理论上的计算精度和收敛阶。另一类是指数型紧致差分格式。由于其系数包含相应的微分算子和网格尺度的指数函数，固有地考虑了迎风效应，因此非常适合解决边界层或大梯度的奇异摄动问题。例如，文献［3］基于截断误差余项修正的思想，构造了求解对流扩散方程的四阶指数型紧致差分格式，对于对流占优问题，该格式在粗网格上能够获得很高的计算精度。文献［4］也发展了四阶指数型有限差分格式，以求解定常对流扩散方程，但文献［3］的格式较文献［4］具有更高的计算稳定性和精度。

依据格式结构的不同，有限差分格式可以分为3大类。第一类是方程型紧致差分格式，这类格式仅包含未知量在离散点处的函数值。例如，文献［3］~［11］中提到的相关格式。第二类是导数型紧致差分格式，这类格式包含未知量及其某阶导数在离散点处的函数值。例如，文献［12］提出的组合紧致差分格式（combined compact difference，CCD），仅使用3个网格点就能达到六阶精度。然而，该格式在计算时需要计算每个节点处的待求量及其一阶、二阶导数值，这导致计算量较大，从而影响计算效率。第三类是混合型紧致差分格式，如文献［13］~［14］中所介绍的格式。与方程型格式相比，混合型紧致差分格式具有更高的计算精度；与导数型格式相比，它具有较小的计算量和更高的计算效率。然而，该类格式通常是基于多项式构造的，对于对流占优问题的计算效果不理想。

本文针对模型方程（1），采用待定系数法和截断误差余项修正法，并结合一阶导数和二阶导数的四阶Padé格式，推导出两种求解对流扩散方程的六阶指数型混合紧致差分格式。该格式具有较高的分辨率，能够有效用于对流占优问题的求解。在计算过程中，待求未知量的导数值采用分部迭代计算，避免了求解大型代数方程组，计算效率更高。

1 对流扩散方程的六阶指数型混合紧致差分格式

将求解区间[X_{\mathrm{1}},X_{\mathrm{2}}]等分为N个子区间[x_i,x_{i+\mathrm{1}}]，x_i=X_{\mathrm{1}}+ih,h=\frac{X_{\mathrm{2}}-X_{\mathrm{1}}}{N},i=\mathrm{0,1},\mathrm{2},\cdots ,N,在x_i点处由泰勒级数展开，得到

\begin{array}{l}{u}_{xi}={\delta }_x{u}_i-\frac{h^{\mathrm{2}}}{\mathrm{6}}{\left(\frac{{\partial }^{\mathrm{3}}u}{\partial x^{\mathrm{3}}}\right)}_i-\frac{h^{\mathrm{4}}}{\mathrm{120}}{\left(\frac{{\partial }^{\mathrm{5}}u}{\partial x^{\mathrm{5}}}\right)}_i+O(h^{\mathrm{6}}),\\ u_{xxi}={\delta }_x^{\mathrm{2}}{u}_i-\frac{h^{\mathrm{2}}}{\mathrm{12}}{\left(\frac{{\partial }^{\mathrm{4}}u}{\partial x^{\mathrm{4}}}\right)}_i-\frac{h^{\mathrm{4}}}{\mathrm{360}}{\left(\frac{{\partial }^{\mathrm{6}}u}{\partial x^{\mathrm{6}}}\right)}_i+O(h^{\mathrm{6}}),\end{array}

其中：

{\delta }_x{u}_i=\frac{u_{i+\mathrm{1}}-u_{i-\mathrm{1}}}{\mathrm{2}h}{,}_{}{\delta }_x^{\mathrm{2}}{u}_i=\frac{u_{i+\mathrm{1}}-\mathrm{2}{u}_i+u_{i-\mathrm{1}}}{h^{\mathrm{2}}}。（2）

假设模型方程（1）有如下形式的差分方程，即

-\alpha {\delta }_x^{\mathrm{2}}{u}_i+\beta {\delta }_x{u}_i=f_i{}_{},（3）

对函数\mathrm{1},x,\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{c}{a}x\right)精确成立，联立方程组得到

\alpha =\left\{\begin{array}{l}\frac{ch}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{ch}{\mathrm{2}a}\right){,}_{}c\ne \mathrm{0}，\beta =c，\\ a{,}_{}c=\mathrm{0}，\beta =c\mathrm{。}\end{array}\right.（4）

通过截断误差分析可知，该差分格式的精度为O(h^{\mathrm{2}})。

在此基础上，假设模型方程（1）的差分格式为

\begin{array}{l}-\alpha {\delta }_x^{\mathrm{2}}{u}_i+c{\delta }_x{u}_i=\\ f_i+c_{\mathrm{1}}{f}_{xi}+c_{\mathrm{2}}{f}_{xxi}+c_{\mathrm{3}}{f}_{xxxi}+c_{\mathrm{4}}{f}_{xxxxi}\mathrm{。}\end{array}（5）

为确定系数c_{\mathrm{1}}{,}_{}{c}_{\mathrm{2}}{,}_{}{c}_{\mathrm{3}}{,}_{}{c}_{\mathrm{4}}，利用原模型方程（1），可将式（5）改写为

\begin{array}{l}-\alpha {\delta }_x^{\mathrm{2}}{u}_i+c{\delta }_x{u}_i=(-a{u}_{xx}+c{u}_x{)}_i+\\ c_{\mathrm{1}}(-a{u}_{xx}+c{u}_x{)}_{xi}+c_{\mathrm{2}}(-a{u}_{xx}+c{u}_x{)}_{xxi}+\\ c_{\mathrm{3}}(-a{u}_{xx}+c{u}_x{)}_{xxxi}+c_{\mathrm{4}}(-a{u}_{xx}+c{u}_x{)}_{xxxxi}，\end{array}（6）

直接计算式（6）右端项的导数，并将式（2）中各式代入，得

\begin{array}{l}{\left(-a{u}_{xx}+c{u}_x\right)}_i=f_i+(\alpha -a+c{c}_{\mathrm{1}})u_{x^{\mathrm{2}}i}+\\ \left(-a{c}_{\mathrm{1}}+c_{\mathrm{2}}c-\frac{c{h}^{\mathrm{2}}}{\mathrm{6}}\right)u_{x^{\mathrm{3}}i}+\left(-a{c}_{\mathrm{2}}+c_{\mathrm{3}}c+\frac{\alpha h^{\mathrm{2}}}{\mathrm{12}}\right)u_{x^{\mathrm{4}}i}+\\ \left(-a{c}_{\mathrm{3}}+c_{\mathrm{4}}c-\frac{c{h}^{\mathrm{4}}}{\mathrm{120}}\right)u_{x^{\mathrm{5}}i}+\left(-a{c}_{\mathrm{4}}+\frac{\alpha h^{\mathrm{4}}}{\mathrm{360}}\right)u_{x^{\mathrm{6}}i}+O\left(h^{\mathrm{6}}\right)\mathrm{。}\end{array}（7）

令\left\{\begin{array}{l}\alpha -a+c{c}_{\mathrm{1}}=\mathrm{0}，\\ -a{c}_{\mathrm{1}}+c_{\mathrm{2}}c-\frac{c{h}^{\mathrm{2}}}{\mathrm{6}}=\mathrm{0}，\\ -a{c}_{\mathrm{2}}+c_{\mathrm{3}}c+\frac{\alpha h^{\mathrm{2}}}{\mathrm{12}}=\mathrm{0}，\\ -a{c}_{\mathrm{3}}+c_{\mathrm{4}}c-\frac{c{h}^{\mathrm{4}}}{\mathrm{120}}=\mathrm{0}，\end{array}\right.（8）

解得

\begin{array}{l}{c}_{\mathrm{1}}=\left\{\begin{array}{l}\frac{a-\alpha }{c}{,}_{}c\ne \mathrm{0}，\\ \mathrm{0}{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{2}}=\left\{\begin{array}{l}\frac{a\left(a-\alpha \right)}{c^{\mathrm{2}}}+\frac{h^{\mathrm{2}}}{\mathrm{6}}{,}_{}c\ne \mathrm{0}，\\ \frac{h^{\mathrm{2}}}{\mathrm{12}}{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{3}}=\left\{\begin{array}{l}\frac{a^{\mathrm{2}}\left(a-\alpha \right)}{c^{\mathrm{3}}}+\frac{\left(\mathrm{2}a-\alpha \right)}{\mathrm{12}c}{h}^{\mathrm{2}}{,}_{}c\ne \mathrm{0}，\\ \mathrm{0}{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{4}}=\left\{\begin{array}{l}\frac{a^{\mathrm{3}}\left(a-\alpha \right)}{c^{\mathrm{4}}}+\frac{a\left(\mathrm{2}a-\alpha \right)}{\mathrm{12}{c}^{\mathrm{2}}}{h}^{\mathrm{2}}+\frac{h^{\mathrm{4}}}{\mathrm{120}}{,}_{}c\ne \mathrm{0}，\\ \frac{h^{\mathrm{4}}}{\mathrm{360}}{,}_{}c=\mathrm{0}\mathrm{。}\end{array}\right.\end{array}（9）

接下来，在式（5）的基础上构造模型方程（1）的六阶指数型混合紧致差分格式。利用泰勒级数展开可得

f_{xi}={\delta }_x{f}_i-\frac{h^{\mathrm{2}}}{\mathrm{6}}{\left(\frac{{\partial }^{\mathrm{3}}f}{\partial x^{\mathrm{3}}}\right)}_i-\frac{h^{\mathrm{4}}}{\mathrm{120}}{\left(\frac{{\partial }^{\mathrm{5}}f}{\partial x^{\mathrm{5}}}\right)}_i+O(h^{\mathrm{6}})，（10）

f_{xxi}={\delta }_x^{\mathrm{2}}{f}_i-\frac{h^{\mathrm{2}}}{\mathrm{12}}{\left(\frac{{\partial }^{\mathrm{4}}f}{\partial x^{\mathrm{4}}}\right)}_i-\frac{h^{\mathrm{4}}}{\mathrm{360}}{\left(\frac{{\partial }^{\mathrm{6}}f}{\partial x^{\mathrm{6}}}\right)}_i+O(h^{\mathrm{6}})。（11）

将式（10）中的f_i用f_{xi}代替，得

f_{xxi}={\delta }_x{f}_{xi}-\frac{h^{\mathrm{2}}}{\mathrm{6}}{\left(\frac{{\partial }^{\mathrm{3}}{f}_x}{\partial x^{\mathrm{3}}}\right)}_i-\frac{h^{\mathrm{4}}}{\mathrm{120}}{\left(\frac{{\partial }^{\mathrm{5}}{f}_x}{\partial x^{\mathrm{5}}}\right)}_i+O(h^{\mathrm{6}})。（12）

由式（10）得

f_{xxxi}=\frac{\mathrm{6}}{h^{\mathrm{2}}}\left({\delta }_x{f}_i-f_{xi}\right)+O(h^{\mathrm{2}})，（13）

而由式（11）和式（12）得

f_{xxxxi}=\frac{\mathrm{12}}{h^{\mathrm{2}}}\left(-{\delta }_x^{\mathrm{2}}{f}_i+{\delta }_x{f}_{xi}\right)+O(h^{\mathrm{2}})\mathrm{。}（14）

将式（13）和式（14）代入式（5）右端并忽略高阶项，得

\begin{array}{l}{F}_i=f_i+c_{\mathrm{1}}{f}_{xi}+c_{\mathrm{2}}{f}_{xxi}+c_{\mathrm{3}}{f}_{xxxi}+c_{\mathrm{4}}{f}_{xxxxi}=\\ \left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x^{\mathrm{2}}\right)f_i+\\ \left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x\right)f_{xi}+c_{\mathrm{2}}{f}_{xxi}\mathrm{。}\end{array}（15）

对式（15）中f_{xi}采用四阶Padé算子离散，其边界采用文献［15］中四阶一致边界格式离散，具体公式为

\left\{\begin{array}{l}\frac{\mathrm{1}}{\mathrm{6}}{\left(f_x\right)}_{i-\mathrm{1}}+\frac{\mathrm{2}}{\mathrm{3}}{\left(f_x\right)}_i+\frac{\mathrm{1}}{\mathrm{6}}{\left(f_x\right)}_{i+\mathrm{1}}=\\ {\delta }_x{f}_i,i=\mathrm{1,2},\cdots ,N-\mathrm{1},\\ f_{\mathrm{0}}^{\text{'}}+\frac{\mathrm{8}}{\mathrm{9}}{f}_{\mathrm{1}}^{\text{'}}=\frac{\mathrm{1}}{h}\left(-\frac{\mathrm{221}}{\mathrm{90}}{f}_{\mathrm{0}}+\frac{\mathrm{433}}{\mathrm{108}}{f}_{\mathrm{1}}-\frac{\mathrm{19}}{\mathrm{6}}{f}_{\mathrm{2}}+\right.\\ \left.\frac{\mathrm{43}}{\mathrm{18}}{f}_{\mathrm{3}}-\frac{\mathrm{25}}{\mathrm{27}}{f}_{\mathrm{4}}+\frac{\mathrm{27}}{\mathrm{180}}{f}_{\mathrm{5}}\right),\\ f_N^{\text{'}}-\frac{\mathrm{8}}{\mathrm{9}}{f}_{N-\mathrm{1}}^{\text{'}}=\frac{\mathrm{1}}{h}\left(\frac{\mathrm{21}}{\mathrm{10}}{f}_N-\frac{\mathrm{641}}{\mathrm{108}}{f}_{N-\mathrm{1}}+\frac{\mathrm{121}}{\mathrm{18}}{f}_{N-\mathrm{2}}-\right.\\ \left.\frac{\mathrm{25}}{\mathrm{6}}{f}_{N-\mathrm{3}}+\frac{\mathrm{41}}{\mathrm{27}}{f}_{N-\mathrm{4}}-\frac{\mathrm{43}}{\mathrm{180}}{f}_{N-\mathrm{5}}\right)\mathrm{。}\end{array}\right.

将计算f(x)一阶导数的四阶Padé算子记为

P_x{f}_i{:}_{}{\delta }_x{f}_i=\frac{\mathrm{1}}{\mathrm{6}}{\left(f_x\right)}_{i-\mathrm{1}}+\frac{\mathrm{2}}{\mathrm{3}}{\left(f_x\right)}_i+\frac{\mathrm{1}}{\mathrm{6}}{\left(f_x\right)}_{i+\mathrm{1}},

则得到

\begin{array}{l}{F}_i=\left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x^{\mathrm{2}}\right)f_i+\\ \left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x\right)P_x{f}_i+c_{\mathrm{2}}{f}_{xxi}\mathrm{。}\end{array}（16）

下面通过对式（16）中f_{xxi}采用两种不同的离散方式，得到模型方程（1）的两种指数型组合紧致差分格式。

1.1 差分格式Ⅰ

对式（16）右端项中f_{xxi}采用四阶Padé算子离散，其边界采用文献［15］所给的四阶一致边界格式离散，具体公式为

\left\{\begin{array}{l}\frac{\mathrm{1}}{\mathrm{12}}{\left(f_{xx}\right)}_{i-\mathrm{1}}+\frac{\mathrm{5}}{\mathrm{6}}{\left(f_{xx}\right)}_i+\frac{\mathrm{1}}{\mathrm{12}}{\left(f_{xx}\right)}_{i+\mathrm{1}}=\\ {\delta }_x^{\mathrm{2}}{f}_i{,}_{}{}_{}i=\mathrm{1,2},\cdots ,N-\mathrm{1},\\ f_{\mathrm{0}}^{″}+\frac{\mathrm{51}}{\mathrm{52}}{f}_{\mathrm{1}}^{″}=\frac{\mathrm{1}}{h^{\mathrm{2}}}\left(\frac{\mathrm{12}\mathrm{293}}{\mathrm{2}\mathrm{340}}{f}_{\mathrm{0}}-\frac{\mathrm{18}\mathrm{903}}{\mathrm{1}\mathrm{040}}{f}_{\mathrm{1}}+\frac{\mathrm{2}\mathrm{891}}{\mathrm{104}}{f}_{\mathrm{2}}-\right.\\ \left.\frac{\mathrm{23}\mathrm{941}}{\mathrm{936}}{f}_{\mathrm{3}}+\frac{\mathrm{387}}{\mathrm{26}}{f}_{\mathrm{4}}-\frac{\mathrm{5}\mathrm{063}}{\mathrm{1}\mathrm{040}}{f}_{\mathrm{5}}+\frac{\mathrm{494}}{\mathrm{720}}{f}_{\mathrm{6}}\right),\\ f_N^{″}-\frac{\mathrm{51}}{\mathrm{52}}{f}_{N-\mathrm{1}}^{″}=\frac{\mathrm{1}}{h^{\mathrm{2}}}\left(\frac{\mathrm{17}\mathrm{599}}{\mathrm{4}\mathrm{680}}{f}_N-\frac{\mathrm{17}\mathrm{237}}{\mathrm{1}\mathrm{040}}{f}_{N-\mathrm{1}}+\right.\\ \frac{\mathrm{795}}{\mathrm{26}}{f}_{N-\mathrm{2}}-\frac{\mathrm{28}\mathrm{735}}{\mathrm{936}}{f}_{N-\mathrm{3}}+\frac{\mathrm{1}\mathrm{871}}{\mathrm{104}}{f}_{N-\mathrm{4}}-\\ \left.\frac{\mathrm{6}\mathrm{117}}{\mathrm{1}\mathrm{040}}{f}_{N-\mathrm{5}}+\frac{\mathrm{596}}{\mathrm{720}}{f}_{N-\mathrm{6}}\right)\mathrm{。}\end{array}\right.

记f(x)二阶导数的四阶Padé算子为

P_x^{\mathrm{2}}{f}_i{:}_{}{\delta }_x^{\mathrm{2}}{f}_i=\frac{\mathrm{1}}{\mathrm{12}}{\left(f_{xx}\right)}_{i-\mathrm{1}}+\frac{\mathrm{5}}{\mathrm{6}}{\left(f_{xx}\right)}_i+\frac{\mathrm{1}}{\mathrm{12}}{\left(f_{xx}\right)}_{i+\mathrm{1}}{}_{},

则式（16）可写成

\begin{array}{l}{F}_i=\left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x^{\mathrm{2}}\right)f_i+\\ \left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x\right)P_x{f}_i+c_{\mathrm{2}}{P}_x^{\mathrm{2}}{f}_i{}_{}\mathrm{。}\end{array}

因此得到差分格式Ⅰ，即

-\alpha {\delta }_x^{\mathrm{2}}{u}_i+c{\delta }_x{u}_i=F_i\mathrm{。}（17）

其中：

\begin{array}{c}\alpha =\left\{\begin{array}{l}\frac{ch}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{ch}{\mathrm{2}a}\right){,}_{}c\ne \mathrm{0},\\ a{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{1}}=\left\{\begin{array}{l}\frac{a-\alpha }{c}{,}_{}c\ne \mathrm{0},\\ \mathrm{0}{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{2}}=\left\{\begin{array}{l}\frac{a\left(a-\alpha \right)}{c^{\mathrm{2}}}+\frac{h^{\mathrm{2}}}{\mathrm{6}}{,}_{}c\ne \mathrm{0},\\ \frac{h^{\mathrm{2}}}{\mathrm{12}}{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{3}}=\left\{\begin{array}{l}\frac{a^{\mathrm{2}}\left(a-\alpha \right)}{c^{\mathrm{3}}}+\frac{\left(\mathrm{2}a-\alpha \right)}{\mathrm{12}c}{h}^{\mathrm{2}}{,}_{}c\ne \mathrm{0},\\ \mathrm{0}{,}_{}c=\mathrm{0}；\end{array}\right.\\ c_{\mathrm{4}}=\left\{\begin{array}{l}\frac{a^{\mathrm{3}}\left(a-\alpha \right)}{c^{\mathrm{4}}}+\frac{a\left(\mathrm{2}a-\alpha \right)}{\mathrm{12}{c}^{\mathrm{2}}}{h}^{\mathrm{2}}+\frac{h^{\mathrm{4}}}{\mathrm{120}}{,}_{}c\ne \mathrm{0},\\ \frac{h^{\mathrm{4}}}{\mathrm{360}}{,}_{}c=\mathrm{0}\mathrm{。}\end{array}\right.\end{array}

1.2 差分格式Ⅱ

将式（11）左、右两端同乘以2，并与式（12）相减，可得

f_{xxi}=\mathrm{2}{\delta }_x^{\mathrm{2}}{f}_i-{\delta }_x{f}_{xi}+O(h^{\mathrm{4}})\mathrm{。}

所以，对式（16）右端项中f_{xxi}采用上述四阶差分算子离散，可得

\begin{array}{l}{F}_i=\left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x+\left(\mathrm{2}{c}_{\mathrm{2}}-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}\right){\delta }_x^{\mathrm{2}}\right)f_i+\\ \left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x-c_{\mathrm{2}}{\delta }_x\right)P_x{f}_i{}_{},\end{array}

因此得到差分格式Ⅱ，即

-\alpha {\delta }_x^{\mathrm{2}}{u}_i+c{\delta }_x{u}_i=F_i{}_{},（18）

其中\alpha，c_{\mathrm{1}}，c_{\mathrm{2}}，c_{\mathrm{3}}，c_{\mathrm{4}}的表示同差分格式Ⅰ。

1.3 截断误差分析

假设Pe=\frac{ch}{a}，则有

\begin{array}{l}{c}_{\mathrm{1}}=h\left(\frac{\mathrm{1}}{Pe}-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}\right),\\ c_{\mathrm{2}}=h^{\mathrm{2}}\left(\frac{\mathrm{1}}{P{e}^{\mathrm{2}}}-\frac{\mathrm{1}}{\mathrm{2}Pe}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{6}}\right),\\ c_{\mathrm{3}}=h^{\mathrm{3}}\left(\frac{\mathrm{1}}{P{e}^{\mathrm{3}}}-\frac{\mathrm{1}}{\mathrm{2}P{e}^{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{6}Pe}-\frac{\mathrm{1}}{\mathrm{24}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}\right),\\ c_{\mathrm{4}}=h^{\mathrm{4}}\left(\frac{\mathrm{1}}{P{e}^{\mathrm{4}}}-\frac{\mathrm{1}}{\mathrm{2}P{e}^{\mathrm{3}}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{6}P{e}^{\mathrm{2}}}-\frac{\mathrm{1}}{\mathrm{24}Pe}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{120}}\right)\mathrm{。}\end{array}（19）

由双曲余切函数的泰勒级数展开为

\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}(x)=\frac{\mathrm{1}}{x}+\frac{x}{\mathrm{3}}-\frac{x^{\mathrm{3}}}{\mathrm{45}}+\frac{\mathrm{2}{x}^{\mathrm{5}}}{\mathrm{945}}-\frac{x^{\mathrm{7}}}{\mathrm{4}\mathrm{725}}+\cdots +\frac{{(-\mathrm{1})}^{n-\mathrm{1}}{\mathrm{2}}^{\mathrm{2}n}{B}_n{x}^{\mathrm{2}n-\mathrm{1}}}{(\mathrm{2}n)!}+\cdots，

可得

c_{\mathrm{1}}=h^{\mathrm{2}}\left[-\frac{c}{\mathrm{12}a}+\frac{c^{\mathrm{3}}{h}^{\mathrm{2}}}{\mathrm{720}{a}^{\mathrm{3}}}-\frac{c^{\mathrm{5}}{h}^{\mathrm{4}}}{\mathrm{30}\mathrm{240}{a}^{\mathrm{5}}}+\frac{c^{\mathrm{7}}{h}^{\mathrm{6}}}{\mathrm{1}\mathrm{209}\mathrm{600}{a}^{\mathrm{7}}}\right]+O(h^{\mathrm{8}}),

c_{\mathrm{2}}=h^{\mathrm{2}}\left[\frac{\mathrm{1}}{\mathrm{12}}+\frac{c^{\mathrm{2}}{h}^{\mathrm{2}}}{\mathrm{720}{a}^{\mathrm{2}}}-\frac{c^{\mathrm{4}}{h}^{\mathrm{4}}}{\mathrm{30}\mathrm{240}{a}^{\mathrm{4}}}+\frac{c^{\mathrm{6}}{h}^{\mathrm{6}}}{\mathrm{1}\mathrm{209}\mathrm{600}{a}^{\mathrm{6}}}\right]+O(h^{\mathrm{8}}),

c_{\mathrm{3}}=h^{\mathrm{4}}\left[-\frac{c}{\mathrm{180}a}+\frac{c^{\mathrm{3}}{h}^{\mathrm{2}}}{\mathrm{12}\mathrm{096}{a}^{\mathrm{3}}}-\frac{c^{\mathrm{5}}{h}^{\mathrm{4}}}{\mathrm{518}\mathrm{400}{a}^{\mathrm{5}}}+\frac{c^{\mathrm{7}}{h}^{\mathrm{6}}}{\mathrm{14}\mathrm{515}\mathrm{200}{a}^{\mathrm{7}}}\right]+O(h^{\mathrm{8}}),

c_{\mathrm{4}}=h^{\mathrm{4}}\left[\frac{\mathrm{1}}{\mathrm{360}}+\frac{c^{\mathrm{2}}{h}^{\mathrm{2}}}{\mathrm{12}\mathrm{096}{a}^{\mathrm{2}}}-\frac{c^{\mathrm{4}}{h}^{\mathrm{4}}}{\mathrm{518}\mathrm{400}{a}^{\mathrm{4}}}+\frac{c^{\mathrm{6}}{h}^{\mathrm{6}}}{\mathrm{14}\mathrm{515}\mathrm{200}{a}^{\mathrm{6}}}\right]+O(h^{\mathrm{8}})\mathrm{。}

此时

记式（16）右端项为

F_i=F_{\mathrm{1}i}+F_{\mathrm{2}i},（20）

其中F_{\mathrm{1}i}=\left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x^{\mathrm{2}}\right)f_i+\left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x\right)P_x{f}_i，F_{\mathrm{2}i}=c_{\mathrm{2}}{f}_{xxi}。对于差分格式Ⅰ和差分格式Ⅱ，式（20）中的第一部分F_{\mathrm{1}i}的截断误差是相同的，第二部分F_{\mathrm{2}i}的截断误差依据f_{xxi}的离散方式不同而不同，记其截断误差为R_i=R_{\mathrm{1}i}+R_{\mathrm{2}i}。

式（16）是对式（5）右端项中f_{xi}采用四阶Padé算子离散，对f_{xxxi}和f_{xxxxi}分别采用式（13）和式（14）离散所得到的只含f_{xi}和f_{xxi}的表达式，因此式（20）中第一部分F_{\mathrm{1}i}的截断误差等价于c_{\mathrm{1}}{f}_{xi}+c_{\mathrm{3}}{f}_{xxxi}+c_{\mathrm{4}}{f}_{xxxxi}的截断误差。然而，当f_{xi}采用四阶Padé算子离散时，其截断误差为

R_{f_{xi}}=\frac{h^{\mathrm{4}}}{\mathrm{180}}{f}_{x^{\mathrm{5}}i}+O\left(h^{\mathrm{5}}\right),

当f_{xxxi}采用式（13）离散时，其截断误差为

R_{f_{xxxi}}=-\frac{h^{\mathrm{2}}}{\mathrm{20}}{f}_{x^{\mathrm{5}}i}+O(h^{\mathrm{4}}),

当f_{xxxxi}采用式（14）离散时，其截断误差为

R_{f_{xxxxi}}=\frac{h^{\mathrm{2}}}{\mathrm{15}}{f}_{x^{\mathrm{6}}i}+O(h^{\mathrm{4}}),

所以F_{\mathrm{1}i}的截断误差为

R_{\mathrm{1}i}=-\frac{h^{\mathrm{6}}}{\mathrm{5}\mathrm{400}}\left(\frac{c}{a}{f}_{x^{\mathrm{5}}i}-f_{x^{\mathrm{6}}i}\right)+O(h^{\mathrm{7}})\mathrm{。}

对于差分格式Ⅰ，f_{xxi}的截断误差为

R_{\mathrm{1}{f}_{xxi}}=\frac{h^{\mathrm{4}}}{\mathrm{240}}{f}_{x^{\mathrm{6}}i}+O(h^{\mathrm{5}}),

所以F_{\mathrm{2}i}的截断误差为

R_{\mathrm{2}i}=\frac{h^{\mathrm{6}}}{\mathrm{2}\mathrm{880}}{f}_{x^{\mathrm{6}}i}+O(h^{\mathrm{7}})\mathrm{。}

对于差分格式Ⅱ，f_{xxi}的截断误差为

R_{\mathrm{2}{f}_{xxi}}=\frac{h^{\mathrm{4}}}{\mathrm{360}}{f}_{x^{\mathrm{6}}i}+O(h^{\mathrm{6}}),

所以F_{\mathrm{2}i}的截断误差为

R_{\mathrm{2}i}=\frac{h^{\mathrm{6}}}{\mathrm{4}\mathrm{320}}{f}_{x^{\mathrm{6}}i}+O(h^{\mathrm{8}})\mathrm{。}

综上所述，本文两种差分格式的截断误差都是O\left(h^{\mathrm{6}}\right)，即本文构造的两种差分格式在时，都是六阶指数型混合紧致差分格式。

1.4 算法描述

将差分格式Ⅰ和差分格式Ⅱ左端的差分算子展开后，两种格式可以统一为

\left(-\frac{\alpha }{h^{\mathrm{2}}}-\frac{c}{\mathrm{2}h}\right)u_{i-\mathrm{1}}+\left(\frac{\mathrm{2}\alpha }{h^{\mathrm{2}}}\right)u_i+\left(-\frac{\alpha }{h^{\mathrm{2}}}+\frac{c}{\mathrm{2}h}\right)u_{i+\mathrm{1}}=F_i\mathrm{。}（21）

对于差分格式Ⅰ，有

\begin{array}{l}{F}_i=\left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x^{\mathrm{2}}\right)f_i+\\ \left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x\right)P_x{f}_i+c_{\mathrm{2}}{P}_x^{\mathrm{2}}{f}_i\end{array}。

对于差分格式Ⅱ，有

\begin{array}{l}{F}_i=\left(\mathrm{1}+\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}{\delta }_x+\left(\mathrm{2}{c}_{\mathrm{2}}-\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}\right){\delta }_x^{\mathrm{2}}\right)f_i+\\ \left(c_{\mathrm{1}}-\frac{\mathrm{6}{c}_{\mathrm{3}}}{h^{\mathrm{2}}}+\frac{\mathrm{12}{c}_{\mathrm{4}}}{h^{\mathrm{2}}}{\delta }_x-c_{\mathrm{2}}{\delta }_x\right)P_x{f}_i\end{array}。

其算法流程图见图1。

2 差分格式的性质

不失一般性，令a=\mathrm{1}，则模型方程可写成

-u_{xx}+c{u}_x=f(x)\mathrm{。}（22）

2.1 人工扩散效应分析

人工扩散技术是在差分方程（显式或隐式）中引入数值扩散项，用于近似微分方程，从而获得对流占优问题的无振荡光滑解。因此，它是高效求解对流占优问题的一种常用且有效的手段。然而，过量添加人工扩散，在处理大梯度变换时会降低格式的计算精度，进而影响计算效率。因此，对于奇异摄动问题，必须谨慎添加人工扩散，其添加量最好小于传统迎风格式中的人工扩散值ch/\mathrm{2}。例如，FOC格式、迎风格式和本文格式的人工扩散项分别表述如下。

● FOC格式^［5］为

-\left(\mathrm{1}+\frac{P{e}^{\mathrm{2}}}{\mathrm{3}}\right){\delta }_x^{\mathrm{2}}{u}_i+c{\delta }_x{u}_i=f_i+\frac{h^{\mathrm{2}}}{\mathrm{12}}{\left(f_{xx}-c{f}_x\right)}_i,

其人工扩散项为\frac{P{e}^{\mathrm{2}}}{\mathrm{3}}；

● 迎风格式的人工扩散项为\frac{Pe}{\mathrm{2}}；

● 本文差分格式的人工扩散项为\frac{Pe}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{Pe}{\mathrm{2}}\right)-\mathrm{1}。

当Pe>\mathrm{0}时，\frac{Pe}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{Pe}{\mathrm{2}}\right)-\mathrm{1}单调递增。因此当Pe=\frac{\mathrm{3}}{\mathrm{2}}时，有

当时，，而

所以有

。

当Pe>\frac{\mathrm{3}}{\mathrm{2}}时，，而

\begin{array}{l}\frac{Pe}{\mathrm{2}}-\left[\frac{Pe}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{Pe}{\mathrm{2}}\right)-\mathrm{1}\right]=\\ \mathrm{1}-\frac{Pe}{{\mathrm{e}}^{Pe}-\mathrm{1}}>\mathrm{1}-\frac{\mathrm{1.5}}{{\mathrm{e}}^{\mathrm{1.5}}-\mathrm{1}}>\mathrm{0}，\end{array}

因此

。

综上所述，当Pe\in (\mathrm{0},+\mathrm{\infty })时，本文差分格式的人工扩散项恒小于其他两种格式的人工扩散项（图2）。

2.2 非振荡行为

对模型方程（22），本文两种差分格式的系数矩阵为

\begin{array}{l}\mathit{A}=\\ \mathrm{t}\mathrm{r}\mathrm{i}[-\gamma (\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\gamma +\mathrm{1}),\mathrm{2}\gamma \mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\gamma ,-\gamma (\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\gamma -\mathrm{1})],\end{array}

其中\gamma =\frac{ch}{\mathrm{2}}=\frac{Pe}{\mathrm{2}}。

引理1^［18］ 若矩阵\mathit{A}严格对角占优或不可约对角占优，则\mathit{A}非奇异。

引理2^［18］ 设\mathit{A}=(a_{ij})为实n\times n矩阵，若矩阵\mathit{A}对角占优或不可约对角占优，且对i\ne j都有a_{ij}\le \mathrm{0}，对所有\mathrm{1}\le i\le n都有a_{ii}>\mathrm{0}，则{\mathit{A}}^{-\mathrm{1}}\ge \mathrm{0}。

引理3^［18］n\times n的复矩阵\mathit{A}为不可约的充要条件是它的方向图G(\mathit{A})强连接。

引理4^［18］ 设\mathit{A}=\left(a_{ij}\right)为实n\times n矩阵，对所有的i\ne j都有a_{ij}\le \mathrm{0}，则当且仅当\mathit{A}是\mathrm{M}矩阵时，\mathit{A}是单调的。

定理1 差分格式（17）和（18）的系数矩阵\mathit{A}是\mathrm{M}矩阵且{\mathit{A}}^{-\mathrm{1}}\ge \mathrm{0}。

证明 （ⅰ） 证明系数矩阵\mathit{A}={\left(a_{ij}\right)}_{n\times n}满足a_{ii}>\mathrm{0}{,}_{}{a}_{ij}\le \mathrm{0}{,}_{}i\ne j,i,j=\mathrm{1,2},\cdots ，n。

对\forall \gamma \ne \mathrm{0}，有\mathrm{2}\gamma \mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}(\gamma )>\mathrm{0}。这是因为当 \gamma >\mathrm{0}时，\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}(\gamma )>\mathrm{1}\Rightarrow \mathrm{2}\gamma \mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}(\gamma )>\mathrm{0}；当时，。

对\forall \gamma \ne \mathrm{0}，有

这是因为当\gamma >\mathrm{0}时，当时，，即对任意\gamma \ne \mathrm{0}，a_{ii}>\mathrm{0}，a_{ij}\le \mathrm{0},i\ne j,i,j=\mathrm{1,2},\cdots ，n。

（ⅱ） 证明系数矩阵\mathit{A}是不可约的。

系数矩阵为

\begin{array}{l}\mathit{A}=\\ \mathrm{t}\mathrm{r}\mathrm{i}[-\gamma (\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\gamma +\mathrm{1}),\mathrm{2}\gamma \mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\gamma ,-\gamma (\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\gamma -\mathrm{1})],\end{array}

其方向图见图3。由于方向图G(\mathit{A})是强连接的，系数矩阵\mathit{A}不可约。

（ⅲ） 证明系数矩阵\mathit{A}对角占优。

当i=\mathrm{2,3},\cdots ,n-\mathrm{1}时，

{\displaystyle \sum _{\begin{array}{l}j=\mathrm{1}\\ j\ne i\end{array}}^n}\left|a_{nj}\right|=\left|-\gamma \left(\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\gamma \right)+\mathrm{1}\right)\right|+\left|-\gamma \left(\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\gamma \right)-\mathrm{1}\right)\right|=\gamma \left(\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\gamma \right)+\mathrm{1}\right)+\gamma \left(\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\gamma \right)-\mathrm{1}\right)=\mathrm{2}\gamma \mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\gamma \right)=\left|a_{ii}\right|,

当i=\mathrm{1}时，有

这是因为当\gamma >\mathrm{0}时，有

当时，有

同理，当i=n时，有

即 \left|a_{ii}\right|\ge {\displaystyle \sum _{\begin{array}{l}j=\mathrm{1}\\ j\ne i\end{array}}^n}\left|a_{ij}\right|。

因此，由（ⅱ）、（ⅲ）及引理1可知系数矩阵\mathit{A}={\left(a_{ij}\right)}_{n\times n}非奇异；由（ⅰ）~（ⅲ）及引理1得{\mathit{A}}^{-\mathrm{1}}\ge \mathrm{0}，故系数矩阵\mathit{A}={\left(a_{ij}\right)}_{n\times n}是\mathrm{M}矩阵。

2.3 收敛性

不妨设\left\{u(x)\mathrm{0}\left|\le x\le L\right.\right\}为定解问题

的解。定义网格函数U=\left\{U_i\left|\mathrm{0}\le i\le m\right.\right\}，其中U_i=u(x_i){,}_{}\mathrm{0}\le i\le m。在网格点上考虑上述定解问题，有

\left\{\begin{array}{l}-u_{xx}(x_i)+c{u}_x(x_i)=f(x_i){,}_{}\mathrm{1}\le i\le m-\mathrm{1},\\ u(\mathrm{0})={\beta }_{\mathrm{1}}{,}_{}u(L)={\beta }_{\mathrm{2}},\end{array}\right.（23）

则式（23）满足差分格式

\left\{\begin{array}{l}-\alpha {\delta }_x^{\mathrm{2}}{U}_i+c{\delta }_x{U}_i=F_i+R_i,\\ U_{\mathrm{0}}={\beta }_{\mathrm{1}}{,}_{}{U}_m={\beta }_{\mathrm{2}},\end{array}\right.

其中\alpha =\left\{\begin{array}{l}\frac{ch}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{ch}{\mathrm{2}}\right){,}_{}c\ne \mathrm{0},\\ \mathrm{1}{,}_{}c=\mathrm{0}\mathrm{。}\end{array}\right.

设\left\{u_i\left|\mathrm{0}\le i\le m\right.\right\}为本文差分格式的解，则

\left\{\begin{array}{l}-\alpha {\delta }_x^{\mathrm{2}}{u}_i+c{\delta }_x{u}_i=F_i,\\ u_{\mathrm{0}}={\beta }_{\mathrm{1}}{,}_{}{u}_m={\beta }_{\mathrm{2}}\mathrm{。}\end{array}\right.

记e_i=u_i-u(x_i){,}_{}\mathrm{0}\le i\le m, 则得到误差方程组

\left\{\begin{array}{l}-\alpha {\delta }_x^{\mathrm{2}}{e}_i+c{\delta }_x{e}_i=R_i{,}_{}\mathrm{1}\le i\le m-\mathrm{1},\\ e_{\mathrm{0}}=\mathrm{0}{,}_{}{e}_m=\mathrm{0},\end{array}\right.

其中\alpha =\left\{\begin{array}{l}\frac{ch}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\left(\frac{ch}{\mathrm{2}}\right){,}_{}c\ne \mathrm{0},\\ \mathrm{1}{,}_{}c=\mathrm{0}\mathrm{。}\end{array}\right.

由1.3节截断误差分析可知，对于差分格式Ⅰ，有

R_i=\left(-\frac{c}{\mathrm{5}\mathrm{400}}{f}_{x^{\mathrm{5}}i}+\frac{\mathrm{23}}{\mathrm{43}\mathrm{200}}{f}_{x^{\mathrm{6}}i}\right)h^{\mathrm{6}},

此时可得

\begin{array}{l}{‖e‖}_{\mathrm{\infty }}\le \frac{L^{\mathrm{2}}}{\mathrm{2}\sqrt[]{\mathrm{6}}}\underset{\mathrm{1}\le i\le m-\mathrm{1}}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|R_i\right|\le \\ \frac{L^{\mathrm{2}}{h}^{\mathrm{6}}}{\mathrm{10}\mathrm{800}\sqrt[]{\mathrm{6}}}\left(-c{M}_{\mathrm{5}}+\frac{\mathrm{23}}{\mathrm{8}}{M}_{\mathrm{6}}\right);\end{array}

对于差分格式Ⅱ，有

R_i=\left(-\frac{c}{\mathrm{5}\mathrm{400}}{f}_{x^{\mathrm{5}}i}+\frac{\mathrm{1}}{\mathrm{2}\mathrm{400}}{f}_{x^{\mathrm{6}}i}\right)h^{\mathrm{6}},

则

\begin{array}{l}{‖e‖}_{\mathrm{\infty }}\le \frac{L^{\mathrm{2}}}{\mathrm{2}\sqrt[]{\mathrm{6}}}\underset{\mathrm{1}\le i\le m-\mathrm{1}}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|R_i\right|\le \\ \frac{L^{\mathrm{2}}{h}^{\mathrm{6}}}{\mathrm{1}\mathrm{200}\sqrt[]{\mathrm{6}}}\left(-\frac{c}{\mathrm{9}}{M}_{\mathrm{5}}+\frac{\mathrm{1}}{\mathrm{4}}{M}_{\mathrm{6}}\right);\end{array}

其中

M_{\mathrm{5}}=\underset{\mathrm{0}\le x\le L}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|f_{x^{\mathrm{5}}}(x)\right|{,}_{}{M}_{\mathrm{6}}=\underset{\mathrm{0}\le x\le L}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|f_{x^{\mathrm{6}}}(x)\right|。

3 色散误差和耗散误差分析

假设\tilde{u}={\mathrm{e}}^{\mathrm{i}kx}，i为虚数单位，则{\delta }_x\tilde{u}=\mathrm{i}\frac{\mathrm{s}\mathrm{i}\mathrm{n}kh}{h}\tilde{u},

{\delta }_x^{\mathrm{2}}\tilde{u}=\frac{\mathrm{2}}{h^{\mathrm{2}}}(\mathrm{c}\mathrm{o}\mathrm{s}\left(kh\right)-\mathrm{1})\tilde{u}。令波数\omega =kh，Pe=ch/a，

则微分方程（1）的精确特征函数值为

\lambda =\frac{a}{h^{\mathrm{2}}}\left({\omega }^{\mathrm{2}}+\mathrm{i}\omega Pe\right)。

又令\alpha =a\gamma ,\mathrm{其中}\gamma =a\frac{Pe}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{h}\frac{Pe}{\mathrm{2}}，则

c_{\mathrm{1}}={\sigma }_{\mathrm{1}}h{,}_{}{c}_{\mathrm{2}}={\sigma }_{\mathrm{2}}{h}^{\mathrm{2}}{,}_{}{c}_{\mathrm{3}}={\sigma }_{\mathrm{3}}{h}^{\mathrm{3}}{,}_{}{c}_{\mathrm{4}}={\sigma }_{\mathrm{4}}{h}^{\mathrm{4}}，

其中：

{\sigma }_{\mathrm{1}}=\frac{\mathrm{1}}{Pe}\left(\mathrm{1}-\gamma \right)，{\sigma }_{\mathrm{2}}=\frac{\mathrm{1}}{P{e}^{\mathrm{2}}}\left(\mathrm{1}-\gamma +\frac{P{e}^{\mathrm{2}}}{\mathrm{6}}\right)，

{\sigma }_{\mathrm{3}}=\frac{\mathrm{1}}{P{e}^{\mathrm{3}}}\left(\mathrm{1}-\gamma +\frac{P{e}^{\mathrm{2}}}{\mathrm{12}}\left(\mathrm{2}-\gamma \right)\right)，

{\sigma }_{\mathrm{4}}=\frac{\mathrm{1}}{P{e}^{\mathrm{4}}}\left(\mathrm{1}-\gamma +\frac{P{e}^{\mathrm{2}}}{\mathrm{12}}\left(\mathrm{2}-\gamma \right)+\frac{P{e}^{\mathrm{4}}}{\mathrm{120}}\right)。

为了进行比较，对四阶Padé（PDE）格式^［14］、古典中心差分（CD）格式和指数四阶紧致（EHOC1^［3］和EHOC2^［4］）格式进行了傅里叶分析，并将结果与本文格式进行了对比。通过计算可 得差分格式Ⅰ、差分格式Ⅱ、PDE格式^［14］、EHOC1格式^［3］、EHOC2格式^［4］和CD格式对应的特征函数可统一表示为

\lambda =\frac{a}{h^{\mathrm{2}}}\frac{T}{T_s}，

其中：T和T_s相应的表达式详见表1。

图4和图5分别给出了在区间［0，π］上，本文两种差分格式及文献中差分格式的耗散误差（Re）和色散误差（Im）比较。图4表明，在网格雷诺数Pe= 0.1， 1， 10时，差分格式Ⅰ和差分格式Ⅱ的耗散误差均优于其他格式；而当Pe增加到100时，差分格式Ⅰ和差分格式Ⅱ的耗散误差显著增大，表现出强耗散性，且差分格式Ⅱ的耗散误差小于差分格式Ⅰ。图5表明，在网格雷诺数Pe= 0.1， 1， 10和100时，差分格式Ⅰ和差分格式Ⅱ的色散误差均小于其他5 种格式。另外，当Pe= 0.1， 1时，差分格式Ⅱ的色散误差小于差分格式Ⅰ，而当Pe增加到10和100时，差分格式Ⅰ的色散误差小于差分格式Ⅱ。

4 数值算例

本节将选取典型算例，编制统一的程序，采用差分格式Ⅰ、差分格式Ⅱ、REC格式^［6］、古典中心差分（CD）格式、FOC格式^［5］、EHOC1格式^［3］和EHOC2^［4］进行计算，通过比较验证本文差分格式的有效性。收敛阶采用的计算式为

\mathrm{R}\mathrm{a}\mathrm{t}\mathrm{e}=\frac{\mathrm{l}\mathrm{n}(E_{\mathrm{1}}/E_{\mathrm{2}})}{\mathrm{l}\mathrm{n}(N_{\mathrm{2}}/N_{\mathrm{1}})}，

其中：E_{\mathrm{1}}{\mathrm{、}}_{}{E}_{\mathrm{2}}分别为网格数取N_{\mathrm{1}}{\mathrm{、}}_{}{N}_{\mathrm{2}}时的最大绝对误差。

4.1 算例1^[4]

考虑常系数非齐次对流扩散方程

边界条件为u(\mathrm{0})=\mathrm{0}{,}_{}u(\mathrm{1})=\mathrm{1}，精确解为

u(x)=\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)+\frac{{\mathrm{e}}^{x/\epsilon }-\mathrm{1}}{{\mathrm{e}}^{\mathrm{1}/\epsilon }-\mathrm{1}}。

当\epsilon很小时，在x=\mathrm{1}处，该解u(x)有一边界层。

表2列出了\epsilon =\mathrm{1.0,0.1,0.01,0.001}时，采用差分格式Ⅰ、差分格式Ⅱ、REC格式^［6］、CD格式、FOC格式^［5］、EHOC1格式^［3］和EHOC2格式^［4］计算的最大误差和收敛阶的比较。从表2中数据可知，差分格式Ⅰ和差分格式Ⅱ均能达到理论上的六阶精度，并且两种格式的最大误差保持在同一量级。随着\epsilon的不断减小，相比文献中差分格式，本文两种差分格式在计算精度上的优势愈发明显。特别地，当\epsilon =\mathrm{0.01}{,}_{}\mathrm{0.001}时，本文差分格式的计算精度远高于REC格式^［6］的计算精度。

4.2 算例2^[4]

考虑常系数非齐次对流扩散方程

其边界条件为u(\mathrm{0})=\mathrm{1}{,}_{}u(\mathrm{1})=\mathrm{l}\mathrm{n}(\mathrm{1}+\epsilon )-\mathrm{1}，精确解为u(x)=\mathrm{l}\mathrm{n}(\mathrm{1}+\epsilon x)+\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{\pi }x)。当\epsilon很小时，该精确解u(x)在x=\mathrm{1}处有一边界层。

表3列出了\epsilon =\mathrm{1.0,0.1,0.01,0.001}时，本文两种差分格式与文献中差分格式的最大误差和收敛阶的比较。从表3中不难发现，本文两种差分格式均能达到理论上的计算精度，并且差分格式Ⅱ的最大误差略小于差分格式Ⅰ。随着\epsilon的不断减小，本文差分格式的优点愈发明显。特别地，当\epsilon =\mathrm{0.01,0.001}时，本文差分格式的计算精度较REC格式^［6］高出2~4个数量级。

5 结论

本文利用待定系数法和截断误差余项修正法，结合一阶和二阶导数的四阶Padé公式，发展了两个求解对流扩散方程的六阶指数型混合紧致差分格式。理论分析和数值算例表明，本文的两种差分格式较文献中的差分格式具有更高的计算精度和稳定性，而且能够有效求解边界层及大梯度问题。本文构造差分格式的思想可以推广到非稳态问题和高维问题的格式构造中。

参考文献

[1] 傅德薰，马延文.计算流体力学［M］.北京：高等教育出版社，2002：191-199.

[2] 陆金甫，关治.偏微分方程数值解法［M］.2版.北京：清华大学出版社， 2004.

[3] TIAN Zhenfu， DAI Shiqiang. High-order compact exponential finite difference methods for convection–diffusion type problems［J］. Journal of Computational Physics， 2007， 220（2）： 952-974.

[4] RADHAKRISHNA PILLAI A C. Fourth-order exponential finite difference methods for boundary value problems of convective diffusion type［J］. International Journal for Numerical Methods in Fluids， 2001， 37（1）： 87-106.

[5] SPOTZ W F. High-order compact finite difference schemes for computational mechanics［D］. The University of Texas at Austin， 1995.

[6] SUN Haiwei， ZHANG Jun. A high-order finite difference discretization strategy based on extrapolation for convection diffusion equations［J］. Numerical Methods for Partial Differential Equations， 2004， 20（1）： 18-32.

[7] TIAN Zhenfu， YU Peixiang. A high-order exponential scheme for solving 1D unsteady convection–diffusion equations［J］. Journal of Computational and Applied Mathematics， 2011， 235（8）： 2477-2491.

[8] 田芳，田振夫.定常对流扩散反应方程非均匀网格上高精度紧致差分格式［J］. 工程数学学报，2009，26（2）： 219-225.

[9] 田芳，葛永斌.求解变系数对流扩散反应方程的指数型高精度紧致差分方法［J］.工程数学学报，2017， 34（3）： 283-296.

[10] TIAN Fang， GE Yongbin， TIAN Zhenfu. Exponential high-order compact finite difference method for convection-dominated diffusion problems on nonuniform grids［J］. Numerical Heat Transfer， Part B： Fundamentals， 2019， 75（3）： 145-177.

[11] YOU D. A high-order Padé ADI method for unsteady convection-diffusion equations［J］. Journal of Computational Physics， 2006， 214（1）： 1-11.

[12] CHU Peter Cheng， FAN Chenwu. A three-point combined compact difference scheme［J］. Journal of Computational Physics， 1998， 140（2）： 370-399.

[13] 梁昌弘，马廷福，葛永斌.两点边值问题的混合型高精度紧致差分格式［J］.宁夏大学学报（自然科学版）， 2017， 38（1）： 1-4.

[14] 马廷福，葛永斌.椭圆型方程两点边值问题的混合型高精度紧致差分格式［J］. 重庆师范大学学报（自然科学版），2019， 36（6）： 87-92.

[15] 王涛，刘铁钢.求解对流扩散方程的一致四阶紧致格式［J］.计算数学，2016， 38（4）： 391-404.

[16] 陆静颖，葛永斌.数值求解一维波动方程的四阶紧致差分方法［J］.宁夏大学学报（自然科学版）， 2020， 41（1）： 17-22.

[17] 袁冬芳，刘文慧，崔桂梅，等.基于人工神经网络的椭圆型微分方程数值求解［J］.宁夏大学学报（自然科学版）， 2021， 43（1）： 6-11.

[18] VARGA R S. Matrix iterative analysis［M］. Englewood Cliffs， NJ： Prentice-Hall， 1962.

基金资助

宁夏自然科学基金资助项目(2020AAC03059)

国家自然科学基金资助项目(11902170)

2024年宁夏回族自治区重点项目引才专项(2024BEH04087)