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

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?

• Can we assume c has at most a finite # of jump discontinuities? – Nawaf Bou-Rabee Oct 18 '16 at 16:46
• Yes, that it what was intended. – P Gibson Oct 18 '16 at 18:34

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)$.
• The condition that $f$ be $C^4$ will be fulfilled with sufficiently smooth initial conditions, but here the initial conditions are distributional. The aim is to directly compute the Green's function if possible, including its purely distributional part. – P Gibson Oct 19 '16 at 13:21
• One can smoothly approximate the initial conditions using a regularization of the $\delta$ function. – Nawaf Bou-Rabee Oct 19 '16 at 14:12