# Solve wave equation with discontinuous coefficients numerically?

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?

• How would you solve it if $C_{1,2}$ were not dependent on $t$? – Maxim Umansky Aug 4 '20 at 16:23
• (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. – cfdlab Aug 5 '20 at 5:24
• If $c$ was not dependent on $t$ I guess we could then use separation of variables to solve it. – Peanutlex Aug 5 '20 at 11:01
• Right! But for a short interval in time [t,t+dt], we could take c(x,t) as constant in time, and then make an update for the solution over this interval. That seems to be a pathway to solving the whole thing. – Maxim Umansky Aug 5 '20 at 15:31
• @MaximUmansky's hint is a great hint. If $c(x,t)$ were continuous in $x$, you would just discretize the wave equations using some FDM, FEM or FVM or any other method, right? So what is stopping you from doing the same thing when it is discontinuous. Of course, there are some considerations; you want mesh to align with the discontinuity, you probably want to use small time steps in case $c_1$ and $c_2$ are highly oscillatory functions. Also are you looking for a numerical solution or an analytical one? If you are looking for an exact solution, probably there is only a weak one. – Abdullah Ali Sivas Aug 6 '20 at 1:44

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.

• Thank you, do you think this may cause problems at $x=0$ where $c$ is discontinuous? – Peanutlex Aug 6 '20 at 10:29
• I have added to the answer a few words on this. There is potentially a concern there which may or may not matter, depending on the kind of solutions you are looking for; but there is a simple fix for it. – Maxim Umansky Aug 6 '20 at 15:17
• One has to be careful in these cases, brute force approaches may not work, and you must test this before accepting an answer, since future readers may rely on this. 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 for discussion on this, and derivation of jump conditions. – cfdlab Aug 9 '20 at 9:44
• @ cfdlab The model that you are advocating $u_{tt} = (c^2 u_x)_x$ is not the equation we are given here. Are you suggesting that the given equation is not physically consistent at the discontinuity? Either way, the approach proposed here, replacing the discontinuous $c$ by a continuous one, would still apply. Physical phenomenal like wave refraction would still be captured with a continuous rapidly changing $c(x)$. This is similar to using artificial dissipation at the shock wave front, which is one of standard techniques. – Maxim Umansky Aug 9 '20 at 17:16

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

• If the wave speed only varies with time, so the proposed PDE is still correct (since $c_x=0$). But it never bad to keep in mind what the correct PDE is in case one modifies the code later on. Additionally, the obtained ODE's from separation of variables are still linear, it is just that the ODE in $T$ is time varying and not nonlinear (this could still pose difficulties though). – fibonatic Oct 29 '20 at 18:47