WENO interpolations and reconstructions using data bounded polynomial approximation

Solution of partial differential equations (PDEs) like hyperbolic conservation laws and convection dominated problems admits strong irregular profiles (large jumps,discontinuities) along with complicated smooth structures which make numerical treatment of such solution challenging. It is worth mentioning that the visualization of the solution of a PDE model is an essential aid to comprehend the physical problem. The goodness of any numerical scheme for these PDEs relies on the accuracy and non-oscillatory nature of approximated numerical fluxes at cell interfaces. Also for numerical stability, scheme should be data-bounded to ensure the bounded growth of the numerical solutions [17]. This consequently depends on boundededness of the polynomial approximation used in underlying numerical scheme. Thus high-order well-behaved polynomial approximation techniques that can interpolate both smooth and discontinuous function data without introducing “Gibbs-like” oscillations are much required. In other words, a high order polynomial approximations should respect the physical properties like boundedness and shape of given data [28]. In fact, data bounded polynomial approximations which are bounded by the minimum and maximum data values may also help to preserve positivity in the solution [29]. Throughout this text, until stated, the term ”data bounded or bounded polynomial approximations” represents that it is bounded by the minimum and maximum of data values used to construct underlying polynomial.

Among high order polynomial approximation techniques essentially non-oscillatory (ENO) and the weighted ENO (WENO) approximations achieved phenomenal popularity due to their ability to approximate the discrete data corresponding to strong discontinuities in an essentially non-oscillatory manner while retaining high-order accuracy in smooth data regions. Though ENO or WENO procedures are generic function approximation techniques, they are mainly developed in the context of numerical schemes for hyperbolic conservation laws (HCL).

Harten et. al [16] first proposed Essentially non-oscillatory(ENO) schemes which later got efficiently implemented by Shu and Osher in [18, 19]. ENO procedure is based on the idea of smoothest stencil selection out of many candidate stencil. In [1], Liu et al., proposed the WENO procedure which uses convex combination of all the candidate stencils unlike choosing the smoothest one in the ENO techniques. Later, Jiang and Shu in [2] introduced the most-celebrated finite-difference WENO framework popularly named as WENO-JS scheme by modifying the smoothness measurement and extended the scheme up to fifth-order accuracy. In [3] Henrick et al., proposed a mapping function to improve optimal order of accuracy at critical points as WENO-JS schemes fail to achieve the formal order of accuracy at critical points. Borges et al., [4] developed new WENO weights by using global smoothness measurement for fifth-order WENO scheme and named it as WENO-Z scheme. Further many versions of WENO schemes are developed in [5, 6, 7, 8, 9, 10, 11, 12, 13, 20]. General information about WENO framework can be found in the surveys by Shu [14] and Zhang et al. [15]. Generally, WENO schemes are based upon WENO reconstruction procedure which generates a piecewise polynomial approximations of a function from given set of its cell averages. WENO reconstruction using polynomials from cell averages is equivalent to the interpolation of point values of the primitive polynomial [14]. Therefore, WENO methods (algorithms) can also be formulated as an interpolation technique [14]. In WENO methods the use of high order polynomials are attractive due to their potential of yielding high accurate approximation of cell interface values.

Despite of the substantial body of work and developments on WENO methods, there are still challenges associated with the stability analysis of these methods. Some of the work carried out in terms of stabilities of WENO schemes are [22, 23, 24] and analysis of a peculiar sign stability property [25] is done for third order WENO reconstruction in [26, 27]. However, to the best of authors knowledge the stability of WENO schemes in terms of boundedness of the solution or associated WENO approximations are not carried out so far.

The main aim of this work is to establish conditions on data-boundedness of three point polynomial approximation and further construct high order data-bounded WENO interpolations and WENO reconstructions at cell interfaces.

The rest of the paper is organized as follows: In Section 2 the data-boundedness conditions for a three-point polynomial approximations is established. This is done by analyzing the conditions on non-linear weights such that the three point polynomial approximation is bounded. These conditions are provided in terms of smoothness parameters which are ratio of consecutive gradients. Section 3 provides suitable non-negative weights for data-bounded WENO approximations along with computational verification for the accuracy and data-boundedness of the WENO approximations. Section 4 briefly discuss the application of data-bounded WENO approximations to construct WENO schemes for hyperbolic conservation laws following author’s work [30], whereas other possible uses are mentioned in 5. Section 6 contains the concluding remarks.

Let the set of discrete data values $\{v_{j}\}_{j=i-1}^{j=i+1}$ be given for a bounded function $v(x)$ at evenly distributed spatial points $\displaystyle\{x_{j}\}_{j=1-1}^{j=i+1}$ with fixed grid spacing $\Delta x$. These data values can be either point values $v_{j}=v(x_{j})$ or cell-average values $\overline{v}_{j}$ of $v(x)$ in cell $I_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ i.e., $\displaystyle\overline{v}_{j}=\frac{1}{\Delta x}\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}v(x)dx$. Consider the following three-point polynomial approximation,

$$\hat{v}(x)=\sum_{p=0}^{1}\alpha_{p}\hat{v}^{p}(x)=\alpha_{0}\hat{v}^{0}(x)+\alpha_{1}\hat{v}^{1}(x)$$

(2.1)

where weights $\alpha_{i}\in\mathbb{R},i=0,1$ are such that $\alpha_{1}=1-\alpha_{0}$ and $\hat{v}^{p}$ are the following linear interpolating polynomial using data $\{v_{i+p-1},v_{i+p}\},\;p=0,1$ in stencils $S_{p}(i):=\{x_{i+p-1},x_{i+p}\}$

	$\displaystyle\hat{v}^{0}(x)$	$\displaystyle=v_{i-1}+\frac{v_{i}-v_{i-1}}{\Delta x}(x-x_{i-1})$		(2.2)
	$\displaystyle\hat{v}^{1}(x)$	$\displaystyle=v_{i}+\frac{v_{i+1}-v_{i}}{\Delta x}(x-x_{i})$

Define the quantities $L_{i}^{+}=\frac{1}{1-r_{i}^{+}},~{}L_{i}^{-}=\frac{1}{1-r_{i}^{-}}$, where $r_{i}^{\pm}$ is the smoothness parameter defined as

$$r_{i}^{\pm}=\frac{\Delta_{\mp}v_{i}}{\Delta_{\pm}v_{i}}\in\mathbb{R}\cup\{\pm\infty\}$$

(2.3)

Also let $m_{i}=\min\{v_{i-1},v_{i},v_{i+1}\}\;\mbox{and}\;M_{i}=\max\{v_{i-1},v_{i},v_{i+1}\}$ then boundedness conditions of the polynomial approximation (2.1) in terms of weight $\alpha_{0}$ and smoothness parameter $r_{i}$, are given by the following theorem.

Proof. Consider the case when smoothness parameter is defined as $r_{i}^{+}$. The proof for $r_{i}^{-}$ follows analogously.
Case I: Monotone Data $r_{i}^{+}\geq 0$: By rearranging the terms, polynomial approximation (2.1) can be written in the following forms,

	$$\hat{v}(x)=v_{i-1}+(1-\alpha_{0})\Delta_{-}v_{i}+\frac{\alpha_{0}}{\Delta x}(x-x_{i-1})\Delta_{-}v_{i}+\frac{\alpha_{1}}{\Delta x}(x-x_{i})\Delta_{+}v_{i}$$		(2.4a)
	$$\hat{v}(x)=v_{i+1}-\alpha_{0}\left(1-\frac{x-x_{i-1}}{\Delta x}\right)\Delta_{-}v_{i}-\left(1-\frac{\alpha_{1}}{\Delta x}(x-x_{i})\right)\Delta_{+}v_{i}$$		(2.4b)
	$$\hat{v}(x)=v_{i}-\alpha_{0}\left(1-\frac{x-x_{i-1}}{\Delta x}\right)\Delta_{-}v_{i}+\frac{\alpha_{1}}{\Delta x}(x-x_{i})\Delta_{+}v_{i}$$		(2.4c)

Let data be monotonically increasing i.e., $m_{i}=v_{i-1}\leq v_{i}\leq v_{i+1}=M_{i}$ then it follows from (2.4a),

$\displaystyle\hat{v}(x)$

$\displaystyle\geq$

$\displaystyle v_{i-1},~{}\text{if}~{}(1-\alpha_{0})\Delta_{-}v_{i}+\frac{\alpha_{0}}{h}\Delta_{-}v_{i}(x-x_{i-1})+\frac{\alpha_{1}}{h}\Delta_{+}v_{i}(x-x_{i})\geq 0$

or

$$\hat{v}(x)\geq m_{i}~{}~{}~{}~{}\text{if}~{}\alpha_{0}(x-x_{i})(1-r_{i}^{+})\leq r_{i}^{+}\Delta x+(x-x_{i})$$

(2.5)

Similarly from (2.4b)

$\displaystyle\hat{v}(x)$

$\displaystyle\leq$

$\displaystyle v_{i+1}~{}(=M_{i}),~{}\text{if}~{}\alpha_{0}(x-x_{i})(1-r_{i}^{+})\geq x-x_{i+1}$

(2.6)

together with (2.5) and (2.6);

$$m_{i}\leq\hat{v}(x)\leq M_{i},\;\mbox{provided}\;(x-x_{i+1})\leq\alpha_{0}(x-x_{i})(1-r_{i}^{+})\leq\left(r_{i}^{+}\Delta x+(x-x_{i})\right).$$

(2.7)

Conditions (a) and (b) of theorem follows from the conditional compound inequality in (2.7). Similarly these conditions for monotonically decreasing data i.e., $M_{i}=u_{i-1}\geq u_{i}\geq u_{i+1}=m_{i}$ can be obtained.
Case 2: Non Monotone data $r_{i}^{+}\leq 0$: It implies that given data consists extrema. Here calculations are given only for the case when $v_{i}$ is minima as similar calculations follows for the case $v_{i}$ is maxima. In case of minima ,

$$v_{i}\leq\min(v_{i-1},v_{i+1}):\Rightarrow m_{i}=v_{i},M_{i}=\max(v_{i+1},v_{i-1}),\mbox{and}\;\Delta_{-}v_{i}\leq 0,~{}\Delta_{+}v_{i}\geq 0.$$

Following sub-cases can occur

• (i)

  $v_{i-1}\leq v_{i+1}\Rightarrow M_{i}=v_{i+1}~{}\mbox{and}~{}v_{i}-v_{i-1}\geq v_{i}-v_{i+1}\Rightarrow r_{i}^{+}\geq-1.$ Therefore $r_{i}^{+}\in[-1,0]$ and $1\leq(1-r_{i}^{+})\leq 2.$ It follows from the equation (2.4c)

  	$\displaystyle\hat{v}(x)$	$\displaystyle\geq$	$\displaystyle v_{i},~{}\text{if}~{}\alpha_{0}\left(1-\frac{x-x_{i-1}}{\Delta x}\right)\Delta_{-}v_{i}-\frac{\alpha_{1}}{\Delta x}(x-x_{i})\Delta_{+}v_{i}\leq 0$		(2.8)
  		$\displaystyle\geq$	$\displaystyle m_{i}~{}\text{if}~{}\alpha_{0}(x-x_{i})\leq\frac{x-x_{i}}{1-r_{i}^{+}}$

  Similarly using (2.4b)

  $\displaystyle\hat{v}(x)$

  $\displaystyle\leq$

  $\displaystyle v_{i+1}~{}(=M_{i}),~{}\text{if}~{}\alpha_{0}(x-x_{i})\geq\frac{x-x_{i+1}}{1-r_{i}^{+}}$

  (2.9)

  Equation (2.8) and (2.9) establish condition (c) of the theorem as,

  $$m_{i}=v_{i}\leq\hat{v}(x)\leq v_{i+1}=M_{i},\;\mbox{provided}\;\frac{x-x_{i+1}}{1-r_{i}^{+}}\leq\alpha_{0}(x-x_{i})\leq\frac{x-x_{i}}{1-r_{i}^{+}}.$$

• (ii)

  $v_{i-1}\geq v_{i+1}\Rightarrow M_{i}=v_{i-1}$ and $v_{i}-v_{i-1}<v_{i}-v_{i+1}\Rightarrow r_{i}^{+}\leq-1$ and $(1-r_{i}^{+})\geq 2.$ Now from (2.4a) we get,

  $\displaystyle\hat{v}(x)$

  $\displaystyle\leq$

  $\displaystyle v_{i-1}~{}(=M_{i}),~{}\text{if}~{}\alpha_{0}(x-x_{i})\leq\frac{r_{i}^{+}\Delta x+(x-x_{i})}{(1-r_{i}^{+})}$

  (2.10)

  Equation (2.8) and (2.10) implies condition (d) of the theroem i.e.,

  $$m_{i}=v_{i}\leq\hat{v}(x)\leq v_{i-1}=M_{i},\;\mbox{provided}\;\frac{r_{i}^{+}\Delta x+(x-x_{i})}{(1-r_{i}^{+})}\leq\alpha_{0}(x-x_{i})\leq\frac{x-x_{i}}{1-r_{i}^{+}}.$$

In this section high order data-bounded WENO approximation of the function $v(x)$ at the cell-interface $x_{i+\frac{1}{2}}$ is discussed. Note that weights $\beta_{0}$ and $\mu_{0}$ defined in (2.13) and (2.16) respectively give WENO approximation. However to ensure high order accuracy in WENO approximations at cell interfaces in sooth solution region, the non-linear weights $\beta_{0}$ and $\mu_{0}$ must attain ideal non-linear weight [30],

The semi-discretized conservative scheme for the one-dimensional hyperbolic conservation laws $u_{t}+f(u)_{x}=0$ is as follows

$$\displaystyle\frac{du_{i}(t)}{dt}=-\frac{1}{\Delta x}\left(\hat{f}_{i+\frac{1}{2}}-\hat{f}_{i-\frac{1}{2}}\right)$$

(4.1)

In the third-order WENO (WENO3) scheme [21] the numerical flux $\hat{f}_{i+\frac{1}{2}}$ is of the form

$$\hat{f}_{i+\frac{1}{2}}=\sum_{l=0}^{1}\omega_{l}\hat{f}^{l}_{i+\frac{1}{2}},$$

(4.2)

where non-linear weights $\omega_{l}$ satisfy the following convexity property

$$\sum_{l=0}^{1}\omega_{l}=1,\;\omega_{l}\geq 0.$$

(4.3)

and the expression of $\hat{f}^{l}_{i+\frac{1}{2}}$ are as follows

	$\displaystyle\hat{f}^{0}_{i+\frac{1}{2}}=\frac{3}{2}f_{i}-\frac{1}{2}f_{i-1}$		(4.4)
	$\displaystyle\hat{f}^{1}_{i+\frac{1}{2}}=\frac{1}{2}f_{i}+\frac{1}{2}f_{i+1}.$

Here the point values of the physical flux of hyperbolic conservation laws i.e., $f(u_{i})=f_{i}$ are considered as the given data values and for maintaining the accuracy of WENO3 schemes, the ideal or linear weights are $d_{0}=\frac{1}{3},~{}~{}d_{1}=\frac{2}{3}.$

• •

  The data-bounded WENO approximations can be applied for several real-life PDE based problems which contains both strong discontinuities and complex smooth solution features. For example, DB-WENO approximations can be applied to the the Convection dominated [14] and Hamilton-Jacobi problems [32] to develop data-bounded schemes for those problems.

• •

  It is important to note that WENO method is actually an approximation procedure, not directly related to PDEs, hence the DB-WENO procedure can also be used in many non-PDE applications, including computer vision and image processing.

Conditions on non-linear weights for constructing data-bounded polynomial approximation are obtained. Based on these bounds, weights are constructed to ensure the data-boundedness of third and fourth order WENO approximations (interpolation and reconstructions) for smooth solution data. Numerical results to show the accuracy and boundedness of the WENO approximations are also given. The fourth order data-bounded WENO approximation which is described in the section 3, can be used to make data-bounded fourth order WENO (WENO4) scheme, but for that detailed charaterization of the structure of WENO weights is required. This work is under investigation and can be reported separately in future.

Acknowledgment: Authors acknowledge the Science and Engineering Board, New Delhi, India for providing necessary financial support through funded projects File No. EMR/2016/000394.

Appendix A: Other Non-linear DB-weights

• $\bullet$

  Some proposed non-linear DB-weights for the region (2.13) in Corrolary 2.1.

  • i.)

    $$\beta_{0}^{\eta}=\eta K$$

    (6.1)

    where $0\leq\eta\leq 1$ and $K$ is defined in equation (2.13).

  • ii.)

    $$\beta_{0}^{2}=\left\{\begin{array}[]{cc}1/4&if~{}r_{i}^{+}\in[-3,5],\\ \frac{8}{3{r_{i}^{+}}^{2}+5}&if~{}r_{i}^{+}<-3,\\ \frac{5}{3{r_{i}^{+}}^{2}-55}&if~{}r_{i}^{+}>5.\end{array}\right.$$

    (6.2)

• $\bullet$

  Some proposed non-linear DB-weights for the region (2.16) in Corrolary 2.2.

  • i.)

    $$\mu_{0}^{\eta}=1-\eta(1-J)$$

    (6.3)

    where $0\leq\eta\leq 1$ and $J$ is defined in equation (2.16).

  • ii.)

    $$\mu_{0}^{2}=\left\{\begin{array}[]{cc}3/4&if~{}r_{i}^{-}\in[-3,5],\\ \frac{3({r_{i}^{-}}^{2}-1)}{3{r_{i}^{-}}^{2}+5}&if~{}r_{i}^{-}<-3,\\ \frac{3({r_{i}^{-}}^{2}-20)}{3{r_{i}^{-}}^{2}-55}&if~{}r_{i}^{-}>5.\end{array}\right.$$

    (6.4)

  Note that by the convexity property of weights, the other non-linear weights $\beta_{1}^{\eta,2}=1-\beta_{0}^{\eta,2},\&~{}\mu_{1}^{\eta,2}=1-\mu_{0}^{\eta,2}.$ From Figure 7 it can be seen that different non-linear DB-weights defined by the above equations are lies inside the data-bound region.