# Exponential approximation space reconstruction WENO scheme for dispersive
PDEs

Lavanya V Salian§ and Samala Rathan†

###### Abstract.

In this work, we construct a fifth-order weighted essentially non-oscillatory
(WENO) scheme with exponential approximation space for solving dispersive
equations. A conservative third-order derivative formulation is developed
directly using WENO spatial reconstruction procedure and third-order TVD
Runge-Kutta scheme is used for the evaluation of time derivative. This
exponential approximation space consists a tension parameter that may be
optimized to fit the specific feature of the charecteristic data, yielding
better results without spurious oscillations compared to the polynomial
approximation space. A detailed formulation is presented for the the
construction of conservative flux approximation, smoothness indicators,
nonlinear weights and verified that the proposed scheme provides the required
fifth convergence order. One and two-dimensional numerical examples are
presented to support the theoretical claims

§Department of Humanities and Sciences, Indian Institute of Petroleum and
Energy-Visakhapatnam, India-530003<EMAIL_ADDRESS>
†Department of Humanities and Sciences, Indian Institute of Petroleum and
Energy-Visakhapatnam, India-530003<EMAIL_ADDRESS>

AMS Classification: 41A10, 65M06

Keywords: Finite difference, WENO scheme, exponential polynomial, dispersion,
order of accuracy.

##  1\. Introduction

In this paper, a conservative finite difference fifth-order weighted
essentially non-oscillatory (WENO) scheme is designed using exponential
polynomials to approximate the third-order derivative term, which evolved in
the nonlinear, possibly prototype dispersion equations of the form,

(1.1)
$\begin{split}u_{t}+f(u)_{x}+g(u)_{xxx}&=0,\quad(x,t)\in\Omega\times(0,T],\\\
u(x,0)&=u_{0}(x),\end{split}$

where $f(u)$ and $g(u)$ are arbitrary (may smooth) nonlinear functions. The
simplest form, where dispersion and nonlinearity occur, is the famous Korteweg
de-Vries (KdV) equation,

(1.2) $u_{t}+auu_{x}+u_{xxx}=0,\quad(x,t)\in\Omega\times(0,T],\\\ $

where $u_{t}$ is the time evolution of the wave propagating in one direction,
$uu_{x}$ is the nonlinearity accounting for the steeping of the wave and
$u_{xxx}$ is the linear dispersion, i.e. spreading of the wave. Solitary waves
are localized traveling waves that maintain a constant speed and shape,
diminishing to asymptotic zero amplitude at a considerable distance [1]. The
solutions for solitary waves are postulated in the form of $u(x,t)=k(x-ct)$,
with $c$ representing the propagation speed of the wave and $k(z)$,
$k^{{}^{\prime}}(z)$, $k^{{}^{\prime\prime}}(z)$ $\to 0$ as $z\to\pm\infty$,
$z=x-ct$. Pioneering work by Zabusky and Kruskal [2] unveiled that solitary
waves engage in nonlinear interactions, resulting in an emerging wave that
conserves both its shape and amplitude. These distinctive solitary waves are
denoted as solitons. There is an increase in research efforts worldwide across
a variety of scientific fields to investigate the concept of solitons. This
concept of solitons has engendered amplified global research efforts spanning
diverse scientific domains. Notably, disciplines encompassing fluid dynamics,
astrophysics, plasma physics, magneto-acoustic waves, and numerous others have
been captivated by the exploration of solitons due to their pronounced
significance. This is owing to the genesis of solitons from the intricate
equilibrium between weak nonlinearity and dispersion.

Later, a class of solitary waves with compact support that is termed as
compactons were discovered by Rosenau and Hyman [3]. The general form of the
Rosenau - Hyman equation, also known as $K(n,n)$ equation, is of the form

(1.3) $u_{t}+(u^{n})_{x}+(u^{n})_{xxx}=0,\quad n>1,\\\ $

where $u(x,t)$ is the wave amplitude as a function of the spatial variable $x$
and time $t$. Compactons represent a distinct subclass of solitons,
distinguished by their finite wavelength, absence of exponential tails, lack
of infinite wings, and resilient adherence to soliton-like behavioral
patterns. Compactons and soliton solutions of the KdV equation exhibit shared
attributes. For instance, the velocity of an individual compacton is directly
proportional to its amplitude. During the movement and interaction of multiple
compactons with varying velocities, nonlinear effects come into play,
resulting in an altered phase upon exit. In contrast to conventional KdV
solitons, the width of a compacton remains constant regardless of its
amplitude. The $K(n,n)$ equation diverges from the typical energy conservation
laws held by the KdV equation and can only be derived from a first-order
Lagrangian for $n=1$. While solitons are analytical solutions, compactons are
non-analytical solutions distinguished by non-analytical points at their
edges. These points of non-analyticity align with genuine nonlinearity points
within the differential equations.

The resolution of the Rosenau-Hyman equation is challenging due to the
concurrent interplay of dispersion effects and nonlinearity. Although the
extensively employed pseudo-spectral method [4] in spatial domains effectively
preserves solution positivity and incorporates high-pass filters to induce
artificial dissipation, the post-compacton collision can lead to sign
alterations within the solution. A variety of alternate methodologies have
been explored in this context. These include the finite difference method with
Pade approximation [5], the local discontinuous Galerkin method [6], the
second-order finite difference approach [7], the adaptive mesh refinement-
based line method [8], and and the direct WENO scheme employing polynomial
bases for dispersion-type equations [9]. The governing equation dictating
dispersive waves, as denoted by (1.1), bears significant resemblances to
hyperbolic conservation laws. Conspicuously, both equation classes are
susceptible to sharp fronts and wavefronts propagating at finite velocities.
Consequently, an avenue of inquiry involves extending numerical techniques
originally devised for resolving hyperbolic conservation laws, such as the
WENO technique, to contend with the intricacies inherent to the dispersive
wave equation. However, this extension mandates a meticulous adaptation of the
WENO procedure to ensure the preservation of conservation, accuracy, and
nonoscillatory behavior.

In 1994, Liu et al.[10] introduced the original version of the WENO scheme
within the context of a finite volume framework, employed for the resolution
of one-dimensional conservation law equations. Subsequently, in 1996, Jiang
and Shu [11] introduced an enhanced rendition of the WENO scheme within the
finite differences framework, exhibiting greater efficiency compared to the
method presented by Liu et al.[10] for addressing both one-dimensional and
two-dimensional conservation law equations. Despite the evident enhancements
introduced by Jiang and Shu [11], their methodology exhibited certain
limitations, primarily manifesting in instances where gradients of higher
order became negligible, thereby resulting in a diminution of accuracy order.
Numerous researchers have endeavored to mitigate the numerical dissipation
near discontinuities and to optimize the computational efficiency of the
conventional WENO-JS scheme through various adaptations. The authors of [16],
Henrick et. al., proposed a fifth-order WENO method named WENO-M. Notably,
WENO-M employs a mapping technique to minimize the deviation of nonlinear
weights utilized in the convex combination of stencils from the optimal
weights, barring instances involving pronounced discontinuities. A distinct
variation of the fifth-order WENO scheme, known as WENO-Z, was proposed by
Borges et al. in their work [17]. The weighting formulation employed in this
context diverges slightly from that of the WENO-JS formulation. In [20], Ha et
al. presented a new smoothness indicator, referred to as WENO-NS, which
utilizes the $L^{1}$-sum of generalized undivided differences to approximate
the derivatives of flux functions. This indicator allows for achieving fifth-
order convergence even in smooth regions, including critical points where the
first derivatives are zero. Further, Rathan and Naga Raju enhanced the
accuracy of WENO scheme at critical points for fifth and seventh order schemes
[21, 22, 23]. Various researchers have proposed different high-order WENO
schemes with the aim of enhancing the efficiency of these numerical schemes
for solving hyperbolic conservation laws [18, 19, 24]. In 2016, Ha et al.[24]
presented a novel WENO scheme utilizing exponential polynomials for solving
hyperbolic conservation laws.

The comparative analysis of numerical outcomes between methodologies based
upon exponential, trigonometric and algebraic polynomial constructions for
hyperbolic conservation laws [24, 25, 26, 27] and non-linear degenerate
parabolic equations [28, 29, 30] has been studied in literature. In instances
involving interpolation of data manifesting rapid gradients or high
oscillations, the utilization of exponential or trigonometric polynomial bases
confers a superior degree of efficiency than algebraic polynomial bases. These
alternative polynomial bases prove to be more suitable for accurately
capturing and representing such complex features in the data. The primary
objective of this study is to propose a novel fifth-order Weighted Essentially
Non-Oscillatory scheme named WENO-E, which employs exponential polynomials for
solving the dispersion equation. The key design strategy of the WENO-E scheme
is to attain the highest possible approximation order in smooth regions while
ensuring accuracy is maintained even at critical points. By utilizing
exponential polynomials, the scheme aims to enhance the accuracy and
robustness of the numerical solution, enabling more effective and reliable
computations for the dispersion equation. A global smoothness indicator based
on generalized undivided differences is introduced to aid in the design of the
nonlinear weights that play a crucial role in WENO reconstructions. Numerical
experiments are conducted and compared with the polynomial-based scheme
(WENO-Z) to demonstrate the WENO-E scheme’s ability to accurately approximate
solutions near singularities. This work is the first of its kind for solving
prototype dispersive equations (1.1) with non polynomial approximation space,
as per the authors’ best knowledge.

The paper’s organization is delineated as follows: Section 2 elaborates on a
broad framework for the finite differences WENO scheme and introduces the
concept of approximating numerical fluxes with exponential polynomials. This
approach is utilized to effectively address the third-order derivative term
inherent in the prototype dispersion equation. In Section 3, we introduce an
innovative technique that leverages exponential polynomials to construct
numerical fluxes, both for the large stencil and its sub-stencils. We also
discuss the associated ideal weights. In Section 4, the paper outlines the use
of $L^{1}$-norm smoothness indicators to construct non-linear weights and
provides an analysis of their accuracy in both smooth regions and at critical
points. To support our theoretical claim, we present few one and two
dimensional numerical results in Section 5. In Section 6, we provide brief
concluding remarks.

##  2\. Finite difference WENO scheme

In this section, we describe a general framework of finite difference WENO
schemes based on exponential polynomials to solve prototype dispersion
equations. Without loss of generality, we shall focus on the one-dimensional
prototype dispersion equations of the form

(2.1)
$\begin{split}u_{t}+f(u)_{x}+g(u)_{xxx}&=0,\quad(x,t)\in\Omega\times(0,T],\\\
u(x,0)&=u_{0}(x),\end{split}$

along with periodic boundary conditions. To extend the algorithm to higher
dimensions, the 1-D algorithm is applied along each coordinate direction.
Assume the uniform spatial mesh as follows:

(2.2) $\left\\{\begin{array}[]{cc}x\in[x_{l},x_{r}],\quad
x_{i}=x_{l}+(i-1)\Delta x,\quad i=1:N,\quad,x_{1}=x_{l},\quad x_{N}=x_{r},\\\
I_{i}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}],\quad
x_{i+\frac{1}{2}}=\frac{(x_{i+1}+x_{i})}{2},\end{array}\right.$

where $\Delta x$ is the spatial step size and $u_{i}^{n}$ is defined as a
nodal point value $u(x_{i},t^{n})$.

###  2.1. Formulation of WENO scheme

We use a conservative finite difference scheme in a method of lines (MOL)
approach to writing

$\frac{du_{i}}{dt}=-\frac{1}{\Delta
x}({\hat{F}_{i+\frac{1}{2}}}-{\hat{F}_{i-\frac{1}{2}}})-\frac{1}{\Delta
x^{3}}({\hat{G}_{i+\frac{1}{2}}}-{\hat{G}_{i-\frac{1}{2}}}),$

where $\hat{F}$ and $\hat{G}$ are the numerical flux for convection and
dispersion, respectively. The detailed construction procedure for convection
flux is given in [11]. For the dispersion term, the $k^{\text{th}}$-order
accurate conservative finite difference scheme is as follows:

(2.3) $\frac{1}{\Delta
x^{3}}({\hat{G}_{i+\frac{1}{2}}}-{\hat{G}_{i-\frac{1}{2}}})\approx
g(u)_{xxx}\big{|}_{x=x_{i}}+\mathcal{O}(\Delta x^{k}),$

where ${\hat{G}_{i+\frac{1}{2}}}$ is a numerical dispersion flux at the cell
boundary $x_{i+\frac{1}{2}}$,
${\hat{G}_{i+\frac{1}{2}}}={\hat{G}(u_{i-r},\dots,u_{i+s})}$ of $(r+s)$
variables. It has to satisfy Lipschitz continuity in each of its arguments and
is consistent with the physical flux ${\hat{G}(u,\dots,u)}=g(u)$. Ahmat et al.
in [9] implicitly considered a function $h(x)$ as follows to guarantee the
conservative property

(2.4) $g(u)=\frac{1}{\Delta x^{3}}\int_{x-\frac{\Delta x}{2}}^{x+\frac{\Delta
x}{2}}\int_{\eta-\frac{\Delta x}{2}}^{\eta+\frac{\Delta
x}{2}}\int_{\zeta-\frac{\Delta x}{2}}^{\zeta+\frac{\Delta
x}{2}}h(\theta)\,d\theta d\zeta d\eta,$

then by triple derivation, we have

(2.5) $g(u)_{xxx}=\frac{h(x+\frac{3}{2}\Delta x)-3h(x+\frac{1}{2}\Delta
x)+3h(x-\frac{1}{2}\Delta x)-h(x-\frac{3}{2}\Delta x)}{\Delta x^{3}}.$

If we define a function $G(x)$, such that

(2.6) $G(x)=h(x+\Delta x)-2h(x)+h(x-\Delta x),$

then we have

(2.7)
$g(u)_{xxx}\big{|}_{x=x_{i}}=\frac{G(x_{i+\frac{1}{2}})-G(x_{i-\frac{1}{2}})}{\Delta
x^{3}}.$

If we have numerical flux ${\hat{G}_{i+\frac{1}{2}}}$ which is an
approximation to $G(x_{i+\frac{1}{2}})$ up to $k^{\text{th}}$-order, then we
have equation (2.3).

We now derive the specific form of the numerical flux ${\hat{G}(u)}$ for
WENO5. We first split the flux into positive and negative parts, that is
$g(u)=g^{+}(u)+g^{-}(u)$ with $\frac{\partial g^{+}(u)}{\partial u}\geq 0$ and
$\frac{\partial g^{-}(u)}{\partial u}\leq 0$, and $G^{+}(x)$ and $G^{-}(x)$
are defined by (2.6) according to $g^{+}(u)$ and $g^{-}(u)$, respectively. The
reconstructed numerical fluxes ${\hat{G}^{+}_{i+\frac{1}{2}}}$ and
${\hat{G}^{-}_{i+\frac{1}{2}}}$ to approximate $G^{+}(x_{i+\frac{1}{2}})$ and
$G^{-}(x_{i+\frac{1}{2}})$ up to $k^{\text{th}}$-order, respectively. Finally,
we define numerical fluxes ${\hat{G}_{i+\frac{1}{2}}}$ as

(2.8)
${\hat{G}_{i+\frac{1}{2}}}={\hat{G}^{+}_{i+\frac{1}{2}}}+{\hat{G}^{-}_{i+\frac{1}{2}}},$

and ${\hat{G}_{i+\frac{1}{2}}}$ is an approximation to $G(x_{i+\frac{1}{2}})$
up to $k^{\text{th}}$-order. The reconstruction procedure for
${\hat{G}^{+}_{i+\frac{1}{2}}}$ is given below and procedure for
${\hat{G}^{-}_{i+\frac{1}{2}}}$ is mirror symmetric respect to
$x_{i+\frac{1}{2}}$.

Figure 1.  The numerical flux ${\hat{G}^{+}_{i+\frac{1}{2}}}$ is constructed
on the stencil $S_{7}=\\{x_{i-2},\dots,x_{i+4}\\}$ with seven uniform points
and three 5-point substencils $S_{0}$, $S_{1}$, $S_{2}$.

###  2.2. Approximation using exponential polynomials

The WENO scheme for solving the dispersion type equation involves
approximating the value of $G(x_{i+\frac{1}{2}})$ using equation (2.7). The
reconstruction process must balance achieving optimal accuracy in smooth
regions while maintaining essential non-oscillatory behaviour in non smooth
regions. Polynomials are widely used to construct numerical fluxes. However,
the limitation of using polynomials is that the approximation space is shift-
and-scale invariant. Thus, it cannot be tailored to suit the specific
characteristics of the given data. This limitation can result in significant
numerical dissipation when interpolating data with rapid gradients, which
hinders the ability to produce sharp edges. To overcome this limitation,
researchers [24, 25, 26, 27, 28] have explored the use of other types of basis
functions, such as exponential polynomials, which have been shown to yield
better results in terms of producing sharp edges and reducing numerical
dissipation.
The general form of exponential polynomials can be written as:

(2.9) $\Phi(x)=x^{n}e^{\lambda x},$

where $n$ is a non-negative integer and $\lambda\in\mathbb{R}$ or
$\lambda\in\iota\mathbb{R}$ $(\iota^{2}=-1)$. Let
$\\{\Phi_{1},\dots,\Phi_{r}\\}$ be a set of exponential polynomials of the
form in (2.9) and let $\Gamma_{r}$ be the space defined by

(2.10) $\Gamma_{r}:=\text{span}\\{\Phi_{1},\dots,\Phi_{r}\\}.$

A necessary condition for the space $\Gamma_{r}$ is that the set
$\\{\Phi_{1},\dots,\Phi_{r}\\}$ is linearly independent so that the
determinants of the Wronskian matrix related to them are non-zero, i.e.,

(2.11) $\det(\Phi_{n}(s_{i}):i,n=1,\dots,r)\neq 0,$

for any $r$-point stencil $\\{s_{i}:i=1,\dots,r\\}$. The space $\Gamma_{r}$
needs to satisfy the following basic requirements for the practical
computation of the proposed interpolation:

  * •

Shift-invariant. The space $\Gamma_{r}$ should be shift-invariant in the sense
that for any $\alpha\in\mathbb{R}$, $f\in\Gamma_{r}$ implies
$f(\cdot-\alpha)\in\Gamma_{r}$. This ensures that the interpolation kernel is
invariant under the shifting of the evaluation location and stencil. It allows
for a set of interpolation kernels to be precomputed for a fixed point and
then applied to every evaluation position based on the chosen stencil at a
given cell boundary. The space $\Gamma_{r}$ defined in (2.11) meets this
requirement since it is shift-invariant.

  * •

Symmetry. The space $\Gamma_{r}$ should be symmetry means that if a function
$f$ is in the space $\Gamma_{r}$, then the reflected function $f(-\cdot)$
should also belong to the same space.

Including the polynomial $\Phi(x)=1$ in the space $\Gamma_{r}$ ensures that
the sum of the interpolation weights over all basis functions is equal to 1,
which is necessary for the interpolation kernel to satisfy the partition of
unity property. In this paper, we choose

(2.12) $\Gamma_{7}:=\text{span}\\{1,x,x^{2},e^{\lambda x},e^{-\lambda
x},\cos{\lambda x},\sin{\lambda x}\\}$

as the basis functions for global stencil $S_{7}:=\\{x_{i-2},\dots,x_{i+4}\\}$
and similarly,

(2.13) $\Gamma_{5}=\text{span}\\{1,x,x^{2},e^{\lambda x},e^{-\lambda x}\\}$

for the five-point substencils $S_{m}:=\\{x_{i-2+m},\dots,x_{i+2+m}\\}$,
$m=0,1,2$. Here, $\Gamma_{7}$ and $\Gamma_{5}$ constitute an extended
Tchebysheff systems on $\mathbb{R}$ so that the non-singularity of the
interpolation matrices in (2.11) is guaranteed as in [31, 24].

##  3\. Fifth-order WENO-E scheme

Conservative numerical schemes are developed by approximating the function
$G(x)$ in (2.7). This approximation, represented as $\hat{G}(x)$, is
constructed using a exponential polynomial form with unspecified coefficients.
When this polynomial is substituted into (2.4), it results in a system of
equations. In this system, the flux is known at the nodes surrounding the
relevant interface, enabling the determination of a distinct set of
coefficients. After obtaining $\hat{G}(x)$, the approximation of the spatial
derivative in (2.3) is as follows

(3.1) $g(u)_{xxx}\big{|}_{x=x_{i}}\approx\frac{1}{\Delta
x^{3}}({\hat{G}_{i+\frac{1}{2}}}-{\hat{G}_{i-\frac{1}{2}}}).$

We will focus on two orders of convergence. The primary concern is the order
at which (3.1) is satisfied, as it dictates the spatial convergence rate of
the overall scheme. Additionally, we consider the order of individual
approximations to the numerical flux, ${\hat{G}_{i\pm\frac{1}{2}}}$, which is
significant in establishing criteria for the acceptance of non-oscillatory
weights.

###  3.1. Constructions of numerical flux ${\hat{G}^{+}_{i+\frac{1}{2}}}$

The approximations of ${\hat{G}^{+}_{i+\frac{1}{2}}}$ are denoted by $p(x)$.
We first consider an exponential polynomial approximation to $h(x)$ on the
7-point stencil $S_{7}$

(3.2) $h(x)\approx q(x)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}e^{\lambda
x}+a_{4}e^{-\lambda x}+a_{5}\cos{\lambda x}+a_{6}\sin{\lambda x},$

with undetermined coefficients $a_{k}$ with $k=0,\dots,6$.
Substituting (3.2) into (2.4)

(3.3) $g^{+}(u(x))=\frac{1}{\Delta x^{3}}\int_{x-\frac{\Delta
x}{2}}^{x+\frac{\Delta x}{2}}\int_{\eta-\frac{\Delta x}{2}}^{\eta+\frac{\Delta
x}{2}}\int_{\zeta-\frac{\Delta x}{2}}^{\zeta+\frac{\Delta
x}{2}}q(\theta)\,d\theta d\zeta d\eta,$

and performing the integration gives

(3.4) $\begin{split}g^{+}(u(x))&=a_{0}+a_{1}x+a_{2}\biggl{(}\frac{\Delta
x^{2}}{4}+x^{2}\biggr{)}+a_{3}\biggl{(}-\frac{4e^{\lambda
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}+\frac{4e^{\lambda x}\cosh[\lambda\Delta
x]\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ &+a_{4}\biggl{(}\frac{e^{\frac{\lambda}{2}(\Delta
x-2x)}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{-3\lambda\Delta x}{2}-\lambda x}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{-\lambda
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ &+a_{5}\biggl{(}\frac{8\cos[\lambda
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}+a_{6}\biggl{(}\frac{8\sin[\lambda
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}.\end{split}$

To determine the coefficients $a_{0},\dots,a_{6}$, one can consider (3.4) as
$g^{+}(u(x_{i}-2\Delta x))=g^{+}(u_{i-2}),\dots,g^{+}(u(x_{i}+4\Delta
x))=g^{+}(u_{i+4})$ with $x_{i}=0$ and solve the resulting $7\times 7$ system
$AX=B$ are specified in Appendix:B. Substituting the coefficients of vector
$X$ into (3.2) and then calculating $p(x)=q(x+\Delta x)-2q(x)+q(x-\Delta x)$
at $x=x_{i+\frac{1}{2}}$, we get

(3.5)
${\hat{G}^{+}_{i+\frac{1}{2}}}=C_{0}g^{+}(u_{i-2})+C_{1}g^{+}(u_{i-1})+C_{2}g^{+}(u_{i})+C_{3}g^{+}(u_{i+1})+C_{4}g^{+}(u_{i+2})+C_{5}g^{+}(u_{i+3})+C_{6}g^{+}(u_{i+4}).$

After applying the Taylor series to the coefficients $C_{j},0\leq j\leq 6$ of
equation (3.5) is given by

(3.6) $\begin{split}C_{0}&=-\frac{1}{15}-\frac{229\lambda^{4}\Delta
x^{4}}{226800}+\mathcal{O}(\Delta x^{8}),\quad
C_{1}=\frac{21}{40}+\frac{1493\lambda^{4}\Delta
x^{4}}{75600}+\mathcal{O}(\Delta x^{8}),\\\
C_{2}&=\frac{1}{8}-\frac{587\lambda^{4}\Delta x^{4}}{15120}+\mathcal{O}(\Delta
x^{8}),\quad\quad C_{3}=-\frac{23}{12}-\frac{65\lambda^{4}\Delta
x^{4}}{9072}+\mathcal{O}(\Delta x^{8}),\\\
C_{4}&=\frac{7}{4}+\frac{11\lambda^{4}\Delta x^{4}}{378}+\mathcal{O}(\Delta
x^{8}),\quad\quad C_{5}=-\frac{19}{40}-\frac{323\lambda^{4}\Delta
x^{4}}{18900}+\mathcal{O}(\Delta x^{8}),\\\
C_{6}&=-\frac{7}{120}+\frac{103\lambda^{4}\Delta
x^{4}}{113400}+\mathcal{O}(\Delta x^{8}).\end{split}$

Thus, the equation (3.5) becomes

$\displaystyle{\hat{G}^{+}_{i+\frac{1}{2}}}=$
$\displaystyle\biggl{(}-\frac{1}{15}-\frac{229\lambda^{4}\Delta
x^{4}}{226800}\biggr{)}g^{+}(u_{i-2})+\biggl{(}\frac{21}{40}+\frac{1493\lambda^{4}\Delta
x^{4}}{75600}\biggr{)}g^{+}(u_{i-1})+\biggl{(}\frac{1}{8}-\frac{587\lambda^{4}\Delta
x^{4}}{15120}\biggr{)}g^{+}(u_{i})$
$\displaystyle+\biggl{(}-\frac{23}{12}-\frac{65\lambda^{4}\Delta
x^{4}}{9072}\biggr{)}g^{+}(u_{i+1})+\biggl{(}\frac{7}{4}+\frac{11\lambda^{4}\Delta
x^{4}}{378}\biggr{)}g^{+}(u_{i+2})$
$\displaystyle+\biggl{(}-\frac{19}{40}-\frac{323\lambda^{4}\Delta
x^{4}}{18900}\biggr{)}g^{+}(u_{i+3})+\biggl{(}-\frac{7}{120}+\frac{103\lambda^{4}\Delta
x^{4}}{113400}\biggr{)}g^{+}(u_{i+4}).$

From the definition of (2.4) we also know that

$\begin{split}{\hat{G}^{+}_{i+\frac{1}{2}}}&=\frac{1}{\Delta
x^{3}}\Biggl{[}\biggl{(}-\frac{1}{15}-\frac{229\lambda^{4}\Delta
x^{4}}{226800}\biggr{)}\biggl{(}\int_{x_{i-\frac{5}{2}}}^{x_{i-\frac{3}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta d\eta\biggr{)}\\\
&+\biggl{(}\frac{21}{40}+\frac{1493\lambda^{4}\Delta
x^{4}}{75600}\biggr{)}\biggl{(}\int_{x_{i-\frac{3}{2}}}^{x_{i-\frac{1}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta d\eta\biggr{)}\\\
&+\biggl{(}\frac{1}{8}-\frac{587\lambda^{4}\Delta
x^{4}}{15120}\biggr{)}\biggl{(}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta d\eta\biggr{)}\\\
&+\biggl{(}-\frac{23}{12}-\frac{65\lambda^{4}\Delta
x^{4}}{9072}\biggr{)}\biggl{(}\int_{x_{i+\frac{1}{2}}}^{x_{i+\frac{3}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta d\eta\biggr{)}\\\
&+\biggl{(}\frac{7}{4}+\frac{11\lambda^{4}\Delta
x^{4}}{378}\biggr{)}\biggl{(}\int_{x_{i+\frac{3}{2}}}^{x_{i+\frac{5}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta d\eta\biggr{)}\\\
&+\biggl{(}-\frac{19}{40}-\frac{323\lambda^{4}\Delta
x^{4}}{18900}\biggr{)}\biggl{(}\int_{x_{i+\frac{5}{2}}}^{x_{i+\frac{7}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta d\eta\biggr{)}\\\
&+\biggl{(}-\frac{7}{120}+\frac{103\lambda^{4}\Delta
x^{4}}{113400}\biggr{)}\biggl{(}\int_{x_{i+\frac{7}{2}}}^{x_{i+\frac{9}{2}}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}h(\theta)\,d\theta d\zeta
d\eta\biggr{)}\Biggr{]}.\end{split}$

To analyze the accuracy of the approximation to the third derivative for
smooth solutions, we assume $h(x)$ has sufficient local regularity.
Substituting the Taylor series expansion at the grid point $x_{i}$,

$h(\theta)=h(x_{i})+\sum_{j=1}^{7}\frac{(\theta-
x_{i})^{j}}{j!}\frac{d^{j}h}{dx^{j}}\bigg{|}_{x=x_{i}}+\mathcal{O}(\Delta
x^{8}),$

and performing the integration, we have

(3.7)
$\begin{split}{\hat{G}^{+}_{i+\frac{1}{2}}}&=\frac{d^{2}h}{dx^{2}}\bigg{|}_{x=x_{i}}\Delta
x^{2}+\frac{1}{2}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}\Delta
x^{3}+\frac{5}{24}\frac{d^{4}h}{dx^{4}}\bigg{|}_{x=x_{i}}\Delta
x^{4}+\frac{1}{16}\frac{d^{5}h}{dx^{5}}\bigg{|}_{x=x_{i}}\Delta
x^{5}+\frac{91}{5760}\frac{d^{6}h}{dx^{6}}\bigg{|}_{x=x_{i}}\Delta x^{6}\\\
&+\biggl{(}-\frac{7}{240}\lambda^{4}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}+\frac{25}{768}\frac{d^{7}h}{dx^{7}}\bigg{|}_{x=x_{i}}\biggr{)}\Delta
x^{7}+\mathcal{O}(\Delta x^{8}).\end{split}$

Comparing (3.7) with Taylor series expansion

(3.8)
$\begin{split}G(x_{i+\frac{1}{2}})=&h(x_{i+\frac{3}{2}})-2h(x_{i+\frac{1}{2}})+h(x_{i-\frac{1}{2}})\\\
=&\frac{d^{2}h}{dx^{2}}\bigg{|}_{x=x_{i}}\Delta
x^{2}+\frac{1}{2}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}\Delta
x^{3}+\frac{5}{24}\frac{d^{4}h}{dx^{4}}\bigg{|}_{x=x_{i}}\Delta
x^{4}+\frac{1}{16}\frac{d^{5}h}{dx^{5}}\bigg{|}_{x=x_{i}}\Delta x^{5}\\\
&+\frac{91}{5760}\frac{d^{6}h}{dx^{6}}\bigg{|}_{x=x_{i}}\Delta
x^{6}+\frac{13}{3840}\frac{d^{7}h}{dx^{7}}\bigg{|}_{x=x_{i}}\Delta
x^{7}+\mathcal{O}(\Delta x^{8}),\end{split}$

we obtain

(3.9)
$\begin{split}{\hat{G}^{+}_{i+\frac{1}{2}}}&=G(x_{i+\frac{1}{2}})+\biggl{(}-\frac{7}{240}\lambda^{4}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}+\frac{7}{240}\frac{d^{7}h}{dx^{7}}\bigg{|}_{x=x_{i}}\biggr{)}\Delta
x^{7}+\mathcal{O}(\Delta x^{8}),\\\
&=G(x_{i+\frac{1}{2}})+\mathcal{A}^{+}\Delta x^{7}+\mathcal{O}(\Delta
x^{8}).\end{split}$

Similarly, we calculate the numerical flux ${\hat{G}^{+}_{i-\frac{1}{2}}}$

(3.10)
$\begin{split}{\hat{G}^{+}_{i-\frac{1}{2}}}&=\biggl{(}-\frac{1}{15}-\frac{229\lambda^{4}\Delta
x^{4}}{226800}\biggr{)}g^{+}(u_{i-3})+\biggl{(}\frac{21}{40}+\frac{1493\lambda^{4}\Delta
x^{4}}{75600}\biggr{)}g^{+}(u_{i-2})+\biggl{(}\frac{1}{8}-\frac{587\lambda^{4}\Delta
x^{4}}{15120}\biggr{)}g^{+}(u_{i-1})\\\
&+\biggl{(}-\frac{23}{12}-\frac{65\lambda^{4}\Delta
x^{4}}{9072}\biggr{)}g^{+}(u_{i})+\biggl{(}\frac{7}{4}+\frac{11\lambda^{4}\Delta
x^{4}}{378}\biggr{)}g^{+}(u_{i+1})+\biggl{(}-\frac{19}{40}-\frac{323\lambda^{4}\Delta
x^{4}}{18900}\biggr{)}g^{+}(u_{i+2})\\\
&+\biggl{(}-\frac{7}{120}+\frac{103\lambda^{4}\Delta
x^{4}}{113400}\biggr{)}g^{+}(u_{i+3})\\\
{\hat{G}^{+}_{i-\frac{1}{2}}}&=G(x_{i-\frac{1}{2}})+\biggl{(}-\frac{7}{240}\lambda^{4}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}+\frac{7}{240}\frac{d^{7}h}{dx^{7}}\bigg{|}_{x=x_{i}}\biggr{)}\Delta
x^{7}+\mathcal{O}(\Delta x^{8})=G(x_{i-\frac{1}{2}})+\mathcal{A}^{-}\Delta
x^{7}+\mathcal{O}(\Delta x^{8}),\end{split}$

where the values $\mathcal{A}^{\pm}$ are independent of $\Delta x$, but
dependent on $\lambda$. Therefore, substituting (3.9) and (3.10) into (2.7),
we have the fifth order approximation

(3.11)
$g(u)_{xxx}\big{|}_{x=x_{i}}=\frac{{\hat{G}^{+}_{i+\frac{1}{2}}}-{\hat{G}^{+}_{i-\frac{1}{2}}}}{\Delta
x^{3}}+\mathcal{O}(\Delta x^{5}).$

###  3.2. Constructions of numerical flux ${\hat{G}^{m}_{i+\frac{1}{2}}}$ for
$m=0,1,2$

When the stencil $S_{7}$ contains a discontinuity, the numerical flux
${\hat{G}^{+}_{i+\frac{1}{2}}}$ formulation on the big stencil $S_{7}$ can
result in oscillations because all seven nodes are affected. To mitigate this
issue, the WENO procedure considers smaller stencils
$S_{m}=\\{x_{i-2+m},\dots,x_{i+2+m}\\}$, $m=0,1,2$, of the big stencil
$S_{7}$. Following a similar argument, on each small stencils $S_{m},$
$m=0,1,2$, by taking approximation to $h(x)$ as

(3.12) $h(x)\approx
q_{m}(x)=a_{0}^{m}+a_{1}^{m}x+a_{2}^{m}x^{2}+a_{3}^{m}e^{\lambda
x}+a_{4}^{m}e^{-\lambda x},$

with undetermined coefficients $a_{k}^{m}$ with $k=0,\dots,4$, $m=0,1,2$.
Substituting (3.12) into (2.4)

(3.13) $\displaystyle\frac{1}{\Delta x^{3}}\int_{x_{j}-\frac{\Delta
x}{2}}^{x_{j}+\frac{\Delta x}{2}}\int_{\eta-\frac{\Delta
x}{2}}^{\eta+\frac{\Delta x}{2}}\int_{\zeta-\frac{\Delta
x}{2}}^{\zeta+\frac{\Delta x}{2}}q_{m}(\theta)\,d\theta d\zeta d\eta$
$\displaystyle=g^{+}(u_{j}),\quad j=i-2+m,\dots,i+2+m,\quad m=0,1,2,$

and performing the integration gives

(3.14)
$\begin{split}g^{+}(u(x))&=a_{0}^{m}+a_{1}^{m}x+a_{2}^{m}\biggl{(}\frac{\Delta
x^{2}}{4}+x^{2}\biggr{)}+a_{3}^{m}\biggl{(}\frac{e^{\frac{-3\lambda\Delta
x}{2}+x\lambda}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}+a_{4}^{m}\biggl{(}\frac{e^{\frac{-3\lambda\Delta x}{2}-\lambda
x}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\biggr{)}.\end{split}$

To determine the coefficients $a_{0}^{m},\dots,a_{4}^{m}$ we solve the
resulting $5\times 5$ systems $A_{m}X_{m}=B_{m}$ is given by,

  1. (1)

For $m=0$,

$\displaystyle A_{0}=\left[\begin{matrix}1&-2\Delta x&\frac{17}{4}\Delta
x^{2}&\frac{e^{-\frac{7\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&-\Delta
x&\frac{5}{4}\Delta x^{2}&\frac{e^{-\frac{5\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}&\frac{e^{-\frac{\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&0&\frac{1}{4}\Delta
x^{2}&\frac{e^{-\frac{3\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{3\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&\Delta
x&\frac{5}{4}\Delta x^{2}&\frac{e^{-\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}&\frac{e^{-\frac{5\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&2\Delta x&\frac{17}{4}\Delta
x^{2}&\frac{e^{\frac{\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{7\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\
\end{matrix}\right],$
$B_{0}=\left[\begin{matrix}g^{+}(u_{i-2})&g^{+}(u_{i-1})&g^{+}(u_{i})&g^{+}(u_{i+1})&g^{+}(u_{i+2})\end{matrix}\right]^{T},$
$X_{0}=\left[\begin{matrix}a_{0}^{0}&a_{1}^{0}&a_{2}^{0}&a_{3}^{0}&a_{4}^{0}\end{matrix}\right]^{T}.$

  2. (2)

For $m=1$,

$\displaystyle A_{1}=\left[\begin{matrix}1&-\Delta x&\frac{5}{4}\Delta
x^{2}&\frac{e^{-\frac{5\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\
1&0&\frac{1}{4}\Delta x^{2}&\frac{e^{-\frac{3\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}&\frac{e^{-\frac{3\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&\Delta x&\frac{5}{4}\Delta
x^{2}&\frac{e^{-\frac{\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{5\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&2\Delta
x&\frac{17}{4}\Delta x^{2}&\frac{e^{\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}&\frac{e^{-\frac{7\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&3\Delta x&\frac{37}{4}\Delta
x^{2}&\frac{e^{\frac{3\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{9\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\
\end{matrix}\right],$
$B_{1}=\left[\begin{matrix}g^{+}(u_{i-1})&g^{+}(u_{i})&g^{+}(u_{i+1})&g^{+}(u_{i+2})&g^{+}(u_{i+3})\end{matrix}\right]^{T},$
$X_{1}=\left[\begin{matrix}a_{0}^{1}&a_{1}^{1}&a_{2}^{1}&a_{3}^{1}&a_{4}^{1}\end{matrix}\right]^{T}.$

  3. (3)

For $m=2$,

$\displaystyle A_{2}=\left[\begin{matrix}1&0&\frac{1}{4}\Delta
x^{2}&\frac{e^{-\frac{3\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{3\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&\Delta
x&\frac{5}{4}\Delta x^{2}&\frac{e^{-\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}&\frac{e^{-\frac{5\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&2\Delta x&\frac{17}{4}\Delta
x^{2}&\frac{e^{\frac{\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{7\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&3\Delta
x&\frac{37}{4}\Delta x^{2}&\frac{e^{\frac{3\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta
x^{3}}&\frac{e^{-\frac{9\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}\\\ 1&4\Delta x&\frac{65}{4}\Delta
x^{2}&\frac{e^{\frac{5\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})^{3}}{\lambda^{3}\Delta x^{3}}&\frac{e^{-\frac{11\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})^{3}}{\lambda^{3}\Delta x^{3}}\\\
\end{matrix}\right],$
$B_{2}=\left[\begin{matrix}g^{+}(u_{i})&g^{+}(u_{i+1})&g^{+}(u_{i+2})&g^{+}(u_{i+3})&g^{+}(u_{i+4})\end{matrix}\right]^{T},$
$X_{2}=\left[\begin{matrix}a_{0}^{2}&a_{1}^{2}&a_{2}^{2}&a_{3}^{2}&a_{4}^{2}\end{matrix}\right]^{T}.$

Substituting the obtained coefficients of $X_{0},$ $X_{1}$ and $X_{2}$ into
(3.12) and then calculating $p_{m}(x)=q_{m}(x+\Delta
x)-2q_{m}(x)+q_{m}(x-\Delta x)$ at $x=x_{i+\frac{1}{2}}$, we get

(3.15) $\displaystyle{\hat{G}^{(0)}_{i+\frac{1}{2}}}$
$\displaystyle=C_{0}^{0}g^{+}(u_{i-2})+C_{1}^{0}g^{+}(u_{i-1})+C_{2}^{0}g^{+}(u_{i})+C_{3}^{0}g^{+}(u_{i+1})+C_{4}^{0}g^{+}(u_{i+2}),$
$\displaystyle{\hat{G}^{(1)}_{i+\frac{1}{2}}}$
$\displaystyle=C_{0}^{1}g^{+}(u_{i-1})+C_{1}^{1}g^{+}(u_{i})+C_{2}^{1}g^{+}(u_{i+1})+C_{3}^{1}g^{+}(u_{i+2})+C_{4}^{1}g^{+}(u_{i+3}),$
$\displaystyle{\hat{G}^{(2)}_{i+\frac{1}{2}}}$
$\displaystyle=C_{0}^{2}g^{+}(u_{i})+C_{1}^{2}g^{+}(u_{i+1})+C_{2}^{2}g^{+}(u_{i+2})+C_{3}^{2}g^{+}(u_{i+3})+C_{4}^{2}g^{+}(u_{i+4}).$

After applying the Taylor series to the coefficients $C_{j}^{m},0\leq j\leq
4,0\leq m\leq 2$ of equation (3.15), we get By shifting each index by -1, we
obtain the flux ${\hat{G}^{(m)}_{i-\frac{1}{2}}}$. Hence the Taylor series
expansions of (LABEL:eqn:22) gives

(3.17) $\displaystyle{\hat{G}^{(0)}_{i\pm\frac{1}{2}}}$
$\displaystyle=G(x_{i\pm\frac{1}{2}})+\biggl{(}-\frac{1}{8}\lambda^{2}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}+\frac{1}{8}\frac{d^{5}h}{dx^{5}}\bigg{|}_{x=x_{i}}\biggr{)}\Delta
x^{5}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle{\hat{G}^{(1)}_{i\pm\frac{1}{2}}}$
$\displaystyle=G(x_{i\pm\frac{1}{2}})+\biggl{(}\frac{1}{8}\lambda^{2}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}-\frac{1}{8}\frac{d^{5}h}{dx^{5}}\bigg{|}_{x=x_{i}}\biggr{)}\Delta
x^{5}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle{\hat{G}^{(2)}_{i\pm\frac{1}{2}}}$
$\displaystyle=G(x_{i\pm\frac{1}{2}})+\biggl{(}-\frac{1}{8}\lambda^{2}\frac{d^{3}h}{dx^{3}}\bigg{|}_{x=x_{i}}+\frac{1}{8}\frac{d^{5}h}{dx^{5}}\bigg{|}_{x=x_{i}}\biggr{)}\Delta
x^{5}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle{\hat{G}^{(m)}_{i\pm\frac{1}{2}}}$
$\displaystyle=G(x_{i\pm\frac{1}{2}})+\mathcal{B}_{m}\Delta
x^{5}+\mathcal{O}(\Delta x^{6}),$

where the values $\mathcal{B}_{m}$ are independent of $\Delta x$, but
dependent on $\lambda$. Hence, we have third order approximations

(3.18) $\displaystyle\begin{aligned}
g(u)_{xxx}\big{|}_{x=x_{i}}&=\frac{{\hat{G}^{(0)}_{i+\frac{1}{2}}}-{\hat{G}^{(0)}_{i-\frac{1}{2}}}}{\Delta
x^{3}}+\mathcal{O}(\Delta x^{3}),\\\
g(u)_{xxx}\big{|}_{x=x_{i}}&=\frac{{\hat{G}^{(1)}_{i+\frac{1}{2}}}-{\hat{G}^{(1)}_{i-\frac{1}{2}}}}{\Delta
x^{3}}+\mathcal{O}(\Delta x^{3}),\\\
g(u)_{xxx}\big{|}_{x=x_{i}}&=\frac{{\hat{G}^{(2)}_{i+\frac{1}{2}}}-{\hat{G}^{(2)}_{i-\frac{1}{2}}}}{\Delta
x^{3}}+\mathcal{O}(\Delta x^{3}).\end{aligned}$

###  3.3. Ideal weights based on exponential polynomials

The final WENO-E approximation is defined by a convex combination of local
fluxes with non-linear weights $\omega_{m}$:

(3.19)
${\hat{G}_{i\pm\frac{1}{2}}}=\sum_{m=0}^{2}\omega_{m}^{\pm}{\hat{G}^{(m)}_{i+\frac{1}{2}}},$

To derive the non-linear weights $\omega_{m}$, we initially determine the
$d_{m}$, referred to as ideal (or optimal) weights. These $d_{m}$ values are
chosen such that their linear combination with
${\hat{G}^{(m)}_{i+\frac{1}{2}}}$ retains the fifth convergence order to
${G(x_{i+\frac{1}{2}})}$. That is,

(3.20)
${\hat{G}^{+}_{i+\frac{1}{2}}}=\sum_{m=0}^{2}d_{m}{\hat{G}^{(m)}_{i+\frac{1}{2}}},$

satisfying $\sum_{m=0}^{2}d_{m}=1$. The ideal weights $d_{m}$, for the
proposed WENO-E scheme can be obtained as

(3.21) $d_{0}=\frac{C_{0}}{C_{0}^{0}},\quad
d_{1}=\frac{C_{1}-d_{0}C_{1}^{0}}{C_{0}^{1}},\quad
d_{2}=\frac{C_{2}-d_{0}C_{2}^{0}-d_{1}C_{1}^{1}}{C_{0}^{2}}.$

In contrast to the classical WENO scheme, the optimal weights $d_{m}$ in the
proposed WENO method may exhibit variation depending on the parameter
$\lambda$’s choice, but they converge towards the original ideal weights as
$\Delta x\to 0$. Adding and subtracting
$\sum_{m=0}^{2}d_{m}{\hat{G}^{(m)}_{i\pm\frac{1}{2}}}$ to (3.19), gives:

(3.22)
${\hat{G}_{i\pm\frac{1}{2}}}=\sum_{m=0}^{2}d_{m}{\hat{G}^{(m)}_{i\pm\frac{1}{2}}}+\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m}){\hat{G}^{(m)}_{i\pm\frac{1}{2}}}=\bigl{[}G(x_{i\pm\frac{1}{2}})++\mathcal{A}^{\pm}\Delta
x^{7}+\mathcal{O}(\Delta
x^{8})\bigr{]}+\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m}){\hat{G}^{(m)}_{i\pm\frac{1}{2}}}.$

(The superscripts $\pm$ corresponds to the $\pm$ in the
${\hat{G}^{(m)}_{i\pm\frac{1}{2}}}$). Expanding the second term with the help
of (3.17) we obtain

(3.23)
$\begin{split}\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m}){\hat{G}^{(m)}_{i\pm\frac{1}{2}}}&=\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m})\biggl{(}G(x_{i\pm\frac{1}{2}})+\mathcal{B}_{m}\Delta
x^{5}+\mathcal{O}(\Delta x^{6})\biggr{)}\\\
&=G(x_{i\pm\frac{1}{2}})\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m})+\Delta
x^{5}\sum_{m=0}^{2}\mathcal{B}_{m}(\omega_{m}^{\pm}-d_{m})+\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m})\mathcal{O}(\Delta
x^{6})\end{split}$

Substituting the result above at a finite difference formula for the
exponential polynomial approximation ${\hat{G}_{i\pm\frac{1}{2}}}$:

(3.24)
$\begin{split}\frac{\hat{G}_{i+\frac{1}{2}}-\hat{G}_{i-\frac{1}{2}}}{\Delta
x^{3}}&=\frac{G(x_{i+\frac{1}{2}})-G(x_{i-\frac{1}{2}})}{\Delta
x^{3}}+\mathcal{O}(\Delta
x^{5})+\Biggl{[}\frac{\sum_{m=0}^{2}(\omega_{m}^{+}-d_{m}){\hat{G}^{(m)}_{i+\frac{1}{2}}}-\sum_{m=0}^{2}(\omega_{m}^{-}-d_{m}){\hat{G}^{(m)}_{i-\frac{1}{2}}}}{\Delta
x^{3}}\Biggr{]},\\\ &=g(u)_{xxx}\big{|}_{x=x_{i}}+\mathcal{O}(\Delta
x^{5})+\Biggl{[}\frac{G(x_{i+\frac{1}{2}})\sum_{m=0}^{2}(\omega_{m}^{+}-d_{m})-G(x_{i-\frac{1}{2}})\sum_{m=0}^{2}(\omega_{m}^{-}-d_{m})}{\Delta
x^{3}}\Biggr{]}\\\ &+\Delta
x^{2}\sum_{m=0}^{2}\mathcal{B}_{m}(\omega_{m}^{+}-\omega_{m}^{-})+\Biggl{[}\sum_{m=0}^{2}(\omega_{m}^{+}-d_{m})-\sum_{m=0}^{2}(\omega_{m}^{-}-d_{m})\Biggr{]}\mathcal{O}(\Delta
x^{3}).\end{split}$

The $\mathcal{O}(\Delta x^{5})$ term remains after division by $\Delta x^{3}$
because $\mathcal{A}^{+}=\mathcal{A}^{-}$ in (3.22). Thus, necessary and
sufficient conditions for fifth-order convergence in (2.7) are given by

(3.25)
$\begin{split}\sum_{m=0}^{2}(\omega_{m}^{\pm}-d_{m})&=\mathcal{O}(\Delta
x^{8}),\\\
\sum_{m=0}^{2}\mathcal{B}_{m}(\omega_{m}^{+}-\omega_{m}^{-})&=\mathcal{O}(\Delta
x^{3}),\\\ \omega_{m}^{\pm}-d_{m}&=\mathcal{O}(\Delta x^{2}).\end{split}$

Note that the first constraint is always satisfied due to normalization
$\sum_{m=0}^{2}\omega_{m}^{\pm}=\sum_{m=0}^{2}d_{m}$ and, from (3.23), we see
that a sufficient condition for fifth-order of convergence is given by

(3.26) $\omega_{m}^{\pm}-d_{m}=\mathcal{O}(\Delta x^{3}).$

##  4\. Smoothness indicators and non-linear weights

The smoothness indicator plays a pivotal role in WENO reconstruction, as it
serves as a fundamental factor in determining the non-linear weights. These
weights are derived by assessing the smoothness of the local solution within
each sub-stencil $S_{m}$. In this section, we introduce a novel set of non-
linear weights that enhances existing fifth-order WENO schemes. We constructs
the local and global smoothness indicators employing an $L^{1}$-norm approach
[26]. This global indicator gauges the approximate magnitude of the
derivatives of the local solution within each sub-stencil. Moreover, it is
demonstrated that within smooth regions, the nonlinear weights closely
approximate the linear weights at a rate of $\mathcal{O}(\Delta x^{3})$, even
near the critical points. The proof presented in this section establishes the
fulfilment of condition (3.26), thereby ensuring that the newly proposed
scheme achieves fifth-order accuracy in smooth areas.

###  4.1. Development of global and local smoothness indicators

For the construction of a smoothness indicator, we use $n^{\text{th}}$-order
generalized undivided differences $\mathcal{D}_{m}^{n}g(u)$ ($n=3,4$) of
$g(u)$ on the stencil $S_{m}$, $m=0,1,2$ is given by

(4.1) $\mathcal{D}_{m}^{n}g(u_{i+\frac{1}{2}}):=\sum_{x_{j}\in
S_{m}}a_{m,j}^{[n]}g(u(x_{j})).$

Let $n_{m}$ denote the number of points inside the stencil $S_{m}$ and define
the coefficient vector $\mathbf{a}_{m}^{[n]}:=(a_{m,j}^{[n]}:x_{j}\in
S_{m})^{T}$ in (4.1) by solving the linear system

$\mathbf{V}\cdot\mathbf{a}^{[n]}_{m}=\mathbf{\delta}^{[n]},$

for the non-singular matrix

$\mathbf{V}:=\biggl{(}\frac{(x_{j}-x_{i+\frac{1}{2}})^{l}}{\Delta
x^{l}l!}:x_{j}\in
S_{m},l=0,\dots,n_{m}-1\biggr{)},\quad\text{and}\quad\mathbf{\delta}^{[n]}:=(\delta_{n,l}:l=0,\dots,n_{m}-1)^{T}.$

Note that the coefficients in (4.1) are independent of $\Delta x$ and
evaluation point $x_{i+\frac{1}{2}}$. The operators $\mathcal{D}_{m}^{3}g(u)$
and $\mathcal{D}_{m}^{4}g(u)$ for $m=0,1,2,$ can be written as

(4.2) $\displaystyle\mathcal{D}_{0}^{3}g(u)$
$\displaystyle=-g(u_{i-1})+3g(u_{i})-3g(u_{i+1})+g(u_{i+2}),$
$\displaystyle\mathcal{D}_{1}^{3}g(u)$
$\displaystyle=-g(u_{i-1})+3g(u_{i})-3g(u_{i+1})+g(u_{i+2}),$
$\displaystyle\mathcal{D}_{2}^{3}g(u)$
$\displaystyle=-2g(u_{i})+7g(u_{i+1})-9g(u_{i+2})+5g(u_{i+3})-g(u_{i+4}),$
$\displaystyle\mathcal{D}_{0}^{4}g(u)$
$\displaystyle=g(u_{i-2})-4g(u_{i-1})+6g(u_{i})-4g(u_{i+1})+g(u_{i+2}),$
$\displaystyle\mathcal{D}_{1}^{4}g(u)$
$\displaystyle=g(u_{i-1})-4g(u_{i})+6g(u_{i+1})-4g(u_{i+2})+g(u_{i+3}),$
$\displaystyle\mathcal{D}_{2}^{4}g(u)$
$\displaystyle=g(u_{i})-4g(u_{i+1})+6g(u_{i+2})-4g(u_{i+3})+g(u_{i+4}).$

Then a simple calculation with Taylor expansion shows that they can
approximate the derivatives with higher accuracy than classical undivided
differences.

######  Theorem 4.1.

Let the stencil $S_{m}$ be a stencil around $x_{i+\frac{1}{2}}$ with
$\\#S_{m}=n_{m}$, and assume that $g\in\mathcal{C}^{n_{m}}(\Omega)$, where
$\Omega$ is an open interval containing $S_{m}$. Then, the functional
$\mathcal{D}_{m}^{n}g(u_{i+\frac{1}{2}})$ in (4.1) has the convergence
property

(4.3)
$\mathcal{D}_{m}^{n}g(u_{i+\frac{1}{2}})=\frac{d^{n}}{dx^{n}}g(u)\bigg{|}_{x=x_{i+\frac{1}{2}}}\Delta
x^{n}+\mathcal{O}(\Delta x^{5}),\quad n=3,4.$

We define the smoothness indicator $\beta_{m}$ in each substencil by

(4.4) $\beta_{m}=|\mathcal{D}_{m}^{3}g(u)|+|\mathcal{D}_{m}^{4}g(u)|,\quad
m=0,1,2,$

and the global smoothness indicator $\zeta$ is simply defined as the absolute
difference between $\beta_{0}$ and $\beta_{2}$, i.e.,

(4.5) $\zeta=|\beta_{0}-\beta_{2}|.$

###  4.2. Construction of non-linear weights and analysis of convergence order

The non-linear weights for the scheme are defined as [32]

(4.6)
$\omega_{m}=\frac{\alpha_{m}}{\sum_{j=0}^{2}\alpha_{j}},\quad\alpha_{m}=d_{m}\biggl{(}1+\frac{\zeta}{\beta_{m}+\Delta
x^{2}}\biggr{)},m=0,1,2.$

To attain the optimal order approximation for $G(x_{i+\frac{1}{2}})$ in smooth
regions, the weights should converge appropriately towards the ideal weights
as $\Delta x$ approaches zero. Conversely, in regions where a discontinuity is
present, the weights should effectively eliminate the contribution of stencils
that contain the discontinuity. Using the concept of Taylor expansion, the
operators $\mathcal{D}_{m}^{3}g(u)$ and $\mathcal{D}_{m}^{4}g(u)$ for
$m=0,1,2$ can be represented as follows

(4.7) $\displaystyle\mathcal{D}_{0}^{3}g(u)$ $\displaystyle=\Delta
x^{3}g^{(3)}_{i+\frac{1}{2}}+\frac{1}{8}\Delta
x^{5}g^{(5)}_{i+\frac{1}{2}}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle\mathcal{D}_{1}^{3}g(u)$ $\displaystyle=\Delta
x^{3}g^{(3)}_{i+\frac{1}{2}}+\frac{1}{8}\Delta
x^{5}g^{(5)}_{i+\frac{1}{2}}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle\mathcal{D}_{2}^{3}g(u)$ $\displaystyle=\Delta
x^{3}g^{(3)}_{i+\frac{1}{2}}-\frac{7}{8}\Delta
x^{5}g^{(5)}_{i+\frac{1}{2}}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle\mathcal{D}_{0}^{4}g(u)$ $\displaystyle=\Delta
x^{4}g^{(4)}_{i+\frac{1}{2}}-\frac{1}{2}\Delta
x^{5}g^{(5)}_{i+\frac{1}{2}}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle\mathcal{D}_{1}^{4}g(u)$ $\displaystyle=\Delta
x^{4}g^{(4)}_{i+\frac{1}{2}}+\frac{1}{2}\Delta
x^{5}g^{(5)}_{i+\frac{1}{2}}+\mathcal{O}(\Delta x^{6}),$
$\displaystyle\mathcal{D}_{2}^{4}g(u)$ $\displaystyle=\Delta
x^{4}g^{(4)}_{i+\frac{1}{2}}+\frac{3}{2}\Delta
x^{5}g^{(5)}_{i+\frac{1}{2}}+\mathcal{O}(\Delta x^{6}).$

By definition of each $\beta_{m}$ with $m=0,1,2,$ it is straightforward that
the truncation error of the local smooth indicators $\beta_{m}$ is of the form

(4.8) $\displaystyle\beta_{0}$ $\displaystyle=|g^{(3)}_{i+\frac{1}{2}}|\Delta
x^{3}+|g^{(4)}_{i+\frac{1}{2}}|\Delta
x^{4}-\frac{3}{8}|g^{(5)}_{i+\frac{1}{2}}|\Delta
x^{5}+\frac{7}{24}|g^{(6)}_{i+\frac{1}{2}}|\Delta
x^{6}-\frac{187}{1920}|g^{(7)}_{i+\frac{1}{2}}|\Delta x^{7}+\mathcal{O}(\Delta
x^{8}),$ $\displaystyle\beta_{1}$
$\displaystyle=|g^{(3)}_{i+\frac{1}{2}}|\Delta
x^{3}+|g^{(4)}_{i+\frac{1}{2}}|\Delta
x^{4}+\frac{5}{8}|g^{(5)}_{i+\frac{1}{2}}|\Delta
x^{5}+\frac{7}{24}|g^{(6)}_{i+\frac{1}{2}}|\Delta
x^{6}+\frac{71}{640}|g^{(7)}_{i+\frac{1}{2}}|\Delta x^{7}+\mathcal{O}(\Delta
x^{8}),$ $\displaystyle\beta_{2}$
$\displaystyle=|g^{(3)}_{i+\frac{1}{2}}|\Delta
x^{3}+|g^{(4)}_{i+\frac{1}{2}}|\Delta
x^{4}+\frac{5}{8}|g^{(5)}_{i+\frac{1}{2}}|\Delta
x^{5}+\frac{7}{24}|g^{(6)}_{i+\frac{1}{2}}|\Delta
x^{6}+\frac{71}{640}|g^{(7)}_{i+\frac{1}{2}}|\Delta x^{7}+\mathcal{O}(\Delta
x^{8}).$ $\displaystyle\zeta$ $\displaystyle=|g^{(5)}_{i+\frac{1}{2}}|\Delta
x^{5}+\frac{5}{24}|g^{(7)}_{i+\frac{1}{2}}|\Delta x^{7}+\mathcal{O}(\Delta
x^{9}).$

It implies that

(4.9) $\alpha_{m}=d_{m}\biggl{(}1+\frac{\zeta}{\beta_{m}+\Delta
x^{2}}\biggr{)}=\begin{cases}d_{m}\bigl{(}1+\mathcal{O}(\Delta
x^{3})\bigr{)},&\text{if}\quad g^{+}_{xxx}\neq 0,\\\
d_{m}\bigl{(}1+\mathcal{O}(\Delta x^{3})\bigr{)},&\text{if}\quad
g^{+}_{xxx}=0,g^{+}_{xxxx}\neq 0.\end{cases}.$

Then, using the fact that $\sum_{m=0}^{2}d_{m}=1$ and
$\omega_{m}=\frac{\alpha_{m}}{\sum_{j=0}^{2}\alpha_{j}}$

(4.10) $d_{m}=\frac{1}{1+\frac{\zeta}{\beta_{m}+\Delta
x^{2}}}\omega_{m}\bigl{(}1+\mathcal{O}(\Delta
x^{3})\bigr{)}\sum_{m=0}^{2}d_{m}=\omega_{m}+\mathcal{O}(\Delta x^{3}).$

Thus, the non-linear weights satisfy the sufficient condition (3.26).

######  Theorem 4.2.

Assume that $d_{m},m=0,1,2$, is the linear weights and that the function
$g(u)_{xxx}$ is smooth around the big stencil $S_{7}$, so the nonlinear
weights (4.6) satisfy the conditions (3.26) even in the presence of critical
points where the third and higher derivatives vanish.

Now we summarise the proposed scheme in algorithmic fashion. We iteratively
solve equation (2.1) until the final time $t=T$. The numerical solution at the
$n^{\text{th}}$ time step $t=t^{n}$ is denoted by
$\\{u^{n}_{i}:=u^{n}(x_{i})\\}$. Commencing with $n=0$ and $\\{u^{0}_{i}\\}$
being a given initial condition, the following steps are taken:

Algorithm 1

  1. (1)

Establishing a uniform mesh distribution.

    1. (a)

Define the spatial domain as $x\in[x_{l},x_{r}]$.

    2. (b)

Set the number of grid points as $N$.

    3. (c)

Calculate the spatial step size as $\Delta x=\frac{(x_{r}-x_{l})}{N}$.

  2. (2)

Finite difference WENO scheme for the third derivatives.

    1. (a)

The conservative form of the equation within the MOL framework is defined,
employing the numerical fluxes $\hat{F}$ and $\hat{G}$ for convection and
dispersion, correspondingly.

    2. (b)

For the convection term, the flux construction procedure outlined in reference
[11] is followed.

    3. (c)

The approximation of fluxes for the dispersion term is achieved by applying
the WENO reconstruction approach as detailed in equation (2.7).

  3. (3)

Construct exponential WENO-E reconstruction at the cell interface
$x_{i+\frac{1}{2}}$.

    1. (a)

Split the numerical flux $\hat{G}_{i+\frac{1}{2}}$ into positive part
$\hat{G}^{+}_{i+\frac{1}{2}}$ and negative part $\hat{G}^{-}_{i+\frac{1}{2}}$
as in (2.8).

    2. (b)

Construct the approximations $\hat{G}^{+}_{i+\frac{1}{2}}$ based on
exponential polynomial (2.12) on stencil $S_{7}$ and local approximations
$\hat{G}^{(m)}_{i+\frac{1}{2}}$ based on exponential polynomial (2.13) on each
substencil $S_{m}$ for $m=0,1,2.$

    3. (c)

Determine the ideal weights $d_{m}$ and form the approximation given by
(3.20).

    4. (d)

Analogously, replication of steps 3(b) and 3(c) to obtain
$\hat{G}^{-}_{i+\frac{1}{2}}$.

  4. (4)

Compute non-linear weights.

    1. (a)

Calculate undivided differences of $g(u)$ on the associated substencils :
$\mathcal{D}^{n}_{S^{+}_{m}}$, $\mathcal{D}^{n}_{S^{-}_{m}}$ for $n=3,4$,
$m=0,1,2.$

    2. (b)

Evaluate the smoothness indicators $\beta_{m}$ on each substencil $S_{m}$ as
in (4.4) and determine the value $\zeta$ as in (4.5).

    3. (c)

Compute the unnormalize weights $\alpha_{m}$ for each $m=0,1,2$ and using
$\zeta$ and $\alpha_{m}$, compute the nonlinear weights $\omega_{m}$ as in
(4.6).

  5. (5)

Repeat step 3 and 4, to construct numerical flux $\hat{G}_{i-\frac{1}{2}}$ by
shifting each index by -1.

  6. (6)

Time discretization.

    1. (a)

Update the time step from $t^{n}$ to $t^{n+1}=t^{n}+\Delta t^{n}$, by applying
third order SSP Runge-Kutta scheme (5.2).

    2. (b)

If $t^{n+1}<T$, set the time step size $\Delta t^{n+1}$ by (5.3) with updated
wave speeds. In case $t^{n+1}+\Delta t^{n+1}>T$, $\Delta t^{n+1}$ is set
$T-t^{n+1}$.

    3. (c)

Repeat steps beginning with step 3 for each time step until time $T$.

######  Remark 4.3.

When $\lambda\Delta x$ is close to zero, the WENO-E scheme reduces to the WENO
scheme based on algebraic polynomials as given in [9]. In the polynomial case,
the numerical flux on the big stencil $S_{7}$ is given by

(4.11)
${\hat{G}^{+}_{i+\frac{1}{2}}}=\biggl{[}\frac{-1}{15}g^{+}(u_{i-2})+\frac{21}{40}g^{+}(u_{i-1})+\frac{1}{8}g^{+}(u_{i})-\frac{23}{12}g^{+}(u_{i+1})+\frac{7}{4}g^{+}(u_{i+2})-\frac{19}{40}g^{+}(u_{i+3})+\frac{7}{120}g^{+}(u_{i+4})\biggr{]}.$

Also the numerical fluxes on the small stencils
$S^{m}=\\{x_{i-2+m},\dots,x_{i+2+m}\\}$ with $m=0,1,2$ are given by

(4.12) $\displaystyle{\hat{G}^{(0)}_{i+\frac{1}{2}}}$
$\displaystyle=-\frac{1}{4}g^{+}(u_{i-2})+\frac{3}{2}g^{+}(u_{i-1})-2g^{+}(u_{i})+\frac{1}{2}g^{+}(u_{i+1})+\frac{1}{4}g^{+}(u_{i+2}),$
$\displaystyle{\hat{G}^{(1)}_{i+\frac{1}{2}}}$
$\displaystyle=\frac{1}{4}g^{+}(u_{i-1})+\frac{1}{2}g^{+}(u_{i})-2g^{+}(u_{i+1})+\frac{3}{2}g^{+}(u_{i+2})-\frac{1}{4}g^{+}(u_{i+3}),$
$\displaystyle{\hat{G}^{(2)}_{i+\frac{1}{2}}}$
$\displaystyle=\frac{7}{4}g^{+}(u_{i})-\frac{9}{2}g^{+}(u_{i+1})+4g^{+}(u_{i+2})-\frac{3}{2}g^{+}(u_{i+3})+\frac{1}{4}g^{+}(u_{i+4}).$

The smoothness indicator $\beta_{m}$ of the interpolation polynomials on the
smaller stencils $S_{m}$ are given by

(4.13) $\beta_{m}^{Z}=\sum_{\kappa=1}^{r}\Delta
x^{2\kappa-1}\int_{I_{j}}\biggl{(}\frac{d^{\kappa}}{dx^{\kappa}}\tilde{p}_{m}(x)\biggr{)}^{2}dx.$

In order to increase the accuracy of the nonlinear weights, an alternative to
mapped nonlinear weights [9], we chose to use the WENO-Z nonlinear weights
introduced in [17]. The nonlinear weights $\omega_{m}^{Z}$ of WENO-Z are
defined by

(4.14)
$\omega_{m}^{Z}=\frac{\alpha_{m}^{Z}}{\sum_{k=0}^{2}\alpha_{k}^{Z}},\quad\alpha_{m}^{Z}=d_{m}\biggl{(}1+\frac{\tau_{5}}{\beta_{m}^{Z}+\epsilon}\biggr{)},\quad\tau_{5}:=|\beta_{0}^{Z}-\beta_{2}^{Z}|,\quad
m=0,1,2,$

where $d_{0}=\frac{4}{15}$, $d_{1}=\frac{1}{2}$, $d_{2}=\frac{7}{30}$ are
linear weights, $\alpha_{m}^{Z}$ stands for unnormalized weights.

######  Remark 4.4.

For the convection term, we use the fifth-order finite difference WENO scheme
with WENO-Z nonlinear weights as (4.14) and its corresponding linear weights
in [17].

######  Remark 4.5.

In later numerical experiments, we use $\epsilon=\Delta x^{2}$ to ensure the
fifth-order convergence for WENO-Z and WENO-E schemes, even at critical
points, while keeping ENO property near discontinuity [32, 33].

######  Remark 4.6.

The extension of the present scheme to the two-dimensional case can be done by
discretizing the spatial dimensions one at a time.

##  5\. Numerical results

In this section, we illustrate several examples to test the proposed
WENO-E-$\lambda\Delta x$ scheme (henceforth referred to as the WENO-E scheme)
derived by the exponential approximation space for linear and non-linear
dispersion-type equations with various initial conditions. We compared the
performance of the proposed scheme with the WENO-Z scheme, which is based on
polynomial approximations. In the WENO-E scheme, the tension parameter
$\lambda$ plays a significant role. The numerical value for $\lambda$ can be
tuned according to the characteristics of the initial data; which however is
time-consuming and resource intense. To simplify, the $\lambda$ value is
chosen based on the initial cell size such that $0\leq\lambda\Delta x\leq
0.1$. When $\lambda\Delta x$ is close to zero, the WENO-E scheme reduces to
the standard WENO-Z scheme based on algebraic polynomials. On the other hand,
when $\lambda\Delta x$ is large, the exponential polynomials can better
capture the rapid variations in the solution near discontinuities, resulting
in improved interpolation results. However, choosing a very large value of
$\lambda\Delta x$ can lead to numerical instability and accuracy issues. Thus,
an optimal value of $\lambda$ is chosen depending on the given data feature
and the required accuracy. For the time integration, we use the third-order
explicit strong stability preserving (SSP) Runge-Kutta method [13]. For the
ordinary differential equation of the form

(5.1) $\frac{du}{dt}=RHS(u),$

the SSP Runge-Kutta method is given by

(5.2) $\begin{split}u^{(1)}&=u^{n}+\Delta tRHS(u^{n}),\\\
u^{(2)}&=\frac{3}{4}u^{n}+\frac{1}{4}u^{(1)}+\frac{1}{4}\Delta
tRHS(u^{(1)}),\\\
u^{n+1}&=\frac{1}{3}u^{n}+\frac{2}{3}u^{(2)}+\frac{2}{3}\Delta
tRHS(u^{(2)}).\\\ \end{split}$

To ensure the $CFL$ stability condition, the time step is given by

(5.3) $\Delta t\leq\min\Biggl{(}\frac{CFL\cdot\Delta
x^{\frac{5}{3}}}{\max(|f^{\prime}(u)|)},\frac{CFL\cdot\Delta
x^{3}}{\max(|g^{\prime}(u)|)}\Biggr{)}.$

We have chosen $CFL=0.3$ as obtained by Fourier analysis for the linear
dispersion equation [9] for further analysis.

######  Example 5.1.

Consider the linear KdV equations (Airy equation),

(5.4) $\left\\{\begin{aligned}
u_{t}+u_{xxx}&=0,\quad(x,t)\in[0,2\pi]\times[0,T],\\\ u_{0}(x)&=\sin(x),\quad
x\in[0,2\pi],\end{aligned}\right.$

with periodic boundary conditions and the exact solution is given by
$u(x,t)=\sin(x+t)$. We solve this linear equation up to the final time $T=1$
by the WENO-E scheme and the WENO-Z scheme. The following error norms are used
to compute the accuracy of schemes.

$L^{\infty}=\max_{0\leq i\leq
N}|u_{e}-u_{a}|,\,\,\,\,\,\,\,\,\,\,L^{1}=\dfrac{1}{N+1}\sum_{i=0}^{N}|u_{e}-u_{a}|,$

where $u_{e}$ and $u_{a}$ denote the exact and approximate solutions of the
PDE. A comprehensive analysis of error norms, specifically $L^{1}$ and
$L^{\infty}$ errors, is conducted across various combinations of the parameter
$\lambda\Delta x$ and grid resolution $N$. The results are detailed in
Appendix 1, where we have considered 15 distinct $\lambda\Delta x$ values
carefully chosen to span the entire error spectrum within the $[0,1]$ range.
$\lambda\Delta x$ can be categorized into three classes based on the observed
outcomes. Firstly, there exists a category where $\lambda\Delta x$
consistently yields precise and accurate results. Secondly, values in this
category demonstrate reasonably accurate results for lower resolution and
remains stagnant over larger values of $N$. Lastly, there is a subset of
$\lambda\Delta x$ values significantly deviating from the optimal
$\lambda\Delta x$, resulting in inconsistent outcomes.

In this example, we observe proximate convergence within the range of
$0.02\leq\lambda\Delta x\leq 0.04$, with $\lambda\Delta x=0.02$ exhibiting the
lowest errors among these values. Additionally, we deliberately included
$\lambda\Delta x$ values of $0.06$ and $0.1$ in our analysis to illustrate the
substantial incongruity in results that arises when $\lambda\Delta x$ values
significantly diverge from the optimal choice. Figure 2 provides a comparison
of errors between the WENO-Z scheme and WENO-E scheme with $\lambda\Delta x$
values of $0.02$, $0.04$, $0.06$, and $0.1$ in terms of the $L^{\infty}$\- and
$L^{1}$-norms, respectively. Furthermore, Table 1 presents a summary of the
$L^{\infty}$-error, $L^{1}$-error, and order of convergence for all schemes
under consideration.

Table 1. Comparison of WENO-Z and WENO-E schemes in terms of $L^{\infty}$\-
and $L^{1}$\- errors along with their convergence rate for Example 5.1 over
the domain $\Omega=[0,2\pi]$ at time $T=1$.
$\bm{N}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
---|---|---|---|---|---
|  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate
10 | 2.6042e-03 | - | 2.5610e-03 | - | 2.5611e-03 | - | 2.5611e-03 | - | 2.5602e-03 | -
20 | 8.704e-05 | 4.9029 | 8.7186e-05 | 4.8765 | 8.7170e-05 | 4.8768 | 8.7085e-05 | 4.8782 | 8.6335e-05 | 4.8902
40 | 2.7752e-06 | 4.9711 | 2.7735e-06 | 4.9743 | 2.7627e-06 | 4.9797 | 2.7154e-06 | 5.0032 | 2.3194e-06 | 5.2181
80 | 8.7052e-08 | 4.9946 | 8.6670e-08 | 5.0001 | 8.1182e-08 | 5.0888 | 5.7395e-08 | 5.5641 | 1.4167e-07 | 4.0331
160 | 2.7262e-09 | 4.9969 | 2.5429e-09 | 5.0910 | 2.0500e-10 | 8.6294 | 1.2112e-08 | 2.2445 | 1.1176e-07 | 3.4211e-01
320 | 1.1102e-10 | 4.6180 | 2.1082e-11 | 6.9143 | 1.3553e-09 | -2.7250 | 7.3103e-09 | 7.2848e-01 | 5.7148e-08 | 9.6769e-01
$\bm{N}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
|  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate
10 | 1.7456e-03 | - | 1.7519e-03 | - | 1.7520e-03 | - | 1.7520e-03 | - | 1.7516e-03 | -
20 | 5.6856e-05 | 4.9403 | 5.7100e-05 | 4.9393 | 5.7089e-05 | 4.9396 | 5.7031e-05 | 4.9411 | 5.6536e-05 | 4.9534
40 | 1.7815e-06 | 4.9962 | 1.7820e-06 | 5.0019 | 1.7751e-06 | 5.0072 | 1.7447e-06 | 5.0307 | 1.4902e-06 | 5.2456
80 | 5.5652e-08 | 5.0005 | 5.5421e-08 | 5.0070 | 5.1911e-08 | 5.0957 | 3.6701e-08 | 5.5710 | 9.0592e-08 | 4.0400
160 | 1.7390e-09 | 5.0001 | 1.6221e-09 | 5.0945 | 1.3076e-10 | 8.6330 | 7.7265e-09 | 2.2479 | 7.1295e-08 | 3.4559e-01
320 | 7.0742e-11 | 4.6195 | 1.3412e-11 | 6.9182 | 8.6372e-10 | -2.7236 | 4.6589e-09 | 7.2983e-01 | 3.6420e-08 | 9.6905e-01

(a) (a) Convergence plot: $L^{\infty}$\- error

(b) (b) Convergence plot: $L^{1}$\- error

Figure 2. Comparison WENO-Z and WENO-E schemes in terms of $L^{\infty}$ and
$L^{1}$ errors (in $\log_{10}$ scale) for Example 5.1 at $T=1$

.

######  Example 5.2.

We compute the solution of the two-dimensional linear dispersion equation

(5.5)
$u_{t}+u_{xxx}+u_{yyy}=0,\quad(x,y,t)\in(0,2\pi)\times(0,2\pi)\times(0,T],$

subject to initial data $u(x,y,0)=\sin(x+y)$ with periodic boundary conditions
in both directions. The exact solution is given by $u(x,y,t)=\sin(x+y+2t)$.
The time step is set to a fixed value of $0.1\Delta x^{3}$, and numerical
solutions are computed at a specific time $T=1$. In Figure 3 and Table 2,
shows a comparison of the errors of the WENO-Z scheme and WENO-E scheme with
$\lambda\Delta x=0.02,0.04,0.06,0.1$ in the $L^{\infty}$\- and $L^{1}$-norms,
respectively. As we refined the grid, WENO-E scheme with $\lambda\Delta
x=0.04$ gives more accurate results and converge to the exact solution with
the order of five. Figure 4 shows the numerical solutions on a grid with
$80\times 80$ grid points for WENO-Z and WENO-E scheme with $\lambda\Delta
x=0.04$.

Table 2. Comparison of WENO-Z and WENO-E schemes in terms of $L^{\infty}$\-
and $L^{1}$\- errors along with their convergence rate for Example 5.2 over
the domain $\Omega=(0,2\pi)\times(0,2\pi)$ at time $T=1$.
$\bm{N_{x}\times N_{y}}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
---|---|---|---|---|---
|  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate
10 $\times$ 10 | 5.4139e-03 | - | 5.3032e-03 | - | 5.3035e-03 | - | 5.3035e-03 | - | 5.3020e-03 | -
20 $\times$ 20 | 1.7482e-04 | 4.9527 | 1.7505e-04 | 4.9210 | 1.7502e-04 | 4.9214 | 1.7484e-04 | 4.9228 | 1.7333e-04 | 4.9349
40 $\times$ 40 | 5.5487e-06 | 4.9776 | 5.5462e-06 | 4.9801 | 5.5245e-06 | 4.9855 | 5.4300e-06 | 5.0090 | 4.6380e-06 | 5.2239
80 $\times$ 80 | 1.7411e-07 | 4.9940 | 1.7336e-07 | 4.9997 | 1.6238e-07 | 5.0884 | 1.1480e-07 | 5.5637 | 2.8337e-07 | 4.0327
160 $\times$ 160 | 5.4551e-09 | 4.9963 | 5.0885e-09 | 5.0904 | 4.0691e-10 | 8.6405 | 2.4219e-08 | 2.2450 | 2.2350e-07 | 0.3425
$\bm{N_{x}\times N_{y}}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
|  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate
10 $\times$ 10 | 3.4560e-03 | - | 3.4272e-03 | - | 3.4273e-03 | - | 3.4273e-03 | - | 3.4261e-03 | -
20 $\times$ 20 | 1.1104e-04 | 4.9599 | 1.1149e-04 | 4.9420 | 1.1147e-04 | 4.9423 | 1.1136e-04 | 4.9438 | 1.1039e-04 | 4.9559
40 $\times$ 40 | 3.5318e-06 | 4.9746 | 3.5329e-06 | 4.9800 | 3.5191e-06 | 4.9853 | 3.4588e-06 | 5.0088 | 2.9543e-06 | 5.2237
80 $\times$ 80 | 1.1080e-07 | 4.9944 | 1.1033e-07 | 5.0009 | 1.0335e-07 | 5.0896 | 7.3069e-08 | 5.5649 | 1.8036e-07 | 4.0339
160 $\times$ 160 | 3.4732e-09 | 4.9955 | 3.2400e-09 | 5.0898 | 2.5903e-10 | 8.6402 | 1.5421e-08 | 2.2443 | 1.4231e-07 | 0.3418

(a) (a) Convergence plot: $L^{\infty}$\- error

(b) (b) Convergence plot: $L^{1}$\- error

Figure 3. Comparison WENO-Z and WENO-E schemes in terms of $L^{1}$ and
$L^{\infty}$ errors (in $\log_{10}$ scale) for Example 5.2 at $T=1$.

(a) (a) Numerical solution: WENO-Z

(b) (b) Numerical solution: WENO-E ($\lambda\Delta x=0.04$)

Figure 4.  Comparison of numerical solutions obtained from using WENO-Z and
WENO-E schemes with $N_{x}\times N_{y}=80\times 80$ at $T=1$ in
$(0,2\pi)\times(0,2\pi)$, for Example, 5.2.

######  Example 5.3.

Consider the nonlinear KdV equation

(5.20) $\left\\{\begin{aligned} u_{t}-3(u^{2})_{x}+u_{xxx}&=0,\quad
x\in[-10,10],\quad t\geq 0,\\\ u_{0}(x)&=-2\operatorname{sech}^{2}(x),\quad
x\in[-10,10].\end{aligned}\right.$

The exact solution is given by $u(x,t)=-2\operatorname{sech}^{2}(x-4t)$. We
embark on the numerical computation of the solution within the spatial domain
$x\in[-10,10]$ at a specific time $T=0.5$ while employing periodic boundary
conditions. The $L^{\infty}$-error, $L^{1}$-error and convergence order are
reported in Table 3 and Figure 5. A comparison of numerical solutions obtained
from WENO-Z and WENO-E schemes with $N=640$ at $T=0.5$ in $x\in[-10,10]$ is
given in Figure 6. The result of WENO-E with $\lambda\Delta x=0.04$ is close
to the exact solution and better compared to the polynomial case. As the
number of grid points increases, the error converges to machine epsilon.

Table 3. Comparison of WENO-Z and WENO-E schemes in terms of $L^{\infty}$\-
and $L^{1}$\- errors along with their convergence rate for Example 5.3 over
the domain $\Omega=[-10,10]$ at time $T=0.5$.
$\bm{N}$ | WENO-Z | WENO-E-0.02 | WENO-E- 0.04 | WENO-E-0.06 | WENO-E-0.1
---|---|---|---|---|---
|  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate
80 | 8.3944e-02 | - | 1.2351e-02 | - | 1.2351e-02 | - | 1.2352e-02 | - | 1.2354e-02 | -
160 | 3.6399e-04 | 7.8494 | 4.1591e-04 | 4.8922 | 4.1591e-04 | 4.8922 | 4.1583e-04 | 4.8926 | 4.1486e-04 | 4.8962
320 | 1.2547e-05 | 4.8586 | 1.3086e-05 | 4.9901 | 1.3072e-05 | 4.9918 | 1.3004e-05 | 4.9990 | 1.2429e-05 | 5.0608
640 | 8.5693e-07 | 3.8720 | 8.5676e-07 | 3.9330 | 8.5721e-07 | 3.9306 | 8.5918e-07 | 3.9199 | 8.7566e-07 | 3.8272
1280 | 8.9345e-07 | -0.0602 | 8.9347e-07 | -0.0605 | 8.9369e-07 | -0.0601 | 8.9467e-07 | -0.0584 | 9.0287e-07 | -0.0441
$\bm{N}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
|  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate
80 | 2.3854e-02 | - | 2.1970e-03 | - | 2.1972e-03 | - | 2.1974e-03 | - | 2.1979e-03 | -
160 | 5.9718e-05 | 8.6419 | 6.8676e-05 | 4.9996 | 6.8676e-05 | 4.9997 | 6.8663e-05 | 5.0001 | 6.8506e-05 | 5.0038
320 | 1.9878e-06 | 4.9089 | 2.0648e-06 | 5.0557 | 2.0624e-06 | 5.0574 | 2.0510e-06 | 5.0651 | 1.9547e-06 | 5.1312
640 | 8.4259e-08 | 4.5602 | 8.4750e-08 | 4.6067 | 8.3420e-08 | 4.6278 | 7.7660e-08 | 4.7230 | 4.4961e-08 | 5.4421
1280 | 2.5115e-08 | 1.7463 | 2.5080e-08 | 1.7567 | 2.4566e-08 | 1.7637 | 2.5510e-08 | 1.6061 | 4.8688e-08 | -0.1149

(a) (a) Convergence plot: $L^{\infty}$\- error

(b) (b) Convergence plot: $L^{1}$\- error

Figure 5. Comparison WENO-Z and WENO-E schemes in terms of $L^{1}$ and
$L^{\infty}$ errors (in $\log_{10}$ scale) for Example 5.3 at $T=0.5$.

(a) Numerical solution

(b) Zoom in near smooth extrema

Figure 6.  Comparison of numerical solutions obtained from using WENO-Z and
WENO-E schemes with $N=640$ at $T=0.5$ in $x\in[-10,10]$, for Example 5.3.

Table 4. Comparison of WENO-Z and WENO-E schemes in terms of $L^{\infty}$\-
and $L^{1}$\- errors along with their convergence rate for Example 5.5 over
the domain $\Omega=[0,2\pi]$ at time $T=\pi/2$.
$\bm{N}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
---|---|---|---|---|---
|  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate |  $\bm{L^{\infty}}$-error | Rate
40 | 6.3724e-04 | - | 6.3724e-04 | - | 6.3724e-04 | - | 6.3718e-04 | - | 6.3673e-04 | -
80 | 1.7932e-05 | 5.1513 | 1.7927e-05 | 5.1516 | 1.7920e-05 | 5.1522 | 1.7890e-05 | 5.1544 | 1.7642e-05 | 5.1736
160 | 7.0911e-07 | 4.6604 | 7.0899e-07 | 4.6602 | 7.0660e-07 | 4.6646 | 6.9624e-07 | 4.6835 | 6.6316e-07 | 4.7336
320 | 4.6784e-08 | 3.9219 | 4.6796e-08 | 3.9213 | 4.6983e-08 | 3.9107 | 4.7792e-08 | 3.8647 | 6.4315e-08 | 3.3661
640 | 5.0636e-09 | 3.2078 | 5.0663e-09 | 3.2074 | 5.1073e-09 | 3.2015 | 5.3168e-09 | 3.1682 | 3.4122e-08 | 0.9144
1280 | 4.1144e-10 | 3.6214 | 4.0924e-10 | 3.6299 | 5.0185e-10 | 3.3472 | 2.3221e-09 | 1.1951 | 1.7924e-08 | 0.9288
$\bm{N}$ | WENO-Z | WENO-E-0.02 | WENO-E-0.04 | WENO-E-0.06 | WENO-E-0.1
|  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate |  $\bm{L^{1}}$-error | Rate
40 | 3.1608e-04 | - | 3.1654e-04 | - | 3.1654e-04 | - | 3.1644e-04 | - | 3.1628e-04 | -
80 | 1.1146e-05 | 4.8256 | 1.1149e-05 | 4.8273 | 1.1146e-05 | 4.8278 | 1.1091e-05 | 4.8345 | 1.1005e-05 | 4.8449
160 | 4.5757e-07 | 4.6064 | 4.5746e-07 | 4.6072 | 4.5581e-07 | 4.6120 | 4.2927e-07 | 4.6914 | 3.8846e-07 | 4.8243
320 | 2.0152e-08 | 4.5050 | 2.0104e-08 | 4.5081 | 1.9391e-08 | 4.5550 | 1.3932e-08 | 4.9454 | 2.5432e-08 | 3.9330
640 | 1.4181e-09 | 3.8290 | 1.4092e-09 | 3.8346 | 1.3336e-09 | 3.8620 | 2.1872e-09 | 2.9288 | 1.7251e-08 | 0.5600
1280 | 1.5958e-10 | 3.1516 | 1.5226e-10 | 3.2102 | 2.1720e-10 | 2.6182 | 1.0973e-09 | 0.9951 | 8.7921e-09 | 0.9724

######  Example 5.4.

Consider the following classical KdV equation with zero dispersion limit of
conservation law

(5.51) $\begin{split}u_{t}+\left(\frac{u^{2}}{2}\right)_{x}+\epsilon
u_{xxx}&=0,\quad x\in[0,1],\quad t\geq 0,\\\ \end{split}$

with continuous initial condition

(5.52) $u(x,0)=2+0.5\sin(2\pi x),\quad x\in[0,1],$

and discontinuous initial condition

(5.53) $u(x,0)=\begin{cases}1,&\text{if}\quad 0.25<x<4,\\\
0,&\text{else}.\end{cases}$

(a) $T=0.5$, $N=200$, $\epsilon=10^{-4}$

(b) $T=0.5$, $N=500$, $\epsilon=10^{-5}$

(c) $T=0.5$, $N=1500$, $\epsilon=10^{-6}$

Figure 7. Numerical solution of KdV equation with $N=200,500,1000$ at $T=0.5$
and $\lambda\Delta x=0.02$ in $x\in[0,1]$, WENO-E for initial condition (5.52)
of Example 5.4.

(a) $T=0.01$, $N=1500$, $\epsilon=10^{-4}$

(b) $T=0.05$, $N=1500$, $\epsilon=10^{-4}$

(c) $T=0.1$, $N=1500$, $\epsilon=10^{-4}$

(d) Wave simulation across various T

(e) Wave simulation in 2D

(f) Contour for wave simulation

Figure 8. Numerical solution of KdV equation with $\epsilon=10^{-4}$, $N=1500$
at $T=0.01,0.05,0.1$ and $\lambda\Delta x=0.02$ in $x\in[0,5]$, WENO-E for
initial condition (5.53) of Example 5.4

The dispersive non-linear KdV equation with initial condition (5.52) with
periodic boundary condition has a dispersive shock wave behavior and produces
continuous wavelets in the vicinity of the discontinuity for small $\epsilon$.
This example demonstrates the effectiveness of our scheme in accurately
resolving high-frequency wavelets even for very small values of $\epsilon$.
The solution is computed with $\epsilon=10^{-4},10^{-5},10^{-6}$, respectively
at $T=0.5$ and plotted in Figure 7. The findings indicate that the simulations
accurately capture physical oscillations, and the solutions are devoid of
noise, especially before and after the occurrence of dispersive shocks. Stable
numerical oscillations can be observed in both downstream and upstream of the
continuous wavelets in solutions obtained with significantly coarser grids.

Upon considering the aforementioned equation with a discontinuous initial
condition and inflow-outflow boundary condition, we can see a dispersion-shock
wave propagating to the left at each discontinuous interface. The amplitude of
the wave increases over time. Figure 8 demonstrates how a top-hat initial
solution transforms into a series of traveling waves. Even with physically
discontinuous initial data, the zero-dispersion limit solutions do not exhibit
discontinuity. Instead, they develop continuous fine-scale wavelets and
eventually break up into solitary waves. The dispersion-shock wave plotted
with the current scheme is more accurate and flexible than those plotted in
previous studies [14, 15, 9] with the same mesh.

######  Example 5.5.

The canonical traveling wave solution for the $K(2,2)$ equation

(5.54) $u_{t}+(u^{2})_{x}+(u^{2})_{xxx}=0,$

is given by the compacton

$u(x,t)=\begin{cases}\frac{4c}{3}\cos^{2}\left(\frac{x-ct}{4}\right),&\text{if}\quad|x-ct|\leq
2\pi,\\\ 0,&\text{else.}\end{cases}$

The accuracy of the scheme is measured away from the interference of the non-
smooth interfaces in the interval $[0,2\pi]$ at $T=\pi/2$ with the periodic
boundary condition. Figure 12 and Table 4 present a comparison of the errors
in the $L^{\infty}$\- and $L^{1}$-norms between the WENO-Z scheme and WENO-E
scheme with different values of $\lambda\Delta x$ $(0.02,0.04,0.06,0.1)$.
WENO-E scheme with $\lambda\Delta x=0.02$ produces better accuracy as plotted
in Figure 11.

(a) $T=0$

(b) $T=10$

(c) $T=50$

(d) $T=120$

(e) Wave simulation across various T

(f) Contour for wave simulation

Figure 9. Numerical solution of K(2,2) equation with $N=400$ at
$T=0,10,50,120$ and $\lambda\Delta x=0.02$ WENO-E in $x\in[-5\pi,25\pi]$ for
initial condition (5.55) of Example 5.5.

(a) $T=10$

(b) $T=25$

(c) $T=50$

(d) Wave simulation across various T

(e) Contour for wave simulation

Figure 10. Numerical solution of K(3,3) equation with $N=600$ at $T=10,25,50$
and $\lambda\Delta x=0.02$ WENO-E in $x\in[0,30\pi]$ for Example 5.6.

(a) Numerical solution

(b) Zoom in near smooth extrema

Figure 11.  Comparison of numerical solutions obtained from using WENO-Z and
WENO-E schemes with $N=320$ at $T=\pi/2$ in $x\in[-4\pi,4\pi]$, for Example
5.5

(a) (a) Convergence plot: $L^{\infty}$\- error

(b) (b) Convergence plot: $L^{1}$\- error

Figure 12. Comparison WENO-Z and WENO-E schemes in terms of $L^{1}$ and
$L^{\infty}$ errors (in $\log_{10}$ scale) for Example 5.5 at $T=\pi/2$.

(a) $T=0$

(b) $T=2$

(c) $T=6$

(d) $T=8$

(e) Wave simulation across various T

(f) Contour for wave simulation

Figure 13. Numerical solution of K(3,3) equation with $N=400$ at $T=0,2,6,8$
and $\lambda\Delta x=0.02$ WENO-E in $x\in[-6\pi,6\pi]$ for initial condition
(5.57) of Example 5.6

To observe the behavior of the decomposition of compacton into compacton-
anticompacton pairs, we take the initial data

(5.55)
$u(x,0)=\begin{cases}\frac{4}{3}\cos^{2}\left(\frac{x}{8}\right),&\text{if}\quad-4\pi\leq
x\leq 4\pi,\\\ 0,&\text{else.}\end{cases}$

In Figure 9, wave motion is simulated at $T=0,10,25,$ and $50$ in the domain
of $[-5\pi,25\pi]$ with $N=400$ cells. The compactons are observed to split
from one initial compacton and move to the right over time. A small residue is
also developed at the left interface, similar to the LDG scheme [6, 9], but no
Gibbs oscillations occur in the non-smooth interface (i.e., the edges of the
compactons). The small residue, which appears to be a compacton-anti-compacton
pair on the left side of the compacton-packet, is believed to be physical and
was also detected by the LDG scheme [6, 9]. Our scheme can differentiate
between numerical and physical oscillation without any additional special
treatment.

######  Example 5.6.

Consider the $K(3,3)$ equation

(5.56) $u_{t}+(u^{3})_{x}+(u^{3})_{xxx}=0.$

In order to observe the interaction between three compactons, we take the
initial data as

$u(x,0)=\begin{cases}\sqrt{3}\cos\left(\frac{x-10}{3}\right),&\text{if}\quad|x-10|\leq\frac{3\pi}{2},\\\
\frac{3}{2}\cos\left(\frac{x-25}{3}\right),&\text{if}\quad|x-25|\leq\frac{3\pi}{2},\\\
\sqrt{\frac{3}{2}}\cos\left(\frac{x-40}{3}\right),&\text{if}\quad|x-40|\leq\frac{3\pi}{2},\\\
0,&\text{otherwise.}\end{cases}$

The exact solution is given by

Figure 10 shows that as the compactons propagate to the right, three
compactons with different speeds ($c=10,25,40$) pass through each other during
non-linear interaction while maintaining their coherent shapes after the
collision. Although the compactons emerge intact from the collision, a minor
residue is reflected back from the collision on the left. Similar findings can
also be found in [6, 9]. These compactons are discovered to be not fully
elastic, as mentioned in [5] and identified in the original compacton study
[7]. The wavelets of the residue become even more apparent with mesh
refinement. Based on earlier compacton studies and numerical observations, we
conclude that this phenomenon is not numerically induced.

Next, we solve (5.6) subject to the initial data

(5.57)
$u(x,0)=\begin{cases}\sqrt{\frac{3}{2}}\cos\left(\frac{x}{6}\right),&\text{if}\quad-3\pi\leq
x\leq 3\pi,\\\ 0,&\text{else.}\end{cases}$

The results of our numerical simulations are shown in Figure 13. As time
evolves, a train of canonical compactons splits from the initial data and
moves to the right. At the same time, a rapid oscillation develops at the left
interface of the initial data. The oscillations that we get on the left side
of the solution indicate that the oscillations are an integral part of the
solution and not a numerical artifact.

After examining the six aforementioned examples, it becomes evident that there
is no universally optimal value for the tension parameter, denoted as
$\lambda\Delta x$. Its effectiveness is contingent upon the specific
characteristics of the equation under consideration. In linear scenarios,
regardless of whether the number of grid points $N$ is small or large, we
observe improved accuracy within the range of $\lambda\Delta x$ between 0.02
and 0.04. Conversely, in non-linear cases, superior accuracy can be achieved
by setting $\lambda\Delta x$ to 0.1 for low values of $N$ and maintaining the
range of $\lambda\Delta x$ between 0.02 and 0.04 for high values of $N$. More
details are provided in Appendix:A. In summary, it is advisable to employ the
range, $\lambda\Delta x\in[0.02,0.04]$, to ensure the optimal performance of
the proposed scheme.

##  6\. Conclusion

In this study, a fifth-order WENO scheme based on exponential polynomials for
solving the non-linear dispersion-type equations has been proposed. The
proposed approximation space utilizes the exponential basis with a tension
parameter that can be optimized to better fit the specific features of the
characteristics of the initial data, resulting in improved results without the
presence of spurious oscillations. Numerical tests were conducted to simulate
several equations with various initial conditions, including the 1-D and 2-D
linear and non-linear KdV equations, $K(2,2)$ equation, and $K(3,3)$ equation.
In direct comparison to the WENO-Z method, our proposed method WENO-E, when
executed with a optimized tension parameter, has superior capacity for
resolving non-linear dispersion-type equations at higher resolutions.

## Acknowledgements

The author Samala Rathan is supported by NBHM, DAE, India (Ref. No.
02011/46/2021 NBHM(R.P.)/R $\&$ D II/14874) and IIPE, Visakhapatnam, India
(IRG Grant No. IIPE/DORD/IRG/001).

## Conflict of interest

The authors declare no potential conflict of interests.

## Data availability

The data that support the findings of this study are available upon reasonable
request.

## Appendix:A

In this section, we provide a detailed tabulated values that enumerates the
$L^{1}$-errors associated with a range of $\lambda\Delta x$ values in Table 5,
6, 7, and 8 for examples 5.1, 5.2, 5.3, and 5.5 respectively. The chosen
spectrum of $\lambda\Delta x$ values is carefully selected to ensure that
there is a balanced and adequate representation of errors. This comprehensive
overview allows us to effectively showcase the impact and significance of
varying $\lambda\Delta x$. The optimal value of $\lambda\Delta x$ are
highlighted for each example in bold. To see a universally applicable
parameter, we’ve categorized $\lambda\Delta x$ values into three groups which
are as follows:

  * •

Category A (Cat-A) is $\lambda\Delta x\in[0.02,0.04]$ which comprises optimal
values that consistently yield accuracy across various grid point counts $N$.

  * •

Category B (Cat-B) is $\lambda\Delta x\in[0.06,0.1]$ which includes values
closely resembling Cat-A at lower grid points $N$, but they neither improve
nor maintain accuracy as $N$ increases, ultimately deteriorating.

  * •

Category C (Cat-C) is $\lambda\Delta x>0.1$ which encompasses values that
exhibit inconsistency and substantial deviation from the actual solution,
numerically distant from Cat-A.

As per these observations, the choice of the tension parameter, $\lambda\Delta
x$, does not possess a universally optimal value, as its effectiveness relies
on the specific characteristics of the equation under consideration. In linear
scenarios, whether the number of grid points $N$ is small or large, we observe
enhanced accuracy within the $\lambda\Delta x$ range of $0.02$ to $0.04$.
Conversely, in nonlinear cases, superior accuracy is achieved by setting
$\lambda\Delta x$ to $0.1$ for low $N$ and adhering to the $\lambda\Delta x$
range of $0.02$ to $0.04$ for high $N$. In overall, one can utilize the Cat-A
i.e., $\lambda\Delta x\in[0.02,0.04]$ to see the good performance of the
proposed scheme.

Table 5. WENO-E schemes in terms of $L^{1}$\- errors for Example 5.1 over the
domain $\Omega=[0,2\pi]$ at time $T=1$.
$\bm{L^{1}}$ | $\bm{N}$
---|---
error | 10 | 20 | 40 | 80 | 160 | 320 | 640
WENO-Z | 1.7456e-03 | 5.6856e-05 | 1.7815e-06 | 5.5652e-08 | 1.7390e-09 | 7.0742e-11 | 1.3677e-10
WENO-E-0.01 | 1.7519e-03 | 5.7101e-05 | 1.7825e-06 | 5.5640e-08 | 1.7317e-09 | 6.7075e-11 | 1.3520e-10
WENO-E-0.02 | 1.7519e-03 | 5.7100e-05 | 1.7821e-06 | 5.5421e-08 | 1.6221e-09 | 1.3413e-11 | 1.0871e-10
WENO-E-0.03 | 1.7519e-03 | 5.7097e-05 | 1.7802e-06 | 5.4471e-08 | 1.1474e-09 | 2.2517e-10 | 3.5957e-11
WENO-E-0.04 | 1.7520e-03 | 5.7089e-05 | 1.7751e-06 | 5.1911e-08 | 1.3076e-10 | 8.6372e-10 | 3.3622e-10
WENO-E-0.05 | 1.7520e-03 | 5.7069e-05 | 1.7643e-06 | 4.6515e-08 | 2.8258e-09 | 2.2103e-09 | 1.0080e-09
WENO-E-0.06 | 1.7520e-03 | 5.7032e-05 | 1.7447e-06 | 3.6701e-08 | 7.7266e-09 | 4.6589e-09 | 2.2314e-09
WENO-E-0.07 | 1.7520e-03 | 5.6970e-05 | 1.7124e-06 | 2.0541e-08 | 1.5797e-08 | 8.6913e-09 | 4.2465e-09
WENO-E-0.08 | 1.7520e-03 | 5.6874e-05 | 1.6629e-06 | 4.2486e-09 | 2.8176e-08 | 1.4877e-08 | 7.3376e-09
WENO-E-0.09 | 1.7519e-03 | 5.6734e-05 | 1.5908e-06 | 4.0299e-08 | 4.6179e-08 | 2.3872e-08 | 1.1833e-08
WENO-E-0.1 | 1.7516e-03 | 5.6537e-05 | 1.4902e-06 | 9.0593e-08 | 7.1295e-08 | 3.6421e-08 | 1.8105e-08
WENO-E-0.2 | 1.7378e-03 | 4.7842e-05 | 2.8973e-06 | 2.2839e-06 | 1.1666e-06 | 5.8367e-07 | 2.9160e-07
WENO-E-0.3 | 1.6711e-03 | 1.0023e-05 | 2.1903e-05 | 1.1783e-05 | 5.9103e-06 | 2.9538e-06 | 1.4761e-06
WENO-E-0.4 | 1.4862e-03 | 9.1817e-05 | 7.3016e-05 | 3.7329e-05 | 1.8667e-05 | 9.3276e-06 | 4.6615e-06
WENO-E-0.5 | 1.0921e-03 | 3.0621e-04 | 1.8055e-04 | 9.1071e-05 | 4.5503e-05 | 2.2736e-05 | 1.4356e-03
WENO-E-0.6 | 3.7362e-04 | 6.9462e-04 | 3.7529e-04 | 1.8838e-04 | 9.4091e-05 | 4.7011e-05 | 1.4620e+32
Table 6. WENO-E schemes in terms of $L^{1}$\- errors for Example 5.2 over the domain $\Omega=(0,2\pi)\times(0,2\pi)$ at time $T=1$. $\bm{L^{1}}$ | $\bm{N_{x}\times N_{y}}$
---|---
error | 10$\times$10 | 20$\times$20 | 40$\times$40 | 80$\times$80 | 160$\times$160
WENO-Z | 3.4560e-03 | 1.1104e-04 | 3.5318e-06 | 1.1080e-07 | 3.4732e-09
WENO-E-0.01 | 3.4272e-03 | 1.1149e-04 | 3.5338e-06 | 1.1077e-07 | 3.4586e-09
WENO-E-0.02 | 3.4272e-03 | 1.1149e-04 | 3.5329e-06 | 1.1034e-07 | 3.2400e-09
WENO-E-0.03 | 3.4273e-03 | 1.1149e-04 | 3.5292e-06 | 1.0845e-07 | 2.2923e-09
WENO-E-0.04 | 3.4273e-03 | 1.1147e-04 | 3.5191e-06 | 1.0335e-07 | 2.5903e-10
WENO-E-0.05 | 3.4273e-03 | 1.1143e-04 | 3.4977e-06 | 9.2607e-08 | 5.6387e-09
WENO-E-0.06 | 3.4273e-03 | 1.1136e-04 | 3.4589e-06 | 7.3069e-08 | 1.5421e-08
WENO-E-0.07 | 3.4272e-03 | 1.1124e-04 | 3.3949e-06 | 4.0895e-08 | 3.1531e-08
WENO-E-0.08 | 3.4270e-03 | 1.1105e-04 | 3.2966e-06 | 8.4582e-09 | 5.6242e-08
WENO-E-0.09 | 3.4267e-03 | 1.1078e-04 | 3.1537e-06 | 8.0231e-08 | 9.2177e-08
WENO-E-0.1 | 3.4262e-03 | 1.1039e-04 | 2.9544e-06 | 1.8036e-07 | 1.4231e-07
WENO-E-0.2 | 3.3973e-03 | 9.3414e-05 | 5.7438e-06 | 4.5470e-06 | 2.3286e-06
WENO-E-0.3 | 3.2631e-03 | 1.9555e-05 | 4.3422e-05 | 2.3459e-05 | 1.1798e-05
WENO-E-0.4 | 2.8944e-03 | 1.7936e-04 | 1.4476e-04 | 7.4321e-05 | 3.7262e-05
WENO-E-0.5 | 2.1115e-03 | 5.9823e-04 | 3.5799e-04 | 1.8133e-04 | 9.0833e-05
WENO-E-0.6 | 6.8564e-04 | 1.3574e-03 | 7.4421e-04 | 3.7510e-04 | 1.8783e-04
Table 7. WENO-E schemes in terms of $L^{1}$\- errors for Example 5.3 over the domain $\Omega=[-10,10]$ at time $T=0.5$. $\bm{L^{1}}$ | $\bm{N}$
---|---
error | 40 | 80 | 160 | 320 | 640 | 1280
WENO-Z | 3.1369e-02 | 2.3854e-02 | 5.9718e-05 | 1.9878e-06 | 8.4259e-08 | 2.5115e-08
WENO-E-0.01 | 3.0119e-02 | 2.1970e-03 | 6.8675e-05 | 2.0650e-06 | 8.4833e-08 | 2.5116e-08
WENO-E-0.02 | 3.0119e-02 | 2.1970e-03 | 6.8676e-05 | 2.0648e-06 | 8.4750e-08 | 2.5080e-08
WENO-E-0.03 | 3.0120e-02 | 2.1971e-03 | 6.8676e-05 | 2.0642e-06 | 8.4390e-08 | 2.4925e-08
WENO-E-0.04 | 3.0120e-02 | 2.1972e-03 | 6.8676e-05 | 2.0624e-06 | 8.3420e-08 | 2.4566e-08
WENO-E-0.05 | 3.0120e-02 | 2.1973e-03 | 6.8672e-05 | 2.0584e-06 | 8.1376e-08 | 2.4272e-08
WENO-E-0.06 | 3.0121e-02 | 2.1974e-03 | 6.8663e-05 | 2.0510e-06 | 7.7660e-08 | 2.5510e-08
WENO-E-0.07 | 3.0122e-02 | 2.1975e-03 | 6.8645e-05 | 2.0389e-06 | 7.1545e-08 | 2.8296e-08
WENO-E-0.08 | 3.0122e-02 | 2.1977e-03 | 6.8615e-05 | 2.0201e-06 | 6.2423e-08 | 3.2753e-08
WENO-E-0.09 | 3.0123e-02 | 2.1978e-03 | 6.8571e-05 | 1.9928e-06 | 5.1088e-08 | 3.9343e-08
WENO-E-0.1 | 3.0124e-02 | 2.1979e-03 | 6.8506e-05 | 1.9547e-06 | 4.4961e-08 | 4.8688e-08
WENO-E-0.2 | 3.0136e-02 | 2.1967e-03 | 6.5444e-05 | 7.2544e-07 | 8.6967e-07 | 4.7576e-07
WENO-E-0.3 | 3.0145e-02 | 2.1807e-03 | 5.1841e-05 | 7.1637e-06 | 4.5905e-06 | 2.3462e-06
WENO-E-0.4 | 3.0135e-02 | 2.1291e-03 | 2.2468e-05 | 2.6852e-05 | 1.4602e-05 | 7.3805e-06
WENO-E-0.5 | 3.0085e-02 | 2.0125e-03 | 7.3678e-05 | 6.8390e-05 | 3.5670e-05 | 1.7974e-05
WENO-E-0.6 | 2.9965e-02 | 1.7948e-03 | 2.1361e-04 | 1.4370e-04 | 7.3836e-05 | 4.3984e-05
Table 8. WENO-E schemes in terms of $L^{1}$\- errors for Example 5.5 over the domain $\Omega=[0,2\pi]$ at time $T=\pi/2$. $\bm{L^{1}}$ | $\bm{N}$
---|---
error | 40 | 80 | 160 | 320 | 640 | 1280
WENO-Z | 3.1608e-04 | 1.1146e-05 | 4.5757e-07 | 2.0152e-08 | 1.4181e-09 | 1.5958e-10
WENO-E-0.01 | 3.1654e-04 | 1.1150e-05 | 4.5757e-07 | 2.0149e-08 | 1.4175e-09 | 1.5902e-10
WENO-E-0.02 | 3.1654e-04 | 1.1150e-05 | 4.5747e-07 | 2.0105e-08 | 1.4092e-09 | 1.5226e-10
WENO-E-0.03 | 3.1654e-04 | 1.1149e-05 | 4.5702e-07 | 1.9910e-08 | 1.3782e-09 | 1.4107e-10
WENO-E-0.04 | 3.1654e-04 | 1.1146e-05 | 4.5581e-07 | 1.9391e-08 | 1.3336e-09 | 2.1720e-10
WENO-E-0.05 | 3.1653e-04 | 1.1141e-05 | 4.5326e-07 | 1.8373e-08 | 1.4339e-09 | 5.1372e-10
WENO-E-0.06 | 3.1652e-04 | 1.1131e-05 | 4.4862e-07 | 1.6655e-08 | 2.1872e-09 | 1.0973e-09
WENO-E-0.07 | 3.1649e-04 | 1.1115e-05 | 4.4098e-07 | 1.4819e-08 | 3.9253e-09 | 2.0677e-09
WENO-E-0.08 | 3.1644e-04 | 1.1091e-05 | 4.2927e-07 | 1.3932e-08 | 6.8062e-09 | 3.5647e-09
WENO-E-0.09 | 3.1638e-04 | 1.1055e-05 | 4.1223e-07 | 1.7225e-08 | 1.1160e-08 | 5.7471e-09
WENO-E-0.1 | 3.1629e-04 | 1.1006e-05 | 3.8847e-07 | 2.5432e-08 | 1.7251e-08 | 8.7921e-09
WENO-E-0.2 | 3.1212e-04 | 8.8400e-06 | 7.5150e-07 | 5.5608e-07 | 2.8347e-07 | 1.4158e-07
WENO-E-0.3 | 2.9387e-04 | 2.2379e-06 | 5.4111e-06 | 2.8729e-06 | 1.4366e-06 | 7.1667e-07
WENO-E-0.4 | 2.4463e-04 | 2.5784e-05 | 1.7949e-05 | 9.1032e-06 | 4.5375e-06 | 2.2632e-06
WENO-E-0.5 | 1.4122e-04 | 7.8889e-05 | 4.4326e-05 | 2.2210e-05 | 1.1061e-05 | 5.5166e-06
WENO-E-0.6 | 5.1057e-05 | 1.7505e-04 | 9.2092e-05 | 4.5940e-05 | 2.2871e-05 | 1.1407e-05

## Appendix:B

The matrix $A$ is

$\displaystyle\left[\begin{smallmatrix}1&-2\Delta x&\frac{17\Delta
x^{2}}{4}&\biggl{(}-\frac{4e^{-2\lambda\Delta x}\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta x^{3}}+\frac{4e^{-2\lambda\Delta
x}\cosh[\lambda\Delta x]\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{5\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{2\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\cos[2\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}-\frac{8\sin[2\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ 1&-\Delta x&\frac{5\Delta
x^{2}}{4}&\biggl{(}-\frac{4e^{-\lambda\Delta x}\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta x^{3}}+\frac{4e^{-\lambda\Delta
x}\cosh[\lambda\Delta x]\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{-\frac{\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{3\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\cos[\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}-\frac{8\sin[\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ 1&0&\frac{\Delta
x^{2}}{4}&\biggl{(}-\frac{4\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta x^{3}}+\frac{4\cosh[\lambda\Delta
x]\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{\frac{-3\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\sin{\left[\frac{\lambda\Delta
x}{2}\right]}^{3}}{\lambda^{3}\Delta x^{3}}\biggr{)}&0\\\ 1&\Delta
x&\frac{5\Delta x^{2}}{4}&\biggl{(}-\frac{4e^{\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}+\frac{4e^{\lambda\Delta x}\cosh[\lambda\Delta
x]\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{\frac{-5\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{-\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{-\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\cos[\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\sin[\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ 1&2\Delta x&\frac{17\Delta
x^{2}}{4}&\biggl{(}-\frac{4e^{2\lambda\Delta x}\sinh\left[\frac{2\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta x^{3}}+\frac{4e^{2\lambda\Delta
x}\cosh[2\lambda\Delta x]\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{\frac{-7\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{-3\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{-2\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\cos[2\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\sin[2\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ 1&3\Delta x&\frac{37\Delta
x^{2}}{4}&\biggl{(}-\frac{4e^{3\lambda\Delta x}\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta x^{3}}+\frac{4e^{3\lambda\Delta
x}\cosh[\lambda\Delta x]\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{\frac{-9\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{-5\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{-3\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\cos[3\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\sin[3\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}\\\ 1&4\Delta x&\frac{65\Delta
x^{2}}{4}&\biggl{(}-\frac{4e^{4\lambda\Delta x}\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta x^{3}}+\frac{4e^{4\lambda\Delta
x}\cosh[\lambda\Delta x]\sinh\left[\frac{\lambda\Delta
x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{e^{\frac{-11\lambda\Delta
x}{2}}(-1+e^{\lambda\Delta x})}{\lambda^{3}\Delta
x^{3}}+\frac{e^{\frac{-7\lambda\Delta x}{2}}(-1+e^{\lambda\Delta
x})}{\lambda^{3}\Delta x^{3}}-\frac{4e^{-4\lambda\Delta
x}\sinh\left[\frac{\lambda\Delta x}{2}\right]}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\cos[4\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}&\biggl{(}\frac{8\sin[4\lambda\Delta
x]\sin{\left[\frac{\lambda\Delta x}{2}\right]}^{3}}{\lambda^{3}\Delta
x^{3}}\biggr{)}\end{smallmatrix}\right],$

the vectors $B$ and $X$ are

$B=\left[g^{+}(u_{i-2}),g^{+}(u_{i-1}),g^{+}(u_{i}),g^{+}(u_{i+1}),g^{+}(u_{i+2}),g^{+}(u_{i+3}),g^{+}(u_{i+4})\right]^{T},$
$X=\left[a_{0},a_{1},a_{2},a_{3},a_{4},a_{5},a_{6}\right]^{T}.$

## References

  * [1] A. M. Wazwaz, Partial differential equations and solitary waves theory, Springer Science and Business Media, 2010.
  * [2] N. J. Zabusky, M. D. Kruskal, Interaction of ” solitons” in a collisionless plasma and the recurrence of initial states, Phys. Rev. Lett., 15(6) (1965), pp. 240.
  * [3] P. Rosenau, J. M. Hyman, Compactons: solitons with finite wavelength, Phys. Rev. Lett., 70(5) (1993), pp. 564.
  * [4] A. Chertock, D. Levy, Particle methods for dispersive equations, J. Comput. Phys., 171(2) (2001), pp. 708-730.
  * [5] B. Mihaila, A. Cardenas, F. Cooper, A. Saxena, Stability and dynamical properties of Rosenau-Hyman compactons using Padé approximants, Phys. Rev. E, 81(5) (2010), pp. 056708.
  * [6] D. Levy, C. W. Shu, J. Yan, Local discontinuous Galerkin methods for nonlinear dispersive equations, J. Comput. Phys., 196 (2) (2004), pp. 751–772.
  * [7] M.S. Ismail, T. R. Taha, A numerical study of compactons, Math. Comput. Simul., 47 (6) (1998), pp. 519–530.
  * [8] P. Saucez, A. V. Wouwer, W. E. Schiesser, P. Zegeling, Method of lines study of nonlinear dispersive waves, J. Comput. Appl. Math. 168 (1–2) (2004), pp. 413–423.
  * [9] M. Ahmat, J. Qiu, Direct WENO scheme for dispersion-type equations, Math. Comput. Simul., 204 (2023), pp. 216-229.
  * [10] X. D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys., 115 (1994), pp. 200-212.
  * [11] G.S. Jiang, C. W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys., 126 (1996), pp. 202–228.
  * [12] C. W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys., 77 (1988), pp. 439-471.
  * [13] S. Gottlieb, D. I. Ketcheson, C. W. Shu, Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations, World Scientific, 2011.
  * [14] A. Mazaheri, H. Nishikawa, Very efficient high-order hyperbolic schemes for time-dependent advection–diffusion problems: Third-, fourth-, and sixth-order, Comput. Fluids, 102 (2014), pp. 131–147.
  * [15] A. Mazaheri, M. Ricchiuto, H. Nishikawa, A first-order hyperbolic system approach for dispersion, J. Comput. Phys., 321 (Supplement C) (2016), pp. 593–605.
  * [16] A. K. Henrick, T. D. Aslam, J. M. Powers, Mapped weighted-essentially-non-oscillatory schemes: Achieving optimal order near critical points, J. Comput. Phys., 207 (2005), pp. 542–567.
  * [17] R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys., 227 (2008), pp. 3191–3211.
  * [18] D. S. Balsara, C. W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. J. Comput. Phys., 160 (2000), pp.405-452.
  * [19] M. Castro, B. Costa, W. S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys., 230 (2011), pp.1766–1792.
  * [20] Y. Ha, C. H. Kim, Y. J. Lee, J. Yoon, An improved weighted essentially non-oscillatory scheme with a new smoothness indicator, J. Comput. Phys., 232 (2013a), pp. 68–86.
  * [21] S. Rathan, G. Naga Raju, A modified fifth-order WENO scheme for hyperbolic conservation laws, Comput. Math. Appl., 75(5), (2018):, pp.1531-1549.
  * [22] S. Rathan, G. Naga Raju, An improved non-linear weights for seventh-order weighted essentially non-oscillatory scheme, Comput. Fluids. 156 (2017), pp. 496-514.
  * [23] S. Rathan, G. Naga Raju, Improved weighted ENO scheme based on parameters involved in nonlinear weights, Appl. Math. Comput. 331 (2018) 120-129
  * [24] Y. Ha, C. Kim, H. Yang, J. Yoon, Sixth-order weighted essentially nonoscillatory schemes based on exponential polynomials, SIAM J. Sci. Comput., 38 (2016), pp. A1984–A2017.
  * [25] Y. Ha, C. H. Kim, H. Yang, J. Yoon, A sixth-order weighted essentially non-oscillatory schemes based on exponential polynomials for Hamilton–Jacobi equations, J. Sci. Comput., 75 (2018), pp. 1675–1700.
  * [26] Y. Ha, Y. J. Lee, J. Yoon, Modified essentially non-oscillatory scheme based on exponential polynomial interpolation for hyperbolic conservation laws, SIAM J. Numer. Anal., 52 (2013b), pp. 864–893.
  * [27] J. Zhu J. Qiu, Trigonometric WENO schemes for hyperbolic conservation laws and highly oscillatory problems, Commun. Comput. Phys., 8 (2010), pp. 1242–1263.
  * [28] R. Abedian, M. Dehghan, A high-order weighted essentially nonoscillatory scheme based on exponential polynomials for nonlinear degenerate parabolic equations, Numer. Methods Partial Differ. Eq., 38 (2022), pp. 970–996.
  * [29] S. Rathan, R. Kumar, A. D. Jagtap, L1-type smoothness indicators based WENO scheme for nonlinear degenerate parabolic equations, Appl. Math. Comput. 375 (2020) 125112.
  * [30] S. Rathan, J. Gu, A sixth-order central WENO scheme for nonlinear degenerate parabolic equations, Comput. Appl. Math. 42 (2023) 182.
  * [31] S. Karlin, W. J. Studden, Tchebycheff systems: With applications in analysis and statistics, Interscience Publishers, (1966).
  * [32] F. Aràndiga, A. Baeza, A. M. Belda, P. Mulet, Analysis of WENO schemes for full and global accuracy. SIAM J. Numer. Anal., 49(2) (2011), pp.893-915.
  * [33] W. S. Don, R. Borges, Accuracy of the weighted essentially non-oscillatory conservative finite difference schemes, J. Comput. Phys., 250 (2013), pp.347-372.