# High-order difference and pseudospectral methods for discontinuous problems

###### Abstract

High order finite-difference or spectral methods are typically problematic in approximating a function with a jump discontinuity. Some common remedies come with a cost in accuracy near discontinuities, or in computational cost, or in complexity of implementation. However, for certain classes of problems involving piecewise analytic functions, the jump in the function and its derivatives are known or easy to compute. We show that high-order or spectral accuracy can then be recovered by simply adding to the Lagrange interpolation formula a linear combination of the jumps. Discretizations developed for smooth problems are thus easily extended to nonsmooth problems. Furthermore, in the context of one-dimensional finite-difference or pseudospectral discretizations, numerical integration and differentiation amount to matrix multiplication. We construct the matrices for such operations, in the presence of known discontinuities, by operating on the corrected Lagrange formula. In a method-of-lines framework, this provides a simple and efficient way to obtain solutions with moving discontinuities to evolution partial differential equations.

^{†}

^{†}thanks: This work was supported by DFG grant SFB/Transregio 7 “Gravitational Wave Astronomy”, STFC grant PP/E001025/1, and (FP7/2007-2013)/ERC grant no. 304978.

## I Introduction

Finite-difference and spectral methods are commonly used for numerically solving ordinary or partial differential equations. Such methods are highly efficient and accurate when the solution is smooth. When discontinuities are present, however, oscillations due to Gibbs phenomena arise that contaminate the solution and cause loss of convergence. Several workarounds exist and are commonly used to suppress or avoid these phenomena. However, they typically come at a cost. For piecewise analytic solutions, domain decomposition avoids the point of discontinuity, but can make implementation complicated and computation expensive, especially in time-dependent problems with moving discontinuities, as remapping of coordinates or regridding is required at each time step. For finite-difference methods, lowering the order of approximation globally or near the discontinuity reduces oscillations but sacrifices accuracy. For spectral methods, applying filtering techniques can improve the convergence properties of the solution away from discontinuities, but the accuracy near discontinuities remains poor Boyd2001. Post-processing of the resulting oscillatory data using direct Gottlieb1992; Gottlieb1994; Gottlieb1996; Gottlieb1995; Gottlieb1995a; Gottlieb1997; GottliebDavidSigal2003 and inverse Jung2004 Gegenbauer reconstruction can recover spectrally convergent non-oscillatory solutions. This approach has been successfully used in several contexts, such as magnetic resonance imaging Archibald2003; Archibald2002; Archibald2004 and evolving nonlinear hyperbolic systems Shu1995. Nevertheless, reconstruction techniques are also computationally expensive and subject to complications. For instance, the Gegenbauer transformation matrices become ill-conditioned for high-degree polynomials Jung2004. In addition, the reconstruction must be constrained to use a sufficiently small ratio of order to truncation, which can impair its effectiveness in certain applications Boyd2005.

In many problems, however, the jumps of the function and its derivatives at the discontinuity are known, or easily computed. One may thus construct scheme that exploit this information. Such ideas have been introduced by Krylov Krylov1907, Lanczos Lanczos1966 and Eckhoff eckhoff1994; eckhoff1995; eckhoff1996; eckhoff1997; eckhoff1998; eckhoff2002 in the context of Fourier methods (see also kantorovichkrylov1958; Barkhudaryan2007; Poghosyan2010; Adcock2010; Batenkov2012; Poghosyan2012) and by Lipman and Levin Lipman2009 in the context of finite-difference methods. In what follows, we consider nodal methods based on Lagrange interpolation. We show that adding a linear combination of the jumps to the Lagrange interpolation formula can account for the presence of the discontinuity and restore the original accuracy. Finite-difference and pseudospectral scheme developed for smooth problems, outlined in §2, can thus be very easily extended to discontinuous or nonsmooth problems, outlined in §3. The merits and possible extensions of this method are summarized in §4.

## Ii Smooth functions

As a primer for the discontinuous methods that follow in §3, we first review the Lagrange interpolation method for sufficiently differentiable functions and its role in numerical integration and differentiation.

### ii.1 Interpolation

Let be a function, with its values known at the ordered, distinct nodes (). The collocation polynomial of degree :

(1) |

that interpolates the given nodes, is unique and determined by solving the linear system of algebraic conditions

(2) |

for the coefficients . The solution may be written in the Lagrange form

(3) |

where

(4) |

are the Lagrange basis polynomials. By construction, the basis polynomials satisfy the conditions , with denoting the Kronecker symbol, so that the polynomial (3) satisfies the conditions (2).

The nodes have so far been left unspecified. For polynomials of low degree, equidistant nodes

(5) |

are commonly used. For polynomials of high degree, however, this choice leads to oscillations due to Runge’s phenomenon and the interpolation process diverges.

The oscillations can be minimized by using nodes that are distributed more densely towards the edges of the interval. For instance, interpolation on Chebyshev-Gauss-Lobatto nodes

(6) |

converges uniformly for every absolutely continuous function Fornberg1998. These nodes are the extrema of a Chebyshev polynomial of degree . If all the function derivatives are bounded on the interval , Chebyshev interpolation always has the property of infinite order (i.e. exponential) convergence.

Spectral methods are often introduced and implemented as expansions in terms of orthogonal basis functions. However, as the use of spectral expansions for most practical purposes requires the introduction of a finite grid, one may question the necessity of special basis functions. Indeed, the construction of a spectral expansion, with coefficients determined via collocation at the nodes , can be replaced by a global interpolating polynomial if the nodes are the abscissas of a Gaussian quadrature (typically the simple roots or extrema of a Jacobi polynomial, of which Chebyshev is a special case) associated with the basis functions Fornberg1998; Boyd2001. This pseudospectral approach allows numerical operations to be performed in physical, rather than spectral, space without compromise in convergence order. This property can be exploited in numerical codes, offering the flexibility of switching between finite-difference and pseudospectral methods by simply shifting the location of the nodes (e.g. changing from Eq. (5) to (6)). We shall make use of this flexibility in §3.

### ii.2 Differentiation

The -th derivative of may be approximated by differentiating the Lagrange interpolating polynomial (3) times. This operation may be written as matrix multiplication:

(7) |

where

(8) |

is the derivative matrix, constructed from the
derivative of the basis polynomials
(4)
with respect to . Note that, for , the above expression gives the identity matrix, , as expected.
For fixed grids, the above matrices can be pre-computed and stored in memory. The floating-point precision of Eq. (7) can be improved using a “negative sum trick” Baltensperger2003
which ensures that the derivative matrices vanish upon multiplying a constant vector. Fast methods for computing the matrices
have been developed
Fornberg1998; Fornberg1988; Fornberg2006; Sadiq2011; Welfert1997 and implemented in computational libraries and software. For instance, the Mathematica command

`n=N`

Solve‘FiniteDifferenceDerivative[

`erivative[n],X,`

ifferenceOrder-¿m]

‘‘DifferentiationMatrix’’ \\ \\uses Fornberg’s algorithm to compute the \verb n -th derivative matrix \verb Dn , given the list of grid-points \verb X={x0,x1,x2,...,xN} , at the desired order of approximation \verb m . The latter is an integer satisfying $n \le m \le N$ that yields a finite-difference error $\mathcal{O}(\Delta x ^{m})$, where $\Delta x$ is the maximum local grid spacing. A choice $m<N$ yields a composite derivative matrix of size $N+1$ constructed from centred $(m+1)$-point stencils. For example, the choice $n=1$ and $m=2$ computes $1^{\rm st}$ derivatives of $2^{\rm nd}$ order polynomials (given by Eq.~\eqref{PolyNLagrange1} for $N=2$) to construct the centred finite-difference approximation %\begin{widetext} \begin{eqnarray} \label{smoothFDf} {{f’}_i} &=& {f_{i + 1}}\frac{{{x_i} - {x_{i - 1}}}}{{({x_{i + 1}} - {x_i})({x_{i + 1}} - {x_{i - 1}})}} \nonumber\\&+& {f_i}\frac{{{x_{i + 1}} + {x_{i - 1}} - 2{x_i}}}{{({x_{i + 1}} - {x_i})({x_i} - {x_{i - 1}})}} \nonumber\\ &-& {f_{i - 1}}\frac{{{x_{i + 1}} - {x_i}}}{{({x_{i + 1}} - {x_{i - 1}})({x_i} - {x_{i - 1}})}}+\mathcal{O}(\Delta x^2). \end{eqnarray} %\end{widetext} The entire domain $[a,b]$ is covered in a composite fashion (that is, the above formula is applied repeatedly at each node $x_i$ using only nearest neighbour points). This gives the centred-derivative matrix \begin{eqnarray} D^{(1)}_{ij} &=&\frac{{{x_i} - {x_{i - 1}}}}{{({x_{i + 1}} - {x_i})({x_{i + 1}} - {x_{i - 1}})}}{\delta _{j,i + 1}} \nonumber \\ &+& \frac{{{x_{i + 1}} + {x_{i - 1}} - 2{x_i}}}{{({x_{i + 1}} - {x_i}) ({x_i} - {x_{i - 1}})}}{\delta _{ij}} \nonumber \\ &-& \frac{{{x_{i + 1}} - {x_i}}}{{({x_{i + 1}} - {x_{i - 1}})({x_i} - {x_{i - 1}})}}{\delta _{j,i - 1}} \end{eqnarray} for all $i=1,\dots,N-1$, while one-sided finite differences are used at the end-points $i=0,N$. A single-domain pseudospectral derivative matrix can be obtained by choosing the maximal value $m=N$ and selecting nodes $x_i$ that coincide with the abscissas of a Gaussian quadrature, %associated with the basis functions such as the Chebyshev-Gauss-Lobatto nodes \eqref{ChebyshevNodes}. It is interesting to note that the pseudospectral derivative matrices obtained from this (single-domain) application of Eq.~\eqref{derivativematricespi} are the \textit{infinite-order limit} of the finite-difference derivative matrices obtained by (composite) application of Eq.~\eqref{derivativematricespi} \cite{Fornberg1998}. This adds further support to the view of pseudospectral methods as a special case of finite-difference methods. %$\sum_jD^{(n)}_{ij}$ \subsection{Integration} Quadrature, or numerical integration, amounts to approximating the definite integral of a given function or, equivalently, integrating an approximation to the function. Integrating the approximation \eqref{PolyNLagrange1} gives a quadrature rule in the form of a weighted sum of function evaluations: %\begin{equation} \label{Integralp1} %\int_{{a}}^{{b}} {f(x)dx}\simeq \int_{{a}}^{{b}} {p(x)dx} = \sum\limits_{j = 0}^N {{D^{(-1)}_{ij}}{f_j}} %\end{equation} \begin{equation} \label{Integralp1} \int_{{a}}^{{b}} {f(x)dx}\simeq \int_{{a}}^{{b}} {p(x)dx} = \sum\limits_{j = 0}^N {{w_{j}}{f_j}}, \end{equation} where %\begin{equation} \label{IntegralMatrixp1} %{D^{(-1)}_{ij}}=\int_{{x_0}}^{{x_i}} {\pi_j(x)dx} %\end{equation} \begin{equation} \label{IntegralMatrixp1} {w_{j}}=\int_{{a}}^{{b}} {\pi_j(x)dx} \end{equation} are the relevant weights of integration. Various well-known quadrature rules stem from the above formulae. For instance, using interpolation of 2 or 3 equidistant nodes yields the composite trapezoidal or Simpson rule respectively. For the Chebyshev-Gauss-Lobatto nodes \eqref{ChebyshevNodes}, the above procedure yields the homonymous quadrature formula, which is exponentialy convergent for smooth integrands. \\ \section{Non-smooth functions} In general, a jump discontinuity in a function or its derivatives cannot be accurately approximated by a single polynomial. %($N \le k-1$). As mentioned earlier, a common workaround is to divide the computational domain so that any discontinuities are located at the boundary between domains. In what follows, we explore a simpler approach that does not require domain decomposition but uses the location and magnitude of the discontinuity if known. %However, if the location and magnitude of the discontinuity is known, one may explore a simpler approach that does not require domain decomposition. We will illustrate that incorporating discontinuities into the approximating function instead of avoiding them offers flexibility and ease of programming without loss of accuracy. \subsection{Interpolation} Let $f:[a,b]\rightarrow \mathbb{R}$ be a $C^{k}$ function, with the values $f_i=f(x_i)$ known at the ordered, distinct nodes $x_i$. The Lagrange interpolation method outlined in \S2 requires sufficient differentiability, that is, $k \ge N+1$. However, non-differentiable solutions appear in multiple physical contexts, and their approximation necessitates modification of the above methods. If $f$ is a piecewise-$C^{N+1}$ function and its jump discontinuities are known, then a simple generalization of the Lagrange method to %any $k \ge -1$ any $k \ge -1$ may be obtained as follows. Let the discontinuity be located at some $\xi \in (a,b)$ and the (non-infinite) jumps in $f$ and its derivatives be given by: \begin{equation} \label{fmJconditions2} {f^{(m)}}(\xi^+ ) - {f^{(m)}}(\xi^- ) = {J_m},\quad m = 0,1, \ldots ,\infty \end{equation} where % $M \le N$ and $\xi^-$ and $\xi^+$ denote the left and right limits to $\xi$. One may then approximate $f$ by a piecewise polynomial function that (i) interpolates the given $N+1$ nodes, and (ii) has the same derivative jumps at the discontinuity. Eq. \eqref{PolyN1} is thus replaced by the ansatz \begin{equation} \label{PolyN2DiscInterpol} p(x) %=\sum\limits_{j = 0}^{N} {{c_j}(x){x^j}} =\theta (x - \xi ){p_ + }(x) + \theta (\xi - x){p_ - }(x) \end{equation} where \begin{equation} \label{PolyN1DiscInterpol} \theta (x)=\left\{ {\begin{array}{*{20}{c}} 1,&{x > 0}\\ {1/2,}&{x = 0}\\ 0,&{x < 0} \end{array}} \right. \end{equation} denotes the Heaviside step function %, %\begin{equation} \label{PolyCoefN2} %{c_j}(x) = \theta (x - \xi )c_j^ + + \theta (\xi - x)c_j^ - %\end{equation} %are piecewise constant coefficients, and \begin{equation} \label{PolyN2} {p_ - }(x) = \sum\limits_{j = 0}^{N} {c_j^-(\xi) {x^j}},\quad{p_ + }(x) = \sum\limits_{j = 0}^{N} {c_j^+(\xi) {x^j}} \end{equation} denote the left and right interpolating polynomials. Note that the polynomials $p_-,p_+$ and the piecewise polynomial function $p$ depend implicitly on the location $\xi$ of the discontinuity. % that interpolate $f$ left and right of the discontinuity at $x=\xi$. Half of the ($2N+2$) polynomial coefficients $c_j^\pm$ %$c_j^\pm$ are to be determined by the collocation conditions \eqref{PolyNConds1} which, given the ansatz \eqref{PolyN2DiscInterpol}, translate to \begin{equation} \label{PolyNConds2} \left\{ \begin{array}{*{20}{c}} p_{+} (x_i) = f_i,&x_i > \xi ,\\ p_{-} (x_i) = f_i,&x_i < \xi . \end{array} \right. \end{equation} %This only determines $N+1$ coefficients. The The remaining coefficients could be determined by the first $N+1$ of the jump conditions \eqref{fmJconditions2}. However, as information on high order jumps may not be beneficial (as discussed below and illustrated in Figs.~1 and 2) or easily accessible, we shall allow the number of jumps enforced to be user-specifiable, by dropping jumps in derivatives higher than order $M\le N$, where $M$ is some chosen integer between $-1$ and $N$. % it tends to minimize Runge-type phenomena at high $N$. The remaining $N+1$ coefficients are then determined by the jump conditions \begin{equation} \label{fmJ0conditions2} p_ + ^{(m)}(\xi ) - p_ - ^{(m)}(\xi ) = \left\{ \begin{array}{l} {J_m},\quad m = 0,1,...,M\\ 0,\quad m = M + 1,...,N \end{array} \right. \end{equation} Setting the higher than $M^{\rm th}$-order derivative jumps to zero means that the polynomial coefficients $c_{M+1},\dots, c_N$ remain constant across $\xi$ as in Lagrange interpolation. This choice is favoured by numerical experiments to be discussed below and allows recovery of Lagrange interpolation as a special case, $M=-1$. It should be noted that the ansatz \eqref{PolyN2} should not be viewed merely as two independent interpolating polynomials matched at the discontinuity. If that were the case, the polynomial degrees would differ left and right of the discontinuity, depending on the number of grid points on each side, and our method would amount to domain decomposition. Instead, our ansatz should be regarded as a single ‘‘polynomial’’ of degree $N$ that covers the whole domain $[a,b]$ and has piecewise constant coefficients with a jump at the discontinuity $\xi$. The left and right coefficients encode information from \textit{all} nodes and are correlated via the jump conditions \eqref{fmJ0conditions2}. As mentioned above, the ‘‘no jump’’ choice $M=-1$ corresponds to Lagrange interpolation, which suffices for smooth functions. %Very small values of $M$ generally yield smoother approximations. A choice of $M$ at least as high as $k+1 = 0,1,2\dots$ is required to approximate $C^{k}$ functions with accuracy better than Lagrange, as this reproduces the first nonvanishing jump, $J_{k+1}$, in the $(k+1)^{\rm th}$ derivative. (For instance, a $C^0$ function can be approximated by specifying at least the first derivative jump $f’(\xi^+)-f’(\xi^-)=J_1$.) Incorporating higher derivative jumps, that is, using $M > k+1 $, increases the approximation accuracy of the function as the jumps $J_{k+2}, J_{k+3},\dots ,J_{M}$ in its $(k+2)^{\rm th},(k+3)^{\rm th}, \dots , M^{\rm th}$ derivative are also reproduced.\footnote{Such behaviour is similar to Hermite \cite{hermite1878} and Birkhoff \cite{Birkhoff1906} interpolation, except the present method only uses derivative \textit{jumps} at $\xi$ rather than derivatives themselves at $x_i$. } Nevertheless, the accuracy degrades for values of $M$ too close to $N$ due to the appearance of spurious oscillations analogous to Runge’s phenomenon. Numerical experiments suggest that intermediate choices of $M$ relative to $N$ (typically $M/N \approx 1/2 $ depending on the function) are optimal and consistently yield highly accurate results. %For small values of $M$, the behaviour of the approximation approaches that of Lagrange interpolation. The latter is recovered with the choice $M=-1$, whence only the second branch of Eq. \eqref{fmJ0conditions2} survives and forces the coefficients of the left and right polynomials \eqref{PolyN2} to coincide. As detailed above, the linear system of algebraic conditions \eqref{PolyNConds2}-\eqref{fmJ0conditions2} can be uniquely solved by substituting the polynomials \eqref{PolyN2} and determining their unknown coefficients. This procedure is simplified if the polynomials are first converted to a Lagrange basis, \begin{equation} \label{ppmCpi} p_\pm(x)=\sum_{j=0}^N C^\pm_j(\xi) \pi_j(x). \end{equation} Substitution in the algebraic conditions \eqref{PolyNConds2}-\eqref{fmJ0conditions2} yields \begin{equation} \label{Cpmxifthetag} C^+_j(\xi) =f_j+\theta(\xi - x_j)g_j(\xi),\quad C^-_j(\xi) =f_j-\theta(x_j-\xi)g_j(\xi), \end{equation} where \begin{equation} \label{gjofx2} {g_j}(\xi) = \sum\limits_{m = 0}^{M } {\frac{{{J_m}}}{{m!}}{{({x_j} - \xi )}^m}} \end{equation} are weights computed from the jump conditions at the discontinuity $\xi$ given the nodes $x_j$. Substituting %When $J_m \rightarrow 0$, these weights vanish and one recovers Lagrange interpolation. Eqs. \eqref{ppmCpi}-\eqref{gjofx2} into the ansatz \eqref{PolyN2DiscInterpol} yields the concise solution \begin{equation} \label{PolyNLagrange2} p(x) = \sum\limits_{j = 0}^N {{[f_j+s_j(x;\xi)]}{\pi _j}(x)} \end{equation} which generalizes the Lagrange formula \eqref{PolyNLagrange1} to a \textit{piecewise} \textit{polynomial} interpolating function. This interpolating function depends implicitly on the location $\xi$ of the discontinuity through the weight functions \begin{equation} \label{sjofx2} {s_j}(x;\xi) = [\theta (x - \xi )\theta (\xi - {x_j}) - \theta (\xi - x)\theta ({x_j} - \xi )]g_j(\xi). \end{equation} When $J_m = 0$, the above quantities vanish and Lagrange interpolation is recovered. When evaluated at a node $x=x_i$, the above expression simplifies to \begin{equation} \label{sjofx2b} {s_j}(x_i;\xi) = [\theta (x_{i} - \xi ) - \theta ({x_j} - \xi )]g_j(\xi) \end{equation} Since the prefactor in ${s_j}(x_i;\xi)$ is antisymmetric in its two indices and $\pi_j(x_i)=\delta_{ij}$ is symmetric, the correction $\sum s_j(x;\xi)\pi_j(x) $ in the interpolating formula \eqref{PolyNLagrange2} vanishes at each node $x=x_i$. This ensures that the collocation conditions \eqref{PolyNConds1} or, equivalently, \eqref{PolyNConds2}, are satisfied by construction. %It is easy to confirm that the piecewise polynomial %\eqref{PolyNLagrange2} satisfies the conditions \eqref{PolyNConds1}, since %the contribution $\sum g_j(x_i)\pi_j(x_i)=\sum g_j(x_i)\delta_{ij}$ %vanishes at the nodes. %\begin{equation} %(1-x^2)\Phi’’(x)-2x \Phi’(x) +l(l+1) \Phi(x)=\delta(x-\xi) %\end{equation} As an example, let us consider the piecewise smooth solutions \begin{equation} \label{PhilLegendrePQ} \Phi_l(x)={P_l}(\xi ){Q_l}(x)\theta (x - \xi ) + {P_l}(x){Q_l}(\xi )\theta (\xi-x) \end{equation} to the inhomogeneous Legendre differential equation with a singular (Dirac delta function) source, where $P_l(x)$ and $Q_l(x)$ denote the Legendre functions of the first and second kind respectively. Such solutions arise in the context of self-forces acting on static scalar and electric test charges in the spacetime of a Schwarzschild black hole \cite{Burko2000}. In this case, the derivative jumps are given by \begin{equation} \label{Phijumps} J_k=\Phi^{(k)}_l(\xi^+)-\Phi^{(k)}_l(\xi^-)={P_l}(\xi )Q^{(k)}_l(\xi) - P^{(k)}_l(\xi){Q_l}(\xi ) \end{equation} In general, although such solutions have finite differentiability at $x=\xi$, their derivative jumps can be obtained analytically from the field equation they satisfy without knowledge of the actual solution. This knowledge of the jumps can be used to inform the method outlined above and obtain numerical solutions to the radial part of the field equations. The accuracy of the numerical solutions depends on how well the underlying interpolating function (used for numerical operations such as integration or differentiation) approximates the actual solution. A comparison of the analytic solution \eqref{PhilLegendrePQ} to interpolating functions based on the Lagrange formula \eqref{PolyNLagrange1} and its generalization, Eq. \eqref{PolyNLagrange2}, is depicted in Fig.~1 for equidistant and Fig.~2 for pseudospectral nodes. Although both formulas exhibit Runge oscillations near the boundaries for equidistant nodes, only the corrected formula faithfully reproduces the discontinuity given sufficient resolution. Pseudospectral nodes eliminate the Runge oscillations near the boundaries, but the corrected formula manages to also eliminate the Gibbs oscillations of the Lagrange formula near the discontinuity. Evidently, \textit{the accuracy lost when applying the Lagrange formula on nonsmooth functions is recovered by simply adding the correction terms} \eqref{sjofx2}. A comparison of convergence rates using pseudospectral nodes and different numbers of derivative jumps is depicted in Fig.~3. As mentioned earlier, the convergence estimates favour the use of an intermediate number of derivative jumps relative to a given number of nodes ($M\approx N/2$). With pseudospectral nodes, the formula \eqref{PolyNLagrange2} is at least $M^{\rm {th}}$-order convergent if derivatives of $f$ higher than $M$ have finite jumps, such as the example given above. For functions whose derivatives higher than order $M$ are continuous, the formula \eqref{PolyNLagrange2} is exponentially convergent. \begin{figure*}[h] \includegraphics[width=0.65\textwidth]{fig1} \includegraphics[width=0.65\textwidth]{fig2} %\vspace{2.5in} \caption{\textbf{Top}: Single domain interpolation of the function ${f}(x) = {P_2}(\xi ){Q_2}(x)\theta (x - \xi ) + {P_2}(x){Q_2}(\xi )\theta (\xi-x)$ at $N+1=13$ equidistant nodes, using the Lagrange formula~\eqref{PolyNLagrange1} or its discontinuous generalization \eqref{PolyNLagrange2} with jump conditions on the first $M=6$ or $12$ derivatives. Although all formulas exhibit Runge oscillations near the boundaries, only Eq.~\eqref{PolyNLagrange2} faithfully reproduces the function near the discontinuity. \textbf{Bottom}: Interpolation error. Runge oscillations are minimized when an intermediate number of jumps ($M\approx N/2$) is used.} \label{diff} \end{figure*} \subsection{Differentiation} A finite-difference or pseudospectral approximation to the $n$-th derivative of $f$ can be obtained by differentiating the piecewise polynomial \eqref{PolyNLagrange2}. When evaluated at a node, $x=x_i$, this yields \begin{equation} \label{fnderivativeapprox2} {f^{(n)}}({x_i}) \simeq {p^{(n)}}({x_i }) = \sum\limits_{j = 0}^N D_{ij}^{(n)}[f_j+s_j(x_{i};\xi)] \end{equation} with the derivative matrices $D_{ij}^{(n)}$ given again by Eq.~\eqref{derivativematricespi}. Again, high-order or exponential convergence can be attained thanks to the correction terms \eqref{sjofx2b}. As an example, let us consider a function with jumps $J_0, J_1$ and $J_2$ in its $0^{\rm{th}}$, $1^{\rm{st}}$ and $2^{\rm{nd}}$ derivative at the $J^{\rm{th}}$ node, $x_J=\xi$. Then, for $3$-point finite-difference stencils, Eq.~\eqref{fnderivativeapprox2} gives the left and right derivatives: \begin{eqnarray} \label{fnderivativeleftright} f’_{J^{\pm}}&=& f_J \\&-& \frac{1}{2}{J_0}\left( {\frac{1}{{{x_J} - {x_{J - 1}}}} + \frac{1}{{{x_{J + 1}} - {x_J}}} - \frac{2}{{{x_{J + 1}} - {x_{J - 1}}}}} \right)\nonumber\\ &+&J_1\frac{{{x_{J \pm 1}} - {x_J}}}{{{x_{J + 1}} - {x_{J - 1}}}} - \frac{1}{2}{J_2} \frac{{({x_{J + 1}} - {x_J})({x_J} - {x_{J - 1}})}}{{{x_{J + 1}} - {x_{J - 1}}}} \nonumber \end{eqnarray} with $f_J$ given by Eq.~\eqref{smoothFDf} for $i= J$. Notice that the above formula satisfies the jump condition $f’_{J^{+}}-f’_{J^{-}}=J_1$ as it ought to. Higher order or pseudospectral formulas can be obtained in a similar fashion. For instance, using the nodes given by Eq.~\eqref{ChebyshevNodes}, the formula \eqref{fnderivativeapprox2} yields a Gauss-Chebyshev-Lobatto scheme that is at least of $M^{\rm {th}}$-order convergence if derivatives of $f$ higher than $M$ have finite jumps, or of infinite-order convergence if derivatives higher than $M$ are continuous. \subsection{Integration} Numerical integration formulas for nonsmooth integrands can be obtained, as in \S2.3, by integrating the interpolating formula \eqref{PolyNLagrange2}. Then, the standard quadrature formula \eqref{Integralp1} generalizes to %\begin{equation} \label{Integralp2} %\int_{{x_0}}^{{x_i}} {f(x)dx}\simeq \int_{{x_0}}^{{x_i}} {p(x)dx} = \sum\limits_{j = 0}^N {{D^{(-1)}_{ij}}[f_j+s_j(x_{i})]} %\end{equation} \begin{equation} \label{Integralp2} \int_{{a}}^{{b}} {f(x)dx}\simeq \int_{{a}}^{{b}} {p(x)dx} = \sum\limits_{j = 0}^N {{[w_{j}}{f_j}+q_j(\xi)]} \end{equation} where the weights $w_j$ are given by Eq. \eqref{IntegralMatrixp1}, and the corrections \begin{align} \label{IntegralMatrixp2} {q_{j}(\xi)}&=\int_{{a}}^{{b}} {s_j(x;\xi)\pi_j(x)dx} \\ &= g_j \theta(\xi-x_j) \int_{{\xi}}^{{b}} {\pi_j(x)dx}-g_j\theta(x_j-\xi)\int_{{a}}^{{\xi}} {\pi_j(x)dx} \nonumber \end{align} account for the discontinuity and vanish when $J_m=0$. Eq.~\eqref{Integralp2} can be used to obtain discontinuous generalizations of standard composite or Gaussian quadrature formulas using, for example, the nodes \eqref{EquidistantNodes} or \eqref{ChebyshevNodes} respectively. This is very useful but relatively straightforward and thus not detailed here. %\begin{widetext}%\centering \begin{figure*}[h] %\begin{center} \includegraphics[width=0.65\textwidth]{fig3} \includegraphics[width=0.65\textwidth]{fig4} %\vspace{2.5in} \caption{Same as Fig. 1 but using Chebyshev-Gauss-Lobatto nodes given by Eq.~\eqref{ChebyshevNodes}, which minimize the maximum error by distributing oscillations more evenly. Eq.~\eqref{PolyNLagrange2} faithfully reproduces the function near and away from the discontinuity, especially when an intermediate number ($M\approx N/2$) of derivative jumps are enforced.} \label{diff2} %\end{center} \end{figure*} %\end{widetext} \clearpage \begin{figure}[h] \includegraphics[width=\columnwidth]{fig5} %\vspace{2.5in} \caption{$l_\infty$ error norm of interpolating the function ${f}(x) = {P_2}(\xi ){Q_2}(x)\theta (x - \xi ) + {P_2}(x){Q_2}(\xi )\theta (\xi - x)$ at $N+1$ Chebyshev-Gauss-Lobatto nodes using the Lagrange formula~\eqref{PolyNLagrange1} or its discontinuous generalization \eqref{PolyNLagrange2} with jump conditions on the first $M=5, 15$ or $30$ derivatives enforced. Having dropped the (nonzero in this case) jumps in derivatives higher than $M$, we observe convergence of order slightly above $M$. For functions with continuity in derivatives higher than $M$, infinite-order (exponential) convergence is attained.} \label{diff3} \end{figure} \section{Concluding remarks} We have shown that, if the derivative jumps across a discontinuity are known, the simple modification \eqref{PolyNLagrange2} to the Lagrange formula \eqref{PolyNLagrange1} can fully account for the presence of a discontinuity and completely avoid Gibbs phenomena. This can be exploited in finite-difference or pseudospectral scheme to recover their original accuracy. The main advantage offered by the present method is that it is efficient and simple to implement, even when the discontinuity moves freely inside the domain, as no requirements for domain decomposition, regridding, or reconstruction arise. Accurate numerical differentiation of piecewise smooth functions is crucial to many problems of physical interest. A common scenario consists of evolving a time-dependent partial differential equation via the \textit{method of lines} \cite{knapp2008method,Schiesser1991}. Spatial operators are then computed via finite- difference or pseudospectral methods described above, and time integration is performed using standard numerical methods for coupled ordinary differential equations, such as Runge-Kutta. The standard methods of spatial discretization in \S2 cease to apply when discontinuities are present; this is when Eq.~\eqref{fnderivativeapprox2} becomes useful. For a moving discontinuity, the weights \eqref{sjofx2b} must by evaluated at each time step; if the jumps $J_m$ at $\xi$ are known, this evaluation is not expensive computationally as it requires a few subtractions and multiplications but no extra function evaluations. The cost and accuracy of obtaining the jumps and their location depends on the problem at hand; in certain problems this information is available ‘‘for free’’. For instance, in calculations of radiation-reaction forces in electrodynamics or gravitational self-force corrections to geodesic motion of small bodies in general relativity \cite{Barack2009,Poisson2011}, the source of the dimensionally reduced field equations is singular, making the solution or its derivative discontinuous across the particle location \cite{Canizares2009,Canizares2011a,Canizares2010c,Jaramillo2011}. However, the jumps \eqref{fmJconditions2} occur at the known particle location and can be computed analytically from the radial part of the field equations. The present method is ideally suited to such problems where information on the jumps is available. Lipman and Levin \cite{Lipman2009} have developed a closely related interpolation method and a moving least squares technique for detecting the location and magnitude of an unknown discontinuity or singularity, with accuracy of higher-order than subcell resolution methods \cite{Harten1989}. This technique would allow an extension of the present method to problems with time-dependent discontinuities that are not known a priori. For pseudospectral discretizations, a mapping devised by Kosloff and Tal-Ezer to improve stability of the Chebyshev collocation scheme \cite{Kosloff1993} and the accuracy of higher derivatives \cite{Costa2000,Don1997}, can also be employed in conjunction with Eq.~\eqref{fnderivativeapprox2} in method-of-lines implementations to evolve discontinuous solutions. In problems where the grid values of the function derivatives are available, in addition to their jumps at the discontinuity, Eqs.~\eqref{PolyNConds2} can be complemented by conditions on the derivatives at the same nodes. Together with Eq.~\eqref{fmJ0conditions2}, this will yield generalizations of Hermite \cite{hermite1878} and Birkhoff \cite{Birkhoff1906} interpolation to discontinuous problems. As a final remark, we note that extension of the present technique to higher dimension is technically possible. However, for separable multi-dimensional problems discretized on tensor-product grids, the present method can be applied directly. \section*{Acknowledgments} We thank Sarp Akcay, Sebastiano Bernuzzi, Bernd Br\"ugmann, Priscilla Ca\~nizares, Kyriaki Dionysopoulou, Ian Hawke, David Hilditch, Jonathan Thornburg, Maarten van de Meent and K\=oji Ury\=u for constructive discussions and comments. \bibliographystyle{siam} \bibliography{library} \end{document}