Fuxing Hu
Department of Mathematics, Huizhou University, Huizhou, China.
DOI: 10.4236/jamp.2026.142044   PDF    HTML   XML   74 Downloads   282 Views   Citations

Abstract

A strategy is devised to make the WENO interpolation in the point values achieve optimal accuracy near the discontinuities. The classical WENO interpolation ensures the optimal accuracy when all stencils are smooth and ENO property when the discontinuity appears. When there exist more than two successive smooth stencils, the maximum theoretical accuracy near discontinuity is also preferred to be obtained. To achieve it, we divide the classical WENO algorithm into several sub-WENO procedures. In each sub-WENO procedure, only two stencils are used and the order of accuracy grows by at most. If both stencils are smooth, then sub-WENO procedure increases the order of accuracy by one. If there is a stencil that is smooth and the left one is non-smooth, then algorithm conserves the order of interpolation by corresponding smooth stencil and keeps the ENO property. If both stencils are non-smooth, then the value constructed by sub-WENO procedure will be ignored in the latter procedures. The whole of new WENO algorithm can be expressed as a tree structure. The indicator of smoothness of every medium stencil in the tree structure is defined by the indicators of smoothness of corresponding stencils on the top of tree. Such definition is proved to be capable of obtaining the optimal accuracy and keeping the ENO property. And the new WENO algorithm has almost the same computational cost as the classical WENO algorithm.

Keywords

Optimal Order of Accuracy, WENO Interpolation, Corner and Jump Discontinuity

Share and Cite:

1. Introduction

The interpolation of piecewise smooth function from the discrete sample points is always an important problem of numerical computation. The simple linear Lagrangian interpolation is suitable for a set of smooth sample points and achieves the optimal accuracy. However, we often need to interpolate a set of discrete points which contains corner (or jump) discontinuities. In these cases, the Gibbs oscillations generated by linear Lagrangian interpolations around the discontinuities will decrease the numerical accuracy. Then some preferred nonlinear algorithms should be used to obtain stable numerical approximation. The nonlinear essentially non-oscillatory (ENO) interpolation, which can tackle this problem well, was first presented in [1] [2] for solving hyperbolic conservation laws. To eliminate the effect of non-smooth parts of sample points, the ENO algorithm adaptively chooses the smoothest stencil from the candidates to reconstruct the piecewise function. The indicator of smoothness of each stencil is defined as the Newton divided difference of sample points on this stencil. More information about ENO algorithm can be found in [3]-[5].

In the procedure of ENO interpolation, while the smoothness of all the candidate stencils is measured, only the smoothest one is finally conserved even though all the stencils are smooth enough to be used. To gather all the potential stencils, in [6], the authors presented the weighted ENO (WENO) algorithm based on ENO. The aim of WENO algorithm is to keep the ENO property in non-smooth regions and meanwhile improve the order of accuracy in smooth regions. To obtain it, in the WENO algorithm, each candidate stencil is assigned a nonlinear weight which measures the contribution of this stencil for the final convex combination. The nonlinear weights are devised to approach linear optimal weights in smooth region to improve the accuracy of interpolation. While there exist stencils across the discontinuities, the weights assigned to these non-smooth stencils can be ignored almost. In [7], the authors showed new indicators of smoothness which emulate the idea of minimizing the total variation of function. The new indicators of smoothness are defined as the $L_2$ norm of the derivatives of the interpolation polynomials. The indicators of smoothness in [7] ensure the optimal order of accuracy in smooth region and ENO property when there exist stencils that contain the discontinuities. About the improvements and applications of WENO algorithm, the interested reader can refer to [8]-[16].

Around the discontinuities, the classical WENO algorithm [7] generally only achieves the same order of accuracy as the corresponding ENO interpolation. It was shown in [7] for reconstructing the numerical flux at cell boundary based on the cell averages. What leads to the degeneration of order of accuracy is the classical WENO cannot distribute reasonable optimal weights to the smooth stencils when the discontinuity appears. In [17], the authors proposed a power WENO algorithm to improve the accuracy near discontinuities. However, it does not obtain the maximum theoretical order of accuracy. Then, the authors in [18] succeed in obtaining the maximum theoretical accuracy close to the discontinuities by an improved WENO algorithm (WENO-AW). Unlike the fixed linear optimal weights proposed in [7] [19], they devise nonlinear optimal weights to tackle this problem. The nonlinear optimal weights are expressed as the nonlinear convex combination of three vectors of linear optimal weights. Each vector of linear optimal weights is appropriate for one special case. The proofs of maximum theoretical accuracy near discontinuities and ENO property were presented in [18]. They also give another algorithm to raise the accuracy of WENO algorithm for the interval which contains the corner discontinuity. In [17] [18] [20], certain modified indicators of smoothness are used to detect the corner and jump discontinuities. At the same time, these indicators of smoothness conserve the optimal accuracy and ENO property. In this paper, we only consider how to improve the order of accuracy near discontinuities by using a more efficient algorithm. If one wants to obtain the optimal accuracy in the interval containing the discontinuity, the algorithm in [18] will be a commendable choice.

As mentioned above, the algorithm in [18] recovers the optimal accuracy near discontinuities by devising a set of nonlinear optimal weights. However, the new nonlinear optimal weights are computed by an extra WENO algorithm. And in this procedure, the indicators of smoothness of the bigger stencils also need to be calculated. As shown in [18], the computational costs of new WENO algorithm are more than double when compared with the classical WENO. In this paper, we present a simple WENO algorithm to recover the optimal accuracy near discontinuities. In order to describe this algorithm clearly, we take the 6th-order WENO algorithm for example. To ensure the optimal accuracy near discontinuities, the classical WENO algorithm is divided into three sub-WENO procedures. We first construct two 5th-order WENO approximations by using the former two 4-points stencils and latter two 4-points stencils, respectively. In the next, we construct the 6th-order WENO approximation by using two 5-points stencils and corresponding 5th-order WENO approximations which have been obtained. From the statement above, we need three sub-WENO procedures in this algorithm and in each sub-WENO procedure we only need information of two stencils. Unlike the classical WENO, which directly constructs the 6th-order approximation by the nonlinear convex combination of three 4th-order approximations, the order of accuracy of new WENO algorithm grows one by one and this nonlinear interpolation truly confirms the optimal accuracy near discontinuities. Since the 5-points stencils are used in this procedure, we have to compute their indicators of smoothness. To control the oscillations and reduce the computational cost, we express the indicators of smoothness of 5-points stencils as the product of two corresponding 4-points substencils. The new WENO algorithm here can be reformulated into the similar compact form as the classical WENO algorithm and not much computation is introduced. Furthermore, it is also easy to extend the algorithm to the cases of higher order.

The organization of this paper is as follows. In Section 2, we review the classical WENO [7] and WENO-AW [18] algorithms. In Section 3, we show the new WENO interpolation and prove the statement of optimal order of accuracy near discontinuity and property of ENO. In Section 5, we test the numerical accuracy and computational costs.

2. Review of WENO Algorithms

In this section, we review the classical WENO [7] and WENO-AW algorithm [18] in point values. The interested reader can also refer to [21] [22] for a full statement of WENO interpolations.

2.1. The Classical WENO Algorithm with Fixed Optimal Weights

Let us consider a set of sample points $\left(x_i,f_i\right)$ , $1\le i\le N$ , where $f_i=f\left(x_i\right)$ and $\text{Δ}x=x_i-x_{i-1}$ . What we want to do is to interpolate the value in middle point $x_{j-\frac{1}{2}}$ of interval $\left(x_{j-1},x_j\right)$ when the discontinuity appears around this interval, but not in this interval. As shown in [23], it is possible to locate the corner discontinuities, but it is no hope to locate the jump discontinuities. Hence, when the corner discontinuity appears in $\left(x_{j-1},x_j\right)$ , the algorithm proposed in [18] can be used to tackle this problem. But when the jump discontinuity appears in $\left(x_{j-1},x_j\right)$ , the order of accuracy of numerical approximation to $x_{j-\frac{1}{2}}$ will be affected inevitably.

Since $x_{j-\frac{1}{2}}$ we want to approximate is in the interval $\left(x_{j-1},x_j\right)$ , each stencil used should contain $\left(x_{j-1},x_j\right)$ . Let us denote by $S_i^m$ the stencil

$\left\{x_{j+i-m+1},\cdots ,x_{j-1},x_j,\cdots ,x_{j+i}\right\},$

where the superscript of $S_i^m$ denotes the number of point of this stencil contains and the subscript denotes the number of point at the right of interval $\left(x_{j-1},x_j\right)$ . The same notation is also used for the interpolation polynomial $p_i^m$ , indicators of smoothness ${\beta }_i^m$ , the weights $w_i^m$ and so on. Let us take the 6th-order WENO for example to approximate the value at $x_{j-\frac{1}{2}}$ in the middle of $\left(x_{j-1},x_j\right)$ . The 6th-order WENO algorithm uses three 4-points stencils,

$\begin{array}{l}{S}_0^4=\left\{x_{j-3},x_{j-2},x_{j-1},x_j\right\},\\ S_1^4=\left\{x_{j-2},x_{j-1},x_j,x_{j+1}\right\},\\ S_2^4=\left\{x_{j-1},x_j,x_{j+1},x_{j+2}\right\}.\end{array}$

On each stencil, we obtain approximation at $x_{j-\frac{1}{2}}$ by a Lagrangian interpolation polynomial of degree 3,

$\begin{array}{l}{p}_0^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(f_{j-3}-5f_{j-2}+15f_{j-1}+5f_j\right),\\ p_1^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(-f_{j-2}+9f_{j-1}+9f_j-f_{j+1}\right),\\ p_2^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(5f_{j-1}+15f_j-5f_{j+1}+f_{j+2}\right).\end{array}$ (1)

In [7], the indicators of smoothness were defined as the $L_2$ norm of the derivatives of the interpolation polynomials. These indicators of smoothness were proposed for upwind methods to solve hyperbolic conservation laws with cell averages. To fit the interpolation in point values, we adopt the one presented in [14] [18] [20],

${\beta }_i^4={\displaystyle \sum }_{l=2}^n{\displaystyle {\int }_{x_{j-1}}^{x_j}{\left(\frac{{\text{d}}^l}{\text{d}{x}^l}{p}_i^4\left(x\right)\right)}^2\text{d}x},i=0,1,2.$ (2)

By removing the first derivative of interpolation polynomial from formula in [7], the modified indicators of smoothness perform well for the corner discontinuities. The explicit forms of (2) of three stencils are

$\begin{array}{l}{\beta }_0^4=\frac{1}{48}{\left(8f_{j-3}-27f_{j-2}+30f_{j-1}-11f_j\right)}^2+\frac{13}{16}{\left(f_{j-2}-2f_{j-1}+f_j\right)}^2,\\ {\beta }_1^4=\frac{1}{48}{\left(8f_{j-2}-21f_{j-1}+18f_j-5f_{j+1}\right)}^2+\frac{13}{16}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2,\\ {\beta }_2^4=\frac{1}{48}{\left(11f_{j-1}-30f_j+27f_{j+1}-8f_{j+2}\right)}^2+\frac{13}{16}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2.\end{array}$ (3)

To ensure the optimal accuracy and ENO property, the nonlinear weights are introduced by

$w_i^4=\frac{{\alpha }_i^4}{{{\displaystyle \sum }}_{l=0}^2{\alpha }_l^4},{\alpha }_i^4=\frac{d_i^4}{{\left(ϵ+{\beta }_i^4\right)}^q},i=0,1,2.$ (4)

The nonlinear weights $w_i^4\ge 0$ and ${\displaystyle {\sum }_{i=0}^2{w}_i^4}=1$ . The linear optimal weights $d_i^4$ are chosen to be $d_0^4=\frac{3}{16}$ , $d_1^4=\frac{10}{16}$ , $d_2^4=\frac{3}{16}$ so that

$p_2^6\left(x_{j-\frac{1}{2}}\right)=d_0^4{p}_0^4\left(x_{j-\frac{1}{2}}\right)+d_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right)+d_2^4{p}_2^4\left(x_{j-\frac{1}{2}}\right).$

The parameter $ϵ$ appears in denominator is used to avoid the division by zero. The fully discusses about $ϵ$ can refer to [8] [24]. The exponent $q$ in (4) is used to increase the difference of scales of weights near the non-smooth region. For high order WENO interpolations, we generally need to choose $q\ge 2$ for control the numerical oscillations around the discontinuities [19] [25]. To imitate the 6th-order linear interpolation at the smooth region and meanwhile compress the interpolation oscillations, the approximation $F_2^6\left(x_{j-\frac{1}{2}}\right)$ of $p_2^6\left(x\right)$ at $x_{j-\frac{1}{2}}$ is chosen to be

$F_2^6\left(x_{j-\frac{1}{2}}\right)=w_0^4{p}_0^4\left(x_{j-\frac{1}{2}}\right)+w_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right)+w_2^4{p}_2^4\left(x_{j-\frac{1}{2}}\right).$ (5)

At smooth regions, Taylor series expansions at $x_{j-\frac{1}{2}}$ of the indicators of smoothness in (3) can be collected to be

${\beta }_i^4={\left(\Delta x^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\Delta x^2\right)\right),i=0,1,2.$ (6)

Replacing ${\beta }_i^4$ in (4) and taking $ϵ$ small enough, if $f^{″}\left(x_{j-\frac{1}{2}}\right)\ne 0$ the nonlinear weights approximate linear ones by

$w_i^4=d_i^4+\mathcal{O}\left(\text{Δ}{x}^2\right),i=0,1,2.$ (7)

Substituting (7) into (5) gives

$\begin{array}{c}{F}_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum }_{i=0}^2{w}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum }_{i=0}^2{d}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i^4-d_i^4\right)p_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i^4-d_i^4\right)p_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum }_{i=0}^2\left(w_i^4-d_i^4\right)f_{j-\frac{1}{2}}+{\displaystyle \sum }_{i=0}^2{d}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i^4-d_i^4\right)\left(p_i^4\left(x_{j-\frac{1}{2}}\right)-f_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ =p_2^6\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\Delta x^6\right)\\ =f_{j-\frac{1}{2}}+\mathcal{O}\left(\Delta x^6\right).\end{array}$ (8)

At the smooth part of discretized data, the classical 6th-order WENO algorithm has optimal accuracy if $f^{″}\left(x_{j-\frac{1}{2}}\right)\ne 0$ .

In the next, we consider the case in which ${\beta }_0^4$ and ${\beta }_1^4$ are smooth, while ${\beta }_2^4$ contains a jump discontinuity. For the corner discontinuity, we can reach the similar conclusions. In the case of jump discontinuity, three indicators of smoothness will take values

${\beta }_0^4=\mathcal{O}\left(\Delta x^4\right),{\beta }_1^4=\mathcal{O}\left(\Delta x^4\right),{\beta }_2^4=\mathcal{O}\left(1\right).$

Then the corresponding weights can be expressed as

$\begin{array}{c}{w}_0^4=\frac{\frac{d_0^4}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{d_0^4}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{d_1^4}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{d_2^4}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{d_0^4}{d_0^4+d_1^4{\left(\frac{ϵ+{\beta }_0^4}{ϵ+{\beta }_1^4}\right)}^q+d_2^4{\left(\frac{ϵ+{\beta }_0^4}{ϵ+{\beta }_2^4}\right)}^q}\\ =\frac{d_0^4}{d_0^4+d_1^4{\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right)}^q+d_2^4\left(\mathcal{O}\left(\text{Δ}{x}^{4q}\right)\right)}=\frac{d_0^4}{d_0^4+d_1^4}+\mathcal{O}\left(\text{Δ}{x}^2\right),\end{array}$

$\begin{array}{c}{w}_1^4=\frac{\frac{d_1^4}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{d_0^4}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{d_1^4}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{d_2^4}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{d_1^4}{d_0^4{\left(\frac{ϵ+{\beta }_1^4}{ϵ+{\beta }_0^4}\right)}^q+d_1^4+d_2^4{\left(\frac{ϵ+{\beta }_1^4}{ϵ+{\beta }_2^4}\right)}^q}\\ =\frac{d_1^4}{d_0^4{\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right)}^q+d_1^4+d_2^4\left(\mathcal{O}\left(\text{Δ}{x}^{4q}\right)\right)}=\frac{d_1^4}{d_0^4+d_1^4}+\mathcal{O}\left(\text{Δ}{x}^2\right)\end{array}$

and

$\begin{array}{c}{w}_2^4=\frac{\frac{d_2^4}{{\left(ϵ+{\beta }_2^4\right)}^q}}{\frac{d_0^4}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{d_1^4}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{d_2^4}{{\left(ϵ+{\beta }_2^4\right)}^q}}\\ =\frac{d_2^4\frac{{\left(ϵ+{\beta }_0^4\right)}^q{\left(ϵ+{\beta }_1^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q}}{d_0^4{\left(ϵ+{\beta }_1^4\right)}^q+d_1^4{\left(ϵ+{\beta }_0^4\right)}^q+d_2^4\frac{{\left(ϵ+{\beta }_0^4\right)}^q{\left(ϵ+{\beta }_1^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\mathcal{O}\left(\text{Δ}{x}^{4q}\right).\end{array}$

So, the contribution of stencil $S_2^4$ which contains discontinuity can be ignored when $\text{Δ}x$ is small enough. It is exactly the ENO property. Generally, since

$\frac{d_0^4}{d_0^4+d_1^4}{p}_0^4\left(x_{j-\frac{1}{2}}\right)+\frac{d_1^4}{d_0^4+d_1^4}{p}_1^4\left(x_{j-\frac{1}{2}}\right)\ne p_1^5\left(x_{j-\frac{1}{2}}\right)$

and using the similar operations in (8), we only obtain

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f_{j-\frac{1}{2}}+\mathcal{O}\left(\text{Δ}{x}^4\right).$ (9)

Finally, we consider the case of ${\beta }_0^4=\mathcal{O}\left(\text{Δ}{x}^4\right)$ , ${\beta }_1^4=\mathcal{O}\left(1\right)$ and ${\beta }_2^4=\mathcal{O}\left(1\right)$ . That is, the discontinuity lies in interval $\left(x_j,x_{j+1}\right)$ and only stencil ${\beta }_0^4$ is smooth. Through the similar analysis as the above cases, the same conclusion as (9) is achieved for WENO approximation at $x_{j-\frac{1}{2}}$ . The idea of classical WENO algorithm is to ensure the optimal accuracy at smooth part of discretized data and ENO property when there exist stencils affected by discontinuity. In the second case, since ${\beta }_0^4$ and ${\beta }_1^4$ are both smooth, the higher 5th-order approximation using the information of stencil $S_1^5$ is preferred to be reached.

2.2. The WENO Algorithm with Adapted Optimal Weights

In the classical WENO algorithm [7], the optimal weights $d_i^4$ ($i=0,1,2$ ) in Equation (4) are fixed. It is considered as the primary reason of degeneration of order of accuracy around discontinuities in [18]. The authors present the adapted optimal weights based on the indicators of smoothness of stencils. For the 6th-order WENO-AW algorithm, if the jump discontinuity lies in $\left(x_{j+1},x_{j+2}\right)$ , i.e., ${\beta }_0^4=\mathcal{O}\left(\text{Δ}{x}^4\right)$ , ${\beta }_1^4=\mathcal{O}\left(\text{Δ}{x}^4\right)$ and ${\beta }_2^4=\mathcal{O}\left(1\right)$ , then optimal weights $\left(2d_0^4,d_1^4,0\right)$ , which satisfy

$2d_0^4{p}_0^4\left(x_{j-\frac{1}{2}}\right)+d_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right)=p_1^5\left(x_{j-\frac{1}{2}}\right),$

are preferred to be chosen to achieve the 5th-order of accuracy. Similarly, if the discontinuity lies in $\left(x_{j-3},x_{j-2}\right)$ , i.e., ${\beta }_0^4=\mathcal{O}\left(1\right)$ , ${\beta }_1^4=\mathcal{O}\left(\text{Δ}{x}^4\right)$ and ${\beta }_2^4=\mathcal{O}\left(\text{Δ}{x}^4\right)$ , then optimal weights $\left(0,d_1^4,2d_2^4\right)$ , which satisfy

$d_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right)+2d_2^4{p}_2^4\left(x_{j-\frac{1}{2}}\right)=p_2^5\left(x_{j-\frac{1}{2}}\right),$

are preferred to be chosen. To improve the accuracy near discontinuities, an adapted strategy to choose the optimal weights among $\left(2d_0^4,d_1^4,0\right)$ , $\left(0,d_1^4,2d_2^4\right)$ and $\left(d_0^4,d_1^4,d_2^4\right)$ is proposed in [18]. In order to avoid the abrupt transition from one vector of optimal weight to another one at the interfaces between smooth and non-smooth regions, the nonlinear adapted optimal weights are defined as a smooth convex combination of three vectors of optimal weights,

$\left({\tilde{d}}_0^4,{\tilde{d}}_1^4,{\tilde{d}}_2^4\right)={\tilde{w}}_1^5\left(2d_0^4,d_1^4,0\right)+{\tilde{w}}_2^6\left(d_0^4,d_1^4,d_2^4\right)+{\tilde{w}}_2^5\left(0,d_1^4,2d_2^4\right).$ (10)

The coefficients of three vectors of linear optimal weights are defined as

${\tilde{w}}_1^5=\frac{{\tilde{\alpha }}_1^5}{{\tilde{\alpha }}_1^5+{\tilde{\alpha }}_2^6+{\tilde{\alpha }}_2^5},{\tilde{w}}_2^6=\frac{{\tilde{\alpha }}_2^6}{{\tilde{\alpha }}_1^5+{\tilde{\alpha }}_2^6+{\tilde{\alpha }}_2^5},{\tilde{w}}_2^5=\frac{{\tilde{\alpha }}_2^5}{{\tilde{\alpha }}_1^5+{\tilde{\alpha }}_2^6+{\tilde{\alpha }}_2^5},$

where

${\tilde{\alpha }}_1^5=\frac{1}{{\left(ϵ+{\beta }_1^5\right)}^q},{\tilde{\alpha }}_2^6=\frac{1}{{\left(ϵ+{\beta }_2^6\right)}^q},{\tilde{\alpha }}_2^5=\frac{1}{{\left(ϵ+{\beta }_2^5\right)}^q}.$

After obtaining the adapted optimal weights $\left({\tilde{d}}_0^4,{\tilde{d}}_1^4,{\tilde{d}}_2^4\right)$ , inserting them into Equation (4) gives

$w_i^4=\frac{{\alpha }_i^4}{{{\displaystyle \sum }}_{l=0}^2{\alpha }_l^4},{\alpha }_i^4=\frac{{\tilde{d}}_i^4}{{\left(ϵ+{\beta }_i^4\right)}^q},i=0,1,2.$(11)

In case 1: all the 4-points stencils are smooth. The nonlinear weights computed by (11) satisfy

$\left(w_0^4,w_1^4,w_2^4\right)=\left(d_0^4,d_1^4,d_2^4\right)+\mathcal{O}\left(\text{Δ}{x}^2\right).$

It is the sufficient condition, as shown by the Theorem 3.2 in [18], to achieve the 6th-order approximation,

$F_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum }_{i=0}^2{w}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)=f_{j-\frac{1}{2}}+\mathcal{O}\left(\text{Δ}{x}^6\right).$

In case 2: the discontinuity lies in the interval $\left(x_{j+1},x_{j+2}\right)$ , i.e., ${\beta }_0^4$ and ${\beta }_1^4$ are both smooth and ${\beta }_2^4$ is non-smooth. If it is a corner discontinuity, then we have

$\left(w_0^4,w_1^4,w_2^4\right)=\left(2d_0^4,d_1^4,0\right)+\left(\mathcal{O}\left(\text{Δ}{x}^2\right),\mathcal{O}\left(\text{Δ}{x}^2\right),\mathcal{O}\left(\text{Δ}{x}^{2p}\right)\right).$

If it is a jump discontinuity, then we have

$\left(w_0^4,w_1^4,w_2^4\right)=\left(2d_0^4,d_1^4,0\right)+\left(\mathcal{O}\left(\text{Δ}{x}^4\right),\mathcal{O}\left(\text{Δ}{x}^4\right),\mathcal{O}\left(\text{Δ}{x}^{4p}\right)\right).$

If the discontinuity lies in the symmetric interval $\left(x_{j-3},x_{j-2}\right)$ of $\left(x_{j+1},x_{j+2}\right)$ , by similar analysis we obtain

$\left(w_0^4,w_1^4,w_2^4\right)=\left(0,d_1^4,2d_2^4\right)+\left(\mathcal{O}\left(\text{Δ}{x}^{2p}\right),\mathcal{O}\left(\text{Δ}{x}^2\right),\mathcal{O}\left(\text{Δ}{x}^2\right)\right)$

for the corner discontinuity and

$\left(w_0^4,w_1^4,w_2^4\right)=\left(0,d_1^4,2d_2^4\right)+\left(\mathcal{O}\left(\text{Δ}{x}^{4p}\right),\mathcal{O}\left(\text{Δ}{x}^4\right),\mathcal{O}\left(\text{Δ}{x}^4\right)\right)$

for the jump discontinuity. Also by the Theorem 3.2 in [18], the 5th-order approximation can be achieved,

$F_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum }_{i=0}^2{w}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)=f_{j-\frac{1}{2}}+\mathcal{O}\left(\text{Δ}{x}^5\right).$

In case 3: the discontinuity lies in the interval $\left(x_j,x_{j+1}\right)$ , i.e., ${\beta }_0^4$ is smooth and ${\beta }_1^4$ and ${\beta }_2^4$ are both non-smooth. If it is a corner discontinuity, then we have

$\left(w_0^4,w_1^4,w_2^4\right)=\left(1,0,0\right)+\mathcal{O}\left(\text{Δ}{x}^{2q}\right).$

If it is a jump discontinuity, then we have

$\left(w_0^4,w_1^4,w_2^4\right)=\left(1,0,0\right)+\mathcal{O}\left(\text{Δ}{x}^{4q}\right).$

If the discontinuity lies in the symmetric interval $\left(x_{j-2},x_{j-1}\right)$ of $\left(x_j,x_{j+1}\right)$ , by similar analysis we obtain

$\left(w_0^4,w_1^4,w_2^4\right)=\left(0,0,1\right)+\mathcal{O}\left(\text{Δ}{x}^{2q}\right)$

for the corner discontinuity and

$\left(w_0^4,w_1^4,w_2^4\right)=\left(0,0,1\right)+\mathcal{O}\left(\text{Δ}{x}^{4q}\right)$

for the jump discontinuity. Also by the Theorem 3.2 in [18], the 4th-order approximation can be achieved,

$F_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum }_{i=0}^2{w}_i^4{p}_i^4\left(x_{j-\frac{1}{2}}\right)=f_{j-\frac{1}{2}}+\mathcal{O}\left(\text{Δ}{x}^4\right).$

When the interval $\left[x_{j-1},x_j\right]$ contains the corner discontinuity, the strategy in [18] is suggested to be used. This strategy ensures the 3rd-order of accuracy for the interval which contains the corner discontinuity.

3. The 6th-Order Sub-WENO Algorithm

In this section, a simple 6th-order sub-WENO algorithm is presented, which has the same aim as the WENO-AW algorithm proposed in [18], to achieve the optimal order near the discontinuities. The second algorithm in [18] is to improve the order of accuracy of when the middle interval $\left(x_{j-1},x_j\right)$ contains the discontinuity. As shown in Section 2.2, the extra three indicators of smoothness of big stencils ${\beta }_1^5$ , ${\beta }_2^5$ and ${\beta }_2^6$ are required to calculate the adapted nonlinear optimal weights. The computational cost of the WENO-AW algorithm in [18] is more than double when compared with the classical WENO algorithm [7]. In the classical 6th-order WENO algorithm, the idea is to combine three 4th-order linear interpolations to achieve 6th-order accuracy at smooth regions and the potential 5th-order interpolations are skipped. It is the probable reason which leads to the degeneration of algorithm near discontinuities. To recover the accuracy near discontinuities, we divide the classical WENO into several sub-WENO procedures. In each sub-WENO procedure, only two stencils are used. The order of accuracy grows by only one at most by combining two stencils and the ENO property is also conserved (Figure 1).

Figure 1. The left is the tree structure of the 6th-order sub-WENO algorithm. The right is the tree structure of the 6th-order classical WENO algorithm.

In the 6th-order sub-WENO algorithm, there are three 4-points stencils can be used and there will be three sub-WENO procedures as shown in Section 3. In the first sub-WENO procedure denoted by $\left\{S_0^4,S_1^4\right\}$ ↪$S_1^5$ , emulating the classical WENO in Section 2, we use the stencils $S_0^4$ and $S_1^4$ to compute the 5th-order approximation at $x_{j-\frac{1}{2}}$ ,

$\begin{array}{l}{F}_1^5\left(x_{j-\frac{1}{2}}\right)=w_0^4{p}_0^4\left(x_{j-\frac{1}{2}}\right)+w_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right).\hfill \end{array}$ (12)

We extract the formulas of interpolation polynomials $p_i^4\left(x\right)\left(i=0,1\right)$ at $x_{j-\frac{1}{2}}$ from (1),

$\begin{array}{l}{p}_0^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(f_{j-3}-5f_{j-2}+15f_{j-1}+5f_j\right),\\ p_1^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(-f_{j-2}+9f_{j-1}+9f_j-f_{j+1}\right).\end{array}$ (13)

The nonlinear weights are computed similarly as (4) by

$\begin{array}{l}{w}_0^4=\frac{{\alpha }_0^4}{{\alpha }_0^4+{\alpha }_1^4},w_1^4=\frac{{\alpha }_1^4}{{\alpha }_0^4+{\alpha }_1^4},\hfill \end{array}$ (14)

and the unnormalized weights are defined as

$\begin{array}{l}{\alpha }_0^4=\frac{d_0^4}{{\left(ϵ+{\beta }_0^4\right)}^q},{\alpha }_1^4=\frac{d_1^4}{{\left(ϵ+{\beta }_1^4\right)}^q}.\hfill \end{array}$ (15)

The linear optimal weights in (15) are chose to be $d_0^4=\frac{3}{8}$ and $d_1^4=\frac{5}{8}$ to satisfy

$p_1^5\left(x_{j-\frac{1}{2}}\right)=d_0^4{p}_0^4\left(x_{j-\frac{1}{2}}\right)+d_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right).$

The indicators of smoothness ${\beta }_i^4\left(i=0,1\right)$ are extracted from (3), and we repeat them for completeness of new algorithm,

$\begin{array}{l}{\beta }_0^4=\frac{1}{48}{\left(8f_{j-3}-27f_{j-2}+30f_{j-1}-11f_j\right)}^2+\frac{13}{16}{\left(f_{j-2}-2f_{j-1}+f_j\right)}^2,\\ {\beta }_1^4=\frac{1}{48}{\left(8f_{j-2}-21f_{j-1}+18f_j-5f_{j+1}\right)}^2+\frac{13}{16}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2.\end{array}$ (16)

Up to now, the first sub-WENO procedure $\left\{S_0^4,S_1^4\right\}$ ↪$S_1^5$ is completed by emulating the classical WENO algorithm. This sub-WENO procedure can achieve 5th-order approximation to $f_{j-\frac{1}{2}}$ at smooth region and keep the ENO property when there is one non-smooth stencil. Actually, when there is one non-smooth stencil, this sub-WENO procedure still ensure the 4th-order of accuracy. The proofs of assertions are ignored here, since the properties of the algorithm after integrating three sub-WENO procedures are our aims.

In the second sub-WENO procedure $\left\{S_1^4,S_2^4\right\}$ ↪$S_2^5$ , we use the stencils $S_1^4$ and $S_2^4$ to obtain another 5th-order approximation to $f_{j-\frac{1}{2}}$ ,

$\begin{array}{l}{F}_2^5\left(x_{j-\frac{1}{2}}\right)=w_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right)+w_2^4{p}_2^4\left(x_{j-\frac{1}{2}}\right),\hfill \end{array}$ (17)

where the values of two 4th-order interpolation polynomials at $x_{j-\frac{1}{2}}$ are

$\begin{array}{l}{p}_1^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(-f_{j-2}+9f_{j-1}+9f_j-f_{j+1}\right),\\ p_2^4\left(x_{j-\frac{1}{2}}\right)=\frac{1}{16}\left(5f_{j-1}+15f_j-5f_{j+1}+f_{j+2}\right).\end{array}$ (18)

The nonlinear weights in (17) are computed as

$\begin{array}{l}{w}_1^4=\frac{{\alpha }_1^4}{{\alpha }_1^4+{\alpha }_2^4},w_2^4=\frac{{\alpha }_2^4}{{\alpha }_1^4+{\alpha }_2^4},\hfill \end{array}$ (19)

and the unnormalized weights are defined as

$\begin{array}{l}{\alpha }_1^4=\frac{d_1^4}{{\left(ϵ+{\beta }_1^4\right)}^q},{\alpha }_2^4=\frac{d_2^4}{{\left(ϵ+{\beta }_2^4\right)}^q}.\hfill \end{array}$ (20)

The linear optimal weights are chose to be $d_1^4=\frac{5}{8}$ and $d_2^4=\frac{3}{8}$ for satisfying

$p_2^5\left(x_{j-\frac{1}{2}}\right)=d_1^4{p}_1^4\left(x_{j-\frac{1}{2}}\right)+d_2^4{p}_2^4\left(x_{j-\frac{1}{2}}\right).$

The indicators of smoothness ${\beta }_i^4\left(i=1,2\right)$ are also chosen from (3),

$\begin{array}{l}{\beta }_1^4=\frac{1}{48}{\left(8f_{j-2}-21f_{j-1}+18f_j-5f_{j+1}\right)}^2+\frac{13}{16}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2,\\ {\beta }_2^4=\frac{1}{48}{\left(11f_{j-1}-30f_j+27f_{j+1}-8f_{j+2}\right)}^2+\frac{13}{16}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2.\end{array}$ (21)

The second sub-WENO procedure $\left\{S_1^4,S_2^4\right\}$ ↪$S_2^5$ is similar to the first one but the nonlinear interpolation is operated on the stencils $S_1^4$ and $S_2^4$ . It is noted that in the two sub-WENO procedures above, we use the several same notations, such as the nonlinear weight $w_1^4$ in (14) and (19). But these same notations are independent in the different sub-WENO procedures.

Through the former sub-WENO procedures, we arrive the 5th-order approximations satisfying

$F_1^5\left(x_{j-\frac{1}{2}}\right)=p_1^5\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^5\right)$

and

$F_2^5\left(x_{j-\frac{1}{2}}\right)=p_2^5\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^5\right)$

on the 5-points stencils $S_1^5$ and $S_2^5$ , respectively, when discretized data is smooth. In the third sub-WENO procedure $\left\{S_1^5,S_2^5\right\}$ ↪$S_2^6$ , we arrange to obtain the 6th-order approximation to $f_{j-\frac{1}{2}}$ by combination of $F_1^5\left(x_{j-\frac{1}{2}}\right)$ and $F_2^5\left(x_{j-\frac{1}{2}}\right)$ ,

$F_2^6\left(x_{j-\frac{1}{2}}\right)=w_1^5{F}_1^5\left(x_{j-\frac{1}{2}}\right)+w_2^5{F}_2^5\left(x_{j-\frac{1}{2}}\right).$ (22)

The nonlinear weights are defined as

$\begin{array}{l}{w}_1^5=\frac{{\alpha }_1^5}{{\alpha }_1^5+{\alpha }_2^5},w_2^5=\frac{{\alpha }_2^5}{{\alpha }_1^5+{\alpha }_2^5},\hfill \end{array}$ (23)

where the unnormalized weights are defined as

$\begin{array}{l}{\alpha }_1^5=\frac{d_1^5}{{\left(ϵ+{\beta }_1^5\right)}^q},{\alpha }_2^5=\frac{d_2^5}{{\left(ϵ+{\beta }_2^5\right)}^q}.\hfill \end{array}$ (24)

The optimal weights are chosen to be $d_1^5=\frac{1}{2}$ and $d_2^5=\frac{1}{2}$ for satisfying

$p_2^6\left(x_{j-\frac{1}{2}}\right)=d_1^5{p}_1^5\left(x_{j-\frac{1}{2}}\right)+d_2^5{p}_2^5\left(x_{j-\frac{1}{2}}\right).$

The only problem is how to choose the indicators of smoothness ${\beta }_1^5$ and ${\beta }_2^5$ on $S_1^5$ and $S_2^5$ , respectively. A natural option of ${\beta }_i^5\left(i=1,2\right)$ is to use the Formula (2)

${\beta }_i^5={\displaystyle \sum }_{l=2}^4{\displaystyle {\int }_{x_{j-1}}^{x_j}{\left(\frac{{\text{d}}^l}{\text{d}{x}^l}{p}_i^5\left(x\right)\right)}^2\text{d}x},i=1,2.$

When the discontinuity lies in $\left(x_{j+1},x_{j+2}\right)$ or $\left(x_{j-3},x_{j-2}\right)$ , the choice of (2) is valid and the final algorithm can achieve 5th-order accuracy and keep the ENO property. However, if there exists a discontinuity appears in $\left(x_j,x_{j+1}\right)$ or $\left(x_{j-2},x_{j-1}\right)$ , then $S_1^5$ and $S_2^5$ are both non-smooth. And in result, both stencils are distributed comparative weights and the oscillation will be inevitable. To approach the optimal accuracy and preserve the ENO property, the effective choices are

$\begin{array}{l}{\beta }_1^5:={\beta }_0^4{\beta }_1^4,{\beta }_2^5:={\beta }_1^4{\beta }_2^4.\hfill \end{array}$ (25)

Taking $S_1^5$ for example, since $S_1^5=S_0^4\cup S_1^4$ , we express the indicator of smoothness parent stencil $S_1^5$ as the product of substencils $S_0^4$ and $S_1^4$ . This option will be extended to the higher order sub-WENO algorithms. If the discontinuity appears in $\left(x_j,x_{j+1}\right)$ or $\left(x_{j+1},x_{j+2}\right)$ , by using the indicators of smoothness in (25), the contribution of stencil $S_2^5$ can be ignored and $S_1^5$ will dominate the final combination. The similar result can be obtained when the discontinuity lies in $\left(x_{j-3},x_{j-2}\right)$ or $\left(x_{j-2},x_{j-1}\right)$ . In addition, the indicators of smoothness ${\beta }_i^4\left(i=0,1,2\right)$ have been obtained in the former sub-WENO procedures and we do not need extra computational cost. For the simplification of the equations, we replace $ϵ+{\beta }_i^5\left(i=1,2\right)$ in denominators of (24) by $\left(ϵ+{\beta }_{i-1}^4\right)\left(ϵ+{\beta }_i^4\right)\left(i=1,2\right)$ , respectively, and (24) becomes

$\begin{array}{l}{\alpha }_1^5=\frac{d_1^5}{{\left(\left(ϵ+{\beta }_0^4\right)\left(ϵ+{\beta }_1^4\right)\right)}^q},{\alpha }_2^5=\frac{d_2^5}{{\left(\left(ϵ+{\beta }_1^4\right)\left(ϵ+{\beta }_2^4\right)\right)}^q}.\hfill \end{array}$ (26)

The common term $ϵ+{\beta }_1^4$ will be cancelled when we substitute the unnormalized weights into (23). So, we actually define the indicators of smoothness ${\beta }_1^5$ and ${\beta }_2^5$ as

$\begin{array}{l}{\beta }_1^5:={\beta }_0^4,{\beta }_2^5:={\beta }_2^4.\hfill \end{array}$ (27)

For the higher order sub-WENO algorithms, the choice of indicators of smoothness is presented in Section 4 in detail. The proofs of achieving optimal accuracy and keeping ENO property will be shown in theorem 1. Before giving this theorem, we first simplify the sub-WENO algorithm which will facilitate the proof of theorem 1.

Inserting the weights (14) into (12) gives

$F_1^5\left(x_{j-\frac{1}{2}}\right)=\frac{3{\left(ϵ+{\beta }_1^4\right)}^q{p}_0^4\left(x_{j-\frac{1}{2}}\right)+5{\left(ϵ+{\beta }_0^4\right)}^q{p}_1^4\left(x_{j-\frac{1}{2}}\right)}{3{\left(ϵ+{\beta }_1^4\right)}^q+5{\left(ϵ+{\beta }_0^4\right)}^q},$ (28)

and similarly inserting the weights (19) into (17) gives

$F_2^5\left(x_{j-\frac{1}{2}}\right)=\frac{5{\left(ϵ+{\beta }_2^4\right)}^q{p}_1^4\left(x_{j-\frac{1}{2}}\right)+3{\left(ϵ+{\beta }_1^4\right)}^q{p}_2^4\left(x_{j-\frac{1}{2}}\right)}{5{\left(ϵ+{\beta }_2^4\right)}^q+3{\left(ϵ+{\beta }_1^4\right)}^q}.$ (29)

Finally, substituting the 5th-order approximations (28) (29) into (22), we obtain

$F_2^6\left(x_{j-\frac{1}{2}}\right)=w_0{p}_0^4\left(x_{j-\frac{1}{2}}\right)+w_1{p}_1^4\left(x_{j-\frac{1}{2}}\right)+w_2{p}_2^4\left(x_{j-\frac{1}{2}}\right),$ (30)

Through reorganization, the nonlinear weights in (30) are

$w_0={\alpha }_0{\alpha }_2,w_1={\alpha }_2\left(1-{\alpha }_0\right)+{\alpha }_1\left(1-{\alpha }_2\right),w_2=\left(1-{\alpha }_1\right)\left(1-{\alpha }_2\right),$ (31)

where the unnormalized weights are

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}},\\ {\alpha }_1=\frac{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_2^4\right)}^q}},\\ {\alpha }_2=\frac{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2^4\right)}^q}}.\end{array}$ (32)

Theorem 1. Suppose that the stencil $S_2^6$ contains a discontinuity at most, the exponent in unnormalized weights $q\ge 1$ and $ϵ\le \mathcal{O}\left(\Delta x^4\right)$ . Then the 6th-order sub-WENO algorithm (30) (31) (32) satisfies the following three cases:

Case 1: if the stencil $S_2^6$ is smooth, then

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^6\right);$

Case 2: if there is a discontinuity lies in $\left(x_{j-3},x_{j-2}\right)$ , then

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^5\right),$

or if there is a discontinuity lies in $\left(x_{j+1},x_{j+2}\right)$ , then

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^5\right);$

Case 3: if there is a discontinuity lies in $\left(x_{j-2},x_{j-1}\right)$ , then

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^4\right),$

or if there is a discontinuity lies in $\left(x_j,x_{j+1}\right)$ , then

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\text{Δ}{x}^4\right).$

Proof 1. For the case 1, since the stencil $S_2^6$ is smooth, we obtain

$\begin{array}{l}{\beta }_i^4={\left(\text{Δ}{x}^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right),i=0,1,2.\hfill \end{array}$ (33)

By Taylor analysis, the unnormalized weights in (32) satisfy

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}=\frac{3{\left(ϵ+{\beta }_1^4\right)}^q}{3{\left(ϵ+{\beta }_1^4\right)}^q+5{\left(ϵ+{\beta }_0^4\right)}^q}=\frac{3}{8}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\alpha }_1=\frac{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{5{\left(ϵ+{\beta }_2^4\right)}^q}{5{\left(ϵ+{\beta }_2^4\right)}^q+3{\left(ϵ+{\beta }_1^4\right)}^q}=\frac{5}{8}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\alpha }_2=\frac{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{{\left(ϵ+{\beta }_2^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q+{\left(ϵ+{\beta }_0^4\right)}^q}=\frac{1}{2}+\mathcal{O}\left(\text{Δ}{x}^2\right).\end{array}$ (34)

Then the weights in (31) satisfy

$\begin{array}{l}{w}_0={\alpha }_0{\alpha }_2=\frac{3}{16}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ w_1={\alpha }_2\left(1-{\alpha }_0\right)+{\alpha }_1\left(1-{\alpha }_2\right)=\frac{10}{16}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ w_2=\left(1-{\alpha }_1\right)\left(1-{\alpha }_2\right)=\frac{3}{16}+\mathcal{O}\left(\text{Δ}{x}^2\right).\end{array}$ (35)

By using the similar discussion as (8) and denoting $\left(d_0,d_1,d_2\right)=\left(\frac{3}{16},\frac{10}{16},\frac{3}{16}\right)$ , we have

$\begin{array}{c}{F}_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum }_{i=0}^2{w}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)p_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)p_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)f_{j-\frac{1}{2}}+{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)\left(p_i^4\left(x_{j-\frac{1}{2}}\right)-f_{j-\frac{1}{2}}\right)+p_2^6\left(x_{j-\frac{1}{2}}\right)\\ =p_2^6\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\Delta x^6\right)\\ =f_{j-\frac{1}{2}}+\mathcal{O}\left(\Delta x^6\right).\end{array}$ (36)

For the case 2, we only discuss the situation in which the discontinuity lies in the interval $\left(x_{j-3},x_{j-2}\right)$ since $\left(x_{j+1},x_{j+2}\right)$ and $\left(x_{j-3},x_{j-2}\right)$ are symmetric. We first analyze the case of corner discontinuity which lies in $\left(x_{j-3},x_{j-2}\right)$ . At this moment three 4-points stencils satisfy

$\begin{array}{l}{\beta }_0^4=\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\beta }_1^4={\left(\text{Δ}{x}^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right),\\ {\beta }_2^4={\left(\text{Δ}{x}^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right).\end{array}$ (37)

The unnormalized weights in (32) satisfy

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}=\frac{3{\left(ϵ+{\beta }_1^4\right)}^q}{3{\left(ϵ+{\beta }_1^4\right)}^q+5{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{2q}\right),\\ {\alpha }_1=\frac{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{5{\left(ϵ+{\beta }_2^4\right)}^q}{5{\left(ϵ+{\beta }_2^4\right)}^q+3{\left(ϵ+{\beta }_1^4\right)}^q}=\frac{5}{8}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\alpha }_2=\frac{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{{\left(ϵ+{\beta }_2^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q+{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{2q}\right).\end{array}$ (38)

The weights in (31) satisfy

$\begin{array}{l}{w}_0={\alpha }_0{\alpha }_2=\mathcal{O}\left(\text{Δ}{x}^{4q}\right),\\ w_1={\alpha }_2\left(1-{\alpha }_0\right)+{\alpha }_1\left(1-{\alpha }_2\right)=\frac{5}{8}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ w_2=\left(1-{\alpha }_1\right)\left(1-{\alpha }_2\right)=\frac{3}{8}+\mathcal{O}\left(\text{Δ}{x}^2\right).\end{array}$ (39)

Denoting $\left(d_0,d_1,d_2\right)=\left(0,\frac{5}{8},\frac{3}{8}\right)$ and one can easily verify ${\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)=p_2^5\left(x_{j-\frac{1}{2}}\right)$ , then we have

$\begin{array}{c}{F}_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum }_{i=0}^2{w}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)p_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)p_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)f_{j-\frac{1}{2}}+{\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum }_{i=0}^2\left(w_i-d_i\right)\left(p_i^4\left(x_{j-\frac{1}{2}}\right)-f_{j-\frac{1}{2}}\right)+p_2^5\left(x_{j-\frac{1}{2}}\right)\\ =p_2^5\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\Delta x^6\right)\\ =f_{j-\frac{1}{2}}+\mathcal{O}\left(\Delta x^5\right).\end{array}$ (40)

If there is a jump discontinuity in the interval $\left(x_{j-3},x_{j-2}\right)$ , then the indicators of smoothness satisfy

$\begin{array}{l}{\beta }_0^4=\mathcal{O}\left(1\right),\\ {\beta }_1^4={\left(\Delta x^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\Delta x^2\right)\right),\\ {\beta }_2^4={\left(\Delta x^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\Delta x^2\right)\right).\end{array}$ (41)

Inserting these indicators of smoothness into the unnormalized weights in (32) gives

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}=\frac{3{\left(ϵ+{\beta }_1^4\right)}^q}{3{\left(ϵ+{\beta }_1^4\right)}^q+5{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{4q}\right),\\ {\alpha }_1=\frac{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{5{\left(ϵ+{\beta }_2^4\right)}^q}{5{\left(ϵ+{\beta }_2^4\right)}^q+3{\left(ϵ+{\beta }_1^4\right)}^q}=\frac{5}{8}+\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\alpha }_2=\frac{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{{\left(ϵ+{\beta }_2^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q+{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{4q}\right).\end{array}$ (42)

Then the nonlinear weights in (31) satisfy

$\begin{array}{l}{w}_0={\alpha }_0{\alpha }_2=\mathcal{O}\left(\Delta x^{8q}\right),\\ w_1={\alpha }_2\left(1-{\alpha }_0\right)+{\alpha }_1\left(1-{\alpha }_2\right)=\frac{5}{8}+\mathcal{O}\left(\Delta x^2\right),\\ w_2=\left(1-{\alpha }_1\right)\left(1-{\alpha }_2\right)=\frac{3}{8}+\mathcal{O}\left(\Delta x^2\right).\end{array}$ (43)

Repeating (40) shows the same result,

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f_{j-\frac{1}{2}}+\mathcal{O}\left(\text{Δ}{x}^5\right).$ (44)

Finally, we consider the case 3, in which the discontinuity lies in interval $\left(x_{j-2},x_{j-1}\right)$ or $\left(x_j,x_{j+1}\right)$ . Again we only analyze the situation of the discontinuity lies in interval $\left(x_{j-2},x_{j-1}\right)$ since $\left(x_j,x_{j+1}\right)$ and $\left(x_{j-2},x_{j-1}\right)$ are symmetric. we take into account the corner discontinuity which lies in $\left(x_{j-2},x_{j-1}\right)$ , then the indicators of smoothness achieve

$\begin{array}{l}{\beta }_0^4=\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\beta }_1^4=\mathcal{O}\left(\text{Δ}{x}^2\right),\\ {\beta }_2^4={\left(\text{Δ}{x}^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right).\end{array}$ (45)

Substituting the indicators of smoothness into the unnormalized weights in (32) shows

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}=\frac{3{\left(ϵ+{\beta }_1^4\right)}^q}{3{\left(ϵ+{\beta }_1^4\right)}^q+5{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(1\right),\\ {\alpha }_1=\frac{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{5{\left(ϵ+{\beta }_2^4\right)}^q}{5{\left(ϵ+{\beta }_2^4\right)}^q+3{\left(ϵ+{\beta }_1^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{2q}\right),\\ {\alpha }_2=\frac{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{{\left(ϵ+{\beta }_2^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q+{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{2q}\right).\end{array}$ (46)

Then the nonlinear weights in (31) satisfy

$\begin{array}{l}{w}_0={\alpha }_0{\alpha }_2=\mathcal{O}\left(\text{Δ}{x}^{2q}\right),\\ w_1={\alpha }_2\left(1-{\alpha }_0\right)+{\alpha }_1\left(1-{\alpha }_2\right)=\mathcal{O}\left(\text{Δ}{x}^{2q}\right),\\ w_2=\left(1-{\alpha }_1\right)\left(1-{\alpha }_2\right)=1+\mathcal{O}\left(\text{Δ}{x}^{2q}\right).\end{array}$ (47)

Denoting $\left(d_0,d_1,d_2\right)=\left(0,0,1\right)$ and one can easily verify ${\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)=p_2^4\left(x_{j-\frac{1}{2}}\right)$ , then we have

$\begin{array}{c}{F}_2^6\left(x_{j-\frac{1}{2}}\right)={\displaystyle \sum _{i=0}^2}{w}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum _{i=0}^2}{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum _{i=0}^2}{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum _{i=0}^2}\left(w_i-d_i\right)p_i^4\left(x_{j-\frac{1}{2}}\right)+{\displaystyle \sum _{i=0}^2}{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum _{i=0}^2}\left(w_i-d_i\right)p_i^4\left(x_{j-\frac{1}{2}}\right)-{\displaystyle \sum _{i=0}^2}\left(w_i-d_i\right)f_{j-\frac{1}{2}}+{\displaystyle \sum _{i=0}^2}{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)\\ ={\displaystyle \sum _{i=0}^2}\left(w_i-d_i\right)\left(p_i^4\left(x_{j-\frac{1}{2}}\right)-f_{j-\frac{1}{2}}\right)+p_2^4\left(x_{j-\frac{1}{2}}\right)\\ =p_2^4\left(x_{j-\frac{1}{2}}\right)+\mathcal{O}\left(\Delta x^{2q+4}\right)\\ =f_{j-\frac{1}{2}}+\mathcal{O}\left(\Delta x^4\right).\end{array}$ (48)

If the jump discontinuity lies in $\left(x_{j-2},x_{j-1}\right)$ , then the indicators of smoothness approach

$\begin{array}{l}{\beta }_0^4=\mathcal{O}\left(1\right),\\ {\beta }_1^4=\mathcal{O}\left(1\right),\\ {\beta }_2^4={\left(\text{Δ}{x}^2{f^{″}}_{j-\frac{1}{2}}\right)}^2\left(1+\mathcal{O}\left(\text{Δ}{x}^2\right)\right).\end{array}$ (49)

Replacing the indicators of smoothness into the unnormalized weights in (32) shows

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}=\frac{3{\left(ϵ+{\beta }_1^4\right)}^q}{3{\left(ϵ+{\beta }_1^4\right)}^q+5{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(1\right),\\ {\alpha }_1=\frac{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_1^4\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{5{\left(ϵ+{\beta }_2^4\right)}^q}{5{\left(ϵ+{\beta }_2^4\right)}^q+3{\left(ϵ+{\beta }_1^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{4q}\right),\\ {\alpha }_2=\frac{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0^4\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2^4\right)}^q}}=\frac{{\left(ϵ+{\beta }_2^4\right)}^q}{{\left(ϵ+{\beta }_2^4\right)}^q+{\left(ϵ+{\beta }_0^4\right)}^q}=\mathcal{O}\left(\text{Δ}{x}^{4q}\right).\end{array}$ (50)

Then the nonlinear weights in (31) satisfy

$\begin{array}{l}{w}_0={\alpha }_0{\alpha }_2=\mathcal{O}\left(\Delta x^{4q}\right),\\ w_1={\alpha }_2\left(1-{\alpha }_0\right)+{\alpha }_1\left(1-{\alpha }_2\right)=\mathcal{O}\left(\Delta x^{4q}\right),\\ w_2=\left(1-{\alpha }_1\right)\left(1-{\alpha }_2\right)=1+\mathcal{O}\left(\Delta x^{4q}\right).\end{array}$ (51)

Denoting $\left(d_0,d_1,d_2\right)=\left(0,0,1\right)$ and one can easily verify ${\displaystyle \sum }_{i=0}^2{d}_i{p}_i^4\left(x_{j-\frac{1}{2}}\right)=p_2^4\left(x_{j-\frac{1}{2}}\right)$ , then we have

$F_2^6\left(x_{j-\frac{1}{2}}\right)=f_{j-\frac{1}{2}}+\mathcal{O}\left(\text{Δ}{x}^4\right).$ (52)

Theorem 1 shows the maximum theoretical accuracy of the 6th-order sub-WENO algorithm near discontinuities. What we left is to ensure the ENO property of new algorithm. The algorithm is deemed to conserve the ENO property if it satisfies the following two conditions,

1) If the stencil $S_i^m$ is smooth, then the nonlinear weight corresponding $S_i^m$ satisfy $w_i^m=\mathcal{O}\left(1\right)$ .

2) If the stencil $S_i^m$ is non-smooth, then the nonlinear weight corresponding $S_i^m$ satisfies $w_i^m\le \mathcal{O}\left(\text{Δ}{x}^m\right)$ .

Clearly, from the nonlinear weights (39) (43) (47) and (51) in theorem 1, we can validate the ENO property of new algorithm if $q\ge 2$ .

4. The Higher Order Sub-WENO Alrorithm

In this section, the 2rth-order sub-WENO algorithms are presented for $r\ge 2$ . In particular, when $r=2$ the 4th-order sub-WENO algorithm is the same as the classical WENO algorithm. The 6th-order sub-WENO algorithm ($r=3$ ) has been shown in Section 3 and we also analyze the optimal accuracy and ENO property. To implement the algorithm clearly, we give the tree structure of 2rth-order sub-WENO algorithm in Section 4 (Figure 2).

Figure 2. The tree structure of 2rth-order sub-WENO algorithm.

In the first level, there are $r$ $\left(r+1\right)$ -points stencils and corresponding indicators of smoothness ${\beta }_i^{r+1}\left(i=0,\cdots ,r-1\right)$ . As the 6th-order sub-WENO algorithm in the Section 3, it will be found that we only need to compute the indicators of smoothness ${\beta }_i^{r+1}\left(i=0,\cdots ,r-1\right)$ of the smallest stencils on the first level. The choice of indicator of smoothness ${\beta }_i^m\left(r+2\le m\le 2r-1,m-r-1\le i\le r-1\right)$ should ensure the optimal accuracy, ENO property and modest computational cost of sub-WENO algorithm. From tree structure, we can find that each nonlinear approximation $F_i^m\left(x_{j-\frac{1}{2}}\right)$ on stencil $S_i^m$ is obtained by the combination of the nonlinear approximations $F_{i-1}^{m-1}\left(x_{j-\frac{1}{2}}\right)$ and $F_i^{m-1}\left(x_{j-\frac{1}{2}}\right)$ on $S_{i-1}^{m-1}$ and $S_i^{m-1}$ , respectively,

$F_i^m\left(x_{j-\frac{1}{2}}\right)=w_{i-1}^{m-1}{F}_{i-1}^{m-1}\left(x_{j-\frac{1}{2}}\right)+w_i^{m-1}{F}_i^{m-1}\left(x_{j-\frac{1}{2}}\right).$ (53)

That is, in each sub-WENO procedure $\left\{S_{i-1}^{m-1},S_i^{m-1}\right\}$ ↪$S_i^m$ , we only use the information from two stencils. The nonlinear weights are defined as

$\begin{array}{l}{w}_{i-1}^{m-1}=\frac{\frac{d_{i-1}^{m-1}}{{\left(ϵ+{\beta }_{i-1}^{m-1}\right)}^q}}{\frac{d_{i-1}^{m-1}}{{\left(ϵ+{\beta }_{i-1}^{m-1}\right)}^q}+\frac{d_i^{m-1}}{{\left(ϵ+{\beta }_i^{m-1}\right)}^q}},\\ w_i^{m-1}=\frac{\frac{d_i^{m-1}}{{\left(ϵ+{\beta }_i^{m-1}\right)}^q}}{\frac{d_{i-1}^{m-1}}{{\left(ϵ+{\beta }_{i-1}^{m-1}\right)}^q}+\frac{d_i^{m-1}}{{\left(ϵ+{\beta }_i^{m-1}\right)}^q}}.\end{array}$ (54)

The linear optimal weights $d_{i-1}^{m-1}$ and $d_i^{m-1}$ in (54) are chosen so that they satisfy

$p_i^m\left(x_{j-\frac{1}{2}}\right)=d_{i-1}^{m-1}{p}_{i-1}^{m-1}\left(x_{j-\frac{1}{2}}\right)+d_i^{m-1}{p}_i^{m-1}\left(x_{j-\frac{1}{2}}\right).$

Since $S_i^m=S_{r+1-m+i}^{r+1}\cup \cdots \cup S_i^{r+1}$ and, as we have done for the 6th-order sub-WENO algorithm, a reasonable choice is

$\begin{array}{l}{\beta }_{i-1}^{m-1}:={\beta }_{r+1-m+i}^{r+1}\cdots {\beta }_{i-1}^{r+1},\\ {\beta }_i^{m-1}:={\beta }_{r+2-m+i}^{r+1}\cdots {\beta }_i^{r+1}.\end{array}$ (55)

Replacing $ϵ+{\beta }_{i-1}^{m-1}$ and $ϵ+{\beta }_i^{m-1}$ in (54) by $\left(ϵ+{\beta }_{r+1-m+i}^{r+1}\right)\cdots \left(ϵ+{\beta }_{i-1}^{r+1}\right)$ and $\left(ϵ+{\beta }_{r+2-m+i}^{r+1}\right)\cdots \left(ϵ+{\beta }_i^{r+1}\right)$ , respectively, and canceling the common terms, then (54) is simplified as

$\begin{array}{l}{w}_{i-1}^{m-1}=\frac{\frac{d_{i-1}^{m-1}}{{\left(ϵ+{\beta }_{r+1-m+i}^{r+1}\right)}^q}}{\frac{d_{i-1}^{m-1}}{{\left(ϵ+{\beta }_{r+1-m+i}^{r+1}\right)}^q}+\frac{d_i^{m-1}}{{\left(ϵ+{\beta }_i^{r+1}\right)}^q}},\\ w_i^{m-1}=\frac{\frac{d_i^{m-1}}{{\left(ϵ+{\beta }_i^{r+1}\right)}^q}}{\frac{d_{i-1}^{m-1}}{{\left(ϵ+{\beta }_{r+1-m+i}^{r+1}\right)}^q}+\frac{d_i^{m-1}}{{\left(ϵ+{\beta }_i^{r+1}\right)}^q}}.\end{array}$ (56)

Section 4 helps us to determine the indicators of smoothness of stencils intuitively. The indicators of smoothness of stencils $S_{r+1-m+i}^{r+1}$ and $S_i^{r+1}$ on the top of “V” corresponding the sub-WENO procedure $\left\{S_{i-1}^{m-1},S_i^{m-1}\right\}$ ↪$S_i^m$ are chosen to be the indicators of smoothness of $S_{i-1}^{m-1}$ and $S_i^{m-1}$ , respectively. It is remarkable that the stencil $S_i^m$ is involved simultaneously in both the sub-WENO procedures $\left\{S_{i-1}^{m-1},S_i^{m-1}\right\}$ ↪$S_i^m$ and $\left\{S_i^{m-1},S_{i+1}^{m-1}\right\}$ ↪$S_{i+1}^m$ . But the indicator of smoothness of $S_i^m$ is different when it belongs to different sub-WENO procedure (Figure 3).

Figure 3. The choice of indicator of smoothness for sub-WENO procedure.

To validate such choice of indicators of smoothness, taking the 8th-order sub-WENO algorithm for example, Section 4 shows the cases in which the discontinuity lies in $\left(x_{j+2},x_{j+3}\right)$ , $\left(x_{j+1},x_{j+2}\right)$ and $\left(x_j,x_{j+1}\right)$ , respectively. The left tree structure of Section 4 shows the situation when the discontinuity lies in the interval $\left(x_{j+2},x_{j+3}\right)$ , i.e., the former three stencils $S_i^5\left(i=0,1,2\right)$ are smooth and the last one $S_3^5$ is non-smooth. From the first level to second level, there are three sub-WENO procedures. The former two sub-WENO procedures both produce the 6th-order approximations and the third sub-WENO procedure only produce a 5th-order approximation since the stencil $S_3^5$ is non-smooth and the nonlinear weight assigned to it can be ignored. From the second level to third level, there are two sub-WENO procedures. In the first sub-WENO procedure $\left\{S_1^6,S_2^6\right\}$ ↪$S_2^7$ , since the indicators of smoothness of stencils $S_1^6$ and $S_2^6$ are defined by

${\beta }_1^6:={\beta }_0^5 and {\beta }_2^6:={\beta }_2^5,$

this sub-WENO procedure generates 7th-order approximation. The second sub-WENO procedure $\left\{S_2^6,S_3^6\right\}$ ↪$S_3^7$ only produces the 6th-order approximation, since ${\beta }_3^6:={\beta }_3^5$ is non-smooth. In the final sub-WENO procedure $\left\{S_2^7,S_3^7\right\}$ ↪$S_3^8$ , the indicators of smoothness are defined by

${\beta }_2^7:={\beta }_0^5\text{and}{\beta }_3^7:={\beta }_3^5.$

Since $S_3^5$ is non-smooth, the contribution of this stencil is ignored and the final approximation $F_3^8\left(x_{j-\frac{1}{2}}\right)$ is only 7th-order. The middle and right tree structures of Section 4 show the situations when the discontinuities lie in $\left(x_{j+1},x_{j+2}\right)$ and $\left(x_j,x_{j+1}\right)$ , respectively. The “$\ast$ ” denotes the situation that the value generated by this sub-WENO procedure is nonsense since both substencils are non-smooth. By the strategy of choosing the indicators of smoothness, the value denoted by “$\ast$ ” will be ignored almost in the latter sub-WENO procedures.

To implement the higher order sub-WENO algorithm easily, we present the necessary formulas explicitly for $r=4$ . Since all the sub-WENO procedures can be integrated, we will only present the compact forms of sub-WENO algorithms as (30) (31) (32). For the 8th-order sub-WENO algorithm, there exist four stencils $S_i^5\left(i=0,1,2,3\right)$ used to approximate $f_{i-\frac{1}{2}}$ . The linear 5th-order interpolations at $x_{j-\frac{1}{2}}$ can be expressed as

$\begin{array}{l}{p}_0^5\left(x_{j-\frac{1}{2}}\right)=\frac{1}{128}\left(-5f_{j-4}+28f_{j-3}-70f_{j-2}+140f_{j-1}+35f_j\right),\\ p_1^5\left(x_{j-\frac{1}{2}}\right)=\frac{1}{128}\left(3f_{j-3}-20f_{j-2}+90f_{j-1}+60f_j-5f_{j+1}\right),\\ p_2^5\left(x_{j-\frac{1}{2}}\right)=\frac{1}{128}\left(-5f_{j-2}+60f_{j-1}+90f_j-20f_{j+1}+3f_{j+2}\right),\\ p_3^5\left(x_{j-\frac{1}{2}}\right)=\frac{1}{128}\left(35f_{j-1}+140f_j-70f_{j+1}+28f_{j+2}-5f_{j+3}\right).\end{array}$ (57)

By combining all the sub-WENO procedures (Figure 4), we can arrive

Figure 4. The choice of indicator of smoothness for sub-WENO procedure keeps the ENO property. The circles “$○$ ” denotes the non-smooth stencils and the stars “$\ast$ ” denotes the values generated by corresponding sub-WENO procedures are nonsense.

$F_3^8\left(x_{j-\frac{1}{2}}\right)=w_0{p}_0^5\left(x_{j-\frac{1}{2}}\right)+w_1{p}_1^5\left(x_{j-\frac{1}{2}}\right)+w_2{p}_2^5\left(x_{j-\frac{1}{2}}\right)+w_3{p}_3^5\left(x_{j-\frac{1}{2}}\right),$ (58)

where

$\begin{array}{l}{w}_0={\alpha }_0{\alpha }_3{\alpha }_5,w_1={\beta }_0{\alpha }_3{\alpha }_5+{\alpha }_1{\beta }_3{\alpha }_5+{\alpha }_1{\alpha }_4{\beta }_5,\\ w_2={\beta }_1{\beta }_3{\alpha }_5+{\beta }_1{\alpha }_4{\beta }_5+{\alpha }_2{\beta }_4{\beta }_5,w_3={\beta }_2{\beta }_4{\beta }_5,\end{array}$ (59)

and ${\beta }_i=1-{\alpha }_i$ , $i=0,\cdots ,5$ . The unnormalized weights are defined as

$\begin{array}{l}{\alpha }_0=\frac{\frac{3}{{\left(ϵ+{\beta }_0\right)}^q}}{\frac{3}{{\left(ϵ+{\beta }_0\right)}^q}+\frac{7}{{\left(ϵ+{\beta }_1\right)}^q}},{\alpha }_1=\frac{\frac{1}{{\left(ϵ+{\beta }_1\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_1\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_2\right)}^q}},\\ {\alpha }_2=\frac{\frac{7}{{\left(ϵ+{\beta }_2\right)}^q}}{\frac{7}{{\left(ϵ+{\beta }_2\right)}^q}+\frac{3}{{\left(ϵ+{\beta }_3\right)}^q}},{\alpha }_3=\frac{\frac{5}{{\left(ϵ+{\beta }_0\right)}^q}}{\frac{5}{{\left(ϵ+{\beta }_0\right)}^q}+\frac{7}{{\left(ϵ+{\beta }_2\right)}^q}},\\ {\alpha }_4=\frac{\frac{7}{{\left(ϵ+{\beta }_1\right)}^q}}{\frac{7}{{\left(ϵ+{\beta }_1\right)}^q}+\frac{5}{{\left(ϵ+{\beta }_3\right)}^q}},{\alpha }_5=\frac{\frac{1}{{\left(ϵ+{\beta }_0\right)}^q}}{\frac{1}{{\left(ϵ+{\beta }_0\right)}^q}+\frac{1}{{\left(ϵ+{\beta }_3\right)}^q}}.\end{array}$ (60)

The indicators of smoothness computed by (2) can be expressed as

$\begin{array}{c}{\beta }_0^5=\frac{1}{36}{\left(2f_{j-4}-11f_{j-3}+27f_{j-2}-29f_{j-1}+11f_j\right)}^2\\ +\frac{39}{36}{\left(f_{j-4}-5f_{j-3}+9f_{j-2}-7f_{j-1}+2f_j\right)}^2\\ +\frac{3124}{2880}{\left(f_{j-4}-4f_{j-3}+6f_{j-2}-4f_{j-1}+f_j\right)}^2,\\ {\beta }_1^5=\frac{1}{36}{\left(f_{j-3}-7f_{j-2}+9f_{j-1}-f_j-2f_{j+1}\right)}^2\\ +\frac{39}{36}{\left(-f_{j-2}+3f_{j-1}-3f_j+f_{j+1}\right)}^2\end{array}$

$\begin{array}{c}+\frac{3124}{2880}{\left(f_{j-3}-4f_{j-2}+6f_{j-1}-4f_j+f_{j+1}\right)}^2,\\ {\beta }_2^5=\frac{1}{36}{\left(2f_{j-2}+f_{j-1}-9f_j+7f_{j+1}-f_{j+2}\right)}^2\\ +\frac{39}{36}{\left(f_{j-2}-3f_{j-1}+3f_j-f_{j+1}\right)}^2\\ +\frac{3124}{2880}{\left(f_{j-2}-4f_{j-1}+6f_j-4f_{j+1}+f_{j+2}\right)}^2,\\ {\beta }_3^5=\frac{1}{36}{\left(11f_{j-1}-29f_j+27f_{j+1}-11f_{j+2}+2f_{j+3}\right)}^2\\ +\frac{39}{36}{\left(2f_{j-1}-7f_j+9f_{j+1}-5f_{j+2}+f_{j+3}\right)}^2\\ +\frac{3124}{2880}{\left(f_{j-1}-4f_j+6f_{j+1}-4f_{j+2}+f_{j+3}\right)}^2.\end{array}$ (61)

5. Numerical Results

In this section, two examples, containing the corner and jump discontinuities respectively, are used to test the order of accuracy of sub-WENO algorithms. We also compare the computational costs between the classical WENO [7], the WENO-AW [18] and sub-WENO algorithms. The 8th-order sub-WENO algorithms is only tested by the first example since we can obtain the similar results for the second example. We choose the parameters $ϵ={10}^{-40}$ and $q=2$ for all the algorithms.

Example 1. Consider the function

$f\left(x\right)=\{\begin{array}{ll}-10{\text{e}}^{-x+3}-3x^2,\hfill & \text{if}-1\le x\le 0,\hfill \\ -10{\text{e}}^{x+3}-3x^2,\hfill & \text{if}0<x\le 1,\hfill \end{array}$ (62)

there is a corner discontinuity lies at $x=0$ . This piecewise function is even and continuous, but the first derivative is not discontinuous. Since the function is symmetrical with $y$ -axis, we only show the absolute errors at $\left\{x_{j-1/2},x_{j+1/2},x_{j+3/2},x_{j+5/2}\right\}$ . While it is only interested in the order of accuracy near discontinuities in this paper, we keep the discontinuities lie in the center of grid $\left[x_{j-1},x_j\right]$ . To get this, the interval $\left[-1,1\right]$ is divided into $2N+1$ grids and then we subdivide each grid into three small uniform grids in every grid refinement. In the numerical experiments, the grid spacing is set to be $\text{Δ}x=\frac{2}{101\cdot 3^i}$ , $i=0,\cdots ,4$ . Table 1 shows the absolute errors and order of accuracy of classical WENO, WENO-AW and sub-WENO algorithms. As presented in Section 2.1, when the discontinuity lies in the interval $\left[x_{j-1},x_j\right]$ , the approximations of classical WENO algorithm at $x_{j+1/2}$ and $x_{j+3/2}$ can only arrive the 4th-order of accuracy. When we predict the value at $x_{j+5/2}$ , since the stencil $S_2^6$ for approximating $x_{j+5/2}$ is smooth, the maximum 6th-order of accuracy is certainly obtained. As shown in [18] and Section 3, the WENO-AW and sub-WENO algorithms are devised to reasonably combine all the smooth stencils to arrive the maximum theoretical accuracy. The WENO-AW and sub-WENO algorithms achieves 5th-order of accuracy when approximate the value at $x_{j+3/2}$ since the 5-points stencil $S_1^5$ is smooth. It is noted that the absolute errors presented by WENO-AW and sub-WENO algorithms are almost identical.

Table 1. The absolute errors and orders of accuracy of the 6th-order classical WENO, WENO-AW and sub-WENO algorithms for the function (62) which contains a corner discontinuity at $x=0$ .

classical WENO
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$
0	1.478e−00 *****	1.080e−06 *****	2.780e−07 *****	4.284e−10 *****
1	4.929e−01 1.000	1.438e−08 3.931	3.439e−09 3.998	5.649e−13 6.034
2	1.643e−01 1.000	1.818e−10 3.979	4.245e−11 4.000	7.646e−16 6.012
3	5.477e−02 1.000	2.261e−12 3.993	5.239e−13 4.000	1.044e−18 6.004
4	1.826e−02 1.000	2.799e−14 3.998	6.467e−15 4.000	1.430e−21 6.001
WENO-AW
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$
0	1.478e−00 *****	1.109e−06 *****	7.869e−09 *****	4.284e−10 *****
1	4.929e−01 1.000	1.450e−08 3.917	3.044e−11 5.056	5.649e−13 6.034
2	1.643e−01 1.000	1.822e−10 3.984	1.227e−13 5.019	7.646e−16 6.012
3	5.477e−02 1.000	2.263e−12 3.995	5.012e−16 5.006	1.044e−18 6.004
4	1.826e−02 1.000	2.800e−14 3.998	2.058e−18 5.002	1.430e−21 6.001
sub-WENO
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$
0	1.478e−00 *****	1.144e−06 *****	7.869e−09 *****	4.283e−10 *****
1	4.929e−01 1.000	1.464e−08 3.967	3.044e−11 5.056	5.649e−13 6.036
2	1.643e−01 1.000	1.828e−10 3.989	1.227e−13 5.019	7.646e−16 6.012
3	5.477e−02 1.000	2.266e−12 3.996	5.012e−16 5.006	1.044e−18 6.004
4	1.826e−02 1.000	2.801e−14 3.999	2.058e−18 5.002	1.430e−21 6.001

Table 2 compares the computational costs of three WENO algorithms. To obtain the reliable CPU costs, we loop the interpolation parts of codes 5 million times. The results are shown in Table 2 and we find that the CPU times of the WENO-AW and sub-WENO algorithms increase approximately by 150% and 5% respectively when compared with the classical WENO algorithm.

Table 2. The CPU costs of the 6th-order classical WENO, WENO-AW and sub-WENO algorithms for the function (62). The interpolation parts of codes are looped 5 million times to obtain the reliable CPU costs.

classical WENO	WENO-AW	sub-WENO
13.8s	35.4s	14.6s

The performance of 8th-order sub-WENO algorithm is presented in Table 3. We can find that the algorithm achieves the optimal order near the corner discontinuity. When we approximate the value at $x_{j+1/2}$ , there exists only one stencil $S_3^5$ is smooth. The nonlinear weight distributed to this stencil is dominant and the contributions of the left stencils are ignored. The prediction at $x_{j+3/2}$ by sub-WENO algorithm gives a reasonable combination of $S_2^5$ and $S_3^5$ since they are both smooth. When approximate the value at $x_{j+5/2}$ , the discontinuity lies in the far left interval $\left[x_{j-1},x_j\right]$ . The stencils $S_1^5$ , $S_2^5$ and $S_3^5$ all are smooth and the sub-WENO algorithm uses them well to obtain 7th-order of accuracy.

Table 3. The absolute errors and orders of accuracy of the 8th-order sub-WENO algorithms for the function (62) which contains a corner discontinuity at $x=0$ .

8th-order sub-WENO
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$	$x_{j+7/2}$
0	1.398e−00 *****	4.656e−08 *****	8.755e−11 *****	5.936e−13 *****	2.865e−14 *****
1	4.659e−01 1.000	1.847e−10 5.033	1.157e−13 6.034	2.690e−16 7.008	4.142e−18 8.048
2	1.553e−01 1.000	7.511e−13 5.011	1.567e−16 6.011	1.226e−19 7.003	6.200e−22 8.016
3	5.177e−02 1.000	3.079e−15 5.004	2.141e−19 6.004	5.600e−23 7.001	9.392e−26 8.006
4	1.726e−02 1.000	1.265e−17 5.001	2.933e−22 6.001	2.560e−26 7.000	1.425e−29 8.004

Example 2. In this example, we slightly modify the piecewise function (62) so that the corner discontinuity in it is changed into jump discontinuity,

$f\left(x\right)=\{\begin{array}{ll}10{\text{e}}^{x+3}+3x^2,\hfill & \text{if}-1\le x\le 0,\hfill \\ -10{\text{e}}^{x+3}-3x^2,\hfill & \text{if}0<x\le 1.\hfill \end{array}$ (63)

We only test the 6th-order WENO algorithms for this example. As shown in Table 4, the classical WENO algorithm only arrives 4th-order of accuracy when we approximate the values at $x_{j+3/2}$ which is near discontinuity. For the WENO-AW and sub-WENO algorithms, they are both obtain the maximum theoretical accuracy. It is noted that, when approximating the value at discontinuity, we will lose the accuracy entirely.

Table 4. The absolute errors and orders of accuracy of classical WENO, WENO-AW and sub-WENO algorithms for the function (63) which contains a jump discontinuity at $x=0$ .

classical WENO
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$
0	197.988 *****	1.250e−06 *****	2.847e−07 *****	4.284e−10 *****
1	199.899−0.009	1.507e−08 4.021	3.464e−09 4.013	5.649e−13 6.034
2	200.537−0.003	1.846e−10 4.007	4.254e−11 4.004	7.646e−16 6.012
3	200.749−0.001	2.273e−12 4.002	5.243e−13 4.002	1.044e−18 6.004
4	200.820−0.000	2.804e−14 4.001	6.469e−15 4.000	1.430e−21 6.001
WENO-AW
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$
0	197.605 *****	1.250e−06 *****	7.869e−09 *****	4.284e−10 *****
1	199.772−0.010	1.507e−08 4.021	3.044e−11 5.056	5.649e−13 6.034
2	200.494−0.003	1.846e−10 4.007	1.227e−13 5.019	7.646e−16 6.012
3	200.735−0.001	2.273e−12 4.002	5.012e−16 5.006	1.044e−18 6.004
4	200.815−0.000	2.804e−14 4.001	2.058e−18 5.002	1.430e−21 6.001
sub-WENO
$i$	$x_{j-1/2}$	$x_{j+1/2}$	$x_{j+3/2}$	$x_{j+5/2}$
0	197.243 *****	1.250e−06 *****	7.869e−09 *****	4.283e−10 *****
1	199.651−0.011	1.507e−08 4.021	3.044e−11 5.056	5.649e−13 6.036
2	200.454−0.004	1.846e−10 4.007	1.227e−13 5.019	7.646e−16 6.012
3	200.722−0.001	2.273e−12 4.002	5.012e−16 5.006	1.044e−18 6.004
4	200.811−0.000	2.804e−14 4.001	2.058e−18 5.002	1.430e−21 6.001

6. Conclusion

In this paper, the sub-WENO algorithm is presented to recover the optimal order of accuracy near the discontinuities. The sub-WENO algorithm is constructed by dividing the classical WENO into several sub-WENO procedures. In each sub-WENO procedure, we only combine two stencils to approximate the value of target points. If the two stencils are both smooth, then sub-WENO procedure increases the order of accuracy by one. If there is a stencil that is smooth and the left one is non-smooth, then algorithm conserves the order of interpolation by corresponding smooth stencil and keeps the ENO property. If both stencils are non-smooth, then the value constructed by sub-WENO procedure will be cut off in the latter procedures. The whole of sub-WENO algorithm can be expressed as tree structure. The choice of smoothness of indicator of stencils in the middle part of tree is also presented. This choice does not increase the computational time much. The proof of order of accuracy and ENO property of the 6th-order sub-WENO algorithm is shown in theorem 1. The numerical tests validate the results we have proved in theorem 1. The sub-WENO algorithm based on the cell averages and its application to hyperbolic conservation laws will be our future work.

Acknowledgements

This work was supported by Natural Science Foundation of Guangdong Province under contract no. 2018A030313879, the Technology Foundation of Huizhou under contract no. 2017C0403019 and Huizhou University under contract no. HZUX1201620.