On the approximation of derivative values using a WENO algorithm with progressive order of accuracy close to discontinuities

In this article, we introduce a new WENO algorithm that aims to calculate an approximation to derivative values of a function in a non-regular grid. We adapt the ideas presented in [Amat et al., SIAM J. Numer. Anal. (2020)] to design the nonlinear weights in a manner such that the order of accuracy is maximum in the intervals close to the discontinuities. Some proofs, remarks on the choice of the stencils and explicit formulas for the weights and smoothness indicators are given. We also present some numerical experiments to confirm the theoretical results.

employed by Levy et al. in Levy et al. (1999), which is designed to approximate solutions of hyperbolic systems of conservation laws, In our case, the motivation would be the possible application of the new algorithm to obtain approximations of the solution of Hamilton-Jacobi type equations, where algorithms that use interpolation to get accurate approximations to the derivative ϕ x have been proved to be very useful. In what follows, we show the construction of the algorithm and provide theoretical and numerical proofs of its accuracy, but we let the solution of the equation in (1) to be future work. Even so, as we show in what follows, the algorithm proposed can be used as a general technique to obtain the derivatives of a function with singularities, providing progressive order of accuracy close to these singularities. In Levy et al. (1999), the authors use a uniform discretization in time but not in space, i.e. x i = h i are the spatial grid points and t n = n t are the time steps. In Levy et al. (1999), two reconstructions are necessary to design the method: firstly, a reconstruction scheme to approximate the average of the function u over a cell at a time t n ; secondly, a reconstruction of f for the derivative of the fluxes, Q i , which satisfies: where r − 1 is the accuracy attained for the reconstruction of the derivative, and supposing that u at a time t n is known. WENO nonlinear interpolation methods can be used for this purpose (see e.g. Levy et al. (1999); Liu et al. (1994); Shu (1999)). In Levy et al. (1999), a symmetric strategy is used, but in this work, we construct a general method that is not necessarily symmetric. Therefore, in this paper, we will construct a new nonlinear method to approximate the derivative point values with progressive order of accuracy close to the discontinuities adapting the ideas presented in Amat et al. (2020). We will present this method for any r . We reformulate this problem as follows: We suppose [a, b] ⊂ R and a grid X = {x j } J j=0 , x 0 = a, x J = b, x j+1 = x j + h j , x j+1/2 = x j + h j /2, x j−1/2 = x j − h j−1 /2, h = max{h j : j = 0, . . . , J − 1}, and we suppose that there exists L > 0 such that For a piecewise smooth function f , we consider a point value discretization f (x j ) = f j , j = 0, . . . , J . To calculate the derivative value at a specific point x i , we take the stencils of r points which contain it, i.e. let I r i,k = {i − (r − 1) + k, . . . , i + k} be the subindex and the stencils S r i,k = {x t : t ∈ I i k } with k = 0, . . . , r −1 and calculate the interpolation polynomial of degree r − 1 such that p r −1 i,k (x j ) = f j , j ∈ I r i,k . The classic WENO interpolator is a convex combination of these polynomials: where ω r −1 i,k are nonlinear weights whose value depends on the existence of a discontinuity affecting the stencil. In our case, we will define It is clear that if the function is smooth at the largest stencil {x i−(r −1) , . . . , x i+(r −1) }, then, if the choice of the nonlinear weights is adequate, we can obtain an approximation of the derivative values of order of accuracy O(h 2r −2 ). However, when a discontinuity crosses the above-mentioned stencil, typically (see, e.g. Liu et al. (1994)) , the order decreases to r − 1.
Note that there are two principal differences between classical WENO interpolation for pointvalue approximation of a function and the WENO algorithm used to obtain an approximation of the derivative values: 1. When we use classical WENO algorithm for point-value interpolation, we take the stencils which contain x i−1 and x i because, typically, we interpolate at x i− 1 2 . 2. The evaluation of the derivative of the polynomial, P i , Eq. (3), is at x i not at x i− 1 2 .
The key of this method is the definition of the nonlinear weights that make use of linear optimal weights and that are defined as follows. Let p 2r −2 i,0 be the polynomial which interpolates f at {x i−(r −1) , . . . , x i+(r −1) }, then there exist C r −1 i,k > 0, k = 0, . . . , r − 1 which satisfy: The following proposition for non-uniform grid is proved in Yáñez (2010) using the ideas proposed in Aràndiga et al. (2010). We reproduce it here for completeness.

Proposition 1.1
The optimal weights are determined by: with 0 ≤ k < r − 1.
Proof Using Lagrange's base, we know that: Therefore: Analogously: We can obtain the optimal weights from the equality: thus, if k = 0, then the unique stencil of r points which contains the point x i−(r −1) is S r i,0 . Then, taking the term of f i−(r −1) , we obtain: Now, if we take k = 1, we have only two stencils: S r i,0 and S r i,1 , which contain the point x i−(r −1)+1 . Then using Eq. (6) at term f i−(r −1)+1 , we get: We suppose that we have obtained the optimal weights C r −1 then  Table 1 Optimal weights for r = 2, 3 and a uniform grid for the term f i−(r −1)+k in Eq. (6) we get: Then, In Tables 1 and 2, we show the values C r −1 i,k , 0 ≤ k ≤ r − 1 for r = 2, 3 and 4. If we take a uniform grid, the following corollary is a direct consequence by Prop. 1.1. We denote as C r −1 i,k = C r −1 k because this value does not depend on the position i. Corollary 1.2 Using a uniform grid, the optimal weights for the approximation of derivative values are given by: Proof From x i = ih, using the Lagrange basis, we have that Analogously,

Table 2
Optimal weights for r = 4 and a non-uniform grid By induction on k, if k = 0 and ξ = s − i, then We suppose the result for k − 1 and calculate C r −1 k , with 0 ≤ t < k ≤ r − 1. Thus, since we obtain: Then, since by induction hypothesis, (9) and Prop. 1.1, we have: In Table 3, some optimal weights are shown. For r = 3, they are the optimal weights shown in Levy et al. (1999).
With these linear weights, we can use the following expressions for the nonlinear weights: with r −1 k=0 ω r −1 i,k = 1. In Eq. (10), the parameter θ is an integer that assures maximum order of accuracy close to the discontinuities. In our case, we will take θ ≥ r and introduce the parameter > 0; to avoid divisions by zero, we will set it to = 10 −16 . The values I r −1 i,k are called smoothness indicators for f (x) on each sub-stencil of r − 1 points. There exist several expressions for I r −1 i,k . For example, the indicators designed in Jiang and Shu (1996) and Levy et al. (1999) arẽ These indicators are suitable for the approximation of the conservation laws (1) with discontinuities in their solutions. The results obtained using them for approximating the derivative values are not satisfactory because they do not correctly detect kink discontinuities. For approximating the derivative values, suitable for the approximation of the Hamilton-Jacobi equations (1), the measurement of the smoothness indicators should start from the second derivative (Jiang and Peng 2000;Amat and Ruiz 2017). In this work, we use the formula given in Amat and Ruiz (2017) adapted to non-uniform grids: Recently, a new WENO-2r algorithm has been introduced in Amat et al. (2020) and Amat et al. (2021). It consists of exploiting a recursive process to calculate the nonlinear weights with the aim of obtaining progressive order of accuracy of the approximation close to discontinuities. In this paper, we adapt these ideas to obtain a new progressive WENO interpolator to approximate the derivative values. The paper is divided into the following sections: in Sect. 2, we show the algorithm and construct the new method for r = 3 and for r = 4. In Sect. 3, we generalize the results for any r and give a general expression for nonlinear optimal weights. In Sect. 4, we present a strategy to compute efficiently the smoothness indicators and study the order of accuracy. Finally, some numerical experiments and some conclusions are presented.

New WENO with progressive adaptation to discontinuities: cases r = 3 and r = 4
The new algorithm designed by Amat et al. in Amat et al. (2020) consists of using the Aitken's interpolation process Abramowitz and Stegun (2010) to calculate progressive linear weights. For simplicity, we present two examples and in the next section we show the generalization for any r .

New WENO for r = 3 with progressive adaptation to discontinuities in non-uniform grids
Therefore, by Aitken's process, we have that Analogously, dp 3 Then, it is clear that: then, we define: where the smoothness indicators I 3 i,0,k 1 , k 1 = 0, 1 in (14) will be defined in Sect. 4 based on those introduced in Levy et al. (1999). Finally, we apply WENO with the new nonlinear weights, i.e. we calculate: with I 2 i,k , k = 0, 1, 2, being the smoothness indicators introduced in Eq. (11). Using the same reasoning, as a corollary, we obtain the formulas for new nonlinear weights in the uniform grid case. Thus, we have that In this case, we will write C k 1 k 2 ,k 3 (x i ) = C k 1 k 2 ,k 3 for any k 1 , k 2 , k 3 . Thus, we get: We can see that the linear optimal weights, in this case, are similar to the weights shown in Levy et al. (1999).
Finally, we repeat the same steps to obtain nonlinear weights.

New WENO with progressive adaptation to discontinuities for r = 4 in non-uniform grids
We start with the polynomial } and apply the same process, i.e. we express the derivative value of this polynomial at x i as a combination of the derivative values of the polynomials p 5 i,0 and respectively; then, again, using Aitken's algorithm, we obtain: Analogously, we represent p 5 i,0 as a combination of p 4 i,0 and p 4 i,1 as and repeating the same steps we get: We repeat this procedure again to obtain: then, the optimal progressive weights are: Thus, we can definẽ and l = 4, 5, k = 0, 1, k 1 = k, k + 1; and with this vector, we apply the classical WENO algorithm.
As a corollary, if the grid is uniform we get: Thus, it is easy to check that Finally, we perform the same algorithm to obtain the nonlinear weights.

General new WENO algorithm for derivative values and general optimal weights
In this section, we will generalize the method for any r . To compute the linear weights for each "level" we can prove the following lemma following the same ideas as in Amat et al. (2020), i.e. using Aitken's process as in Sects. 2.1 and 2.2. then To obtain the interpolators p l+1 i,k , p l i,k and p l i,k+1 , the stencils used are {x i−(r −1)+k , . . . , x i−(r −1)+k+l+1 }, {x i−(r −1)+k , . . . , x i−(r −1)+k+l } and {x i−(r −1)+k+1 , . . . , x i−(r −1)+k+l+1 }, respectively. Then, using Aitken's interpolation process (Abramowitz and Stegun 2010), we obtain: If we differentiate the previous expression: It is trivial to check that for all i, we have: As a corollary, we can calculate the optimal weights if the grid is uniform.
Corollary 3.2 Let 0 < r − 1 ≤ l ≤ 2r − 3 and 0 ≤ k ≤ (2r − 3) − l, if the grid is uniform, i.e. there exists h = 1/J such that x j = a + j · h, 0 ≤ j ≤ J and we denote as C l k,k and C l k,k+1 the values which satisfy: dp l+1 Proof It is direct by Lemma 3.1 and We apply Eq. (15) for each level and get: Thus, if we define the weights and the vector C r−1 i,k with 0 ≤ k ≤ r − 2 as where C r −1 i,k,k , C r −1 i,k,k+1 are defined in Eq. (16), then we get that: with C r −1 i,k , 0 ≤ k ≤ r − 1 the optimal weights obtained in Prop. 1.1, Tables 1 and 2 for r = 2, 3 and 4. It is analogous in the uniform grid case.
We have a tree scheme, where each branch produces a polynomial of a determined degree. Now, the idea is to use all the points which are not contaminated by a discontinuity. To follow this idea, we reduce this branch to O(h 2 ) using nonlinear weights as follows: We substitute in Eq. (21) the linear weights by nonlinear weights where θ, are the parameters above-mentioned and I l i,k,k 1 are smoothness indicators defined at level l = r , . . . , 2r − 3. Therefore, the last ingredient of this scheme is to define the indicators to "remove" (to obtain O(h 2 )) the non-suitable branch. We will use the strategy used in Amat et al. (2020) explained in detail in Sect. 4.
Finally, we define the new weights as: UsingC r −1 i , we apply classical WENO and calculateP i (

Smoothness indicators and analysis of the accuracy
Let us start with the analysis of the smoothness indicators presented by Amat and Ruiz in Amat and Ruiz (2017), Eq. (12): . . , r − 1, in non-uniform grids. We use the same ideas introduced in Aràndiga et al. (2010). First of all, at the smooth zones we obtain I r −1 i,k = O(h 4 ) using the following adapted result proved in Amat and Ruiz (2017).
Lemma 4.1 Let 0 ≤ k ≤ r − 1 and p r −1 i,k the interpolating polynomial of degree r − 1 ≥ 2 of f that uses the nodes of the stencil S r i,k . Then, the smoothness indicator obtained through (11) satisfy We will prove the following auxiliary lemmas using the ideas presented in Aràndiga et al. (2010).

Lemma 4.2 Let
Proof Let p r −1 i,k , p r −1 i,k 1 be the two interpolating polynomials of f of degree r − 1 ≥ 2 at nodes in the stencil S r i,n n = k, k 1 . As k = k 1 and 0 ≤ k, (25) Therefore, we get Thus, we obtain: Lemma 4.3 Let 0 ≤ k, k 1 ≤ r − 1, 1 ≤ θ and let I r −1 n , n = k, k 1 be smoothness indicators of f on the stencil S r i,n = {x i−(r −1)+n , · · · , x i+n }, and p r −1 i,k 1 be the interpolating polynomial of f at nodes in the stencil S r and Proof The proof of Eq. (28) is direct by Lemmas 4.1 and 4.2. To prove Eq. (29), we use the following algebraic manipulation (see Aràndiga et al. (2010)): Notice that if I r −1 i,k 1 = O(h m ) with m > 4 and 0 < < h m is a fixed value, then Lemma 4.3 is not satisfied and the optimal order is not obtained. This is explained in detail in Aràndiga et al. (2010). In all our experiments, as we have mentioned in Section 1, we have set = 10 −16 . This way, the order of accuracy of the smoothness indicators will be O(h 4 ) at the smooth parts of the data. In the rest of the paper, we always consider these conditions. Proposition 4.4 Let 0 ≤ k ≤ r − 1, 1 ≤ θ and ω r −1 i,k the nonlinear weights defined in Eq. (10), then with m = 2 if the discontinuity is in the function and m = 1 if the discontinuity is in the first derivative.
Proof Let 0 ≤ k, n < r − 1. If f is smooth in S r n , then using (10) and Lemma 4.1: with m = 2 being the discontinuity in the function and m = 1 the discontinuity in the first derivative. Finally, if f is smooth in {x i−(r −1) , . . . , x i+(r −1) }, we fix a value 0 ≤ k 1 ≤ r − 1, and using (29) in Lemma 4.3, we get for all 0 ≤ k ≤ r − 1: Therefore, To construct the smoothness indicators for each "level", we consider where the discontinuities are placed. To clarify the explanation, we start with an example, for r = 4. In Sect. 2.2, we have seen that we should define I l i,k,k 1 with l = 4, 5, k = 0, 1, k 1 = k, k + 1, then we have the following scheme: We suppose that the discontinuities are far enough from each other, i.e. there only exists one discontinuity at an interval. Thus, we have the following cases: • There does not exist any discontinuity: we use all the weights.  Fig. 3, we write in red the branch that is not used because a discontinuity is contained in the stencil; in green, the stencil that is used although it contains a discontinuity, and in blue the stencils which do not contain any discontinuity. Therefore, we define I 6 i,0,0 = I 3 i,0,0 and I 5 i,1,1 = I 3 i,1,1 .
Diagram showing the structure of the optimal weights needed to obtain optimal order of accuracy for r = 4 Using this example, we introduce the following definition.
Second case: discontinuity at (x i+1 , x i+2 ). Red: polynomials that are not used. Blue: polynomials used. Green: polynomial affected by a discontinuity, but not eliminated branch Third case: discontinuity at (x i , x i+1 ). Red: polynomials that are not used. Blue: polynomials used. Green: polynomial affected by a discontinuity but not eliminated branch Definition 1 General smoothness indicators Let I r −1 i,k , with k = 0, . . . , r − 1 be the smoothness indicators shown in Eq. (11), then we define the smoothness indicators for any r as: with r ≤ l ≤ 2r − 3.
The following lemma is similar to Prop. 4.4 for nonlinear weights at each level.

If I l i,k,k is affected by a singularity, then ω
Here, m = 2 if the discontinuity is in the function, and m = 1 if the discontinuity is in the first derivative.
Proof It is similar to Prop. 4.4: from C l i,k,k + C l i,k,k+1 = 1, using the definition given in (22), we have that: , by Eq. (29) in Lemma 4.3, we get: Analogously for ω l i,k,k+1 . 2. If I l i,k,k = O(h 4 ) and I l i,k,k+1 = O(1), if the discontinuity is in the function or I l i,k,k+1 = O(h 2 ) if the discontinuity is in the first derivative: If we analyse the examples presented, we can determine a rule for the weights and try to prove it. In the first case, Fig. 2, the discontinuity is in the interval (x i+(l 0 −1) , x i+(l 0 −1) ) = (x i+2 , x i+3 ), (i.e. l 0 = 3), and we can see that the branch marked in red is the one corresponding to p 5 i,1 ; thus, ω l i,0,1 = ω (r −2)+l 0 i,0,1 = ω 5 i,0,1 should be O(h 2 ). In the second case, Fig.  3, the discontinuity is in the interval (x i+1 , x i+2 ), (l 0 = 2), and the red branches in those corresponding to p 4 i,1 , then ω (r −2)+l 0 i,0,1 = ω 4 i,0,1 and ω 5 i,0,1 should be O(h 2 ). We prove that the weights of the different branches which contain a discontinuity go to 0 as h → 0, in the following lemma.
, then for all (r − 2) + l 0 ≤ l ≤ 2r − 3 the nonlinear weights defined in Eq. (22) satisfy: , with m = 2 if the discontinuity is in the function and m = 1 if the discontinuity is in the first derivative. −1) , . . . , x i } and l 0 ≥ 1. Subsequently, the stencil used to calculate if the function contains a kink discontinuity. By Def. 1, we have that I l i,0,1 = I r −1 i,l+2−r = O(hm),m = 0, 2. Using Lemma 4.5, we finish the proof.
, then for all r ≤ l ≤ l 0 + (r − 2) − 1 the nonlinear weights defined in Eq. (22) satisfy: Proof Let 2 ≤ l 0 ≤ r − 1 and r ≤ l ≤ l 0 + (r − 2) − 1. We analyse the stencils used to calculate I l i,k,k and I l i,k,k+1 with 0 ≤ k ≤ l 0 + (r − 2) − 1 − l. First of all, we take into account two previous considerations: From r ≤ l, we get: and from k ≤ l 0 + (r − 2) − 1 − l: By Def. 1, we have that: Using (34), the stencil does not cross the discontinuity and . By Lemma 4.5, we get the result.
With the ingredients presented in the previous sections, we can prove the following lemma. We suppose that the isolated discontinuity is to the right of the point x i . By symmetry, the analysis for the left side would be similar.
with m = 2 if the discontinuity is in the function and m = 1 if the discontinuity is in the first derivative, and dp r +l 0 −2 Proof Applying Eq. (19) to the interpolatory polynomial p r +l 0 −2 i,0 and Lemmas 4.6 and 4.7, we obtain that: with m = 2 if the discontinuity is in the function and m = 1 if the discontinuity is in the first derivative, with Therefore, we have proved that, when the isolated discontinuities affect the stencil, we get maximal order in a neighbourhood of the interval where it is located. In the next section, we show some experiments that confirm this theoretical result.

Numerical experiments and conclusions
In this section, we present some numerical experiments to verify the theoretical results shown in previous sections. We will use the following function: discretized in the interval [− π 6 , 1 − π 6 ] using J q = 2 q + 1 uniform spaced points where η = 0, 10. In the first case, the function f 0 is continuous, but it presents a discontinuity in the first derivative, and in the second case, f 10 has a jump discontinuity. We denote the interval where the discontinuity is contained as (x i−1 , x i ) and compute the error as with ap being the approximation of the derivative values using the following methods: • [lin-(2r − 2)] Linear Lagrange method using a centred stencil of 2r − 1 points.
Finally, we compute the numerical order of accuracy using the formula: We start using a large stencil of five points, i.e. r = 3, and, using the linear method, we expect to obtain four consecutive points where the order is lost. Using classical WENO method, we will typically obtain order 3 at the points where the order is lost using the linear method. Using the new algorithm, we expect to obtain progressive order of accuracy (r − 1, r , 2r − 2) = (2, 3, 4). We show the results for r = 3 in Tables 4, 5 and 6 for f 0 and  Tables 10, 11 and 12 for f 10 . For r = 4, we can see the results in Tables 7, 8 and 9 for the first  experiment and Tables 13, 14 and 15 for the second one. It is clear that the new algorithm produces better approximation close to the discontinuity. We obtain the same results to the left of the discontinuity. Finally, it is important to remark that if the function is continuous, but presents a discontinuity in the first derivative, i.e. a kink, we cannot use the smoothness indicators introduced in Eq. (11).  Table 7 Grid refinement analysis for the linear Lagrange method for r = 4 algorithm for the function in (40) Table 8 Grid refinement analysis for the classical WENO method for r = 4 algorithm for the function in (40) Table 9 Grid refinement analysis for the new WENO method for r = 4 algorithm for the function in (40)

Table 10
Grid refinement analysis for the linear Lagrange method for r = 3 algorithm for the function in (40)

Table 11
Grid refinement analysis for the classical WENO method for r = 3 algorithm for the function in (40)

Table 12
Grid refinement analysis for the new WENO method for r = 3 algorithm for the function in (40)

Table 13
Grid refinement analysis for the linear Lagrange method for r = 4 algorithm for the function in (40)

Table 14
Grid refinement analysis for the classical WENO method for r = 4 algorithm for the function in (40)

Table 15
Grid refinement analysis for the new WENO method for r = 4 algorithm for the function in (40)

Conclusions
In this work, a new WENO algorithm with progressive order of accuracy close to discontinuities has been introduced. It allows to calculate approximations of derivative values of a function using regular or non-regular grids. It is based on the same ideas used in Amat et al. (2020) which consist in using Aitken process to calculate the nonlinear weights. The explicit formulas for the optimal weights have been showed and the order of accuracy in each interval has been proved. Finally, some experiments have been presented that confirm the theoretical results obtained.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.