A family of non-oscillatory 6-point interpolatory subdivision schemes

In this paper we propose and analyze a new family of nonlinear subdivision schemes which can be considered non-oscillatory versions of the 6-point Deslauries-Dubuc (DD) interpolatory scheme, just as the Power p schemes are considered nonlinear non-oscillatory versions of the 4-point DD interpolatory scheme. Their design principle may be related to that of the Power p schemes and it is based on a weighted analog of the Power p mean. We prove that the new schemes reproduce exactly polynomials of degree three and stay ’close’ to the 6-point DD scheme in smooth regions. In addition, we prove that the first and second difference schemes are well defined for each member of the family, which allows us to give a simple proof of the uniform convergence of these schemes and also to study their stability as in [19, 22]. However our theoretical study of stability is not conclusive and we perform a series of numerical experiments that seem to point out that only a few members of the new family of schemes are stable. On the other hand, extensive numerical testing reveals that, for smooth data, the approximation order and the regularity of the limit function may be similar to that of the 6-point DD scheme and larger than what is obtained with the Power p schemes.


Introduction
Subdivision schemes are recursive processes used for the fast generation of curves and surfaces in computer-aided geometric design, as well as an essential ingredient in many multiscale algorithms used in data compression. In some applications, the given data need to be retained at each step of the refinement process, which requires the use of interpolatory subdivision schemes. The so-called Deslauries-Dubuc (DD henceforth) subdivision schemes [6] are a well known family of interpolatory subdivision schemes which can be interpreted as a recursive application of a piecewise polynomial interpolatory tool [10,11]. A general setting by which a piecewise polynomial interpolation technique can be used to provide the set of local rules that defines a subdivision scheme has been described in [10]: Assuming that χ l ⊂ χ l+1 are two nested grids on R m , f l is a set of known data associated to the grid χ l and I[x, ·] is a piecewise polynomial reconstruction technique, new data associated to the grid χ l+1 can be generated as follows This process allows to define a recursive subdivision scheme where sequences of values on denser and denser meshes are obtained according to the set of local rules derived from I[x, ·]. Clearly (1) leads to an interpolatory subdivision scheme if I is an interpolatory reconstruction, i.e. I[x l i , f l ] = f l i , ∀x l i ∈ χ l . For m = 1 and a binary refinement strategy, i.e. x l+1 2i = x l i and χ l+1 \χ l ≡ {x l+1 2j +1 } j ∈Z , we have so that values on a given mesh are 'copied' at the same location on higher resolution levels, while the interpolatory technique I[·, ·] specifies the local rules used for the generation of new data values. It is well known (see e.g. [9,10]) that the DD subdivision schemes can be written in the form (2) with I(·, ·) a Lagrange interpolatory reconstruction that considers an interpolatory stencil centered around the evaluation point. In general, the use of piecewise polynomial Lagrange interpolation based on a stencil that uses l points to the left and r points to the right of the evaluation point leads to a linear (i.e. data-independent) subdivision scheme that can be written as with mask coefficients a l,r j that can be computed from the interpolatory rule I[·, ·] (see [9,10]). It is well known that these schemes lead to a Gibbs-like oscillatory behavior when applied to discrete data with large gradients (see Fig. 1) and several nonlinear piecewise polynomial interpolatory techniques have been considered in the literature, within this same framework, in an attempt to construct interpolatory subdivision schemes that avoid undesired oscillations. Examples of such schemes are the ENO-WENO subdivision schemes [10,11,14,16] or the PPH scheme [2,4].
The oscillatory behavior displayed in Fig. 1 is typical of data-independent subdivision schemes based on Lagrange interpolation, which does not preserve the shape properties of data with large gradients when the degree of the polynomial pieces is larger than 1. ENO-WENO subdivision schemes manage to avoid the Gibbslike oscillatory behavior by selecting an appropriate stencil for the interpolatory reconstruction [10,11,14,16]. Other nonlinear interpolatory subdivision schemes, like the Power p schemes [17] or the shape-preserving schemes in [8], owe their non-oscillatory character to the judicious use of certain nonlinear averages.
The aim of this paper is to design, and analyze, non-oscillatory 6-point schemes that can be considered nonlinear analogs of the 6-point DD linear scheme, S 3,3 , in the same sense as the Power p schemes are considered nonlinear, non-oscillatory versions of S 2,2 , the 4 point DD subdivision scheme. For this, we shall use a weighted nonlinear average defined in [23] which generalizes the Power p mean defined in [15] and used in the design of the Power p schemes [17]. The new schemes proposed in this paper can be written in the following general form where F : l ∞ (Z) → l ∞ (Z) is a nonlinear operator, δ : l ∞ (Z) → l ∞ (Z) is linear and continuous and S L is a linear and convergent subdivision scheme. The Power p schemes and other related subdivision schemes considered in [1,3,22] can also be written in the form (4), which allowed the authors to study their convergence and stability by using the following results [1].
Theorem 1 Let S N be a nonlinear subdivision scheme which can be written in the form (4).
The scheme S N is uniformly convergent provided that F and δ satisfy the following conditions: C2. ∃L > 0, 0 < T < 1 : ||δS L N (f )|| ∞ ≤ T ||δf || ∞ ∀f ∈ l ∞ (Z) The scheme S N is (Lipschitz) stable, i.e. ∃C > 0 such that provided that F and δ satisfy the following conditions: Remark 2 If a scheme of the form (4) is convergent, the smoothness of the limit functions S ∞ N f , f ∈ l ∞ (Z), may be related to the smoothness of S ∞ L f . In particular it can be shown (see [12,17] The new schemes proposed in this paper share other features with the Power p schemes. For both families of schemes the linear operator in Eq. 4 can be considered as δ = ∇ 2 , where and the subdivision schemes are defined by piecewise smooth functions that are globally Lipschitz. The paper is organized as follows: in Section 2 we provide an explanation of the non-oscillatory character of the Power p schemes which can be used as a design tool to obtain new families of non-oscillatory 6-point interpolatory subdivision schemes. These shall require a nonlinear analog of the Power p mean, the Weighted Power p , proposed in [23]. The new families of 6-point schemes are defined and analyzed in Section 3. In this section we examine the polynomial reproduction properties and the existence of difference schemes, as well as the convergence and approximation properties of the new families of schemes. Section 4 is devoted to the issue of the stability of the new schemes. In Section 4.1 we study the Weighted Power p mean, and its Generalized Gradient, an essential ingredient in the application of the theory developed in [19] for the study of the stability of a nonlinear scheme. The limitations of this theory in the study of the stability of the proposed schemes are analyzed in Section 4.2. In Section 5 we study the stability issue from a computational point of view, and also present several numerical examples that illustrate our theoretical results. We close in Section 6 with some conclusions.
In addition, the relation between the Generalized Gradients of the piecewise smooth functions that define a nonlinear scheme, and the contractivity of such scheme is carefully explained in an Appendix to this paper.

Nonlinear averages and Non-oscillatory schemes
The Power p interpolatory subdivision schemes [17,19] are binary interpolatory subdivision schemes for which the generation of new data values (at odd points) is given by the following rule where is the so-called Power p mean [15], a nonlinear function that satisfies (see [15,17] for details) It is straightforward to see that ψ 2,2 in Eq. 3 can be written as The obvious similarity between Eq. 9 and Eq. 6 makes the Power p schemes nonlinear versions of the 4-point DD scheme. In addition, if f = (f i ) i∈Z , f i = F (x i ) with F (x) a smooth convex function, and X = {x i } an h-uniform grid, it can be proven that This property is obtained from the following relation (see e.g. [17]), which holds for x · y > 0, p ≥ 1 On the other hand, as shown in Fig. 2, the behavior of S H p when refining discrete data with large gradients is quite different from that of S 2,2 . In what follows, we give an explanation of the lack of oscillations observed in Fig. 2 which may be used to design 6-point nonlinear analogs of S 3,3 , the 6 point DD scheme. The starting point  Fig. 2 Non-oscillatory behavior of Power p schemes, S Hp , for p = 2, 3, compared to the 4-point DD scheme, S 2,2 in our construction is the following relation, which follows from Neville's Algorithm for Lagrange interpolation (see e.g. [13]). Moreover, it is not difficult to see that we can write, for l + r > 1 where L l,r is a linear operator such that (L l,r f ) 2n = 0, due to the interpolatory property. For the two schemes below hence From Eq. 12, for l = 2, r = 2, we get while for S H p in Eq. 6 we can write Taking into account (15) and Eq. 16, the behavior of the S 2,2 and S H p schemes may be explained in terms of the interpolatory stencils associated to the schemes S 2,1 and S 1,2 , shown in Fig. 4.
It is a well known fact that any Lagrange-type interpolatory technique suffers an O(1) accuracy loss as soon as the interpolatory stencil crosses a discontinuity. The data in Figs. 1-2 correspond to f i = F (x i ), i ∈ Z, F (x) smooth except for an isolated discontinuity θ ∈ (x m , x m+1 ). For these data Since we get (see Fig. 4) The values (S 2,2 f ) 2j +1 are displayed as × in Figs. 1-2, and the O(1) perturbations in Eq. 19 are clearly visible in Fig. 2 at the intervals adjacent to the one containing the discontinuity. After repeated application of the subdivision scheme, they cause the oscillations observed in the limit function, S ∞ 2,2 f . On the other hand, since min{|x|, |y|} ≤ |H p (x, y)| ≤ p min{|x|, |y|} for xy > 0, we may write Thus, the behavior of the S H p schemes at the intervals neighboring the singularity is closer to the behavior of the S l,r scheme, (l, r) = {(1, 2), (2, 1)}, which uses a stencil that does not cross the singularity, see Figs. 2, 3 and 4, which is, ultimately, the reason for the lack of oscillatory behavior. We would like to proceed in a similar manner in order to limit the influence of schemes with singularity-crossing stencils for the 6-point DD linear scheme. For S 3,3 we may write The stencils associated to the schemes S 2,3 and S 3,2 would not allow to avoid oscillations at all intervals neighboring a singularity in the data. On the other hand, given the distribution of the stencils for S 3,1 , S 1,3 , S 2,2 shown in Fig. 5, we conclude that it might be possible to design non-oscillatory versions of the S 3,3 scheme by considering nonlinear analogs of the linear averages involved in Eq. 20 that allow us to remain closer to the S l,r scheme whose interpolatory stencil does not cross the singularity. For the general weighted average expression we may consider the Weighted Power p mean proposed in [23].
It is proven in [23] that W p,a,b (x, y) generalizes the H p mean. Indeed, it can be easily checked that We recall next some of the properties of W p,a,b (x, y) in Eq. 22. The reader is referred to [23] for the proofs.

6-point Nonlinear, Non-oscillatory schemes
Taking into account Eqs. 20 and 13, we can write where it can easily be shown that We may obtain two families of nonlinear 6-point schemes simply by replacing each linear average by the appropriate nonlinear mean (recall that W p, 1 Because of Eq. 24-(b), these schemes remain closer to the subdivision scheme with the least oscillatory behavior, hence they are expected to display a non-oscillatory behavior similar to that of the Power p schemes, see Fig. 6. In addition, since these subdivision schemes can be written as a nonlinear perturbation of the monotone S 1,1 linear scheme, many of their properties can be analyzed with the same tools used in [19] for the Power p schemes. We examine next the polynomial reproduction properties of these families of nonlinear subdivision schemes and the existence of difference schemes.

Polynomial reproduction properties and difference schemes
Throughout this section, we denote by k the set of polynomials of degree ≤ k, and by 1 the constant sequence given by 1 i = 1, i ∈ Z.

Proposition 5
The schemes SWH p,q , SHW q,p reproduce exactly 3 .
Proof Since S 3,1 , S 2,2 and S 1,3 reproduce 3 exactly, we have that for P ∈ 3 and In the linear case, the relation between exact polynomial reproduction and the existence of the associated difference schemes S [l] with S [0] = S, and ∇ l S = S [l] ∇ l is well known [20]. As observed in [19], offset invariance [18] is the right concept to characterize the existence of difference schemes in the nonlinear case.

Definition 6 [19] A binary subdivision operator S is offset invariant (OSI) for k
if for each f ∈ l ∞ (Z) and any polynomial P (x) ∈ m , m ≤ k there exists a polynomial, Q, of degree < m such that Schemes of the form (4) with δ = ∇ k are offset invariant for k−1 . To check this, It is proven in [19] that offset invariance for k guarantees the existence of the difference schemes S [l] for l ≤ k + 1. Thus the new families of nonlinear subdivision schemes are offset invariant for 1 , which guarantees the existence of the first and second difference schemes. For the families in Eqs. 27-28, these schemes can be easily computed by elementary means. Introducing the restriction operator (m < n) and the functions L l,r : so that From Eq. 29, it can be easily deduced that ⎧ ⎨ ⎩ (SWH [2] p,q w) 2n = 2G p,q • χ n−2,n+1 w, We obtain a similar result for SHW q,p , substituting G p,q by R q,p We remark that the new schemes reproduce exactly 3 , however they are only offset invariant for 1 and, hence, difference schemes S [k] do not exist for k > 2.

Convergence
As observed in [19,22], the existence of difference schemes may be very helpful in proving the convergence of a nonlinear scheme. Since δ = ∇ 2 , so that C2 in Theorem 1 is equivalent to the following condition

Theorem 7
The schemes SWH p,q and SHW q,p are uniformly convergent, for all p, q ≥ 1.

Order of accuracy
The order of approximation of a subdivision scheme measures the approximation properties of the recursive process when applied to discrete data coming from smooth functions.

Definition 8 A convergent subdivision scheme S has approximation order r if
For a given subdivision scheme, the order of approximation after one iteration is usually much easier to obtain.

Definition 9
Let S be a subdivision scheme that satisfies Then r is called the order of approximation after one iteration of S.
Obviously, the order of approximation after one iteration of interpolatory subdivision schemes based on Lagrange interpolation is at least as high as that of the interpolatory reconstruction used in its design. We notice that Eq. 10 implies that the order of approximation after one iteration of the Power p schemes is at least 4, when refining smooth convex functions and p ≥ 2, since For the schemes defined in this paper, we can also measure how close the new schemes are to the 6-point DD scheme for smooth convex functions.

Remark 10
We know that if f i = F (ih) and F (x) is smooth and convex, (∇ 2 f ) i do not change sign. We can show, by straightforward Taylor expansions, that hence we also expect that, for smooth convex functions and h small enough, (L l,r ∇ 2 f ) i will not change sign either.
Proof Notice that Since F is a smooth function, the Taylor expansions in Eq. 39 show that x, y, z are O(h 2 ) and non-zero, provided that h is sufficiently small. Since (24) Notice that (using Taylor expansions, when necessary) Hence, assuming without loss of generality that x, z ≥ 0, Eq. 7 leads to For Z 2 (x, y, z), denoting s = H p (x, z) and using Eq. 22, we get that (s, y > 0), Moreover, Since 4 < 3q + 2 for q ≥ 1, we may conclude that from which we deduce the desired result for the schemes SWH p,q .
For the SHW q,p family, we proceed in a similar way. Assume x, y, z > 0 and write As before, it is easy to deduce that Hence, for x, y, z ≥ 0, using Eq. 22 and proceeding as in Eq. 40 we get, Thus Y 1 (x, y, z) = O(h 2p+2 ) (notice that the two terms of the 1 2 , 1 2 average in Y 1 have the same sign). For from which we deduce the desired result also for the schemes SHW q,p .

Corollary 12 Under the same conditions as in the previous proposition
As shown in [5], the order of approximation of a stable subdivision scheme can be deduced from its order of approximation after one iteration. The following result from [5] holds for linear as well as nonlinear subdivision schemes and it serves also as a motivation to study the stability of the nonlinear schemes under consideration.

Theorem 13
Let S be a convergent subdivision scheme whose approximation order after one iteration is r ≥ 1. Then if S is stable, it has approximation order r.
In [19] it is shown that the Power p schemes are stable subdivision schemes for p < 3 and unstable for p ≥ 4, hence (38) and Theorem 13 ensure that the order of approximation of the Power p schemes is 4, when refining smooth convex functions, for p < 3.
In the following section we examine the question of stability for the families of schemes (27) and Eq. 28, in order to check whether or not similar conclusions can be extracted for the new families of schemes presented in this paper.

Stability of the 6-point nonlinear schemes
Theorem 1 establishes that stability follows from two facts: Lipschitz-continuity of the operator F and contractivity of δS L N , for some L > 0. When δ = ∇ 2 and S N is offset invariant for 1 (the case of the new families of 6-point schemes), condition S2 can be equivalently expressed as follows ∃ L > 0, 0 < μ < 1 : ||(S [2] N We shall see next that the second difference schemes in Eq. 31 are defined by nonlinear functions that admit uniformly bounded Generalized Gradients. In [22], this fact was used to show the stability of a nonlinear, monotonicity preserving, scheme by expressing the contractivity condition (42) in terms of the Generalized Jacobian of the scheme and using Corollary 24 in the Appendix (or see [19]). However, we shall see that this technique does not seem to be as useful for the schemes considered in this paper.
The first step is to show that the Weighted-Power p mean (22) belongs to a special class of continuous, piecewise smooth functions: the class of C 1 pw (R 2 ) functions. Functions in this class are continuous, piecewise smooth and have uniformly bounded directional derivatives except (maybe) at 0 ∈ R m and across certain hyperplanes separating regions of C 1 smoothness. Directional derivatives along the separating hyperplanes do, also, exist. For this class of functions it is possible to define a Generalized Gradient using only the gradients on smooth regions. As the classical gradient for smooth functions, the linear map associated to any Generalized Gradient recovers all directional derivatives that 'make sense', and satisfies a chain rule property for the composition with Lipschitz curves. We refer the reader to [19] or the Appendix to this paper for the definition, and the main properties, of this class of functions.

The Weighted-Harmonic mean: Generalized Gradients
Property (24)-(a) implies that W p,a,b (x, y) is a continuous function in R 2 . It is obviously differentiable in the interior of the sectors in R 2 separated by the three hyperplanes H 1 = {x = 0}, H 2 = {y = 0}, H 3 = {x = y}. As observed in [19] (see also the Appendix) a Generalized Gradient for W p,a,b (x, y) can be defined using only the gradients in smoothness regions (see Eq. A.3), provided that certain compatibility conditions are satisfied over the separating hyperplanes. Since W p,a,b (x, y) = −W p,a,b (−x, −y), it is enough to consider the half plane y ≥ 0. Then, the compatibility conditions (A.2) amount to showing that To check these conditions, we first notice that W p,b,a (y, x) = W p,a,b (x, y), hence it is enough to consider the gradients of the 1-homogeneous 2 function (see Fig. 7) A straightforward computation leads to Since ∇φ a,b is 0-homogeneous, for y > 0, ∇φ a,b (x, y) = ∇φ a,b (x/y, y/y) = ∇φ a,b (t, 1), t = x/y, hence with 2 A function F is n-homogeneous if F (λx) = λ n F (x), being x a vector and λ a scalar. Fig. 7 W p,a,b (x, y) in its smoothness sectors Clearly, the functions μ p,a,b (t), tμ p,a,b (t), are continuous in [0, 1]. Moreover, Taking into account Fig. 7, we observe that the compatibility conditions (43), Eq. 44, Eq. 45 follow from the fact that ∀p ≥ 1, ∀a, b ≥ 0, a + b = 1, In addition, ∇φ a,b (x, y) is uniformly bounded in {(x, y) : y > x > 0} because ∇φ a,b (t, 1) is bounded for 0 ≤ t ≤ 1. Hence, ∇W p,a,b (x, y) is also uniformly bounded for (x, y) in any region of smoothness and, thus, W p,a,b (x, y) in Eq. 22 belongs to the space C 1 pw (R 2 ). From Corollary 19, W p,a,b (x, y) is a Lipschitz function. However, uniform bounds for ||DW p,a,b (x, y)|| depend on the parameters p, a, b in a more involved way than for the Power p averages. We illustrate the required computations by examining the cases p = 1, 2.
. This is an increasing function in [0, 1], hence which proves the result.

The Generalized Jacobian of the second difference schemes
Let us consider the family of schemes SWH p,q in Eq. 31. As specified in the Appendix, in order to define a Generalized Jacobian of SWH [2] p,q we need to justify the existence of uniformly bounded Generalized Gradients for the function G p,q in Eq. 30, which satisfy the chain rule (A.10).
Proof Notice that G p,q = − 1 16 W p, 3 8 Collecting all of the above we havẽ a.e.(0, 1) so that Eq. 52 provides an adequate definition of a Generalized Jacobian of G p,q .
Taking into account the bounds obtained previously for the generalized gradient of the function W p,a,b (x, y) for p = 1, 2 (these bounds are optimal for p = 2) and the known (optimal) bounds for H q (x, y), we display in Table 1 the values of the bound in Eq. 54 for 1 ≤ p, q ≤ 2.
Since (S represents any scheme in the family) ||DS [2] || ∞ = max n {||(DS [2] ) [2n,:] || 1 , ||(DS [2] ) [2n+1,:] || 1 }, from Eq. 54, Eq. 55 and the results in Table 1, we cannot ensure ||DS [2] || ∞ < 1 for any member of our family of schemes, or, equivalently, we cannot ensure the contractivity of the second difference scheme. As observed in [22], the technique described in the Appendix will be successful when the 1-norm of some row of the Generalized Jacobian is strictly uniformly bounded by 1. In this case, and by carefully considering the form of the matrix products, it may be possible to arrive at products of Generalized Jacobians whose norms are strictly bounded by 1. Taking into account the bounds in Table 1, it seems that this strategy might only be feasible for p = q = 1, because such case is the only one where Eq. 54 is strictly bounded by 1. Since the task of obtaining theoretically bounds for products of Generalized Jacobian is very involved, we examine the issue numerically in the following section.

Numerical experiments
In this section we carry out a series of numerical experiments that illustrate the theoretical developments of the previous sections. We consider first the issue of the stability of the new schemes, from a numerical perspective. In addition, we also consider the smoothness of limit functions, as well as the approximation order of the non-oscillatory 6-point schemes, comparing the numerical results with those obtained for the Power p schemes.

Stability
To examine the question of stability for each chosen subdivision scheme, S, we compute the quantities C j S (h) for each j ≥ 1, 0 < h < 1, For the computation we consider a sufficiently large set of sequences f 0 = {f 0 i } and perturbation sequences θ = {θ i }, with components randomly chosen from the set {−1, 0, 1} (hence θ ∞ = 1). We notice that if S is Lipschitz stable, C j S (h) ≤ C, Table 1 Bounds of ||D(SWH [2] p,q ) [  ∀j ≥ 1, ∀h > 0, so that any deviation with respect to this behavior may be considered a sign of the instability of the scheme. Figure 8 displays C 10 S (h) as a function of h, for S = SWH p,q and various values of p, q. The plots in Fig. 8 seem to indicate that the scheme SWH p,q is not stable for p + q > 3.
As observed in Theorem 1, Lipschitz stability follows from the contractivity of an appropriate power of the second difference scheme, which is a condition that can also be examined numerically. For h = 10 −7 , the smallest value of h considered in Fig. 8, we compute in order to check if the hypothesis S2, in its equivalent formulation (42), is fulfilled. In Fig. 9 we display the values of T j S (h) for 1 ≤ j ≤ 6. We clearly notice that ∃L ≥ 1 such that T L S (h) < 1 for S = SWH p,q , (p, q) = (1, 1), (1, 2), a behavior that would be obtained if these schemes were stable. On the other hand, T j S (h) appears to grow with j for p + q ≥ 4 indicating that condition (42) is not fulfilled.
We also observe that the T j S does not grow for S = SWH 2,1 , but it does not appear to become smaller than one, an indication that condition (42) is only a sufficient condition for stability. Table 2 summarizes the observations that might be deduced from our testing process. The same experiments were performed for the SHW q,p family, with similar results. Plot of (T [2] S (h)) j with respect to j for S = SWH p,q , h = 10 −7 and several values of p, q

Smoothness of the limit function
In Section 3.2 we have shown that the schemes proposed in this paper are convergent, i.e. S ∞ f 0 is a continuous function ∀f 0 ∈ l ∞ (Z) for S = SWH p,q , SHW q,p . The regularity of the limit function obtained from Theorem 1 (see Remark 10) seems to be much smaller than what is observed in practice. In this section we perform a numerical study of the regularity of the proposed schemes based on the following observations (see [5]): Let us assume that f (x) = S ∞ f 0 ∈ C r− , S an interpolatory subdivision scheme, and r = l + β, l ∈ N, 0 ≤ β < 1. Then, Therefore, in order to estimate the numerically regularity of a subdivision scheme, S, in a given region, [a, b], we compute (for several values of l) For our numerical testing process we consider f 0 Table 3 we display the results corresponding to an initial mesh with N = 18 points in [−6, 6], so that x = 0 does not belong to the initial mesh, and in Table 4 we display the results obtained when a uniform mesh with N = 21 points in [−6, 6], so that x = 0 belongs to the initial grid, is used to compute f 0 (see also Fig. 10).
The tables show that, for these examples, the global regularity of the limit function is at least r = 1. Moreover, when x = 0 (the abscissa of the maximum of F (x)) is included in the initial grid, the global regularity of the limit function obtained with the Power p schemes is smaller than that obtained with the new schemes when max{p, q} > 1. According to Table 3, in this case the limit function seems to be globally C 1 for max{p, q} > 1. In addition, we also observe in both tables that for (p, q) = (2, 2) we get the same smoothness as for the S 3,3 scheme, around x = 0, much higher than that of the Power 2 scheme.  3]. Initial data and limit functions displayed in Fig. 10, left column

Approximation order
The order of approximation of a subdivision scheme measures its ability to reconstruct smooth functions from relatively coarse samples.
is a smooth function, we study the difference between the limit function f ∞ (x) = S ∞ f 0 (x) and F (x) in a given region by considering S ∞ f 0 ≈ S L f 0 (with L = 7 in all test cases) and measuring In Tables 5-8

we display E S,[a,b] (h) for different values of h, different regions [a, b]
and different functions F (x). The tables also show the numerical order of accuracy, r n , obtained by a least squares fit of the data (log 2 (h l ), log 2 (E S (h l )) for h l = h 0 /2 l , and a given value of h 0 . The ultimate purpose of the numerical testing process is twofold. On one hand the Tables show that, by choosing p, q appropriately, it is possible to obtain the same order of approximation (as well as errors of a similar magnitude) as that of the S 3,3 scheme. In addition, the Tables also show that, for each scheme in the family, the order of approximation of the limit function is the same as the theoretical order of approximation after one iteration, r t in the Tables. For the new schemes proposed in the paper, r t is obtained in Corollary 12 for convex regions, but it can also be Table 4 Same as Table 3. Initial data and limit functions displayed in Fig. 10, right column which holds for smooth functions, provided that (L l,r •∇ 2 f ) 2j +1 have the same sign. We notice further that if F is a smooth function and f j = F (x j ), straightforward Taylor expansions lead to Hence, at smooth convex regions (at least for h small enough) L l,r • ∇ 2 does not change sign. Moreover, if x n = nh is an inflection point and F (x n ) = 0, the formulas above show that, for h small enough and (l, r) = (1, 3), (3,1), (2,2), (L l,r ∇ 2 f ) 2n+m do not change sign for m ≥ 0 or m < 0, and the calculation of r t by Taylor expansions is feasible. This is the case of the two functions considered in the numerical tests. The testing process below is carried out for the family of schemes SWH p,q . We have also performed the same study for the SHW q,p family. As expected, the resulting tables are similar and the conclusions are also the same, hence we omit them. We include the errors corresponding to the S 3,3 and S H p schemes for the sake of comparison. An * in the Table 6 means that it is not possible to find r t by Taylor expansions.

Gaussian data
, h 0 = 0.1. Table 5  We observe that the estimated order of approximation of SWH p,q coincides with r t , the order of approximation after one iteration in the convex region.
In Table 6 we display the corresponding results for [a, b] = [−1, −0.3], which contains the inflection point x = 0.5. The Table indicates that the computed r n coincides with r t for all the nonlinear 6-point schemes. As observed above, r t can still be computed by Taylor expansions for the 6-points nonlinear schemes. It should be noticed that for F (x) = e −2x 2 which has a vertical asymptote in x = ±0.5 and hence it is unbounded around the inflection points. Since F ((2n + 1)h/2) = F (nh)h/2 + O(h 2 ), Eq. 60 leads to around the inflection point, i.e., the order of approximation after one iteration is 5 in the non-convex region. Performing the analogous computation for (p, q) = (3, 1),  (3, 2) we find that the 'theoretical' approximation order after one iteration is 5 and 6, respectively, also in the non-convex region.
We also remark that, in all the tables displayed, the magnitude of the errors corresponding to the 6-point nonlinear schemes whose order of accuracy is 4, 5 or 6 are similar to those of S 3,3 and better than that of the Power p schemes.

Tangent data
We repeat the previous study for F (x) = tan(π x). In this case, the function is convex in the interval [0.1, 0.3], and changes convexity at [−0.25, 0.25]. We consider h 0 = 0.1. The conclusions are similar. In particular, the magnitude of the errors is similar to those of S 3,3 , and better than those obtained with the Power p schemes, for SWH p,q schemes whose order of accuracy is 4, 5 or 6. We also remark that order of approximation of the scheme coincides with r t (not displayed in the tables). Notice that for F (x) = tan(π x) F I V (x)/F (x) = 8π 6 (cos(2πx) − 5) 2 tan(π x) sec 6 (π x) which is a bounded function around x = 0. Hence, according to Eq. 60, the order of approximation after one iteration of SWH p,q dd is 6 also in the non-convex region.

Conclusions
We have constructed two families of non-oscillatory subdivision schemes that can be considered nonlinear versions of the 6-point Deslauries-Dubuc interpolatory subdivision scheme. We have studied their convergence by exploiting the (piecewise) smoothness properties of the functions that define these subdivision schemes, following a novel technique developed in [19]. The stability of these schemes turns out to be more difficult to study with the techniques employed in [22] and we have explored this issue computationally. The numerical results reveal indeed that the techniques based on finding appropriate bounds for the Generalized Jacobian of the second difference scheme, as in [22], have no chance to succeed, except for (p, q) = (1, 1) (where contractivity might be proven for L = 2) and (p, q) = (2, 1) (for L = 4).
We have also performed several numerical experiments that suggest that the approximation properties of the new schemes can be as good as those of the 6-point linear scheme when reconstructing smooth functions. In addition, numerical experiments show that the smoothness of the limit functions obtained from convex data may be larger than the smoothness of the limit functions obtained with Power p schemes. which we shall denote by j . Therefore R m = ( j j ) ∪ ( i H i ), the sets j j and i H i are disjoint, and ∂ j , the boundary of j , is always included in the union of the separating hyperplanes.
In addition, functions in C 1 pw (R m ) have uniformly bounded gradients in smooth regions, i.e.
and the smooth gradients satisfy the following compatibility condition on the separating hyperplanes: Let 0 = x ∈ ∂ ⊂ H, where is a smoothness region for ψ and H is one of the hyperplanes separating the regions of smoothness of ψ. Then where D w ψ(x) is the derivative of ψ at x in the direction of w. Functions in C 1 pw (R m ) admit Generalized Gradients. As the standard gradient of a smooth function, the main property of a Generalized Gradient is that the associated linear map recovers all directional derivatives that 'make sense' at a given point.
provides an adequate definition of a Generalized Gradient of ψ ∈ C 1 pw (R m ), since for each 0 = x ∈ R m , Dψ(x) defines a linear map that satisfies for any 0 = v ∈ R m when x belongs to a smoothness region, and also for any Notice that Dψ(x) in Eq. A.3 might not be uniquely defined when x belongs to a hyperplane separating two or more regions of smoothness, if the limit in Eq. A.3 is different for different smoothness regions with a common boundary. However the compatibility condition (A.2) ensures (A.4) for all directional derivatives that make sense, independently of the chosen definition for the vector Dψ(x).
The following results generalize some of the properties satisfied by the gradient of a smooth function.
where Dψ is a generalized gradient of ψ. Then, hence, under the hypothesis of the Lemma, γ ((0, 1)) =: ⊂ or ⊂ H − {0}, and for any Generalized Gradient Dψ and for any t ∈ (0, 1) (notice that if x, y ∈ H, a separating hyperplane, then x − y ∈ H). Since g(t) is differentiable in (0, 1), by the classical Mean Value Theorem (MTV) Proof We consider again the straight line γ (t) = y + t (x − y) and the function g(t) = ψ(γ (t)) in Eq. A.5. Notice that γ (t) can either cut the separating hyperplanes at a finite number of points, or belong entirely to one of them. We prove (A.9) in each case.
The following result establishes that the chain rule holds for the composition of C 1 pw (R m ) functions with Lipschitz curves.
The chain rule (A.10) turns out to be the key property for our purposes. The next result shows that this property is also satisfied by a possibly wider class of functions. Notice that ||D(ψ • M)(x)|| is also uniformly bounded. Hence, we may also associate the notion of Generalized Gradients to certain functions, related (by composition) to C 1 pw (R m ) functions.

Generalized Jacobians
The notion of Generalized Gradient leads, in a rather natural way, to that of Generalized Jacobian for functions ψ : R m → R p , ψ = (ψ 1 , . . . , ψ p ), ψ i ∈ C 1 pw (R m ). As in the smooth case, . . .

Dψ p (x)
⎞ ⎟ ⎠ provides the definition of a Generalized Jacobian of ψ at x ∈ R m . Obviously, the definition may not be unique, but, as stated below, the linear maps associated to such matrices also satisfy the chain rule for the composition with Lipschitz curves.
Since the countable union of sets of zero measure has zero measure, Eq. A.18 follows as in Theorem 22.
Obviously, the last two theorems apply as long as the functions ψ i admit uniformly bounded Generalized Gradients, Dψ i , satisfying the chain rule (A.10).
Hence, since γ j (1) = S j f , γ j (0) = S j g, ∀ j ≥ 0 we can write From Eq. A.20, we easily deduce the following contractivity result, which has been used in [22] to prove the stability of a monotone nonlinear scheme.