求解Hamilton-Jacobi方程的四阶WENO格式
2022-08-18程晓晗封建湖
程晓晗, 封建湖
(长安大学 理学院 陕西 西安 710064)
0 引言
1 数值方法
考虑一维H-J方程初值问题
(1)
为简单起见,将待求解区域均匀划分为N个网格单元Ii=[xi-1/2,xi+1/2],其中h=xi+1/2-xi-1/2为网格步长。求解式(1)的半离散形式为
(2)
(3)
图1 重构模板
步骤1 给定子模板S0={xi-2,xi-1,xi},S1={xi-1,xi,xi+1},S2={xi,xi+1,xi+2}(如图1所示),构造线性多项式p0(x),p1(x),p2(x)满足
i-1+k,k=0,1,2,
(4)
则可得φx(xi)的三种二阶近似为
(5)
步骤2 若在模板S={S0,S1,S2}上构造三次多项式p(x)满足
(6)
则可得φx(xi)的四阶近似为
(7)
(8)
其中线性权为d0=1/6,d1=2/3,d2=1/6。
步骤4 类似于Borges等的思想[8],选取非线性权为
(9)
ε的引入是为了防止分母为零,本文中选取ε=10-6。β0、β1为模板S0、S1的光滑因子,即
(10)
β2为模板S的光滑因子,即
1 922Δφi-494Δφi+1)+Δφi-1(3 443Δφi-1-
5 966Δφi+1 602Δφi+1)+Δφi(2 843Δφi-
(11)
选取参数τ为
(12)
将β0、β1、β2在点xi处进行泰勒展开有
(13)
将式(13)代入式(12)可得τ=O(h6)。当ε≪βk时有
(14)
(15)
步骤6 对半离散格式(2)在时间方向上的离散采用三阶Runge-Kutta方法[9]
(16)
就可得到下一时间层的解。
以上方法可按照逐维计算的方法直接推广至多维情形, 在此不再赘述。
2 数值试验
下面给出几个一维和二维经典算例来检验本文的四阶WENO格式(WENO4)的性能。
算例1考虑线性对流方程
该方程的精确解是φ(x,t)=-cos(π(x-t)),可用来测试格式的数值精度。为了使时间方向也能达到四阶精度,我们取时间步长Δt=Δx4/3。采用不同网格点数计算到t=2,在区间[-1,1]上的L1和L∞误差如表1所示。可以看出,WENO4格式达到了理论上的四阶精度。
表1 格式的数值精度(一维情形)
注:“-”表示无数值。
算例2考虑非线性方程
其中H分别取一个凸算子H(u)=(u+1)2/2和非凸算子H(u)=-cos(u+α)。这两种情形都会在t=1/π2时产生奇性。采用40个网格点计算,图2给出了t=2/π2时刻的数值结果,其中实线是参考解,由本文格式采用6 400个网格点计算得到。从图2中可以看出,WENO4格式成功地捕捉到了“尖角”处。
算例3考虑二维线性对流方程
该方程的精确解是φ(x,y,t)=-cos(π(x+y-2t)),可以用来测试二维情况下格式的数值精度。采用不同网格点数计算到t=1,选取区域[-1,1]×[-1,1]上的L1和L∞误差结果如表2所示。可以看出,二维情形下WENO4格式也可以达到理论上的四阶精度。
表2 格式的数值精度(二维情形)
注:“-”表示无数值。
算例4考虑二维Burgers方程
该问题在t=0.5/π2时刻解是光滑的,在t=1.5/π2时刻解会产生间断的一阶导数。采用40×40网格点计算,其结果如图3所示。可以看出,数值解与经典文献[4, 6, 10]的数值解相吻合。
图3 算例4的数值结果
算例5考虑几何光学中的Eikonal方程[6]
采用40×40网格点计算到t=0.6,其等势线与曲面图如图4所示。可以看出,WENO4格式成功捕捉到了本问题发展成的“尖角”部分。
图4 算例5的数值结果
算例6考虑优化控制问题[6]
与前面算例的通量仅依赖于▽φ不同,本问题的通量还依赖于(x,y)。经过有限时刻,解的导数会产生间断,最感兴趣的量是ω=sign(φy)。采用40×40网格点计算到t=1,图5给出了解φ和ω=sign(φy)的数值计算结果,可以看出WENO4格式能较好分辨出解的间断导数。
图5 算例6的数值结果