求解Hamilton-Jacobi方程的四阶WENO格式

程晓晗， 封建湖

(长安大学 理学院 陕西 西安 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的数值结果

3 结论