Algebraic Fourier reconstruction of piecewise smooth functions

Accurate reconstruction of piecewise-smooth functions from a finite number of Fourier coefficients is an important problem in various applications. The inherent inaccuracy, in particular the Gibbs phenomenon, is being intensively investigated during the last decades. Several nonlinear reconstruction methods have been proposed, and it is by now well-established that the"classical"convergence order can be completely restored up to the discontinuities. Still, the maximal accuracy of determining the positions of these discontinuities remains an open question. In this paper we prove that the locations of the jumps (and subsequently the pointwise values of the function) can be reconstructed with at least"half the classical accuracy". In particular, we develop a constructive approximation procedure which, given the first $k$ Fourier coefficients of a piecewise-$C^{2d+1}$ function, recovers the locations of the jumps with accuracy $\sim k^{-(d+2)}$, and the values of the function between the jumps with accuracy $\sim k^{-(d+1)}$ (similar estimates are obtained for the associated jump magnitudes). A key ingredient of the algorithm is to start with the case of a single discontinuity, where a modified version of one of the existing algebraic methods (due to K.Eckhoff) may be applied. It turns out that the additional orders of smoothness produce a highly correlated error terms in the Fourier coefficients, which eventually cancel out in the corresponding algebraic equations. To handle more than one jump, we propose to apply a localization procedure via a convolution in the Fourier domain.


Introduction
Consider the problem of reconstructing a function ℄ R from a finite number of its Fourier coefficients It is well-known that for periodic smooth functions, the truncated Fourier series ´ µ ß Ü converges to very fast, subsequently making Fourier analysis very attractive in a vast number of applications.We have by the classical Lebesgue lemma (see e.g.[30]) that where Å ´ µ is the error of the best uniform approximation to by trigonometric polynomials of degree at most Å.This number, in turn, depends on the smoothness of the function.In particular: (1) If is -times continuously differentiable (including at the endpoints) and ¬ ¬ ´ µ ´Üµ ¬ ¬ Ê, then (see [39, Vol.I, Chapter 3, Theorem 13.6]) Å ´ µ (2) If is analytic, then by classical results of S.Bernstein (see e.g.[30, Chapter IX]) there exist constants and Õ ½ such that 2) Yet many realistic phenomena exhibit discontinuities, in which case the unknown function is only piecewise-smooth.As a result, the trigonometric polynomial F Å ´ µ no longer provides a good approximation to due to the slow convergence of the Fourier series (one of the manifestations of this fact is commonly known as the "Gibbs phenomenon").It has very serious implications, for example when using spectral methods to calculate solutions of PDEs with shocks.Therefore an important question arises: "Can such piecewise-smooth functions be reconstructed from their Fourier measurements, with accuracy which is comparable to the 'classical' one (such as (1.1)This problem has received much attention, especially in the last few decades ( [23,11,9,5,20,8,29,14,15,26,3,22,13,10,4,37] would be only a partial list).It has long been known that the key problem for Fourier series acceleration is the detection of the shock locations.By now it is well-established that classical convergence rates can be restored uniformly up to the discontinuities (see e.g.[22]), but the corresponding question for the jump locations themselves is still open.Notice that any linear approximation procedure with free (a-priori unknown) jump locations will not be able to achieve accuracy higher than ½ Ô Å -see [17].
Several partial results and conjectures in this direction are known, in particular the following.The concentration method of Gelb&Tadmor [19,20] recovers the jumps with first order accuracy, and it can be extended to higher orders.Kvernadze [26] proves that his method can recover jumps of a ¿ function with second order accuracy.In [17,7] we have conjectured that the locations of the jumps of a piecewise function can be recovered with accuracy from its first Fourier coefficients (a similar conjecture is made in [36]).Both Eckhoff [14] and Banerjee&Geer [3] made the same conjectures with respect to their particular reconstruction methods.We would also like to mention a related but different problem: reconstruction of piecewise-smooth functions from point measurements.There, adaptive approximations can achieve asymptotic accuracy for piecewise functions [2,31,27].
With this motivation, our main goal in this paper is to arrive at a better understanding of the "optimal", or the "best possible" accuracy of reconstruction, especially with respect to the locations and the magnitudes of the jumps.As a means to achieve this goal, we develop a reconstruction method which allows for explicit accuracy analysis.Our method is a "hybrid" between a Fourier filtering technique which is first applied to localize the jumps, and the algebraic approach of Eckhoff/Kvernadze which is used in order to resolve each discontinuity one at a time to a high order of accuracy.It is precisely this "localization" which makes the subsequent analysis tractable.
Our accuracy analysis is "asymptotic" in nature, although we provide the explicit constants at every step.These constants in general depend upon various a-priori estimates (such as the minimal distance between the jumps, or the upper bound on the jump magnitudes), which are presumably available.See discussion in Section 2 below, in particular (2.4).
Let us now give a brief summary of the main results.
(1) If a function with a single jump has at least ¾ • ½ continuous derivatives everywhere except the jump, then the jump location can be recovered from the first Å Fourier coefficients with error at most Å ¾ (Theorem 4.13).
In addition, a jump in the Ð-th derivative can be recovered with error at most Å Ð ½ (Theorem 4.21).A key observation in the analysis is that the additional orders of smoothness produce highly correlated error terms in the Fourier coefficients, which eventually cancel out in the corresponding algebraic equations.
(2) The localization step does not "destroy" the above accuracy estimates (Theorem 5.2).Thus, the pointwise values of are recovered with the accuracy Å ½ (Theorem 6.1) up to the jumps.
(3) Numerical simulations are consistent with the theoretical accuracy predictions ( Section 7).
By means of this constructive approximation procedure with provable asymptotic convergence properties, we therefore demonstrate that the algebraic reconstruction methods for piecewise-smooth data can be at least "half accurate" compared to the classical approximation theory for smooth data.
We provide an overview of the reconstruction procedure in Section 2. For expository reasons, the details of the localization step and the analysis of its accuracy are postponed until Section 5.The resolution method of a single jump is presented in Section 3, while Section 4 is devoted to proving its asymptotic convergence order.Finally, the accuracy of the whole reconstruction is analyzed in Section 6.Some common notations used throughout the paper are summarized below.
¯N denotes the natural numbers, R -the real numbers, C -the complex numbers.¯ denotes the class of smooth functions which are continuously differentiable times everywhere.½ is the class of smooth functions having continuous derivatives of all orders.
¯ Ö ´Þµ is the ball of radius Ö centered at Þ, and Ö ´Þµ is the boundary of such a ball.

The algebraic reconstruction method
Let us assume that has Ã ¼ jump discontinuities Ã ½ .Furthermore, we assume that ¾ in every segment ´ ½ µ, and we denote the associated jump magnitudes at by We write the piecewise smooth as the sum ©•¨, where ©´Üµ is smooth and periodic and ¨´Üµ is a piecewise polynomial of degree , uniquely determined by such that it "absorbs" all the discontinuities of and its first derivatives.
This idea is very old and goes back at least to A.N.Krylov ( [25,4]).Eckhoff derives the following explicit representation for ¨´Üµ: where ÎÒ ´Ü µ is understood to be periodically extended to ℄ and Ò´Üµ is the Ò-th Bernoulli polynomial.For completeness, let us dervie the formula for the Fourier coefficients of ¨´Üµ (it can also be found in [14]).Lemma 2.1.Let ¨´Üµ be a piecewise polynomial of degree , with jump discontinuities Ã ½ and the associated jump Proof.One integration by parts yields for and so after • ½-fold repetition we obtain (recall that ¨´ •½µ ¼): A key observation is that if © is sufficiently smooth, then the contribution of ´©µ to ´ µ is negligible for large .Therefore, for some large enough Å one can build from the equations (2.2) an approximate system Here and in the rest of the paper we use the notation def ß .
In fact, this system (up to a change of variables and the number of equations) lies at the heart of the algebraic reconstruction methods of Eckhoff [14], Banerjee&Geer [3] and Kvernadze [26].Banerjee&Geer solve it for all the parameters at once by least squares minimization.Eckhoff and Kvernadze eliminate all the first, resulting in a system of polynomial equations for the , whose coefficients have nonlinear dependence on the initial data.
In contrast, we propose to solve this system separately for each , beacuse this case reduces to a single polynomial equation with respect to .We achieve this "separation" by filtering the original Fourier coefficients such that only the part related to a particular remains.This step requires some a-priori knowledge of the approximate locations of the jumps.Fortunately, such an information can easily be obtained by a variety of methods -see Section 5.
Let us finish this section by presenting the main steps of the reconstruction.We denote the approximately reconstructed parameters with a tilde sign.If not stated otherwise, it is understood that these approximations depend on the index Å.It is important to note that we distinguish between the actual smoothness of the function and the reconstruction order.
Algorithm 2.2.Let be a piecewise-smooth function with jumps at Ã ½ , continuously differentiable ½ times between the jumps.Fix a reconstruction order to be some nonnegative integer

½.
Let there be given the first Å • • ¾ Fourier coefficients of .
( ´ ©µ def ´ µ ´ ¨µ Å Take the final approximation to be 3) The rest of this paper is devoted to providing all the details of the above algorithm and analyzing its accuracy.In particular, we will seek estimates of the form where ¯ £ ££ £££ are some absolute constants depending only on the "size" of the problem; ¯G G ´ ½ Ã µ represents the geometry of the jump points (such as minimal distance between two adjacent jumps); ¯A A ´ ¼ ½ ½ Ã µ represents some a-priori bounds on the jump magnitudes, such as lower and upper bounds; ¯Ê is an absolute bound for the Fourier coefficients of the smooth component ©: ¯ ½ ¾ ¿ and « ¬ are some "simple" functions.
In the course of our investigation we shall be defining more specific bounds, but it will always be assumed that those can be expressed in terms of the above quantities.
Since we are interested in "asymptotic" estimates, we will in general allow the inequalities (2.4) to hold for all Å starting from some index Ã £ which may be large and depend on the parameters of the problem.However, if a particular bound holds for all Ã £ then it will in general hold for ½ ¾ Ã £ as well, with some larger multiplicative constants £ , but which are harder to compute explicitly.

Resolving a single jump
Let def ß .The goal is to recover ¨from the approximate system of equations To find , we eliminate ¼ from the equations.The result is a single polynomial equation having the exact value as one of its solutions.In Eckhoff's paper, this elimination is described in great detail, while here we present only the end result. Let Lemma 3.1.The point satisfies: Proof.The proof is an immediate consequence of Lemma A.4 (see Appendix A).
Since the exact coefficients Ñ (and as a result the polynomials Ô ´Þµ) are unknown, we approximate these with the known quantities Now we are ready to formulate the procedure of recovering the parameters of a single jump.
Algorithm 3.2.Let us be given the first Å • • ¾ Fourier coefficients of the function which has a single jump ¾ ℄.
(1) Solve the polynomial equation Õ Å ´Þµ ¼ and take to be the root which is closest to the unit circle.In Section 4 below, we shall provide the justification for this choice.
(2) The jump magnitudes ¼ are reconstructed as follows.By (3.2), the exact values of satisfy We use the approximations Ö Ñ , and solve the system of linear equations with respect to the unknowns

¨ Ð ©
by any one of the standard methods.

Accuracy analysis: a single jump
Our goal in this section is to analyze Algorithm 3.2 and calculate its accuracy.We shall express all our estimates in terms of the index , keeping in mind that it should be replaced with Å to be consistent with the definitions of the previous sections.
4.1.Accuracy analysis: jump location.We start with the determination of the jump point .Our strategy will be to investigate the polynomial Õ ´Þµ, and determine the bounds on locations of its roots.We can informally summarize the main results as follows: (1) Starting from some , the roots of Õ ´Þµ are "separated" from each other by at least ½ .
(2) If the function is continuously differentiable at least ½ ¾ • ½ times everywhere except at , then one of those roots deviates from the "true" value by at most ¾ .
We regard Õ ´Þµ as a perturbation of Ô ´Þµ.With this point of view, we shall first describe the roots of Ô ´Þµ, and then calculate the "deviations" due to the difference ´Þµ Õ ´Þµ Ô ´Þµ In the subsequent analysis we denote the roots of Ô ´Þµ by Þ ´ µ for ¼ ½ , with the convention that Þ ´ µ ¼ . Also, we denote the roots of Õ ´Þµ by ´ µ .It will be convenient to study Ô ´Þµ in a different coordinate system.For this purpose, consider the following transformation of the punctured Þ-plane: Then the inverse map is given by Now we translate the problem into the Ù-plane.
Therefore it makes sense to study the roots of × ´Ùµ.We denote these roots by ´ µ ¼ . The observation below is an immediate consequence of Lemma 3.1.
Claim 4.3.× ´¼µ ¼ Therefore we will always take ´ µ ¼ ¼ In what follows, we shall break × ´Ùµ into a sum of terms and subsequently apply a perturbation analysis to determine its roots.We begin with some simplifications: Now we make a change in indexing according to the following scheme: With this definition, we can write We will need a technical result.
The polynomial × ´Ùµ must therefore be of the form where in particular We break up the polynomial × ´Ùµ into a "dominant" and a "perturbation" part: × ´Ùµ × ´Ùµ • ´Ùµ where Next we shall see that the dominant component × ´Ùµ determines the locations of the roots of × ´Ùµ up to the first order accuracy, while the other component ´Ùµ is responsible for second-order perturbations of these roots.
We denote the roots of × ´Ùµ by ´ µ ¼ with ´ µ ¼ ¼.It turns out that × ´Ùµ can be completely characterized.

Now we show that
´Ùµ perturbes the zeros of × ´Ùµ by at most ¾ .Since the coefficients ´Ùµ depend linearly on

¼
, we can expect that the bound will depend on the quantity È Ð ¼ Ð .For convenience, let us therefore define There exist constants ¿ Ã½ such that for all Ã½ £ and for all Proof.Our method of proof is based on Rouche's theorem (Theorem B.1).We shall define a sequence ´ µ ¿ £ ¾ (where ¿ is to be determined) and consider disks of radius ´ µ around each one of the roots ´ µ .Our goal is to find ¿ so that ´ µ • ´ µ ß on the boundaries of these disks.
¯In order to bound ¬ from below, we shall use Lemma B.2.We need to bound from below the first derivative at ´ µ , as well as to bound from above the second derivative in the disk
(1) We always have

Now
´ µ is always a root of × ½ ´Ùµ, therefore the value of Ù × ´Ùµ ¬ ¬ ¬ Ù ´ µ is independent of and thus we can write Since all the roots are simple, this is guaranteed to be a strictly positive bound. ( ´ µ ¬ ¬ ¬ ½ and therefore Ù £ ¾ ½ ´½ µ .Using (4.5) and differentiating twice, we get ´ µ • ´ µ ß ¬ ¬ ¬ ´ µ where are some linear functions of ¼ Let ´ µ N R be any function satisfying ¼ ´ µ ½ , and consider Ù ´ µ • ´ µ ß .By Lemma 4.8, for some constant .
We set and let ´ µ ¿ £ ¾ .We need the inequality ´ µ to be satisfied, and this is obviously possible if In this case we have shown that ´ µ ¾ £ ¾ and also which completes the proof.
Remark 4.10.We have in fact shown that for each Ã½ the polynomial × ´Ùµ has precisely • ½ distinct roots.
Now we can go back to the original polynomial Ô ´Þµ and accurately describe the location of its roots ´ µ Ì ½ ´ µ .Being careful to avoid the singularity ´ µ ½ (by choosing large enough ), we now show that the geometry of the roots ´ µ is preserved under Ì ½ .In particular, the numbers Þ ´ µ remain separated from each other (following (4.6)), each of them being close (following Lemma 4.9) to one of the numbers The only thing which is different are the constants.

If in addition
and on the other hand Remaining in the Þ-plane, we now turn to investigate Õ ´Þµ and its roots Ò ´ µ Ó . Recall that we consider Õ ´Þµ to be a "perturbation" of Ô ´Þµ by another polynomial ´Þµ, i.e.

Õ ´Þµ Ô ´Þµ • ´Þµ
The coefficients of ´Þµ depend on the Fourier coefficients of the "smooth part" of our pieceiwise-smooth function .It turns out that in the general setting, the coefficients of ´Þµ are large compared to those of Ô ´Þµ and therefore the perturbations of the roots are large too.If, however, there is enough structure in those coefficients due to additional orders of smoothness, then the perturbation of the roots is small.This is the essense of the key Lemma 4.12 below.
Recall that has in fact ½ continuous derivatives everywhere in ℄Ò , and denote the additional jump magnitudes at by

•½
½ .For every Ð ½, let ¨Ð denote the piecewise polynomial of degree Ð with jump point and jump where © £ is ½-times smooth everywhere in ℄.Thus there exists a constant Ê £ such that ´©£ µ Ê £ ½ ¾ (4.10) Let us also denote The idea of the proof is the same as in Lemma 4.9.Namely, we shall seek the constants ½ ½ and Ã such that if (1) In order to show that ´ µ everywhere on ´ µ ´ µ .Furthermore, in this case This is almost what we want -we would like to have such a bound on the boundary of a neighborhood of ´ µ instead of ´ µ .If ¼ then these values coincide, and so we're done.Otherwise, recall that for Ã½ £ we also have ¿ £ ¾ .So in order to make sure that ´ µ belongs to ´ µ ´ µ , we just require that ´ µ ¿ £ ¾ .
We have thus shown the following: (a) For every function ¼ ´ µ and for every Ã½ £ , the following bound holds for every Ù on the boundary Ã½ £ and ´ µ as above, which additionally satisfies ´ µ ¿ £ ¾ , the above bound holds for every Ù on the boundary of a neighborhood of ´ µ of the same diameter ¾ ´ µ, for every ¼ ½ .
(2) We can now show that similar bounds hold for Ô ´Þµ.Again, only the constants will be different.The map Ì ½ (4.1), being a Möbius transformation, maps ´ µ ´ µ to a circular neighborhood of Þ ´ µ (which is not necessarily ¾ and so Ù £ • ½ ½ ¾ .On the other hand, in this case Ù £ • ½ ¾.
Ã £ and ´ µ as above, which additionally satisfies ´ µ ¿ £ ¾ , the above bound holds for every Þ on the boundary of a neighborhood of Þ ´ µ of the same diameter as above, for every ¼ ½ • ´ µ ß where ´ µ is some function satisfying ¼ ´ µ ½ .Our goal now is to find a uniform upper bound for ¬ ¬ ´Þ µ ¬ ¬ .
(a) By (4.9) we have ¾ and therefore where « ´ µ ¾¼ ¾ ´ µ for some constant ¾¼.Furthermore, using the estimate of Lemma A.3 we have Ã where Ã is an explicit constant (see Lemma A.3). Therefore Combining all the above estimates we therefore have for where ½ is to be determined, and suppose also that We have shown above that there exists a neighborhood ´ µ containing Þ ´ µ of diameter at most ´ µ ½ ¡ À ¡ ¾ such that for every Þ £ ¾ ´ µ we have On the other hand, for every such Þ £ we have by (4.13) ´ ££ • Ê £ µ ¾ ¾ Therefore we must choose ½ and for which the condition is satisfied, together with (4.14).For example: In this case, Õ ´Þµ has a simple zero ´ µ in ´ µ so that •¾ where ½ is to be determined.We again require that ´ µ . We have shown that whenever Ã £ , there exists a neighborhood ´ µ ¼ containing such that for every Þ £ ¾ ´ µ •¾ .On the other hand, by (4.13) for every such Þ £ . This is possible for example when Thus we have completed the proof of Lemma 4.12.
We can finally combine everything and prove the main result of this section.
½ denote the roots of the Laguerre polynomial Then there exist constants ½ ¾ and Ã such that for every Ã À the following statemenets are true: (1) The numbers

Ò
Ý ´ µ Ó lie on the ray Ç , so that ¬ ½, and: ) Algorithm 3.2 provides an approximation for which is accurate up to order ´ •¾µ .
Proof.We have already proved (1) (for Ã¾) and (3) (for Ã¿ £ ) -see Lemma 4.12.(2) follows from Lemma 4.12 and Lemma 4.11 by choosing ¾ ½ • ½½ and Ã À.In order to prove (4), we need to show that no root ´ µ is closer to the unit circle than ´ µ ¼ .From geometric considerations (see Figure 4.1 on page 15), it is sufficient to require that which is true whenever Accuracy analysis: jump magnitudes.Suppose that ½ ¾ • ½ and let Ã À so that our algorithm gives an approximation ´ µ with error at most ½ À ¡ ¾ , in accordance with Theorem 4.13.Our goal is to analyze the accuracy of calculating the approximate jump magnitudes ´ µ Ð , given by the solution of the linear system (3.5).For convenience, we denote The superscripts ´ µ are omitted.The picture on the right shrinks towards the unit circle as ½.
We consider only the case of exactly • ½ equations.Thus we can write this system in the following form: where Î is the ´ • ½µ ¢ ´ • ½µ system matrix

¿ (4.16)
Our goal is to estimate the error ´ µ def ´ µ for ¼ ½ . Let ´ µ def Ö • ´ µ Ñ • Then subtracting (4.16) from (4.15) gives This is the key relation of this section.In order to estimate the magnitude of ´ µ , we shall first write out explicit expansions for the quantities ´ µ , and then investigate how these expansions are transformed when multiplied by the matrix Î ½ .Our analysis will show that the special combination of the structures of both this matrix and the expansion coefficients results in remarkable cancellations.
Let us start by investigating the structure of the matrix Î .
Therefore we can write In light of (4.17), we would now like to examine the action of the matrix Î ½ ¼ on the vectors in the right hand side of (4.20).
Proof.First recall the well-known 1 power series expansion (1) On one hand, the ´Ð • ½µ-st entry in the product on the right-hand side of (4.23) equals to (2) On the other hand, by Proposition B.3 we have for some bounded function Ê ¼ ½ It is now easily seen that the multiplication by Î ½ ¼ "orders up" the vectors in (4.20) by decreasing powers of .Further multiplication by Ë from the left preserves this structure, as is evident from the following calculation.
Lemma 4.20.Let be arbitrary constants.Then there exist constants such that Proof.Let ½ • ½ and consider the -th entry of the product, say Ý : We can now prove the main result of this section.
Theorem 4.21.Assume that ½ ¾ • ½ and Ã À, so that by Theorem 4.13 we have Then there exist constants ¿¿ Ã½¼ such that for every Proof.Combine (4.17

Localizing the discontinuities
As we have seen, both the location and the magnitudes of the jump can be reconstructed with high accuracy.The remaining ingredient in our method is to divide the initial function into regions containing a single jump, and subsequently apply the reconstruction algorithm in each region.
1 It can be proven by induction on , using the identity Our approach is to multiply the initial function by a "bump" which vanishes outside some neighborhood of the -th jump.
This step requires a-priori estimates of the jump positions, which can fortunately be obtained by a variety of methods, for example: (1) The concentration method of Gelb&Tadmor [20]; (2) The method of partial sums due to Banerjee&Geer [3]; (3) Eckhoff's method with order zero (the Prony method).
All the above methods provide accurate estimates of up to first order.For definiteness, we present the description of the last method and a rigorous proof of its convergence in Appendix D. Now, the multiplication is implemented as Fourier domain convolution.Because of the Fourier uncertainty principle, the Fourier series of our bump will have inifinite support and therefore every practically computable convolution will always be an approximation to the exact one.Nevertheless, an error of order at most ½ ¾ in the Fourier coefficients will be "absorbed" in the constant Ê £ (4.10) and therefore we will still have accurate estimates for the reconstruction of each separate jump.This will require us to use bump functions which are ½ .An explicit construction of such a function is provided in Appendix C.
We assume that the following quantities are known a-priori: ¯the lower and upper bounds for the jump magnitudes of order zero: Â½ ¼ Â¾; ¯the minimal distance between any two jumps Â¿ ¼; Our localization algorithm can be summarized as follows.¿ , and this will be possible if Å is not smaller than required by Theorem D.9.
(2) For each : Proof.It is clear that the exact function ¡ has exactly one jump at and jump magnitudes ¼ . In order to prove that Theorems 4.13 and 4.21 can be applied, it is sufficient to show that the error ´ ½•¾µ .By the Fourier convolution theorem the exact Fourier coefficients of are equal to: ´ µ ½ ½ ´ µ ´ µ while our algorithm approximates these by the truncated convolution (5.2).Let us estimate the convolution tail ´ µ ´ µ On one hand, the Fourier coefficients of can be bounded using (5.1): ´ µ ¿ ´Â¾ • Ì µ ½ On the other hand, taking « ½ • ½ we have by Theorem C. 1 ´ µ where ´× Õµ is the Hurwitz zeta function.Therefore Algorithm 3.
In this section we are going to calculate the overall accuracy of approximation.Let us therefore suppose that ½ ¾ • ½, and so using the Fourier coefficients ¾Å ´ µ ¾Å ´ µ we have reconstructed the singular part ¨´Üµ with accuracy specified by Theorem 5.2.Recall from (2.3) that our final approximation is defined by

Å
´ µ ´ ¨µ¡ ß Ü • ẅhere Intuitively, the approximation error function will look as depicted in Figure 6.1 on page 21 -very small almost everywhere except in some shrinking neighborhoods of the jump points.Let Ý ¾ ℄ Ò Ã ½ .
If we take Å large enough so that the error estimate of Theorem 4.13 will be less than the distance to the nearest jump Ý , then Ý will lie in the "flat" region of Figure 6.1 on page 21 and the error ¬ ¬ ´Ýµ ´Ýµ ¬ ¬ will be small.This is precisely the content of our final theorem.
Let us denote the "jump-free" region by , and let it be ½-times continuously differentiable between the jumps.Let Ö ¼.Then for every integer satisfying ¾ •½ ½, there exist explicit functions depending on all the a-priori bounds such that for all Å Algorithm 2.2 reconstructs the locations and the magnitudes of the jumps with accuracy provided by Theorem 5.2, and with the pointwise accuracy ´ ´ µ ´¨µµ ß Ü •Ẅ e write the overall approximation error as Let us examine the two terms on the right-hand side separately.
According to our previous notation, © ¨is a -times continuously differentiable everywhere function.The term Å is easily seen to be the usual Fourier truncation error of ©, since Let ¢ def ¨ ¨denote the "singular error function" (see Figure 6.1 on page 21).Then the second term can be written as: where « and ¬ are provided by Theorem 5.2.For every ¯ Ö, we define Using the formula (2.1) we therefore have The functions Î Ð belong to Ð , and therefore by the well-known estimate (see also (1.1)), there exist constants Ë Ð such that Let us now investigate Ï ´Ýµ.Let Ä denote an upper bound for the magnitudes of the jumps: Clearly the functions Í Ð ¯´Ýµ satisfy: (1) Since Ý ¾ Dr, then Í Ð ¯´Ýµ ¼¯for some absolute constant ¼.This bound can be obtained by just Taylor- expanding the functions Î Ð ´Ý • ¯µ at ¯ ¼.In particular for ¯ « ´Åµ we have (2) In the "no man's land" of length « ´Åµ between and , Í Ð ¯is bounded by ¾ ¡ Ä.Furthermore, as we have just seen, in the flat regions Í Ð «´Åµ is bounded by ¼ «Å ¾ .Therefore the Fourier coefficients of Ï are certainly bounded by and so Combining (6.1), (6.2), (6.3), (6.4), (6.5)and (6.6) completes the proof.

Numerical results
In this section we present results of various numerical simulations whose primary goal is to validate the asymptotic accuracy predictions for large Å.We have used a straightforward implementation and made no attempt to optimize it further.In particular, the Fourier coefficients are assumed to be known with arbitrary precision (this is important for the localization, see below).
7.1.Recovery of a single jump.Given ½ and Å, a piecewise function with one discontinuity is generated according to the formula where the numbers ¾ ℄ Ð ¾ R ½ Ð ¼ and ¾ C Å Å are chosen at random, such that ½ ¾ and .
The Fourier coefficients are calculated with the exact formula These coefficients are then passed to the reconstruction routine for a single jump, of order .This routine implements Algorithm 3.2 in a standard MATLAB environment with double-precision calculations.
The following experiments were carried out: (1) Keeping and ½ fixed, compare the accuracy of recovering the jump location and all the jump magnitudes for different values of Å.The results can be seen in Figure 7.1 on page 23.We also plot the distribution of roots of the corresponding polynomials Õ Å ´Þµ -compare with Figure 4.1 on page 15.
(2) Keeping ½ fixed, compare the accuracy for different reconstruction orders ½ ½.The results are presented in The optimality of ½ ¾ ½, as well as the asymptotic order of convergence, are clearly seen to fit the theoretical predictions.
The instability and eventual breakup of the measured accuracy for large values of Å is due to the finite-precision calculations.Accuracy of reconstructing the jump point, d=4 7.2.Localization.We have restricted ourselves to the following simplified setting: the function has two jumps at ½ ¼ and ¾ ¿, and we localize the jump at the origin by a bump having width ¿ around the initial approximation ½ ½ ¼ .The explicit formulas for the Fourier coefficients of the bump are derived in Appendix C. We have used Mathematica in order to carry out the computations with arbitrary precision.
The results can be seen in Figure 7.4 on page 24.Localization convergence can clearly be seen here, although it starts from very large coefficients.

Discussion
In this paper we have demonstrated that nonlinear Fourier reconstruction of piecewise-smooth functions can achieve accuracy with asymptotic order of at least half the order of smoothness.As indicated by our theoretical results as well as the numerical simulations, a reconstruction method whose order is more than half the order of smoothness becomes less accurate.So it appears that the algebraic approach has certain limitations, and the interesting question is whether these limitations are inherent or superficial.We hope that our results may provide a clue towards obtaining sharp upper bounds.
In addition, it seems that Eckhoff's conjecture is false as stated in [14], namely that the jumps of a piecewise-smooth function can be reconstructed with accuracy ¾ .Using a method of highest possible order doesn't take into account the stiffness of the problem.In fact, it can be shown that the Lipschitz constant of the solution map ´ µ Å• •½ Å Ð of order is proportional to Å .We plan to present these results elsewhere.
Hopefully, our analysis can be related to the algebraic reconstruction schemes of Kvernadze and Banerjee&Geer as well.
Now assume that ´ µ R • R is a given function.Let us perform a change of variable Ý ½ , and define ´Ýµ def ½ Ý ¡ ´ µ.With this notation, we have We subsequently define the "dual" operator as: The operator has an interesting property of "killing" the lowest-order Taylor coefficient at 0.
Proposition A.2. Let À´Ýµ be analytic at Ý ¼, such that À´Ýµ ÑÝ Ñ • Ñ•½Ý Ñ•½ • .Then for Ò ¾ N, the function Ò À ´Ýµ is analytic at ¼ with Taylor expansion The proof is by induction on Ò.The basis Ò ¼ is given.Assuming that Í ´Ýµ Ò ½ À´Ýµ is analytic at 0 and Making the analytic change of coordinates Ý Þ ´Ýµ we conclude that the Taylor expansion of Í ´Þ ´Ýµµ at the origin is This expansion holds in some neighborhood of the origin.
Lemma A.3.Let Ð ¾ N. Then there exist positive constants The proof of the claim is in two steps.First, we shall develop the expression Ð ´ µ into power series in ½ converging for sufficiently large .Then, based on this representation we shall establish the required bound.

´ • ½µ
But then for arbitrary Then we rewrite the given expression as The following statements are true for all × ¼ ½ • ½: ( Now let ´×µ def × ¬ be a polynomial of degree ¬.By (A.1), the above expression can be rewritten as (1) Assume Ø ×.Then « • ½ ¬ • ½, and so by Lemma A.1 we have ¡ «•½ ´×µ ¼.This completes the proof of the first part.
(2) Let Ø × ½.By Lemma A.1 we have ¡ «•½ ´×µ ´« • ½µ and so by (A.3) ´ for all ¾ N and some constant independent of ; (3) For every fixed the following inequality holds for all Þ ¾ ½ ´Þ¼µ Then there exists a constant independent of ´ µ such that for all , the following holds for every Þ ¾ ´ µ ´Þ¼µ: ´ µ Proof.Let us write the truncated Taylor expansion of È around Þ¼ with remainder in Lagrange form: ´ µ ¾ ´ µ.
In general, for approximation of order , we have for Ü Standard majorization of the Taylor series tail by a geometric series.

Appendix C. Explicit construction of a bump
In this appendix we present an explicit construction of the bump function which we used in our numerical simulations.We also derive an explicit bound for the size of its Fourier coefficients, to be used in the proof of localization accuracy.
Let there be given two parameters Ø and with ¾ Ø, together with the point ¾ R. Our goal is to build a function Ø ´Ü ) which satisfies the following conditions: (G1) (G4) the Fourier coefficients of decay as rapidly as possible.
The idea is to take a standard ½ mollifier, scale it and convolve with a box function.
We define two new parameters: the scaling factor × and the width of the box Ö.Note that our construction implies Ö ¾×, because otherwise the result of the convolution will not have a flat segment in the middle.
Let us therefore take the standard ½ mollifier ©´Üµ ´ ½ ´½ Ü ¾ µ for Ü ½ ¼ otherwise and scale it between × and × for some × ¼: Now we take a box function centered at , having width Ö: Finally we convolve the two and get a smooth bump: The new parameters × Ö should be compatible with the original Ø.In particular, we want to have a strip of width Ø in the center, and the extent of the whole bump should not exceed .Therefore we have the following compatibility conditions: The function so constructed clearly satisfies conditions ´ ½µ ´ ¿µ above.Let us now maximize the decay of its Fourier coefficients.By definition: ´ Notice first that ©´Þµ is zero outside the region ½ Þ ½, therefore we can make the change of variables Þ Ø Ü Ø Ø and rewrite the integral as Now we scale back: Þ ×Ý and obtain the explicit formula ´ µ Finally we would like to determine the optimal values for × and Ö so that ´ µ decrease as rapidly as possible with ½.
First note that since © ¾ ½ , then for every « ½ there exists a constant ¼ ´«µ such that ´©µ

¼ ¡ «
The formula (C.2) suggests that we should take × to be as large as possible.Applying the conditions (C.1) we get that the following values maximize ×: We have thus proved the following result: Theorem C.1.Given Ø with ¾ Ø and a point , let Ø ´Ü µ be the bump constructed above.Then it satisfies the conditions ´ ½µ ´ µ such that for every « ½ there exists a constant ½ ½ ´«µ such that for all ¾ N ´ Ø µ ½ ¡ ´¾ Øµ « ½ « (D.1) is a well-known system of equations which is sometimes called the Prony system ( [33]) or Sylvester-Ramanujan system ( [28]).The original method of solution (due to Baron de Prony, [32]) is to exploit the following fact.After this preparation, we can now describe the algorithm for obtaining initial estimates using Prony's method.Recall that our a-priori bounds are given by Â½ Â¾ Â¿ and Ì -see Section 5.
Algorithm D.3.Let us be given the first Å •¾Ã ½ Fourier coefficients ´ µ of a function with Ã unknown discontinuities Ã ½ , which is continuously differentiable between these discontinuities.Denote the magnitudes of the jumps by Ã ½ .We use standard result from numerical linear algebra.
Consider (D.3).The error in the right-hand side is given by (5.1).Therefore we now need to estimate the condition number of the matrix À .Although all the entries are bounded, it may still happen 2 that ´À µ is unbounded.Fortunately, this is not the case.To see this, we are going to factorize À into a component which depends on , and a component which doesn't.Corollary D.6.For all ¾ N ´À µ Â¾ Â½ ´Î µ Remark D.7. ´Î µ is well-studied in e.g.[18].It essentially depends on the minimal distance between the nodes.In particular: Proof.In the context of Lemma D.4, our original system is À q m (D.2) and the perturbed system is À q m (D.3).Note that Ñ Â½ ¡ ¿ for some ¿.From previous considerations we therefore have ´«µ ¡ ¡ « ½

)
Solve the system (3.1) by localization (Algorithm 5.1) and resolution (Algorithm 3.2).This will provide us with approximate values for the parameters the approximate Fourier coefficients of the smooth part:

Lemma 4 . 8 .
The numbersÒ ´ µ Ó satisfy the following properties: 3) The constants and therefore satisfy the assumptions of Lemma B.2.We definedef Ñ Ò ½ ¡ .The conclusion is that there exists a constant such that for every function ´ µ N R satisfying ¼ ´ µ

. ( 3 )
Once we have the lower bound for ¬ ¬ Ô ´Þµ ¬ ¬ on circles of diameter at most ´ µ containing Þ ´ µ , let us now establish an upper bound for ¬ ¬ ´Þµ ¬ ¬ on these circles.Let Þ£ belong to such a circle.On one hand, its distance from Þ ´ µ is at most .On the other hand, by Lemma 4.11 On one hand, we have the bound (4.10).On the other hand, Þ

( a )Theorem 5 . 2 .
Construct the bump centered at with parameters Ø ¾¡ Â¿ ¿ and Â¿, according to Appendix C. Calculate its Fourier coefficients in the range ¿Å ¿Å according to (C.2).(b) Now let ¡ .For each ¼ ½ Å • • ½ calculate Use the above approximate Fourier coefficients ´Åµ ´ µ as the input to Algorithm 3.2 for reconstructing all the parameters of a single jump.Algorithm 5.1 will produce estimates of all these parameters with the accuracy as stated in Theorems 4.13 and 4.21, Ê £ being replaced with some other constant Ê £ Ê £ ´Ê£ Ì Â¾ Â¿µ.

Figure 7 .2 on page 23 .( 3 )
Figure 7.2 on page 23.(3) Keeping the reconstruction order fixed, compare the accuracy for different smoothness values ½.The results are presented in Figure 7.3 on page 23.

Figure C. 1 .
Figure C.1.The construction of the bump

( 1 )( 2 Ö
Calculate the sequenceÖ ¾ ´ß µ ´ µNow we would like to analyze the accuracy of Algorithm D.3.First, we need to estimate the error in solving the system (D.3).
and let Ü¼ • ¡Ü denote the exact solution of this perturbed system.Denote norm ¡ and the induced matrix norm.Then AEÜ ½ ¡ AE ´AE • AE µ (D.4)

8 .
There exist constants ¾ Ã½¾ such that for all ¼ ½

2 9 .
We would like to estimate AEq according to (D.4).If Consider for instance À ½ ½ • ½ ½ ½ ½ We can finally estimate the accuracy of Algorithm D.3.To shorten notation, let ´A G Ì µ def Ì For every ¼ « ½ there exist constants ´«µ Ã½¿ ´« µ such that for every Ã½¿ Algorithm D.3 reconstructs the locations of the jumps with accuracy We have shown that the perturbation É É has coefficients of magnitude ½ .We will use the same reasoning as in Section 4 in order to estimate the quantity ¬ ¬ ¬ ¬ -which is the perturbation of the roots of É ´Þµ.Let Þ ¾.Since the coefficients of É ´Þµ do not depend on , we can obviously find constants and such that (1)É ¼ ´ µ (2) É ¼¼ ´ÞµBy a reasoning similar to Lemma B.2 we conclude that there exist constants ½ ¼ such that for every function ¼ then ´ µ and so (D.5) holds.On the other hand, by Lemma D.8 we have that and therefore É´Þµ has a simple zero inside ´ µ ´ µ. µ « ½ where ¬ ´ µ .Then by Taylor expansion of the logarithm we will have (recall ½) for large enough Ã½¿ ´«µ Substituting ´Ýµ Ý Ð , we conclude that ´Ýµ is a real analytic function of Ý in the disk Ý ½ , and so it can be written as Applying Proposition A.2 to ´Ýµ we conclude that ¼ ¡¡ ¡ Ð• ½ ¼.Therefore the expansion This completes the proof of the second part.Let the polynomial Õ´Þµ ¾ C Þ℄ be a sum Õ´Þµ Ô´Þµ • ´Þµ.Let Þ¼ be a simple zero of Ô´Þµ.If there exists ¾ R • such that Õ´Þµ has a simple zero inside ´Þ¼µ.Lemma B.2.Let there be given a sequence of polynomials È ´Þµ C C and a point Þ¼ ¾ C such that (1) È ´Þ¼µ ¼ for all ¾ N; then