A non-separable progressive multivariate WENO-$2r$ point value

The weighted essentially non-oscillatory {technique} using a stencil of $2r$ points (WENO-$2r$) is an interpolatory method that consists in obtaining a higher approximation order from the non-linear combination of interpolants of $r+1$ nodes. The result is an interpolant of order $2r$ at the smooth parts and order $r+1$ when an isolated discontinuity falls at any grid interval of the large stencil except at the central one. Recently, a new WENO method based on Aitken-Neville's algorithm has been designed for interpolation of equally spaced data at the mid-points and presents progressive order of accuracy close to discontinuities. This paper is devoted to constructing a general progressive WENO method for non-necessarily uniformly spaced data and several variables interpolating in any point of the central interval. Also, we provide explicit formulas for linear and non-linear weights and prove the order obtained. Finally, some numerical experiments are presented to check the theoretical results.


Introduction and review
In the last years, WENO methods have been developed and used in several applications, mainly to obtain numerical solutions of partial differential equations (PDEs) but also in other fields, such as image processing or computer-aided design (see e.g.[6,14,15]).The idea is to compute a non-linear combination of interpolations through polynomials of degree r, aiming to obtain the maximum possible order 2r when the data is free of discontinuities, and order r + 1 at the non-smooth parts.In what follows, we briefly review the method.Let us denote as X a uniform partition of the interval [a, b] in J subintervals: and consider the point-value discretization of a piecewise smooth function f at the nodes x i , . In this setting, the WENO method with 2r nodes interpolates at the mid-point of the interval (x i−1 , x i ), denoted by x i− 1  2 , using the stencil S 2r 0 = {x i−r , • • • , x i+r−1 }.We construct this interpolation by means of the following convex combination: where ω r k ≥ 0, k = 0, • • • , r − 1 are non-linear (data-dependent) positive weights such that The values of the weights ω r k are designed to obtain an order of accuracy 2r at x i− 1 2 when the function is smooth in the large stencil, as follows: There are optimal weights C r k ≥ 0, with k = 0, . . ., r − 1 that satisfy the following equality: A formula for these values is obtained in [6]: As stated in [15], the weights ω r k satisfy at smooth zones, with κ ≥ r − 1, thus assuring that the interpolation in (1) attains order of accuracy 2r In [12,15], (3) is achieved via the following expressions: In the previous expression, the parameter t is an integer that assures maximum order of accuracy close to the discontinuities.The parameter ǫ > 0 is introduced to avoid divisions by zero, and is usually forced to take the size of the smoothness indicators at smooth zones.In our numerical tests, we will set it to ǫ = h 2 .The values I r k are called smoothness indicators for f (x) on each sub-stencil of r points.
Since it was first introduced in [15], several successive improvements have been proposed for WENO algorithms.They have been focused on enhancing their accuracy, efficiency, and robustness in approximating solutions to hyperbolic conservation laws, but also on extending their applicability to other contexts such as approximation and interpolation of data.Initially introduced by Liu, Osher, and Chan in 1994 [15], WENO schemes have gone through several advancements.The WENO schemes proposed by Jiang and Shu in 1996 [12], improved the original idea in [15], by proposing new smoothness indicators inspired by the measure of the total variation.These smoothness indicators were more capable of detecting discontinuities and allowed to extended the idea of WENO to higher orders of accuracy.However, the classical WENO scheme might experience a loss of accuracy when encountering critical points in the solution: for instance, a fifth-order WENO scheme may only attain third-order accuracy near smooth extrema [16].In fact, in [16], the authors proposed the WENO-M scheme, that not only addressed the accuracy issue but also marked the first significant improvement in the solution quality near shocks and high gradients.However, the introduced mapping method proved to be computationally expensive.This study led to the publication of the article [9] where the authors proposed the WENO-Z scheme, that uses a new set of weights, derived from previously unused information within the classical WENO scheme: a higher-order global smoothness indicator constructed through a linear combination of the original smoothness indicators.This scheme achieved superior results with almost the same computational effort as the classical WENO method.After that, many variants of WENO schemes were proposed, starting from these two methodologies (see, e.g.[10,13,17]).In the context of data approximation, we proposed a first improvement of WENO algorithm in 2018 [2], where the algorithm attained a progressive order of accuracy close to the discontinuities using a recursive formulation of the WENO weights.Using WENO interpolation, maximum order of accuracy is obtained in smooth parts of the data but the accuracy is reduced to order r + 1 when, at least, a discontinuity crosses a stencil.In [3], Amat, Ruiz, Shu and Yáñez present a new WENO method using the Aitken-Neville algorithm to obtain progressive order of approximation.This method is introduced for the point-value discretization in a uniform grid, and to approximate at the mid-points of the intervals.In [4], the algorithm is extended to calculate the derivative value of a function knowing its evaluation in a non-uniform grid and for any point of the considered interval.
In this paper, we extend this method to n dimensions and prove its theoretical properties.We start with dimension 1 in Section 2, following the ideas presented in [3,4].We provide a new recursive algorithm to obtain a non-linear scalar approximation for any point in the central interval in Section 2.Then, in Section 3, we review the multidimensional Lagrange interpolation using tensor product, and the WENO version designed by Aràndiga et al. in [7].Subsequently, we introduce, in Section 4, our general method for any dimension n starting with n = 2, and give an explicit expression for the optimal weights.Afterwards, we generalize the method to any number of dimension n, and prove that the approximation attains the maximum order of accuracy when the nodes are in a region where f is smooth, and has an increasing order of accuracy depending on where the isolated discontinuity is.In Section 6, we present the smoothness indicators.The theoretical results are confirmed by the numerical examples, that are presented in Section 7. Finally, some conclusions are drawn in the last section.

A new univariate WENO-2r algorithm with progressive order of accuracy close to discontinuities
Recently, Amat, Ruiz, Shu and Yáñez in [3] introduced a new centered WENO-2r method, which consists in using all the points free of discontinuities to interpolate a value at the mid-point of an interval.In this section, we generalize this interpolation to non-uniform grids.Let [a, b] ⊂ R be an interval, we consider a non-uniform grid , and construct the polynomial which interpolates at these nodes, that we denote by p 2r−1 0 .Using the Aitken-Neville formula, we express this polynomial using the following expression which involves the two interpolatory polynomials, p 2r−2 j0 , j 0 = 0, 1 with respective stencils: where C 2r−2 0,0 (x * ) and C 2r−2 0,1 (x * ) are the optimal weights.We repeat this process with each polynomial up to degree r.A representation of this process can be seen in Figure 1 (where we have removed the dependence on the value x * ), [3].Thus, we obtain the general explicit expression: We calculate these values through the following lemma (the version for uniform grids, and when the interpolation is at mid-points is proved in [3]). where Proof.The proof is direct by taking into account that the stencil for constructing p l+1 j is {x i−r+j , . . ., x i−r+j+l+1 } and for p l j and p l j+1 are {x i−r+j , . . ., x i−r+j+l } and {x i−r+j+1 , . . ., x i−r+j+l+1 } respectively.As Corollary, we obtain the values showed in [3] for mid-point interpolation. then From Eq. ( 5), we can deduce the following recursive method: with the non-linear weights defined as: where C r j,j1 (x * ), j 1 = j, j + 1 are determined in Eq. ( 7).The values I l j,j1 , explained in subsection 2.1, are smoothness indicators.
With this formulation, the approximation is defined by being p2r−1 0 (x * ) the result of the recursive process, Eq. (10).

On smoothness indicators and the accuracy of the new univariate progressive WENO method
The choice of the smoothness indicators is crucial.To design the new non-linear algorithm, we only calculate the smoothness indicators at level l = r, and we use them in all the levels, as we can see in this subsection.In [6] it is proved that if the smoothness indicators satisfy the following conditions: P1 The order of a smoothness indicator that is free of discontinuities is h 2 , i.e.
P2 The distance between two smoothness indicators free of discontinuities is h r+1 , i.e. let k, k ′ ∈ {0, 1, . . ., r − 1} be such that there does not exist any discontinuity in S r k and S r k ′ , then P3 When a discontinuity crosses the stencil S r k then Then, the optimal weights satisfy Eq. (3) for determined parameters t and ǫ (Proposition 2 in [6]).We exploit this result to construct the smoothness indicators in each step of the Aitken-Neville's algorithm.We will use the following definition (Definition 4.1 in [3]).
Definition 1.Let l = r, . . ., 2r − 2, and I r k , with k = 0, . . ., r − 1, be the smoothness indicators with properties P1,P2 and P3.Then, we define the smoothness indicators at level l as: Therefore, in each step, we will remove the part which is "contaminated" by a discontinuity.In particular, we can prove the following proposition adapting the proof of Proposition 4.4 in [3] to any point x * ∈ (x i−1 , x i ).
and let I l k,k , and I l k,k+1 be smoothness indicators defined in Definition 1.The following weights: with, and with C l k,k (x * ) + C l k,k+1 (x * ) = 1, will fall under one of the following cases: 1.If neither I l k,k nor I l k,k+1 are affected by a discontinuity, then ωl We have the following main result: , and Ĩ2r−1 (x * ; f ) the result of the recursive process, Eq. (10 By symmetry, we only analyze when there exists an isolated discontinuity at an interval [x i−1+l0 , x i+l0 ], l 0 = 1, . . ., r−1 (analogously, we obtain the equivalent symmetric results for [x i−r+l0 , x i−r+l0+1 ], l 0 = 0, . . ., r − 2).
In the next sections, we generalize this method to multi-dimensions.We also introduce some possible smoothness indicators which satisfy P1,P2, and P3 in Section 6.

Comparison of multivariate linear Lagrange interpolation with the non-linear WENO method in Cartesian grids
In this section, we briefly review the Lagrange interpolation problem in several variables when the data are located in Cartesian grids (non-necessarily equally-spaced) and we construct a multivariate WENO method following the ideas presented in [7] for two variables.We present the necessary ingredients to extend the progressive WENO-2r algorithm introduced in Section 2 to several variables.

Multivariate linear interpolation
Let us start supposing a j , b j ∈ R with a j < b j , j = 1, . . ., n, and an hypercube denoted as: We call the points of a grid for each interval as: and define We suppose an unknown function f : n j=1 [a j , b j ] → R, and consider our data as the evaluation of this function in the points of the Cartesian grid, i.e.
Let the polynomials of n variables be denoted by: when τ 1 = . . .= τ n = τ we call it Π τ n := Π τ1,...,τn n .We consider (i 1 , . . ., i n ) ∈ N n to be the reference index where the approximation is centered at, and r ∈ N such that 0 ≤ i j − r ≤ i j + r − 1 ≤ J j , j = 1, . . ., n, and the centered stencil: then, the problem consists in calculating a polynomial of degree 2r − 1 such that: To do so, we will use the Lagrange base of polynomials: then, the unique polynomial is It can be checked (see [7]) that if j ], and f ∈ C 2nr (R) then the error satisfies: and if m = (m 1 , . . ., m n ) ∈ N n with m j ≤ 2r − 1, j = 1, . . ., n, we get where

Multivariate WENO interpolation in Cartesian grids
Using the same notation as in the previous section, the goal is to construct a non-linear interpolant in multidimensions in the same way as we reviewed in 1d, i.e., an interpolant with maximum order 2r in the smooth parts and with order r + 1 when there exists a discontinuity which crosses, at least, one small stencil.
In this case, we introduce the method supposing that we want to obtain an approximation of the function f evaluated at any point of the hypercube The construction is similar to the 1d case: Firstly, we make the linear combination of interpolants of lower degree, r, if k = (k 1 , . . ., k n ) then where C r k (x * ) are the optimal weights, and it is not difficult to prove that they have the explicit form: being C r kj (x * j ), k j = 0, . . ., r − 1, j = 1, . . ., n the optimal weights in 1d.In the case that we want to approximate at the mid-point of the hypercube, i.e.
, then they are defined in Eq. ( 2).The polynomials p r k interpolate at the nodes } × {x Then, we replace the optimal weights by non-linear ones using the following formula being I r k the smoothness indicators with the same requirements pointed out in Section 2, i.e.: P1 The order of the smoothness indicator free of discontinuities is h 2 , i.e.
P2 The distance between two smoothness indicators free of discontinuities is h r+1 , i.e. let k = (k 1 , . . ., k n ), and P3 When a discontinuity crosses the stencil S r k then: An example of smoothness indicators will be introduced in Section 6.The parameters ǫ and t are chosen to obtain maximum order.In our case, following [6] and [7] we will take ǫ = h 2 and t = 1 2 (r + 1).Therefore, we state the next result, which is similar to Proposition 2 in [6] and Theorem 1 in [7].
at smooth regions, O(h r+1 ), if, at least, one stencil lies in a smooth region.

A new progressive bivariate, n = 2, WENO method
The idea of this new method is to reach the maximum possible order of accuracy when one discontinuity crosses the largest stencil.For this purpose, we use the Aitken-Neville-based algorithm developed in multi-dimensions, [11].We start with a polynomial of degree 2r − 1 and decompose it in 2 2 polynomials of degree 2r − 2 obtaining as weights polynomials of degree 1 that will be replaced by non-linear weights, depending on the location of the discontinuity.Afterwards, we continue decomposing the 2 2 polynomials of degree 2r − 2 in polynomials of degree 2r − 3, and so on.In each step, the non-linear weights determine the stencils free of discontinuity, which will be used to approximate the value.The procedure is similar to the method expounded in Section 2. For ease of reading, we start with n = 2 in this section, and then we extend our results to any dimension n in Section 5.

A first example: A progressive bivariate WENO method with r = 2
Let us start by designing a bivariate WENO method with r = 2.We suppose a non necessarily uniform grid in 2 ), l j = 0, . . ., J j , with j = 1, 2. Let (i 1 , i 2 ) ∈ N 2 such that 0 ≤ i j − 2 and i j + 1 ≤ J j , with j = 1, 2. We determine the largest stencil and we want to interpolate at any point 2 ].We compute the polynomial such that p 3 (0,0) (x 2 2 ) ∈ S 4 (0,0) .Again, we can express it as the sum of 2 2 polynomials that interpolate in the stencils: } × {x (1,1) using, as we mentioned before, the Aitken-Neville-type formula given in [11]: with Finally, we define the non-linear weights as: , where α2 (0,0),k defined as: , where I 2 k , k ∈ {0, 1} 2 are smoothness indicators which satisfy the properties P1, P2 and P3.The parameters t and ǫ are defined in Section 3.2 as h 2 and r + 1 = 3 respectively.Thus, we define the new interpolant as: If we apply it at mid-points, the interpolant presented in [7] and the method showed in Section 3.2 for n = 2 are similar if the same smoothness indicators are used.Therefore, when r = 2 there are not differences between the progressive WENO and the classical one.We show the construction of the new method for r = 3.

The case r = 3
In this case, we use all the points of the following largest stencil: } × {x }.

The general case
Now, we construct the new bivariate WENO-2r algorithm using a larger stencil of (2r) 2 points S 2r (0,0) = {x The goal is to obtain maximum order of approximation if there exists an isolated discontinuity.We start representing the polynomial of degree 2r − 1 as follows: , and repeat the process up to degree r: being Ω l = j l + {0, 1} 2 with 1 ≤ l ≤ r − 3. Consequently, we determine the values C l j,j1 with r ≤ l ≤ 2r − 2, and j 1 ∈ j + {0, 1} 2 , proving the following lemma.

A new progressive multivariate WENO method
In this section we generalize the method designed in Section 2 for 1d, and in Section 4.3 for 2d.We follow the same steps: First of all, we consider the data as in Section 3.1, i.e., if n j=1 [a j , b j ] is an hypercube and are the points of a non-regular grid, we suppose our data as the evaluation of a unknown function at these points Let (i 1 , . . ., i n ) ∈ N n and r ∈ N be such that 0 j ], and the centered stencil: with 0 = (0, . . ., 0).We compute the interpolatory polynomial of degree 2r − 1 with nodes S 2r 0 , p 2r−1 0 (x * ), and we express it as the combination of 2 n polynomials of degree 2r − 2. Thus, we get Now, we repeat this process up to polynomials of degree r (for simplicity remove the dependence of x * ), we denote as Ω l = j l + {0, 1} n : with being C l ji,ji , and C l ji,ji+1 , i = 1, . . ., n the weights defined in Eq. ( 9).We establish the recursive process as Eq. ( 10): being for l = r, . . ., 2r − 2, and k ∈ {0, . . ., 2r − 2 − l}, where I l k,k1 are the smoothness indicators determined by the following formula: being I r k , with 0 ≤ k j ≤ r − 1, j = 1, . . ., n, smoothness indicators satisfying the properties P1, P2, and P3.Therefore, the approximation will be: Finally, we calculate the order of accuracy stating the next theorem.

Smoothness indicators
In this section, we present some smoothness indicators which satisfy the above mentioned properties.We generalize the smoothness indicators introduced by Aràndiga et al. in [7] which are an adaptation of the ones presented in [5].The idea is to design some functionals that fulfill P1, P2, and P3. Given Hence, we define where Γ = j ], J = {0, 1, . . ., r} n \ {0} and smoothness indicators to get an expression which is scale independent.
We suppose that p is a polynomial interpolant at a stencil of (r + 1) n points, then where E is defined in Eq. ( 18), and ||l|| ∞ = max j=1,...,n |l j |.Then, if k, k ′ ∈ {0, . . ., r − 1} n , we get: Finally, if a discontinuity crosses the stencil S r k then: since some of the quadratic terms of I r k will not converge to 0.
Remark 6.1.An evaluation of these smoothness indicators for 2d in gridded data is shown in [7].Also, some adaptation to smoothness indicators with better capabilities, and computationally more efficient than those introduced in [1], or in [8] can be performed, but the process requires an study on the order of accuracy, which exceeds the scope of this paper.

Numerical experiments
In this section, we check our theoretical results through some numerical examples.We divide it into four parts, starting with 1d experiments for uniform and non-uniform grids and performing some tests for the multivariable case.In both cases, the order of accuracy is analyzed measuring the error in a set of points with a determined grid and refining it.For this purpose we define the error in a finite set Υ ⊂ R n as: Secondly, in all cases, we study the behavior of the resulting interpolator in the zones close to the discontinuities and conclude that by employing our new non-linear algorithm, if the discontinuity is isolated, the Gibbs phenomenon is avoided.

Examples in 1D for uniform grids
We perform the first experiment using the same function studied in [3]: discretizing it with To analyze the order close to the discontinuity we locate the point 0 in each level ℓ, compute the errors in the adjacent intervals, and estimate the numerical order of accuracy.Therefore, let j ℓ 0 ∈ N such that 0 ∈ [x j ℓ 0 −1 , x j ℓ 0 ], then we interpolate the value of the function f at 10000 points equally spaced in the interval [x j ℓ 0 −5 , x j ℓ 0 +4 ], denoted by J ℓ , determine the error using Eq. ( 42), E ℓ s = E(J ℓ ∩ [x j ℓ 0 +s−1 , x j ℓ 0 +s ]) and approximate the numerical order in each interval s as: This allows us to analyze the progressively increasing order in the neighboring intervals to the one containing the discontinuity.For r = 3, we can see in Table 1 that the order of accuracy increases as we proved in Theorem 5.1.Also, for r = 4, Table 3, we observe the same behavior of the numerical order.However, when we apply the classical WENO algorithm, in both cases r = 3, 4, Tables 2 and 4, the numerical order is reduced to r + 1, in all the intervals where one of the small stencils is contaminated by the singularity.Note that when we use our algorithm, we interpolate in any point of the interval, but when using the classical WENO method, we can only interpolate at the mid-point.

Examples in 1D for non-uniform grids
In this subsection, we perform two experiments: one with the function defined in (43) and, also, with the function studied in [6]: We construct a non-uniform grid {x i } 32 i=0 with x i ∼ U [0, 1], being U [0, 1] a uniform distribution in [0, 1] and interpolate the values of the function in a uniform grid of 10000 points.The result is showed in Figure 5, we can see that the two interpolants avoid Gibbs-phenomenon.We compare our method with classical WENO adapted to non-uniform grids, we change the optimal weights in (2) using the formulas presented in Lemma 2.1.In Tables 6 and 8 we can observe similar results to those obtained for uniform grids: the order of accuracy in the adjacent cells to the isolated discontinuity, o l 2 and o l −2 , is smaller than the one obtained using the new WENO algorithm.These numerical results are consistent with the theoretical ones.To analyze the numerical order we take the non-uniform mesh, {x 5 i } 32 i=0 , we divide it in each level using the following formula: x ℓ+1 2j+1 = 1 2 (x ℓ j + x ℓ j+1 ), j = 0, . . ., 2 ℓ − 1, for ℓ ≥ 5 and compute the numerical order following Eq.(44).Again, we see, in Tables 5 and 7, that the behaviour of the numerical order is the expected one, and proved in Theorem 5.1.

Examples in 2d for uniform grids
We start with a smooth function to check numerically that the order is 2r.Thus, we consider the function     and interpolate it using a Cartesian grid in the square [−1, 1] 2 : at the points applying the new WENO algorithm for r = 3.Note that the set of points B ℓ is conformed by a convex combination of points of the grid.We have chosen this collection of points, but we could have selected any other.The result is shown in Fig. 6, and we observe that it is similar to the original function f 3 .In order to compute the numerical order, we use the same strategy presented in Section 7.1, we calculate E ℓ = E(B ℓ ) and approximate the order as in Eq. ( 44).
Regarding the table presented in Fig. 6 it is clear that the order of accuracy is the expected one.
1.6882e-08 5 1.3037e-08 5.94 6.4910e-11 8.02 6 2.0584e-10 5.98 2.5102e-13 8.01 7 3.2245e-12 6.00 1.4433e-15 7.44 In order to analyze the order close to the discontinuities, we perform an experiment with the function: We take the mid-point ) = (0, 0) of the cartesian grid X ℓ ; the set (see Figure 7.(a) red big points), and calculate an approximation to the function in the points (see Figure 7.(a) gray, green, yellow and blue small points).We chose this set for the convenience of analyzing the order, but any other random sample of points could have been selected.A region between four data points is defined as: , and the errors and the numerical orders in each region as: Now, we compare with a modification of the classical WENO method adapted to work using tensor products.Note that the WENO-2d method defined in [7] is constructed for uniform-grids.In this case, we reformulate the classical WENO to approximate at any point in the considered interval and for non-uniform-grids.Theoretically, we have proved that the order of accuracy obtained at the green points in Figure 7.(a) is equal to 4, as there is one unique square stencil that is not contaminated by the discontinuity.For the yellow ones, 5 is the theoretical order of accuracy, and 6 for the rest of the points.Using the classical WENO algorithm, we only distinguish two zones, order 4 (yellow and green zones) and order 6 for the rest of the points.We can check this fact in Tables 9 and 10, where we determine the order, o 7 (s1,s2) in each region.In this case, in the smooth part, we obtain numerical order 7.The interpolation is shown in Figure 7.(c).We can see again that the resulting intepolator avoids Gibbs phenomenon (Figure 7  Then we calculate the errors and numerical orders in Ω ℓ = X ℓ ∩ [−5.5h ℓ , 5.5h ℓ ] 2 , as in previous subsections, with E ℓ = E(Ω ℓ ).We can observe again in Table 11, that the results obtained in the numerical experiments satisfy the theoretical ones.To finish these numerical tests, we perform an interpolation of the function  being X7 a non-uniform grid constructed employing the formula described in Eq. (50) (see Figure 8(a)).We show the result in Figure 8(c) and (d) confirming that Gibbs phenomenon does not appear.

Conclusions and future work
In this work, we have presented a new WENO method with a progressive order of accuracy close to the discontinuities.This new algorithm is simpler computationally than the one presented in [3] and more general.In particular, our main contributions can be summarised in the following: we design a new recursive method to compute the interpolation; we construct non-linear explicit formulas to calculate the weights to interpolate any point in the central interval; we propose an algorithm for non-uniform grids, and finally, we extend the method to n dimensions.This strategy presents a drawback, the stencils used to interpolate are square and the number of points is larger than necessary to get the optimal order of accuracy.To improve this fact and to apply this technique in partial differential equations contexts are the principal lines of our future work.

r− 1 k=0 ω r k = 1 ,
and p r k are the Lagrange interpolants with nodes S r k = {x i−r+k , • • • , x i+k }.

Figure 1 :
Figure1: Diagram showing the structure of the optimal weights needed to obtain optimal order of accuracy,[3].

Table 4 :
Grid refinement analysis for the classical WENO-8 algorithm for the function in (43).

Table 5 :
Non-uniform grid refinement analysis for the new WENO-6 algorithm for the function(43)

Table 7 :
Non-uniform grid refinement analysis for the new WENO-6 algorithm for the function (45).

Table 8 :
Non-uniform grid refinement analysis for the classical WENO-6 algorithm for the function (45)