Analysis of a New Nonlinear Interpolatory Subdivision Scheme on σ Quasi-Uniform Grids

: In this paper, we introduce and analyze the behavior of a nonlinear subdivision operator called PPH, which comes from its associated PPH nonlinear reconstruction operator on nonuniform grids. The acronym PPH stands for Piecewise Polynomial Harmonic, since the reconstruction is built by using piecewise polynomials deﬁned by means of an adaption based on the use of the weighted Harmonic mean. The novelty of this work lies in the generalization of the already existing PPH subdivision scheme to the nonuniform case. We deﬁne the corresponding subdivision scheme and study some important issues related to subdivision schemes such as convergence, smoothness of the limit function, and preservation of convexity. In order to obtain general results, we consider σ quasi-uniform grids. We also perform some numerical experiments to reinforce the theoretical results.


Introduction
Subdivision schemes are closely related to reconstruction operators. They have been used in the last few decades in many applications ranging from the numerical solution of partial differential equations to image processing and computer-aided geometric design. Subdivision schemes give simple and fast algorithms to approximate the limit function from a set of initial data at a coarse resolution level. There is an method of immediately generating subdivision schemes from reconstruction operators and more specifically from prediction operators [1][2][3]. Due to this connection, subdivision schemes inherit many of the properties of their associated reconstruction operators. In particular, the subdivision scheme is nonlinear if the reconstruction operator is nonlinear and it is said interpolatory if it comes from a reconstruction operator that is an interpolation.
Nonlinear subdivision schemes have emerged as good candidates to adapt to the concrete data in use. The research in this field has expanded, with new contributions each year, and has received the attention of many researchers; see for example [4][5][6][7][8]. Nonlinearity means data-dependent subdivision schemes, which may also involve nonlinear operations in their definition. Then, by definition, they are designed to overcome certain drawbacks that appear when dealing with their linear counterparts, such as bad behavior in the presence of isolated discontinuities for instance. An example of these types of operators was defined in [1] and was named PPH (Piecewise Polynomial Harmonic). This scheme basically consists on a clever modification of the classical four-point Lagrange subdivision scheme. Several studies have been carried out about their properties and performance in different applications, see for example [1,9,10]. Two main purposes of this subdivision scheme are related to dealing with data containing isolated discontinuities, reducing the undesirable effects, and preserving the convexity of the initial data while maintaining a centered support based on four points.
In [11], the authors extended the definition of the PPH reconstruction operator to nonuniform grids. In turn, this fact allows us to extend the PPH subdivision scheme to nonuniform grids and to carry out a parallel study in this new scenario. In order to overcome some technical difficulties in the theoretical proofs, we restricted some results to σ quasi-uniform grids. The resultant scheme is quite interesting in terms of applications due to the almost C 1 smoothness of the limit function, allowing us to approximate accurately continuous functions with corners, and due to appropriate properties regarding convexity preservation of the initial data; see [9]. In this paper, we focus on proving the convergence of the scheme towards an almost C 1 limit function and we address numerically the issue of stability, which is a central issue in order to be useful for applications.
The paper is organized as follows: Section 2 is devoted to a review of the PPH reconstruction operator over nonuniform grids. Section 3 presents a short review about Harten's interpolatory multiresolution setting, which is closely connected to interpolatory subdivision schemes. In Section 4, we define the associated subdivision scheme, which we show amounts to the PPH subdivision scheme when use a restriction 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 5, we analyze the main issues about subdivision schemes. In particular, we prove some results about convergence, smoothness of the limit function, and convexity preservation. In Section 6, we carry out some numerical tests to check the theoretical smoothness of the limit function and the performance of the nonlinear subdivision scheme. Finally, we give some conclusions in Section 7.

A Nonlinear PPH Reconstruction Operator on Nonuniform Grids
In this section, we review the definition of the nonlinear reconstruction that gives rise to the nonlinear subdivision scheme under study in this article. More information about this reconstruction operator can be found in [1,11,12].
Let us define a nonuniform grid X = (x i ) i ∈ Z . Let us also denote h i := x i − x i−1 , the nonuniform spacing between abscissae. Let us consider the set of values { f j−1 , f j , f j+1 , f j+2 } for some j ∈ Z corresponding to the abscissae {x j−1 , x j , x j+1 , x j+2 } of the nonuniform grid X.
We need to introduce the definition of the second-order divided differences and the weighted arithmetic mean of D j and D j+1 defined as with the weights We require also some definitions and lemmas that appear in [11]. 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 V as the function Lemma 1. If x ≥ 0 and y ≥ 0, the harmonic mean is bounded as follows: Lemma 2. Let a > 0 be 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, We review the following definition for the PPH reconstruction on nonuniform meshes. The details and main properties of this reconstruction operator can be found in [11,12].
Definition 2 (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 by and V j = V(D j , D j+1 ), with V being the weighted harmonic mean defined in (4) with the weights w j,0 and w j,1 in (3). We define R(x) as the PPH nonlinear reconstruction operator given by where R j (x) is the unique interpolation polynomial that satisfies We can write the PPH reconstruction by using the middle point where the the coefficients a j,i , i = 0, . . . , 3 are calculated by imposing conditions (11). Depending on the local case, Case 1 or Case 2, the coefficients have different expressions. Case 1. |D j | ≤ |D j+1 |, In this case, the coefficients of the polynomial (12) take the form Case 2. |D j | > |D j+1 |, In this case, we obtain the following coefficients for the polynomial (12) With the previous definitions and lemmas, we are now ready to introduce the PPH subdivision scheme. However, before doing this, we also review some basic concepts of Harten's interpolatory multiresolution setting and its connection with subdivision schemes.

Harten's Interpolatory Multiresolution Setting
Let us consider a set of nested grids in R, and the point-value discretization where V k is the space of real sequences related to the resolution of X k and C B (R) is the set of bounded continuous functions on R. A reconstruction operator R k associated with this discretization is any right inverse of D k , which means that for all f k ∈ V k , R k f k ∈ C B (R), and D k R k = I, that is The sequences {D k } k∈N and {R k } k∈N define a multiresolution transform [2]. The prediction operator, i.e., D k+1 R k : V k → V k+1 , defines a subdivision scheme. Relation (16) implies that the subdivision scheme is interpolatory. If R k is a nonlinear reconstruction operator, then the corresponding subdivision scheme S := D k+1 R k becomes also nonlinear.

A Nonlinear PPH Subdivision Scheme on Nonuniform Grids
Let us consider a particular set of nonuniform nested grids X k = (x k i ) i∈Z , k ≥ 0, generated from an initial grid X 0 . Definition 3. Given X 0 = {x i } i∈Z as a nonuniform grid in R, we define, for k ∈ N (the larger the k, the larger the resolution), as the set of nested grids given by Let us also consider h k i = x k i − x k i−1 , the nonuniform spacing between abscissae. Given a set of control points f k = ( f k i ) i∈Z , we define the nonlinear PPH subdivision scheme as where (4) and it is computed with the weights (w k i,0 ) i∈Z , (w k i,1 ) i∈Z is given in (3), and the second-order divided differences D k i and D k i+1 are defined in (1).
Notice that the expression of the subdivision scheme at odd indexes coincides with the coefficientã 0 of the PPH reconstruction operator in (13) or (14) due to the fact that the defined subdivision scheme satisfies S = D k+1 R k . This means that the expression of the subdivision scheme is symmetric, even if the modification of the data has been carried out to the left (14) or to the right (13) for the concrete piece of the underlying reconstruction operator.
Suppose that the initial data come from a convex smooth function; then by definition through its associated reconstruction operator, we get a fourth-order accurate subdivision scheme. When the data come from an underlying smooth function with inflexion points, the order is reduced around these inflexion points [12]. The use of the weighted harmonic mean in Definition 17 guarantees certain adaptation near jump discontinuities. In the presence of an isolated singularity, we have two adjacent intervals where D i = O(1) and For these cases, the harmonic mean maintains an order of V i = O(1). If both D i and D i+1 are affected by discontinuity, then no adaption takes place. However, this situation occurs only in the prediction of one value per scale and per discontinuity. It is also interesting to remark that, for uniform meshes, i.e., h i = h ∀i, then all of the given expressions reduce to equivalent expressions in [1] valid only for the uniform case.
Notice that Definition 17 of the PPH subdivision schemes was introduced for general nonuniform meshes. From now on, one needs to take into account that some results are true for general grids while others require restriction to a particular type of nonlinear mesh that is the most common in practice.
In next the section, we study some main issues about the defined subdivision scheme. In particular we prove convergence, almost C 1 smoothness in the limit function, and we provide a result concerning convexity preservation.

Main Properties of the PPH Subdivision Scheme in Nonuniform Meshes
We start the section with some definitions taken from [1] that are used in the rest of the article.

Definition 4.
A nonlinear subdivision scheme is called uniformly convergent if for every set of initial data f 0 ∈ l ∞ (Z), there exists a continuous function S ∞ f 0 ∈ C(R), such that

Definition 5.
A convergent nonlinear subdivision scheme is called stable if there exists a constant C such that for every pair of initial data f 0 ,f 0 ∈ l ∞ (Z),

Definition 6.
Let N ≥ 0 be a fixed integer. A nonlinear interpolatory subdivision scheme has the property of polynomial reproduction of order N if for all P ∈ Π N , where Π N stands for the vector space of polynomials of degree less or equal to N, we have S p =p, where p andp are defined by p k = P(2 −k ·) andp k = P(2 −(k+1) ·).

Definition 8.
A nonlinear subdivision scheme is called Lipschitz continuous if there exists a constant C > 0 such that, for every f , g ∈ l ∞ (Z), the following is verified: We can now provide some basic results before addressing the convergence of the scheme. In order to prove the coming theoretical results, we work with σ quasi-uniform grids, according to the following definition: Proposition 1. The nonlinear subdivision scheme associated wtih the PPH reconstruction: (1) reproduces polynomials of degree N ≤ 2,
(1) If f is a polynomial of degree less or equal to 2, therefore, the proposed scheme reproduces polynomials of degree 2.
(2) By definition of the PPH subdivision scheme for a given j ∈ Z, we have that Thus, and therefore the nonlinear subdivision scheme is bounded. .
and in the second case, we similarly obtain then using the same arguments as in case (a), we obtain If D j ( f )D j+1 ( f ) < 0, we consider the function Z(x, y) = xy w j,0 y+w j,1 x defined for all xy > 0. It is easy to check that the Jacobian of the function Z verifies Thus, the mean value theorem easily leads to Clearly, C = 1 + max{2σ 3 , σ 4 } is a convenient constant that completes the proof.
The next lemma and proposition allow us to prove the existence of a contractive scheme S 1 for the differences Lemma 3. Let D be the set defined by D := j ∈ Z : D j D j+1 > 0 , and let the expressions E 1j , E 2j and M be defined as follows: Then, the following inequalities are satisfied Proof. (1) and (2) are trivial. Let us break down (3). Given j ∈ D, we have Proposition 2. Associated with the PPH nonlinear reconstruction, on nonuniform grids, there exists a nonlinear subdivision scheme S 1 for the differences. If the grid is σ-quasy uniform, then S 1 is bounded, i.e., satisfies , then λ 1 < 1 and S 1 is contractive.
Proof. (a) Existence of S 1 . S 1 δ k has the following expressions for even and odd indexes. (a.1) Even indexes and depending on the value of V j we differentiate two cases, Then and again by proceeding in a similar way, we obtain the following for the two different cases, We consider again even and odd indexes. Thus, The subdivision scheme S 1 is contractive if

Corollary 1.
For σ-quasy uniform grids where σ < 2, the scheme S 1 is contractive, since We now provide a simple and technical lemma to support the proof of next lemma. .
Proof. By definition of γ j,2 in (9) together with property (5) of the harmonic mean , and the expression of D k The case of γ j,−1 can be derived analogously.
We need two more lemmas that are used in the proof of Theorem 1.

Lemma 5.
Let {R k } be the sequence of nonlinear PPH reconstruction operators associated with a sequence of nested σ quasi-uniform grids {X k } satisfying Definition 3 and S as the PPH interpolatory subdivision scheme. There exists C ∈ R such that, if f k+1 = S f k , then ∀k, Proof. Let f k ∈ l ∞ (Z), and x ∈ R. Let j be such that x ∈ [x k j , x k j+1 ], and assume that where R L k stands for the centered Lagrange reconstruction operators of the same order. (1) We prove first the bound for the second term on the right-hand side. Since According to Lemma 4 in [11] |A Here, we review that the Lagrange polynomial bases sum to one: From now on, we drop the explicit dependence on x for the sake of clarity and write simply A m and B m when referring to these quantities. Since f k+1 = S f k and S is interpolatory, we have Taking into account property (5) of the harmonic mean | and similarly for the term in f k+1 2j+1 , Taking again into account the property (5) of the harmonic mean, we have | V k j | ≤ |D k j | w k j,0 and we can write Using (27) and (28), we get The modulus of the first term at the right-hand side can be rewritten as Then, using (25) and (26) |R (2) Let us now estimate |R L Let us suppose, without loss of generality, that we are in the first case, i.e., |Dj k | ≤ |D k j+1 | .
Using Definition (2) and applying the triangular inequality, we get Taking now into account that ∑ 2 s=−1 γ j,s = 0, we can rewrite the first term and we can bound it as follows: The second term on the right-hand side of (31) can be bounded using Lemma 4. Considering (30), (32), and Lemma 4, we have For the other case, |D k j | > |D k j+1 | using the same ideas, we also get the same bound. (3) Let us study now |R k+1 ( f k+1 )(x) − R L k+1 ( f k+1 )(x)|. Inequality (33) allows us to write Since by Proposition 2 the operator S 1 is bounded by σ 2 , we get that with C 3 = σ 2 C 2 . Finally, combining the results in (29), (33), and (34), we obtain the following which completes the proof.
The following theorem uses standard arguments and previous lemmas to prove the convergence of the nonlinear PPH subdivision scheme.
Theorem 1 (Convergence). Let {R k } be the sequence of nonlinear PPH reconstruction operators associated with a sequence of nested σ quasi-uniform grids {X k } with σ < 1 + 1 2M satisfying Definition 3. Then, the associated PPH interpolatory subdivision scheme S is uniformly convergent.
Proof. The basis of the proof is to observe that {R k ( f k )} k∈N is a Cauchy sequence in C B (R), the space of continuous and bounded functions in R. Let and from Proposition 2 ∃ C 2 ∈ R such that As σ < 1 + 1 2M , S 1 is contractive, which means (C 2 < 1) and lim k→∞ (C 2 ) k = 0.
Thus, given i.e., given > 0, In order to continue studying the degree of smoothness of the limit function, we need one more lemma. Lemma 6. Let {R k } be the sequence of nonlinear PPH reconstruction operators associated with a sequence of nested σ quasi-uniform grids {X k } satisfying Definition 3. The interpolatory PPH reconstruction operators R k have the following properties: (2) For each level k ≥ 1, for all x, y such that |x − y| < λ k−1 1 h 0 min , with λ 1 = 1 2 + (σ − 1)M < 1, the contractivity constant of the scheme S 1 of the differences and h 0 min = min j∈Z h k j , there exist a constant C such that Proof.
(1) The proof of this point can be found in Proposition 3 in [11].
(2) We prove property (2). We write According to Expression (33) inside the proof of Lemma 5, we have that Then, we focus now on the second term. Let us suppose x ∈ [x k j , x k j+1 ], and let us see that y ∈ [x k s , x k s+1 ] with |s − j| ≤ 4. Let us take the integer number which implies |s − j| ≤ 4. We now write Regrouping terms we can plug f k s into the previous formula as follows: Now taking into account that |B i | ≤ σ, |D i | ≤ σ, i = −1, 0, 1, 2 according to Lemma 4 in [11], and that |s − j| ≤ 4, we get which finishes the proof.
With all of these requisites, the limit function turns out to be Hölder continuous with α = 1.

Theorem 2 (Smoothness).
Let {R k } be the sequence of nonlinear PPH reconstruction operators associated with a sequence of nested σ quasi-uniform grids {X k } with σ < 1 + 1 2M satisfying Definition 3. Then, the associated PPH interpolatory subdivision scheme S is Hölder continuous with α = 1.

Proof.
In order to prove a Lipschitz condition for the limit function, we have that By using Lemma 5 and Proposition 2, we get If |x − y| ≥ h 0 min , then using the boundedness of the limit function S ∞ ( f ) derived from Theorem 1, we get If |x − y| < h 0 min , then there exists k ∈ N such that λ k 1 h 0 min < |x − y| < λ k−1 1 h 0 min . Thus, from point 2 of Lemma 6, we obtain and therefore, Then, from Proposition 2, Finally, from (38) and (40), we deduce Lipschitz condition, which completes the proof.
We complete our theoretical study with the important issue of preservation of convexity of the initial data. In order to address this question, we introduce two definitions.
Definition 11. An interpolatory subdivision scheme is said to be convexity preserving if and only if the data set {(x k i , f k i )} is strictly convex for all level k of subdivision.
Using these definitions, we can give the following theorem.
Theorem 3 (Convexity). Let {R k } be the sequence of nonlinear PPH reconstruction operators associated with a sequence of nested σ quasi-uniform grids {X k } with σ < 1 + 1 2M satisfying Definition 3. Then, the associated PPH interpolatory subdivision scheme S is convexity preserving if and only if Proof. The proof is based on the fact that, if D k i > 0, ∀i ∈ Z, at a given scale k ∈ N, then we have that the interpolatory subdivision scheme is convexity preserving if D k+1 2i+1 > 0 and D k+1 2i+2 > 0, ∀i ∈ Z. We start by computing D k+1 2i+1 > 0, .
Taking into account the relations between the scales k and k + 1, we get Using that the odd points at the scale k + 1 are predicted by (17), we obtain due to the fact that D k i > 0 and D k i+1 > 0.
Plugging the corresponding values for f k+1 2i+1 and f k+1 2i+3 according to (17) into last expression now, we arrive at .
After simple algebraical manipulations, we reach Considering that we already proved D k+1 2i+1 > 0 in (42), in order for the interpolatory subdivision scheme to be convexity preserving it remains only to evaluate D k+1 2i+2 > 0, that is which concludes the proof.

Corollary 2.
Let {R k } be the sequence of nonlinear PPH reconstruction operators associated with a sequence of nested σ quasi-uniform grids {X k } with σ < 1 + 1 2M satisfying Definition 3. If | V k i | < 2 min{D k i , D k i+1 }, ∀i ∈ Z, ∀k ∈ N, then the associated PPH interpolatory subdivision scheme S is convexity preserving.
Proof. Let us consider D k i > 0, ∀i ∈ Z, at a given scale k ∈ N. We get the following chain of inequalities: which proves the property of convexity preservation by applying Theorem 3. Proof. Let us consider D k i > 0, ∀i ∈ Z, at a given scale k ∈ N. Taking into account that V k i is a mean, we have and therefore, we can apply Corollary 2.

Remark 1.
If the initial data f 0 i , i ∈ Z come from a smooth function, we have the hypothesis of Corollary 3 satisfied for h 0 = max i∈Z h 0 i sufficiently small since due to the fact that µ 1 − µ 0 = O(h 0 ), where µ 0 , µ 1 , c are intermediate points between x k i−1 and x k i+2 .

Remark 2.
In the case of dealing with uniform grids, we have that V k i coincides with the classical harmonic mean, and therefore, | V k i | < 2 min{D k i , D k i+1 }, and Corollary 2 applies. Thus, we have a convexity preserving interpolatory subdivision scheme for any initial data.

Remark 3.
If instead of the weighted harmonic mean V k i we use the classical harmonic mean in the definition of the subdivision scheme given in (17), we immediately obtain a convexity preserving subdivision scheme because the hypothesis of Corollary 2 is met. However, we reduce the approximation order to second order in this case, while the original scheme comes from a reconstruction that is fourth-order accurate for strictly convex functions (see [11]).

Numerical Experiments
In this section, we carry out some numerical experiments to analyze the outputs obtained and to compare them with the expected theoretical results. Our first experiment is focused on the presented result about the smoothness of the limit function. We estimate the exponent α of the Hölder continuity of the limit function. In order to do so, we considered the following functions f (x) and g(x) given by x(cos (2πx) + 1), 0.3 < x ≤ 0.7, We also consider the point-value discretization f 0 given by the function values at a nonuniform grid X 1 with 30 points in the interval [0, 1]. Then, we carry out an estimation of the quotient for different levels k of refinement, k = 10, k = 15, and k = 17, and for different values of α, α = 0.75, α = 0.99, α = 1, α = 1.1, and α = 1.25.
In Figure 1, we show the original function considered in solid blue and the subdivision curve after k = 5 subdivision levels in dash-dotted black. In Table 1, we can observe that the constant C converges with the resolution levels to a fix value for α = 1. For values of α smaller than 1, the estimated value of C decreases with the number of resolution levels k, which means that the Hölder exponent of the subdivision scheme is higher. In turn, for values of α larger than 1, the estimated value of C increases with the number of resolution levels k, which means that the Hölder exponent of the subdivision scheme must be lower. Notice that the constant C in the definition of Hölder continuity depends on f (x) but must become stable as we approach the limit function with larger and larger k. We also carried out the same experiment, varying the number and position of the grid points, and for both functions given in (43). We used a nonuniform grid X 2 with 20 non-equally spaced abscissae. As can be seen in Tables 2 and 3, the results are consistent with our previous observations, obtaining in all cases an estimation for the Hölder exponent α = 1. However, we can appreciate that the value of the constant C depends not only on the function from which the point values are taken but also on the starting grid, since the limit functions for different grids are quite similar in the sense that they approximate the underlying function with fourth order, but they are not the same. In our second experiment, we only performed a comparison between the presented PPH subdivision scheme and the classical linear scheme with four points based on Lagrange interpolation. We plotted in Figure 1 the subdivision curve obtained for both methods, PPH and Lagrange, after k = 5 levels of subdivision and starting from the nonuniform grid X 1 with 30 initial points used in the first numerical experiment and the associated point values of the function f (x). The original function is also provided to compare the approximation capabilities of both subdivision methods. We see the original function in a solid blue line, the Lagrange subdivision scheme in a dashed red line, and the PPH subdivision scheme in a dash-dotted black line. As it can be appreciated in Figure 1, the Gibbs effects and undesirable oscillations due to the presence of a jump discontinuity are highly reduced with the PPH scheme in contrast with what happens with the linear scheme. Notice the high oscillations that appear near the jump discontinuities when using the linear scheme, which is known to happen when one implements any linear scheme. In Figure 1, to the right, we show a zoom of the area around the first jump discontinuity to observe more clearly the behavior of the nonlinear scheme. In Figure 2, we also plot the results obtained with the nonuniform grid X 2 with 20 non-equally spaced points considered in the previous experiment. We considered both functions f (x) and g(x). Again, the same type of Gibbs effects appear around the jump discontinuity for the linear method. The corner is not so problematic.
In Table 4, we see the errors || f k − S k f 0 || p , p = 1, 2, ∞, committed by approximating the original data f k , i.e., the right point-values of the function f (x) at the corresponding abscissas with S k f 0 for k = 5 subdivision levels, where S k f 0 stands for the iterative application k times of the analyzed subdivision schemes, namely PPH and Lagrange, starting from the initial function point values f 0 at the given grid X 1 with 30 abscissae. In Table 5, we give the corresponding results for the grid X 2 with 20 abscissae. In Table 6, we consider this time to be the errors ||g k − S k g 0 || p , p = 1, 2, ∞, for the function g(x) using the grid X 2 with 20 abscissae.   Table 4. Subdivision errors || f k − S k f 0 || p , p = 1, 2, ∞, committed by approximating the original data f k with S k f 0 for k = 5 subdivision levels starting from the initial function point values f 0 at the given grid X 1 with 30 points.

Method
|| f k − S k f 0 || 1 || f k − S k f 0 || 2 || f k − S k f 0 || ∞ PPH 0.0098 0.0454 0.4239 Lagrange 0.0388 0.1355 0.9767 Table 5. Subdivision errors || f k − S k f 0 || p , p = 1, 2, ∞, committed by approximating the original data f k with S k f 0 for k = 5 subdivision levels starting from the initial function point values f 0 at the given grid X 2 with 20 abscissae.  Table 6. Subdivision errors ||g k − S k g 0 || p , p = 1, 2, ∞, committed by approximating the original data g k with S k g 0 for k = 5 subdivision levels starting from the initial function point values g 0 at the given grid X 2 with 20 abscissae.

Conclusions
We defined and analyzed the PPH subdivision scheme on nonuniform grids, which were derived from its associated reconstruction operator. We paid special attention to the case of σ quasi-uniform grids and initial data coming from strictly convex (concave) smooth functions.
We theoretically proved some crucial issues when dealing with subdivision schemes, such as the existence of a contractive scheme for the first differences, convergence, smoothness of the limit function, and preservation of the convexity properties of the initial data.
In the numerical experiments section, we carried out some experiments that reinforce the theoretical results, in particular, we observed the Hölder continuity of the limit function, giving a numerical estimation of the exponent α that coincides with the result in Theorem 2. We also carried out another experiment to analyze the performance of the subdivision scheme with initial data that contain a numerical jump discontinuity, observing that Gibbs effects and oscillations are negligible. Finally, a potential real application in 2D was given by zooming in on some coarse data from geological areas corresponding to unaccessible seabeds.  Acknowledgments: We hank the anonymous referees and our colleagues S.Amat and D. Yañez for their help improving the article.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: