求解奇异摄动两点边值问题的奇性分离法

The Singularity-Separated Method for the Singularly Perturbed Two-Point Boundary Value Problem

摘要: 本文使用奇性分离法求解奇异摄动两点边值问题。首先通过修改边界条件得到弱奇性的第三边值辅助问题，将其解记为w(x)，其次利用特征值构造一个奇异函数v(x)，最后将原两点边值问题的解u(x)表示为u(x)=w(x)-v(x)。由于将解的奇性进行了分离，数值求解时不必使用局部加密网格。数值实验中边界层仅用1个单元的稀疏网格就能得到高精度的有限元解。

Abstract: The singularity-separated method is used to solve the singularly perturbed two-point boundary value problem. The third boundary value problem whose solution is w(x) is constructed by modifying the boundary-value condition and a singular function v(x) is constructed by the eigenvalues. Then the solution u(x) of the two-point boundary value problem can be expressed as w(x)-v(x). Numerical results show that the FE-solutions have the high accuracy instead of local refinement meshgrid.

1. 引言

奇异摄动问题是科学研究和工程实践中的常见问题。由于微分方程中的二阶导数项含有小参数 $\epsilon >0$，方程的解将在很窄的区间上发生剧烈变化，出现边界层或内部层，导致方程求解困难。

我们考虑一维奇异摄动两点边值问题

$\{\begin{array}{ll}Au=-\epsilon u^{″}+b{u}^{\prime }+cu=f\left(x\right)\hfill & \text{in}\left[0,1\right]\hfill \\ u\left(0\right)=0,u\left(1\right)=0\hfill & \hfill \end{array}$ (1)

其中小参数 $\epsilon ={10}^{-3}\sim {10}^{-10}$，常数 $b>0$， $c>0$。该问题的解在右端点 $x=1$ 附近出现一个边界层，记其宽度为 $\tau$。早期在均匀网格上采用有限差分法进行求解时会产生强烈的非物理振荡 [1]，而使用Shishkin网格后，差分法 [2]、有限元法 [3] 以及DG方法 [4] 等各种数值方法都能取得高精度的数值解，所以Shishkin网格能有效求解奇异摄动问题。

然而在处理多维问题时，Shishkin网格的使用会出现一些不足。例如二维区域上使用Shishkin网格时，狭窄的边界层内集中了大量节点，使得边界附近很多剖分单元的长宽比很大，而区域角点附近的剖分单元非常小，在使用有限元法进行计算时，这种网格的划分会对稳定性和精确度等方面带来许多影响。若能将奇异解的奇性进行分离，那么计算时就可能不必使用Shishkin网格也能得到高精度的数值解。

2. 奇性分离法

奇异摄动两点边值问题(1)的通解为

$u\left(x\right)=u_0\left(x\right)-C_1{\varphi }_1\left(x\right)-C_2{\varphi }_2\left(x\right),{\varphi }_1\left(x\right)={\text{e}}^{{\lambda }_{1}x},{\varphi }_2\left(x\right)={\text{e}}^{{\lambda }_2\left(x-1\right)},$

其中 $u_0\left(x\right)$ 为问题(1)的一个特解， $C_1{\varphi }_1\left(x\right)+C_2{\varphi }_2\left(x\right)$ 为对应的齐次方程 $Au=0$ 的通解，

${\lambda }_1=\frac{-2c}{b+\sqrt{b^2+4c\epsilon }}\approx -\frac{c}{b}$， ${\lambda }_2=\frac{b+\sqrt{b^2+4c\epsilon }}{2\epsilon }\approx \frac{b}{\epsilon }$。

由边界条件得 $u\left(0\right)=u_0\left(0\right)-C_1-C_2{\text{e}}^{-{\lambda }_2}=0$， $u\left(1\right)=u_0\left(1\right)-C_1{\text{e}}^{{\lambda }_1}-C_2=0$，由于 ${\text{e}}^{-{\lambda }_2}$ 非常接近0，可取 $C_1=u_0\left(0\right)$， $C_2=u_0\left(1\right)-u_0\left(0\right){\text{e}}^{{\lambda }_1}$。我们定义正则函数 $w\left(x\right)=u_0\left(x\right)-C_1{\varphi }_1\left(x\right)$ 以及奇异函数 $v\left(x\right)=C_2{\varphi }_2\left(x\right)$，于是(1)的解 $u\left(x\right)$ 可以表示为 $u\left(x\right)=w\left(x\right)-v\left(x\right)$，这样就实现了对解中奇性的分离。

故我们提出如下奇性分离法 [5]：

1) 构造第三边值辅助问题

$\{\begin{array}{ll}Aw=-\epsilon w^{″}+b{w}^{\prime }+cw=f\left(x\right)\hfill & \text{in}\left[0,1\right]\hfill \\ w\left(0\right)=0,A_{1}w\left(1\right)=b{w}^{\prime }\left(1\right)+cw\left(1\right)=f\left(1\right)\hfill & \hfill \end{array}$ (2)

这样有 $\epsilon w^{″}\left(1\right)=0$，辅助问题解的奇性大大减弱。

2) 构造奇异函数

$v\left(x\right)=w\left(1\right){\varphi }_2\left(x\right),{\varphi }_2\left(x\right)={\text{e}}^{{\lambda }_2\left(x-1\right)},$

它满足 $Av=0$，且 $v\left(1\right)=w\left(1\right)$， $v\left(0\right)=w\left(1\right){\text{e}}^{-{\lambda }_2}\approx 0$。

3) 记 $z=u-\left(w-v\right)$，满足

$\{\begin{array}{l}Az=-\epsilon z^{″}+b{z}^{\prime }+cz=0\\ z\left(0\right)=v\left(0\right)=w\left(1\right){\text{e}}^{-{\lambda }_2}\approx 0,z\left(1\right)=0.\end{array}$ (3)

我们将证明

定理1问题(3)的解 $z=u-\left(w-v\right)$ 有以下估计

$\left|z\left(x\right)\right|\le \left|v\left(0\right)\right|,\left|z^{\prime }\left(x\right)\right|\le C\left|v\left(0\right)\right|/\epsilon ,v\left(0\right)=o\left({\text{e}}^{-b/\epsilon }\right).$

根据此定理，可以认为 $z=0$，于是可以把问题(1)的解表示为

$u\left(x\right)=w\left(x\right)-v\left(x\right)+o\left({\text{e}}^{-b/\epsilon }\right),u^{\prime }\left(x\right)=w^{\prime }\left(x\right)-v^{\prime }\left(x\right)+o\left({\text{e}}^{-b/\epsilon }\right)/\epsilon .$

3. 定理的证明

定理1问题

$\{\begin{array}{ll}Az=-\epsilon z^{″}+b{z}^{\prime }+cz=0\hfill & \text{in}\left[0,1\right]\hfill \\ z\left(0\right)=v\left(0\right),z\left(1\right)=0\hfill & \hfill \end{array}$ (4)

的解z有以下估计

$\left|z\left(x\right)\right|\le \left|v\left(0\right)\right|,\left|z^{\prime }\left(x\right)\right|\le C\left|v\left(0\right)\right|/\epsilon .$

证明：记 $M=\max \left|z\left(x\right)\right|$。

假设z在区间[0, 1]内一点 $x_0$ 处取得正的最大值 $z\left(x_0\right)$，那么 $z^{\prime }\left(x_0\right)=0$， $z^{″}\left(x_0\right)<0$，则 $Az\left(x_0\right)>0$，与 $Az=0$ 矛盾，故z若有正的最大值则应该发生在端点处，于是有

$z\left(x\right)\le \max \left\{\max \left\{z\left(0\right),z\left(1\right)\right\},0\right\}\le \left|v\left(0\right)\right|.$

同理，若z在区间[0, 1]内一点 $x_1$ 处取得负的最小值 $z\left(x_1\right)$，那么 $z^{\prime }\left(x_1\right)=0$， $z^{″}\left(x_1\right)>0$，则 $Az\left(x_1\right)<0$，与 $Az=0$ 矛盾，可得

$z\left(x\right)\ge \min \left\{\min \left\{z\left(0\right),z\left(1\right)\right\},0\right\}\ge -\left|v\left(0\right)\right|.$

下面证明第二个估计。

在方程 $Az=0$ 的两边同乘x并在区间 $\left(0,x\right)$ 上分部积分得

$-\epsilon x{z}^{\prime }+\epsilon {\displaystyle {\int }_0^x{z}^{\prime }\text{d}x}+bxz\left(x\right)-{\displaystyle {\int }_0^{x}bz\text{d}x}+{\displaystyle {\int }_0^{x}cxz\text{d}x}=0,$

$-\epsilon x{z}^{\prime }+bxz+\epsilon z\left(x\right)-\epsilon z\left(0\right)+{\displaystyle {\int }_0^x\left(-bz+cxz\right)\text{d}x}=0,$

因此 $\epsilon x\left|z^{\prime }\right|\le CM$，当 $x\ge \frac{1}{2}$ 时有 $\epsilon \left|z^{\prime }\right|\le CM$。

类似地，在方程 $Az=0$ 的两边同乘 $x-1$ 并在 $\left(x,1\right)$ 上分部积分得

$\epsilon \left(x-1\right)z^{\prime }+\epsilon {\displaystyle {\int }_x^1{z}^{\prime }\text{d}x}-b\left(x-1\right)z\left(x\right)-b{\displaystyle {\int }_x^{1}z\text{d}x}+c{\displaystyle {\int }_x^1\left(x-1\right)z\text{d}x}=0,$

$\epsilon \left(x-1\right)z^{\prime }-b\left(x-1\right)z+\epsilon z\left(1\right)-\epsilon z\left(x\right)-b{\displaystyle {\int }_x^{1}z\text{d}x}+c{\displaystyle {\int }_x^1\left(x-1\right)z\text{d}x}=0,$

当 $x\le \frac{1}{2}$ 时，依然有 $\epsilon \left|z^{\prime }\right|\le CM$。

定理1得证。

4. 数值实验

本节基于奇性分离法求解奇异摄动两点边值问题(1)，取 $b=0$， $c=1$， $f=1-x$。先求解弱奇性的第三边值辅助问题(2)，再构造奇异函数得到原问题(1)的解。使用有限元方法计算时需将(2)化为一阶方程组

$\{\begin{array}{l}-\epsilon p^{\prime }+b{w}^{\prime }+cw=f\\ w^{\prime }-p=0\\ w\left(0\right)=0,\left(b{p}^{\prime }+cw\right)\left(1\right)=f(\; 1\; )\end{array}$

数值实验时不使用Shishkin网格，我们在边界层 $J_0=\left[0,\tau \right]$ 上取1个单元，正则区间 $J_1=\left[\tau ,1\right]$ 上取10个单元。计算采用二次平均间断有限元(ADG)方法。表1中列出的是正则区间和边界层上的有限元解误差。从表中可以看出，函数u和导数q在正则区间上的有限元解误差与 $\epsilon$ 无关，而边界层内， $\epsilon$ 越小，u的有限元解误差反而越小，且有 ${\hat{e}}_q\approx {\hat{e}}_u/\epsilon$，符合定理1的结论。

Table 1. The error of the flux of ADG under 10 + 1 meshgrid

表1. 10 + 1网格下二次平均间断元通量误差

基金项目

湖南省教育厅科学研究项目(18C0137)。