程晓晗， 封建湖
(长安大学 理学院 陕西 西安 710064)
考虑一维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)
就可得到下一时间层的解。
以上方法可按照逐维计算的方法直接推广至多维情形, 在此不再赘述。
下面给出几个一维和二维经典算例来检验本文的四阶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的数值结果