APP下载

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

2022-08-18程晓晗封建湖

郑州大学学报(理学版) 2022年5期
关键词:算例线性数值

程晓晗, 封建湖

(长安大学 理学院 陕西 西安 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 结论

猜你喜欢

算例线性数值
体积占比不同的组合式石蜡相变传热数值模拟
数值大小比较“招招鲜”
舰船测风传感器安装位置数值仿真
铝合金加筋板焊接温度场和残余应力数值模拟
关于非齐次线性微分方程的一个证明
非齐次线性微分方程的常数变易法
提高小学低年级数学计算能力的方法
线性耳饰
论怎样提高低年级学生的计算能力
试论在小学数学教学中如何提高学生的计算能力