Solve wave equation with discontinuous coefficients numerically?

Question

I would like to solve the following equation $$\frac{\partial^2 y}{\partial t^2} - c^2(x,t)\frac{\partial^2 y}{\partial x^2}=0,$$ for $y=y(x,t)$ numerically. The wave speed, $c(x,t)$, is of the form $$c(x,t) = \begin{cases} c_1(t), & x<0 \\ c_2(t), & x\ge 0. \end{cases}$$ The initial conditions are $$y(x,0)=f(x),$$ $$\left.\frac{\partial y}{\partial t}\right|_{t=0}=g(x).$$ The boundary conditions are that the solution is periodic in $x$, with period $2l$ such that $$y(-l,t)=y(l,t).$$ Note that if you can find a solution with a more convenient set of boundary conditions please let me know. We need $y(x,t)$ and $\partial y / \partial x$ to be continuous.

I am not sure how to solve this. My first thoughts are to solve the problem like this. Let $$y(x,t)=\begin{cases} y_1(x,t), & x<0 \\ y_2(x,t), & x\ge0, \end{cases}$$ where $$\frac{\partial^2 y_1}{\partial t^2} - c_1^2(t)\frac{\partial^2 y_1}{\partial x^2}=0,$$ $$\frac{\partial^2 y_2}{\partial t^2} - c_2^2(t)\frac{\partial^2 y_2}{\partial x^2}=0,$$ where the boundary conditions are now $$y_1(-l,t)=y_2(l,t),$$ $$y_1(0,t)=y_2(0,t),$$ but how do I proceed from here?

Answer 1

Here is a brute force solution that would work no matter what is the discontinuity and nonlinearity in $c(x,t)$.

Write your PDE as a system of two:

$ \dot{y}=z\\ \dot{z}=c^2(x,t) y_{xx} $

Now, discretize it on a uniform spatial grid in x:

$ \vec{x}= [x_0, x_1,..., x_{n-1}] \\ \vec{y}= [y_0, y_1,..., y_{n-1}] \\ \vec{z}= [z_0, z_1,..., z_{n-1}] \\ $

Now the PDE becomes a set of 2n ODEs,

$ \dot{y}_k = z_k \\ \dot{z}_k = c^2(x_k,t) (y_{k+1}+y_{k-1}-2 y_k)/{h}^2 $

where $h$ is the grid spacing and $k \in$ {0,1,...,n-1}, except to account for the periodic boundary condition at k=0 use

$ \dot{z}_k = c^2(x_k,t) (y_{k+1}+y_{n-1}-2 y_k)/{h}^2 $

and at k=n-1 use

$ \dot{z}_k = c^2(x_k,t) (y_{0}+y_{k-1}-2 y_k)/{h}^2 $

Put this system of 2n ODEs for the state vector $[\vec{y},\vec{z}]$ into your favorite ODE solver, with initial conditions $[f(\vec{x}),g(\vec{x})]$, and that's it.

A caveat in this approach is that the solution is treated as smooth everywhere although it would not be actually smooth at the $c(x)$ discontinuity; there the solution would have discontinuous first derivatives. This may or may not cause a problem for the numerical solution, depending on how strong the discontinuity is and what kind of solution is sought. But this can affect the grid convergence rate of the solution, and, more importantly, the accuracy of the wavefront refraction angle at the discontinuity. However, there is a simple remedy for all those problems: replace the discontinuity by a smooth resolvable transition layer, e.g., near the discontinuity at x=0 use

$ c(x) = \frac{1}{2}(c_1+c_2) + \frac{1}{2}(c_2-c_1) \frac{x}{\delta} $

where $\delta$ is the transition layer width.

Answer 2

$c$ depending on time is not the issue. You will use an RK scheme which takes care of this. The issue is $c$ is discontinuous in $x$. I recommend SBP-SAT schemes for this.

(1) Derive an energy equation at PDE level. (2) Search literature for SBP-SAT schemes which enforce interface conditions via SAT penalty terms, which are designed to mimic the energy equation at numerical level. E.g., see doi.org/10.1007/s10915-016-0297-3 and in particular, see section "The One Dimensional Wave Equation with a Grid Interface" which is similar to your setting.

I think you can also find exact solution using separation of variables but I have not worked it out. Solve the wave equation in each sub-domain $$ y_i(x,t) = T_i(t) X_i(x), \qquad i=1,2 $$ $$ \frac{T''_i(t)}{c^2_i(t)T_i(t)} = \frac{X''_i(x)}{X_i(x)} = \textrm{constant} $$ Then match the two solutions to have continuity of solution and derivative at $x=0$ and at periodic boundaries. The difficulty is solving $T_i(t)$ since it could be nonlinear equation if $c_i(t)$ is not constant.

Update: If speed is discontinuous, one has to be careful about the form of the pde model. The correct model should be of the form $$ u_{tt} = (c^2 u_x)_x $$ See doi:10.1016/j.jcp.2008.06.023