Non-linear WENO B-spline based approximation method

In this work we present a new WENO b-spline based quasi-interpolation algorithm. The novelty of this construction resides in the application of the WENO weights to the b-spline functions, that are a partition of unity, instead to the coefficients that multiply the b-spline functions of the spline. The result obtained conserves the smoothness of the original spline and presents adaption to discontinuities in the function. Another new idea that we introduce in this work is the use of different base weight functions from those proposed in classical WENO algorithms. Apart from introducing the construction of the new algorithms, we present theoretical results regarding the order of accuracy obtained at smooth zones and close to the discontinuity, as well as theoretical considerations about how to design the new weight functions. Through a tensor product strategy, we extend our results to several dimensions. In order to check the theoretical results obtained, we present an extended battery of numerical experiments in one, two and tree dimensions that support our conclussions.


Introduction
The WENO (weighted essentially non-oscillatory) method is an approximation strategy used for dealing with jump discontinuities that may occur in the solutions of hyperbolic PDEs.The main idea of the WENO method is to build approximations using information from only one side of the discontinuity.This is achieved by utilizing certain smoothness indicators in order to construct a weighted approximation that uses data from smooth zones mainly.We note that in the development of numerical schemes for hyperbolic PDEs, the smoothness of the local approximation is less important and is less addressed.
A very brief and incomplete historical revision about the WENO algorithm can be found in what follows.Liu, Osher and Chan introduced for the first time the WENO algorithm in [22], and it was proposed as a local nonlinear convex combination of different interpolants, covering adjacent stencils.The values of the weights of this combination were functions of some smoothness indicators, calculated over the same stencil of the interpolants.These weights are designed so that the stencils affected by a discontinuity have a negligible contribution to the convex combination.Using a measure inspired by the total variation, Jiang and Shu introduced more suitable and efficient smoothness indicators in [23].The final objective of this technique was to provide high order of accuracy at smooth zones of the function while providing adaption to singularities using the same information as ENO (essentially non-oscillatory), [12,13,14,15,16,17,18,19,20,21].In the literature can be found many references about WENO algorithm.An incomplete list that the reader can visit to obtain more information about the state of the art of this algorithm can be: [15,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42], with a special emphasis in [43,44] and the references therein.
More recently, some investigation has been done about how to design WENO algorithms for splines [1,5].In this paper, we address the problem of approximating functions with jump discontinuities using their data on a uniform grid.To do so, we develop a new class of WENO-like methods which are based upon B-spline quasi-interpolation operators.There are several new ideas as well as new features presented within this development.The main idea is to apply the smoothness-based weights to the B-spline basis functions rather than to their coefficients.The second novelty is the introduction of a new weighting strategy challenging the classical one.The approximation method is defined with splines of any degree, and the resulting approximations are as smooth as the B-splines involved.
The method is first presented for the univariate case, whilst the extension to higher dimensions is straightforward using the tensor product strategy.An analysis of the approximation error is presented, and extensive numerical tests and experiments are performed.Numerical results are also available for the approximation of functions with jump discontinuities in 2D and in 3D.

Preliminaries
In this section, we present the idea to construct a new non-linear operator which avoids the Gibbs-type phenomenon in the adjacent intervals close to a singularity and conserves the order of accuracy in the smooth zones.Using the same notation as [2], we consider a function f , a constant h > 0, and the equidistant data {f n := f (nh)}.For the approximation by a spline of degree p we denote ), and define the quasi-interpolation operator where the different components are: , with equidistant knots S p = − p+1 2 , . . ., p+1
• The linear operator L p , that is defined using the ideas presented in [10] as: with the coefficients c p,j , j = − p 2 , . . ., p 2 : where δ i,j is the Kronecker delta function and t(i, j) are the central factorial numbers of the first kind (see [10] and [6]), which can be computed recursively as: and t(0, 0 The main power of the quasi-interpolator Q p lies in its capacity of obtaining explicit smooth, high-order, approximations to the function, as described in the following propositions (see for example [10]): However, when the data present singularities or strong gradients, the reconstruction presents some artifacts close to the discontinuities (see e.g.[1,5]).This Gibbs-type phenomenon can be avoided by constructing a new non-linear operator, as described in the next section.

Non-linear WENO-based B-spline method
This section aims to obtain a non-linear approximation operator by applying the arguments introduced for the WENO method to the operator Q p .If we analyze Eq. ( 1) for a fixed point x * ∈ R, due to the compact support of the B p function, we can rewrite it as where We call the values C p k (x * ) the B-optimal weights for x * .Note that Eq. ( 4) is a convex combination since k∈J (x * ) C p k (x * ) = 1.We now substitute these linear weights with non-linear ones following the WENO strategy (see, e.g.[23,22]).Thus, being where Ψ is a function that satisfies certain conditions to ensure maximum order at the smooth parts of the data.We will introduce these conditions within the next results.In the WENO algorithm context, the function Ψ has been typically defined as: where ǫ h is a parameter to avoid zero values in the denominator, for example, ǫ h = h 2 , and the parameter t is an integer that assures maximum order of accuracy close to the discontinuities.Finally, the values I k,p are the smoothness indicators (see, e.g.[23,22]) which are designed to detect which stencils are contaminated by a discontinuity.In our case, we choose the following: The properties of these operators are well-known, and they are essential to get optimal order of approximation at the smooth parts of the function, and to avoid the Gibbs phenomenon near singularities.We summarize the principal characteristics in the next proposition.
Proof.The proof is straightforward using Taylor's expansion.
For convenience, in the rest of the paper, we denote I k,p = I k,p (f ).
In order to guarantee that the behavior of the reconstruction obtained using the new operator is similar to the linear one at the smooth zones, we impose some conditions on the function Ψ.In the following proposition, we present the definition of the new non-linear approximation operator, together with the first requirement on the function Ψ. Proposition 3.2 Let h > 0, x * ∈ R and let I k,p be the smoothness indicators of f , with k ∈ J (x * ).Let Ψ be a function such that with K k,l ∈ R independent of h.Then, Proof.We fix a value l ∈ J (x * ), we have that for any k ∈ J (x * ) there exist K k,l ∈ R such that Finally, With these ingredients, we have the following main result concerning the approximation order of the non-linear operator at the smooth regions of f .
], for all k ∈ J (x * ) then: Proof.The proof is straightforward, by Prop.3.2: and from ], for all k ∈ J (x * ), using Proposition 2.2, we get: With the condition given in Proposition 3.3 the accuracy should be optimal, at least, when the data do not present singularities.In the following lemma, we present sufficient conditions for constructing a function Ψ with the required properties.
as h → 0, with s 1 , s ∈ N, s 1 ≥ s and q 1 , q 2 ∈ Z, s 1 + q 2 ≥ q 1 .Then, the function Ψ(x) = 1 ψ(x) satisfies: Proof.As ψ ∈ C 2 then Ψ ∈ C 2 and, by Taylor's expansions, there exists ξ ∈ (min(I l,p , I k,p ), max(I l,p , I k,p )) such that Note that, if we determine Ψ(x) = C = 0, for all x ∈ R we get good behaviour when the data do not present any singularity.This is clear because, in this case, C p k = ω p k , for all k and, hence, Q p = Q NL p .Therefore, we impose different conditions to adapt our new operator to the discontinuities in the following result.

Proposition 3.5 Let us consider
, and let I n,p be smoothness indicators of f .If the function f presents a jump or a kink (a jump in the first order derivative) in at least one interval, n ∈ J (x * ), and, as h → 0, with s ∈ N, q ∈ Z, s > q, then: Proof.By (15), we have that and if f presents a jump or a kink in [h(k − p 2 ), h(k + p 2 )] we get Finally, we prove that under these conditions we have approximation order O(h) close to the singularities.
Proof.We suppose that p is an odd number, the proof is similar for an even p. then From the definition of the operator L p , Eq. ( 2), we have that we have that

This expression holds since, if k
(17) Finally, by ( 16) and ( 17), we have that

The design of the function Ψ
In this subsection, we present some possible options to obtain a function Ψ which satisfies the abovementioned conditions.We know, by Eq. (11), that s 1 = 2p + 1 if p is even, s 1 = 2p − 1 if it is odd, and that Ψ = 1/ψ with ψ ∈ C 2 (R) and ψ(x) = 0, for all x ∈ R + .Therefore, we only summarize the conditions of the function ψ: ].Note that the first three conditions are devoted to obtaining optimal order, so, it is necessary that: The fourth condition is to adapt the operator to deal with discontinuities.In this case, the constraint is q < s.We suppose that if the function is not continuous then the function presents a jump discontinuity in at least one stencil, i.e., I k,p = O(1) as h → 0.
We start with the classical smoothness indicators designed by Jiang and Shu in [23], that where analyzed and modified in [4].Thereby, we consider the function ψ s (x) = h 2 +x.In this case, we have s = 2, q 1 = q 2 = 0 and q = 0. Therefore, all the conditions are satisfied.The behavior of the resulting approximation using the new operator is adequate.
Our second example is the function ψ c (x) = C + x/h.It is clear that s = 1, its derivative is ψ ′ (x) = 1/h then q 1 = −1 and q 2 = 0. Finally, in the presence of a jump discontinuity, we have that q = −1.Again, all the conditions are fulfilled.For this function, we should pre-determine a constant C. If it is very big, the efficiency of the algorithm in the smooth parts increases, however, the capacity of detection of discontinuities is reduced.If C is chosen to be small, the absolute error in the smooth part is higher than the one obtained when C is bigger.We choose C = 1 and we perform some numerical experiments in the next section.We introduce the last example to avoid determining constants.We define the function ψ d (x) = e x/h and calculate the constants: s = 0, q ∈ Z − .A numerical accuracy comparison with these different options of ψ is presented in section 5, where some advantage is observed when using the function ψ d .Also, ψ d seems to provide more sharpness in the approximations close to discontinuities.

Generalization to several dimensions
One of the main advantages of the new method is its ease of generalization to higher dimensions.In this section, we present the general method in k-dimensions.For this, we use the notation presented in [2] (also see [10]).Let h > 0, we consider a real function f : R k → R, the vectors p = (p 1 , . . ., p k ), n = (n 1 , . . ., n k ) and the values We denote as: ).Thus, the quasi-linear operator in n dimensions for a point x * = (x * 1 , . . ., x * k ) ∈ R k is: where L p is the tensor product of L p , Equation (2), i.e.
with c p,j defined in Equation ( 3), and Using the same notation as in Eq. ( 5), we have that: To generalize the operator Q NL p , Eq. ( 6), we replace each value C p l n l (x * l ) by its corresponding ω p l n l (x * l ) in the same way as we have developed above.If we denote: the new non-linear operator is: As we can see, a tensor product strategy has been used.Therefore, if the function is sufficiently smooth the order of accuracy is 1 + min 1≤l≤k {p l }, and O(h) when the function presents some singularity.

Numerical experiments
In this section, we analyze numerically some of the properties of the algorithms proposed in previous sections.More specifically, we check the behavior of the proposed WENO B-spline-based algorithms at smooth zones, and close to the discontinuities.We are interested in studying the numerical order and the absence of the Gibbs-type phenomenon.We also present results for bivariate and trivariate data, obtained by applying a tensor product strategy, as explained in Section 4.

Accuracy at smooth zones
In this section, we analyze the behavior of the new algorithm and the classical one at regions of smoothness in order to check the numerical accuracy reached.To do so, we consider the smooth function and perform a grid refinement analysis.The error in the approximated solution is denoted by E l and it is computed using a grid size h l using the infinity norm inside the considered interval, We can see that l represents the step of the grid refinement analysis, which consists in reducing the grid size h to h/2 when going from l to l + 1.The order of accuracy is obtained in general as, The results for the quadratic spline are presented in table 2, for the cubic in table 3, for the quartic in table 4 and for the quintic in table 5.All the tables show the results for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1, and the classical spline.We can see how all the versions of the algorithm attain the correct numerical accuracy when applied to a smooth function: order O(h 3 ) for the quadratic spline, order O(h 4 ) for the cubic, and so on.

Accuracy close to the discontinuity
It is also interesting to check the numerical accuracy of the new method close to the discontinuity, as stated in Preposition 3.6.We consider the function in ( 25) Table 2: Grid refinement analysis for the accuracy of the quadratic WENO spline and the classical quadratic spline using the infinity norm.The original data has been sampled from the function in (22).We show the results for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1, and the classical spline.The low-resolution nodes have been sampled with m = 2 l nodes and the high-resolution data with 11(m − 1) + m points.Table 3: Grid refinement analysis for the accuracy of the cubic WENO spline and the classical cubic spline using the infinity norm.The original data has been sampled from the function in (22).We show the results for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1, and the classical spline.The low-resolution nodes have been sampled with m = 2 l nodes and the high-resolution data with 10(m − 1) + m points.Table 4: Grid refinement analysis for the accuracy of the quartic WENO spline and the classical quartic spline using the infinity norm.The original data has been sampled from the function in (22).We show the results for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1, and the classical spline.The low-resolution nodes have been sampled with m = 2 l nodes and the high-resolution data with 11(m − 1) + m points.
We check the accuracy in (0.5, 1] but avoiding the interval that contains the discontinuity close to x = 0.5. The reason is that we know that the approximation in the interval that contains the discontinuity presents the order of accuracy is O(h 0 ).To do so, as before, we use the infinity norm of the error in (23) and the formula for the numerical order in (24).The results for the quadratic spline are presented in table 6, for the cubic in table 7, for the quartic in table 8 and for the quintic in table 9.As before, all the tables show the results for the functions ψ s , ψ c , ψ d .All of them show order O(h) in the infinity norm, as proved in Preposition 3.6.This order of approximation corresponds to the intervals that are adjacent to the one that contains the discontinuity.
2 10 Error WENO-ψ s quintic spline (E l ) 3.0918e-04 9.0968e-07 3.0077e-09 9.5293e-12 3.9302e-14 Table 5: Grid refinement analysis for the accuracy of the quintic WENO spline and the classical quintic spline using the infinity norm.The original data has been sampled from the function in (22).We show the results for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1, and the classical spline.The low-resolution nodes have been sampled with m = 2 l nodes and the high-resolution data with 10(m − 1) + m points.(25).In this case, we present the results of the quartic WENO algorithm for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1.

Adaption to the presence of discontinuities: Gibbs phenomenon
In this subsection, we analyze the behavior of the algorithms in terms of their adaption close the discontinuity.
Figure 1 shows the results of approximating the function in (25) close to the discontinuity.To obtain these results, we have started from a sampling of the function in (25)  Table 9: Grid refinement analysis in the infinity norm close to the discontinuity of the function shown in (25).In this case, we present the results of the quintic WENO algorithm for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1.
quartic spline, or 10 intermediate points for the cubic and quintic spline.As before, we present the results for the quadratic, cubic, quartic, and quintic WENO B-spline-based algorithm the functions ψ s , ψ c , ψ d , presented in Subsection 3.1, and the respective classical splines.It is easy to appreciate how the three versions of the WENO B-spline-based algorithm presented managed to reduce the oscillation close to the discontinuity that the classical spline presents.We can also see that each function ψ s , ψ c or ψ d allows the WENO algorithm to present different degrees of sharpness close to the discontinuity, ψ d being the sharpest one.For the cubic spline, we can see that the three functions ψ s , ψ c , and ψ d behave very similarly close to the discontinuity.For the quintic spline ψ s and ψ c also behave very similarly.It is important to mention that at the interval that contains the discontinuity all the algorithms present smearing, so it is not possible to control the error.

Experiments for cubic splines and bi-variate data using tensor product
In this section, we present the results that we obtain when extending the one-dimensional algorithms presented in previous sections to several dimensions.The process is described in Section 4. Figure 2 shows the results obtained when approximating the function f (x, y) = cos(xy), (x − 0.5) 2 + (y − 0.5) 2 ≤ r 2 , sin(xy), (x − 0.5) 2 + (y − 0.5 with r = 1 4 .The results are shown in Figure 26.We can see how the nonlinear algorithms manage to reduce the oscillations of the classical approach close to the discontinuity.

Experiments for cubic splines and three-variate data using tensor product
In this section, we present similar results to the ones obtained in the previous subsection, but applied to trivariate data.We have used again the tensor product approach described in Section 4. Figure 3 shows the results obtained when approximating the function with r = 0.4.In this case, we have chosen the function for the weights ψ d , presented in Subsection 3.1.
To obtain the results shown in Figure 27, the function in (27)   The column to the right presents the error obtained by each of them.The first row presents the result of the quadratic algorithm, the second presents the results of the cubic, and so on.
in the 3D reconstruction of the linear algorithm can be inferred from the high-frequency oscillations in the contour plots close to the discontinuity surface and from the presence of many small closed contours also close to the discontinuity surface.

Conclusions
In this article, we have introduced a new WENO B-spline-based algorithm, where the adaption properties are reached through the nonlinear modification of the b-spline functions, which are known to be a partition of unity.These B-splines are used to construct nonlinear weight functions, that must satisfy certain properties in order to assure that the spline can avoid Gibbs-like oscillations close to discontinuities in the function.This approach is completely new and, up to our knowledge, has not been proposed in the literature before.Apart from the adaption to the discontinuities, the new algorithm conserves the smoothness of the classical spline.We have provided theoretical proofs for the accuracy of this new construction at smooth zones and close to the discontinuities, as well as conditions for the weight functions to conserve the accuracy of the data at smooth zones.Through a tensor product strategy, we are capable to extend our results to data of any number of dimensions.All the numerical experiments that we have presented satisfy the theoretical conclusions that we have reached.

Proposition 2 . 1
Consider a p degree B-spline, B p , with knots S p , then the local operator Q p defined in Equation (1) is p − 1 continuous.Proposition 2.2 Consider a p degree B-spline, B p , with knots S p , then the local operator Q p defined in Equation (1) reproduces Π p (R) and achieves O(h p+1 ) approximation order to C p+1 functions.

Figure 1 :
Figure1: The column to the left presents the limit function obtained for ψ s , ψ c , ψ d and the linear algorithm.The column to the right presents the error obtained by each of them.The first row presents the result of the quadratic algorithm, the second presents the results of the cubic, and so on.

Figure 2 :
Figure 2: Left column, approximation of the function in (26) obtained for ψ s (first row), ψ c (second row), ψ d (third row) and the linear algorithm (fourth row).Right column, the distribution of the error obtained by each of them.

Figure 3 :
Figure 3: Left column, approximation of volumetric data obtained from the trivariate function in (27).The approximation shown has been obtained for ψ d (first row), and the linear algorithm (second row).Right column, the same figure but projected over the plane xy.

Table 1 :
Table 1 we show some values of c p,j .

Table 6 :
(25) refinement analysis in the infinity norm close to the discontinuity of the function shown in(25).In this case, we present the results of the quadratic WENO algorithm for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1.

Table 7 :
(25) refinement analysis in the infinity norm close to the discontinuity of the function shown in(25).In this case, we present the results of the cubic WENO algorithm for the functions ψ s , ψ c , ψ d , presented in Subsection 3.1.

Table 8 :
Grid refinement analysis in the infinity norm close to the discontinuity of the function shown in of 400 initial points.Then we obtained approximations at the position of these initial points, plus 11 intermediate points for the quadratic and has been sampled with 200 × 200 × 200 points.The reconstruction has been obtained at the position of the initial values plus at two equidistant intermediate values between the original ones in each direction.In this case, we present the approximation of the trivariate function in (27) using a contour slice routine: we have obtained slices of the volumetric data in the z direction and we have represented contour plots of the bivariate data obtained on each slice.We can see how the nonlinear algorithm manages to reduce the oscillations of the classical approach close to the discontinuity, that in this case is a surface, obtaining a sharp reconstruction close to it.The oscillations