奇异摄动问题不仅在应用数学、物理学领域有广泛的应用(例如,对流扩散和反应扩散方程、著名的Navier-Stokes方程、半导体设备建模的漂移扩散方程、量子力学等),还广泛应用于与人们的生产生活密切相关的工程技术领域(例如,在开采石油、天然气的研究中所涉及的渗流问题,在天气预报中为描述大气运动的流体力学和热力学方程等)[1-2].一般情况下,这类问题的高阶导数项包含正的小参数,称为摄动参数.当该小参数无限趋向于零时,问题的解在某些很小的局部区域内变化非常剧烈,即解存在边界层或内点层.因而很多的问题很难求出精确解,尤其是非线性的问题.因此,伴随着计算机技术的高速发展,这类问题的数值解法越来越受到广大科研人员的关注[3-4].
基于此,本文考虑如下半线性的奇异摄动反应扩散方程:
(1)
其中ε∈(0,1]为小参数,a(x,u(x))∈C1([0,1]×R)且满足
(2)
当参数ε→0时, 问题(1)的解在区间[0,1]的两端存在边界层.由于边界层出现, 标准的数值方法在均匀网格下很难获得理想的数值结果.目前非均匀网格方法主要包括层适应网格方法和自适应移动网格方法.层适应网格方法由于构造简单,且容易分析其收敛性而受到许多学者的青睐[5-9].然而,层适应网格方法必须事先给出边界层的位置及宽度.相比之下,自适应移动网格方法更加适合于奇异摄动问题的数值求解[10-11].
Chadha和Kopteva[12]针对奇异摄动半线性反应扩散方程,构造出一种具有二阶精度的自适应移动网格算法.文献[13-14]研究了一维奇异摄动反应扩散方程的二次C1-样条配置法,分别给出了算法的后验误差估计和自适应网格生成算法.Linβ[15]讨论了奇异摄动一维反应扩散方程的有限元方法,证明了最大范数的后验误差估计,并设计了相应的自适应网格算法.Kopteva[16]在任意网格上构造了奇异摄动半线性反应扩散方程的有限差分方法,利用Green函数的相关性质推出了数值格式的后验误差估计.Kopteva等[17]考虑了一维线性奇异摄动反应扩散方程的自适应网格算法,并利用数值格式的局部截断误差估计设计了不同的网格控制函数,同时指出基于经典的弧长控制函数的自适应网格算法只能达到一阶收敛.值得一提的是,文献[16]给出的后验误差估计比较复杂,文献[17]没有给出后验误差估计的具体证明过程.因此,针对奇异摄动反应扩散方程的自适应移动网格算法,给出简单有效的后验误差估计显得尤为必要.
本文在文献[16-17]的基础上,进一步考虑问题(1)的自适应移动网格算法.首先,在任意非均匀网格上构造迎风有限差分格式,并利用多项式插值技术和精确解u(x)的稳定性结果,给出最大范数的后验误差估计和相应的网格生成算法.最后的数值实验表明本文提出的自适应网格算法具有二阶收敛.另外,本文构造的后验误差估计只需a(x,u(x))∈C1([0,1]×R),而文献[16-17]中提出的后验误差估计要求a(x,u(x))∈C2([0,1]×R).
注1 在本文中,C表示与ε和网格大小N无关的任意常数,在不同的位置可以表示不同的数值.在我们的估计中,对于任意函数
定义连续最大范数为
对于实值网格函数φ
{φ(xi)
,定义离散最大范数为‖φ‖∞=max0≤i≤N|φ(xi)|.
为了得出连续问题(1)的稳定性结果,我们首先给出方程(1)的线性化形式:
![]()
-ε2u″(x)+b(x)u(x)=-a(x,0),
(3)
其中b(x)=au(x,λu(x)),0<λ<1.显然,线性差分算子
满足极大值原理.此外,对于精确解u(x)有如下的稳定性结论.
引理1[16] 设u(x)为问题(1)的解,则有如下的稳定性估计:
(4)
推论1 对于任意两个函数v(x)与w(x),满足v(0)=w(0),v(1)=w(1),且
Lv(x)-Lw(x)=F(x),
其中F(x)为分段连续函数,则有
(5)
设
为任意的非均匀网格,其中N为正整数,对于1≤i≤n,取网格步长hi=xi-xi-1.在网格
上,使用如下有限差分格式离散问题(1):
(6)
其中Ui为u(xi)的数值解,
记U=(U0,U1,…,UN)T,从文献[16]的引理2.4可以看出,差分算子LN满足离散极大值原理,故满足如下的稳定性结果:
(7)
在本节中,为了得到离散格式(6)的后验误差估计,在任意区间Ji=[xi-1,xi],i=1,2,…,N上构造如下分段二次函数![]()
(8)
显然,由式(8)可得
(9)
定理1 设u(x)为问题(1)的精确解,{Ui
为通过离散格式(6)计算出的数值解,
为由式(8)所定义的分段二次函数,则
(10)
其中
证明 由离散格式(6),对任意x∈(xi-1,xi)有
(11)
进一步可得
再由推论1可以证得式(10).
推论2 设在问题(1)中,a(x,u(x))∈C1([0,1]×R)且满足式(2),对于任意x∈Ji,1≤i≤N,有
(12)
其中
证明 对任意x∈(xi-1,xi),对
在xi点Taylor展开有
由定理1得到
再结合式(8)有
则显然有式(12)得证.
在本节中,基于等分布原理[12],根据定理1中的后验误差估计来设计自适应网格算法.令
(13)
根据网格等分布原理,关键的问题是找到节点{(xi,Ui)
,使其满足
(14)
其中![]()
进而选择如下的控制函数:
(15)
为了获得满足式(14)的网格节点
构造如下的网格生成算法.
Step 1 网格初始化:确定区间数N,取初始网格
为均匀网格,其中
取定常数γ*>1.
Step 2 在当前网格
上使用差分格式(6)计算数值解
其中
分别在区间
上使用式(13)计算
进一步计算控制函数![]()
Step 3 取
计算![]()
Step 4 设Yi=iLN/N,i=0,1,…,N,由点
构造分段线性插值函数φ,用φ(Yi)估计新网格![]()
Step 5 网格判断:如果γ(k)≤γ*,则将
作为最终网格,计算数值解
作为最终解,退出;否则转step 2继续迭代.
为了验证本文所给出的后验误差估计及相应的自适应移动网格算法的有效性,考虑如下的奇异摄动反应扩散方程[18]:
-ε2u″(x)+(1+u)(1+(1+u)2)=0, x∈(0,1), u(0)=u(1)=0.
(16)
由于离散格式(6)是一个非线性系统,我们采用经典的Newton迭代格式求解,其中迭代初始解向量选择0,迭代终止条件设为‖UN,(k+1)-UN,(k)‖∞≤N-2,其中UN,(k)为第k次的迭代解向量,k=0,1,….
由于方程(16)的精确解是未知的,为计算最大误差和收敛阶,使用双重网格的方法来估计误差.设
为离散格式(6)在网格{xi
上计算得到的数值解,
为离散格式(6)在网格
上计算得到数值解.记U2N(x)为{xi
在节点
上的分段线性插值函数,则可以如下定义数值解的最大误差:
(17)
对于每一个N及所有的ε,令
进一步地,收敛阶可定义为
(18)
由定理1的后验误差估计,可选择如下控制函数:
(19)
类似地,由推论2,可构造如下控制函数:
(20)
在具体计算中,依次选择了ε=10-j, j=0,1,…,7,在N分别取2k,k=6,7,…,12时,利用控制函数(19)和(20)构造的自适应移动网格方法计算得到的误差和收敛阶分别见表1、表2,其中对于每一个ε,第一行表示对应于每个N的最大误差,第二行表示收敛阶.表1中,对于参数ε=10-j, j=0,1,…,6,上述优化算法中的终止常数γ*取3.0,但对于ε=10-7,γ*取8.0.
表1 基于控制函数(19)的自适应移动网格算法的计算结果
Table 1 Numerical results of the presented adaptive grid method based on monitor function(19)
εN=26N=27N=28N=29N=210N=211N=212Rε1007.25E-61.81E-64.53E-71.13E-72.83E-87.08E-91.77E-9-2.002.002.002.002.002.00-2.0010-12.69E-47.10E-52.01E-56.36E-62.06E-66.41E-71.99E-7-1.921.821.661.631.681.69-1.7310-25.23E-41.15E-42.61E-56.21E-61.50E-63.65E-79.01E-8-2.192.142.072.052.042.02-2.0810-34.89E-41.82E-41.21E-41.93E-42.16E-64.61E-71.09E-7-1.430.59-0.686.482.232.08-2.0210-47.68E-34.00E-33.07E-56.48E-67.08E-68.31E-66.65E-6-0.947.032.24-0.13-0.230.32-1.7010-58.71E-31.20E-49.19E-41.26E-41.79E-63.92E-79.18E-8-6.18-2.942.866.142.192.09-2.7610-61.22E-22.28E-31.51E-38.13E-41.35E-41.50E-51.06E-7-2.430.590.902.593.177.14-2.8010-73.15E-28.39E-31.65E-34.55E-41.87E-61.33E-41.04E-6-1.912.351.867.92-6.147.00-2.48EN3.15E-28.39E-31.65E-38.13E-41.35E-41.33E-46.65E-6-1.912.351.022.590.034.32-2.03
表2中,对于参数ε=10-j, j=0,1,3,4,7,γ*取3.0,对于参数ε=10-j, j=2,5,6,γ*取4.5.我们还将每一个ε对应的误差阶的平均值记为Rε,见表1、表2的最后一列.
表2 基于控制函数(20)的自适应移动网格算法的计算结果
Table 2 Numerical results of the presented adaptive grid method based on monitor function(20)
εN=26N=27N=28N=29N=210N=211N=212Rε1007.25E-61.81E-64.53E-71.13E-72.83E-87.08E-91.77E-9-2.002.002.002.002.002.00-2.0010-18.14E-42.08E-45.22E-51.31E-53.27E-68.17E-72.04E-7-1.971.992.002.002.002.00-1.9910-21.44E-32.78E-45.27E-51.19E-52.73E-66.13E-71.41E-7-2.382.402.142.132.152.12-2.2210-31.34E-31.74E-47.13E-52.26E-54.74E-61.09E-62.50E-7-2.941.291.662.252.132.12-2.0610-42.28E-31.74E-45.46E-51.52E-52.46E-64.88E-78.58E-8-3.721.671.842.632.342.51-2.4510-52.29E-27.46E-42.50E-41.67E-51.07E-52.85E-61.68E-7-4.941.583.910.641.914.09-2.8410-61.34E-22.14E-23.66E-44.05E-56.13E-63.99E-62.97E-7--0.685.873.172.720.623.75-2.5810-72.23E-27.82E-31.59E-33.07E-45.25E-51.33E-62.95E-7-1.512.302.372.555.302.18-2.70EN2.29E-22.14E-21.59E-33.07E-45.25E-53.99E-62.97E-7-0.093.752.372.553.723.75-2.70
表3 均匀网格下的计算结果
Table 3 Numerical results calculated on a uniform mesh
εN=26N=27N=28N=29N=210N=211N=2121007.25E-61.81E-64.53E-71.13E-72.83E-87.08E-91.77E-92.002.002.002.002.002.00-10-18.14E-42.08E-45.22E-51.31E-53.27E-68.17E-72.04E-71.971.992.002.002.002.00-10-22.76E-21.43E-24.47E-31.26E-33.23E-48.15E-52.04E-50.951.681.821.971.992.00-10-33.81E-31.25E-22.71E-23.04E-21.88E-26.77E-31.93E-3-1.71-1.12-0.160.691.471.81-10-44.09E-51.63E-46.48E-42.50E-38.77E-32.26E-23.15E-2-2.00-1.99-1.95-1.81-1.37-0.48-10-54.10E-71.64E-66.55E-62.62E-51.05E-44.16E-41.63E-3-2.00-2.00-2.00-2.00-1.99-1.97-10-64.10E-91.64E-86.55E-82.62E-71.05E-64.19E-61.68E-5-2.00-2.00-2.00-2.00-2.00-2.0010-74.10E-111.64E-106.55E-102.62E-91.05E-84.19E-81.68E-7-2.00-2.00-2.00-2.00-2.00-2.00-EN2.76E-21.43E-22.71E-23.04E-21.88E-22.26E-23.15E-20.95-0.92-0.160.69-0.27-0.48-
同时,为了进行比较,使用均匀网格上的有限差分格式(6)来计算该问题,相应的数值结果见表3.
从表1~3可以看出,对于ε=10-j, j=0,1,2,在均匀网格上,使用有限差分格式(6)也可以获得最佳的收敛阶.但当ε=10-j, j=3,4,5,6,7时,在自适应网格上获得的数值结果明显要比在均匀网格上获得的数值结果更好.此外,从表1、表2中每个ε对应的平均收敛阶可以看出,所提出的自适应网格方法是二阶收敛的,且独立于摄动参数ε.
为了演示网格生成算法中网格的移动过程,当N=64,ε=10-5,10-2时,图1和图2分别给出了基于控制函数(19)和(20)的自适应网格算法的网格迭代过程.同时,图3给出了当N=64,ε=10-2时,自适应网格算法计算的数值解.从图2、图3可以看出,在x=0和x=1附近的点比较密集,也表明问题(16)在这两个点处表现为两个边界层.最后,图4给出了N分别取2k,k=8,9,10,ε依次取10-j, j=7,6,…,1,0时,误差
随参数ε的变化曲线图.显然,
随着ε的增大而减少.
图1 网格的演化过程(ε=10-5, N=64)
Fig. 1 Evolution of the mesh(ε=10-5, N=64)
图2 网格的演化过程(ε=10-2, N=64)
Fig. 2 Evolution of the mesh(ε=10-2,N=64)
图3 数值解(ε=10-2, N=64)
Fig. 3 Numerical solutions(ε=10-2, N=64)
图4 误差
随参数ε的变化曲线
Fig. 4 The variation curve of
with ε
本文主要针对奇摄动反应扩散方程,设计了一种自适应网格方法.基于二次多项式插值,得出后验误差估计,用于构造自适应网格算法.数值实验表明,使用本文构造的方法计算的数值解具有二阶精度.值得指出的是,所提出的自适应网格方法还可以推广到其他的奇异摄动反应扩散问题中.
致谢 本文作者衷心感谢池州学院自然科学重点项目(CZ2018ZR06)对本文的资助.
[1] 史娟荣, 莫嘉琪. 一类奇异摄动燃烧模型的渐近解[J]. 应用数学和力学, 2016, 37(7): 691-698.(SHI Juanrong, MO Jiaqi. Asymptotic solutions to a class of singular perturbation burning models[J]. Applied Mathematics and Mechanics, 2016, 37(7): 691-698.(in Chinese))
[2]
R, TEOFANOV L. A uniform numerical method for semilinear reaction-diffusion problems with a boundary turning point[J]. Numerical Algorithms, 2010, 54: 431-444.
[3] GENG F, QIAN S, LI S. Numerical solutions of singularly perturbed convection-diffusion problems[J]. International Journal of Numerical Methods for Heat & Fluid Flow, 2014, 24(6): 1268-1274.
[4] RAMOS H, VIGO-AGUIAR J, NATESAN S, et al. Numerical solution of nonlinear singularly perturbed problems on nonuniform meshes by using a non-standard algorithm[J]. Journal of Mathematical Chemistry, 2010, 48(1): 38-54.
[5] ADLER J, MACLACHLAN S, MADDEN N. A first-order system Petrov-Galerkin discretization for a reaction-diffusion problem on a fitted mesh[J]. IMA Journal of Numerical Analysis, 2015, 36(3): 1-29.
[6] BRDAR M, ZARIN H. A singularly perturbed problem with two parameters on a Bakhvalov-type mesh[J]. Journal of Computational and Applied Mathematics, 2016, 292: 307-319.
[7] 蔡新, 蔡丹琳, 吴瑞潜, 等. 奇异摄动反应扩散问题的高阶不等距计算方法[J]. 应用数学和力学, 2009, 30(2): 171-178.(CAI Xin, CAI Danlin, WU Ruiqian, et al. High accurate non-equidistant method for singular perturbation reaction-diffusion problem[J]. Applied Mathematics and Mechanics, 2009, 30(2): 171-178.(in Chinese))
[8] 邵文婷, 郑烁宇. 奇异摄动内层问题的宽度估计及其数值求解[J]. 计算机工程与应用, 2020, 56(4): 44-49.(SHAO Wenting, ZHENG Shuoyu. Width estimation of singular perturbed interior layer problem and its numerical solution[J]. Computer Engineering and Applications, 2020, 56(4): 44-49.(in Chinese))
[9] LIU L B, CHEN Y P. A posteriori error estimation in maximum norm for a strongly coupled system of two singularly perturbed convection diffusion problems[J]. Journal of Computational and Applied Mathematics, 2017, 313: 152-167.
[10] LIU L B, CHEN Y P. A robust adaptive grid method for a system of two singularly perturbed convection-diffusion equations with weak coupling[J]. Journal of Scientific Computing, 2014, 61(1): 1-16.
[11] 毛志, 刘利斌. 一类奇异摄动问题的自适应移动网格算法[J]. 湘潭大学学报(自然科学版), 2019, 41(1): 64-71.(MAO Zhi, LIU Libin. An adaptive moving grid algorithm for a system of singularly perturbed problems[J]. Journal of Xiangtan University (Natural Science Edition), 2019, 41(1): 64-71.(in Chinese))
[12] CHADHA N M, KOPTEVA N. A robust grid equidistribution method for a one-dimensional singularly perturbed semilinear reaction-diffusion problem[J]. IMA Journal of Numerical Analysis, 2011, 31(1): 188-211.
[13] LINB T, RADOJEV G, ZARIN H. Approximation of singularly perturbed reaction-diffusion problems by quadratic C1-splines[J]. Numerical Algorithms, 2012, 61: 35-55.
[14] LINB T, RADOJEV G. Robust a posteriori error bounds for spline collocation applied to singularly perturbed reaction-diffusion problems[J]. Electronic Transactions on Numerical Analysis, 2016, 45: 342-353.
[15] LINB T. A posteriori error estimation for arbitrary order fem applied to singularly perturbed one-dimensional reaction-diffusion problems[J]. Applications of Mathematics, 2014, 59(3): 241-256.
[16] KOPTEVA N. Maximum norm a posteriori error estimates for a 1D singularly perturbed semilinear reaction-diffusion problem[J]. IMA Journal of Numerical Analysis, 2006, 27(3): 1-14.
[17] KOPTEVA N, MADDENB N, STYNES M. Grid equidistribution for reaction-diffusion problems in one dimension[J]. Numerical Algorithms, 2005, 40(3): 305-322.
[18] RAO S C S, KUMAR S, KUMAR M. A parameter-uniform B-spline collocation method for singularly perturbed semilinear reaction-diffusion problems[J]. Journal of Optimization Theory and Applications, 2010, 146: 795-809.