On a regularization-correction approach for the approximation of piecewise smooth functions
CCorrected subdivision approximation ofpiecewise smooth functions
Sergio Amat, ∗ David Levin, † Juan Ruiz- ´Alvarez ‡ Abstract
Subdivision schemes are useful mathematical tools for the gene-ration of curves and surfaces. The linear approaches suffer from Gibbsoscillations when approximating functions with singularities. On theother hand, when we analyze the convergence of nonlinear subdivisionschemes the regularity of the limit function is smaller than in the linearcase. The goal of this paper is to introduce a corrected implementationof linear interpolatory subdivision schemes addressing both propertiesat the same time: regularity and adaption to singularities with highorder of accuracy. In both cases of point-value data and of cell-averagedata, we are able to construct a subdivision-based algorithm produc-ing approximations with high precision and limit functions with highpiecewise regularity and without diffusion nor oscillations in the pres-ence of discontinuities.
Key Words.
Non-smooth approximation, subdivision schemes, regu-larity, Gibbs phenomenon, diffusion. This work was funded by project 20928/PI/18 (Proyecto financiado por la ComunidadAut´onoma de la Regi´on de Murcia a trav´es de la convocatoria de Ayudas a proyectos parael desarrollo de investigaci´on cient´ıfica y t´ecnica por grupos competitivos, incluida en elPrograma Regional de Fomento de la Investigaci´on Cient´ıfica y T´ecnica (Plan de Actuaci´on2018) de la Fundaci´on S´eneca-Agencia de Ciencia y Tecnolog´ıa de la Regi´on de Murcia)and by the Spanish national research project PID2019-108336GB-I00. ∗ Departamento de Matem´atica Aplicada y Estad´ıstica. Universidad Polit´ecnica deCartagena (Spain). e-mail: [email protected] † School of Mathematical Sciences. Tel-Aviv University, Tel-Aviv (Israel). e-mail: [email protected] ‡ Departamento de Matem´atica Aplicada y Estad´ıstica. Universidad Polit´ecnica deCartagena (Spain). e-mail: [email protected] a r X i v : . [ m a t h . NA ] D ec Introduction
Subdivision schemes are very well known mathematical tools for the ge-neration of curves and surfaces. In the past years many articles have beenpublished about linear subdivision schemes. An incomplete list of referencesis, for example, [12, 16, 19, 20, 21, 24, 32, 33, 34, 35]. Non linear subdivisionschemes have also attracted much attention. See for example, [5, 7, 8, 13, 27,28, 29, 30] and the references therein for an incomplete list of publicationsabout the subject.Several properties of subdivision schemes are desirable for practical a-pplications: convergence, regularity of the limit function, stability, order ofapproximation and adaptation to singularities. The last property usually re-quires the implementation of nonlinear strategies, resulting in the so callednonlinear subdivision schemes. These schemes present a regularity of thelimit functions that is smaller than the one obtained by linear schemes. Inpractice it is usually desired that the schemes have at least C regularity, asthis kind of schemes are associated to improved results in the aerodynam-ics or hydrodynamics of the generated surface. When the data is computedwith high accuracy, interpolatory subdivision schemes are preferred over ap-proximating schemes, since using the former, the original data is preserved.These schemes are usually based on Lagrange interpolation using a centeredstencil.The ENO (Essential Non-Oscillatory) reconstruction is an approxima-tion based on a nonlinear operator with the same order of accuracy as thelinear scheme, but capable of reducing the effect of isolated singularities. Ituses a selection of the stencil that aims to avoid zones affected by singu-larities. Thus, the precision is only reduced at the intervals which containa singularity. The aim of the ENO subcell resolution (ENO-SR) techniqueis to improve the approximation properties of the reconstruction even atintervals that contain a singularity. Therefore, ENO-SR reconstruction hasan improved order of accuracy over the linear reconstruction for piecewisesmooth functions. In [9], the authors proposed a rigorous analysis of theENO-SR procedure (both for the point-value and the cell-average sampling).In this article we give the details that are necessary to describe the strategyintroduced in [9]. In [9], for a piecewise C function f with a jump [ f (cid:48) ] ofthe first derivative at x ∗ , the authors proved that the singularity is alwaysdetected for a grid-spacing discretization h smaller than the critical scale h c being, h c := | [ f (cid:48) ] | t ∈ R \{ x ∗ } | f (cid:48)(cid:48) ( t ) | . (1)2his critical scale represents the minimal level of resolution needed to dis-tinguish between the singularity and a smooth region.The detection algorithm for the point-values sampling case in [9] candetect corner and jump discontinuities. Although, only corner singularitiescan be located using this discretization, as the location of jumps in thefunction is lost during the discretization process. Jumps in the function canbe located if we use other kinds of discretization that accumulate informationover the whole interval of discretization instead of sampling information atisolate locations. One example is the cell-average discretization, that weintroduce later on in Section 2. In [9], the authors also perform an analysisfor jump discontinuities and data discretized by cell-averages, obtaining thecritical scale, h c := | [ f ] | t ∈ R \{ x ∗ } | f (cid:48) ( t ) | , (2)where [ f ] is the jump in the function f .The algorithm for the point-values sampling is based on the comparisonof second order differences for the detection of intervals that are potentialcandidates of containing a singularity:∆ f j = f ( x j +1 ) − f ( x j ) + f ( x j − ) . (3)Then, in [9], the intervals are labeled as B if they are suspicious of containinga singularity or G if they are not. The rules used for this process are:1. If | ∆ h f ( x j − ) | > | ∆ h f ( x j − ± nh ) | , n = 1 , · · · , m, the intervals I j − = [ x j − , x j − ] and I j = [ x j − , x j ] are labeled as B .In this case we do not know which of the two intervals contain thesingularity, so both are labeled as B .2. If | ∆ h f ( x j ) | > | ∆ h f ( x j + nh ) | , n = 1 , · · · , m − , and | ∆ h f ( x j − ) | > | ∆ h f ( x j − − nh ) | , n = 1 , · · · , m − , then the interval I j is labeled as B . In this case, the two largest divideddifferences contain the interval I j , so it is a candidate of containingthe singularity. 3t is assumed that the rest of the intervals do not contain a singularity. In[9] it is also proved that under the critical scale h c presented in (1) for thepoint-values sampling and (2) for the cell-averages sampling, any intervalthat contains a singularity is labeled as B and all the intervals labeled as G are at smooth zones of the function. Even though, the algorithm canproduce bad detections, i.e. intervals labeled as B at smooth regions of thefunction. Following [9], we can assume that given a continuous piecewisefunctions there exits a critical scale h c depending only on derivatives of f at smooth regions such that all the intervals with a singularity are labeledas B or BB .Once all the intervals have been labeled as G or B , in order to locatethe position of corner singularities, we build cubic polynomials to the leftand to the right of the suspicious interval, labeled as B or BB . Then, wedefine a function H as the difference between the two polynomials, and wesuppose that there is a unique root of this function inside the suspiciousinterval. Note that for a sufficiently small grid-spacing, only one root of H can fall inside the considered interval. In Lemma 3, statement 2 of [9]the authors prove this fact. Thus, a good approximation for the positionof the discontinuity can be easily obtained finding the roots of H . Theinterested reader can refer to [10] for a nice discussion about the process.The accuracy of the results depends on the order of the polynomials, as theauthors prove in Lemma 3, statement 3 of [9]. In the cell-average setting,jumps in the function transform into corner singularities in the primitiveand their position can be easily obtained following the process describedabove.In this paper we include one step more for the location algorithm whena BB cell is detected. In this case, once we have located the discontinuityusing the process described, we know which of the two contiguous B cellscontains the discontinuity. Then we treat this cell as B cell and repeatthe process of location in order to assure the best possible accuracy in thelocation of the singularity.In [11], a specific prediction operator in the interpolatory framework wasproposed and analyzed. It is oriented towards the representation of piece-wise smooth continuous functions. Given a family of singularity points, aposition-dependent prediction operator is defined. Its convergence was es-tablished using the matrix formalism. The algorithm lead to the successfulcontrol of the classical Gibbs phenomenon associated to the approximation oflocally discontinuous functions. This position-dependent algorithm is equiv-alent to the ENO-SR subdivision scheme assuming the previous knowledgeof the singularity positions, in particular improving the accuracy of the linear4pproach.Despite their good accuracy properties, the theoretical regularities ofENO and ENO-SR schemes are C − in the point-values, that is smaller thanthe regularity of the 4-point linear subdivision scheme, that is C − , as wecan see in [13]. Other adaptive interpolation schemes like the PPH algorithm[2] manage to eliminate the Gibbs phenomenon close to the discontinuities,but introducing diffusion and, thus, reducing the order of accuracy but alsothe regularity (that is C − numerically) close to the discontinuities.The goal of this paper is to develop a corrected implementation of thelinear subdivision scheme that combines adaption to singularities and, atthe same time, high order of regularity avoiding diffusion and Gibbs phe-nomenon. Indeed, we present a procedure to develop corrected subdivisionapproximations for functions with discontinuitis which have the five impor-tant properties:1. Interpolation.2. High precision.3. High piecewise regularity.4. No diffusion.5. No oscilationsAs far as we know, this is the first time that subdivision schemes that ownall these properties at the same time appear in the literature. We present theparticular case of the 4-point Dubuc-Deslauriers interpolatory subdivisionscheme [15] through which we obtain a C − piecewise regular limit functionand that is capable of reproducing piecewise cubic polynomials. Indeed, itis straightforward to obtain higher regularity and reproduction of piecewisepolynomials of higher degree, just using the same technique with largerstencils. It is also possible and straightforward to extend the correctionconstruction presented in this article to other subdivisions schemes aimedto deal with data that contain singularities.The paper is organized as follows: in Section 2 we introduce the dis-cretizations that is being used in the paper, Section 3 is dedicated to intro-duce our corrected approximation approach, and the analysis of the approx-imation order of the resulting approximants. Section 4 is dedicated to checkthat the numerical regularity of the limit function in the point-values andthe cell-average discretization corresponds to the one found theoretically. InSection 5 we present some experiments for univariate and bivariate functionsdiscretized by point-values’ sampling, in Section 6 we show some experimentsfor univariate and bivariate functions discretized by cell-averages’ samplingand, finally, in Section 7 we expose the conclusions.5 The discretizations: point-values and cell averages
In the case of a point-values discretization we can just consider that we aregiven a vector of values on a grid. This means that our data is interpretedas the sampling of a function at grid-points { x j = jh } . In this case, as it hasbeen mentioned above, the position of the discontinuities in the function islost during the discretization process and there is no hope to recover theirexact position. Other kind of singularities, such as discontinuities in the firstderivative, can be located using the point-values discretization. Anotheroption is to consider the original data as averages of a certain piecewisecontinuous function over the intervals of a grid. This is the cell-averagesetting, which also allows to locate jump discontinuities in the function. Inboth cases we consider the grid points in [0 , X = { x j } Nj =0 , x j = jh, h = 1 N .
For the point-values case we use the discretization { f j = f ( x j ) } Nj =0 at thedata points { x j = jh } Nj =0 .On the other hand, for the cell-averages sampling we are given the localaverages values, ¯ f j = 1 h (cid:90) x j x j − f ( x ) dx, j = 1 , · · · , N. (4)Also in this case we aim at approximating the underlying function f .Let us define the sequence { F j } as, F j = h j (cid:88) i =1 ¯ f i = (cid:90) x j f ( y ) dy, j = 1 , · · · , N, (5)taking F = 0. Denoting by F the primitive function of f , i.e. F ( x ) = (cid:82) x f ( y ) dy , the values { F j } are the point-values discretization of F . Nowwe are back in the case of point-value data, for F , and after finding anapproximation G ( x ) to F ( x ), g ( x ) = G (cid:48) ( x ) would be the approximation to f ( x ) = F (cid:48) ( x ), such that, ¯ g j = F j − F j − h . (6)In the following sections we suggest the framework of a corrected subdi-vision algorithm that is adapted to the discontinuities, that avoids the Gibbsphenomenon, that attains the same regularity at smooth zones as the equiv-alent linear algorithm, without diffusion and that reproduces polynomialsof the degree associated to the linear scheme used.6 A corrected implementation of linear subdivi-sion schemes
We present our approach for the approximation of a function with one sin-gular point. Later on we explain how to use it for the case of several singularpoints.Let f be a piecewise C -smooth function on [0 , x ∗ , and assume we are given the vector of values { f ( x i ) } Ni =0 at the datapoints { x i = ih } Ni =0 , N = 1 + 1 /h . We denote by f − ( x ) and f + ( x ) thefunctions to the left and to the right of the discontinuity respectively.In our framework we are going to use the 4-point Dubuc-Deslaurierssubdivision scheme: (cid:26) ( Sf k ) j = f kj , ( Sf k ) j +1 = − f kj − + f kj + f kj +1 − f kj +2 . (7)This scheme has fourth order of accuracy for C functions, and it has a C − regularity [15], [17]. We aim at retrieving these properties for functions withjump singularities in the function and the first derivative for data discretizedby point-values or by cell-averages. In this subsection we introduce the corrected approximation algorithm. Let T − ( x ) and T +3 ( x ) be the third order Taylor approximations of f − and of f + at x ∗ respectively. Consider the one sided cubic polynomial (cid:26) T + ( x ) = 0 , x < x ∗ ,T + ( x ) = T +3 ( x ) − T − ( x ) , x ≥ x ∗ . (8)For x ≥ x ∗ T + ( x ) = [ f ] + [ f (cid:48) ]( x − x ∗ ) + 12 [ f (cid:48)(cid:48) ]( x − x ∗ ) + 16 [ f (cid:48)(cid:48)(cid:48) ]( x − x ∗ ) , (9)where [ f ] , [ f (cid:48) ] , [ f (cid:48)(cid:48) ] , [ f (cid:48)(cid:48)(cid:48) ] denote the jumps in the derivatives of f at x ∗ :[ f ] = f + ( x ∗ ) − f − ( x ∗ ) , [ f (cid:48) ] = ( f + ) (cid:48) ( x ∗ ) − ( f − ) (cid:48) ( x ∗ ) , [ f (cid:48)(cid:48) ] = ( f + ) (cid:48)(cid:48) ( x ∗ ) − ( f − ) (cid:48)(cid:48) ( x ∗ ) and [ f (cid:48)(cid:48)(cid:48) ] = ( f + ) (cid:48)(cid:48)(cid:48) ( x ∗ ) − ( f − ) (cid:48)(cid:48)(cid:48) ( x ∗ ) . It follows that g ≡ f − T + ∈ C [0 , , f at x ∗ , we will use insteada fourth order approximation of T + ( x ): For x ≥ x ∗ , (cid:101) T + ( x ) = (cid:102) [ f ] + (cid:103) [ f (cid:48) ]( x − x ∗ ) + 12 (cid:103) [ f (cid:48)(cid:48) ]( x − x ∗ ) + 16 (cid:103) [ f (cid:48)(cid:48)(cid:48) ]( x − x ∗ ) , (10)where (cid:102) [ f ] (cid:102) [ f (cid:48) ] (cid:103) [ f (cid:48)(cid:48) ] (cid:103) [ f (cid:48)(cid:48)(cid:48) ] = [ f ][ f (cid:48) ][ f (cid:48)(cid:48) ][ f (cid:48)(cid:48)(cid:48) ] + O ( h ) O ( h ) O ( h ) O ( h ) . (11)Note that near x = x ∗ , i.e., if | x − x ∗ | = O ( h ), | T + ( x ) − ˜ T + ( x ) | = O ( h ) , as h → . (12)Our algorithm has three steps: • Step : Smoothing the data.We compute the new data using ˜ T + ( x ) in (10),˜ g ( x i ) = f ( x i ) − ˜ T + ( x i ) , i = 0 , ..., N. (13)It is clear that g ( x ) does not present singularities up to the thirdderivative. As we show below, using a fourth oder approximation ˜ g to g , implies a truncation error that is O ( h ). • Step : Subdivision.We apply the 4-point interpolatory subdivision scheme to the data˜ g ( x j ) and we denote the limit function by ˜ g ∞ . • Step : Correcting the approximation.˜ f ∞ ( x ) = ˜ g ∞ ( x ) + (cid:101) T + ( x ) (14)In Figure 1 we present an example with one of the functions that weuse in the numerical experiments (more specifically, the one in (25)). Thecircles represent a function discretized by point-values with a discontinuityin the first derivative. The dots represent the smoothed data in (13), and8e can see that it owns better regularity properties than the original data.Using the strategy presented in this section, we apply a subdivision schemeto the corrected data and then we compute the corrected approximation in(14). Original dataOriginal data plus correction term
Figure 1: The circles in this graphs correspond to the original data obtainedfrom the discretization of the function in (25). The dots correspond to thecorrected data in (13).
Remark 1.
It is important to note here that the algorithm can be appliedto functions with m > discontinuities. In this case the correction can beapplied iteratively from left to right as the discontinuities are found. Thetotal correction term would be, ˜ T total + ( x ) = m (cid:88) n =1 ˜ T n + ( x ) , being (cid:101) T n + ( x ) = (cid:103) [ f ] n + (cid:103) [ f (cid:48) ] n ( x − x ∗ n ) + 12 (cid:103) [ f (cid:48)(cid:48) ] n ( x − x ∗ n ) + 16 (cid:103) [ f (cid:48)(cid:48)(cid:48) ] n ( x − x ∗ n ) , (15) the correction term (10) at each of the m discontinuities placed at x ∗ n , n =1 , · · · , m , found in the data. In this case, ˜ T total + ( x ) should be used insteadon ˜ T + ( x ) in steps 1 and 3 of the algorithm. In order to work with data discretized by cell-averages, { ¯ f j } , we firstlyobtain the point-values { F j } of the primitive function F (5). Then we apply9ur strategy for the point-values discretization described above to obtain anapproximation G ∼ F and then approximate f by G (cid:48) . The main differencefrom the case of point-values data is that the function F is a continuousfunction, with a jump in its first derivative. Thus, an additional step inthe algorithm for cell-averages data is to define x ∗ as the intersection pointof the two local polynomials derived from the data { F j } . In Figure 2 werepresent this process. Original cell average data
PrimitivePrimitive plus correction term
Subdivided cell average data
Figure 2: Left, primitive of the function presented in (30). Center, primitiveof the data in the plot to the left, and the corrected primitive. Right, resultof the subdivision of the function discretized by cell-averages.In the next section we present the construction of a fourth order accurateapproximation of the jumps in the function and its derivatives for datadiscretized by point-values.
The approximation of the jumps with the desired accuracy can be obtainedusing the strategy derived by us in [1]. In order to obtain the approximatejumps in the derivatives of f we need to know the position of the discon-tinuity x ∗ up to certain accuracy. In the introduction we mentioned howto obtain an approximation of the position of the discontinuity x ∗ with thedesired accuracy following [9]. If we work with stencils of four points, weneed O ( h ) accuracy in the detection. Let’s suppose that we know that thediscontinuity is placed at a distance α from x j in the interval { x j , x j +1 } .If we want O ( h ) accuracy for the approximation of the jump relations, weneed to use four points to each side of the discontinuity. Let’s suppose thatwe use the stencil { f + j − , f + j − , f + j − , f + j , f − j +1 , f − j +2 , f − j +3 , f − j +4 } placed at the10ositions { x j − , x j − , x j − , x j , x j +1 , x j +2 , x j +3 , x j +4 } . Then, we can appro-ximate the values of the derivatives of f from both sides of the discontinuityjust using third order Taylor expansion around x ∗ and setting the followingsystem of equations for the − side, f − j = f − ( x ∗ ) − f − x ( x ∗ ) α + 12 f − xx ( x ∗ ) α − f − xxx ( x ∗ ) α ,f − j − = f − ( x ∗ ) − f − x ( x ∗ )( h + α ) + 12 f − xx ( x ∗ )( h + α ) − f − xxx ( x ∗ )( h + α ) ,f − j − = f − ( x ∗ ) − f − x ( x ∗ )(2 h + α ) + 12 f − xx ( x ∗ )(2 h + α ) − f − xxx ( x ∗ )(2 h + α ) ,f − j − = f − ( x ∗ ) − f − x ( x ∗ )(3 h + α ) + 12 f − xx ( x ∗ )(3 h + α ) − f − xxx ( x ∗ )(3 h + α ) , (16)and for the + side, f + j +1 = f + ( x ∗ ) + f + x ( x ∗ )( h − α ) + 12 f + xx ( x ∗ )( h − α ) + 13! f + xxx ( x ∗ )( h − α ) ,f + j +2 = f + ( x ∗ ) + f + x ( x ∗ )(2 h − α ) + 12 f + xx ( x ∗ )(2 h − α ) + 13! f + xxx ( x ∗ )(2 h − α ) ,f + j +3 = f + ( x ∗ ) + f − x ( x ∗ )(3 h − α ) + 12 f + xx ( x ∗ )(3 h − α ) + 13! f + xxx ( x ∗ )(3 h − α ) ,f + j +4 = f + ( x ∗ ) + f + x ( x ∗ )(4 h − α ) + 12 f + xx ( x ∗ )(4 h − α ) + 13! f + xxx ( x ∗ )(4 h − α ) , (17)being α = x ∗ − x j . It is clear that the previous systems can be writtenin matrix form where the system matrix is a Vandermonde matrix, that isalways invertible.Solving the two systems (16) and (17), where f − ( x ∗ ) , f − x ( x ∗ ) , f − xx ( x ∗ ) , f − xxx ( x ∗ )and f + ( x ∗ ) , f + x ( x ∗ ) , f + xx ( x ∗ ) , f + xxx ( x ∗ ) are the unknowns, we can obtain appro-ximations for the jumps in the derivatives of f .It is easy to prove [1] that (cid:102) [ f ] (cid:102) [ f (cid:48) ] (cid:103) [ f (cid:48)(cid:48) ] (cid:103) [ f (cid:48)(cid:48)(cid:48) ] = f + ( x ∗ ) − f − ( x ∗ ) f + x ( x ∗ ) − f − x ( x ∗ ) f + xx ( x ∗ ) − f − xx ( x ∗ ) f + xxx ( x ∗ ) − f − xxx ( x ∗ ) + O ( h ) O ( h ) O ( h ) O ( h ) . (18) Let us analyze the approximation properties of the corrected scheme for thecase of point-values data and for cell-averages data.11t is important to note that the technique described in this and previoussections can also be applied to treat the approximation near the boundarypoints 0 and 1. Outside the boundary we use a zero padding strategy, andwe treat each boundary as a discontinuity which position is known. Thus,treating the boundaries is just considering more than one discontinuity inthe data, as explained in Remark 1. Using this strategy enables us to extendthe approximation results discussed below to the entire interval [0 , x j , x j +1 ). For the point-values discretization, we will considertwo cases: when we find jumps in the function or in the first derivative.Let’s start by the former.We are working with jump discontinuities in the point-values so, strictly,we can not know the exact location of the discontinuity. Thus, we can onlyhope to obtain a limit function that contains a jump discontinuity in thesame interval as the original function. As the exact location of the discon-tinuity is an information that can not be obtained by any means from thedata, we can just assume that the approximated location of the discontinuitycorresponds with the exact location of the discontinuity. Theorem 3.1. - The point-values’ case: jump discontinuities.
Let f be a piecewise C -smooth function on [0 , , with a jump discontinuity,and assume we are given the vector of values { f j } Nj =0 at the data points { x j = jh } Nj =0 , N = 1 + 1 /h . Let’s choose the singular point to be at x ∗ ∈ ( x j , x j +1 ) . Then the corrected approximation obtained by the correc-ted algorithm interpolates the data points, it is C − { [0 , x ∗ ) ∪ ( x ∗ , } , and || f − ˜ f ∞ || ∞ , [0 , \{ x ∗ } = O ( h ) as h → , for any piecewise C -smooth f thathas a singular point at x ∗ and that has as point-values { f ( x j ) } Nj =0 = { f j } Nj =0 .Proof. The 4-point subdivision scheme is an interpolatory scheme and thelimit function it defines is C − [15]. Hence,˜ g ∞ ( x j ) = f ( x j ) − ˜ T + ( x j ) , and ˜ g ∞ ∈ C − [0 , f and the smoothness property.To prove the approximation order, we use the fact that the subdivisionscheme is a local procedure and that it reproduces cubic polynomials. Fur-thermore, by [15], [17], for f ∈ C ( R ) the 4-point subdivision scheme gives O ( h ) approximation order to f . Therefore, away from the singularity, for x < x j − and for x > x j +3 , | ˜ g ∞ − ( f − ˜ T + ) | = O ( h ) , h →
0. Correcting the approximation by ˜ T + yields the approximationresult away from x ∗ .To show the approximation order near x ∗ , we rewrite the data for thesubdivision process, ˜ g ( x j ) = f ( x j ) − ˜ T + ( x j ), as˜ g ( x j ) = f ( x j ) − T + ( x j ) − T − ( x j ) + ( T + ( x j ) − ˜ T + ( x j )) + T − ( x j ) . Denoting q = f − T + − T − , we observe that q = f − T − for x < x ∗ and q = f − T +3 for x > x ∗ . In both cases q ( x ) = O ( h ) near x ∗ . Next we notethat by (11) ( T + ( x i ) − ˜ T + ( x i )) = O ( h ) , for i − ≤ j ≤ i + 5. Using the above, plus the locality of the subdivisionscheme and the reproduction of cubic polynomials, it follows that˜ g ∞ ( x ) = T − ( x ) + O ( h ) , x j − ≤ x ≤ x j +3 . Finally, in view of (8) and (12) we obtain˜ f ∞ ( x ) = ˜ g ∞ ( x ) + ˜ T + ( x ) = (cid:26) T − ( x ) + O ( h ) , x j − ≤ x ≤ x ∗ ,T +3 ( x ) + O ( h ) , x ∗ ≤ x ≤ x j +3 , (19)implying that˜ f ∞ ( x ) = f ( x ) + O ( h ) , x j − ≤ x ≤ x j +3 , x (cid:54) = x ∗ , which completes the proof. Remark 2.
In order to include the point x ∗ in Theorem 3.1, we need tosuppose that the initial data comes from a function that is continuous fromthe left or from the right of x ∗ , as this information, as well as the exactposition of the discontinuity is lost in the discretization process. Now we can continue with the case when we find corner singularities, i.e.discontinuities in the first derivative.
Theorem 3.2. - The point-values’ case: discontinuities in the firstderivative.
Let f be a piecewise C -smooth function on [0 , , with a jumpdiscontinuity in the first derivative, and assume we are given the vectorof values { f j } Nj =0 at the data points { x j = jh } Nj =0 , N = 1 + 1 /h . Let’ssuppose that the singular point is placed at s ∗ ∈ ( x j , x j +1 ) and that we havelocated the singularity at the approximated location x ∗ ∈ ( x j , x j +1 ) . Then,the corrected approximation obtained by the corrected algorithm interpolatesthe data points, it is C − { [0 , x ∗ ) ∪ ( x ∗ , } , and || f − ˜ f ∞ || ∞ = O ( h ) as h → . roof. Outside the interval ( x ∗ , s ∗ ) (if the approximated location x ∗ is tothe right of s ∗ , the analysis is equivalent), the proof is similar to the onefollowed in Theorem 3.1.Let’s now analyze the case when we subdivide in the interval ( x ∗ , s ∗ ).We can just follow the proof of Theorem 1 in [9]. We suppose that thediscontinuity is placed in the interval [ x j , x j +1 ], being h = x j +1 − x j auniform grid-spacing. A graph of this case is plotted in Figure, 3. As we areusing cubic polynomials for the location of the discontinuity, it is clear thatthe distance between x ∗ and s ∗ must be O ( h ), as proved in statement 3 ofLemma 3 of [9]. Let us analyze the accuracy attained by the subdivisionscheme in the point-values discretization in the interval ( x ∗ , s ∗ ). The left sideof the discontinuity placed at s ∗ is labeled as the − side and the right side of s ∗ as the + side. The approximated location of the discontinuity is labeled as x ∗ . Following the process described in Section 3.3 we obtain the jumps in thefunction and its derivatives at x ∗ . It is clear that if the piecewise function isnot composed of polynomials of degree smaller or equal than 3, the locationof the discontinuity will be approximated. In this case, obtaining the jumpin the function and its derivatives at x ∗ is just as extending f + from s ∗ upto x ∗ , as shown in Figure 3. Let’s write the limit function obtained by thecorrected algorithm as,˜ f ∞ ( x ) = (cid:26) ( ˜ f ∞ ) − ( x ) x < x ∗ , ( ˜ f ∞ ) + ( x ) x > x ∗ . (20)The error obtained at any point x ∈ ( x ∗ , s ∗ ) will be, | f ( x ) − ˜ f ∞ ( x ) | = | f − ( x ) − ( ˜ f ∞ ) + ( x ) | ≤ | f − ( x ) − f + ( x ) | + | f + ( x ) − ( ˜ f ∞ ) + ( x ) | . (21)Following similar arguments to the ones used in the first part of the proof,the second term is bounded and is O ( h ). For the first term, we can usesecond order Taylor expansion around s ∗ to write, | f − ( x ) − f + ( x ) | ≤ | [ f (cid:48) ] | ( s ∗ − x ) + sup R \ s ∗ | f (cid:48)(cid:48) | ( s ∗ − x ) ≤ ( s ∗ − x ∗ ) (cid:32) | [ f (cid:48) ] | + sup R \ s ∗ | f (cid:48)(cid:48) | h (cid:33) . (22)As h < h c (see (1)), then | f − ( x ) − f + ( x ) | ≤ | [ f (cid:48) ] | ( s ∗ − x ∗ ) . (23)Now, we can use point 3 of Lemma 2 of [9], that establishes that ( s ∗ − x ∗ ) ≤ C h m sup x ∈ R \{ s ∗} | f ( m ) ( x ) || [ f (cid:48) ] | (with m = 4 in our case), to obtain that the first term14f the right hand side of (21) is also bounded and is O ( h ). This concludesthe proof. · · · f j − f j − f j f j +1 f j +2 f j +3 · · · x j − x j − x ∗ s ∗ x j x j +1 x j +2 x j +3 (cid:100) (cid:100) (cid:100) (cid:100) (cid:100) (cid:100)(cid:105) − (cid:105) + Figure 3: In this figure we can see an example of a corner singularity placedin the interval ( x j , x j +1 ), being h = x j +1 − x j a uniform grid-spacing. Theleft side of the discontinuity placed at s ∗ is labeled as the − side and theright side of s ∗ as the + side. The approximated location of the discontinuityis labeled as x ∗ .Due to the order of accuracy attained by the corrected algorithm close tothe discontinuity and to the regularity and convergence of the linear schemewe can state the following corollary: Corollary 3.3.
The corrected subdivision scheme does not introduce Gibbsphenomenon nor diffusion close to the discontinuities.Proof.
Given the fact that the correction introduced in step 1 of the algo-rithm eliminates the presence of the discontinuity up to the order of accuracyof the linear scheme, we can consider that the linear scheme in step 2 is a-pplied to smooth data after the correction, so no Gibbs effect can appear.Step 3 consists in correcting back the subdivided data, so no diffusion canappear.
Corollary 3.4.
The corrected subdivision scheme reproduces piecewise cubicpolynomials.Proof.
For piecewise cubic polynomials, the location of the discontinuity isexact. Remind that we are using cubic interpolating polynomials for thelocation algorithm. Following Theorem 3.1 we get the proof. Mind thatif the singularity is in the first derivative, the proof given for Theorem 3.1would be equivalent for this case, as now we know that the exact positionof the singularity is x ∗ . 15n what follows, we enumerate the main properties of the corrected al-gorithm: • The corrected subdivision scheme has the same (piecewise) regularityas the linear subdivision used. • The corrected subdivision scheme does not produce any oscillations inthe presence of singularities. • The algorithm has high accuracy without diffusion. • The algorithm is interpolatory. • The corrected subdivision scheme is exact for piecewise polynomialfunctions of the same degree as the one used in the construction of thescheme. • The corrected subdivision scheme is stable due to the stability of thelinear scheme used.Let’s consider now the cell-averages discretization. In this case, we willonly analyze the detection of jumps in the function.
Theorem 3.5. - The cell-averages’ case.
Let f be a piecewise C -smoothfunction on [0 , , with a jump discontinuity at s ∗ , and assume we are giventhe cell-averages ¯ f j = 1 h (cid:90) x j x j − f ( x ) dx, j = 1 , · · · , N. Using the algorithm in Section 3.2, the following results hold:1. The approximate singular point x ∗ satisfies | s ∗ − x ∗ | ≤ Ch .
2. The approximation G ( x ) to the primitive function F ( x ) interpolatesthe data { F j } Nj =1 , and | F ( x ) − G ( x ) | ≤ Ch , x ∈ [0 , .3. The approximation g ( x ) = G (cid:48) ( x ) to f ( x ) satisfies | f ( x ) − g ( x ) | = O ( h ) as h → for x ∈ [0 , , x not in the closed interval between s ∗ and x ∗ . 4. Denoting δ = x ∗ − s ∗ , | f ( x ) − g ( x + δ ) | ≤ Ch ∀ x ∈ [0 , \ { s ∗ } . Proof.
Figure 4 represents a possible scenario for this theorem. Let us proveevery point of the theorem: 16. First we note that the primitive function F is piecewise C . Afterlocating the interval [ x j , x j +1 ] containing the singular point, the al-gorithm finds the two cubic polynomials, ˜ T − and ˜ T + , respectivelyinterpolating the data at four points to the left and to the right of thesingularity interval. As we are using cubic polynomials for the locationof the singularity, the intersection point x ∗ of the two polynomials in[ x j , x j +1 ] satisfies | s ∗ − x ∗ | ≤ Ch , as proved in point 3 of Lemma 3of [9].2. Now we can follow the steps of Theorem 3.1 to deduce that | F ( x ) − G ( x ) | ≤ Ch , x ∈ [0 , F ∈ C ( R ) the 4-pointsubdivision scheme gives O ( h ) approximation order to the derivativeof F . This, together with the observation that both F and G aredifferentiable outside the interval between s ∗ and x ∗ , yield the resultin item 3.4. To prove item 4 let us assume, w.l.o.g., that x ∗ > s ∗ . Both f and g ( · + δ ) have their jump discontinuity at s ∗ . For x < s ∗ we have x + δ < x ∗ . Since | δ | = O ( h ), and g is continuous for x < x ∗ , itfollows that g ( x + δ ) = g ( x ) + O ( h ) . Hence, | f ( x ) − g ( x + δ ) | = | f ( x ) − g ( x ) | + O ( h ) , and by item 3, | f ( x ) − g ( x ) | = O ( h ). Similarly, for x > s ∗ , x + δ > x ∗ and | f ( x + δ ) − f ( x ) | = O ( h ). Hence, | f ( x ) − g ( x + δ ) | = | f ( x + δ ) − g ( x + δ ) | + O ( h ) = O ( h ) . Remark 3.
Inside the interval between s ∗ and x ∗ , as mentioned in Section7 of [9], we can not hope to obtain a better accuracy than O (1) in the infinitynorm. Even though, we can provide a result in the L norm (in fact, thisresult is provided for the L p norm in Section 7 of [9]) having into accountthat, || f ∞ − f || L ≤ || f ∞ − g || L + || g − f || L , where f ∞ is the limit function in the cell-averages and g is the smoothfunction obtained by our algorithm plus the final translation, which has the · · · · · · · · · · · F j F j +1 x j x j +1 s ∗ x ∗ (cid:124) (cid:123)(cid:122) (cid:125) δ s ∗ x ∗ (cid:124) (cid:123)(cid:122) (cid:125) δ x j x j +1 (cid:100) (cid:100) f ( x ) g ( x ) F ( x ) G ( x ) Figure 4: In this figure we represent a jump discontinuity in the interval( x j , x j +1 ), that transforms into a discontinuity in the first derivative for theprimitive. The approximated location of the discontinuity is labeled as x ∗ and the exact location is labeled as s ∗ . The exact function is represented by f ( x ) and the approximation is represented by g ( x ). The exact primitive isrepresented by F ( x ) and the approximation is represented by G ( x ). same cell-averages as the function f and has the discontinuity at x ∗ . Thefirst term of the right hand side of the inequality is controlled by the orderof the linear subdivision scheme (as the correction do not perturb the order)and the second term is controlled using the results in Section 7 of [9]. We would like to check numerically the regularity of the data provided by thecorrected subdivision scheme, that is inherited from the linear subdivisionscheme used. We also want to compare it with the numerical regularityprovided by the ENO-SR method. Thus, an appropriate estimate for theregularity of a subdivision scheme in the point-values or the cell-averagesshould be used. We start by providing the definition of H¨older regularity:
Definition 4.1. An l times continuously differentiable function f : Ω ⊂ R → R is said to have H¨older regularity R H = l + α , if ∃ C < ∞ such that (cid:12)(cid:12)(cid:12)(cid:12) ∂ l f ( y ) ∂x l − ∂ l f ( z ) ∂x l (cid:12)(cid:12)(cid:12)(cid:12) ≤ C | y − z | α , ∀ y, z ∈ Ω Remark 4.
It is clear that if we apply this definition to the primitive F = (cid:82) xa f ( x ) dx of f , the value of l will be reduced by one and the regularity of f will result one order lower than the regularity of F .
18n [23] it is given a definition for the regularity of a subdivision scheme inthe point-value case, considering that ∆ x i are divided differences and ∆ x i are second divided differences: Definition 4.2. A C l subdivision scheme for the point-values is said to haveH¨older regularity l + α l , if ∃ C < ∞ such that lim l →∞ l ! (cid:12)(cid:12)(cid:12) ∆ l f kj +1 − ∆ l f kj (cid:12)(cid:12)(cid:12) ≤ C (2 − k h ) α l , If α l = 1, this definition implies that the subdivision scheme is almost C l +1 . In Section 2 we explained that working with the cell-averages of a func-tion, it is possible to obtain the exact point-values of the primitive. Then,we can just measure the regularity of the function f discretized throughcell-averages using the regularity of the primitive. Definition 4.3.
The subdivision scheme in the point-values for the primitive F is said to have H¨older regularity l + α l if ∃ C < ∞ such that lim l →∞ l ! (cid:12)(cid:12)(cid:12) ∆ l F kj +1 − ∆ l F kj (cid:12)(cid:12)(cid:12) ≤ C (2 − k h ) α l , and this means that the subdivision scheme for the cell-averages of f hasH¨older regularity l − α l , as (cid:12)(cid:12)(cid:12) ∆ l F kj +1 − ∆ l F kj (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) ∆ l − f kj +1 − ∆ l − f kj (cid:12)(cid:12)(cid:12) , Then, following [23], the regularity of a limit function obtained throughpoint-values or cell-averages can be evaluated numerically. Using S and S ,the subdivision schemes for the (undivided) differences of order l = 1 and l = 2 associated to the subdivision scheme S , the following quantities areestimated for l = 1 , j : α l = − log (cid:32) l || ( S k +1 l f ) j +1 − ( S k +1 l f ) j || ∞ || ( S kl f ) j +1 − ( S kl f ) j || ∞ (cid:33) . (24)This expression provides an estimate for α and α such that the limitfunctions belong to C β − and C β − if the data is discretized in thepoint-values. On the other hand, for the case of cell-averages data, we justhave to apply the expression in (24) to the values of the primitive F . In thiscase, the limit function belongs to C β − and C β − . Mind that we haveonly deducted one from the regularity of the limit function for the primitive.19 Numerical results for the case of point-valuessampling
In the experiments presented in this section we have calculated a high accu-racy approximation of the position of the discontinuity using the techniquedescribed in the introduction. We have also computed high order approxi-mations of the jumps in the function and its derivatives using the processdescribed in Section 3.3.Let’s start by the function, f ( x ) = (cid:26) a + (cid:0) x − π (cid:1) (cid:0) x − π − (cid:1) + x + sin(10 x ) , if x < π x + sin(10 x ) , if x ≥ π , (25)with a = 0 and x ∈ [0 , x = π is [ f ( x = π/ f (cid:48) ( x = π/ f (cid:48)(cid:48) ( x = π/ − f ( x ) = (cid:0) x − π (cid:1) (cid:0) x − π − (cid:1) + x + sin(10 x ) , if x < π x + sin(10 x ) , if π ≤ x < π , ( x − π )( x − π −
5) + x + sin(10 x ) , if x ≥ π , (26)that presents two discontinuities. In Figure 7 we present the result obtainedwhen applying the corrected algorithm to the function in (26). We can seethat the algorithm can deal with more than one discontinuity with no prob-lem, just applying steps 1 and 3 of Section 3 as many times as discontinuitiesare detected. 20 FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
Figure 5: Limit function obtained by the linear algorithm (left), the quasi-linear algorithm (center) and the corrected algorithm (right). In order toobtain these graphs we have started from 16 initial data points of the func-tion in (25).
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
Figure 6: Zoom of the limit functions shown in Figure 5 obtained by thelinear algorithm (left), the quasi-linear algorithm (center) and the correctedalgorithm (right).
We have previously mentioned that it is not possible to locate the positionof a discontinuity in the function using a discretization by point-values.Although this is true, it is indeed possible to locate the interval of length h that contains the discontinuity. Our objective is to obtain subdivided datathat keeps the regularity of the linear subdivision scheme applied, that doesnot present diffusion nor Gibbs phenomenon and that have optimal accuracyclose to the discontinuities. Then, as argued in Theorem 3.1, we just assumethat the discontinuity is placed at the middle of the interval. Let’s applythe corrected subdivision algorithm to data obtained from the sampling of21 FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
Figure 7: Limit function obtained by the corrected algorithm after five levelsof subdivision using as initial data the piecewise smooth function in (26) thatpresents two singularities.the function, f ( x ) = (cid:0) x − π (cid:1) (cid:0) x − π − (cid:1) + x + sin(10 x ) + 1 , if x < π ,x + sin(10 x ) , if π ≤ x < π , ( x − π )( x − π −
5) + x + sin(10 x ) + 2 , if x ≥ π . (27)Figure 8 presents the result obtained after five levels of subdivision using100 initial points. As expected, we do not obtain Gibbs phenomenon nordiffusion. It is important to mention that in Figure 8 we have plotted thelimit function and the original data with continuous lines to point out theposition of the jump in the function, but no data or subdivided data isplaced in the middle of the jump.In the following subsections we will check the regularity and the accuracyobtained using the corrected algorithm. Let’s consider the regularity of the limit function obtained when subdividingdata acquired through a point-values discretization of the function (25). InTable 1 we present numerical estimations of the regularity constant of thelinear algorithm, the quasi-linear algorithm and the corrected algorithm. Toobtain this table, we start from 100 initial data points and we subdivide from L = 5 to L = 10 levels of subdivision in order to obtain an approximation ofthe limit function. We measure the numerical regularity for x < π assuringthat the singularity is not contained in the data. From this table we can seethat the numerical estimate of the regularity for the corrected algorithm is22 FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
Figure 8: Limit function obtained after 5 levels of subdivision using thecorrected algorithm for a piecewise continuous function. The initial datahas 100 points.very close to the one obtained by the linear scheme. The quasi-linear schemeclearly is less regular. L β Linear 0.8117 0.8335 0.8507 0.8647 0.8763
Quasi-linear -6.0754 -0.9273 0.03565 0.0723 0.1556
Corrected 0.9967 0.9983 0.9992 0.9996 0.9998 β Linear 0.0167 0.0084 0.0042 0.0021 0.0011
Quasi-linear -9.5160 -1.9331 -0.9654 -0.9281 -0.8446 -0.6277
Corrected 0.5414 0.2706 0.1156 0.0491 0.0227
Table 1: Numerical estimation of the limit functions regularity C β − and C β − for the different schemes presented and the function in (25). In this subsection we present an experiment oriented to check the order ofaccuracy of the schemes presented. In order to do this, we check the errorof interpolation in the infinity norm obtained in the whole domain and thenwe perform a grid refinement analysis. We define the order of accuracy ofthe reconstruction as, order k +1 = log (cid:18) E k E k +1 (cid:19) , k being the error obtained with a grid spacing h k and E k +1 the errorobtained with a grid spacing h k / l ∞ norm, || f k || ∞ = max j =0 , ··· ,N k {| f kj |} . Taking a fine mesh { x j } j =0 , ··· ,N k , we compute E k ∞ = || f L − S Lk ( f k ) || ∞ , (28)where f L is the original data at the resolution L and S Lk ( f k ) is the resultobtained by the subdivision scheme when applied L − k times to the dataat the resolution k . Table 2 presents the results obtained by the three algo-rithms. In order to obtain this table we have applied the subdivision scheme10 times to the original data, thus, L − k = 10. The data at the resolution k has been obtained through the sampling of the function in (25). We cansee how the linear algorithm losses the accuracy due to the presence of thediscontinuity, while the quasi-linear algorithm and the corrected approachkeep high order of accuracy in the whole domain. Corrected Linear Quasi-linear N k E k ∞ order k E k ∞ order k E k ∞ order k
17 2.3041e-02 - 1.1052e-01 - 2.6140e-02 -35 5.3611e-03 2.1036 4.0630e-02 1.4437 5.3611e-03 2.285665 1.6162e-04 5.0518 2.9258e-02 0.4737 1.6164e-04 5.0516129 2.7694e-05 2.5450 5.2048e-03 2.4909 2.7694e-05 2.5452257 1.7574e-06 3.9780 2.5469e-03 1.0311 1.7574e-06 3.9780513 1.0309e-07 4.0916 1.2169e-03 1.0655 1.0309e-07 4.09161025 5.3956e-09 4.2559 8.9145e-04 0.4490 5.3966e-09 4.25572049 2.2313e-10
Table 2: Grid refinement analysis in the l ∞ norm for the function in (25)for the three subdivision schemes analyzed.It is important to remark that for any initial data obtained from a piece-wise polynomial function of degree smaller or equal than three, the correctedsubdivision scheme and the quasi-linear method will obtain an error equalto the machine precision.We can also try to check the accuracy attained by the corrected sub-division scheme scheme when working with piecewise continuous functionsas the ones analyzed in Subsection 5.1. For the next experiment, the dataat the resolution k has been obtained through the sampling of the function24n (25) with a = 10, that is a piecewise continuous function with a jumpdiscontinuity of size a at x = π . The philosophy of this experiment is theone explained in Subsection 5.1: as the exact position of the discontinuityis lost when discretizing a function by point-values, we consider that thediscontinuity is placed at the middle of the suspicious intervals. Then, inorder to obtain the error and the order of accuracy, we compare with theoriginal function but with the discontinuity placed in the middle of the sus-picious interval. The results are presented in Table 3 and we can see that thecorrected algorithm attains the maximum possible accuracy in the infinitynorm. Corrected N k E k ∞ order k
17 3.6320e-02 -35 2.5607e-03 3.826265 1.5596e-04 4.0373129 9.1954e-06 4.0841257 5.6303e-07 4.0296513 3.4794e-08 4.01631025 2.1618e-09 4.00852049 1.3470e-10
Table 3: Grid refinement analysis in the l ∞ norm for the corrected subdivi-sion scheme for the function in (25), that in this case is a piecewise contin-uous function with a jump discontinuity in the function equal to a = 10. Inthis case, as the exact position of the discontinuity is lost when discretizinga function through the point-values, we consider that the discontinuity ofthe function in (25) is not placed at π/ Taking into account what has been explained in previous subsections, wecan try to subdivide two dimensional data. Let’s consider for example thenext bivariate function, f ( x ) = (cid:26) cos( πx ) cos( πy ) , if ( x + ) + ( y − ) < − cos( πx ) sin( πy ) , if ( x + ) + ( y − ) ≥ , (29)that has been represented in Figure 9. In this case we need to know theposition of the discontinuity, for example through a level set function.We can proceed as follows: 25 Detect and locate the possible discontinuities using the rows and/orthe columns of the data. Fit a curve or construct a level set surfacethat allows to locate the position of the discontinuities in the y and x direction with at least O ( h ) accuracy. • Obtain and store the one dimensional correction term for all the rows. • Add the correction term to the rows and subdivide the rows using thelinear algorithm. • Subdivide the columns of the data that have been corrected in the x direction. • Subdivide the columns of the correction term, that is a smooth func-tion in the whole domain. • Use the level set function to set to zero the high resolution correctionterm where we do not need to correct the subdivided data. • Subtract the masked correction term from the subdivided data.Of course, in order to apply this technique successfully, it is necessary thatdiscontinuities are far enough from each other and from the boundaries.Figure 9 top to the left presents the original bivariate data sampled fromthe function in (29). Top to the right we present the resultant subdivideddata obtained following the process described before. Bottom left we canobserve the subdivided correction term, that is clearly smooth. Bottom tothe right we can see the subdivided and masked correction term.
In this section we work with piecewise continuous functions supposing thatthe data is discretized by cell-averages, so that we are able to localize theposition of the discontinuity up to the accuracy needed.In all the experiments presented in this section we will use the followingfunction discretized by cell-averages, f ( x ) = (cid:26)
10 + (cid:0) x − π (cid:1) (cid:0) x − π − (cid:1) + x + sin(10 x ) , if x < π ,x + sin(10 x ) , if x ≥ π , (30)with x ∈ [0 , F of f , the jump inthe function at x = π is [ F ( x = π/
11 0.8 10.6 0.8 y x y x -11 0.8 10.6 0.8 y x y x Figure 9: Top to the left, plot of the function in (29). Top to the right,subdivided data. Bottom to the left subdivided correction term. Bottom tothe right masked subdivided correction termis [ F (cid:48) ( x = π/ − [ f ( x = π/ −
10 and in the second derivative is[ F (cid:48)(cid:48) ( x = π/ f (cid:48) ( x = π/ f j atthe positions x j − + h . Figure 10 shows that the linear algorithm producesoscillations close to the discontinuity. The corrected and the quasi-linearalgorithms do not produce oscillations and attain a very good approximationclose to the discontinuities. Figure 11 shows a zoom around the discon-tinuity. It is important to mention here that the point attained in themiddle of the jump by the quasi-linear method and the corrected algorithmis not due to diffusion introduced by the algorithms, it is due to the kindof discretization used. Mind that, if the function presents a discontinuityin the interval ( x kj − , x kj ), then the discretization in (4) will always returna cell value at some point in the middle of the jump, that can be observed27n the graphs of Figures 10 and 11. This value simply corresponds to themean of the function in the interval that contains the discontinuity, i.e. thecell-average value in (4). FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
Figure 10: Limit function obtained by the linear algorithm (left), the quasi-linear algorithm (center) and the corrected algorithm (right). In order toobtain these graphs, we have started from 20 initial cell-averages of thefunction in (25).
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
FunctionDiscretized dataLimit function
Figure 11: Zoom of the limit functions shown in Figure 10 obtained by thelinear algorithm (left), the quasi-linear algorithm (center) and the correctedalgorithm (right).
In this subsection we analyze the numerical regularity attained by each of thealgorithms for data discretized by cell-averages. In Table 4 we present somenumerical estimations of the regularity constant of the different algorithmsanalyzed. The data has been obtained from the discretization through cell-averages of the function in (30). Following what was done in the point-values28iscretization, in order to obtain this table we start from 100 initial datacells and we subdivide from L = 5 to L = 10 levels of subdivision in orderto obtain an approximation of the limit function. We measure the numer-ical regularity for x < π , assuring that the discontinuity is not containedin the data. From this table we can see that the numerical estimate of theregularity for the corrected algorithm is close to the one obtained by thelinear scheme. The quasi-linear scheme seems to be less regular. L β Linear 0.2240 0.8134 0.8348 0.8518 0.8656
Quasi-linear -8.1221 0.5087 -0.7584 0.1382 0.3238 -9.3917e-03
Corrected 0.9981 0.9991 0.9995 0.9998 0.9999 β Linear -0.7327 2.3471e-03 1.175e-03 5.8784e-04 2.9401e-04
Quasi-linear -12.2278 -0.4911 -1.7577 -0.8619 -0.67621 -1.0095
Corrected 0.2886 0.1270 5.9099e-02 2.9719e-02 9.8887e-03
Table 4: Numerical estimation of the limit functions regularity C β − and C β − for the different schemes presented and the function in (30) dis-cretized by cell-averages. In this section we reproduce the grid refinement analysis that we performedfor the point-values’ case, using the infinity norm in (28), but we also usethe l norm. Taking the mesh { x kj } j =0 , ··· ,N k , we compute E k = h k N k (cid:88) j =1 | ¯ f kj − S ( ¯ f k − ) j | , where, now ¯ f kj = h k (cid:82) x kj +1 x kj − f ( x ) dx is a cell-average value of the function atthe resolution k , (4). In this case S ( ¯ f k − ) j is the result obtained by thesubdivision scheme applied to the point-values of the primitive (5) of ¯ f k − and then returned to the cell-average values at the level k using (6). N k isthe length of ¯ f k and S ( ¯ f k − ) j at the resolution k .In this case we have used the function presented in (30) discretized bycell-averages. Table 5 presents the errors and orders of approximation ob-tained by the three algorithms in the infinity norm. For our approach weshow the error outside the interval [ x ∗ , s ∗ ], in order to check the results sta-blished in Theorem 3.5. For the linear and quasi-linear subdivision schemes29e present the infinity norm in the whole domain. We can see how thelinear algorithm losses the accuracy close to the discontinuity, while thequasi-linear algorithm and the corrected approach keep high order of accu-racy. In particular, the corrected algorithm attains O ( h ) accuracy, that isin accordance with the results of Theorem 3.5. We can also perform thesame grid refinement analysis but using the l norm instead. The resultsare presented in Table 6. We can see in this table that we attain the orderof accuracy expected. Mind that the accuracy has been reduced by one forall the algorithms as we are using a subdivision algorithm with a stencil ofthree cells. The corrected approach attaches the same order of accuracy asthe quasi-linear algorithm. Corrected Linear Quasi-linear N k E k ∞ order k E k ∞ order k E k ∞ order k
64 1.2739e-02 - 5.1079 0.3483 6.0545e-01 2.3782128 2.3556e-03 2.4350 5.5371 -0.1164 2.9149e-01 1.0546256 5.9829e-04 1.9772 5.8782 -0.0862 3.1212e-02 3.22325012 6.5693e-05 3.1870 6.2890 -0.0975 3.2997e-03 3.24171024 7.3102e-06 3.1678 6.5697 -0.0630 3.1356e-04 3.39552048
Table 5: Grid refinement analysis after ten levels of subdivision in the l ∞ norm for the function in (30) discretized by cell-averages and for the threesubdivision schemes presented. Corrected Linear Quasi-linear N k E k order k E k order k E k order k
64 1.2052e-03 - 1.2716e-01 - 2.0014e-03 -128 1.4370e-04 3.0681 8.0777e-02 0.6546 2.7861e-04 2.8447256 1.9401e-05 2.8889 2.2254e-02 1.8599 3.8808e-05 2.8438512 2.0882e-06 3.2158 1.0952e-02 1.0229 4.7043e-06 3.04431024 2.4270e-07 3.1050 5.5043e-03 0.9925 5.8186e-07 3.01522048 2.9298e-08
Table 6: Grid refinement analysis in the l norm for the function in (30)discretized by cell-averages and for the three subdivision schemes presented. In this case we have applied the algorithms analyzed in previous sectionsto two-dimensional data using a tensor product approach, i.e. we directly30rocess one-dimensional data by rows and then by columns.In this section we will work with the bivariate function discretized bycell-averages, f ( x ) = cos( πx ) cos( πy ) , if 0 ≤ x < . , ≤ y < . , − cos( πx ) cos( πy ) + 2 , if 0 . < x ≤ , ≤ y < . , − cos( πx ) cos( πy ) + 2 , if 0 < x ≤ . , . ≤ y < , − cos( πx ) cos( πy ) + 4 , if 0 . < x ≤ , . ≤ y < . (31)Figure 12 shows the result of one step of subdivision using the tensor productapproach for the data presented in Figure 12 top to the left. This meansthat we apply the one dimensional subdivision scheme by rows and then bycolumns. The result presented in Figure 12 top to the right correspondsto the linear algorithm. We can observe that the effect of the discontinuityappears in the subdivided data in the form of diffusion and Gibbs effect. Theresult of the quasi-linear algorithm and the corrected approach is presentedin Figure 12 bottom to the left and to the right respectively. We can seethat both results are very similar. As shown in previous sections, the maindifference is the regularity of the data close to the discontinuity, that ishigher for the corrected approach. In terms of accuracy, both algorithmsperform similar. In this paper we have introduced a corrected implementation of a linearsubdivision scheme in order to adapt the procedure for the approximation ofnon-smooth data. The first advantage of this approach is that the resultantscheme presents the same regularity as the linear scheme. Moreover, thetheoretical analysis of the convergence and regularity is derived directly fromthe properties of the linear scheme. The accuracy of the corrected approachis obtained from the accuracy reached in the location of the discontinuitiesand in the accuracy of the approximation of the jump in the function and itsderivatives, that is done through Taylor’s expansions. Thus, the correctedalgorithm presents the same regularity, at each regularity zone, as the linearscheme [14, 17, 18] plus the same accuracy as the quasi-linear scheme [9].We present the results for the 4-point linear subdivision scheme, but theapproach is indeed applicable to any other subdivision scheme aimed to dealwith the approximation of functions with singularities. By construction, thecorrected subdivision algorithm does not present Gibbs phenomenon anddoes not introduce diffusion. As far as we know, this is the first time that31 y x y x y x y x Figure 12: Result of one step of subdivision using the tensor product ap-proach for the data presented at the top to the left. The figure at the top tothe right corresponds to the linear algorithm. The result of the quasi-linearalgorithm and the corrected approach is presented at the bottom to the leftand to the right respectively.subdivision schemes that own all these properties at the same time appearin the literature. The numerical results confirm our theoretical analysis.
References [1] Amat, S.; Li, Z.; Ruiz, J. On an new algorithm for function approxima-tion with full accuracy in the presence of discontinuities based on theimmersed interface method. J. Sci. Comput. 75 (2018), no. 3, 1500-1534.[2] Amat, S.; Liandrat, J. On the stability of PPH nonlinear multiresolu-tion. Appl. Comput. Harmon. Anal. 18 (2005), no. 2, 198-206.323] Amat, S.; Liandrat, J.; Ruiz, J.; Trillo, J. C. On a power WENO schemewith improved accuracy near discontinuities. SIAM J. Sci. Comput. 39(2017), no. 6, 2472-2507.[4] Amat, S.; Ruiz, J.; Trillo, J. C. On an algorithm to adapt spline ap-proximations to the presence of discontinuities. Numer. Algorithms 80(2019), no. 3, 903-936.[5] Amat, S.; Dadourian, K.; Liandrat, J. On a nonlinear 4-point ternaryand interpolatory multiresolution scheme eliminating the Gibbs phe-nomenon. Int. J. Numer. Anal. Model. 7 (2), (2010) 261–280.[6] Amat, S.; Liandrat, J. On the stability of the PPH nonlinear multireso-lution,
Appl. Comp. Harm. Anal. , 18 (2), (2005) 198–206.[7] Amat, S.; Ruiz, J.; Trillo, J. C.; Y´a˜nez, D. F. Analysis of the Gibbsphenomenon in stationary subdivision schemes. Appl. Math. Lett., 76,(2018) 157–163.[8] Amir, A.; Levin, D. High order approximation to non-smooth multi-variate functions. Comput. Aided Geom. Design 63 (2018), 31-65.[9] Ar`andiga, F.; Cohen, A.; Donat, R.; Dyn, N. Interpolation and ap-proximation of piecewise smooth functions. SIAM J. Numer. Anal. 43(2005), no. 1, 41-57.[10] Ar`andiga, F.; Donat, R.; Mulet, P. Adaptive interpolation of images,Signal Process. (83) (2003), no. 2, 459–464.[11] Baccou, J.; Liandrat, J. Position-dependent Lagrange interpolatingmultirresolutions. International Journal of Wavelets, Multiresolutionand Information Processing 5 (2007) no. 4, 513-539.[12] Beccari, C.; Casciola, G.; Romani, L. An interpolating 4-point C ternary non-stationary subdivision scheme with tension control. Com-put. Aided Geom. Design, 24 (4) (2007) 210–219.[13] Cohen, A.; Dyn, N.; Matei B. Quasilinear subdivision schemes with ap-plications to ENO interpolation. Applied and Computational HarmonicAnalysis 15 (2003) no. 2, 89-116.[14] Daubechies, I.; Lagarias, J. Two scale differences equations: I. Exis-tence and global regularity of solutions, SIAM J. Math. Anal. 22 (1991)1388-1410. 3315] Deslauriers, G.; Dubuc, S. Symmetric iterative interpolation processes.Constr. Approx 5 (1989), 49-68 .[16] Dyn, N.; Levin, D. Subdivision Schemes in Geometric Modelling. Act.Num. 11, (2002) 73–144.[17] Dyn, N.; Levin, D.; Gregory, J. A. A 4-point interpolatory subdivisionscheme for curve design. Comput. Aided Geom. Design, (4) (1987), no.4, 257-268.[18] Dyn, N.; Gregory, J.; Levin, D. Analysis of uniform binary subdivisionschemes for curve design, Constr. Approx. 7 (1991)127-147.[19] Hassan, M.F.; Ivrissimtzis, I.P.; Dodgson, N.A. Sabin, M.A.; An inter-polating 4-point ternary stationary subdivision scheme, Comput. AidedGeom. Design, 19, (2002) 1–18.[20] Hassan, M.F.; Dodgson, N.A. Ternary and three-point univariate sub-division schemes, in: A. Cohen, J.-L. Merrien, L.L. Schumaker (Eds.),Curve and Surface Fitting: Sant-Malo 2002, Nashboro Press, Brent-wood, (2003) 199–208.[21] Jeon, M.; Han, D.; Park, K.; Choi, G. Ternary univariate curvature-preserving subdivision. J. Appl. Math. Comput., 18 (1-2), (2005) 235–246.[22] Kagan, Yael; Levin, D. High order reconstruction from cross-sections.Curves and surfaces, 289-303, Lecture Notes in Comput. Sci., 9213,Springer, Cham, 2015.[23] Kuijt, F. Convexity Preserving Interpolation: Nonlinear Subdivisionand Splines. PhD thesis, University of Twente, (1998).[24] Kwan, P.K.; Byung-Gook, L.; Gang, J.Y. A ternary 4-point approxi-mating subdivision scheme. App. Math. Comp. 190, (2007) 1563–1573.[25] Lipman, Y.; Levin, D. Approximating piecewise-smooth functions. IMAJ. Numer. Anal. 30 (2010), no. 4, 1159-1183.[26] Levin, D. Between moving least-squares and moving least- l1