## 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 :

(1) |

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

(2) |

and

(3) |

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 .^{5}^{5}5The
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 points^{6}^{6}6The optimal choice of collocation
points depends on the spectral basis chosen. and use the Lagrange
polynomials

(4) |

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

(5) |

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

(6) |

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

(7) |

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

(8) |

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

(9) |

### i.2 Influence of discontinuities

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

(10) |

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:

(11) |

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

(12) |

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 :

(13) |

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:

(14) |

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:

(15) |

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:

(16) |

He then takes their inverse Fourier transform

(17) |

to define a global mollifier :

(18) |

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

(19) |

where is given by

(20) |

for a free parameter .^{7}^{7}7 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
.^{8}^{8}8A 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:

(21) |

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

(22) |

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

(23) |

for all points in the domain Tadmor2007 .^{9}^{9}9Within 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:

(24) |

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

(25) |

where

(26) |

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.#### 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:

(27) |

where

(28) |

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.

## 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

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

(29) |

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:

(30) |

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,

(31) |

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

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:

(32) |

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

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:

(33) |

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

(34) |

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(35) |

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

(36) |

(37) |

and

(38) |

where

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.

### iii.2 Efficacy

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.

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.

## Iv Mollification

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 ^{10}^{10}10 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

(39) |

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.^{11}^{11}11We 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.

### 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

(40) |

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 .^{12}^{12}12Gelb2001EnhancedSph
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.^{13}^{13}13We 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.

## 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.^{14}^{14}14And 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.^{15}^{15}15We 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 .

## References

- (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, https://doi.org/10.1175/1520-0493(2001)129¡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–. http://www.scipy.org/.
- (21) L. Kidder, H. Pfeiffer, M. Scheel, et al. Spectral Einstein Code, 2000. http://www.black-holes.org/SpEC.html.
- (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. https://bitbucket.org/Yurlungur/spectral-shock-capturing.
- (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. www.classicalrealanalysis.com., 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.