边界层问题广泛存在于流体力学、弹性力学、声光学、化学反应、最优控制等领域, 如高雷诺数的流体流动问题、热量传输问题、多孔介质中的渗流问题等等.这类问题描述的现象往往在局部区域具有奇异性, 对应的微分方程中含有反映介质性质的某个摄动小参数, 方程的解或其导数在边界层内变化非常剧烈.

利用数值方法求解边界层问题时, 若在一致网格上求解, 其数值解在边界层急剧变化区域可能会产生非物理振荡, 无法得到满意的结果. Lin$ \beta $ T等学者在非一致网格上研究各类边界层问题, 取得了一些研究成果, 如文献[1-9], 其中文献[1]针对一类边界层问题, 基于多种数值格式, 在非一致网格上得到了数值解收敛于真解的理论结果.文献[2-4]在非一致网格上研究了带两个小参数的奇异摄动问题, 得到了较好的理论结果.文献[5]在自适应移动网格上求解一类二阶双曲型奇异摄动问题, 改进了一致网格上的计算结果.

Bakhvalov-Shishkin网格^[10-14](简记为B-S网格)是能有效求解边界层问题的一类非一致层适应网格, 该类网格具有在边界层局部加密的特性.文献[11-13]在B-S网格上研究几类奇异摄动边界层问题, 得到数值解一致收敛于真解的结论.文献[14]在该类网格上用多尺度有限元逼近法来实现边界层问题的高效逼近.在B-S网格的构造中, 需要选择网格过渡点来确定网格分布.对于网格过渡点中的参数选择, 大部分文献都是人为选取的参数.而不同的网格参数会使得网格点发生变化, 从而影响求解的精度.因此, 如何选取网格过渡点中的参数, 使得数值解具有最好的精度, 是非常有意义的一项工作.

差分进化算法(DE)是Storn和Price^[15]首次提出的一类基于群体的全局优化算法, 目前广泛应用于数据挖掘、模式识别、人工神经网络等领域.文献[16]提出了一种基于泛化反向学习的多目标约束差分进化算法.文献[17]采用改进的差分进化算法, 研究了云制造资源的优化组合问题.文献[18]将差分进化算法与谱配点法结合, 有效求解了一类奇异摄动问题.

本文研究如下的边界层问题

(1.1) $ \begin{equation} \left\{ \begin{array}{ll} -\varepsilon u''(x)+p(x)u'(x)+q(x)u(x) = f(x), \; x \in (0, 1), \\ u(0) = u_0, u(1) = u_1, \end{array} \right. \end{equation} $

其中扩散项系数$ \varepsilon $是摄动小参数, 满足$ 0<\varepsilon \ll 1 $, 对流项系数为$ p(x) $.当$ \varepsilon $极小时, 该问题为稳态的对流占优问题, 可以描述稳态的热传导、粒子扩散、化学反应等现象.假设函数$ p(x), q(x), f(x) $满足充分光滑, 且$ q(x)\geq 0 $.

当对流项系数$ p(x)> 0 $时, 设存在常数$ \beta_1 $使得$ p(x) \geq \beta_1>0 $, 此时问题的解在右边界有厚度为$ O(\varepsilon) $的边界层.而当$ p(x)< 0 $时, 设存在常数$ \beta_2 $使得$ \mid p(x) \mid \geq \beta_2>0 $, 此时问题的解在左边界有厚度为$ O(\varepsilon) $的边界层.我们将对含左边界层和右边界层的这两类问题构造相应的Bakhvalov-Shishkin网格, 采用差分进化算法对网格参数进行优化, 获得具有最优精度的数值解.

我们在求解区间$ [0, 1] $上考虑模型问题(1.1).在$ [0, 1] $上构造网格点$ \{x_i\}(i = 0, 1, \cdots , N) $将问题(1.1)进行离散, 其中$ x_0 = 0, x_N = 1 $, 网格步长$ h_i = x_i-x_{i-1} $.记$ u_i^{N} $为$ x_i $处对应的数值解, 并记

$ D^+u_i^{N} = \frac{u_{i+1}^{N}-u_{i}^{N}}{h_{i+1}}, \; D^-u_i^{N} = \frac{u_{i}^{N}-u_{i-1}^{N}}{h_{i}}, $

$ {\hbar}_{i} = \frac{h_{i}+h_{i+1}}{2}, \; \delta^2u_{i}^{N} = \frac{D^+u_i^{N}-D^-u_i^{N}}{\hbar_i}. $

设网格剖分数$ N $为偶数, 当问题的解具有右边界层时, 定义网格过渡点$ \tau_1 = \frac{2\varepsilon }{\beta_1}\ln N $, 其中常数$ \beta_1 $满足$ p(x) \geq \beta_1>0 $. B-S网格的构造为将求解区间[0, 1]分解为两个子区间$ [0, 1-\tau_1] $和$ [1-\tau_1, 1] $, 再将子区间$ [0, 1-\tau_1] $作$ N/2 $等分, 将子区间$ [1-\tau_1, 1] $也分为$ N/2 $分, 但要使得$ e^{-\beta_1(1-x)/{2\varepsilon }} $等分布在网格节点$ \{x_i\}(i = N/2, N/2+1, \cdots , N) $上, 得到如下网格分布

(2.1) $ \begin{equation} x_i = \left\{ \begin{array}{ll} { } \frac{2i}{N}(1-\frac{2\varepsilon }{\beta_1}\ln N), & i = 0, 1, \cdots , N/2, \\ { } 1+\frac{2\varepsilon }{\beta_1}\ln [1-2(1-\frac{1}{N})(1-\frac{i}{N})], \ &i = N/2+1, \cdots , N. \end{array} \right. \end{equation} $

在该网格上构造问题(1.1)的差分格式为

(2.2) $ \begin{equation} \left\{ \begin{array}{ll} { } -\varepsilon \delta^2u_{i}^{N}+p_i D^-u_i^{N}+q_iu_i = f_i, \; 1\leq i\leq N-1, \\ { } u_0^N = u_0, u_N^N = u_1.\end{array} \right. \end{equation} $

当问题的解具有左边界层时, 定义网格过渡点$ \tau_2 = \frac{2\varepsilon }{\beta_2}\ln N $, 其中常数$ \beta_2 $满足$ \mid p(x) \mid \geq \beta_2>0 $.将区间[0, 1]分解为两个子区间$ [0, \tau_2] $和$ [\tau_2, 1] $, 将子区间$ [\tau_2, 1] $作$ N/2 $等分, 将子区间$ [0, \tau_2] $也分为$ N/2 $分, 但要使得$ e^{-\beta_2 x/{2\varepsilon }} $等分布在网格节点$ \{x_i\}(i = 0, 1, \cdots , N/2) $上, 得到如下B-S网格分布

(2.3) $ \begin{equation} x_i = \left\{ \begin{array}{ll} { } -\frac{2\varepsilon}{\beta_2}\ln(1-2(1-\frac{1}{N})\frac{i}{N}), & i = 0, 1, \cdots , N/2, \\ { } 1-(1-\frac{2\varepsilon }{\beta_2}\ln N)\frac{2(N-i)}{N}, \ & i = N/2+1, \cdots , N. \end{array} \right. \end{equation} $

在该网格上构造问题(1.1)的差分格式为

(2.4) $ \begin{equation} \left\{ \begin{array}{ll} -\varepsilon \delta^2u_{i}^{N}+p_i D^+u_i^{N}+q_iu_i = f_i, \; 1\leq i\leq N-1, \\ u_0^N = u_0, u_N^N = u_1.\end{array} \right. \end{equation} $

显然, 当小参数$ \varepsilon $和网格剖分数$ N $确定时, 两种情况下的网格过渡点分别由参数$ \beta_1 $和$ \beta_2 $决定.选取不同的参数会使得网格点发生变化, 从而影响求解的精度.

引理2.1^[13]  当$ p(x)> 0 $时, 设$ u(x_i) $为问题(1.1)在网格(2.1)上的精确解, $ u_i^N $为在网格(2.1)上求解差分格式(2.2)的数值解, 则$ u_i^N $关于$ \varepsilon $一致收敛于$ u(x_i) $, 且有

$ \max\limits_{i}|u(x_i)-u_i^N|\leq CN^{-1}. $

引理2.2^[12]  当$ p(x)< 0 $时, 设$ \tilde{u}(x_i) $为问题(1.1)在网格(2.3)上的精确解, $ \tilde{u}_i^N $为在网格(2.3)上求解差分格式(2.4)的数值解, 则$ \tilde{u}_i^N $关于$ \varepsilon $一致收敛于$ \tilde{u}(x_i) $, 且有

$ \max\limits_{i}|\tilde{u}(x_i)-\tilde{u}_i^N|\leq CN^{-1}. $

由引理可知, 在B-S网格上求解边界层问题(1.1)的差分格式时, 不论小参数$ \varepsilon $如何取值, 数值解均一阶收敛于真解.而选取不同的网格过渡点构造的B-S网格不会影响数值解的收敛阶.

将B-S网格(2.1)和(2.3)中的网格参数$ \beta_1 $和$ \beta_2 $统一记为$ \beta $.若$ u^N $是网格剖分为N个子区间时在B-S网格上求解边界层问题(1.1)的数值解, 它可视为网格参数$ \beta $的函数$ u^N(\beta) $.由于很多实际问题中很难得到边界层问题的精确解, 因此使用双重网格的数值解来估计数值解误差, 记

(3.1) $ \begin{equation} E(\beta) = \parallel u^N(\beta)-u^{2N}(\beta) \parallel_\infty, \end{equation} $

其中$ u^{2N}(\beta) $为将网格剖分为$ 2N $个子区间时的数值解, 此时的网格点由网格剖分为$ N $个子区间时的B-S网格点及两个相邻网格点的中点组成.

为了使误差尽可能的小, 构造如下的目标函数

(3.2) $ \begin{equation} \min E(\beta) = \parallel u^N(\beta)-u^{2N}(\beta) \parallel_\infty. \end{equation} $

以误差范数(3.2)为目标函数寻找最优的网格参数$ \beta $.

差分进化算法将参数$ \beta $初始化为随机种群, 通过把父代种群中任意两个个体的向量差加权后与第三个个体求和来产生新个体, 然后将新个体与当代种群中相应的个体相比较, 如果新个体的适应度优于当前个体的适应度, 则在下一代中用新个体取代旧个体, 否则仍保存旧个体.通过不断地进化, 保留优良个体, 实现种群的进化, 引导搜索向最优解逼近.

下面给出采用差分进化算法结合差分格式在B-S网格上求解边界层问题(1.1)的算法步骤.

第一步  对参数$ \beta $进行种群随机初始化, 记$ \{\beta_i^t\}, i = 1, 2, \cdots , NP $为第$ t $代种群的个体, 其中$ NP $为种群规模.此时$ t = 0 $.

第二步  将$ \{\beta_i^t\}, i = 1, 2, \cdots , NP $代入(2.1)或(2.3)式计算对应的B-S网格.对B-S网格上的差分格式进行求解, 得到B-S网格上$ N+1 $个网格点上的数值解$ u^N $.在原B-S网格上添加相邻网格点的中点, 得到$ 2N+1 $个网格点, 再次求解差分格式, 得到数值解$ u^{2N} $, 进而计算出种群中每个个体对应的目标函数(3.2)的适应值.

第三步(变异操作)  对每个个体$ \beta_i^t $, 令新的变异个体$ \nu_i^t = \beta_{\gamma_1}^t+F(\beta_{\gamma_2}^t-\beta_{\gamma_3}^t) $, 其中$ \gamma_1, \gamma_2, \gamma_3 $是随机下标, 且$ \gamma_1 \neq \gamma_2 \neq \gamma_3 \neq i $, $ F $为变异因子.

第四步(交叉操作)  对变异种群$ \{\nu_i^t\} $和父代种群$ \{\beta_i^t\} $通过交叉操作得到新的种群$ \{z_i^t\} $:

(3.3) $ \begin{equation} z_i^t = \left\{ \begin{array}{l} \nu_i^t, rand\leq CR, \\ \beta_i^t, rand>CR, \end{array} \right. \end{equation} $

其中$ rand $是区间$ (0, 1) $的随机数, $ CR $为范围在$ (0, 1) $间的交叉因子.

第五步  对个体$ z_i^t $和父代种群中个体$ \beta_i^t $对应的误差范数(3.1)进行比较, 选择函数值较小的作为下一代的个体

(3.4) $ \begin{equation} \beta_i^{t+1} = \left\{ \begin{array}{l} z_i^t, E(z_i^t)<E(\beta_i^t), \\ \beta_i^t, E(z_i^t)\geq E(\beta_i^t), \end{array} \right. \end{equation} $

第六步  若$ t $小于设定的最大迭代次数$ tmax $, 则令$ t = t+1 $, 转第二步.否则, 停止迭代, 将上一步得到的最优个体作为网格参数$ \beta $的值, 相应的误差范数作为最优误差范数.

根据上节给出的算法步骤, 采用差分进化算法在B-S网格上对问题(1.1)进行数值实验.实验时取种群规模$ NP $为50, 最大迭代次数$ tmax $为50次, 变异因子$ F $为0.5, 交叉因子$ CR $为0.9.通过以下三个算例来进行数值模拟.

算例1  考虑具有右边界层的问题

$ \begin{eqnarray*} \left\{ \begin{array}{ll} { } -\varepsilon u''(x)+\frac{2}{x+2}u'(x)+\frac{1}{x+2}u(x) = f(x), \; x \in (0, 1), \\ { } u(0) = (\frac{2}{3})^{\frac{1}{\varepsilon}+1}, u(1) = 1, \end{array} \right. \end{eqnarray*} $

这里右端项

$ f(x) = \frac{1}{3}{(\frac{x+2}{3})}^{1/\varepsilon}(\frac{x+3+1/\varepsilon}{x+2}). $

此例中$ p(x) = \frac{2}{x+2}\geq \frac{2}{3}>0 $, 故网格参数$ \beta_1 $需要满足$ 0<\beta_1\leq \frac{2}{3} $.在实验时取参数$ \beta_1 $的搜索范围为$ [\frac{1}{100}, \frac{2}{3}] $.

算例2  考虑具有右边界层的问题

$ \begin{eqnarray*} \left\{ \begin{array}{ll} -\varepsilon u''(x)+u'(x)+xu(x) = f(x), \; x \in (0, 1), \\ u(0) = 0, u(1) = 1, \end{array} \right. \end{eqnarray*} $

这里右端项

$ f(x) = \frac{x(e^{-(1-x)/\varepsilon}-e^{-1/\varepsilon})}{1-e^{-1/\varepsilon}}. $

此例中$ p(x) = 1>0 $, 故网格参数$ \beta_1 $需要满足$ 0<\beta_1\leq 1 $.在实验时取参数$ \beta_1 $的搜索范围为$ [0.01, 1] $.

算例3  考虑具有左边界层的问题

$ \begin{eqnarray*} \left\{ \begin{array}{ll} -\varepsilon u''(x)-u'(x)+2u(x) = e^{x-1}, \; x \in (0, 1), \\ u(0) = 0, u(1) = 0. \end{array} \right. \end{eqnarray*} $

此例中$ p(x) = -1<0 $, 故网格参数$ \beta_2 $需要满足$ 0<\beta_2\leq \mid p(x)\mid = 1 $.在实验时取参数$ \beta_2 $的搜索范围为$ [0.01, 1] $.

取网格剖分数$ N = 32, 64, 128, 256 $, 小参数$ \varepsilon = 10^{-2}, 10^{-3}, 10^{-4}, 10^{-5}, 10^{-6} $, 对三个算例进行数值模拟. 表 1–表 3列出了以(3.2)式为目标函数时, 采用差分进化算法获取的最优网格参数及目标函数值.可见在目标函数最优时, 该方法具有收敛性.

当问题的解具有右边界层时, 网格过渡点为$ \tau_1 = \frac{2\varepsilon }{\beta_1}\ln N $.当问题的解具有左边界层时, 网格过渡点为$ \tau_2 = \frac{2\varepsilon }{\beta_2}\ln N $.由网格过渡点的定义可知, 当固定小参数$ \varepsilon $不变时, 网格参数$ \beta_1 $或$ \beta_2 $应随网格剖分数$ N $的增加而增大.从表 1–$ \! $表 3数据来看, 当摄动参数$ \varepsilon $取较小值$ 10^{-3}, 10^{-4}, 10^{-5}, 10^{-6} $时, 网格参数确实存在这样的变化规律.

图 1作出了算例1和算例3中$ N = 64, \varepsilon = 10^{-5} $时, 目标函数随迭代次数增加的变化.可见经过较少次数的迭代, 目标函数可接近最优值, 说明差分进化算法求解此问题具有较快的收敛速度. 图 1(a)显示, 在$ N = 64, \varepsilon = 10^{-5} $时, 约经过10次的迭代, 算例1的目标函数接近最优值.取迭代次数为10次和50次, 差分进化算法计算算例1的计算时间分别为0.796秒和3.853秒.可见, 利用差分进化算法通过迭代寻找最优网格参数具有较高的计算效率.

取摄动小参数$ \varepsilon = 10^{-2} $和$ \varepsilon = 10^{-5} $, 表 4–表 6分别列出了3个算例中, 采用差分进化算法选取网格参数$ \beta $和取参数$ \beta $为固定值(包括搜索上限)时的误差.表中数据反映出, 采用差分进化算法确定的最优参数$ \beta $对应的误差明显优于$ \beta $取其它固定值时的误差, 且具有收敛性.比如表 4中, 当小参数$ \varepsilon $极小, 取$ \varepsilon = 10^{-5}, N = 128 $时, 若网格参数取搜索上限$ \beta_1 = 2/3 $, 误差数量级为E+0, 而采用差分进化算法求解得到的误差数量级为E-03. 表 5和表 6中, 网格参数取0.1和取差分进化算法确定的最优参数相比, 误差区别也很明显.

我们将网格参数$ \beta $取差分进化算法中最优解与其他固定值时的误差范数目标函数(3.2)图像作比较, 发现误差的区别主要体现在边界层.算例1和算例2的解在$ x = 1 $附近的边界层变化剧烈, 而算例3的解在$ x = 0 $附近的边界层变化剧烈.为了便于观察边界层误差情况, 取$ N = 64, \; \varepsilon = 10^{-5} $, 图 2和图 3分别作出了算例1和算例2中网格参数$ \beta_1 $取最优解与固定值时右边界层(分别为区间[0.999, 1]和[0.9995, 1])的误差曲线, 图 4作出了算例3中网格参数$ \beta_2 $取最优解与固定值时左边界层(区间[0, 0.0004])的误差曲线.误差曲线反映出, 在边界层, 采用差分进化算法得到的数值解精度比人为选取网格参数要高出许多.

因此, 与人为选择B-S网格中的网格参数相比, 利用差分进化算法优化选择网格参数进而求解边界层问题能得到更好的数值结果.该方法经过少量次数的迭代可得到最优网格参数, 在计算精度、收敛性方面均有优势.

本文采用差分进化算法优化选择Bakhvalov-Shishkin网格中参数, 提出了求解具有左、右边界层的奇异摄动问题的差分进化算法.实验结果表明, 与选择固定的网格参数相比, 采用差分进化算法计算能极大提高计算精度, 特别是边界层的计算精度, 并且数值结果具有收敛性.本文的方法也可以推广到Bakhvalov-Shishkin网格上其它问题的计算或高维问题的计算中.