Finite difference methods for 2D and 3D wave equations

A natural next step is to consider extensions of the methods for various variants of the one-dimensional wave equation to two-dimensional (2D) and three-dimensional (3D) versions of the wave equation.

Multi-dimensional wave equations

The general wave equation in \(d\) space dimensions, with constant wave velocity \(c\), can be written in the compact form

(1)\[ \frac{\partial^2 u}{\partial t^2} = c^2\nabla^2 u\hbox{ for }\boldsymbol{x}\in\Omega\subset\mathbb{R}^d,\ t\in (0,T] {\thinspace .}\]

In a 2D problem (\(d=2\)),

\[\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} ,\]

while in three space dimensions (\(d=3\)),

\[\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} {\thinspace .}\]

Many applications involve variable coefficients, and the general wave equation in \(d\) dimensions is in this case written as

(2)\[ \varrho\frac{\partial^2 u}{\partial t^2} = \nabla\cdot (q\nabla u) + f\hbox{ for }\boldsymbol{x}\in\Omega\subset\mathbb{R}^d,\ t\in (0,T],\]

which in 2D becomes

\[\varrho(x,y) \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x}\left( q(x,y) \frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left( q(x,y) \frac{\partial u}{\partial y}\right) + f(x,y,t) {\thinspace .}\]

To save some writing and space we may use the index notation, where subscript \(t\), \(x\), \(y\), or \(z\) means differentiation with respect to that coordinate. For example,

\[\begin{split}\frac{\partial^2 u}{\partial t^2} &= u_{tt},\\ \frac{\partial}{\partial y}\left( q(x,y) \frac{\partial u}{\partial y}\right) &= (q u_y)_y {\thinspace .}\end{split}\]

The 3D versions of the two model PDEs, with and without variable coefficients, can with now with the aid of the index notation for differentiation be stated as

(3)\[ u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f,\]

(4)\[ \varrho u_{tt} = (q u_x)_x + (q u_z)_z + (q u_z)_z + f\]\[ {\thinspace .}\]

At each point of the boundary \(\partial\Omega\) of \(\Omega\) we need one boundary condition involving the unknown \(u\). The boundary conditions are of three principal types:

1. \(u\) is prescribed (\(u=0\) or a known time variation for an incoming wave),
2. \(\partial u/\partial n = \boldsymbol{n}\cdot\nabla u\) prescribed (zero for reflecting boundaries),
3. an open boundary condition (also called radiation condition) is specified to let waves travel undisturbed out of the domain, see wave:app:exer:tsunami1D:radiation for details.

All the listed wave equations with second-order derivatives in time need two initial conditions:

1. \(u=I\),
2. \(u_t = V\).

Mesh

We introduce a mesh in time and in space. The mesh in time consists of time points

\[\begin{split}t_0=0 < t_1 <\cdots < t_{N_t},\end{split}\]

often with a constant spacing \(\Delta t= t_{n+1}-t_{n}\), \(n\in{{\mathcal{I^-}_t}}\).

Finite difference methods are easy to implement on simple rectangle- or box-shaped domains. More complicated shapes of the domain require substantially more advanced techniques and implementational efforts. On a rectangle- or box-shaped domain mesh points are introduced separately in the various space directions:

\[\begin{split}&x_0 < x_1 <\cdots < x_{N_x} \hbox{ in }x \hbox{ direction},\\ &y_0 < y_1 <\cdots < y_{N_y} \hbox{ in }y \hbox{ direction},\\ &z_0 < z_1 <\cdots < z_{N_z} \hbox{ in }z \hbox{ direction}{\thinspace .}\end{split}\]

We can write a general mesh point as \((x_i,y_j,z_k,t_n)\), with \(i\in{\mathcal{I}_x}\), \(j\in{\mathcal{I}_y}\), \(k\in{\mathcal{I}_z}\), and \(n\in{\mathcal{I}_t}\).

It is a very common choice to use constant mesh spacings: \(\Delta x = x_{i+1}-x_{i}\), \(i\in{{\mathcal{I^-}_x}}\), \(\Delta y = y_{j+1}-y_{j}\), \(j\in{{\mathcal{I^-}_y}}\), and \(\Delta z = z_{k+1}-z_{k}\), \(k\in{{\mathcal{I^-}_z}}\). With equal mesh spacings one often introduces \(h = \Delta x = \Delta y =\Delta z\).

The unknown \(u\) at mesh point \((x_i,y_j,z_k,t_n)\) is denoted by \(u^{n}_{i,j,k}\). In 2D problems we just skip the \(z\) coordinate (by assuming no variation in that direction: \(\partial/\partial z=0\)) and write \(u^n_{i,j}\).

Discretization

Two- and three-dimensional wave equations are easily discretized by assembling building blocks for discretization of 1D wave equations, because the multi-dimensional versions just contain terms of the same type that occurs in 1D.

Discretizing the PDEs

Equation (3) can be discretized as

\[[D_tD_t u = c^2(D_xD_x u + D_yD_yu + D_zD_z u) + f]^n_{i,j,k} {\thinspace .}\]

A 2D version might be instructive to write out in detail:

\[[D_tD_t u = c^2(D_xD_x u + D_yD_yu) + f]^n_{i,j,k},\]

which becomes

\[\frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2} = c^2 \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2} + c^2 \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2} + f^n_{i,j},\]

Assuming as usual that all values at the time levels \(n\) and \(n-1\) are known, we can solve for the only unknown \(u^{n+1}_{i,j}\). The result can be compactly written as

(5)\[ u^{n+1}_{i,j} = 2u^n_{i,j} + u^{n-1}_{i,j} + c^2\Delta t^2[D_xD_x u + D_yD_y u]^n_{i,j}{\thinspace .}\]

As in the 1D case, we need to develop a special formula for \(u^1_{i,j}\) where we combine the general scheme for \(u^{n+1}_{i,j}\), when \(n=0\), with the discretization of the initial condition:

\[[D_{2t}u = V]^0_{i,j}\quad\Rightarrow\quad u^{-1}_{i,j} = u^1_{i,j} - 2\Delta t V_{i,j} {\thinspace .}\]

The result becomes, in compact form,

(6)\[ u^{n+1}_{i,j} = u^n_{i,j} -2\Delta V_{i,j} + {\frac{1}{2}} c^2\Delta t^2[D_xD_x u + D_yD_y u]^n_{i,j}{\thinspace .}\]

The PDE (4) with variable coefficients is discretized term by term using the corresponding elements from the 1D case:

\[[\varrho D_tD_t u = (D_x\overline{q}^x D_x u + D_y\overline{q}^y D_yu + D_z\overline{q}^z D_z u) + f]^n_{i,j,k} {\thinspace .}\]

When written out and solved for the unknown \(u^{n+1}_{i,j,k}\), one gets the scheme

\[\begin{split}u^{n+1}_{i,j,k} &= - u^{n-1}_{i,j,k} + 2u^{n}_{i,j,k} + \\ &= \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \frac{1}{2}(q_{i,j,k} + q_{i+1,j,k})(u^{n}_{i+1,j,k} - u^{n}_{i,j,k}) - \\ &\qquad\quad \frac{1}{2}(q_{i-1,j,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i-1,j,k})) + \\ &= \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \frac{1}{2}(q_{i,j,k} + q_{i,j+1,k})(u^{n}_{i,j+1,k} - u^{n}_{i,j,k}) - \\ &\qquad\quad\frac{1}{2}(q_{i,j-1,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j-1,k})) + \\ &= \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \frac{1}{2}(q_{i,j,k} + q_{i,j,k+1})(u^{n}_{i,j,k+1} - u^{n}_{i,j,k}) -\\ &\qquad\quad \frac{1}{2}(q_{i,j,k-1} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j,k-1})) + \\ + &\qquad \Delta t^2 f^n_{i,j,k} {\thinspace .}\end{split}\]

Also here we need to develop a special formula for \(u^1_{i,j,k}\) by combining the scheme for \(n=0\) with the discrete initial condition, which is just a matter of inserting \(u^{-1}_{i,j,k}=u^1_{i,j,k} - 2\Delta tV_{i,j,k}\) in the scheme and solving for \(u^1_{i,j,k}\).

Handling boundary conditions where is \(u\) known

The schemes listed above are valid for the internal points in the mesh. After updating these, we need to visit all the mesh points at the boundaries and set the prescribed \(u\) value.

Discretizing the Neumann condition

The condition \(\partial u/\partial n = 0\) was implemented in 1D by discretizing it with a \(D_{2x}u\) centered difference, and thereafter eliminating the fictitious \(u\) point outside the mesh by using the general scheme at the boundary point. Alternatively, one can introduce ghost cells and update a ghost value to for use in the Neumann condition. Exactly the same ideas are reused in multi dimensions.

Consider \(\partial u/\partial n = 0\) at a boundary \(y=0\). The normal direction is then in \(-y\) direction, so

\[\frac{\partial u}{\partial n} = -\frac{\partial u}{\partial y},\]

and we set

\[[-D_{2y} u = 0]^n_{i,0}\quad\Rightarrow\quad \frac{u^n_{i,1}-u^n_{i,-1}}{2\Delta y} = 0 {\thinspace .}\]

From this it follows that \(u^n_{i,-1}=u^n_{i,1}\). The discretized PDE at the boundary point \((i,0)\) reads

\[\frac{u^{n+1}_{i,0} - 2u^{n}_{i,0} + u^{n-1}_{i,0}}{\Delta t^2} = c^2 \frac{u^{n}_{i+1,0} - 2u^{n}_{i,0} + u^{n}_{i-1,0}}{\Delta x^2} + c^2 \frac{u^{n}_{i,1} - 2u^{n}_{i,0} + u^{n}_{i,-1}}{\Delta y^2} + f^n_{i,j},\]

We can then just insert \(u^1_{i,1}\) for \(u^n_{i,-1}\) in this equation and then solve for the boundary value \(u^{n+1}_{i,0}\) as done in 1D.

From these calculations, we see a pattern: the general scheme applies at the boundary \(j=0\) too if we just replace \(j-1\) by \(j+1\). Such a pattern is particularly useful for implementations. The details follow from the explained 1D case in the section Implementation of Neumann conditions.

The alternative approach to eliminating fictitious values outside the mesh is to have \(u^n_{i,-1}\) available as a ghost value. The mesh is extended with one extra line (2D) or plane (3D) of ghost cells at a Neumann boundary. In the present example it means that we need a line ghost cells below the \(y\) axis. The ghost values must be updated according to \(u^{n+1}_{i,-1}=u^{n+1}_{i,1}\).