A class of C2 quasi-interpolating splines free of Gibbs phenomenon

In many applications, it is useful to use piecewise polynomials that satisfy certain regularity conditions at the joint points. Cubic spline functions emerge as good candidates having C2 regularity. On the other hand, if the data points present discontinuities, the classical spline approximations produce Gibbs oscillations. In a recent paper, we have introduced a new nonlinear spline approximation avoiding the presence of these oscillations. Unfortunately, this new reconstruction loses the C2 regularity. This paper introduces a new nonlinear spline that preserves the regularity at all the joint points except at the end points of an interval containing a discontinuity, and that avoids the Gibbs oscillations.


Introduction
A spline can be defined as any function constructed from one or more polynomial pieces that are joined together satisfying given differentiability requirements. There exist different types of splines (see, e.g., [9]) which can be classified depending on the restrictions introduced (see, e.g., [21]). C 2 -splines are used in the industry in very different kind of applications. The most familiar is computer-aided design, where splines are used to represent geometric entities (see, e.g., [7,8,12,23]) or image processing (see, e.g., [10,18,27]). Another example is computer-aided manufacturing, where splines are involved in the modeling of geometric parts or used in the representation paths for machine tool cutters (see, e.g., [24,28]). They are also used in the reconstruction of field data in many areas of engineering. Other kinds of splines are used to create approximations to the solution of ordinary differential equations, partial differential equations, integral equations, etc. using the finite element method (see, e.g., [6,7,17,22]).
However, when the data contains real or numerical discontinuities, classical interpolation algorithms are not satisfactory. Recently, some new techniques have been developed in order to recover discontinuous functions, for example using Radial Bases Functions (RBF), see [5,20]. In these cases, standard splines suffer from Gibbs oscillations, this phenomenon is first detected in the context of truncated Fourier expansions, [13], and it has been widely analyzed (see, e.g., [11,14,19]). In our paper [3], we introduced a nonlinear procedure that modifies a given spline in order to avoid these oscillations. However, this recent reconstruction loses the C 2 regularity. The goal of the present work is to present a new spline combining both properties: the C 2 regularity and the adaptation to discontinuities.
This article is organized as follows: Section 2 describes the framework that we will use, explains how classical splines are obtained and how to adapt them to the presence of discontinuities. Section 3 introduces and analyzes the properties of a new nonlinear weighted mean that is used to construct a new class of splines adapted to the presence of jump discontinuities having C 2 regularity. Section 4 studies the approximation properties of the new class of splines close to discontinuities and a specific example is considered in Section 5. In Section 6, we analyze how the new spline avoids the Gibbs oscillations. Section 7 presents some numerical experiments where it is shown how the new nonlinear splines perform close to jump discontinuities with uniform grid spacing. Finally, Section 8 presents some conclusions.

Classical cubic splines versus nonlinear cubic splines
In this section, we recall the classical cubic interpolation and analyze the nonlinear modification introduced by Amat, Ruiz, Shu and Trillo in [3]. We will prove that this new type of spline is not necessarily C 2 , i.e., the continuity is sacrificed to avoid the presence of Gibbs oscillations. In order to obtain C 2 regularity adapted to discontinuities, we slightly change the points of interpolation, i.e., we construct a quasi-interpolatory spline approximations using the same technique developed in [3].
Following the notation described in [3], let us consider the set of piecewise continuous functions in the interval [a, b], and let X be a non-uniform partition of the interval [a, b] into m subintervals, Then, we can define the point-values discretization as the values of a function at these points, y i = y(x i ), i = 0, · · · , m.
The spline g(x) in the interval [a, b] depends of the m + 1 pairs of values (x i , y i ), i = 0, · · · , m and in each subinterval [x i , x i+1 ] it is a cubic polynomial g i (x). More precisely, subjected to the following conditions (2) and g 0 (x 0 ) = y 0 , g m−1 (x m ) = y m . This way, a cubic spline that is at least C 2 is obtained. If we solve the previous equations, we get that, where D i stands for the first derivative at every interior knot. We obtain a system of equations for the D i s replacing the conditions for the spline (1) in (2), If we take the natural boundary conditions, i.e., g 0 (x 0 ) = 0 and g m−1 (x m ) = 0, from (1), the left boundary condition implies the equation, and for the right boundary we obtain, Denoting we get the system: with . . , m. Denoting: the previous system can be also written as AD = 3f, with · · · · · · · · · · · · · · · · · · · · · 0 0 0 0 · · · 2 α m 0 0 0 0 · · · 1 2 where M α,β (x, y) is a weighted mean depending on two constants 0 < α, β < 1 such that α + β = 1 and are functions M α,β : R × R → R defined as M α,β (x, y) = αx + βy.
As mentioned in the introduction, it is not difficult to check that linear C 2 splines produce Gibbs phenomenon at zones close to the discontinuities. For this reason, in [3], the authors construct a new nonlinear spline that we present in the following subsection.
Remark 2.1 Throughout, the paper we define h = max(h) m i=1 , and we consider approximation orders in powers of h as h 0. By the notation e = O(h α ) we mean |e| ≤ Ch α where C is independent of h.

A class of nonlinear cubic splines
In order to design the nonlinear spline, we defined a nonlinear weighted mean. Thus, we present the following definition.
Definition 1 (Non-linear weighted mean) Let 0 < α, β < 1, α + β = 1 be two constants and R + = {x ∈ R : x > 0}. A nonlinear weighted mean,H + α,β : R + × R + → R + , will be a function which satisfies the following five properties: 5) Adaptation to discontinuities (AD) property. There exists p ≥ 1 such that for all x, y ∈ R: We extend this definition for all (x, y) ∈ R 2 as: However, we will lose order of accuracy at critical points, i.e., when In order to solve this problem Amat et al. propose in [3] the following nonlinear weighted mean: where T ε : R 2 → R is a function defined as: being ε(x, y) > 0 a function defined for all (x, y) ∈ R 2 . This function depends on the data at each interval. For simplicity, in the rest of the paper, we denote H α,β = H T ε α,β . If we choose a convenient function ε(x, y) (see Section 5 or [3]), it is possible to prove that the (C), (OA) and (AD) properties are satisfied by H α,β for all the points in the plane R 2 .
Using the same reasoning showed in [3], denoting as H α i ,β i (δ i , δ i−1 ) = H i , i = 2, ..., m, we redefine the system (8) changing the vector f bỹ Thus, the new system is AD = 3f where A is the matrix defined in (8). Then, the coefficients of this nonlinear spline are given by: By a direct computation it is easy to check thatg i−1 (x i ) =g i (x i ) =D i , so the nonlinear spline has at least C 1 -regularity. However, the valuesD = (D i ) m i=0 obtained solving this system are not necessarily equal to the values D obtained in (8). This fact causes the spline to lose its C 2 regularity because the new approximation to the derivatives, (D i ), does not satisfy (4). In particular, in [3], this system is obtained for a specific nonlinear mean.
In order to generate a spline that is C 2 continuous, we propose to change the values that we use to interpolate. Thus, the resulting approximation is not interpolatory.

Construction of C 2 -cubic splines adapted to discontinuities
The idea of this new type of spline is to obtain new data (ŷ i ) m i=0 such that the associated linear system (7) agrees with the nonlinear splines defined using the modified right hand side in (12). Thus, by construction the spline will have C 2 regularity at the modified points. In order to obtain values close to the original ones, we cannot introduce modifications at the slopes defined by the points whose intervals contain discontinuities. In particular, at these points we cannot expect C 2 regularity.
We define the new values,ŷ i , implicitly bŷ such that their weighted arithmetic means coincide with the weighted nonlinear means used at the right hand side of the system (12). Thus, we define: where being 0 < η < 1 a constant. The function G is defined taking into account that in smooth regions |δ i − δ i−1 | = O(h) and close to the discontinuity |δ i − δ i−1 | = O (1). With this modification, the slope generated by the pairs of values defining an interval that contains a discontinuity is not modified and the other values remain at a distance of O(h p ) from the original data values. In all the paper, we assume that the parameter h is sufficiently small.
The associated spline will be C 2 except between the pairs of values defining an interval that contains a discontinuity, but the spline becomes only quasi-interpolatory. In any case, we will prove that the order of approximation is the same in both nonlinear schemes.
Using these results, if we take j 0 = 1 and We do not change the values δ i at the interval where there exists a discontinuity since (15). Then, by (15), we get: Yet, the values y i are changed at an interval which is close to a discontinuity since: Analogously, for i 0 + 1, as then: by (19). We have the hypotheses of Corollary 3.2 with j 0 = i 0 + 2 and j 1 = m With these new values, it is clear thatδ m is not necessarily equal to δ m . In order to obtain exactly the system (12), we will suppose as boundary condition that satisfieŝ We define then the piecewise splineĝ(x) bŷ whereD i are the values obtained solving the system (12) except at the interval where the discontinuity is placed.
We can summarize this result in the following theorem. (15) being H α,β (x, y) a nonlinear weighted mean satisfying the (OA) property with p ≥ 1. If we define:

y, z, h) is the function introduced in
then, the associated linear system used to calculateD i , which satisfies the conditions in (2) witĥ Finally, mixing the conditions of the linear and nonlinear splines, (13) and (22), we define a new approximating spline, u(x) :=g(x), as: (25) whereD i are the values that approximate the derivatives obtained solving the system (12) and theŷ i are defined using (14). We can prove the following theorem.
Theorem 3.4 Using identical notation as before, let 2 < i 0 < m − 1 be an integer such that there exists a discontinuity of (24) and (25) is Proof Let 2 ≤ i < i 0 , theD i are the solutions of the system (12), and β i = λ i−1 λ i +λ i−1 , then we get: Using this relation we can see that:ĝ Analogously

On the order of accuracy
Let us start with the special case when the function is sufficiently smooth. In this case, we expect to obtain similar approximations as the original linear spline, i.e., that the distance between the new spline and the original linear spline is of the order introduced by the (OA) property, (9). We prove two auxiliary lemmas that we will use throughout the section. They consist of proving that having an estimate of the order of the approximation for the derivative values, we can estimate the order of approximation of the associated spline.
be the linear spline defined by (1) and (3), beingg(x) defined by (13) and u(x) the new spline defined by (24) and (25) with H α,β (x, y) a nonlinear weighted mean which satisfies the (OA) property with p ≥ 1 in both cases. Letŷ i be defined by (14), and letD i =g (x i ). If there exists j 0 such that then Proof By (3), (13), and (28), we get Analogously, by (3), (25), and (28), we have be the linear spline defined by (1) with (3),g(x) defined by (13) and u(x) the new spline defined by (25) with H α,β (x, y) a nonlinear weighted mean which satisfies the (OA) property with p ≥ 1 in both cases. If there exists j 0 such that and q ≥ 1 such that |y( Proof The proof is straightforward using Lemma 4. p+1,q) ). Now, with these lemmas, we can calculate the order of approximation.

Proposition 4.3 Let (x i , y i ) m i=0 be a set of points with
i=0 < 1 and let g(x) be the linear cubic spline defined by (1) with (3),g(x) defined by (13) and u(x) the new spline defined by (25) with H α,β (x, y) a nonlinear weighted mean which satisfies the (OA) property with p ≥ 1 in both cases, then for all x ∈ [a, b]: Proof The proof follows by verifying the conditions of Lemma 4.1 for all i = 0, . . . , m − 1. Using systems (8) and (12), we obtain: In the following corollary, we prove that if the data are values that come from the discretization of a smooth function by point values, then the new spline has the same order of accuracy as the linear one. Therefore, at the smooth zones, since the classical cubic spline is of order q = 4, if we take a nonlinear average satisfying the (OA) property with p = 3, then the new spline will reach maximum order of accuracy. Now, we will study the case with jump discontinuities. For this case, we calculate the distance between D andD. We bound the norm of the inverse matrix, A −1 using the following lemma proved in [15] and [16] and used, for example, in [4].
Lemma 4.5 [15,16] If the elements of A −1 are denoted by (a −1 j,k ) m j,k=0 then If there exists a discontinuity at the interval [ thus, if j = 0, . . . , m, we can calculate that since m−1 k=1 2 3 ·γ −|j −k| ≤ 2, for all j = 1, . . . , m−1. In order to obtain γ −|j −i 0 +1| ≤ h p , we have to impose that supposing that h < 1. We can summarize this result in Proposition 4.6. (36) depending on values p, h and γ , the number of points m, and the location of the discontinuity i 0 have to satisfy:

Remark 4.1 It is clear that by
Thus, we assume that the previous expression is satisfied in the following results.

Proposition 4.6 Let
, let AD = f and AD =f be the systems defined in (8) and (12) respectively with H α,β (x, y) satisfying the (OA) property with p ≥ 1, then: Proof We obtain the result because if then γ −|j −i 0 +1| ≤ h p and by (35).
Finally, we prove that at intervals placed sufficiently far away from the discontinuities, the order of approximation obtained by the nonlinear splines is maximal.

Corollary 4.7 Let
Proof The proof directly follows since by Lemma 4.1 and Prop. 4.6, we have

On a specific example of a nonlinear weighted mean
In order to define a particular nonlinear weighted mean that satisfies the (C), (OA) and (AD) properties for any p, we need two principal components: a new limiter and a nonlinear translation. We take the new mean introduced in [3]. Thus, let 0 < α, β < 1 be two constants with α + β = 1, and p ≥ 1, we definẽ (37) These means are a generalization of the means introduced in [25] and reduce to them for uniform grid spacing. The new limiter satisfies the properties presented in the following proposition. The proof is trivial and is not included for brevity.
In order to obtain the maximal approximation order, as we have seen in Section 4, we choose p = 3, i.e.,H W 3 . In particular, if we are working with a uniform grid spacing, h i = h i+1 for all i = 1, . . . , m − 1 then α i = β i = 1/2 and the nonlinear weighted means introduced in the previous Subsection transform into de p-power means introduced in [25]. The p-power mean has the general expression, The most important properties that make the power means appropriate are (see [1,25] for more details):

|H p (x, y)| ≤ max (|x|, |y|).
Finally, to define T ε we only have to determine the function ε : R 2 → R + . In [3], it is suggested a nonlinear ε(x, y) based on the smoothness indicators proposed in [26], that is defined as: where I S j is a smoothness indicator that satisfies that is O(1) close to a jump discontinuity and O(h 2 ) at smooth parts. It is clear that at a jump discontinuity ε j is O(h 4 ) and if |x| < |y|, we get then the translation does not affect the adaption attained at jump discontinuities unless a change of sign in the first derivative is placed exactly at the discontinuity but it is adapted to the presence of critical points in the sense that is O(1) at critical points (and around them) and that goes to zero away from them. In order to measure the smoothness of our function, we will use the smoothness indicator designed by Jiang and Shu in [26]. The proposed expression for a non-uniform grid spacing is, where r is the degree of the polynomial p j (x). We have chosen the smoothness indicators in the point values for r = 2 that result from (42). They have the following expression for a uniform grid-spacing, If we choose a smoothness indicator with a wider stencil the results are similar. With this choice, it is proved that smoothness indicators are O(1) at a jump discontinuity and O(h 2 ) at smooth parts (see [2] for more details).

Analysis of the Gibbs phenomenon for nonlinear cubic splines
Let us consider the elimination of the Gibbs phenomenon for jump discontinuities. First of all, let's remember the definition of the Gibbs phenomenon introduced by D. Gottlieb and C.W. Shu in [11]. Given a punctually discontinuous function f and its sampling f h defined by f h n = f (nh), the Gibbs phenomenon deals with the convergence of the approximation g h based on f h towards f when h goes to 0. It can be characterized by two features [11]: 1. Away from the discontinuity the convergence is rather slow and for any point x, 2. There is an overshoot, close to the discontinuity, that does not diminish with the reduction of h. Thus, max |f (x) − g h (x)| does not tend to zero with h.
Following the ideas introduced in [3] it is easy to check that the new nonlinear spline does not suffer from the Gibbs phenomenon. Theorem 6.1 The nonlinear spline obtained through (24) and (25) does not present Gibbs oscillations.
Proof Let us analyze the right hand side of the system in (12). If we remember now that δ i and δ i−1 are divided differences, we know that they are O(1/h) in the presence of a jump discontinuity. Taking into account property 5 of Proposition 5.2, we know that the right hand side of (12) will be, Then, the vector of derivativesD that results from solving (12) will be, This means that the perturbation introduced by the nonlinear spline is O(h) except at the interval that contains the discontinuity. Now we can try to prove that the nonlinear spline provides a prediction that is in the interval [ŷ i ,ŷ i+1 ] when h goes to zero. In order to do so, we express the equation of the spline (1), using the change of variables s = x−x i h , aŝ where the first two terms of (46) amount to a dilation and a translation of the element s 2 (3 − 2s), that has a minimum at s = 0 and a maximum at s = 1, so it can not introduce Gibbs phenomenon. Let us analyze the third and fourth terms of (46). We can see that s 3 − 2s 2 + s and s 3 − s 2 are oscillating functions so the apparition of Gibbs phenomenon can be explained due to the presence of large coefficients accompanying these two elements. In the case of the nonlinear spline in (12) we have already analyzed thatD i = O (1). If this is the case, the two last terms in (46) go to zero when h → 0. Thus, there can not be Gibbs oscillation.

Numerical experiments
The splines proposed are a modification of those introduced in [3] in order to maintain a regularity C 2 . In this numerical section, we are interested in checking that the new splines own the good properties of the splines presented in [3] concerning order of approximation and elimination of Gibbs-like oscillations. Moreover, we are interested in proving that the algorithm is robust with respect to the presence of perturbations in the data, for this reason we will present experiments where we perturb the functions with additive white Gaussian noise. In addition, we would like to confirm that the intervals affected by a discontinuity verify Proposition 4.6. In particular, that we always maintain at least order 1, which is not verified in the case of linear splines causing Gibbs-like oscillations, and that the order increases as we move away from the discontinuity until we recover the maximum order of accuracy.
We perform some numerical tests in order to validate the theoretical results. In particular, we divide this section in three parts: Firstly, we analyze the distance obtained when we apply the function G, (15) with η = 1/8 (we take this constant for all the examples) on the points {y i } m i=0 to get {ŷ i } m i=0 , system given in (14). We will use a uniform grid depending on the parameter h and we will approximate the order obtained. In the second part of this section, we will analyze the Gibbs phenomenon and we will prove that the resulting spline has similar behaviour as that of the splines designed in [3]. Finally, in third part, we study the order of accuracy checking the theoretical results obtained in Section 4. We use the following two functions: and presented in Fig. 1.

Distance between the original and the new data
For our first experiment, we discretize the functions f (x) and l(x) in the interval [−1, 1] using 2 k points, h k = 2 2 k −1 to get the data (y k i ) 2 k −1 i=0 and we apply the algorithm introduced in (14) to obtain (ŷ k i ) 2 k −1 i=0 . We compute ∞ -norm to the left of the discontinuity, i.e., and the numerical approximation order which is: As we can see in Table 1, the numerical order of accuracy converges to 3 when h k decreases.

Numerical experiments using a uniform-grid
In this subsection, we reproduce the experiments presented in [3] but introducing the modification proposed in previous sections. We start from a sampling of the functions in (47) and (48) with 8192 points in [−1, 1]. We subsample the original data to 512 points and subdivide this data to obtain 8192 samples using the proposed spline (modifying the low resolution original data as a preprocessing through (14) and (15)) and the linear spline (7) (without modifying the low resolution original data) over data with and without noise. Then we compare the result with the original data at the highest resolution. In Figs. 2 and 3, the original subsampled data has been represented with red filled circles, the original high resolution data with blue crosses and the reconstruction with black points. All the numerical experiments presented use a   (7) and a zoom around the discontinuity. We can observe the Gibbs phenomenon that appears close to the discontinuity. At the bottom, we can see the nonlinear reconstruction but in this case we have added to the low resolution data white gaussian noise of amplitude 0.5 uniform grid spacing. The results for the non linear spline have been done using the nonlinear translated mean (10) with the translation defined in (11) and the defined in (40). For the experiments with noise, we have chosen additive white gaussian noise of amplitude 0.1 for the function in (48) and 0.5 for the function in (47) (Fig. 3).
We can see that the results are similar to the ones obtained in [3] from a quantitative point of view, but in this case the new spline has been proved theoretically to Fig. 3 Top left, function in (48) reconstructed with an adapted cubic spline. Top right, a zoom around the discontinuity. The low resolution discretization (red filled circles) has 512 points and the reconstruction (black points) 8192 points.The original high resolution data has been represented with blue crosses. At the middle row, we can see the reconstruction obtained through the linear spline (7) and a zoom around the discontinuity. We can observe the Gibbs phenomenon that appears close to the discontinuity. At the bottom, we can see the nonlinear reconstruction but in this case we have added to the low resolution data white gaussian noise of amplitude 0.1 maintain C 2 regularity, an important qualitative property for applications as we mentioned in Section 1.

Accuracy close to the discontinuity
In this subsection, we check the spatial distribution of the accuracy attained by the spline close to the discontinuity. In order to do so, we measure the error in the infinity  norm between the nodes of the splines obtained for the function in (48) in an experiment of subdivision similar to the ones presented in Section 7.2. We will start from a low resolution data of m = 513 nodes, i.e., h = 2 512 and then construct a linear and a nonlinear spline, computed at 15 interior points between the original nodes. Then we divide h by two in order to perform a grid refinement analysis. In particular, we center our attention at the 17 intervals (between the original nodes) to the left of x = 0 that are closer to the discontinuity, we denote the error at the interval (−nh, −(n − 1)h) as e −nh , n = 1, . . . , 17. Dividing the mesh size by two and applying the order formula, (49), we obtain an approximation of the order of accuracy attained in each interval denoting the orders by o −nh . Since [ is the interval where the discontinuity is contained, then Therefore, we can use Prop. 4.6 to calculate the theoretical interval where the order p + 1 is recovered, which is: where γ = 2 + √ 3 because in this example the grid is uniform. In Table 2, we show these constants for every h.
It is interesting to notice that the theoretical values are slightly worse than the numerical ones. In Table 2 we can see that for h =1.95e-03 the value −n 0 = −14.21 and in Table 3, order of accuracy p = 4 is recovered in the interval [−8h, −7h], (o −8h = 4.09), i.e., there exists an error of 6 intervals where the optimal order is reached. Analogously, for all h, the optimal order 4 is reached in previous intervals (marked in bold in Table 3) to the theoretical interval achieved using Prop. 4.6 ( Table 2, column 2). This phenomenon occurs because in Lemma 4.5 and in (35), we use some manipulations in the inequalities that may be unnecessary in some cases. Also, in the results shown in Fig. 4, we can see how the order grows as we get away from the discontinuity.  Table 4 In this table, we present a grid refinement analysis for the approximation order of the linear spline in the intervals from (−5h,  Table 5 In this table, we present a grid refinement analysis for the approximation order of the linear spline in the intervals from (−17h, Finally, we compare the order of accuracy obtained using linear splines. Tables 4  and 5 show the order of accuracy obtained for the linear splines. Comparing with the results shown in Tables 2 and 3, we can observe that in the nonlinear case, at least order 1 is maintained in all the intervals.

Conclusions
In this paper, we have included a new nonlinear spline avoiding the Gibbs oscillations near discontinuities and preserving the C 2 regularity crucial for some applications. Some theoretical results related to the order of approximation and to the elimination of the Gibbs oscillations have been proved. The numerical results presented confirm the good properties of this new spline.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was funded by the Programa de Apoyo a la investigación de la fundación Séneca-Agencia de Ciencia y Tecnología de la Región de Murcia 20928/PI/18, by the national research project MMTM2015-64382-P and PID2019-108336GB-I00 (MINECO/FEDER), by grant MTM2017-83942 funded by Spanish MINECO and by grant PID2020-117211GB-I00 funded by MCIN/AEI/10.13039/501100011033. Data availability Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Conflict of interest The authors declare no competing interests.
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. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.