基于分级网格的有限体积元方法求解奇异摄动两点边值问题

The Finite Volume Element Method Based Improved Grade Mesh for the Singularly Perturbed Two-Point Boundary Value Problems

摘要: 针对一类奇异摄动两点边值问题，本文讨论了基于分级网格及其对偶网格的有限体积元法及其收敛性。在基于分级网格及其对偶网格的基函数及空间的基础上，推导了有限体积元方法的计算格式，所给例子表明所研究计算方法具有稳定和有效性。

Abstract: For a class of singularly perturbed two-point boundary value problems, the finite volume element method based on hierarchical meshes and their dual meshes and its convergence are discussed. Based on the basis function and space of the hierarchical grids and its dual grids, the calculation format of the finite volume element method is deduced. The examples given show that the proposed method is stable and effective.

1. 分级网格及其对偶网格的构造

对奇异摄动两点边值问题，由于小参数 $\epsilon$ 的影响会产生边界层，普通的数值方法在均匀Shishkin网格上很难在靠近边界层处得到理想的数值解。为了得到奇异摄动问题稳定可靠的数值解，应要在其边界层区域放置比边界层外区域更多的网格点，以适应问题的奇异摄动特性。我们采用分级网格来解决边界层问题 [1] [2]。

根据奇异摄动边值问题的边界层，我们可分为以下三类，(1) 边界层存在于 $\left[a,b\right]$ 两端点处，(2) 边界层存在于区间 $\left[a,b\right]$ 的左端点，(3) 边界层存在于区间 $\left[a,b\right]$ 的右端点。本文将着重分析边界层存在于区间 $\left[a,b\right]$ 的左端点。

当边界层存在于区间 $\left[a,b\right]$ 的左端点，将区间 $\left[a,b\right]$ 划分为两个子区间 $\left[a,a+\tau \right]$ 、 $\left[a+\tau ,b\right]$，每个子区间等分为 $\frac{N}{2}$ 个网格点。调用参数 $\tau$ 定义为：

$\tau =\min \left\{\frac{b-a}{2},\lambda \epsilon \ln N\right\}$.

这里 $\lambda$ 由边界层内网格的疏密程度决定， $\lambda$ 是正的常数。比如可取 $\lambda =2.0$。N表示网格节点个数，且N是2的倍数。 $\epsilon$ 是很小的正参数 $0<\epsilon \ll 1$。

其中网格节点定义为

$x_i=\{\begin{array}{l}a+\tau {\left(2i/N\right)}^{\mu },i=0,1,2,\cdots ,N/2,\\ a+\tau +2\left(i-\frac{N}{2}\right)\left(b-a-\tau \right)/N,i=N/2+1,\cdots ,N.\end{array}$

网格的步长定义为：

$h_i=\{\begin{array}{l}\tau {\left[2i/N\right]}^{\mu }-\tau {\left(2\left(i-1\right)/N\right)}^{\mu },i=1,2,\cdots ,N/2,\\ 2\left(b-a-\tau \right)/N,i=N/2+1,\cdots ,N.\end{array}$

我们对应的对偶分级网格点可以表示为：

$J_h^{*}=\left\{x_0,x_{i-1/2}\left(i=1,2,\cdots ,N\right),x_N\right\}.$

其中对偶分级网格节点定义为：

$\begin{array}{l}{x}_{i-1/2}=\{\begin{array}{l}a+\tau {\left[2\left(i-1/2\right)/N\right]}^{\mu },i=1,2,\cdots ,N/2,\\ a+\tau +2\left(i-1/2-\frac{N}{2}\right)\left(b-a-\tau \right)/N,i=N/2+1,\cdots ,N-1,\end{array}\\ x_0=a,x_N=b.\end{array}$

对偶分级网格的步长记为 $h_i^{*}$.

2. 基函数及空间构造

我们针对分级网格，可以得到相应的节点(不包括边界点)基函数如下：

${\phi }_i\left(x\right)=\{\begin{array}{l}\frac{x-x_{i-1}}{h_i},x\in \left[x_{i-1},x_i\right],\\ \frac{x_{i+1}-x}{h_{i\text{+1}}},x\in \left[x_i,x_{i+1}\right],\\ 0,其他.\end{array}i=1,\cdots ,N-1.$

同时，在对偶网格中对应内节点的对偶单元构造相应的分段常数基函数为：

${\psi }_j\left(x\right)=\{\begin{array}{l}1,x\in \left[x_{j-1/2},x_{j+1/2}\right],\\ 0,其他,\end{array}j=1,2,\cdots ,N-1.$

定义线性有限元空间作为测试函数空间 $S_h=\text{span}\left\{{\phi }_1,{\phi }_2.\cdots ,{\phi }_{N-1}\right\}$，分段常数函数空间 $V_h=\text{span}\left\{{\psi }_{\text{1}},{\psi }_{\text{2}}.\cdots ,{\psi }_{N-1}\right\}$ 作为检验函数空间，显然 $\dim S_h=\dim V_h=N-1$。

3. 有限体积元方法及其格式

考虑如下线性奇异摄动两点边值问题 [3]：

$\{\begin{array}{l}-\epsilon \widetilde{y^{″}}\left(x\right)+p\left(x\right)\widetilde{y^{\prime }}\left(x\right)+q\left(x\right)\tilde{y}\left(x\right)=\tilde{f}\left(x\right),x\in \left(a,b\right),\\ \tilde{y}\left(a\right)=\alpha ,\tilde{y}\left(b\right)=\beta ,\end{array}$ (1)

其中 $\epsilon$ 是很小的正参数 $0<\epsilon \ll 1$，和 $\beta$ 是已知量。假设在 $\left[a,b\right]$ 上 $p\left(x\right),q\left(x\right),f\left(x\right)$ 是充分光滑函数。

作如下变量变换，我们得到：

$y\left(x\right)=\tilde{y}\left(x\right)-\left[\alpha +\frac{\beta -\alpha }{b-a}\left(x-a\right)\right],$

就可以将(1)化为如下齐次边界条件问题：

$\{\begin{array}{l}-\epsilon y^{″}\left(x\right)+p\left(x\right)y^{\prime }\left(x\right)+q\left(x\right)y\left(x\right)=f\left(x\right),x\in \left(a,b\right),\\ y\left(a\right)=0,y\left(b\right)=0,\end{array}$ (2)

其中 $f\left(x\right)=\tilde{f}\left(x\right)-\left(p\left(x\right)+\left(x-a\right)q\left(x\right)\right)\left(\beta -\alpha \right)/\left(b-a\right)$。

对以上奇异摄动方程(2)对偶剖分中任意单元 $I_i^{*}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ 上积分

$-\epsilon {\displaystyle {\int }_{I_i^{*}}{y}^{″}\left(x\right)\text{d}x}+{\displaystyle {\int }_{I_i^{*}}p\left(x\right)y^{\prime }\left(x\right)\text{d}x}+{\displaystyle {\int }_{I_i^{*}}q\left(x\right)y\left(x\right)\text{d}x}={\displaystyle {\int }_{I_i^{*}}f\left(x\right)\text{d}x}.$ (3)

上式具体可以写成

(4)

取相应的有限元 $y_h\in S_h$ 作y近似

$\begin{array}{l}\epsilon \left[{y^{\prime }}_h\left(x_{i-1/2}\right)-{y^{\prime }}_h\left(x_{i+1/2}\right)\right]+{\displaystyle {\int }_{x_{i-1/2}}^{x_{i+1/2}}p\left(x\right){y^{\prime }}_h\left(x\right)\text{d}x}+{\displaystyle {\int }_{x_{i-1/2}}^{x_{i+1/2}}q\left(x\right)y_h\left(x\right)\text{d}x}\\ ={\displaystyle {\int }_{x_{i-1/2}}^{x_{i+1/2}}f\left(x\right)\text{d}x},i=1,2,\cdots ,N-1.\end{array}$ (5)

由于 $y_h\in S_h$，所以满足

$y_h\left(x\right)={\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{y}_i{\phi }_i\left(x\right)}.$ (6)

在区间 $\left[x_{i-1},x_i\right]$ 中，有

将表达式(6)代入(5)，我们就可以得到一个关于 $y_1,y_2,\cdots ,y_{N-1}$ 的线性代数方程组

$\begin{array}{l}\epsilon \left(-\frac{1}{h_i}{y}_{i-1}+\left(\frac{1}{h_i}+\frac{1}{h_{i+1}}\right)y_i-\frac{1}{h_{i+1}}{y}_{i+1}\right)+[-\frac{1}{h_i}\left({\displaystyle {\int }_{x_{i-1/2}}^{x_i}p\left(x\right)\text{d}x}\right)y_{i-1}\\ +\left(\frac{1}{h_i}{\displaystyle {\int }_{x_{i-1/2}}^{x_i}p\left(x\right)\text{d}x}-\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+1/2}}p\left(x\right)\text{d}x}\right)y_i+\frac{1}{h_{i+1}}\left({\displaystyle {\int }_{x_i}^{x_{i+1/2}}p\left(x\right)\text{d}x}\right)y_{i+1}]\\ +[\frac{1}{h_i}\left({\displaystyle {\int }_{x_{i-1/2}}^{x_i}q\left(x\right)\left(x_i-x\right)\text{d}x}\right)y_{i-1}+\left(\frac{1}{h_i}{\displaystyle {\int }_{x_{i-1/2}}^{x_i}q\left(x\right)\left(x-x_{i-1}\right)\text{d}x}+\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+1/2}}q\left(x\right)\left(x_{i+1}-x\right)\text{d}x}\right)y_i\\ +\frac{1}{h_{i+1}}\left({\displaystyle {\int }_{x_i}^{x_{i+1/2}}q\left(x\right)\left(x-x_i\right)\text{d}x}\right)y_{i+1}]={\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f\left(x\right)\text{d}x},i=1,2,\cdots ,N-1.\end{array}$ (7)

注意，在边界上： $y_0=y_N=0$，令上式中的矩阵

$A={\left(\begin{array}{cccccccc}\frac{1}{h_1}+\frac{1}{h_2}& -\frac{1}{h_2}& 0& 0& \cdots & 0& 0& 0\\ -\frac{1}{h_2}& \frac{1}{h_2}+\frac{1}{h_3}& -\frac{1}{h_3}& 0& \cdots & 0& 0& 0\\ 0& -\frac{1}{h_3}& \frac{1}{h_3}+\frac{1}{h_4}& -\frac{1}{h_4}& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -\frac{1}{h_{N-2}}& \frac{1}{h_{N-2}}+\frac{1}{h_{N-1}}& -\frac{1}{h_{N-1}}\\ 0& 0& 0& 0& \cdots & 0& -\frac{1}{h_{N-1}}& \frac{1}{h_{N-1}}+\frac{1}{h_N}\end{array}\right)}_{\left(N-1\right)\times \left(N-1\right)},$

$B={\left(\begin{array}{cccccccc}{P}_1& {\tilde{P}}_1& 0& 0& \cdots & 0& 0& 0\\ {\bar{P}}_2& P_2& {\tilde{P}}_2& 0& \cdots & 0& 0& 0\\ 0& {\bar{P}}_3& P_3& {\tilde{P}}_3& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & {\bar{P}}_{N-2}& P_{N-2}& {\tilde{P}}_{N-2}\\ 0& 0& 0& 0& \cdots & 0& {\bar{P}}_{N-1}& P_{N-1}\end{array}\right)}_{\left(N-1\right)\times \left(N-1\right)},$

其中

$\begin{array}{l}{P}_i=\frac{1}{h_i}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_i}p\left(x\right)\text{d}x}-\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+\frac{1}{2}}}p\left(x\right)\text{d}x}\text{,}i=1,2,\cdots N-1.\\ {\tilde{P}}_i=\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+\frac{1}{2}}}p\left(x\right)\text{d}x},i=1,2,\cdots ,N-1.\\ {\bar{P}}_i=-\frac{1}{h_i}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_i}p\left(x\right)\text{d}x},i=1,2,\cdots ,N-1.\end{array}$

$C={\left(\begin{array}{cccccccc}{Q}_1& {\tilde{Q}}_1& 0& 0& \cdots & 0& 0& 0\\ {\bar{Q}}_2& Q_2& {\tilde{Q}}_2& 0& \cdots & 0& 0& 0\\ 0& {\bar{Q}}_3& Q_3& {\tilde{Q}}_3& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & {\bar{Q}}_{N-2}& Q_{N-2}& {\tilde{Q}}_{N-2}\\ 0& 0& 0& 0& \cdots & 0& {\bar{Q}}_{N-1}& Q_{N-1}\end{array}\right)}_{\left(N-1\right)\times \left(N-1\right)},$

其中

$\begin{array}{l}{Q}_i=\frac{1}{h_i}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_i}q\left(x\right)\left(x-x_{i-1}\right)\text{d}x}+\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+\frac{1}{2}}}q\left(x\right)\left(x_{i+1}-x\right)\text{d}x}\text{,}i=1,2,\cdots ,N-1.\\ {\tilde{Q}}_i=\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+\frac{1}{2}}}q\left(x\right)\left(x-x_i\right)\text{d}x},i=1,2,\cdots ,N-1.\\ {\bar{Q}}_i=\frac{1}{h_i}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_i}q\left(x\right)\left(x_i-x\right)\text{d}x},i=1,2,\cdots ,N-1.\end{array}$

$Y={\left(y_1,y_2,\cdots ,y_{N-1}\right)}^{\text{T}},$

$F={\left({\displaystyle {\int }_a^{b}f\left(x\right){\psi }_{\text{1}}\left(x\right)\text{d}x},{\displaystyle {\int }_a^{b}f\left(x\right){\psi }_2\left(x\right)\text{d}x},\cdots ,{\displaystyle {\int }_a^{b}f\left(x\right){\psi }_{N-1}\left(x\right)\text{d}x}\right)}^{\text{T}}.$

则(7)式化为如下线性矩阵形式代数方程组

$\left(\epsilon A+B+C\right)Y=F.$

我们对以上的有限体积元方法的研究，有文献 [4] 可以得到了如下的结论：

定理3.1 $y_h\in S_h$ 是方程(2)的有限体积元解，在分级网格 $J_h$ 的所有节点上

$\left(y-y_h\right)\left(x_j\right)={\rm O}\left(N^{-2}\right),j=1,2,\cdots ,N-1.$

在分级网格 $J_h$ 上，其超收敛性估计：

$\left|{\left(y-y_h\right)}^{\prime }\left(z\right)\right|\le C{N}^{-1}\text{.}$

4. 数值例子

本小节给出如下数值例子来验证有限体积元法的有效性。

例题 考虑如下奇异摄动两点边值问题 [5]：

$\{\begin{array}{l}-\epsilon y^{″}\left(x\right)-y^{\prime }\left(x\right)+y\left(x\right)=-1,x\in \left(0,1\right),\\ y\left(0\right)=0,y\left(1\right)=0.\end{array}$ (8)

其精确解为

$y\left(x\right)=\frac{\left({\text{e}}^{m_2}-1\right){\text{e}}^{m_{1}x}+\left(1-{\text{e}}^{m_1}\right){\text{e}}^{m_{2}x}}{{\text{e}}^{m_2}-{\text{e}}^{m_1}}-1,$

其中 $m_1=\frac{1}{2\epsilon }\left(-1+\sqrt{1+4\epsilon }\right),m_2=\frac{1}{2\epsilon }\left(-1-\sqrt{1+4\epsilon }\right)$。

微分方程(8)中 $p\left(x\right)=-1,q\left(x\right)=1,f\left(x\right)=-1$，对应的有限体积元格式为：

$\begin{array}{l}\epsilon \left(-\frac{1}{h_i}{y}_{i-1}+\left(\frac{1}{h_i}+\frac{1}{h_{i+1}}\right)y_i-\frac{1}{h_{i+1}}{y}_{i+1}\right)+[-\frac{1}{h_i}\left({\displaystyle {\int }_{x_{i-1/2}}^{x_i}\left(-1\right)\text{d}x}\right)y_{i-1}+\left(\frac{1}{h_i}{\displaystyle {\int }_{x_{i-1/2}}^{x_i}\left(-1\right)\text{d}x}-\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+1/2}}\left(-1\right)\text{d}x}\right)y_i\\ +\frac{1}{h_{i+1}}\left({\displaystyle {\int }_{x_i}^{x_{i+1/2}}\left(-1\right)\text{d}x}\right)y_{i+1}]+[\frac{1}{h_i}\left({\displaystyle {\int }_{x_{i-1/2}}^{x_i}\left(x_i-x\right)\text{d}x}\right)y_{i-1}+\left(\frac{1}{h_i}{\displaystyle {\int }_{x_{i-1/2}}^{x_i}\left(x-x_{i-1}\right)\text{d}x}+\frac{1}{h_{i+1}}{\displaystyle {\int }_{x_i}^{x_{i+1/2}}\left(x_{i+1}-x\right)\text{d}x}\right)y_i\\ +\frac{1}{h_{i+1}}\left({\displaystyle {\int }_{x_i}^{x_{i+1/2}}\left(x-x_i\right)\text{d}x}\right)y_{i+1}]={\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\text{d}x},i=1,2,\cdots ,N-1.\end{array}$

本例是左边界层，取 $\tau =\min \left\{0.5,\lambda \epsilon \ln N\right\}$ 将区间 $\left[0,1\right]$ 划分为两个子区间 $\left[0,\tau \right]$ 、 $\left[\tau ,1\right]$，每个子区间等分为 $\frac{N}{2}$ 个网格点。在分级网格下，网格节点为：

区间步长为：

$h_i=\{\begin{array}{l}\tau \left[{\left(2i/N\right)}^{\mu }-{\left(2\left(i-1\right)/N\right)}^{\mu }\right],i=1,2,\cdots ,\frac{N}{2},\\ 2\left(1-\tau \right)/N,i=\frac{N}{2}+1,\cdots ,N.\end{array}$

对于相同的小参数 $\epsilon$ 和剖分数N，针对不同的 $\mu$ 值，得到分级网格的数值解，我们比较其与精确解的误差，如表1和表2。

通过表1和表2，我们可以看出，在分级网格的剖分下，数值解与精确解的误差相对来说较小。且当 $\mu$ 越大，剖分精度越高，与精确解的误差越小。通过以上数据，我们可以得到，分级网格可以更好地计算边界层的数值解，靠近边界层处节点误差与其他节点的误差同阶。

Table 1. Finite volume element error when ε = 10 − 3 , N = 1000 , λ = 2

表1. 当 $\epsilon ={10}^{-3},N=1000,\lambda =2$ 时的有限体积元误差

Table 2. Finite volume element error when ε = 10 − 3 , N = 1000 , λ = 2

表2. 当 $\epsilon ={10}^{-3},N=1000,\lambda =2$ 时的有限体积元误差

基金项目

本文为国家自然科学基金项目(11571102)资助课题。