Discrete wave simulation - absorbing boundaries?

Question

I wrote a simple 2D wave simulation using the following equations: $$\frac{\partial^2 u}{\partial t^2}=c^2\nabla^2u$$ Where $\nabla^2$ is the discrete laplace operator using a Von Neumann neighborhood lattice.
$$\mathbf {D}^2_{xy} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix}$$ This gives the discrete formula: $$\Delta{}^t_{t+1}u_{x,y} - \Delta^{t}_{t-1}u_{x,y} = c^2(u_{x+1,y,t}+u_{x-1,y,t}+u_{x,y+1,t}+u_{x,y-1,t}-4u_{x,y,t})$$ Rearranging the terms gives the solution for the next timestep: $$u_{x,y,t+1} = c^2(u_{x+1,y,t}+u_{x-1,y,t}+u_{x,y+1,t}+u_{x,y-1,t}-4u_{x,y,t}) + \Delta^{t}_{t-1}u_{x,y} + u_{x,y,t}$$

This seems to work well, but I am currently taking $u=0$ when out of boundary. (which is apparently called the Dirichlet condition)

$$A_{m\times n}\\ \begin{cases} 0, & \text{if}\ (i<0) \lor (j<0) \lor (i\geq m) \lor (j\geq n) \\ A_{ij}, & \text{otherwise} \end{cases}$$

But according to the wave equation, that will make the boundaries act like walls, which waves can bounce on, like this:

Animated Example

But I cannot figure out how to make the boundaries "pass through" as if the grid were infinite. Other simulators do it, but I could not figure out the equations for it.

Example of a simulation that has absorbing boundaries, the waves do not reflect on the borders.

Is there a simple method that I can apply to my currently used equation? I don't need anything fancy, only something simple to apply and fast to compute.

Answer 1

Overview

First a few comments on absorbing boundaries: your second picture seems to imply that it's really the edges of your square which absorb the wave, but this is not quite true. Rather, one can consider your picture as in the following:

That is, you see only the excerpt of the full space in which the wave is considered as undisturbed. And outside of what you see is something happening. Most straightforwardly, that would be a grid large enough so that no reflections can occur. But depending on the parameters, this is numerically too costly.

For this case, one can apply "absorbing" techniques whose general idea is to make the solution vanish outside the region of interest, without changing the behavior inside (the outside boundary is here sketched as a circle, which might be more appropriate due to symmetry, but depending on your wave this doesn't have to be). In order to do so, a balance must be found between (i) having as little impact on the inner region as possible and (ii) an as small outer region as possible.

There are many different techniques to implement the absorbing which try to find such an appropriate balance, among them (complex) absorbing potential, complex exterior scaling, masking functions, etc. (see also the references in the comments).

Masking functions / absorbing potentials

As the OP asked for a simple method, I would suggest to try a masking function: for this you apply a given function $\text{mask}(x,y)$ to your solution at each time step which is smaller than in the outside region $R_\text{out}$. A suitable function could be

$$ f_\text{mask}(x,y) = \tanh^2\left(\lambda (R - \sqrt{x^2+y^2})\right) $$

where $R$ denotes the maximum radius on the grid and the parameter $\lambda$ determines the width of the mask which is to be tuned appropriately. Many other choices are possible here as well, e.g. a suitably streched quarter wavelength of a cosine, etc.

There are basically two alternatives how the masking function is applied. The simple version is to simply multiply the solution $u(x,y)$ by $\text{mask}(x,y)$ at each timestep. More sophisticated, one can include the logarithm of the masking function to the partial differential equation, $$\frac{\partial u}{\partial t}=v + \log(f_\text{mask}(x,y)) u\\ \frac{\partial v}{\partial t} = c^2\nabla^2u $$ Why the log? Because in the operator one needs something which is zero inside, and smaller than zero outside. The latter version is very similar to the (complex) absorbing potential technique.

Although there are more efficient methods available, imo nothing beats the conceptual simplicity of the masking approach.

---

EDIT: the previous applies to first-order-in-time PDEs, the kind of which I'm usually dealing with. For second order equations, represented as two first-order equations for the function and the temporal derivative, the absorbing potential must be applied only to the first equation corresponding to the function. The reason is that one does not want the derivative to be damped out, but the function value.