本文考虑如下带间断系数的奇异摄动对流扩散方程

其中 $0<\varepsilon\ll 1$ 是摄动参数, $\Omega^{-}=(0,d)$, $\Omega^{+}=(d,1)$, $b,c,f\in C^{2k}(\Omega^{-}\cup \Omega^{+})$, 且对于给定的常数 $\gamma$, 满足

此外, 假设存在与 $\varepsilon$ 无关的常数 $\beta_1,\beta_{1}^{ \star},\beta_2,\beta_{2}^{ \star}$, $C$, 使得

在这里, 对于任意函数 $w(x)$, 定义 $[w(d)]=w(d^{+})-w(d^{-})$.方程 (1.1) 的解 $u(x)$ 通常在不连续点 $x=d$ 处产生一个宽度为 $O(\varepsilon \ln(1/\varepsilon))$ 的内点层^[1].

具有非光滑系数的奇异摄动微分方程在流体力学, 量子力学, 弹性力学, 气体多孔电极理论, 气象学, 海洋学等领域有着广泛的应用. 这类问题的一个显著特点是高阶导数项包含一个摄动参数 $\varepsilon$, 且对流项和右边源项在区域内部存在间断点. 当 $\varepsilon\to 0$ 时, 这类问题的解在区域的边界和内部某个区域内会产生剧烈的振荡, 即存在所谓的边界层和内点层现象^[2]. 正因为边界层或内点层的出现, 在很多情况下很难给出这类问题的精确解.

因此, 目前常用的数值方法是在层适应网格 (Shishkin 网格和 Bakhvalov 网格) 下构造出相应的迎风有限差分格式, 并证明离散格式的一致收敛性^[3,4]. 然而, 据我们所知, 这类方法的收敛性一般只有一阶或者二阶精度.

与有限差分方法相比, 有限元方法更方便处理复杂区域, 且易构造更高阶格式. 因此, 少数学者研究了带间断系数的奇异摄动对流扩散方程的有限元方法. 例如,Zarin 等^[5]在 Shishkin 网格下构造了带间断系数的奇异摄动对流扩散方程的标准有限元方法, 并在能量范数下证明了数值格式是几乎一阶一致收敛的. 在此基础上, 张进等^[1] 给出了能量范数下几乎 $k$ 阶一致收敛的证明. Babu 等^[6]在 Shishkin 网格下使用流线扩散有限元方法证明了数值格式在最大范数下是几乎二阶收敛的. 据我们查阅相关资料发现, 在 Bakhvalov 网格下, 目前还没有相关文献报道问题 (1.1) 的高阶有限元方法, 其困难主要源于标准的 Lagrange 插值不适用于 Bakhvalov 网格下有限元方法的一致收敛分析. 更具体地说, Lagrange 插值不能在网格过渡点附近提供足够的 $L^2$ 范数稳定性, 所以有限元方法在 Bakhvalov 网格上的研究进展缓慢.

Linß 等^[7]指出, 虽然 Galerkin 有限元在层适应网格上能获得一致收敛的结果, 但是数值结果对层适应网格中参数的选择较敏感, 所以仍然有必要在层适应网格上使用稳定性较强的数值方法.自 1970 年以来, 非对称内惩罚伽辽金 (NIPG) 方法逐渐成为一种流行的稳定技术. 由于该方法采用了内部惩罚项来约束单元边界的不连续, 因此它具有协调有限元所没有的灵活性和优点. 此外, 与 DG 方法的其他变体—例如对称内罚伽辽金方法和不完全内罚伽辽金方法相比, 该方法对任何非负惩罚参数和任何网格都具有稳定性和收敛性. 鉴于其自身独特的优势, NIPG 方法引起了越来越多研究者的关注^[8⇓⇓-11].

基于此, 本文针对带间断系数的奇异摄动对流扩散方程 (1.1), 构造了一种 Bakhvalov-type 网格. 通过设计一种新的复合型插值方法, 它由层外的 Gauß Radau 投影和层内的 Lagrange 插值组成, 从而推导出 Bakhvalov-type 网格下 NIPG 方法的最优一致收敛性. 最后的数值实验支持了我们的理论结果.

假设 1.1 我们将在整篇文章中假设 $\varepsilon \leq CN^{-1}$, 这在实践中并不是一种限制.

为了后续误差分析的需要, 在这部分给出连续解的存在性、稳定性及导数估计, 见如下引理 2.1-2.3

引理 2.1^[1] 若 $u(x)$ 是问题 (1.1) 的解, 则 $u(x)\in C^{0}(\bar{\Omega})\cap C^{1}(\Omega)\cap C^{2}(\Omega^{-}\cup \Omega^{+})$.

引理 2.2^[1] 若问题 (1.1) 的解 $u(x)\in C^{0}(\bar{\Omega})\cap C^{1}(\Omega)\cap C^{2}(\Omega^{-}\cup \Omega^{+})$, 则

其中 $\lambda=\min\{\beta_1/d,\beta_2/(1-d)\}$.

引理 2.3^[1] 假设问题 (1.1) 的解可分解为 $u=S+E$, 其中光滑部分 $S$ 和层部分 $E$ 分别满足 $LS=f$ 和 $LE=0$. 那么对于任意 $0\leq i \leq k+1$, 有

由引理 2.3 可知, 问题 (1.1) 的解在 $x = d$ 处存在内点层. 因此, 为了构造出有效的 Bakhvalov-type 网格, 我们首先将 $\bar{\Omega}$ 分为四个子区间

其中 $\tau_1=d+\frac{\sigma \varepsilon }{\beta_1}\ln \varepsilon$, $\tau_2=d-\frac{\sigma \varepsilon}{\beta_2}\ln \varepsilon$, 且满足 $d-\tau_1\leq d/2$ 和$\tau_2-d\leq (1-d)/2$.假设 $N\in \mathbb{N},N\geq 16$ 是一个能被 4 整除的整数并且 $\sigma\geq k+1$.我们将上述四个子区间分别分成 $N/4$ 个小区间, 则每一个网格节点 $x_j$ 可以表示为

Bakhvalov-type 网格具有如下性质

引理 3.1 令 $h_{j}=x_{j}-x_{j-1}$, 则有

进一步, 对于 $0 \leq \nu \leq \sigma$, 有

证 证明见文献 [12,引理 3].

令 $\mathcal{T}_{N}=\{ I_{j}=[x_{j-1},x_{j}]:j=1,\cdots,N \}$ 为区域 $\Omega$ 的一个网格剖分, 则对任意单元 $I_{j}\in \mathcal{T}_{N}$, 定义 $s$ 阶的分片 Sobolev 空间

和相应的分片 Sobolev 范数和半范数

其中 $\left\|\cdot\right\|_{s,I_{j}}$ 和 $\left|\cdot\right| _{s,I_{j}}$ 是 $H^{s}(I_{j})$ 中常见的 Sobolev 范数和半范数. 在这里, $\left\|\cdot\right\|_{I}$ 和 $(\cdot,\cdot)_{I}$ 分别表示 $L^{2}(I)$ 范数和 $L^{2}(I)$ 内积.

为了构造出问题 (1.1) 的 NIPG 有限元方法, 在 Bakhvalov-type 网格下, 定义如下有限元空间

其中 $\mathbb{P}_{k}(I_{j})$ 是 $I_{j}$ 上次数不超过 $k$ 次的多项式空间. $V_{N}^{k}$ 空间中的函数在单元边界处是完全间断的.

对于给定的函数 $u(x)\in V_{N}^{k}$, 令$u(x_{j}^{+})$ 和 $u(x_{j}^{-})$ 分别表示 $u$ 在 $x_{j}$ 的右侧和左侧取值, 则函数 $u$ 在内部节点 $x_{j}$ 处的跳跃和平均为

为表述方便, 我们将边界节点 $x_{0}$ 和 $x_{N}$ 处的跳跃和平均定义为

基于以上的定义, 我们给出问题 (1.1) 的 NIPG 有限元离散格式: 寻找 $u_{N}\in V_{N}^{k}$ 使得

其中, 对任意 $u,v \in H^{1}(\Omega,\mathcal{T}_{N})$, 算子 $a_j(\cdot,\cdot)(j=1,2)$, $L(\cdot)$ 的定义如下

这里 $\lambda_{j}(j=0,\cdots,N)$ 是与节点 $x_{j}$ 相关的非负惩罚参数, 本文取 $\lambda_{j}=N, j=0,\cdots,N$.

引理 3.2 若 $u$ 和 $u_{N}$ 分别是问题 (1.1) 和 (3.10) 的精确解, 则双线性形式 $a(\cdot,\cdot)$ 满足如下 Galerkin 正交性

证 类似于文献 [11,引理 2.2].

为了讨论双线性形式 $a(\cdot,\cdot)$ 的稳定性, 我们定义如下能量范数

通过简单计算, 我们可以得到

基于 Lax-Milgram 引理, 问题 (3.10) 存在唯一解 $u_N$.

引入 $u$ 的一个复合型插值 $\Pi u$, 即

其中 $L_{k}u$ 是 $u$ 的 $k$ 次 Lagrange 插值, 而 $P_{h}u$ 是 $u$ 的 Gauß Radau 投影. 同时给出 Gauß Radau 投影 $(P_{h}w)|_{I_{j}}$ 的定义^[13]: 对于 $j=1,2,\cdots,N$,

根据文献 [13] 中的投影结果, 对所有 $w\in H^{k+1}(I_j),j=1,2,\cdots,N$, 有

由 Lagrange 插值算子的定义以及 Sobolev 空间中的插值理论^[14], 对所有 $w\in W^{k+1,p}(I_{j}),$$j=1,2,\cdots,N$, 有

其中 $l=0,1$ 且 $1\leq p,q \leq \infty$.

引理 4.1 假设 $\sigma \geq k+1$, 在 Bakhvalov-type 网格上有

证 利用文献 [15,引理 2] 以及 (4.4) 和 (4.6) 式可得出该引理.

情况1 通过使用引理 4.1 中的方法, 我们可以很容易地在 $[d,1]$ 上得到类似的结果, 即

引理 4.2^[16] 令 $w\in H^{1}(I_j)$, 则

引理 4.3 令 $\eta=u-\Pi u$, 在 Bakhvalov-type 网格上有

证 类似地利用文献 [16,引理 3.6] 提出的方法, 可以很容易得到该引理.

由引理 4.1 和情况 1 可以推出下面插值误差的估计结果.

定理 4.1 令 $\sigma \geq k+1$, 那么在 Bakhvalov-type 网格上, 有

引入 $\xi=\Pi u-u_N$. 由 Galerkin 正交性, (3.13) 式和分部积分可得

接下来, 我们将分别分析 (5.1) 式右端的每一项 $\Gamma_i(i=1,\cdots,7)$.

首先, 利用 Hölder 不等式和 (4.13) 式可得

对于 $\Gamma_2$, 根据 (4.12) 式和 $\lambda_j=N$, 可得

为了估计 $\Gamma_3$, 我们首先将 $\Gamma_3$ 分解成

根据 (4.1) 式和 Lagrange 插值的连续性质, 可知 $[\eta(x_j)]=0,j=\frac{N}{4}+2,\cdots,\frac{3N}{4}-2$. 同时, 利用逆不等式, (4.7), (3.2), (3.3) 式和 $\varepsilon\leq N^{-1}$, 可以得到

同理可得

综上可得

对于 $\Gamma_4$ 的分析可以类似于 $\Gamma_3$. 利用 (4.7) 式, 可得

现在我们来考虑 $\Gamma_5$ 和 $\Gamma_6$, 将其分成 $1\leq j \leq \frac{N}{4}+1$, $\frac{N}{4}+2 \leq j \leq \frac{3N}{4}-1$ 和 $\frac{3N}{4}\leq j \leq N$ 三部分进行分析. 对于 $1\leq j \leq \frac{N}{4}+1$, 利用 (4.2) 和 (4.3) 式可以得到

其中 $b(x_{j-\frac{1}{2}})$ 是 $b(x)$ 在区间 $[x_{j-1},x_j]$ 的中点 $x_{j-\frac{1}{2}}$ 处的取值. 由于 $b(x)$ 是一个充分光滑的函数, 那么利用 Lagrange 中值定理可知在 $x$ 和 $x_{j}$ 之间存在一点 $\chi$ 使得 $b(x)-b(x_{j-\frac{1}{2}})=b'(\chi)(x-x_{j-\frac{1}{2}})$ 成立. 因此, 使用逆不等式和 (4.7) 式, 可得

对于 $\frac{3N}{4}\leq j \leq N$, 同理可得

对于 $\frac{N}{4}+2 \leq j \leq \frac{3N}{4}-1$, 我们利用逆不等式, (3.8), (3.9), (4.4), (4.6), (4.7) 式, 通过简单推导可得

综上所述, 我们可以得到以下估计

对于 $\Gamma_7$, 利用 Hölder 不等式和 (4.13) 式可得

最后, 将 (5.2)-(5.7) 式代入到 (4.13) 式, 有

现在我们给出本文的主要结论.

定理 5.1 设 $u$ 和 $u_N$ 分别是问题 (1.1) 和问题 (3.10) 的解, 并令 $\lambda_j=N(j=0,1,\cdots,N)$. 则在 Bakhvalov-type 网格上有

证 由三角不等式, (3.10) 和 (5.8) 式可得

证毕.

为了验证本文所提出的 NIPG 有限元方法的有效性, 我们用两个具体例子来进行数值实验. 考虑到精确解难以求出, 我们使用如下双网格原理来计算数值解的能量范数

基于该误差估计 (6.1) 式, 定义收敛阶的计算公式如下

在以下所有的计算中, 取 $\sigma=k+1$.

例 1考虑如下常系数对流扩散方程

其中

当 $\varepsilon=10^{-j}(j=3,\cdots,8)$ 和 $N=16,32,\cdots,1024$ 时, 表1 和表2 分别给出了在 $k=1$ 和 $k=2$ 时, NIPG 方法计算得到的能量误差和相应的收敛阶, 其中每一个 $\varepsilon$ 对应的第一行表示能量误差 $E^{N}$, 第二行表示相应的收敛阶 $\rho^{N}$. 图1 给出了$N=128,\varepsilon=10^{-4}$ 时, 数值解的变化曲线图. 进一步验证了数值解在不连续点 $x=0.5$ 附近存在内点层.

例 2 考虑如下变系数对流扩散方程

其中

对于不同的 $\varepsilon$ 和 $N$, 表3 和表4 分别列出了 $k=1,2$ 时, NIPG 方法计算问题 (6.4) 得到的能量误差和收敛阶. 当 $N=128,\varepsilon=10^{-4}$ 时, 图2 画出了数值解的曲线图. 该图验证了数值解在不连续点 $x=0.7$ 附近存在内点层.

综合上述两例, 由表1-4 中的数值结果可以看出, 本文提出的 NIPG 方法在 Bakhvalov-type 网格下能达到 $k$ 阶最优一致收敛.