Methods to compute the Green's function for the 1D wave equation with nonsmooth coefficient?

Question

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?

Answer 1

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)$.