拟线性奇异摄动微分方程的自适应差分方法

Adaptive Difference Method for Quasilinear Singularly Perturbed Differential Equations

摘要: 利用简单迎风差分格式和自适应网格方法相结合求解拟线性奇异摄动方程两点边值问题，通过等分布弧长控制函数而产生自适应网格，得到自适应差分格式的关于摄动参数ε一致的一阶误差估计，最后数值算例验证数值解的一阶一致收敛性，表明该方法的可靠性与准确性。

Abstract: A simple upwind difference scheme and an adaptive grid method are combined to solve the quasilinear singularly perturbed two-point boundary value problem. The adaptive mesh is generated by equally distributing the arc-length control function. The first-order error estimate of the adaptive difference scheme uniform with respect to the perturbation parameter ε is obtained. Finally, numerical examples are given to verify the uniform first-order convergence and show the reliability and accuracy of the method.

1. 引言

奇异摄动问题在流体力学和地球物理学等领域中应用广泛。考虑拟线性守恒型奇异摄动方程的两点边值问题：

$Tu\left(x\right):=-\epsilon u^{″}\left(x\right)+b{\left(x,u\left(x\right)\right)}^{\prime }=f\left(x\right),x\in \left(0,1\right)$ (1.1)

$u\left(0\right)=u\left(1\right)=0$ (1.2)

其中 $\epsilon$ 是一个很小的摄动参数， $0<\epsilon \ll 1$。假设 $b\in C^1\left(\left[0,1\right]\ast ℝ\right)$， $f\in C\left[0,1\right]$，且存在常数 $\beta$ 和 $\bar{\beta }$，使得 $0<\beta <b_u\left(x,u\right)\le \bar{\beta }$， $\forall x\in \left[0,1\right]$， $\forall u\in ℝ$ 成立。则问题(1.1) (1.2)存在唯一解 $u\in C^2\left[0,1\right]$。解在边界 $x=1$ 处有一个变化很剧烈的指数边界层，用许多数值方法都不能有效逼近问题的准确解。因此我们使用自适应网格方法，即采用等分布控制函数来自适应地生成网格，进而在边界层区域处放置相对较多的网格点，来保证计算出来的解是更好的近似解(见 [1] )。

文献 [2] 通过等分布数值解的弧长实现自适应移动网格，证明了问题(1.1)在这种自适应移动网格下的简单迎风格式求解结果是一阶收敛的。文献 [3] 讨论了在Sakhvalov-Shishkin网格上求解比问题(1.1)更一般的拟线性奇异摄动方程，建立带权的混合差分格式，估计解误差和导数误差，最终得到二阶收敛的解和导数及其误差估计。文献 [4] 采用一阶和二阶差分格式离散了问题(1.1)，得到依赖于数值解差商的后验误差估计。文献 [5] 在此基础上提出了一种可调节控制函数的自适应方法，数值实验结果验证了方法的可行性。对于非守恒型线性奇异摄动方程，文献 [6] 中利用误差校正方法将高阶不稳定格式和标准迎风格式差分格式联合起来，改进了其精度。文献 [7] 又对文献 [6] 进行改进，利用Richardson外推技术得到了一个更精确的逼近解。文献 [8] 讨论了自适应网格方法求解齐次线性非守恒型奇异摄动方程。

考虑拟线性守恒型奇异摄动方程的两点边值问题(1.1)，鉴于自适应差分方法克服了层适应差分方法需要预知准确解情况的缺点，本文通过简单迎风格式与等分布数值解弧长相结合求得数值解，证明该方法是一阶收敛的且关于摄动参数是一致的，数值算例验证一阶一致收敛性，表明该方法的可靠性与准确性。

2.自适应差分方法及其误差估计

2.1.自适应简单迎风格式

为了离散求解拟线性两点边值问题(1.1)，设

$Av=-\epsilon v^{\prime }\left(x\right)+b\left(x,v\left(x\right)\right),v\in C^2\left[0,1\right]$ (2.1)

则(1.1)可写成

${\left(Au\right)}^{\prime }=-\epsilon u^{″}+b{\left(x,u\left(x\right)\right)}^{\prime }=f$ (2.2)

定义任意网格函数 $\left\{v_i\right\}$ 的离散算子 $A^N$，将(2.1)离散化，可以得到

$A^N{v}_i=-\epsilon D^{+}{v}_i+b\left(x_i,v_i\right)$ (2.3)

其中 $i=1,2,3,\cdots ,N-1$ 且 $0=x_0<x_1<\cdots <x_N=1$

将原方程(1.1)离散化可以得到

$T^N{u}_i^N:=D^{-}\left(A^N{u}_i^N\right)=f_i,i=1,2,\cdots ,N-1$ (2.4)

$u_0^N=u_N^N=0$

其中 $\left\{u_i^N\right\}$ 是在网格 $\left\{x_i\right\}$ 上计算的数值解。

在自适应网格时，除了 $x_0$ 和 $x_N$ 以外的N − 1个节点，每一个点处都可以列出两个式子：

i) 此节点上的差分格式；

ii) 此节点左右两个弧长相等的式子；

联立这 $2N-2$ 个式子，就可以得到自适应的 $x_i$ 和 $u_i^N$。

2.2. 线性守恒型两点边值问题的稳定性

考虑(1.1)的线性情况

$Lu=-\epsilon u^{″}\left(x\right)+{\left(pu\left(x\right)\right)}^{\prime }=f\left(x\right)$ (2.5)

$u\left(0\right)=u\left(1\right)=0$

其中

(i) $0<\beta \le p\left(x\right)\le \bar{\beta }$ (2.6a)

(ii) $f\in C\left[0,1\right]$ (2.6b)

引理2.1设(2.5)中的 $f\left(x\right)=F^{\prime }\left(x\right)\in C\left[0,1\right]$，则方程中存在一个唯一解 $u\in C^2\left[0,1\right]$，使得

${\Vert u\Vert }_{\infty }\le \frac{2}{\beta }{\Vert Lu\Vert }_{\text{*}}$

其中， ${\Vert \text{v}\Vert }_{\text{*}}$ 表示 ${\Vert \text{v}\Vert }_{\text{*}}={\min }_{\text{<msup> <mi> V </mi> <mo> ′ </mo> </msup>}=\text{v}}{\Vert \text{V}\Vert }_{\infty }$。

证明 (2.5)两边同时对 积分可得

$u^{\prime }-\frac{pu}{\epsilon }=-\frac{F\left(x\right)+C_1}{\epsilon }$

可以解出

$u\left(x\right)=C{\text{e}}^{{\displaystyle {\int }_x^1\frac{p\left(t\right)}{\epsilon }\text{d}t}}+{\text{e}}^{-{\displaystyle {\int }_x^1\frac{p\left(t\right)}{\epsilon }\text{d}t}}{\displaystyle {\int }_x^1\frac{F\left(s\right)+C_1}{\epsilon }{\text{e}}^{{\displaystyle {\int }_s^1\frac{p\left(t\right)}{\epsilon }\text{d}t}}\text{d}s}$

由于 $u\left(1\right)=0$，可以解出 $\text{C}=0$，则有

$u\left(x\right)={\displaystyle {\int }_x^1\frac{F\left(s\right)}{\epsilon }{\text{e}}^{-{\displaystyle {\int }_x^s\frac{p\left(t\right)}{\epsilon }\text{d}t}}\text{d}s}+C_1{\displaystyle {\int }_x^1\frac{1}{\epsilon }{\text{e}}^{-{\displaystyle {\int }_x^s\frac{p\left(t\right)}{\epsilon }\text{d}t}}\text{d}s}$ (2.7)

设

$W\left(x\right)={\displaystyle {\int }_x^1\frac{F\left(s\right)}{\epsilon }{\text{e}}^{-{\displaystyle {\int }_x^s\frac{p\left(t\right)}{\epsilon }\text{d}t}}\text{d}s}$ (2.8)

$V\left(x\right)={\displaystyle {\int }_x^1\frac{1}{\epsilon }{\text{e}}^{-{\displaystyle {\int }_x^s\frac{p\left(t\right)}{\epsilon }\text{d}t}}\text{d}s}$ (2.9)

则有

$u\left(x\right)=W\left(x\right)+C_{1}V(\; x\; )$

将 $u\left(0\right)=0$ 代入，可以解出 $C_1=-\frac{W\left(0\right)}{V\left(0\right)}$，则有

$u\left(x\right)=W\left(x\right)-\frac{W\left(0\right)}{V\left(0\right)}V(\; x\; )$

由(2.7)和(2.8)可得

$\left|w\left(x\right)\right|\le {\Vert F\left(x\right)\Vert }_{\infty }V\left(x\right)$ (2.10)

$0\le V\left(x\right)\le {\displaystyle {\int }_x^1\frac{1}{\epsilon }{\text{e}}^{-{\displaystyle {\int }_x^s\frac{\beta }{2}\text{d}t}}\text{d}s}\le {\displaystyle {\int }_x^1\frac{1}{\beta }{\text{e}}^{-\frac{\beta }{\epsilon }\left(s-x\right)}\text{d}\frac{\left(s-x\right)\beta }{\epsilon }}\le \frac{1}{\beta }$

由(2.10)有

$\left|W\left(0\right)\right|\le {\Vert F\left(x\right)\Vert }_{\infty }V\left(0\right)$ (2.11)

则有

$\begin{array}{c}\left|u\left(x\right)\right|\le \left|W\left(x\right)\right|+\frac{\left|W\left(0\right)\right|}{V\left(0\right)}V\left(x\right)\le \left|W\left(x\right)\right|+{\Vert F\left(x\right)\Vert }_{\infty }V\left(x\right)\\ \le 2{\Vert F\left(x\right)\Vert }_{\infty }V\left(x\right)\le \frac{2}{\beta }{\Vert F\left(x\right)\Vert }_{\infty }\end{array}$ (2.12)

由(2.12)和(2.5)可得引理2.1，证毕。

2.3. 离散格林函数与弧长有界性

定义算子 $T^N$ 的标准线性化 ${\tilde{T}}^N$，其中 ${\tilde{A}}^N{u}_i^N=\epsilon D^{-}{u}_i^N+p^i{u}_i^N$， $p^i={\displaystyle {\int }_{s=0}^1{b}_u\left(x_i,s{y}_i\right)\text{d}s}$。则对于 $i=1,2,\cdots ,N-1$ 有

${\tilde{T}}^N{u}_i^N={\tilde{T}}^N\left(u_i^N-\stackrel{\to }{0}\right)=T^N{u}_i^N-T^N\stackrel{\to }{0}=f_i-{\left(T^N\stackrel{\to }{0}\right)}_i,u_0^N=u_N^N=0$ (2.13)

由(1.1)满足 $0<\beta <b_u\left(x,u\right)\le \bar{\beta }$ 知 $p^i\ge \beta$，容易验证 ${\tilde{T}}^N$ 的系数矩阵是一个M矩阵。

现在定义关于点 $x_j$ 的离散格林函数 $G\left(x_i,x_j\right)$

${\tilde{T}}^{N}G\left(x_i,x_j\right)=\frac{{\delta }_{ij}}{h_{i+1}},G\left(0,x_j\right)=G\left(1,x_j\right)=0,i=1,2,\cdots ,N-1$ (2.14)

其中 ${\tilde{T}}^N$ 作用于变量i，克罗内克函数 ${\delta }_{ij}$ 的值当 $i=j$ 时得1，当 $i\ne j$ 时得0。则对于每个i，我们都有

$u_i^N={\displaystyle \underset{j=1}{\overset{N-1}{\sum }}{h}_{j+1}G\left(x_i,x_j\right)\left[f_i-{\left(T^N\stackrel{\to }{0}\right)}_j\right]}$ (2.15)

现在引入格林函数的一些性质(见 [2] )：

引理2.2对于 $i=1,2,\cdots ,N$ 且 $j=1,2,\cdots ,N-1$ 有

$0\le G\left(x_i,x_j\right)\le \frac{1}{\beta }$

引理2.3对于每个 $j\in \left\{1,2,\cdots ,N-1\right\}$，有

$\underset{i=1}{\overset{N-1}{{\displaystyle \sum }}}\left|G\left(x_i,x_j\right)-G\left(x_{i-1},x_j\right)\right|\le \frac{2}{\beta }$

引理2.4令 $\left\{u_i^N\right\}$ 为任意网格 $\left\{x_i\right\}$ 上关于(2.4)的解， $L^N$ 为 $\left\{u_i^N\right\}$ 函数曲线的总弧长，则有

$1\le L^N\le C_2$

其中， $C_2=1+\frac{2}{\beta }{\Vert f-T^N\stackrel{\to }{0}\Vert }_{\infty }$。

证明 首先， $L^N\ge {\displaystyle {\sum }_{i=1}^N{h}_i}=1$. 其次，由(2.15)和引理2.3有

$L^N={\displaystyle \underset{i=1}{\overset{N}{\sum }}{h}_i\sqrt{1+{\left(D^{-}{u}_i^N\right)}^2}}\le {\displaystyle \underset{i=1}{\overset{N}{\sum }}{h}_i\left[1+\left|D^{-}{u}_i^N\right|\right]}=1+{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left|u_i^N-u_{i-1}^N\right|}$ (2.16)

因此可得

$\begin{array}{c}{L}^N\le 1+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}{h}_{j+1}\left|f_j-{\left(T^N\stackrel{\to }{0}\right)}_j\right|\cdot \left|G\left(x_i,x_j\right)-G\left(x_{i-1},x_j\right)\right|}}\\ \le 1+{\Vert f-\left(T^N\stackrel{\to }{0}\right)\Vert }_{\infty }{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}{h}_{j+1}}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left|G\left(x_i,x_j\right)-G\left(x_{i-1},x_j\right)\right|}\\ \le 1+\frac{2}{\beta }{\Vert f-T^N\stackrel{\to }{0}\Vert }_{\infty }\end{array}$ (2.17)

证毕。

2.4. 后验误差估计和等分布弧长及误差估计

定理2.5 (后验误差估计)设 $f\left(x\right)\in C^1\left[0,1\right]$，则存在一个常数C，使后验界满足

${\Vert u^N-u\Vert }_{\infty }\le \frac{2}{\beta }\underset{i}{\max }\left\{\bar{\beta }\left|u_i^N-u_i^{N-1}\right|+C{h}_i\right\}$

证明 因为 $T{u}^N-Tu=L\left(u^N-u\right)$，故由引理2.1可知

${\Vert u^N-u\Vert }_{\infty }\le \frac{2}{\beta }{\Vert T{u}^N-Tu\Vert }_{*}$ (2.18)

因此，只需证

$\frac{2}{\beta }{\Vert T{u}^N-Tu\Vert }_{*}\le \frac{2}{\beta }\underset{i}{\max }\left\{\bar{\beta }\left|u_i^N-u_i^{N-1}\right|+C{h}_i\right\}$ (2.19)

即可。

又因为

${\Vert T{u}^N-Tu\Vert }_{*}={\Vert T{u}^N-f\Vert }_{*}={\Vert {\left(A{u}^N-F-a\right)}^{\prime }\Vert }_{*}$ (2.20)

其中a为常数。

令

$F_i^N={\displaystyle \underset{j=1}{\overset{i}{\sum }}{h}_j{f}_j}={\displaystyle \underset{j=1}{\overset{i}{\sum }}\left(A^N{u}_j^N-A^N{u}_{j-1}^N\right)}=A^N{u}_i^N-A^N{u}_0^N$ (2.21)

$a=A^N{u}_0^N=A^N{u}_i^N-F_i^N$

则有

${\Vert T{u}^N-Tu\Vert }_{*}={\Vert A{u}^N-A^N{u}_i^N+F_i^N-F\Vert }_{*}\le {\Vert A{u}^N-A^N{u}_i^N+F_i^N-F\Vert }_{\infty }$

令 $\bar{\eta }=A{u}^N-A^N{u}_i^N,\tilde{\eta }=F_i^N-F,\eta =\bar{\eta }+\tilde{\eta }$，对于 $x\in \left(x_{i-1},x_i\right)$，有

$\left|\tilde{\eta }\right|\le C_1{h}_i$ (2.22)

还有

$\begin{array}{c}\underset{x_{i-1}\le x\le x_i}{\sup }\left|\bar{\eta }\right|=\underset{x_{i-1}\le x\le x_i}{\sup }\left|\left(-\epsilon D^{+}{u}_i+b\left(x,u_i^N\right)\right)-\left(-\epsilon D^{+}{u}_i+b\left(x_i,u_i^N\right)\right)\right|\\ =\underset{x_{i-1}\le x\le x_i}{\sup }\left|b\left(x,u_i^N\right)-b\left(x_i,u_i^N\right)\right|=\underset{x_{i-1}\le x\le x_i}{\sup }\left|{\displaystyle {\int }_x^{x_i}\frac{\text{d}\left[b\left(x,u^N\left(x\right)\right)\right]}{\text{d}x}\text{d}x}\right|\\ \le \underset{x_{i-1}\le x\le x_i}{\sup }{\displaystyle {\int }_x^{x_i}\left(\left|b_x\left(x,u^N\left(x\right)\right)\right|+\left|b_u\left(x,u{\left(x\right)}^N\right)\right|\left|D^{-}{u}_i\right|\right)\text{d}x}\\ \le h_i\underset{x_{i-1}\le x\le x_i}{\sup }\left|b_x\left(x,u^N\left(x\right)\right)\right|+\bar{\beta }\left|u_i^N-u_{i-1}^N\right|\\ \le C_2{h}_i+\bar{\beta }\left|u_i^N-u_{i-1}^N\right|\end{array}$ (2.23)

由(2.22)和(2.23)就有

${\Vert \eta \Vert }_{\infty }\le \underset{i}{\max }\left\{\bar{\beta }\left|u_i^N-u_{i-1}^N\right|+C{h}_i\right\}$

则

$\frac{2}{\beta }{\Vert T{u}^N-Tu\Vert }_{*}\le \frac{2}{\beta }\underset{i}{\max }\left\{\bar{\beta }\left|u_i^N-u_{i-1}^N\right|+C{h}_i\right\}$

证毕。

又因为

$\begin{array}{l}\frac{2}{\beta }\underset{i}{\max }\left\{\bar{\beta }\left|u_i^N-u_{i-1}^N\right|+C{h}_i\right\}\\ \le C\underset{i}{\max }\sqrt{{\left(u_i^N-u_{i-1}^N\right)}^2+h_i^2}\\ =C\underset{i}{\max }{h}_i\sqrt{1+{\left(D^{-}{u}_i^N\right)}^2}\end{array}$ (2.24)

故由引理2.5和(2.24)有

${\Vert u^N-u\Vert }_{\infty }\le C\underset{i}{\max }{h}_i\sqrt{1+{\left(D^{-}{u}_i^N\right)}^2}$ (2.25)

其中不等式右侧为分段线性函数曲线从 $\left(x_{i-1},u_{i-1}^N\right)$ 到 $\left(x_i,u_i^N\right)$ 的最长线段长度，为了得到更精确的近似解，应使其值尽可能小。由此引入分段线性插值 $u^N$ 的等分布弧长控制函数

$M^N\left(x\right):=\sqrt{1+{\left({\left(u^N\right)}^{\prime }\left(x\right)\right)}^2}$ (2.26)

令

$M_i=\sqrt{1+{\left(D^{-}{u}_i^N\right)}^2},i=1,2,\cdots ,N$ (2.27)

用等分布控制函数来自适应生成网格，设 $h_i$ 为网格的宽度， $i=1,2,\cdots ,N$，则有

$h_i{M}_i=\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{h}_j{M}_j},i=1,2,\cdots ,N$ (2.28)

定理2.6 (误差估计)设 $\left\{u_i^N\right\}$ 是在满足(2.28)的自适应网格上通过方程(1.1)的离散化形式(2.4)计算而得的解，则对某常数C，有

${\Vert u-u^N\Vert }_{\infty }\le C{N}^{-1}$

证明 由 $l_i=\frac{L^N}{N},i=1,2,\cdots ,N$，其中局部弧长 $l_i=\sqrt{{\left(x_i-x_{i-1}\right)}^2+{\left(u_i{}^N-u_{i-1}^N\right)}^2}$，总弧长 $L^N={\displaystyle {\sum }_{i=1}^N{l}_i}$。又由引理2.4，有 $l_i=\frac{L^N}{N}\le \frac{C_2}{N}$，再结合(2.25)，就有 ${\Vert u-u^N\Vert }_{\infty }\le C{N}^{-1}$ 成立。证毕。

3. 数值算例

例1. 考虑拟线性奇异摄动问题

$\{\begin{array}{l}Tu\left(x\right):=-\epsilon u^{″}\left(x\right)-{\left(e^u\right)}^{\prime }+\frac{\pi }{2}\sin \frac{\pi x}{2}{e}^{2u}=0,x\in \left(0,1\right),\\ u\left(0\right)=u\left(1\right)=0.\end{array}$

此方程的准确解为 $u\left(x\right)=u_A\left(x\right)+O\left(\epsilon \right)$，其中

$u_A\left(x\right)=-\ln \left[\left(1+\cos \frac{\pi x}{2}\right)\left(1-\frac{1}{2}{e}^{-\frac{x}{2\epsilon }}\right)\right]$.

对于不同的 $\epsilon$ 和N，得到相应的误差 $e:=\Vert u_A-u^N\Vert$ 和收敛阶 $r:={\log }_2\frac{\max \left|u_i-u_i^N\right|}{\max \left|u_i-u_i^{2N}\right|}$，如表1。

例2. 考虑线性奇异摄动问题

$\{\begin{array}{l}Tu\left(x\right):=-\epsilon u^{″}\left(x\right)-{\left(u+\frac{x}{4}{u}^2\right)}^{\prime }=f\left(x,\epsilon \right),x\in \left(0,1\right)\\ u\left(0\right)=u\left(1\right)=0,\end{array}$

等式右边的 $f\left(x,\epsilon \right)$ 是由下面的准确解u所决定的：

$u\left(x\right)=\frac{\left(e-1\right)e^{-\frac{x}{\epsilon }}-e+e^{-\frac{1}{\epsilon }}}{1-e^{-\frac{1}{\epsilon }}}+e^x$.

通过表1和表2中的数值结果可以发现，收敛阶r随着N的增大在逐渐趋于1，所以选取适当的初值，在自适应网格上简单迎风差分格式求得的近似解和准确解之间的误差满足一阶收敛的估计，并且与 $\epsilon$ 的取值无关。

Table 1. Numerical results of Example 1 using the adaptive simple upwind difference method

表1. 用自适应简单迎风差分方法求解例1的数值结果

Table 2. Numerical results of Example 2 using the adaptive simple upwind difference method

表2. 用自适应简单迎风差分方法求解例2的数值结果

致谢

本文得到北京市大学生创新创业训练计划项目支持，感谢项目指导教师郑权教授的悉心指导。

参考文献