1、2023 年 6 月第 39 卷 第 2 期纯粹数学与应用数学Pure and Applied MathematicsJun.2023Vol.39 No.2二维奇异摄动对流扩散方程的自适应移动网格算法刘利斌1,2,徐磊1,包小兵3(1.南宁师范大学数学与统计学院,广西 南宁530100;2.南宁师范大学广西人机交互与智能决策重点实验室,广西 南宁530100;3.池州学院大数据与人工智能学院,安徽 池州247000)摘要：针对二维奇异摄动对流扩散方程,在任意网格下给出了经典的迎风有限差分格式.利用二元多项式插值技术,推导出一阶最大范数的后验误差估计,并以此设计了一个自适应网格生成算法.数值实验
2、表明本文构造的自适应移动网格算法是有效的.关键词：奇异摄动;后验误差估计;迎风有限差分格式;自适应移动网格算法中图分类号：O241.81文献标识码：A文章编号：1008-5513(2023)02-0199-15DOI：10.3969/j.issn.1008-5513.2023.02.0031 引引引言言言本文考虑如下二维奇异摄动对流扩散方程Lu:=u b1ux b2uy+cu=f,(x,y)=(0,1)2,(1)u=0,(x,y)=,(2)其中 0 1 0,b2(x,y)2 0,c(x,y)0,(x,y).(3)收稿日期：2022-03-23.接收日期：2022-05-11.基金项目：广西省自
3、然科学基金(2020GXNSFAA159010);安徽省优秀青年人才项目(gxyq2021225).作者简介：刘利斌(1982-),博士,教授,研究方向：奇异摄动问题的数值方法.200纯粹数学与应用数学第 39 卷在条件(3)及函数 f 满足一定的相容性条件下(见文献 1),问题(1)-(2)存在唯一解 u(x,y).当 0,问题(1)-(2)的解 u 在边界 x=0 和 y=0 处存在边界层2.众所周知,奇异摄动对流扩散方程在流体力学,地下油藏,化学反应等领域有着广泛的应用,其中含高雷诺数的 Navier-Stokes 方程即为奇异摄动对流扩散方程3.这类问题的一个显著特点是高阶导数项含有一
4、个小参数,该参数 的出现会使解在某些区域变化非常剧烈,即存在边界层现象.由于边界层的出现,常用的数值方法(例如有限差分和有限元方法)在均匀网格下得到数值解的误差往往与 1/有关,进而当 0 时,其数值方法无法获得一致收敛(即数值解的误差与 无关)的数值解.正因为如此,许多学者一直在不断尝试构造一致收敛的数值方法.迄今为止,常见的一致收敛的数值方法有层适应网格方法(Shishkin网格和 Bakhvalov 网格)和自适应移动网格算法,见专著 2-4.层适应网格是一种特殊网格,其构造方法比较简单,因而受到了众多学者的广泛关注,其发表的文献也有许多.例如,文献 5-6 的作者针对二维奇异摄动对流扩
5、散方程(1)-(2),在 Shishkin 网格下分别构造了迎风有限差分格式和混合差分格式,并证明了数值格式的一致收敛性.Nhan和 Vulanovi7-8在 Bakhvalov 网格下分别讨论了一维和二维奇异摄动对流扩散方程的参数一致的数值方法.与层适应网格不同的是,自适应移动网格方法是一种后处理方法,该方法不要精确解的任何先验信息.因而在微分方程数值方法研究中越来越受到学者们的重视9-11.在文献 9-11 的基础上,周琴和杨银12-13分别讨论了含源项非定常对流扩散方程和二阶双曲型方程的移动网格方法.与此同时,基于迎风有限差分格式,文献 14-19 的作者讨论了若干一维奇异摄动问题的后验
6、误差估计,并利用控制函数和网格等分布原理设计了相应的自适应网格生成算法,其中 Kopteva 和 Stynes15中指出尽管一维奇异摄动对流扩散方程的自适应移动网格算法不能直接推广到二维问题,但是其它相关的移动网格算法已经用于某些二维问题的数值求解20-22,并取得了较好的数值结果.值得一提的是虽然文献 20-22 设计了二维问题的移动网格方法,但是没有给出误差估计.由此可见,研究二维奇异摄动问题的后验误差估计,并设计相应的自适应移动网格算法具有一定的意义.基于此,本文针对二维奇异摄动对流扩散方程(1)-(2),在任意网格下给出了经典的迎风有限差分格式.并且基于二元多项式插值函数,构造了数值解
7、的后验误差估计,利用控制函数和网格等分布原理,在 x 轴方向和 y 轴方向分别构造了类似的自适应移动网格生成算法.最后的数值实验验证了本文提出的自适应移动网格算法的有效性.数值结果也表明该自适应移动网格算法能达到最优的一阶一致收敛.注注注1.1在整篇论文中,C 表示一个与 和 N 无关的常数,且在不同的位置取不同的值.第 2 期刘利斌 等:二维奇异摄动对流扩散方程的自适应移动网格算法2012 预预预备备备知知知识识识为了推导出数值解的后验误差估计,首先考虑如下非齐次边界的奇异摄动对流扩散方程Lv:=v b1vx b2vy+cv=f,(x,y)=(0,1)2,(4)v=g,(x,y)=,(5)其
8、中系数 bi,i=1,2,c 和源项 f 与问题(1)-(2)中的系数一致.显然,问题(4)-(5)中的微分算子 L 满足极大值原理2,则问题(4)-(5)存在唯一解 v,且满足如下稳定性估计:引引引理理理2.123令 v(x,y)是问题(4)-(5)的精确解,则在假设条件(3)下,有v C(f+g,),(6)其中 f=max(x,y)|f|,g,=max(x,y)|g|.推推推论论论2.1令函数,满足如下条件L L=F,(x,y)=(0,1)2,=G,(x,y)=,其中 F 为分片光滑的连续函数,则 C(L L+G,).令N=(xi,yj):i,j=0,1,N 为任意的非均匀网格,其中节点
9、xi和 yj满足 0=x0 x1 xN=1,0=y0 y1 yN=1.对于任意的 1 i,j N,记 hx,i=xixi1,x,i=(hx,i+1+hx,i)/2,且 hy,j=yjyj1,y,j=(hy,j+1+hy,j)/2.同时,对于 1 i,j N 1 及任意的网格函数 i,j,定义如下差分算子D2xi,j=1x,i(D+xi,j Dxi,j),D2yi,j=1y,j(D+yi,j Dyi,j),Dxi,j=i,j i1,jhx,i,D+xi,j=i+1,j i,jhx,i+1,Dyi,j=i,j i,j1hy,j,D+yi,j=i,j+1 i,jhy,j+1.在上述网格下,可构造问题
10、(1)-(2)的迎风有限差分格式如下:LNuNi,j:=(D2x+D2y)b1,ijD+x b2,ijD+y+ci,juNi,j=fi,j,1 i,j N 1,(7)uNi,j=0,(xi,yj)N,(8)202纯粹数学与应用数学第 39 卷其中 uNi,j是 u(x,y)在网格点(xi,yj)处的数值解,N=N.接下来,基于离散的极大值原理2,差分格式(7)-(8)的解满足如下稳定性结果.引引引理理理2.2令 uNi,j为离散格式(7)-(8)在任意网格下的解,则max(i,j)?uNi,j?C f,N,其中 f,N=max(i,j)NN|fi,j|.3 后后后验验验误误误差差差估估估计计计
11、及及及自自自适适适应应应网网网格格格算算算法法法在这一部分,首先给出数值解的后验误差估计.然后,根据后验误差估计和网格等分布原理设计出相应的自适应网格生成算法.首先,结合边界条件(8),有D2yuN0,j=D2yuNN,j=0,j=1,N 1,D+yuN0,j=D+yuNN,j=0,j=0,1,N 1,D2xuNi,0=D2xuNi,N=0,i=1,N 1,D+xuNi,0=D+xuNi,N=0,i=0,1,N 1.同时,假设差商 D+xuNN,j,j=0,1,N 和 D+yuNi,N,i=0,1,N 存在,则可将离散格式(7)延拓到网格节点 i=0 和 i=N,得D2xuN0,j:=1f(0
12、,yj)b1(0,yj)D+xuN0,j,j=0,1,N,(9)D2xuNN,j:=1f(1,yj)b1(1,yj)D+xuNN,j,j=0,1,N.(10)类似地,D2yuNi,0:=1f(xi,0)b2(xi,0)D+yuNi,0,i=0,1,N,(11)D2yuNi,N:=1f(xi,1)b2(xi,1)D+yuNi,N,i=0,1,N.(12)显然,由上面的延拓关系(9)-(12)可看出,当(xi,yj)分别取(0,0),(0,1),(1,0),(1,1)时,D2xuNi,j=D2yuNi,j=0.进一步可得出离散格式(7)对于 i,j=0,1,N 都成立.第 2 期刘利斌 等:二维奇
13、异摄动对流扩散方程的自适应移动网格算法203接下来,对于(x,y)Ii,j=xi,xi+1 yj,yj+1,i,j=0,1,N 1,定义如下分片二次插值函数:uN(x,y)=12D2xuNi,j(x xi)(x xi+1)+12D2yuNi,j(y yj)(y yj+1)+1hx,i+1hy,j+1uNi+1,j+1(x xi)(y yj)+uNi,j+1(xi+1 x)(y yj)+1hx,i+1hy,j+1uNi+1,j(x xi)(yj+1 y)+uNi,j(xi+1 x)(yj+1 y).(13)显然,对于(x,y)Ii,j,有uN(x,y)x=D2xuNi,j(x xi+1/2)+1
14、hy,j+1D+xuNi,j+1(y yj)+D+xuNi,j(yj+1 y),(14)uN(x,y)y=D2yuNi,j(y yj+1/2)+1hx,i+1D+yuNi+1,j(x xi)+D+yuNi,j(xi+1 x),(15)2uN(x,y)x2=D2xuNi,j,2uN(x,y)y2=D2yuNi,j,uN(xi,yj)=uNi,j,(16)uN(xi,yj+1)=uNi,j+1,uN(xi+1,yj)=uNi+1,j,uN(xi+1,yj+1)=uNi+1,j+1,(17)其中 xi+1/2=(xi+1+xi)/2,yj+1/2=(yj+1+yj)/2.最后,基于上述插值函数(13
15、),可得到如下后验误差估计:定定定理理理3.1令 u(x,y)为方程(1)-(2)的精确解,uNi,j为离散格式(7)-(8)的解,且 uN(x,y)为(13)式中定义的分片二次插值函数.则?u(x,y)uN(x,y)?Cmax0i,jN1hi+1,j+1(1+i,j+i,j+?D+yD+xuNi,j?),(18)其中hi+1,j+1=maxhx,i+1,hy,j+1,i,j=?D2xuNi,j?+max?D+xuNi,j+1?,?D+xuNi,j?,i,j=?D2yuNi,j?+max?D+yuNi+1,j?,?D+yuNi,j?.证证证明明明对于(x,y)Ii,j,i,j=0,1,N 1,
16、由方程(1),离散格式(7)及延拓关系(9)-(12),有Lu(x,y)LuN(x,y)=f(x,y)+(D2xuNi,j+D2yuNi,j)+b1uN(x,y)x+b2uN(x,y)y cuN(x,y)=R1+R2+R3+R4,(19)204纯粹数学与应用数学第 39 卷其中R1=f(x,y)f(xi,yj),(20)R2=b1uN(x,y)x b1,ijD+xuNi,j,(21)R3=b2uN(x,y)y b2,ijD+yuNi,j,(22)R4=ci,juNi,j cuN(x,y).(23)因为 f(x,y)为充分光滑的函数,则(x,y)Iij,i,j=0,1,N 1,由(20)式可得|R1|=|f(x,y)f(xi,yj)|C(hx,i+1+hy,j+1)Chi+1,j+1.(24)令(x,y)=b1(x,y)uN(x,y)x,由(14)式可得(xi,yj)=b1,ijD2xuNi,jhx,i+1/2+b1,ijD+xuNi,j,且?(x,y)x?C(?D2xuNi,j?hx,i+1+?D+xuNi,j+1?+?D+xuNi,j?+?D2xuNi,j?)C(?D2xuNi,j?