求解非线性退化抛物方程的HWENO-LW格式

HWENO-LW Scheme for Solving Nonlinear Degenerate Parabolic Equations

摘要: 本文构造了基于Lax-Wendroff时间离散的有限差分HWENO (Hermite加权本质非振荡)格式，用于求解非线性退化抛物方程。与传统的Runge-Kutta时间离散方法相比，Lax-Wendroff方法在提高计算效率的同时，能够在解的光滑区域实现时空一致的高阶精度。通过数值算例，验证了该方法的有效性。

Abstract: This paper constructs a finite difference HWENO (Hermite Weighted Essentially Non-Oscillatory) scheme based on Lax-Wendroff time discretization for solving nonlinear degenerate parabolic equations. Compared to traditional Runge-Kutta time discretization methods, the Lax-Wendroff method improves computational efficiency while achieving high-order accuracy in both space and time in smooth regions of the solution. The effectiveness of the method is validated through numerical examples.

1. 引言

对流扩散方程是粘性流体力学及其他实际问题中常见的一类基本运动方程，在环境科学，流体力学，气体动力学等领域有着广泛应用。非线性抛物方程是其中一类重要的方程，常用于描述非线性质量传递、多孔介质的毛吸附、热辐射扩散、不可压缩流体运动力学等问题。非线性抛物方程解具有双曲守恒律方程解的类似性质，即真解会随着时间的推移，出现间断和陡峭的过渡层，使得精确解很难获得，同样给数值求解带来了很多困难，比如陡峭界面和紧支集界面以有限的速度传播。

近年来，流体力学领域的高精度数值计算方法取得了快速发展，并广泛应用于工程实践。常见的数值方法包括TVD格式[1]，ENO格式[2]以及WENO格式[3]-[8]，2005年Qiu等人提出了HWENO格式[9] [10]，并应用于Euler方程和Hamilton-Jacobi方程的求解。该方法通过Hermite插值多项式提高数值模板的紧凑性，同时保持数值格式的稳定性。此外，针对Euler方程，Hamilton-Jacobi方程和浅水波方程等问题，需选择适当的时间离散方法，如具有TVD性质的Runge-Kutta法、线性多步法[11]或类似于Lax-Wendroff格式[12]的Lax-Wendroff时间离散[13]-[15]方法，以实现时空全离散的数值求解。

本文结合Lax-Wendroff时间离散方法与HWENO格式来求解非线性退化抛物方程。Lax-Wendroff方法能实现时空一致的高阶精度。与同阶WENO格式相比，在相同模板条件下，除原函数点信息外，HWENO格式还具有导函数点的信息，从而使得在相同精度的格式条件下，HWENO格式所需模板点数更少。例如，六阶WENO格式需要六个模板点，而六阶HWENO仅需四个模板点。本文设计的六阶HWENO格式基于四点模板，通过三个三点四次函数的非线性组合来近似节点处的二阶导数。为了确保格式的稳定性，对HWENO格式构造过程中出现的负线性权重进行了特殊处理。数值结果验证了该格式的六阶精度及无振荡特性。

2. 非线性退化抛物方程的空间离散

2.1. 大模板上的空间重构

考虑一维情况下的非线性退化抛物方程，

$\{\begin{array}{l}{u}_t=l{\left(u\right)}_{xx},\\ u\left(x,0\right)=u_0\left(x\right),\end{array}$ (1)

为了简化计算，采用均匀网格$\left\{x_i\right\},1\le i\le N$ 进行离散化，$I_i=\left[x_i,x_{i+1}\right]$ 表示计算单元，且$\Delta x=x_{i+1}-x_i$ 是常数。为了构造HWENO格式，将方程(1)两边对x求导，从而得到下列方程组：

$\{\begin{array}{ll}{u}_t=l{\left(u\right)}_{xx},\hfill & u\left(x,0\right)=u_0\left(x\right),\hfill \\ v_t=z{\left(u,v\right)}_{xx},\hfill & v\left(x,0\right)=v_0\left(x\right).\hfill \end{array}$ (2)

方程组(2)的守恒有限差分半离散形式为：

$\{\begin{array}{l}\frac{\text{d}{u}_i\left(t\right)}{\text{d}t}=\frac{1}{\Delta x^2}\left({\widehat{l}}_{i+\frac{1}{2}}-{\widehat{l}}_{i-\frac{1}{2}}\right),\\ \frac{\text{d}{v}_i\left(t\right)}{\text{d}t}=\frac{1}{\Delta x^2}\left({\widehat{z}}_{i+\frac{1}{2}}-{\widehat{z}}_{i-\frac{1}{2}}\right).\end{array}$ (3)

定义函数$\phi \left(x\right)$ 满足

$\{\begin{array}{l}l{\left(u\left(x\right)\right)}_{xx}=\frac{1}{\Delta x^2}{\displaystyle {\int }_{x-\frac{\Delta x}{2}}^{x+\frac{\Delta x}{2}}{\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}\phi \left(\xi \right)\text{d}\xi \text{d}\eta ,}}\hfill \\ z{\left(u\left(x\right),v\left(x\right)\right)}_{xx}=\frac{1}{\Delta x^2}{\displaystyle {\int }_{x-\frac{\Delta x}{2}}^{x+\frac{\Delta x}{2}}{\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{\phi }^{\prime }\left(\xi \right)\text{d}\xi \text{d}\eta ,}}\hfill \end{array}$

在模板$S=\left\{x_{i-1},x_i,x_{i+1},x_{i+2}\right\}$ 上，利用函数点及导函数点得到埃尔米特插值多项式，得到数值通量：

$\begin{array}{l}{\widehat{l}}_{i+\frac{1}{2}}=p\left(x_{i+\frac{1}{2}}\right)=\frac{1}{54}\left(-135l\left(u_i\right)+135l\left(u_{i+1}\right)-7l\left(u_{i-1}\right)+7l\left(u_{i+2}\right)\right)\\ +\frac{\Delta x}{36}\left(-z\left(u_{i-1},v_{i-1}\right)-33z\left(u_i,v_i\right)-33z\left(u_{i+1},v_{i+1}\right)-z\left(u_{i+2},v_{i+2}\right)\right),\\ {\widehat{z}}_{i+\frac{1}{2}}=p^{\prime }\left(x_{i+\frac{1}{2}}\right)=\frac{1}{\Delta x}\left(l\left(u_i\right)+l\left(u_{i+1}\right)-l\left(u_{i-1}\right)+l\left(u_{i+2}\right)\right)\\ +\frac{1}{4}\left(-z\left(u_{i-1},v_{i-1}\right)-9z\left(u_i,v_i\right)+9z\left(u_{i+1},v_{i+1}\right)+z\left(u_{i+2},v_{i+2}\right)\right).\end{array}$

2.2. HWENO格式对原函数方程的重构

在小模板点$s_0=\left\{x_{i-1},x_i,x_{i+1}\right\},s_1=\left\{x_i,x_{i+1},x_{i+2}\right\},s_2=\left\{x_{i-1},x_i,x_{i+1},x_{i+2}\right\}$ 上利用函数点值及导函数点值构造埃尔米特多项式，满足：

$\begin{array}{l}\frac{1}{\Delta x^2}{\displaystyle {\int }_{I_{i+j}}\left({\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{h}_0\left(\xi \right)\text{d}\xi }\right)\text{d}\eta }=l\left(u_j\right),\frac{1}{\Delta x^2}{\displaystyle {\int }_{I_{i+j}}\left({\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{h^{\prime }}_0\left(\xi \right)\text{d}\xi }\right)\text{d}\eta }=z\left(u_j,v_j\right),j=-1,0,1;\\ \frac{1}{\Delta x^2}{\displaystyle {\int }_{I_{i+j}}\left({\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{h}_1\left(\xi \right)\text{d}\xi }\right)\text{d}\eta }=l\left(u_j\right),\frac{1}{\Delta x^2}{\displaystyle {\int }_{I_{i+j}}\left({\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{h^{\prime }}_1\left(\xi \right)\text{d}\xi }\right)\text{d}\eta }=z\left(u_j,v_j\right),j=0,1,2;\\ \frac{1}{\Delta x^2}{\displaystyle {\int }_{I_{i+j}}\left({\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{h}_2\left(\xi \right)\text{d}\xi }\right)\text{d}\eta }=l\left(u_j\right),j=-1,0,1,2;\\ \frac{1}{\Delta x^2}{\displaystyle {\int }_{I_{i+j}}\left({\displaystyle {\int }_{\eta -\frac{\Delta x}{2}}^{\eta +\frac{\Delta x}{2}}{h^{\prime }}_2\left(\xi \right)\text{d}\xi }\right)\text{d}\eta }=z\left(u_j,v_j\right),j=-1,2.\end{array}$

由$p^{\left(m\right)}\left(x\right)=h_m\left(x+\frac{\Delta x}{2}\right)-h_m\left(x-\frac{\Delta x}{2}\right)$ 得到四阶多项式

$p^{\left(0\right)}\left(x_{i+\frac{1}{2}}\right),p^{\left(1\right)}\left(x_{i+\frac{1}{2}}\right),p^{\left(2\right)}\left(x_{i+\frac{1}{2}}\right),$

利用$p\left(x_{i+\frac{1}{2}}\right)={\displaystyle \sum _{m=0}^2{d}_m{p}^{\left(m\right)}\left(x_{i+\frac{1}{2}}\right)}$ ，得到线性权$d_0=\frac{11}{7},d_1=\frac{11}{7},d_2=-\frac{15}{7}$ ，从而求得光滑指示子

${\beta }_m={\displaystyle \sum _{l=1}^2\Delta x^{2l-1}{\displaystyle {\int }_{I_i}{\left(\frac{{\text{d}}^l}{\text{d}{x}^l}{p}^{\left(m\right)}\left(x\right)\right)}^2\text{d}x}}.$ (4)

为确保格式的稳定性，采用特殊方法[16]将负线性权重转换为正线性权重

${\tilde{\sigma }}_m^{+}=\frac{1}{2}\left(d_m+\alpha \left|d_m\right|\right),{\tilde{\sigma }}_m^{-}={\tilde{\sigma }}_m^{+}-d_m,m=0,1,2,$

得到最终线性权${\sigma }_m^{\pm }={\tilde{\sigma }}_m^{\pm }/{\displaystyle \sum _{m=0}^2{\tilde{\sigma }}_m^{\pm }}$ 。基于光滑指示子定义非线性权

${\omega }_m^{\pm }=\frac{{\tilde{\omega }}_m^{\pm }}{{\displaystyle {\sum }_{l=0}^2{\tilde{\omega }}_l^{\pm }}},{\omega }_l^{\pm }=\frac{{\sigma }_l^{\pm }}{{\left(\epsilon +{\beta }_l\right)}^2},m,l=0,1,2.$

本文取$\epsilon ={10}^{-6}$ 。将正负项线性组合，得到最终的非线性权

${\omega }_m={\gamma }^{+}{\omega }_m^{+}-{\gamma }^{-}{\omega }_m^{-},m=0,1,2.$

从而得

${\widehat{l}}_{i+\frac{1}{2}}={\displaystyle \sum _{m=0}^2{\omega }_m{p}_{i+\frac{1}{2}}^{\left(m\right)}}.$

2.3. HWENO格式对导函数方程的重构

类似地，在与2.2节相同的小模板点上使用函数点值及导函数点值构造埃尔米特四阶多项式，得到${p^{\prime }}^{\left(0\right)}\left(x_{i+\frac{1}{2}}\right),{p^{\prime }}^{\left(1\right)}\left(x_{i+\frac{1}{2}}\right),{p^{\prime }}^{\left(2\right)}\left(x_{i+\frac{1}{2}}\right)$ ，此时线性权为${d^{\prime }}_0=\frac{3}{7},{d^{\prime }}_1=\frac{3}{7},{d^{\prime }}_2=\frac{1}{7}$ ，从而得到最终估计

${\widehat{z}}_{i+\frac{1}{2}}={\displaystyle \sum _{m=0}^2{{\omega }^{\prime }}_m{p^{\prime }}_{i+\frac{1}{2}}^{\left(m\right)}}.$

2.4. Mapped非线性权重

对上述非线性权进行精度分析，

$\begin{array}{l}{\displaystyle \sum _{m=0}^2\left({\omega }_m-d_m\right)p_{i\pm \frac{1}{2}}^{\left(m\right)}}={\theta }_{i\pm \frac{1}{2}}{\displaystyle \sum _{m=0}^2\left({\omega }_m-d_m\right)}+\Delta x^6{\displaystyle \sum _{m=0}^2{A}_m\left({\omega }_m-d_m\right)}+{\displaystyle \sum _{m=0}^2\left({\omega }_m-d_m\right){\rm O}}\left(\Delta x^7\right),\\ {\displaystyle \sum _{m=0}^2\left({{\omega }^{\prime }}_m-{d^{\prime }}_m\right){p^{\prime }}_{i\pm \frac{1}{2}}^{\left(m\right)}}={{\theta }^{\prime }}_{i\pm \frac{1}{2}}{\displaystyle \sum _{m=0}^2\left({{\omega }^{\prime }}_m-{d^{\prime }}_m\right)}+\Delta x^5{\displaystyle \sum _{m=0}^2{B}_m\left({{\omega }^{\prime }}_m-{d^{\prime }}_m\right)}+{\displaystyle \sum _{m=0}^2\left({{\omega }^{\prime }}_m-{d^{\prime }}_m\right){\rm O}}\left(\Delta x^6\right),\end{array}$

由于${\displaystyle \sum _{m=0}^2{\omega }_m}={\displaystyle \sum _{m=0}^2{d}_m}=1$ ，${\displaystyle \sum _{m=0}^2{{\omega }^{\prime }}_m}={\displaystyle \sum _{m=0}^2{d^{\prime }}_m}=1$ ，为了在光滑区域确保原函数和导函数均达到六阶精度，根据精度分析，非线性权需满足以下条件：

$\begin{array}{l}{\omega }_m-d_m={\rm O}\left(\Delta x^2\right),\\ {{\omega }^{\prime }}_m-{d^{\prime }}_m={\rm O}\left(\Delta x^3\right),\end{array}$ (5)

对光滑指示子进行精度分析，将其带入非线性权，得到

$\begin{array}{c}{\sigma }_m^{\pm }={\omega }_m^{\pm }\left({\displaystyle \sum _{k=0}^2\frac{{\sigma }_m^{\pm }}{{\left(\epsilon +{\beta }_k\right)}^2}}\right){\left(\epsilon +{\beta }_k\right)}^2\\ ={\omega }_m^{\pm }\left(\frac{1}{D^2}\left(1+{\rm O}\left(\Delta x^2\right)\right)\right){\left(D\left(1+{\rm O}\left(\Delta x^2\right)\right)\right)}^2\\ ={\omega }_m^{\pm }+{\rm O}\left(\Delta x^2\right),m=0,1,2.\end{array}$

即

$\begin{array}{l}{\omega }_m-d_m={\rm O}\left(\Delta x^2\right),\\ {{\omega }^{\prime }}_m-{d^{\prime }}_m={\rm O}\left(\Delta x^2\right).\end{array}$

发现导函数方程的非线性权不满足精度条件(5)，利用特殊处理方法[17]对非线性权进行处理

$g_m\left({\omega }^{\prime }\right)=\frac{{\omega }^{\prime }\left({d^{\prime }}_m+{d^{\prime }}_m^2-3{d^{\prime }}_m{\omega }^{\prime }+{{\omega }^{\prime }}^2\right)}{{d^{\prime }}_m^2+{\omega }^{\prime }\left(1-2{d^{\prime }}_m\right)},m=0,1,2.$

取${{\delta }^{\prime }}_m=g_m\left({{\omega }^{\prime }}_m\right)$ ，定义最终非线性权

${{\omega }^{\prime }}_m^{(M)}=\frac{{{\delta }^{\prime }}_m}{{\displaystyle {\sum }_{m=0}^2{{\delta }^{\prime }}_m}},m=0,1,2.$ (6)

导函数估计的最终形式为

${\widehat{z}}_{i+\frac{1}{2}}={\displaystyle \sum _{m=0}^2{{\omega }^{\prime }}_m^{\left(M\right)}{p^{\prime }}_{i+\frac{1}{2}}^{\left(m\right)}}$

对映射后的非线性权进行精度分析

$\begin{array}{c}{{\delta }^{\prime }}_m=g_m\left({d^{\prime }}_m\right)+{g^{\prime }}_m\left({d^{\prime }}_m\right)\left({{\omega }^{\prime }}_m-{d^{\prime }}_m\right)+\frac{{g^{″}}_m\left({d^{\prime }}_m\right)}{2}{\left({{\omega }^{\prime }}_m-{d^{\prime }}_m\right)}^2+\cdots \\ ={d^{\prime }}_m+{\rm O}\left(\Delta x^3\right)+\frac{{g^{‴}}_m\left({d^{\prime }}_m\right)}{24}{\rm O}\left(\Delta x^4\right),\end{array}$

根据$\frac{1}{{\displaystyle {\sum }_{m=0}^2{{\delta }^{\prime }}_m}}=\frac{1}{1+{\rm O}\left(\Delta x^3\right)}=1+{\rm O}\left(\Delta x^3\right)$ ，将其代入(6)得到

$\begin{array}{c}{{\omega }^{\prime }}_m^{\left(M\right)}=\left({d^{\prime }}_m+{\rm O}\left(\Delta x^3\right)+\frac{{g^{‴}}_m\left({d^{\prime }}_m\right)}{24}{\rm O}\left(\Delta x^4\right)\right)\left(1+{\rm O}\left(\Delta x^3\right)\right)\\ ={d^{\prime }}_m+{\rm O}\left(\Delta x^3\right)+\frac{{g^{‴}}_m\left({d^{\prime }}_m\right)}{24}{\rm O}\left(\Delta x^4\right).\end{array}$

通过分析，映射后的非线性权能够满足精度要求(5)，因此可以保持HWENO格式在不连续区域的非振荡性，并在连续区域达到六阶精度，相关验证结果将在后续的数值试验中给出。

3. Lax-Wendroff时间离散

设时间步长为$\Delta t$ ，$u^n$ 表示第n时间层，$u^{\left(r\right)}$ 表示u的r阶时间导数，则$u,v$ 在时间方向上的Taylor展开式为

$\begin{array}{l}{u}^{n+1}=u^n+\Delta t{u}^{\prime }+\frac{\Delta t^2}{2}{u}^{″}+\cdots +\frac{\Delta t^r}{r!}{u}^{\left(r\right)}+{\rm O}\left(\Delta t^{r+1}\right),\\ v^{n+1}=v^n+\Delta t{v}^{\prime }+\frac{\Delta t^2}{2}{v}^{″}+\cdots +\frac{\Delta t^r}{r!}{v}^{\left(r\right)}+{\rm O}\left(\Delta t^{r+1}\right).\end{array}$

上述展开式的精度可以进行选择，本文取$r=4$ ，可参考文献[14]。Lax-Wendroff时间离散方法的构造过程如下：

第一步，利用HWENO格式求解得到一阶导数$u^{\prime },v^{\prime }$ 。

第二步，先通过偏微分方程求解二阶导数，再通过插值进行近似，得到：

$\begin{array}{l}{u^{″}}_i\approx -\frac{1}{12\Delta x^2}\left({\varphi }_{i-2}-16{\varphi }_{i-1}+30{\varphi }_i-16{\varphi }_{i+1}+{\varphi }_{i+2}\right)\\ {v^{″}}_i\approx -\frac{1}{12\Delta x^2}\left({\phi }_{i-2}-16{\phi }_{i-1}+30{\phi }_i-16{\phi }_{i+1}+{\phi }_{i+2}\right).\end{array}$

其中，

$\begin{array}{l}\varphi =l^{\prime }\left(u\right)u^{\prime },\\ \phi =l^{″}\left(u\right)v{u}^{\prime }+l^{\prime }\left(u,v\right)v^{\prime }.\end{array}$

第三步，类似的，对于三阶导数：

$\begin{array}{l}{u^{‴}}_i\approx -\frac{1}{12\Delta x^2}\left({\varphi }_{i-2}-16{\varphi }_{i-1}+30{\varphi }_i-16{\varphi }_{i+1}+{\varphi }_{i+2}\right)\\ {v^{‴}}_i\approx -\frac{1}{12\Delta x^2}\left({\phi }_{i-2}-16{\phi }_{i-1}+30{\phi }_i-16{\phi }_{i+1}+{\phi }_{i+2}\right)\end{array}$

其中，

$\begin{array}{l}\varphi =l^{\prime }\left(u\right)u^{″}+l^{″}\left(u\right){\left(u^{\prime }\right)}^2,\\ \phi =l^{‴}\left(u\right)v{\left(u^{\prime }\right)}^2+l^{″}\left(u\right)\left(2v^{\prime }{u}^{\prime }+v{u}^{″}\right)+l^{\prime }\left(u\right)v^{″}.\end{array}$

第四步，对于四阶导数：

$\begin{array}{l}{u}_i^{\left(4\right)}\approx -\frac{1}{12\Delta x^2}\left({\varphi }_{i-2}-16{\varphi }_{i-1}+30{\varphi }_i-16{\varphi }_{i+1}+{\varphi }_{i+2}\right)\\ v_i^{\left(4\right)}\approx -\frac{1}{12\Delta x^2}\left({\phi }_{i-2}-16{\phi }_{i-1}+30{\phi }_i-16{\phi }_{i+1}+{\phi }_{i+2}\right)\end{array}$

其中，

$\begin{array}{l}\varphi =l^{\prime }\left(u\right)u^{‴}+3l^{″}\left(u\right)u^{″}{u}^{\prime }+l^{‴}\left(u\right){\left(u^{\prime }\right)}^3,\\ \phi =l^{‴}\left(u\right)v^{\prime }{\left(u^{\prime }\right)}^2+2l^{‴}\left(u\right)v^{\prime }{u}^{\prime }{u}^{″}+l^{(4)}\left(u\right){\left(u^{\prime }\right)}^{3}v+2\left(l^{″}\left(u\right)v^{\prime }{u}^{″}+l^{″}\left(u\right)v^{″}{u}^{\prime }+l^{‴}\left(u\right){\left(u^{\prime }\right)}^2{v}^{\prime }\right)\\ +l^{″}\left(u\right)v{u}^{‴}+l^{″}\left(u\right)v^{\prime }{u}^{″}+l^{‴}\left(u\right)v^{\prime }{u}^{\prime }{u}^{″}+l^{\prime }\left(u\right)v^{‴}+l^{″}\left(u\right)v^{″}{u}^{\prime }.\end{array}$

利用此方法，可以得到时空同步的数值离散格式。其原理与Lax-Wendroff格式类似，即通过空间导数替代泰勒展开中的时间导数，因此称为Lax-Wendroff时间离散。理论上，这种离散格式能够实现任意阶精度，但随着精度的提升，计算量也会相应增加。

4. 数值试验

本节利用大量数值算例验证格式的精度及稳定性。为了确保在精度测试中空间误差占主导地位，我们选取适合的CFL数为0.35，进行数值试验使用的计算机硬件配置：联想AMD Ryzen 5 2500U。

4.1. 精度验证

我们首先验证六阶有限差分HWENO格式在求解周期边界热传导方程初值问题时的精度。

$\{\begin{array}{l}{u}_t=u_{xx},\hfill \\ u\left(x,0\right)=\sin \left(x\right),-\pi \le x<\pi \hfill \end{array}$ (7)

方程(7)的精确解为

$u^e\left(x,t\right)={\text{e}}^{-t}\sin \left(x\right).$

在表1中，我们展示了$T=2$ 时使用HWENO-LW与HWENO-RK格式求解的数值结果，表2则给出了两种格式的计算时间。结果表明，两种格式都能达到预期的精度阶，且计算精度相近。然而，根据表2中的结果，HWENO-LW格式的计算时间比HWENO-RK节省了约$\frac{1}{3}$ ，说明Lax-Wendroff格式在计算效率上更具优势。

Table 1. The order of accuracy of the HWENO scheme

表1. HWENO格式的精度阶

	HWENO-LW				HWENO-RK
N	$L^1$ error	order	$L^{\infty }$ error	order	$L^1$ error	order	$L^{\infty }$ error	order
10	1.79E−05		2.75E−05		3.36E−06		5.18E−06
20	3.03E−07	5.89	4.79E−07	5.84	5.25E−08	6.00	8.31E−08	5.96
40	4.79E−09	5.98	7.55E−09	5.99	8.09E−10	6.02	1.27E−09	6.03
80	7.51E−11	6.00	1.18E−10	6.00	1.26E−11	6.00	1.98E−11	6.01
160	1.16E−12	6.02	1.82E−12	6.02	1.98E−13	5.99	3.11E−13	5.99

Table 2. The computational time of the HWENO scheme

表2. HWENO格式的计算时间

	N = 10	N = 20	N = 40	N = 80	N = 160
HWENO-LW	0.015625	0.046875	0.781250	2.046875	10.546875
HWENO-RK	0.046875	0.125000	1.625000	4.734375	13.062500

4.2. 对流扩散Buckley-Leverett方程

研究标量对流扩散方程即Buckley-Leverett方程[18]

$u_t+F{\left(u\right)}_x=\epsilon {\left(\tau \left(u\right)u_x\right)}_x,\epsilon \tau \left(u\right)\ge 0.$

该模型常用来模拟地下石油勘探中的两相流。

这里我们取参数$\epsilon =0.01$ ，非凸的流通量函数$F\left(u\right)$

$F\left(u\right)=\frac{u^2}{u^2+{\left(1-u\right)}^2},$

和

$\tau \left(u\right)=\{\begin{array}{ll}4u\left(1-u\right),\hfill & 0\le u\le 1,\hfill \\ 0,\hfill & 其他.\hfill \end{array}$

初值条件为

$u\left(x,0\right)=\{\begin{array}{cc}0,& 0\le x\le 1-\frac{1}{\sqrt{2}},\\ 1,& 1-\frac{1}{\sqrt{2}}\le x\le 1,\end{array}$

$u\left(0,t\right)=1$ 为边界条件。

对于对流项，采用五阶有限差分HWENO格式和Lax-Friedrichs数值流通量分裂。图1给出了$T=0.2$ ，计算区域$[0,1]$ 被分成$N=100$ 个一致网络时的数值解，并与[18]的数值结果进行对比，结果表明HWENO-LW与HWENO-RK格式均能精确地收敛到正确的熵解，并且在间断处保持了无振荡特性。从表3的运行时间可知，HWENO-LW格式的计算效率更高。

Table 3. The computational time of the HWENO scheme for solving the Buckley-Leverett equation

表3. HWENO格式求解Buckley-Leverett方程的解的运行时间

时间	HWENO-LW	HWENO-RK
N = 100	15.254 s	20.568 s

Figure 1. The solution of the Buckley-Leverett equation

图1. Buckley-Leverett方程的解

4.3. 强退化对流扩散方程

考虑强退化对流扩散方程[18]

$u_t+F{\left(u\right)}_x=\epsilon {\left(\tau \left(u\right)u_x\right)}_x,\epsilon \tau \left(u\right)\ge 0.$

参数$\epsilon =0.1$ ，$F\left(u\right)=u^2$ ，及

$\tau \left(u\right)=\{\begin{array}{cc}0,& \left|u\right|\le 0.25,\\ 1,& \left|u\right|>0.25,\end{array}$

初值条件为

$u\left(x,0\right)=\{\begin{array}{ll}1,\hfill & -\frac{1}{\sqrt{2}}-0.4<x<-\frac{1}{\sqrt{2}}+0.4,\hfill \\ -1,\hfill & \frac{1}{\sqrt{2}}-0.4<x<\frac{1}{\sqrt{2}}+0.4,\hfill \\ 0,\hfill & 其他.\hfill \end{array}$

计算区域为$\left[-2,2\right]$ 。图2给出了$T=0.7$ ，网格数取$N=100$ 时的数值计算结果。£表示HWENO-LW格式数值计算结果，×表示HWENO-RK格式的计算结果，与数值结果[18]一致，并且在间断处是本质无振荡的。

Figure 2. The solution of the Riemann problem

图2. 黎曼问题的解

4.4. 二维多孔介质方程

对于高维问题在实际实现中，可以逐维分别进行HWENO离散。这里我们数值模拟二维多孔介质方程[19]

$u_t={\left(u^2\right)}_{xx}+{\left(u^2\right)}_{yy},$

其中初值条件$u\left(x,y,0\right)$ 为两个凸包

$u\left(x,y,0\right)=\{\begin{array}{ll}{\text{e}}^{\frac{-1}{6-{\left(x-2\right)}^2-{\left(y+2\right)}^2},}\hfill & {\left(x-2\right)}^2+{\left(y+2\right)}^2<6,\hfill \\ {\text{e}}^{\frac{-1}{6-{\left(x+2\right)}^2-{\left(y-2\right)}^2},}\hfill & {\left(x+2\right)}^2+{\left(y-2\right)}^2<6,\hfill \\ 0,\hfill & 其他.\hfill \end{array}$

计算区域为$\left[-10,10\right]\times \left[-10,10\right]$ 。图3中从左上到右下分别给出了$T=0,0.5,1.0,4.0$ 时的数值解，其中计算区域分成80 × 80的均匀网络，数值结果同[19]中一致。

Figure 3. The solution of the two-dimensional porous media equation

图3. 二维多孔介质方程的解

5. 结论

本文采用Lax-Wendroff时间离散的有限差分HWENO格式来求解非线性退化抛物方程。具体来说，HWENO方法用于处理空间导数项，而通过将空间导数替换为Taylor展开中的时间导数项，构造了Lax-Wendroff时间离散方法，从而得到了HWENO-LW格式。数值实验表明，与传统的Runge-Kutta时间离散方法相比，Lax-Wendroff方法在保持精度的同时显著提高了计算效率，证明了该方法的有效性和工程应用前景。

基金项目

该工作得到了中国山西省自然科学基金(第202103021224041号)的部分支持。