A Piecewise Polynomial Harmonic Nonlinear Interpolatory Reconstruction Operator on Non Uniform Grids—Adaptation around Jump Discontinuities and Elimination of Gibbs Phenomenon

: In this paper, we analyze the behavior of a nonlinear reconstruction operator called PPH around discontinuities. The acronym PPH stands for Piecewise Polynomial Harmonic, since it uses piecewise polynomials deﬁned by means of an adaption based on the use of the weighted Harmonic mean. This study is carried out in the general case of nonuniform grids, although for some results we restrict to σ quasi-uniform grids. In particular we analyze the numerical order of approximation close to jump discontinuities and the elimination of the Gibbs effects. We show, both theoretically and with numerical examples, that the numerical order is reduced but not completely lost as it is the case in their linear counterparts. Moreover we observe that the reconstruction is free of any Gibbs effects for sufﬁciently small grid sizes.


Introduction
Due to the extended use of reconstruction operators in many fields of application, ranging from hyperbolic conservation laws to computer aided geometric design, it is of great importance to dispose of efficient methods to build them for different situations. In general, and for the sake of simplicity, the considered functions are polynomials. High degree polynomials are, however, usually avoided because they are known to generate oscillations and undesirable effects.
Linear operators behave improperly in presence of jump discontinuities, so that different nonlinear operators have emerged to deal with this problematic. Recent approaches to deal with similar problems of functions affected by discontinuities can be found for example in [1][2][3][4][5]. And these nonlinear methods also give rise to interesting applications. To mention some of them one can refer to [6][7][8][9][10][11].
In this article we pay attention to one of these operators that was defined in [12] under the name PPH (Piecewise Polynomial Harmonic). This operator can be seen as a nonlinear counterpart of the classical four points piecewise Lagrange interpolation. The theoretical analysis as much as the practical applications were developed in uniform grids in previous articles (see, for example, [12][13][14][15][16][17][18]). In turn these reconstruction operators are the heart of the definition of associated subdivision and multiresolution schemes [5,[19][20][21].
In this paper, we extend the definition of the PPH reconstruction operator to data over non uniform grids and we study some properties of this operator. In particular, we analyze the behavior of the operator in presence of jump discontinuities. We prove adaption to the jump discontinuity in the sense that some order of approximation is maintained in the area close to the discontinuity, on the contrary to what happens with linear operators that lose completely the approximation order. We also prove, as much theoretically as in numerical experiments, the absence of any Gibss phenomena.
The paper is organized as follows-in Section 2 we remind the nonlinear PPH reconstruction operator [22] on nonuniform grids. Section 3 is dedicated to study the adaption of the operator to the presence of jump discontinuities, making some emphasis in the order of approximation. In Section 4 we analyze the behavior of the operator with respect to the Gibbs phenomena. In Section 5 we present some numerical tests. Finally, some conclusions are given in Section 6.

A Nonlinear PPH Reconstruction Operator on Non Uniform Grids
In this section we recall the definition of the nonlinear PPH reconstruction operator on nonuniform grids, see [22]. We include the necessary elements for the rest of the article. In [22] the reconstruction operator is designed to deal with strictly convex functions, albeit it is also of interest in the case of working with piecewise smooth functions affected by isolated jump discontinuities. This will be our case of interest in this section and in the rest of the article.
Let us define a nonuniform grid X = (x i ) i ∈ Z in R. Let us also denote h i := x i − x i−1 , the nonuniform spacing between abscissae. We consider underlying piecewise continuous functions f (x) with at most a finite set of isolated corner or jump discontinuities, and let us call f i := f (x i ) the ordinates corresponding to the point values of the function at the given abscissae. We also introduce the following notations. In first place, the second order divided differences in second place a weighted arithmetic mean of D j and D j+1 defined as with the weights Given these ingredients in [22] we can find the following definitions, and results that we will use later.

Lemma 1. Let us consider the set of ordinates
Then the values f j−1 and f j+2 at the extremes can be expressed as with the constants γ j,i , i = −1, 0, 1, 2 given by Definition 1. Given x, y ∈ R, and w x , w y ∈ R such that w x > 0, w y > 0, and w x + w y = 1, we denote as V the function Lemma 2. If x ≥ 0 and y ≥ 0, the harmonic mean is bounded as follows Next definition, which is commonly used in numerical analysis, is going to be essential through the rest of the article.

Definition 2.
An expression e(h) = O(h r ), r ∈ Z means that there exist h 0 > 0 and M > 0 Lemma 3. Let a > 0 a fixed positive real number, and let x ≥ a and y ≥ a. If |x − y| = O(h), and xy > 0, then the weighted harmonic mean is also close to the weighted arithmetic mean M(x, y) = w x x + w y y, Definition 3 (PPH reconstruction). Let X = (x i ) i∈Z be a nonuniform mesh. Let f = ( f i ) i∈Z a sequence in l ∞ (Z). Let D j and D j+1 be the second order divided differences, and for each j ∈ Z let us consider the modified values { f j−1 , f j , f j+1 , f j+2 } built according to the following rule • Case 2: where γ j,i , i = −1, 0, 1, 2 are given in (5) and V j = V(D j , D j+1 ), with V the weighted harmonic mean defined in (6) with the weights w j,0 and w j,1 in (3). We define the PPH nonlinear reconstruction operator as where PPH j (x) is the unique interpolation polynomial which satisfies According to Definition 3, it is possible to establish a parallelism with Lagrange interpolation, in fact we can write the PPH reconstruction as where the the coefficients a j,i , i = 0, 1, 2, 3 are calculated by imposing conditions (12). We explain each one of the two possible local cases, Case 1 or Case 2. The coefficients will have symmetrical expressions. Case 1. |D j | ≤ |D j+1 |, which means that a potential singularity may lay in [x j+1 , x j+2 ]. It has been proposed to replace f j+2 with f j+2 in Equation (9) by changing the weighted arithmetic mean in Equation (4b) for its corresponding weighted harmonic mean. This replacement has been performed to carry out a witty modification of the value f j+2 in such a way that its difference with respect to the original f j+2 is large in presence of a discontinuity, but remains sufficiently small in smooth areas maintaining the approximation order. Lemma 2 is crucial for the adaption in case of dealing with the presence of a jump discontinuity, while Lemma 3 plays a fundamental part in proving fourth approximation order for smooth areas of an underlying function. In this case the coefficients a j,i , i = 0, 1, 2, 3 of the PPH polynomial read For our purposes, in the next sections we need to examine deeper the relation with Lagrange interpolation. In particular we get that and considering the Lagrange interpolation polynomial written in the same form as in (13), that is we get that the difference of these coefficients with the ones of PPH j (x) is given by Case 2. |D j | > |D j+1 |, which means that a possible singularity lies in [x j−1 , x j ]. In this case, in Definition 3, the value f j−1 is replaced with f j−1 by using expression (10). Similar comments apply in this case due to symmetry considerations. The coefficients for the polynomial (13) now read The expressions relating the coefficients of the PPH polynomial with the Lagrange interpolation polynomial now write In next section, we will study the approximation order of the PPH reconstruction operator in presence of isolated jump discontinuities.

Approximation Order around Jump Discontinuities
We are going to study the approximation order of the given reconstruction for functions of class C 4 (R) with an isolated jump discontinuity at a given point µ. We consider only the case of working with σ quasi-uniform grids, according with the following definition.
In what follows we give a proposition proving full order accuracy for convex regions of the function, that is fourth order accuracy, and observing that the approximation order is reduced to second order close to the singularities and to third order close to inflection points. We would like to focuss especial attention to the intervals around the discontinuity where the order is reduced, but not completely lost.
Theorem 1. Let f (x) be a function of class C 4 (R\{µ}), with a jump discontinuity at the point µ.
], a > 0, a fixed positive real number, Ω the set of all inflexion points of f (x), and d(x, Ω) the distance function defined by Then, the reconstruction PPH(x) satisfies

3.
In Proof. We do the proof point by point. (2) and (6) we can write From hypothesis we have that the initial data are strictly convex in the considered area [x i−1 , x i+2 ] (for a concave function the arguments remain the same) and therefore they satisfy Since second order divided differences amount to second derivatives at an intermediate point divided by two, i.e and from (21) we get that Plugging this information into (17) Thus where p L i (x) is the Lagrange interpolatory polynomial. Taking into account again the triangular inequality using that Lagrange interpolation also attains fourth order accuracy.

2.
We now prove Point 2. Since d( , Ω) < a, and depending on the exact distance to the inflection point we encounter , and the rest of the proof follows exactly the same track as in Point 1, giving the enunciated result.

For proving Point 3, we observe that in this case
, and again following the same track as in previous points we get and therefore in this case the accuracy is reduced to third order.

4.
In order to prove Point 4, let us suppose without lost of generalization that The other case it is proven analogously. Since by hypothesis the function f (x) is smooth in [x j−2 , x j ] , and it presents a jump discontinuity at the interval [x j , The difference between these coefficients and the ones of PPH j−1 (x) shown in Equation (14) is given by At this stage we distinguish two cases: Taking into account Equations (6), (7) and (25) and the triangular inequality we obtain Equations (6) and (25) and the triangular inequality lead us to And these last chains of equations finish the proof.
We observe that close to the jump discontinuity, that is, in the intervals [x j−1 , x j ] and [x j+1 , x j+2 ], we do not lose all accuracy, but we maintain at least second order accuracy. Unfortunately, in the central interval [x j , x j+1 ] containing the singularity this approach does not allow us to obtain any gain with respect to other reconstruction operators.

Remark 1.
Notice that linear reconstruction operators based on an stencil of four points typically lose the approximation order in three intervals around discontinuities, while the introduced nonlinear reconstruction operator only loses completely the aproximation order in the interval containing the jump discontinuity and maintains at least second order accuracy, that is, O(h 2 ), in the adjacent intervals. In the interval containing the jump discontinuity the approximation order is lost also in the nonlinear reconstruction strategy, since with point values of the function it is impossible to detect the exact position of the jump discontinuity. Remark 2. The order reduction due to inflection points can be tackled using a translation strategy in the definition of the Harmonic mean, to avoid arguments of different signs. This strategy complicates the definition of the operator, but it has been satisfactorily introduced on various occasions [12,23]. In practice the translation is needed not only at the interval containing the inflection point, but also in adjacent intervals.

Analysis of Gibbs Phenomena around Jump Discontinuities
In this section we are going to give a result analyzing the behavior of the proposed nonlinear reconstruction with respect to the generation of possible Gibbs effects due to the presence of jump discontinuities in the underlying function. In particular we prove the following proposition.
Before enunciating the theorem we introduce some definitions.
Definition 5. Given X 0 = {x i } i∈Z a σ quasi-uniform grid in R, we define, for k ∈ N (the larger the k the larger the resolution), the set of nested grids given by Let us also denote [ f ] the size of a jump discontinuity, r k j (x) the straight line joining the points (x k j , f k j ) and (x k j+1 , f k j+1 ), d k i (x), i = j − 1, j, j + 1 the vertical distance from the reconstruction PPH k j (x) to the horizontal line passing through the middle point of (x k j , f k j ) and (x k j+1 , f k j+1 ). The respective expressions come given by We will also use r k max as the maximum distance between PPH k j (x) and r k j (x) measured perpendicularly to r k j (x).
Theorem 2. Let X k = {x k i } i∈Z , k ∈ N ∪ {0} be a set of nested σ quasi-uniform grids in R. Let f ∈ C 4 (R) be a function with four continuous derivatives in all the real line with an isolated jump discontinuity at the abscissa µ located at a certain [x k j , x k j+1 ] for each k, where j depends on k. Then, ∃k 0 : ∀k ≥ k 0 the reconstruction PPH k (x) associated to the data f k := ( f (x k i )) i∈Z does not generate Gibbs phenomena. In particular, the following statements hold: Proof. Let us consider k large enough, k ≥ k 0 , such that Then and from (28), (29) and (6) we get We carry out the rest of the proof addressing point after point.

1.
Since only three intervals are affected by the jump discontinuity for construction, then ∀k

2.
We are going to show now that the oscillations due to the presence of the discontinuity diminish at the interval [x k j−1 , x k j ] with k increasing. In [x k j−1 , x k j ] the PPH reconstruction amounts to As |D k j−1 | ≤ |D k j |, the coefficients are given by (14) adapted to the interval j − 1 Taking into account property (7) of the harmonic mean, we can write Considering (31), the distance d k j−1 (x) can be bounded by , applying arguments based on symmetry and taking into account that |D k j+1 | ≥ |D k j+2 | we also get that d k , as V k j = 0 due to (30), the expression of PPH k j (x) according to (13), (14), At this point we consider two subcases depending on |D k j | and |D k j+1 | where The maximum value of the function d k j (x) in the interval [x k j , x k j+1 ] is either at the extremes of the interval or among any possible critical point x c verifying (d k j ) (x c ) = 0. At the extremes of the interval we have d k j (x k j ) = d k j (x k j+1 ) = 1 2 f k j+1 − f k j , and the condition is satisfied. We are going to prove that the local reconstruction For this purpose, we shall prove that (PPH k j ) (x) = 0 ∀x ∈ [x k j , x k j+1 ]. We start computing (PPH k j ) (x) and (PPH k j ) (x), From (37) and (30) we get To analyze the sign of (PPH k j ) (x j ) we replace in (36) D k j by its expression (28) and we consider two subcases depending on the sign of [ f ], , what amounts to say that there is not local maximum value of PPH k j (x) inside the interval.
Following a similar path to case 4.1 we arrive to and therefore PPH k

5.
We start computing the points where the slope of the tangent of PPH k j (x) equals to the slope of the straight line r k j (x). We consider two subcases, 5.1 |D k j | ≤ |D k j+1 | In this case, the above mentioned points where the tangent of p k j (x) is parallel to r k j (x) are given by: .
The largest distance from these points to r k j (x) is the maximum distance between PPH k j (x) and r k j (x) measured perpendicularly to r k j (x). For both points this distance coincides with The required points P 1 and P 2 in this case take the form: and r k max is given by Remark 3. The hypothesis in Theorem 2 concerning the use of a nested set of σ quasi-uniform grids amounts in practice to build the reconstruction with a small enough maximum grid size.
In the next section we carry out some numerical experiments to check that the practical observations coincide with the theoretical results.

Numerical Experiment
In this section we present a simple numerical test to validate the theoretical results. Our experiment computes the approximation order of the considered reconstruction in several areas corresponding with the different points in Theorem 1. In particular we measure the approximation order in the following areas, identified with the given acronyms: A 0 : In the subinterval containing the discontinuity. A 1 : In a region where the function is smooth without inflexion points. A 2 : In a region where the function is smooth but contains a inflexion point. Let X 0 = (0, 3, 8, 11, 17, 23, 25, 27, 31, 32, 36, 37.5, 38, 39.3, 40) π 20 be a non uniform grid in [0, 2π] and f (x) the following smooth function with a jump discontinuity at x = 1.2π, and an inflexion point at x = 3π 2 , Given the initial abscissas x i , i ∈ I = {0, . . . , 14}, we consider the set of nested , and I k = {x k 0 , . . . , x k n k }, with n k = 2n k−1 − 1, n 0 = 14, k = 0, 1, . . . , 7. For each level of resolution k we build the PPH reconstruction using the data (x k i , f (x k i )), i ∈ I k computing the approximation errors in infinity norm with respect to the original function using a denser set of abscissas, that is, we compute a numerical approximation of Then, we compute the numerical approximation order as Notice that due to Theorem 1 we can assume that for fine enough grids In Tables 1 and 2 we present the errors committed by Lagrange and PPH reconstructions respectively when using as initial nodes the defined nested grids X k . The errors appear separately for each kind of region A 0 , A 1 , A 2 , A 3 and A 4 . The largest error comes near the jump discontinuity for Lagrange reconstruction, as it can be observed in the column corresponding with A 0 .  In Table 3 we present the obtained approximation orders for the studied PPH reconstruction and just for the sake of comparison we also add the approximation orders for the classical four points piecewise Lagrange polynomial interpolation. We have computed the approximation order in the specified different regions A 0 , A 1 , A 2 , A 3 and A 4 . More in concrete, in the case of region A 1 we use the interval [2,3] for the x variable, in the case of region A 2 the interval [4,5], and in the case of region A 3 the intervals [x d k +k , 2π], where k indicates the resolution level and the index d k is such that the inflexion point falls into the interval [x d k −1 , x d k ] for each k. We can observe that in the region A 0 both reconstructions are affected by the jump discontinuity and they lose the approximation order due mainly to the subinterval containing the discontinuity. In the region of type A 1 both reconstructions attain fourth order accuracy as expected. In the case A 2 the PPH reconstruction reduces the approximation order to third order due to the presence of the inflexion point. Similarly in the vicinity of the inflexion point, region A 3 , the PPH reconstruction stays between p = 3 and p = 4. In the adjacent intervals to the singularity, case A 4 we clearly observe an improvement with respect to Lagrange interpolation, since we obtain order p = 2 while Lagrange completely loses the approximation order. Notice that the order reduction produced in the regions A 2 and A 3 occurs in very limited areas and it can be corrected using a translation strategy (see [12,23]) that we have not implemented in this experiment with the aim of studying the original reconstruction operator. Table 3. Approximation orders obtained at iteration k, k = 1, .., 7 for the considered cases A 0 , A 1 , A 2 , A 3 and A 4 using the PPH and Lagrange reconstructions. In Figure 1 we plot the function f (x) and the Lagrange and PPH reconstructions obtained from the initial grids X k , k = 0, 1, 2. We can see that around the singularity, Lagrange reconstruction looses the approximation order and the Gibss phenomena appears. In this zone, PPH reconstruction performs in a more proper way, avoiding any Gibbs effects. We can see that no oscillations appear in the PPH reconstruction even for the coarsest grid. These observations can be seen more clearly in Figure 2 where we have plotted a zoom of this region for k = 3 for both reconstruction operators Lagrange and PPH. We also point out that the oscillations due to the jump discontinuity in Lagrange reconstruction do not diminish to zero with the subdivision level. In fact, from k = 2 we have check out that the reconstruction values at the local maxima and minima of the oscillations remain almost constant.
In the jump interval the distance r k max decreases as k increases, since r k max = O(h k ). In Table 4 the values for k = 0, 1, . . . , 7 are shown. We can see that at a certain subdivision level the given values are approximately decreasing with the ratio 1 2 . Therefore, PPH reconstruction approaches to the straight line r k j (x) as k increases.

Conclusions
We have studied the behavior of the PPH reconstruction operator in presence of jump discontinuities for the case of working with σ quasi-uniform grids. For this purpose, the arithmetic and harmonic means used in the uniform case are changed for weighted means with concrete weights, so that the main properties that allow for maintaining order of approximation in smooth areas and adaptation near singularities continue being true.
A explicit result concerning the approximation order, Theorem 1, has been proved, showing at least second order of approximation for the adjacent intervals to the one containing the jump discontinuity, and ensuring fourth order of approximation in convex (concave) parts of the function far from inflexion points. At a interval containing a inflexion point we get third order of approximation and in the vicinity the order grows progressively till fourth order.
A main result of this article is Theorem 2 in Section 4 proving that the presented reconstruction operator does not generate any Gibbs phenomena in the concrete sense indicated in the enunciate for σ quasi-uniform grids where the maximum space between nodes of the grid is small enough.
Finally we have carried out some numerical experiments to reinforce the theoretical results proven as much in Proposition 1 as in Theorem 1.