On the Convexity Preservation of a Quasi C 3 Nonlinear Interpolatory Reconstruction Operator on σ Quasi-Uniform Grids

: This paper is devoted to introducing a nonlinear reconstruction operator, the piecewise polynomial harmonic (PPH), on nonuniform grids. We deﬁne this operator and we study its main properties, such as its reproduction of second-degree polynomials, approximation order, and conditions for convexity preservation. In particular, for σ quasi-uniform grids with σ ≤ 4, we get a quasi C 3 reconstruction that maintains the convexity properties of the initial data. We give some numerical experiments regarding the approximation order and the convexity


Introduction
Reconstruction operators are widely used in computer-aided geometric design. For simplicity, the functions that are typically used as operators are polynomials. In order to avoid undesirable phenomena generated by high-degree polynomials, reconstructions are usually built piecewise. Due to the bad behavior of linear operators in the presence of discontinuities, it has become necessary to design nonlinear operators to overcome this drawback. One of these operators was defined in [1] and was called the piecewise polynomial harmonic (PPH). This operator essentially consists of a clever modification of the classical four-point piecewise Lagrange interpolation. The initial purpose of this definition was to deal with discontinuities, reducing the undesirable effects to only one interval instead of the three intervals affected in a reconstruction built with a four-point stencil. In addition to that, as we will see throughout this paper, the reconstruction may also play an important role in design purposes, since it keeps the convexity properties of the given starting data.
For the sake of simplicity, as much in the theoretical analysis as in the practical implementation and computational time, studies usually start with data given in uniform grids. Nevertheless, some applications require dealing with data over nonuniform grids. At times, it is not trivial to adapt operators defined over uniform grids to the nonuniform case. The above-mentioned PPH operator was defined over a uniform grid and some of its properties were studied in [1]. These reconstruction operators are the basis for the definition of associated subdivision and multi-resolution schemes. In this paper, we use the definition that we made of the PPH reconstruction operator for data over nonuniform grids in [2], and we study some properties of this operator in greater depth. In particular, we focus on the smoothness of the reconstruction and the convexity-preserving properties of the initial data. We show that PPH reconstruction gives a C ∞ function, except for the knots where the function remains C 0 and the differences between the first, second, and third one-sided derivatives are of the third, second, and first order, respectively (see Definition 5).
In [3], the authors proved that the related subdivision scheme in uniform meshes preserves the convexity of the control points. In this article, we attempt to determine if this result about preserving convexity can be extended for the reconstruction operator and not only in uniform meshes, but also in σ quasi-uniform meshes with σ ≤ 4.
The paper is organized as follows: Section 2 is devoted to defining the PPH reconstruction operator over nonuniform grids. For this purpose, we will use the weighted harmonic mean with appropriate weights. Then, we show that the new reconstruction operator amounts to the original PPH reconstruction operator when we restrict to uniform grids. The definition is given for general nonuniform meshes, although in order to establish some theoretical results, we consider σ quasi-uniform meshes. In Section 3, we study some basic properties of PPH reconstruction, such as the reproduction of polynomials of the second degree, approximation order, smoothness, boundedness of the operator, Lipschitz continuity, and convexity preservation. In Section 4, we analyze the reconstruction when dealing with strictly convex (or concave) initial data. In Section 5, we present some numerical tests. Finally, some conclusions are included in Section 6.

A Nonlinear PPH Interpolation Procedure on Nonuniform Grids
Let us define the nonuniform grid X = (x i ) i ∈ Z . Let us also denote the nonuniform spacing between abscissae as h i := x i − x i−1 . We will work with continuous piecewise reconstructions of a given underlying continuous function f (x) with, at most, a finite set of isolated corner or jump discontinuities, that is, where R j (x) is a third-degree polynomial satisfying From now on, we will use the notation f i := f (x i ). Taking (1) into account, this implies that we are interested in building the appropriate polynomial piece R j (x) in the interval [x j , x j+1 ]. Let us consider the set of values { f j−1 , f j , f j+1 , f j+2 } for some j ∈ Z corresponding to subsequent ordinates of a function f (x) at the abscissae {x j−1 , x j , x j+1 , x j+2 } of a nonuniform grid X, and p L j (x) is the Lagrange interpolatory polynomial built with the points (x i , f i ), i = j − 1, j, j + 1, j + 2, that is, the unique polynomial of degree less or equal 3 satisfying The polynomial p L j (x) can be expressed as where It is well known that from conditions (3), one obtains the following linear system of equations, where the coefficient matrix is a Vandermonde matrix with different nodes and is, therefore, invertible: In order to express the solution of system (5) in a form that is convenient for our purposes, we introduce the definition of the second-order divided differences, and a weighted arithmetic mean of D j and D j+1 , defined as with the weights With these definitions, after solving the system (5), we get the following expressions for the coefficients of the polynomial (4): which can also be expressed as At this point, we give some more definitions and lemmas that we will need later.

Lemma 1.
Let us consider the set of ordinates { f j−1 , f j , f j+1 , f j+2 } for some j ∈ Z at the abscissae {x j−1 , x j , x j+1 , x j+2 } of a nonuniform grid X = (x i ) i ∈ Z . 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 Proof. This proof is an immediate calculation just by expanding the expression of the weighted arithmetic mean in (7) in terms of f i , i = j − 1, j, j + 1, j + 2, that is, In the presence of isolated singularities, predictions made using Lagrange reconstruction operators lose their accuracy in the vicinity of the discontinuity; in fact, three intervals are expected to be affected, since we are considering a stencil of four points. In order to reduce the affected intervals to only one, the one containing the singularity, we introduce a weighted harmonic mean over nonuniform grids, which will be used in the general definition of the PPH reconstruction operator. Notice that it is not possible to recover the exact position of a jump discontinuity inside an interval by using point value discretization of an underlying function. For the case of a corner discontinuity, a strategy such as the subcell resolution technique [4] could be used to detect its position. This harmonic mean is built as the inverse of the weighted arithmetic mean of the inverses of the given values. We define the following function. Definition 2. 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: Before giving another important lemma for our purposes, we will introduce a definition about a basic concept that will be used throughout the rest of the article.
Lemma 3. Let a > 0 be a fixed positive real number and let x ≥ a, y ≥ a. If |x − y| = O(h), then the weighted harmonic mean is also close to the weighted arithmetic mean M(x, y) = w x x + w y y, Remark 1. The smaller the value of a > 0 in Lemma 3, the smaller the h 0 in Definition 3 required to attain the expected theoretical order.
It is well known that the divided differences (6) work as smoothness indicators [5][6][7][8][9][10][11][12][13][14]. If a potential singularity appears at the interval [x j+1 , x j+2 ], we propose that the data (x j+2 , f j+2 ) are not interpolated, and that the ordinate f j+2 is exchanged for another value that is more convenient for what happens in the target interval [x j , x j+1 ], where we want to implement the local polynomial piece according to (1). In the same manner, if a potential singularity lies in the interval [x j−1 , x j ], a symmetrical modification is carried out. According to these observations, we can give the following definition for the PPH reconstruction on nonuniform meshes.

Definition 4 (PPH reconstruction)
. Let X = (x i ) i∈Z be a nonuniform mesh. Let f = ( f i ) i∈Z be 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: where γ j,i , i = −1, 0, 1, 2 are given in (12) and (14) with the weights w j,0 and w j,1 in (8). We define the PPH nonlinear reconstruction operator as where PPH j (x) is the unique interpolation polynomial that satisfies According to Definition 4, and establishing a parallelism with expression (4), we can write the PPH reconstruction as where the the coefficients a j,i , i = 0, . . . , 3 are calculated by imposing conditions (20). Depending on the local case, Case 1 or Case 2, the coefficients will have different expressions.
The replacement of f j+2 with f j+2 by exchanging the weighted arithmetic mean in Equation (11b) for its corresponding weighted harmonic mean has been proposed. It is also important to point out that Equation (17) shows that f j+2 is not significantly affected by a potential singularity at the interval [x j+1 , x j+2 ] , since, by property (15) in Lemma 2, | V j | ≤ 1 w j,0 |D j |, and in turn, D j is not affected by this discontinuity. Therefore, the influence of f j+2 on the values of the reconstruction in the interval [x j , x j+1 ] will be limited. In this case, the coefficients of the new polynomial (21) come from solving the linear system  Thus, the coefficients a j,i , i = 0, . . . , 3 take the form Case 2: |D j | > |D j+1 |, i.e., the possible singularity lies in [x j−1 , x j ]. In this case, in Definition 4, the value f j−1 is replaced with f j−1 by using expression (18). The net effect is again the exchange of the weighted arithmetic mean in Equation (11a) for the corresponding weighted harmonic mean. On this occasion, we get an adaption of the reconstruction to a potential singularity in [x j−1 , x j ], since the effect of the value f j−1 is largely reduced. In fact, by property (15) in Lemma 2, | V j | ≤ 1 w j,1 |D j+1 |, and D j+1 is not affected by any discontinuity. By solving the system (22), we obtain the following coefficients for the polynomial ( 21): Remark 2. The replacement of the weighted arithmetic mean for the corresponding harmonic mean in Definition 4 does not only guarantee adaptation near singularities, but also enlarges the region where the reconstruction preserves convexity according to expressions (40) and (43), as we will see in the next section.

Remark 3.
In both cases, the value of the PPH reconstruction at the midpoint x j+ 1 2 of x j , x j+1 gets the value PPH j (x j+ 1 2 ) = a j,0 . This expression directly defines an associated subdivision scheme and, consequently, also an associated multi-resolution scheme in nonuniform meshes. The interested reader is referred to [1,5] for more details in the context of uniform meshes.

Remark 4.
Notice that, considering uniform meshes, i.e., h i = h ∀i, all the given expressions reduce to the equivalent expressions in [1], which are valid only for the uniform case.

Remark 5.
Notice that Definition 4 of the PPH reconstruction operator has been given for general nonuniform meshes. From now on, one needs to take into account that some results are true for general grids, while others need the restriction to σ quasi-uniform meshes.

Main Properties of the PPH Reconstruction Operator in Nonuniform Meshes
In this section, we study some interesting properties of the new reconstruction operator. More precisely, we study the reproduction of polynomials, accuracy of the reconstruction, smoothness, boundedness, Lipschitz continuity, and convexity preservation. We start with the reproduction of polynomials up to degree 2.

Reproduction of Polynomials up to Degree 2
If the underlying function f (x) is a polynomial of degree 2, then D j = D j+1 = D is constant and D j D j+1 = D 2 ≥ 0. Using Equations (7), (10), (14), (23), and (24), we get So, PPH j (x) = p L j (x), i.e., PPH j (x), reproduces polynomials of a degree less than or equal 2, since p L j (x) does this.

Approximation Order for Strictly Convex (Concave) Functions
We will prove full-order accuracy, that is, fourth order, for a reconstruction that locally uses four centered points to get the approximation at a given interval [x j , x j+1 ] for any j ∈ Z. In particular, we can enunciate the following proposition.
Proposition 1. Let f (x) be a strictly convex (concave) function of class C 4 (R) and let a ∈ R, a > 0 be such that f (x) ≥ a > 0, ∀x ∈ R (let a ∈ R, a < 0 be such that f (x) ≤ a < 0, ∀x ∈ R). Let X = (x i ) i∈Z be a σ quasi-uniform mesh in R, with h i = x i − x i−1 , ∀i ∈ Z, and f = ( f i ) i∈Z , the sequence of point values of the function f (x), ]. This implies that PPH(x) = PPH j (x). Now, let us suppose that the initial data f = ( f i ) i∈Z come from a strictly convex function (for a concave function, the arguments remain the same) satisfying the given hypothesis f (x) ≥ a > 0, ∀x ∈ R for some a > 0. Then, D j D j+1 > 0, since second-order divided differences amount to second derivatives at an intermediate point divided by two, i.e., with µ 1 ∈ (x j−1 , x j+1 ) and µ 2 ∈ (x j , x j+2 ). Moreover, we have where M is a bound of the third derivative of f (x) in the compact interval [x j−1 , x j+2 ].
Since from Equations (7) and (14), we can write Putting this information into (9) and (23), if |D j | ≤ |D j+1 |, or into (10) and (24), if |D j | > |D j+1 |, we get that Thus, where p L j (x) is the Lagrange interpolatory polynomial. Taking the triangular inequality into account again, and using that Lagrange interpolation also attains fourth-order accuracy.

Smoothness
In this part, we study the smoothness of the resulting reconstruction, and for this purpose, we give the following definition. C s function). A function f : R → R is said to be quasi C s (R) if it satisfies: (a) f (x) belongs to class C s (R) except for a numerable set of points X = (

Definition 5 (Quasi
(b) There exist one-sided derivatives until order s, f (m) (x + i ) and f (m) (x − i ), m = 0, . . . , s, and Before giving the main result regarding smoothness, we will prove an auxiliary lemma that we need. in (a, b), and let us suppose that there exist h 0 , M > 0 and r ≥ 1 such that ∀ 0 < h < h 0 : Proof. From the fact that f is derivable in (a, b), we have that given x ∈ (a, b) for all ε > 0, there exists h ε > 0 such that, ∀h : 0 <h < h ε , Let h ∈ (0, h 0 ); then, we take ε h := h r−1 , and there exists h ε h > 0 such that, ∀h : We now define h 1 = min{h, h ε h }. Then, for allh withh ∈ (0, h 1 ), we get: We are now ready to present the following proposition with respect to the PPH reconstruction given in Definition 4.

Proposition 2. Let f (x) be a strictly convex (concave) function of class C 4 (R) and a
Proof. By construction, the PPH reconstruction is C ∞ ((x i , x i+1 )) for all i ∈ Z, since it is nothing else but a piecewise polynomial. Let us study the situation at a grid point x j where two polynomial pieces join. Again, by construction, PPH j−1 (x j ) = PPH j (x j ), and therefore, the reconstruction is a continuous function. Using the proof of Proposition 1, we know that From (28), we get that Thus, from Lemma 4, we get that In particular, Equations (29) and (30) are true for the abscissa x j , which proves the property of quasi C 3 at the grid points.

Boundedness and Lipschitz Continuity
We start by giving the exact definitions of the concepts treated in this section.

Definition 6.
A nonlinear reconstruction operator R : l ∞ (Z) → C(R) is called bounded if there exists a constant C > 0 such that

Definition 7.
A nonlinear reconstruction operator R : l ∞ (Z) → C(R) is called Lipschitz continuous if there exists a constant C > 0 such that ∀ x, y ∈ R, it is verified that Before addressing these properties, we need to prove some lemmas. Lemma 5. Let X = (x i ) i∈Z be a σ quasi-uniform mesh in R, with h i = x i − x i−1 ∀i ∈ Z, and let L m (x) m = −1, 0, 1, 2 be the Lagrange basis for a four-point stencil {x j−1 , x j , x j+1 , x j+2 }. Then, Proof. As is well known, the Lagrange bases are given by In order to obtain the bound for L 1 (x), we distinguish two cases.
, working in a similar way, we also get Finally, Lemma 6. Let X = (x i ) i∈Z be a σ quasi-uniform mesh in R, with h i = x i − x i−1 ∀i ∈ Z, and let us consider the expressions f j+2 in (17) and f j−1 in (18). Then, we have the following bounds: (12), we get

Proof. From Equations
According to property (15) of the harmonic mean | V j | ≤ |D j | w j,0 , , we also get .
Following a similar path for f j−1 , Using the property (15) of the harmonic mean | V j | ≤ |D j+1 | w j,1 , we get , which leads us to . Proposition 3. The nonlinear PPH reconstruction operator is a bounded operator over σ quasiuniform meshes.
Proof. Let x ∈ R and j ∈ Z such that x ∈ [x j , x j+1 ]. Depending on the relative size of D j and D j+1 , the PPH reconstruction operator replaces the value f j+2 with f j+2 or f j−1 by f j−1 as follows: where B m = L m (x), m = −1, 0, 1, 2, stand for the Lagrange polynomials. Applying the triangular inequality for each case, we get According to Lemmas 5 and 6, we obtain the following bound for both cases: The following lemma will be used for proving the Lipschitz continuity.
Lemma 7. Let X = (x i ) i∈Z be a σ quasi-uniform mesh in R, with h i = x i − x i−1 ∀i ∈ Z, and let f = ( f i ) i∈Z be a sequence in l ∞ (Z). Then, the nonlinear reconstruction operator defined in (4) satisfies that ∀j ∈ Z, The reconstruction PPH j (x) and its derivative read Without lost of generalization, we will suppose that |D j | ≤ |D j+1 |. The case |D j | > |D j+1 | can be carried out similarly. First, we prove the following inequalities: where h min = min i∈Z h i depends on the particular σ quasi-uniform mesh.
Using the expressions (23) for the coefficients of the polynomial derivative in (35), we have Thus, where C = 1 h min 2 + 2σ + 4 3 σ 2 depends on the σ quasi-uniform mesh.

Proposition 4.
The nonlinear PPH reconstruction operator is Lipschitz continuous over σ quasiuniform meshes.
Proof. Let us suppose first that there exists j ∈ Z such that x, y ∈ [x j , x j+1 ]. Using the Lagrange mean value theorem, ∃ θ ∈ (x j , x j+1 ), such as Thus, using Lemma 7 now, we get In the general case, we can suppose that x < y, x ∈ [x j 1 , x j 1 +1 ], y ∈ [x j 2 , x j 2 +1 ] with j 1 ≤ j 2 . If j 1 = j 2 ,, we have already proved the result. For j 1 < j 2 ,

Convexity Preservation
We first introduce a definition concerning what we call strictly convex data and a strictly convexity-preserving reconstruction operator. Definition 8. Let X = (x i ) i∈Z be a nonuniform mesh in R, with h i = x i − x i−1 ∀i ∈ Z, and let f = ( f i ) i∈Z be a sequence in l ∞ (Z). We say that the data are strictly convex (concave) if, for all i ∈ Z, it is satisfied that D i > 0 (D i < 0), where D i stands for the second-order divided differences. Definition 9. Let X = (x i ) i∈Z be a nonuniform mesh in R with h i = x i − x i−1 ∀i ∈ Z, and let f = ( f i ) i∈Z be a strictly convex (concave) sequence. We say that an operator R : Next, we give a proposition that introduces sufficient conditions on the grid for convexity preservation of the proposed reconstruction.
Proposition 5. Let X = (x i ) i∈Z be a σ quasi-uniform mesh in R, h i = x i − x i−1 , ∀i ∈ Z, and σ ≤ 4. Let f = ( f i ) i∈Z be a sequence of strictly convex data. Then, the reconstruction PPH(x) is strictly convexity preserving in each (x j , x j+1 ), that is, it is a piecewise convex function satisfying Proof. Let x inR and j ∈ Z such that x ∈ (x j , x j+1 ). Let us also consider that D i > 0∀i ∈ Z. The case D i < 0∀i ∈ Z is proved in the same way.
Computing derivatives in Equation (21), we get PPH j (x) = 2 a j,2 + 6 a j,3 x − x j+ 1 2 . (38) In order to analyze the sign of PPH j (x) we need to consider two cases due to the fact that the expression of PPH j (x) is different for |D j | ≤ |D j+1 | than for |D j | > |D j+1 |.
Replacing coefficients a j,2 , a j,3 coming from Equation (23) in expression (38) results in Taking into account that V j − D j ≥ 0, from (39), we get that proving PPH (x) > 0 is trivial if V j = D j . Otherwise, the inequality PPH (x) > 0 reads Replacing V j with its expression in Equation (14), we obtain Evaluating the previous expression at x j , we obtain the condition for convexity preservation in (x j , x j+1 ). This condition reads Since X is a σ quasi-uniform mesh with σ ≤ 4, we have h j+1 ≤ 2(h j + h j+2 ) , and therefore, the condition (42) is immediately satisfied. This proves the proposition in this case.
This time, by replacing the coefficients a j,2 , a j,3 coming from Equation (24) in expression (38) and following a similar track to that in Case 1, we obtain expressions similar to ( 40) and (41) for the abscissae verifying PPH j (x) > 0: x < x j+ 1 Now, evaluating at x j+1 , we get Thus, since X is a σ quasi-uniform mesh with σ ≤ 4, we get the result.

Remark 6.
As can be observed in expressions (41) and (44), the conditions that assure the strictly convexity-preserving property depend on the second-order divided differences of the initial data. The hypotheses of Proposition 5 are only sufficient conditions, but not necessary conditions.

Remark 7.
Working in a similar way with the Lagrange reconstruction operator p L j (x), we obtain the following expression that is analogue to (41) for the abscissa-fulfilling condition p L j (x) > 0: Then, if we are under the supposition that D j < D j+1 , calling x PPH and x p L to the second members of inequalities (41) and (46), respectively, we get i.e., PPH reconstruction operator preserves the strict convexity in a wider interval than the Lagrange reconstruction operator does. A similar conclusion can be reached under the supposition that D j > D j+1 .

PPH Reconstruction Operator over σ Quasi-Uniform Meshes for Strictly Convex (Concave) Initial Data
In this section, we gather the most important properties of the presented PPH reconstruction for strictly convex (concave) starting input data, and we give them in a unifying theorem. We want to emphasize the potential practical importance of the studied technique for designing processes.
Theorem 1. Let f (x) be a strictly convex (concave) function of class C 4 (R) and let a ∈ R, a > 0 such that f (x) ≥ a > 0, ∀x ∈ R (a ∈ R, a < 0 such that f (x) ≤ a < 0, ∀x ∈ R). Let X = (x i ) i∈Z be a σ quasi-uniform mesh in R with h i = x i − x i−1 , ∀i ∈ Z, and let f = ( f i ) i∈Z be the sequence of point values of the function f (x), f i = f (x i ). Then, the reconstruction PPH(x) satisfies 1.
Reproduction of polynomials up to the second degree.

3.
It presents a quasi C 3 smoothness.

4.
It is bounded and Lipschitz continuous.

5.
It is strictly convexity preserving in each (x j , x j+1 ).
Proof. Taking the previous results into account, the proof of this theorem is now immediate. In fact, the first affirmation is proven in Section 3.1, and the rest of the affirmations are proven in Propositions 1-5, respectively.

Numerical Experiments
In this section, we present three simple numerical experiments. The first one is dedicated to comparing the convexity preservation between the Lagrange and PPH reconstructions. Let us consider the initial convex set of points, (0, 10), (8,9), (25, 12), and (30, 30) , that is, D j > 0, D j+1 > 0 . In Figure 1, we have depicted the reconstruction operators corresponding to Lagrange and PPH, and we have marked with triangles the inflection points for each reconstruction (5.66 and 10.16, respectively). We observe that PPH preserves convexity in a wider range h j+1 + 2h j+2 6 = 4.5 than the Lagrange reconstruction does (see expression (47)). In fact, according to Theorem 1, PPH reconstruction is strictly convexity preserving for the abscissae corresponding with the central interval (8,25), while Lagrange reconstruction is not. to f (x) − PPH(x) ∞ , thus evaluating a much denser set of points. Both the errors and the approximation orders p for each iteration are shown in Table 1, where we can see that PPH reconstruction tends to fourth-order accuracy with this smooth concave function, as is expected according to Proposition 1.
The appropriate behavior of the reconstruction can be checked in Figure 2, where the preservation of the concavity and the accuracy of the approximation can be observed. Table 1. Approximation errors E k in the l ∞ norm and corresponding approximation orders p obtained after k iterations for the PPH reconstruction with f (x) = sin(x), k = 0, 1, ..., 7.

Conclusions
We have defined and studied the PPH reconstruction operator over nonuniform grids, paying special attention to the case of σ quasi-uniform grids and initial data coming from strictly convex (concave) underlying functions.
We have theoretically proven some very interesting properties of the new reconstruction operator from the point of view of a potential use in graphical design applications. These properties include the reproduction of polynomials up to the second degree, approximation order, smoothness, boundedness of the operator, Lipschitz continuity, and convexity preservation. In particular, we would like to emphasize the quasi C 3 smoothness of the operator and the preservation of strict convexity according to the result contained in Theorem 1.
In the section on the numerical experiments, we checked that the behavior corresponded to the developed theory-in particular, the reconstruction attained fourth-order accuracy and preserved the convexity of the initial data. The results clearly show that the reconstruction introduces improvements in comparison with the Lagrange reconstruction. Therefore, the numerical experiments that we carried out reinforce the theoretical results.