考虑二维线性传输方程
$\left\{ \begin{array}{l} {u_t} + a{u_x} + b{u_y} = 0,\\ u\left( {x,y,0} \right) = {u_0}\left( {x,y} \right),\forall x \in \mathit{\Omega }\mathit{.} \end{array} \right.$ | (1) |
这里a和b是常数.一维线性传输方程为
$\left\{ \begin{array}{l} {u_t} + a{u_x} = 0,\\ u\left( {x,0} \right) = {u_0}\left( x \right),\forall x \in \mathit{\Omega }\mathit{.} \end{array} \right.$ | (2) |
线性传输方程是一类很重要的双曲型守恒型方程,很多数值方法的建立都是从该方程开始.一些计算工作者针对双曲守恒型方程提出了大量的数值格式,如Godunov[1],全变差下降(TVD)格式[2],间断伽略金(DG)格式[3],本质非振荡(ENO)格式[4],加权本质非振荡(WENO)格式[5]等.所有这些格式为了格式的稳定性都适当地引入了数值粘性,但是数值粘性会带来数值磨损.如何在保证格式的稳定性的条件下引入最少的数值粘性一直是非常有挑战性的课题.1987年,Tadmor[6]发展了一类二阶的熵守恒格式,并得到了一个对数值粘性更精确的量化方法,得出一个三点格式只需含有比熵守恒格式更多的数值粘性则是熵稳定的.2003年,Tadmor[7]进一步发展了该方法,使得该方法适用于各种守恒系统.随后,Tadmor和Zhong[8]将该方法应用到了Navier-Stokes方程和浅水波方程.2006年,Roe将熵守恒格式和著名的Roe格式结合起来得到的格式[9]不仅满足熵稳定条件还能很好地计算激波,该格式只有一阶精度.国内封建湖团队发展了大量的高分辨率的熵相容格式,刘友琼等[10]针对浅水波方程发展了熵相容格式, 程晓晗等[11]发展了双曲守恒律方程的WENO型熵相容格式.
在理论上,标量双曲型守恒型方程的解满足最大值原理,原始的ENO格式,WENO格式以及DG格式都不满足最大值原理.针对一维标量双曲型方程满足最大值原理的三阶非振荡格式最早由Liu提出[12],随后Zhang等[13]将它推广到高阶,发展了高阶精度满足最大值原理的DG格式和WENO格式.
近几年来,茅德康等[14]针对双曲守恒性方程提出了一种新的Godunov型的熵格式.熵格式的基本思想是在数值离散中同时离散熵等多个物理量的方程,由此就同时从多个侧面来数值地刻画物理过程,实验证明该方法可以有效地控制数值粘性和长时间的计算.与其他方法不同的是,熵格式不仅计算了数值解还计算了数值熵,计算中每一步有数值解和数值熵两个实体,该格式通过数值熵来控制数值粘性.非常巧妙的是,即使是对线性方程,该格式也是非线性的.崔艳芬和茅德康[15]对熵格式进行了数值分析,揭示了这种格式在计算中各步误差相互抵消这一性质.但是熵格式不满足最大值原理,文献[16]的数值实验可以看出数值解会出现过高过低现象,长时间计算后还会出现非物理振荡.
本文主要是针对一维和二维线性传输方程提出满足最大值原理的熵格式,该格式可以看作是满足最大值原理的格式和熵格式的一个结合.具体做法是,通过对熵的多项式重构引入最大值原理限制器.数值实验表明满足最大值原理的熵格式集中了熵格式和最大值原理格式的优势,不仅适合长时间计算而且不会出现过高过低的现象,能消除熵格式产生的非物理振荡.
1 一维线性传输方程满足最大值原理熵格式的描述任意给定的非线性函数U(u),一维线性传输方程(2) 的光滑解式(3) 也成立:
$\left\{ \begin{array}{l} U{\left( u \right)_t} + aU{\left( u \right)_x} = 0,\\ u\left( {x,0} \right) = {u_0}\left( x \right),\forall x \in \mathit{\Omega }\mathit{.} \end{array} \right.$ | (3) |
如果取U(u)为凸函数,那么它就是熵.为了方便描述,取U(u) = u2.初值问题(2) 有通解如下:
$u\left( {x,t} \right) = {u_0}\left( {x - at} \right)$ | (4) |
记
$M = \mathop {\max }\limits_x {u_0}\left( x \right),$ | (5) |
$m = \mathop {{\rm{Min}}}\limits_x {u_0}\left( x \right),$ | (6) |
分别是初始值的最大值和最小值.如果对于任意的x和t,有u(x; t)∈[m; M],那么u(x, t)满足最大值原理.
首先将区间Ω剖分为单元集合, Ω=∪jIj,其中
为了叙述的简单,网格步长取为固定值h.引入数值解为对tn时刻精确解网格平均的近似,
${u_{j,n}} \approx \frac{1}{h}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {u\left( {x,{t_n}} \right){\rm{d}}x} $ | (7) |
另外引入了一个数值熵,它是对精确解的熵的网格平均的逼近,
${U_{j,n}} \approx \frac{1}{h}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {U\left( {u\left( {x,{t_n}} \right)} \right){\rm{d}}x} $ | (8) |
熵格式包含数值解uj, n和数值熵Uj, n两个数值实体.重构函数和un和Un都有关.
熵格式是Godunov型的格式,按重构,发展和网格平均三步来进行.
(1) 重构.首先给出tn时刻每一个网格Ij内熵重构解,它是一个线性函数,
${R_{\rm{E}}}\left( {x;{u_n},{U_n}} \right) = {u_{j,n}} + {s_{j,n,e}}\left( {x - {x_j}} \right),$ | (9) |
这里sj, n, e称为重构斜率.熵重构(8) 自动满足:
$\frac{1}{h}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {{R_{\rm{E}}}\left( {x;{u_n},{U_n}} \right){\rm{d}}x} = {u_{j,n}},$ | (10) |
为了计算sj, n, e,要求
$\frac{1}{h}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {U\left( {{R_{\rm{E}}}\left( {x;{u_n},{U_n}} \right)} \right){\rm{d}}x} = {U_{j,n}}$ | (11) |
即重构解的熵的网格平均等于该网格的数值熵.将式(8) 代入式(10),不难得到:
${\left( {{s_{j,n,{\rm{e}}}}} \right)^2} = \frac{{12\left( {{U_{j,n}} - {{\left( {{u_{j,n}}} \right)}^2}} \right)}}{{{h^2}}}.$ | (12) |
可以求出sj, n, e,符号取为和(uj+1,n-uj-1,n)一致.为了使得格式满足最大值原理,得对熵重构函数RE(x; un, Un), 加上限制器,最终的重构为
$\begin{array}{l} {R_{{\rm{ME}}}}\left( {x;{u_n},{U_n}} \right) = \theta \left( {{R_{{\rm{ME}}}}\left( {x;{u_n},{U_n}} \right) - {u_{j,n}}} \right) + {u_{j,n}},\\ \theta = {\rm{min}}\left\{ {\left| {\frac{{{M_0} - {u_{j,n}}}}{{M' - {u_{j,n}}}}} \right|,\left| {\frac{{{m_0} - {u_{j,n}}}}{{m' - {u_{j,n}}}}} \right|,1} \right\}, \end{array}$ | (13) |
这里M′=maxx∈IjRE(x; un, Un), m′=minx∈IjRE(x; un, Un).详细参见文献[13].
(2) 发展.以重构函数RME(x; un, Un)作为tn层的初值,求解以下初值问题:
$\left\{ \begin{array}{l} {v_t} + a{v_x} = 0, - \infty < x < \infty ,{t_n} < t < {t_{n + 1}},\\ v\left( {x,{t_n}} \right) = {R_{{\rm{ME}}}}\left( {x;{u_n},{U_n}} \right), - \infty < x < \infty . \end{array} \right.$ | (14) |
定解问题(13) 的精确解为v(x, t)=RME(x-at; un, Un).
(3) 网格平均.在t=tn+1时的数值解和数值熵分别计算为
$u_j^{n + 1} = \frac{1}{h}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {v\left( {x,{t_{n + 1}}} \right){\rm{d}}x} ,$ | (15) |
$U_j^{n + 1} = \frac{1}{h}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {U\left( {v\left( {x,{t_{n + 1}}} \right)} \right){\rm{d}}x} .$ | (16) |
二维的满足最大值原理的熵格式是将文献[16]中的重构函数引入一个最大值原理限制器,下面的格式描述沿用文献[14]的基本思路.和一维情形类似,任意给定的非线性函数U(u),二维线性传输方程(1) 的光滑解式(17) 也成立:
$\left\{ \begin{array}{l} U{\left( u \right)_t} + aU{\left( u \right)_x} + bU{\left( u \right)_y} = 0,\\ u\left( {x,y,0} \right) = {u_0}\left( {x,y} \right),\forall x \in \mathit{\Omega }\mathit{.} \end{array} \right.$ | (17) |
为了描述的简单,取U(u)=u2.初始函数u0(x, y)的最大值和最小值分别记为M和m,即:
$M = \mathop {{\rm{Max}}}\limits_{x,y} {u_0}\left( {x,y} \right),$ | (18) |
$m = \mathop {{\rm{Min}}}\limits_{x,y} {u_0}\left( {x,y} \right),$ | (19) |
如果对于任意的x,y和t,有u(x, y, t)∈[m; M],那么u(x, y, t)满足最大值原理.
对区域Ω作四边形剖分,x方向和y方向的网格步长都取为h,网格记为
${u_{j,n,k}} \approx \frac{1}{{{h^2}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {u\left( {x,y,{t_n}} \right){\rm{d}}x} } {\rm{d}}\mathit{y,}$ | (20) |
${U_{j,n,k}} \approx \frac{1}{{{h^2}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {U\left( {u\left( {x,y,{t_n}} \right)} \right){\rm{d}}x} } {\rm{d}}\mathit{y,}$ | (21) |
这里u(x, y, tn)表示方程在tn时的的精确解.
如同一维一样,熵格式是Godunov型的格式,按重构,发展和网格平均三步来进行.
(1) 重构.在tn层上熵重构函数RE(x, y; un, Un)在网格Ij, k上取为x,y的一次函数,
${R_{\rm{E}}}\left( {x,\mathit{y};{u_n},{U_n}} \right) = {u_{j,n,k}} + \alpha \left( {{{\left( {{s_x}} \right)}_{j,n,\mathit{k}}}\left( {x - {x_j}} \right) + {{\left( {{s_y}} \right)}_{j,n,\mathit{k}}}\left( {y - {y_k}} \right)} \right),$ | (22) |
其中
$\frac{1}{{{h^2}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {{R_{\rm{E}}}\left( {x,y;{\mathit{u}_n},{\mathit{U}_n}} \right){\rm{d}}x} } {\rm{d}}\mathit{y,} = {u_{j,n,k}},$ | (23) |
为了确定α,要求重构熵函数还满足如下守恒关系:
$\frac{1}{{{h^2}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {U\left( {{R_{\rm{E}}}\left( {x,y;{\mathit{u}_n},{\mathit{U}_n}} \right)} \right){\rm{d}}x} } {\rm{d}}\mathit{y,} = {U_{j,n,k}}$ | (24) |
由RE(x, y; un, Un)自动满足式(23).可通过式(24) 解出,
${\alpha ^2} = \frac{{12\left( {{U_{j,n}} - {{\left( {{u_{j,n}}} \right)}^2}} \right)}}{{{h^2}\left( {{{\left( {{{\left( {{s_x}} \right)}_{j,n,k}}} \right)}^2} + {{\left( {{{\left( {{s_y}} \right)}_{j,n,k}}} \right)}^2}} \right)}}$ | (25) |
此时可以开根号求得α,符号取正的,为1附近的一个正数.为了使得格式满足最大值原理,得对熵重构函数RE(x, y; un, Un)加上限制器,最终的重构为
$\theta \left( {{R_{{\rm{ME}}}}\left( {x,y;{u_n},{U_n}} \right) - {u_{j,n,k}}} \right) + {u_{j,n,k}},$ | (26) |
$\theta = {\rm{min}}\left\{ {\left| {\frac{{{M_0} - {u_{j,n,k}}}}{{M' - {u_{j,n,k}}}}} \right|,\left| {\frac{{{m_0} - {u_{j,n,k}}}}{{m' - {u_{j,n,k}}}}} \right|,1} \right\}$ | (27) |
这里M′=max(x, y)∈Ij, kRE(x, y; un, Un), m′=min(x, y)∈Ij, kRE(x, y; un, Un).详细参见文献[13].
(2) 发展.求解初值问题:
$\left\{ \begin{array}{l} {v_t} + a{v_x} + b{v_y} = 0,\\ - \infty < x,y < \infty ,{t_n} < t < {t_{n + 1}},\\ v\left( {x,y,{t_n}} \right) = {R_{{\rm{ME}}}}\left( {x,y;{u_n},{U_n}} \right), - \infty < x,y < \infty . \end{array} \right.$ | (28) |
它的精确解为v(x, y, t)=RME(x-at, y-bt; un, Un).
(3) 网格平均. t=tn+1时刻的数值解和数值熵的计算分别为
${u_{j,n + 1,k}} = \frac{1}{{{h^2}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {v\left( {x,y,{t_{n + 1}}} \right){\rm{d}}x{\rm{d}}y,} } $ | (29) |
${U_{j,n + 1,k}} = \frac{1}{{{h^2}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {{{\left( {v\left( {x,y,{t_{n + 1}}} \right)} \right)}^2}{\rm{d}}x{\rm{d}}y.} } $ | (30) |
在实际计算中,采用通量的形式
${u_{j,n + 1,k}} = {u_{j,n,k}} - \lambda \left( {{{\hat f}_{j + 1/2,n,k}} + {{\hat f}_{j - 1/2,n,k}}} \right) - \lambda \left( {{{\hat g}_{j,n,k + 1/2}} + {{\hat g}_{j,n,k - 1/2}}} \right),$ | (31) |
${U_{j,n + 1,k}} = {U_{j,n,k}} - \lambda \left( {{{\hat F}_{j + 1/2,n,k}} + {{\hat F}_{j - 1/2,n,k}}} \right) - \lambda \left( {{{\hat G}_{j,n,k + 1/2}} + {{\hat G}_{j,n,k - 1/2}}} \right),$ | (32) |
这里τ是时间步长,
${{\hat f}_{j \pm 1/2,n,k}} = \frac{1}{{h\tau }}\int_{{t_n}}^{{t_{n + 1}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {av\left( {{x_{j \pm 1/2}},y,t} \right){\rm{d}}y{\rm{d}}t,} } $ | (33) |
${{\hat g}_{j,n,k \pm 1/2}} = \frac{1}{{h\tau }}\int_{{t_n}}^{{t_{n + 1}}} {\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {bv\left( {x,{y_{k \pm 1/2}},t} \right){\rm{d}}\mathit{x}{\rm{d}}t,} } $ | (34) |
${{\hat F}_{j \pm 1/2,n,k}} = \frac{1}{{h\tau }}\int_{{t_n}}^{{t_{n + 1}}} {\int_{{y_{k - 1/2}}}^{{y_{k + 1/2}}} {a{{\left( {v\left( {{x_{j \pm 1/2}},y,t} \right)} \right)}^2}{\rm{d}}y{\rm{d}}t,} } $ | (35) |
${{\hat G}_{j,n,k \pm 1/2}} = \frac{1}{{h\tau }}\int_{{t_n}}^{{t_{n + 1}}} {\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {b{{\left( {v\left( {x,{y_{k \pm 1/2}},t} \right)} \right)}^2}{\rm{d}}\mathit{x}{\rm{d}}t,} } $ | (36) |
式(31)~(34) 的积分采用Gauss积分公式进行计算.
3 数值实验计算了文献[14]4个经典算例,对文献[14]中的熵格式和本文发展的满足最大值原理的熵格式的计算结果进行比较.库朗数(CFL)取为0.5,“+”表示熵格式,“○”表示满足最大值原理的熵格式.
算例1 考虑初值问题(2),
$u\left( {x,0} \right) = \left\{ \begin{array}{l} \exp \left\{ { - \frac{1}{{1 - 16{{\left( {x - \frac{1}{2}} \right)}^2}}}} \right\},x \in \left[ {\frac{1}{4},\frac{3}{4}} \right),\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x \notin \left[ {\frac{1}{4},\frac{3}{4}} \right). \end{array} \right.$ | (37) |
边界条件取周期边界条件.图 1为网格数为200,计算时间分别为t=20和t=200熵格式和满足最大値原理的熵格式与精确解的比较图.从图 1可以看出熵格式长时间计算会出现非物理振荡,而且会出现过高过低现象,而满足最大值原理的熵格式不会产生非物理振荡也没有出现过高过低现象.
算例2 考虑初值问题(2),
$u\left( {x,0} \right) = \exp \left\{ { - 100 \cdot {{\left( {x - \frac{1}{2}} \right)}^2}} \right\}\sin \left( {80x} \right),0 \le x < 1.$ | (38) |
周期边界条件.这个算例比算例1的机构要复杂得多.图 2为网格数为200,计算时间分别为t=20和t=200熵格式和满足最大値原理的熵格式与精确解的比较图.从图 2不难看出熵格式的数值解不能保证在初值的最大值和最小值之间,而满足最大值原理的熵格式的数值解一直保持在最大值和最小值之间.
算例3 考虑初值问题(1),
$u\left( {x,y,0} \right) = \left\{ \begin{array}{l} \exp \left\{ { - \frac{1}{{1 - 16{{\left( {x - \frac{1}{2}} \right)}^2}}}} \right\} \times \exp \left\{ { - \frac{1}{{1 - 16{{\left( {y - \frac{1}{2}} \right)}^2}}}} \right\},\\ x \in \left[ {\frac{1}{4},\frac{3}{4}} \right) \times \left[ {\frac{1}{4},\frac{3}{4}} \right),\\ 0,x \notin \left[ {\frac{1}{4},\frac{3}{4}} \right) \times \left[ {\frac{1}{4},\frac{3}{4}} \right). \end{array} \right.$ | (39) |
周期边界条件.图 3的左边是网格数为200,计算时间t=20时熵格式和满足最大値原理格式以及精确解与平面y=0.505的截面图的比较图.图 3的右边是数值解的三维立体图.从图 3可以看出熵格式产生了过高过低现象而且有非物理振荡产生,而满足最大值原理的熵格式没有出现非物理振荡而且其数值解一直保持在初值的最大值和最小值之间.
算例4 考虑方程(1),初值为:
$\begin{array}{l} u\left( {x,y,0} \right) = \exp \left\{ { - 1000 \cdot {{\left( {x - \frac{1}{2}} \right)}^2}} \right\}\\ \begin{array}{*{20}{l}} {\sin \left( {80x} \right) \cdot \exp \left\{ { - 100 \times {{\left( {y - \frac{1}{2}} \right)}^2}} \right\}\sin \left( {80y} \right),}\\ {\left( {x,y} \right) \in \left[ {0,1} \right) \times \left[ {0,1} \right).} \end{array} \end{array}$ | (40) |
周期性边界条件.图 4的左边是网格数为200,计算时间t=10时熵格式和满足最大値原理格式以及精确解与平面y=0.505的截面图的比较图.图 4的右边是数值解的三维立体图.从图 4还可以看出,熵格式的数值解磨损得非常厉害,而满足最大值原理的熵格式保持结构较好.
本文针对线性传输方程发展了满足最大值原理的熵格式.利用熵格式和满足最大值原理的熵格式对一维和二维线性传输方程进行了数值实验,数值实验表明满足最大值原理的熵格式不会出现非物理振荡,数值解保持在初值的最大值和最小值之间,比熵格式具有更好的保结构性质.
[1] |
GODUNOV S K. A finite difference method for computation of discontinuous solutions of the equations of fluid dynamics[J]. Matematicheskii Sbornik, 1959, 47(11): 357 |
[2] |
HARTEN A, OSHER S. Uniformly high-order accurate non-oscillatory scheme Ⅱ[J]. SIAM Journal on Numerical Analysis, 1987, 24(2): 279 DOI:10.1137/0724022 |
[3] |
COCKBURN B, SHU C W. TVB Runge-Kutta local projecting discontinuous Galerkin finite element methods for conservation laws Ⅱ: general frame work[J]. Mathematics of Computation, 1989, 52(186): 411 |
[4] |
SHU C W, OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes Ⅱ[J]. Journal of Computational Physics, 1989, 83(2): 32 |
[5] |
JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1): 202 DOI:10.1006/jcph.1996.0130 |
[6] |
TADMOR E. Numerical viscosity of entropy stable schemes for systems of conservation laws[J]. Mathematics of Computation, 1987, 49(179): 91 DOI:10.1090/S0025-5718-1987-0890255-3 |
[7] |
TADMOR E. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[J]. Acta Numerica, 2003, 12(12): 451 |
[8] |
TADMOR E, ZHONG W. Entropy stable approximations of navier-stokes equations with no artificial numerical viscosity[J]. Journal of Hyperbolic Differential Equations, 2006, 3(3): 529 DOI:10.1142/S0219891606000896 |
[9] |
CHENG X H, NIE Y F, FENG J H, et al. Self-adjusting entropy-stable scheme for compressible Euler equations[J]. Chinese Physics B, 2015, 24(2): 16 |
[10] |
刘友琼, 封建湖, 梁楠, 等. 求解浅水波方程的熵相容格式[J]. 应用数学与力学, 2013, 34(12): 1247 LIU Youqiong, FENG Jianhu, LIANG Nan, et al. An entropy-consistent flux scheme for shallow water equations[J]. Applied Mathematics and Mechanics, 2013, 34(12): 1247 |
[11] |
程晓晗, 封建湖, 聂玉峰. 求解双曲守恒律方程的WENO型熵相容格式[J]. 爆炸与冲击, 2014, 34(4): 501 CHENG Xiaohan, FENG Jianhu, NIE Yufeng. WENO type entropy consistent scheme for hyperbolic conservation laws[J]. Explosion & Shock Waves, 2014, 34(4): 501 DOI:10.11883/1001-1455(2014)04-0501-07 |
[12] |
OSHER S, LIU X D. Non-oscillatory high order accurate self similar maximum principle satisfying shock capturing schemes[J]. SIAM Journal on Numerical Analysis, 1996, 33(2): 760 DOI:10.1137/0733038 |
[13] |
ZHANG X X, SHU C W. On maximum-principle-satisfying high order schemes for scalar conservation laws[J]. Journal of Computational Physics, 2010, 229(9): 3091 DOI:10.1016/j.jcp.2009.12.030 |
[14] |
LI H X, WANG Z G, MAO D K. Numerically neither dissipative nor compressive scheme for linear advection equation and its application to the Euler system[J]. Journal of Scientific Computing, 2008, 36(3): 285 DOI:10.1007/s10915-008-9192-x |
[15] |
CUI Y F, MAO D K. Error self-canceling of a difference scheme maintaining two conservation laws for linear advection equation[J]. Mathematics of Computation, 2012, 81(278): 715 |
[16] |
CHEN R S, ZOU M, XIAO L. Entropy-TVD scheme for the shallow water equations in one dimension[J]. Journal of Scientific Computing, 2017, 71(2): 822 DOI:10.1007/s10915-016-0322-6 |