I am seeking advice on the best available numerical methods to compute the Green's function for a 1D wave equation with rough coefficient. 


Suppose that the coefficient $c(x)$ in the 1D wave equation $u\_{tt}-c(x)^2u\_{xx}=0$ has constant values $c\_0$ and $c\_1$ to the left and right respectively of a bounded interval $I=(x\_0,x\_1)$, but is variable and nonsmooth on $I$. More precisely, suppose that $c$ is piecewise $C^1$ on $I$ (so that it may have a finite number of jump discontinuities, for example). Consider initial conditions corresponding to a left or right-moving unit impulse at a source point $x\_s$: $u(x,0)=\delta(x-x\_s)$, $u\_t(x,0)=\pm\delta^\prime(x-x\_s)$. One wants to compute the (distributional) solution $u(x\_r,t)$ for $0<t<T$ at some other point $x\_r$. (It may be assumed that $x\_s$ is a point of continuity of $c$.) 


1) What established numerical techniques are applicable to this problem? 


2) What would qualify as the "gold standard" against which new methods should be judged?



Methods to compute the Green's function for the 1D wave equation with nonsmooth coefficient?
As a benchmark for comparison, one can use an explicit variable step size finite difference scheme. This is a well established numerical technique for 1D PDE problems even with unbounded domains. 


I will focus on the spatial discretization of the term $c(x) \partial\_{xx} f(x)$, because the temporal discretization is standard. Even though $c(x)$ may have jump discontinuities, this term is benign because $c(x)$ is not differentiated. Note that numerical stability may require that the local time step size is adjusted in order to ensure that the discrete cone of influence emanating backwards from a given point always contains the continuous one (aka a CFL condition). This is the main point of introducing a variable step size approximation. 


Let $S = \{ x\_i \}$ be a collection of grid points on $\mathbb{R}$ with forward, backward, and average spatial grid sizes given by:
$$
\delta x\_i^+ = x\_{i+1} - x\_i \;, \quad \delta x\_i^- = x\_i - x\_{i-1}\;, \quad \delta x\_i = \frac{\delta x\_i^+ + \delta x\_i^-}{2}
$$
Note that we don't assume that $S$ contains the jump discontinuities of $c(x)$. It might be advantageous to include these points in order to satisfy the above mentioned CFL condition. Also, in order to accurately resolve the effect of the jumps in $c(x)$, one needs to pick the local grid size sufficiently small.


At any grid point (including a jump discontinuity of $c(x)$) use a standard variable grid size, central scheme
$$
(c(x) \partial\_{xx} f(x))\_i \approx \frac{c\_i}{ \delta x\_i} \left( \frac{f\_{i+1} - f\_i}{\delta x\_i^+} - \frac{f\_{i} - f\_{i-1}}{\delta x\_i^-} \right)
$$
If $f$ is $C^4$ (in space), then it is straightforward to verify that this approximation is locally second-order accurate at every grid point. My guess is that the overall accuracy of this scheme in the $\ell\_1$ norm is also second-order. To prove this, one needs to understand the properties of the exact solution in a neighborhood of the jump discontinuities of $c(x)$.


