Log In Sign Up

Spectral Methods in the Presence of Discontinuities

Spectral methods provide an elegant and efficient way of numerically solving differential equations of all kinds. For smooth problems, truncation error for spectral methods vanishes exponentially in the infinity norm and L_2-norm. However, for non-smooth problems, convergence is significantly worse---the L_2-norm of the error for a discontinuous problem will converge at a sub-linear rate and the infinity norm will not converge at all. We explore and improve upon a post-processing technique---optimally convergent mollifiers---to recover exponential convergence from a poorly-converging spectral reconstruction of non-smooth data. This is an important first step towards using these techniques for simulations of realistic systems.


page 1

page 2

page 3

page 4


Spectral shock detection for dynamically developing discontinuities

Pseudospectral schemes are a class of numerical methods capable of solvi...

A spectral approach to non-linear weakly singular fractional integro-differential equations

In this work, a class of non-linear weakly singular fractional integro-d...

A New Parareal Algorithm for Problems with Discontinuous Sources

The Parareal algorithm allows to solve evolution problems exploiting par...

A non-local gradient based approach of infinity Laplacian with Γ-convergence

We propose an infinity Laplacian method to address the problem of interp...

I Spectral Methods

i.1 Pseudospectral projection

The principle idea behind spectral methods (and indeed most numerical methods) is to represent a function of interest as a linear combination of a set of suitably chosen basis functions :


where is the partial sum of , the expansion coefficients are given by




is the inner product between and over some domain with weight function .

For smooth functions, this representation suffers an error whose infinity norm decays faster than any power of and which vanishes in the limit of GrandclementSpectral ; GottliebSpectralForHyperbolic .555The details of convergence can vary with spectral representation and with both the smoothness and analyticity of the function being approximated. For example, the spectral representation of an analytic function will converge more rapidly than a representation of a smooth one, although both will converge exponentially. For a formal overview of these issues, see boyd2013chebyshev . Handling nonlinear differential equations in the basis-coefficient (or modal) representation requires regular integration over the whole domain, which is unfavorable. Instead, one can choose a set of collocation points666The optimal choice of collocation points depends on the spectral basis chosen. and use the Lagrange polynomials


as a second set of basis functions. In this case, integration can be approximated via Gauss quadrature:


where are chosen collocation points and are their associated weights. In this approach, interpolates between the collocation points . This family of methods are known as pseudospectral or collocation-spectral methods.

One can easily transform between the collocation representation of a function (which is just the Lagrange basis in disguise) and the spectral coefficients of its interpolating partial sum via the matrix operation


where is the inverse of the Vandermonde matrix. Similarly, one can approximately map the partial sum of to the partial sum of its derivative via the modal differentiation matrix


Moreover, one can approximately map the restrictions to derivative of restricted to via the nodal differentiation matrix


which is nothing more than the modal differentiation matrix appropriately transformed into the collocation representation GrandclementSpectral ; GottliebSpectralForHyperbolic . This results in the following convenient approximations:


i.2 Influence of discontinuities

Figure 2: Illustration of spectral reconstruction and the influence of discontinuities: LABEL: exponentially convergent expansion of a smooth Gaussian. At expansion order the reconstruction already resembles the true solution. LABEL: Gibbs phenomenon on a simple example of a Top Hat function. High frequency oscillations near discontinuities do not vanish with increasing number of modes , as they retain constant amplitude and get ‘pushed closer’ towards the discontinuity.

Unfortunately, the spectral expansion is only so extraordinarily accurate for smooth functions. For Chebyshev polynomials, the -norm of the difference between an arbitrary function and its partial sum is approximately given by


for some positive constant GrandclementSpectral ; GottliebSpectralForHyperbolic . And one can draw similar bounds for other sets of orthogonal basis functions. This means that, in presence of discontinuities, convergence onto the function of interest is non-uniform and pointwise slow—the error decays as away from the discontinuity and the norm exhibits at best convergence thomson2008elementary . Worse, the infinity norm of the error does not converge at all! It is the local behaviour in the vicinity of the function’s discontinuity which spoils the global approximation within the domain—this is the famous Gibbs phenomenon wilbraham1848certain ; gibbs1898fourier ; MichelsonGibbsReply ; boyd2013chebyshev .

Close to the discontinuity, the reconstruction is polluted by high-frequency oscillations of non-decreasing amplitude, irrespective of the number of modes used, as shown in Figure 2. This behaviour prevents the recovery of true solution, making the method unsuitable for discontinuous problems of interest, such as fluid shocks or stellar surfaces.

i.3 Recovering the underlying solution free from Gibbs phenomenon

Before we introduce methods designed for eliminating the Gibbs phenomenon, we mention out an important property of the Chebyshev polynomials:


which implies that for any function, a change of variables implies that


Equations (11) and (12) tell us that Chebyshev expansion is equivalent to a Fourier cosine series with a change of variables. This, in turn, allows us to study Fourier series and apply what we learn to their Chebyshev equivalents. In this subsection we shall refer to the analysis of Fourier expansion, in order to remain consistent with the literature.

i.3.1 Filtering

Vandeven Vandeven1991a suggests that convergence of a Fourier series of a discontinuous function :


can be accelerated by introducing spectral filters , which are smooth functions of compact support , characterised by . Upon application of , one arrives at filtered spectral expansion of function , defined below:


which converges to faster than the Gibbs-affected . This is achieved by decreasing the importance of higher order terms in the spectral expansion, which smooths out Gibbs oscillations, while preserving the low-frequency bulk approximation to

. However, this operation is equivalent to real space convolution with a Fourier transform of the

filter, which acts like a smoothing kernel and results in global smoothing of the function and subsequent smoothing out of the discontinuity.

i.3.2 Mollification

In gottlieb1985recovering , Gottlieb and Tadmor use the idea of spectral filters to introduce mollification, which is a physical space equivalent of filtering carried out in Fourier space:


where is equivalent to the in equation (14) above. They define mollifiers as unit mass, non-negative kernels of compact support. These are then adapted at every point within the domain, such that the essential support—i.e., the majority of the nonzero portion—of does not cross the discontinuities in . To avoid the discontinuities, the method must identify where they are. Therefore, recovery of the Gibbs-free solution requires two basic steps: edge detection, followed by real-space mollification.

In Tadmor2007 Tadmor defines mollifiers by starting from a family of filters characterised by optimal exponential decay in both physical and Fourier space:


He then takes their inverse Fourier transform


to define a global mollifier :


where is the Hermite polynomial of order and the number of frequencies used in the inverse Fourier transform is decided by compactness concerns as described below.

Finally, to ensure appropriate essential support, is dilated as


where is given by


for a free parameter .777 In numerical investigation performed in Tadmor2007 , Tadmor uses the value of , which we implement in our tests as well. The parameter , on the other hand, introduces information about the discontinuities within the domain. The interval is the largest interval of smoothness enclosing and centered on .888A larger region of continuity could be constructed which encloses but is not centered on it. This could, in principle, improve the effectiveness of the mollifier and is worthy of future study. The number of terms in the sum in equation (18) is then given by:


Using this recipe, Tadmor arrives at the exponentially accurate mollifier given by:


and the Fourier-space filtering operation (15) becomes the real-space mollification operation given by


for all points in the domain Tadmor2007 .999Within an interval , the mollifier (22) is symmetric in and . However, when and are sufficiently far apart, this symmetry is broken. The pointwise convergence of spectral reconstruction is sub-linear near a discontinuity, but improves to exponential convergence far from it, in accordance with Equation (24) below:


where the constant is dictated by the specific piecewise analyticity properties of . Close to the discontinuity is very small and hence the error is large. Far away from the discontinuity, however, increases, forcing the error to become small.

Note that there is a different integral and a different mollification operation for each point . In this way, mollification is analogous to solving an inverse-problem via Green’s functions methods. We can see this by explicitly rewriting equation (23) as




For a full formal derivation, we refer the interested Reader to Ref. Tadmor2007

and the references therein. For the estimates of convergence of

see Tadmor2007 p. 373-374.

Figure 3: Illustration of shape evolution of adaptive single-sided mollifiers as they approach a boundary for piecewise mollification. Figure presents 500 mollifiers constructed for a grid of 500 points in the domain. Each point on the -axis represents the central point of a mollifier for a function with discontinuities at and . The -axis is a copy of the -axis and the color is the amplitude of the mollifier. Thus the extent of the color in the -axis is the extent of the mollifier in the physical domain. The mollifier’s support evolves in space, becoming most compact close to the function discontinuity. In order to avoid combining information from two regions of smoothness across the discontinuity, the mollifying kernel is non-zero only on the side of discontinuity of interest.

i.3.3 Preserving discontinuities

In the limit of the dilation operation defined in Equation (19) becomes a Dirac delta function. However, for finite , the mollifier is supported on the whole real line. This means that, even though it becomes narrower near discontinuities, it still admixes information across them, which leads to and overall smoothing of the edge. This is a milder version of the effect discussed in Section I.3.1.

We would like our method to preserve the discontinuous character of , at the same time making use of the excellent smoothing properties of adaptive mollifiers. For this reason we utilize knowledge of jump locations and force the mollifier amplitude to vanish outside the region of smoothness so that:




for a mollifier centered at point surrounded by edge positions at and . The parameter ensures the mollifier is appropriately normalized to have unit mass.

This approach is ad hoc—we do not have a formal justification for truncating the mollifiers in this way. However, we believe that the results presented below justify our choice. A formal analysis of mollifiers which vanish outside the region of smoothness may reveal a more optimal construction. This is a subject for future study.

This method results in recovery of truly discontinuous functions, cured of spurious oscillations. Illustration of the single-sided mollifiers is presented in Figures 1 and 3, which show how are forced to zero immediately beyond the theoretical position of the discontinuity.

Ii Chebyshev Polynomials of First Kind

In the remainder of this work, we focus on Chebyshev pseudospectral methods. We therefore spend a few moments to describe them.

ii.1 Initial decomposition

Figure 4: -norm difference between the spectral reconstruction and true solution as a function of highest order of expansion . The norm has been computed using high precision quadrature integration. LABEL: shows the smooth function case, where the difference decreases exponentially with increasing , leveling off at machine precision. LABEL: exhibits slower convergence due to presence of discontinuities, which retain error at the discontinuity for all and allow for convergence away from it. The computed slope of linear fit to the plot provides the approximate order of convergence.

Chebyshev polynomials of first kind form an orthogonal set on under the weight such that:


Under restriction to the domain and with the application of Equation (29), Chebyshev polynomials can be used to approximately represent systems characterized by non-periodic boundary conditions. (To this purpose the domain must be appropriately shifted and rescaled). As discussed in Section I.1, such a representation (for smooth problems) is subject to an error which decays faster than any power of , as illustrated in Figure 4. The plot shows the norm of the difference between the spectral reconstruction and true solution as a function of the highest order of expansion for a Gaussian function defined below:


where to ensure that the relevant piece of lies well within the domain. The norm decreases exponentially with increasing , leveling off at , reaching machine precision. In the case of a discontinuous function,


the norm of error behaves differently, decreasing non-uniformly at a rate , where , as shown in Figure 4. This trend is consistent with the recognized influence of discontinuities mentioned in Section I.2. The Gibbs oscillations present around discontinuities do not decrease in amplitude, as is increased. They persist for any choice of , conserving their amplitude as they move closer towards the point of jump in the function.

ii.2 Gauss quadrature and projection

Figure 5: Spectral reconstructions obtained both through projection and using Gaussian quadrature. Projection coefficients were calculated using high accuracy quadrature method in Mathematica Mathematica , by computing the solution to equation (2). In the smooth function case results differ for low orders of , however both reconstructions converge onto each other and the true solution exponentially with increasing .
Figure 6: Spectral reconstructions obtained both through projection and using Gaussian quadrature. Projection coefficients were calculated using high accuracy quadrature method in Mathematica Mathematica , by computing the solution to equation (2). In the case of a discontinuous function reconstructions converge onto each other slower than in 5. The high frequency oscillations close to the discontinuity never vanish due to the Gibbs phenomenon.

The Chebyshev polynomials posses another favorable property, which facilitates their application in pseudospectral methods introduced in I.1. The weights and collocation points associated with all types of Chebyshev-Gauss quadratures are completely analytic and can be straightforwardly computed. In this work, we are using the Chebyshev-Gauss-Lobatto quadrature, which includes collocation points at the domain boundaries. Their locations and associated weights are given by:


As described in Section I.1, precise calculation of integrals in order to compute the expansion coefficients is computationally expensive. The quadrature solution provides convenient means of handling the operation using Equation (6), at the cost of introducing truncation error in the approximation of the integral. In the case of smooth functions the pseudospectral and spectral reconstructions converge to the true solution at an exponential rate. In other words, they converge to each other, as illustrated in Figure 5.

Figure 5 presents spectral reconstructions of a smooth Gaussian function using two sets of expansion coefficients. The quadrature ones are computed using the RHS of Equation (5), while the projected ones are computed in Mathematica Mathematica as an inner product between the Gaussian and a given Chebyshev polynomial, using its high accuracy quadrature method. The same procedure is repeated in the case of a discontinuous function in Figure 6. Because of the Gibbs phenomenon, neither the spectral nor the pseudospectral projection converges onto the original function across the entire domain and the projections differ from each other. This effect needs to be recognized as an additional source of error in any method applied to discontinuous functions.

ii.3 Spatial derivatives

Figure 7: Visual comparison between analytic derivative of function defined in Equation (30) and one computed using Equation (7) for a smooth and discontinuous function case. LABEL: shows perfect overlap between the analytic and computed result for highest order of expansion. LABEL: illustrates the influence of the Gibbs phenomenon on the derivative of a discontinuous function. Instead of obtaining two delta functions at the locations of the discontinuities, one observes an oscillating pattern with peaks of highest amplitude located close to the discontinuity. Note that the peaks are not located exactly at the discontinuity. There is some “phase error.”

Spectral reconstructions of derivatives of smooth functions inherit the supreme error convergence properties, visualized in Figure 7. The result computed using Equation (7) completely overlaps with the analytic result, with the -norm of the difference of order , as expected given Figure 4.

With discontinuous functions, however, the Gibbs phenomenon influences the derivative across the entire domain, producing spurious oscillations. This effect is presented in Figure 7, where the computed spectral reconstruction differs significantly from the analytic solution. It is also important to note that the information about the exact position of the discontinuity is subject to uncertainty due to the significant size of spacing between collocation points. This effect is less pronounced closer to the domain boundaries, where collocation points are more closely spaced.

Iii Edge Detection

As we saw in Section II.3, the spectral derivative of a discontinuous function is not well behaved. This is unfortunate, since the (infinitely) large gradient at a discontinuity provides an excellent means of detecting one. Here we describe the technique, first developed by Gelb Gelb1999a ; Gelb2005a ; gelb2006robust , for regularizing these spectral derivatives and using them as a way to localize discontinuities in spectral data.

iii.1 Method

We define the jump function of a piecewise smooth function to be:


In Gelb2001a , Gelb showed that for admissible concentration factors and spectral coefficients , in the basis of polynomial derivatives , the expression below converges pointwise to :


The basic idea here is to filter (à la Vandeven Vandeven1991a ) the spectral expansion of the gradient of so that the spurious oscillations in the gradient disappear.

As discussed in Gelb1999a ; Gelb2001a ; Gelb2008a there exist multiple concentration factors with convergence rate of , which exhibit different behaviour close to the discontinuity. Each

produces its own oscillatory pattern near the jump discontinuities, often trading exponential convergence away from the jump for decreased oscillations close to it. In order to exploit both good convergence and no-oscillation features of individual factors, Gelb and collaborators

Gelb2008a designed a minmod function as


where each corresponds to the jump function approximation computed through Equation (34) with a different concentration factor . Note that the minmod function is a pointwise object, since the chosen can be different at each point . The estimated convergence rate of given in Equation (35) to is given by the lowest formal bound on the individual concentration methods and evaluates to (for a formal proof, see Gelb2008a and references therein).

In our investigation, we follow Gelb2008a and make use of Equation (35) by applying concentration factors first defined in Gelb1999a as





Certainly other concentration factors are valid and the application of additional factors is worth investigating. As in Gelb2008a , we set , , and in our tests. Results of our computations are discussed in Section III.2 below.

Figure 8: Jump functions computed using Equations (36), (37), and (38). LABEL: presents the result of applying , while LABEL: and LABEL: correspond to and respectively. Each jump function is characterised by different oscillation patterns, which allow for the powerful performance of the minmod recipe in (35).

iii.2 Efficacy

Figure 9: Minmod function computed following (35), using concentration factors in (36), (37), (38). LABEL: presents the case of discontinuities located well within the domain, where identification of minmod function extrema allows for straightforward computation of the discontinuity location. Errors are quantified in Figure 10. LABEL: illustrates the inability to easily identify the true discontinuity location because of the minmod function being polluted by spurious extrema. Errors associated with this more difficult case are quantified in Figure 11.

The concentration factors given in Equations (36), (37), and (38) produce jump function approximations akin to those presented in Figure 8. The output is then fed into Equation (35), which results in the functions visualised in Figure 9. Location of discontinuities is then found by searching for local extrema in the minmod result.

In simple cases, when the discontinuities are located well within the domain and have comparable magnitudes, the method proves successful, returning discontinuity position with an error of for highest order of expansion. To illustrate this claim, we present the error in discontinuity location and its behaviour with increasing in Figure 10, where the jumps occur at and . Oscillations in the error behaviour are due to changing spatial distribution of collocation points with increasing number of modes. The envelope, however, seems to obey the error convergence predicted by theoretical arguments in Gelb2008a . As shown in Figure 9, the minmod function cleanly predicts two peaks, identifying the edges without issue.

Figure 10: Difference between the true discontinuity location and the one computed using edge detection recipe provided by Gelb2008a for a discontinuity located well within the domain, as a function of order of expansion . There were two discontinuities of the same jump magnitude, present at and on the domain . The minmod function is sampled on a real space grid of 500 points. The edge detection algorithm performs well, with the error convergence upper bound following relation derived in Gelb2008a .

Discontinuity location accuracy is not only a function of the order of expansion , but also the spatial position of the jump within the domain. In the center of the domain, where the collocation point density is smallest, the accuracy is naturally lowest. As the discontinuity moves closer to the domain boundary, the magnitude of difference between detected and true jump location decreases due to higher density of collocation points—compare the absolute values of the error in Figure 11 with those in Figure 10. However, multiple peaks begin to appear in the minmod function, due to poor behaviour of jump function approximations close to the domain boundary. These spurious peaks can make edge detection unreliable near domain boundaries. This issue is well illustrated in Figure 9, where clear identification of the correct maximum is not possible. Knowledge of spectral expansion coefficients at this point is not enough to identify a single discontinuity in this regime.

In the attempt to mitigate this effect and allow for correct discontinuity location up to the domain boundaries, we solved the problem heuristically. Our approach utilized spectral reconstruction of the derivative of the function of interest, which contains information about the expected sign of the jump. Despite the presence of Gibbs phenomenon, we were able to discard spurious minmod peaks by comparing their sign against neighboring extrema in the derivative. In the case of fluid shocks, one might be able to use shock indicators for further guidance. This approach is not yet robust and requires more experimentation in more realistic settings. Another way to mitigate this difficulty might be via nonlinear enhancement, as described in

Gelb2001a .

These spurious peaks provide further difficulties in the presence of multiple discontinuities of varying magnitudes. It is not always clear how to distinguish a smaller, physical peak in the minmod function from a spurious one. Moreover, a small discontinuity may be washed out entirely if it is accompanied by a larger one.

Figure 11: Difference between the true discontinuity location and the one computed using edge detection recipe provided by Gelb2008a for a discontinuity located close to the domain boundary, as a function of order of expansion . There were two discontinuities of the same jump magnitude, present at and on the domain . The minmod function was sampled on a real space grid of 500 points. The edge detection algorithm performs well for high enough expansion order with the absolute error lower than in the case of discontinuity location well within the domain. However, it finds spurious peaks and fails to find the physical discontinuity for .

Iv Mollification

Figure 12: Effect of piecewise mollification applied to a top-hat function projected onto a basis of Chebyshev polynomials. Piecewise mollification (thick solid) smooths out the Gibbs oscillations in the reconstruction (dashed) and preserves the discontinuous properties of the original function (fine solid), irrespective of a shallow slope of the reconstruction. LABEL: Two discontinuities located well within the domain. LABEL: A single discontinuity well within the domain.
Figure 13: Pointwise difference between the analytic and computed result for mollified (left) and unmollified (right) spectral reconstruction sampled on an equally spaced grid of 500 points in the domain . The panels show the case of a single discontinuity within the domain, as presented in Figure 12. The error in mollified reconstruction converges exponentially away from the discontinuity, as shown in Figure 15, while the unmollified reconstruction suffers poor convergence. Performance of mollification with modes involved is comparable to that of the unmollified spectral reconstruction.
Figure 14: Pointwise difference between the analytic result and mollified spectral reconstruction sampled on an equally spaced grid of 500 points in the domain . Here we have two discontinuities within the domain, as seen in Figure 12. The error exhibits decay between the discontinuities and decay at domain boundaries, which is further demonstrated in Figures 15 and 16.
Figure 15: Difference between the analytic function and mollified spectral reconstruction at as a function of expansion order . The point was chosen carefully, to avoid probing the spectral reconstruction at a collocation point. Cases under consideration are the same as in Figures 14 and 13—two discontinuities within the domain (left) and a single discontinuity (right). Piecewise mollification seems to prove its exponential accuracy away from a single discontinuity (right), however suffers worse performance between the two of them (left). The spectral reconstruction behaves poorly in both cases, not exhibiting convergence at all.
Figure 16: Difference between the analytic function and mollified spectral reconstruction at x=-1.00 domain boundary as a function of expansion order . The case under consideration is the same as in Figure 14—two discontinuities within the domain. Piecewise mollification seems to provide approximately 3rd order convergence at the domain boundary. The spectral reconstruction case is not included, as the domain boundary is one of the collocation points for Chebyshev-Gauss-Lobatto quadrature. By construction, this point is required to agree with the function exactly.

Mollification, as discussed in Section I.3.2 is an operation carried out on the spectral reconstruction in real space. In our tests, we approximate continuum mollification by discrete mollification on a very fine, evenly-spaced grid of 500 points ( values) within the domain . For each , we compute an appropriate adaptive mollifier, based on the point position within the domain and the location of the closest detected discontinuity. We then multiply the Gibbs-polluted reconstruction by the mollifier at all

and perform discrete integration over the whole domain (e.g., via Chebyshev integration or via a quadrature rule). This means we compute 500 integrals, which is potentially very expensive. For realistic problems, a more efficient integration scheme is probably desirable.

iv.1 The Handling of Boundaries

Non-periodic boundary conditions require careful handling of domain boundaries. Lack of information outside of the domain influences the integration result, with the effect most pronounced for grid points closest to the edges of the domain at and . Every grid point outside of the domain effectively gives zero contribution to the discrete integral, creating an effect akin to a smoothed discontinuity, discussed in Section I.3.1.

In an attempt to improve convergence at the boundaries, we considered different approaches such as the introduction of narrow buffer zones outside of the domain and treating the boundaries as discontinuities to renormalize the mollifiers. We first extended the information beyond domain boundaries by adding a mirror reflection of the values within from the boundary 101010 is defined by (20) with and for the respective boundaries. to make sure we cover the width of essential support of the mollifier and minimize the influence of missing values on the discrete integral. We have chosen to mirror the values in order not to introduce bias and only make use of the reconstructed function. This approach yielded promising results, producing errors of order for highest order in Chebyshev expansion.

We then tried a different approach, using ideas derived from mollifier modification across the discontinuity as described in Section I.3.3. We required the mollifiers to reject information beyond the domain by forcing their amplitude to zero across the boundary and renormalizing them to have unit mass. This solution resulted in a better error behavior, with values of order for highest order in Chebyshev expansion. It also had the advantage of being entirely contained within the domain, therefore we decided to implement it in our method.

iv.2 Results and error convergence

Piecewise mollification with the single-sided mollifiers illustrated in Figures 1 and 3 has proved effective in reducing Gibbs oscillations, while still preserving the discontinuous character of the functions under investigation. Figure 12 shows the mollified result plotted against the unprocessed spectral reconstruction and the original top hat function computed for modes in the Chebyshev polynomial expansion. Figure 12 presents the top hat positioned well within the domain, while Figure 12 illustrates a possible result of advection of the top hat waveform such that majority of it has already managed to leave the domain. In both cases, the mollified reconstruction removes unwanted oscillations, and retains a (slightly smoothed out) discontinuity.

In order to quantify the method’s performance we compute a pointwise error across the domain by taking the difference between mollified result and analytic function value at our 500 discrete . We show this error for the single discontinuity case in Figure 13 and the two discontinuity case in Figure 14. The pointwise error depends strongly on the position and number of discontinuities. Nevertheless performance is quite satisfactory in all cases we investigated. In the single discontinuity case, we also compare the pointwise errors to the unmollified case. The improvement via mollification is dramatic.

These relationships are further investigated in Figures 15 and 16, where we plot the pointwise difference between the mollified reconstruction and the analytic function for and respectively. We find that for a point away from a single discontinuity the rate of error convergence exhibits exponential behaviour, while it seems to be of order between the two discontinuities. At the domain boundaries, regardless of the solution character within it, the error seems to converge at approximately 3rd order.

Our results are in qualitative agreement with Tadmor2007 , where the author demonstrates exponential error convergence away from discontinuities, nearly error at the discontinuity itself and similar magnitude of error present at the domain boundaries. However, it is important to recognize that our numerical tests used different trial functions, spectral projection and discrete integration methods, which do not allow for a direct comparison between the two results.

iv.3 Stability of the method

In some cases, the edge detection algorithm may fail to detect an edge or locate it with significant error. (See for example Figure 9.) In both cases the mollified solution still retained its good behaviour away from the discontinuity, with the increase in error influencing only the local environment of the underlying real jump. In case of no edge detection close to the boundary, the function would simply lose its discontinuous character there and be smoothed out, as if we applied filtering to its spectral coefficients. When the discontinuity was misplaced w.r.t to the original jump, it would still remain discontinuous and without Gibbs oscillations, however the error would be locally extended in . Further discussion of the robust character of the method is presented in Section V.

iv.4 Computational Cost

Throughout our investigation, we performed the discrete integration at every point using simple trapezoidal rule. To ensure that this second-order error does not contribute to our analysis, we chose a very fine grid for integration. Since one needs to perform one integral per sampling point, this procedure introduces significant computational cost—roughly the number of sampling points squared. Reducing this high cost is an impediment to the successful application of mollifiers to realistic problems and it requires careful attention in the future.

V Hyperbolic PDEs

As a final test and as a proof-of-concept, we solve the one-dimensional linear advection Equation


on the domain with periodic boundary conditions. We use a discontinuous tophat function like Equation (31) given as initial data. We solve Equation (39) using the method of lines and the pseudospectral methods described in Section I.1. We then post-process the solution via the discontinuous mollifiers defined in Section I.3.3. We detect edges using Gelb’s minmod method as described in Section III. We recommend this same basic approach for future applications of this method. By performing the pseudospectral simulation, one can leverage the efficiency of spectral methods. Then, when one wants to analyze the output, one can choose specific times and snapshots to post-process and recover spectral accuracy via mollification.

v.1 Results for Discontinuous Mollification

We present snapshots of the mollified solution to Equation (39) in Figure 17. When the discontinuities are far from domain boundaries, the solution behaves as described in Section IV.2. There is a slight smearing of the discontinuity (providing error), but otherwise the mollified solution is extremely good. When a discontinuity is very near domain boundaries, however, this smoothing increases and can spoil the mollified solution near the discontinuity.111111We emphasize that a failure of the mollifier does not prevent the successful evolution of the system. All mollification is performed in post-processing. Fortunately, this smearing is local, and the global solution is still well-behaved. We believe this robustness of the mollifier is critical for its successful application in the future.

Figure 17: Snapshots of the mollified reconstruction of an advected square wave on a periodic domain. The reconstruction suffers from significant errors when the wave reaches the domain boundary, when edge detection method fails. This, however does not entirely inhibit the recovery of the underlying solution. Read from left to right, top to bottom: LABEL: the square wave begins well within the domain and the mollified reconstruction provides 4th order convergence presented in Figure 15; LABEL: as the wave approaches domain boundary the reconstruction is not affected; LABEL: the mollified reconstruction suffers significant error at the domain boundary, however does not prevent recovery of information between the detected discontinuity and the boundary; LABEL: the reconstruction remains stable as the wave crossed the periodic boundary; LABEL: the reconstruction is only affected in the direct vicinity of x=1.0 boundary; LABEL: the reconstruction behaves just like in LABEL:

v.2 Comparison to Gegenbauer Polynomials

The Gegenbauer reconstruction, first proposed by Gottlieb and Shu Gottlieb1992a has received much attention as a post-processing technique that enables spectral methods to resolve shocks and discontinuities. The Gegenbauer polynomials are those Jacobi polynomials which are orthogonal under the norm


In this approach, smooth regions of a discontinuous spectral solution are reprojected onto a finite-sized basis of a subset of the Gegenbauer polynomials with fixed .  is chosen such that the new basis is sufficiently “different” from the original spectral basis. In this way, Gibbs oscillations can be removed up to the discontinuity and a spectrally accurate, discontinuous solution can be constructed. For more details on the Gegenbauer reconstruction, see Gottlieb2011a and references therein. For some examples of “realistic” one-dimensional applications of the Gegenbauer approach, see Gelb2000Enhanced ; Gelb2001EnhancedSph ; MeisterFilter .121212Gelb2001EnhancedSph studies a two-dimensional system, but the Gegenbauer reconstruction is only applied to a one-dimensional test case.

The Gegenbauer reconstruction thus fills the same role as our discontinuous mollifiers. It is a post-processing technique that removes the Gibbs oscillations from a spectral solution of a PDE system and it requires some method of finding the discontinuities, such as the one described in Section III. It is therefore worth comparing the Gegenbauer reconstruction to our mollifiers.

Figure 18 shows the Gegenbauer reconstruction for the advection equation (39) using the techniques described in Section I.1 and the method of lines. We use Equation (31) as initial data. Our Gegenbauer code is open source and the interested Reader can find it in MillerGegenbauerCode . The figure shows that the reconstruction behaves well when the discontinuities are in the center of the domain. However, when the discontinuities are near domain boundaries, the global solution breaks down. In this way, the Gegenbauer reconstruction is not robust. As we discuss in Section V, mollified reconstructions are not vulnerable to this instability.131313We note that there are techniques designed to improve the robustness of the Gegenbauer reconstruction, which we did not explore. For example, see Gelb2005a , gelb2006robust and references therein.

In realistic settings, shocks and discontinuities may well be present near domain boundaries. Moreover, in multiple dimensions, it may be much more difficult to localize a discontinuity. We believe the robustness of the mollifiers, compared to the Gegenbauer reconstruction may prove critical in handling these difficulties.

Figure 18: Relevant snapshots of a Gegenbauer reconstruction of an advected square wave on a periodic domain. The reconstruction fails when the wave reaches the periodic domain boundary. Read from left to right, top to bottom: LABEL: the square wave begins near the center of the domain and the Gegenbauer reconstruction provides root-exponential convergence up to the point of discontinuity; LABEL: as the square wave approaches the edge, the reconstruction closest to the edge fails; LABEL: the Gegenbauer reconstruction fails entirely for the portion of the solution near the domain boundary; LABEL: the reconstruction stabilizes while the square wave crosses the domain boundary; LABEL: the reconstruction completely breaks down; LABEL: the reconstruction recovers and converges exponentially far from the domain boundaries.

Vi Concluding Thoughts

In the course of this work, we have investigated the efficacy of several techniques that detect discontinuities in spectral data and that help reduce the error introduced by the associated Gibbs phenomenon. We have found that, although these techniques are very promising, they require careful implementation for practical applications. The spectral edge detection developed by Gelb is a good example of this. Theoretical accuracy is attained in ideal situations, such as for a discontinuity in the center of the domain, but spurious oscillations spoil the solution when the discontinuity is near a domain boundary. (See, e.g., figures 10 and 11.)

Moreover, the best technique in theory may be inferior in practice. In ideal cases, the Gegenbauer reconstruction provides uniform root-exponential convergence (in the number of modes of the original system) to a true solution, up to and including discontinuities Gottlieb1992a . However, Boyd demonstrated that as the number of Gegenbauer modes increases, convergence can easily be destroyed BoydTroubleWithGegenbauer . In our own experiments, we found the Gegenbauer reconstruction to perform very poorly when the discontinuity was poorly localized (as is the case with Gelb’s edge detection near domain boundaries), while Tadmor’s mollifiers performed significantly better. (See, e.g., figures 18 and 17.) We believe our experience aligns well with Boyd’s argument.

As a result of our experiments, we have improved upon both the edge detection developed by Gelb and collaborators Gelb1999a ; Gelb2001a ; Gelb2008a and the mollifiers developed by Tadmor and collaborators gottlieb1985recovering , Tadmor2002 , Tanner2006 . In particular, we have developed several techniques for handling erroneous behaviour at the boundaries of the domain and, more importantly, we have introduced mollifiers that vanish outside the region of smoothness. These one-sided mollifiers allow for the recovery of truly discontinuous solutions in a way that is more resistant to perturbation than the Gegenbauer reconstruction. We believe our improved versions of these techniques are a promising path towards using high-order spectral methods in production simulations of non-smooth problems. However, there are many obstacles that must be overcome. The study presented here is preliminary—it is only in 1D and covers relatively simple, linear systems.

In higher dimensions, the increased dimensionality of the discontinuity itself poses a problem. In one dimension, a discontinuity is localized to a point. However, in two and three dimensions it is a line and a surface respectively. Worse, the extra dimensionality allows discontinuities to have complex geometric and topological structure.141414And of course, in many dimensions, as in one, there may be multiple disjoint discontinuities. The obvious extension of our one-dimensional results to this setting is via a Cartesian product. However, it is not clear how effective this approach will be.

Nonlinear systems provide another challenge. In the linear case, the error introduced by the Gibbs phenomenon simply advects across the grid. However, in nonlinear systems, aliasing error will translate this error into “physical” source terms and can drive an instability. Tadmor’s spectral viscosity method provides stability and convergence criteria for spectral solutions of these nonlinear systems Tadmor1990 . However, it is not clear how effective this technique will be when multiple kinds of physics and solution methods, for example gravity and radiation transport, are added to the calculation.

Another issue with nonlinear systems is that they have regimes of validity. For example, the Gibbs phenomenon may drive fluid density or pressure to a non-positive value. One may need to play tricks (such as an artificial atmosphere) to preserve stability. And these tricks may damage global convergence.

On the other hand, nonlinear systems also provide physical measures of success and physical tools which can be used to detect shocks and discontinuities, perhaps augmenting the edge detection discussed here. One could, for example, look at generalized Riemann conditions or entropy production for contact discontinuities and shock discontinuities respectively.

Finally, mollification efficiency is a challenge. Mollifying a spectral reconstruction requires one integral per sample point. If the quadrature points of numerical integration are the sample points, this translates into a computational cost quadratic in the number of sample points. This cost is somewhat mitigated by the fact that mollification is performed as a post-processing step and does not need to be applied at every time step. Nevertheless, an efficient and adaptive integration scheme is required.151515We note that the fact that the mollifier depends on space, and in fact changes dramatically when one moves across a discontinuity, precludes its evaluation in the spectral domain. Therefore, we can not take advantage of more efficient fft methods to improve efficiency.

We do not believe these difficulties are insurmountable. However, they require careful study. Indeed, we believe they are promising avenues of future research. The extreme efficiency of spectral methods mean that, if these difficulties are resolved, calculations that are currently performed on a supercomputer could be performed on a desktop. These potential gains make this extra care worthwhile.

Vii Acknowledgments

We thank Roland Haas and Oleg Korobkin for their insightful comments on the draft. JMM thanks David Radice for initially guiding him to the appropriate literature.

The authors acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). The research was also supported by the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. JMM acknowledges the support of the USA’s National Nuclear Security Administration and the US Department of Energy’s Scientific Discovery Through Advanced Computing program.

We are grateful to the countless developers contributing to open source projects on which we relied in this work, including Python rossumPythonWhitePaper , numpy and scipy numpy ; scipyLib , and Matplotlib hunterMatplotlib .


  • (1) J. Boyd. Chebyshev and Fourier Spectral Methods: Second Revised Edition. Dover Books on Mathematics. Dover Publications, 2013.
  • (2) J. P. Boyd. Trouble with gegenbauer reconstruction for defeating gibbs’ phenomenon: Runge phenomenon in the diagonal limit of gegenbauer polynomial approximations. Journal of Computational Physics, 204(1):253 – 264, 2005.
  • (3) L. Boyle, M. Kesden, and S. Nissanke. Binary black-hole merger: Symmetry and the spin expansion. Phys. Rev. Lett., 100:151101, Apr 2008.
  • (4) A. Gelb and D. Cates. Detection of edges in spectral data III-refinement of the concentration method. Journal of Scientific Computing, 36(1):1–43, 2008.
  • (5) A. Gelb and J. P. Gleeson. Spectral viscosity for shallow water equations in spherical geometry. Monthly Weather Review, 129(9):2346–2360, 2001,¡2346:SVFSWE¿2.0.CO;2.
  • (6) A. Gelb and Z. Jackiewicz. Determining Analyticity for Parameter Optimization of the Gegenbauer Reconstruction Method. SIAM Journal on Scientific Computing, 27(3):1014–1031, 2005.
  • (7) A. Gelb and E. Tadmor. Detection of Edges in Spectral Data. Applied and Computational Harmonic Analysis, 7:101–135, 1999, 0112016.
  • (8) A. Gelb and E. Tadmor. Enhanced spectral viscosity approximations for conservation laws. Applied Numerical Mathematics, 33(1):3 – 21, 2000.
  • (9) A. Gelb and E. Tadmor. Detection of Edges in Spectral Data II . Nonlinear Enhancement. Society for Industrial and Applied Mathematics Journal of Numerical Analysis, 38(4):1389–1408, 2001.
  • (10) A. Gelb and J. Tanner. Robust reprojection methods for the resolution of the gibbs phenomenon. Applied and Computational Harmonic Analysis, 20(1):3–25, 2006.
  • (11) J. W. Gibbs. Fourier’s series. Nature, 59(1522):200, 1898.
  • (12) D. Gottlieb and J. Hesthaven. Spectral methods for hyperbolic problems. Journal of Computational and Applied Mathematics, 128(1):83 – 131, 2001. Numerical Analysis 2000. Vol. VII: Partial Differential Equations.
  • (13) D. Gottlieb, C.-W. Shu, A. Solomonoff, and H. Vandeven. On the Gibbs phenomenon 1: Recovering exponential accuracy from the Fourier partial sum of a non-periodic analytic function. 43:81–98, 1992.
  • (14) D. Gottlieb and E. Tadmor. Recovering pointwise values of discontinuous data within spectral accuracy. In E. M. Murman and S. S. Abarbanel, editors, Progress and Supercomputing in Computational Fluid Dynamics, Progress in Scientific Computing. Springer, Birkhäuser Boston, 1985.
  • (15) S. Gottlieb, J. H. Jung, and S. Kim. A review of David Gottlieb’s work on the resolution of the Gibbs phenomenon. Communications in Computational Physics, 9(3):497–519, 2011.
  • (16) P. Grandclément and J. Novak. Spectral methods for numerical relativity. Living Reviews in Relativity, 12(1):1, Jan 2009.
  • (17) E. Hewitt and R. E. Hewitt. The gibbs-wilbraham phenomenon: An episode in fourier analysis. Archive for History of Exact Sciences, 21(2):129–160, Jun 1979.
  • (18) J. D. Hunter. Matplotlib: A 2D graphics environment. Computing In Science & Engineering, 9(3):90–95, 2007.
  • (19) W. R. Inc. Mathematica 9.0, 2012.
  • (20) E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001–.
  • (21) L. Kidder, H. Pfeiffer, M. Scheel, et al. Spectral Einstein Code, 2000.
  • (22) A. Meister, S. Ortleb, and T. Sonar. Application of spectral filtering to discontinuous galerkin methods on triangulations. Numerical Methods for Partial Differential Equations, 28(6):1840–1868, 2012.
  • (23) A. A. Michelson. Fourier’s Series. Nature (London), 59:200, Dec. 1898.
  • (24) J. M. Miller. Spectral shock capturing in python.
  • (25) G. Rossum. Python reference manual. Technical report, Amsterdam, The Netherlands, The Netherlands, 1995.
  • (26) M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews, and H. P. Pfeiffer. High-accuracy waveforms for binary black hole inspiral, merger, and ringdown. Phys. Rev. D, 79:024003, Jan 2009.
  • (27) E. Tadmor. Shock capturing by the spectral viscosity method. Computer Methods in Applied Mechanics and Engineering, 80(1-3):197–208, 1990.
  • (28) E. Tadmor. Filters, mollifiers and the computation of the gibbs phenomenon. Acta Numerica, pages 305–379, 2007.
  • (29) E. Tadmor and J. Tanner. Adaptive mollifiers for high resolution recovery of piecewise smooth data from its spectral information. Foundations of Computational Mathematics, 2(2):155–189, Jan 2002.
  • (30) J. Tanner. Optimal filter and mollifier for piecewise smooth spectral data. 75:767–790, 04 2006.
  • (31) B. Thomson, J. Bruckner, and A. Bruckner. Elementary Real Analysis, Second Edition., 2008.
  • (32) S. van der Walt, S. C. Colbert, and G. Varoquaux. The numpy array: A structure for efficient numerical computation. Computing in Science Engineering, 13(2):22–30, March 2011.
  • (33) H. Vandeven. Family of spectral filters for discontinuous problems. Journal of Scientific Computing, 6(2):159–192, 1991.
  • (34) H. Wilbraham. On a certain periodic function. Cambridge and Dublin Math. J, 3(198):1848, 1848.