Home / Finite Difference Methods

The Wave Equation in Time Domain

The wave equation in time domain can be stated as follows, \begin{equation} \frac{\partial^2 u}{\partial t^2}=c^2 \nabla^2 u+f,\qquad x\in \Omega \subset \mathbb{R}^d,~t\in(T_1,T_2] \end{equation} where $\Delta$ is the Laplacian, $f$ is a forcing function (for example our source) and $c$ is the wave velocity at which the time and spatially varying wave $u$ propagates.

A Fourth-Order Finite-Difference Method for
the Three Dimensional Wave Equation in Time Domain

Let $I=(a,b)\times(c,d)\times(e,f)$ and $J=[a,b]\times[c,d]\times[e,f]$. A three-dimensional form of the wave equation can be found as follows, \begin{align} \begin{cases} \dfrac{\partial^2 u}{\partial t^2} = c^2(x,y,z) \Delta u +f(x,y,z,t),&(x,y,z)\in I,~t\in(T_1,T_2],\\ u(x,y,z,T_1)=g(x,y,z),&(x,y,z)\in J,\\ \dfrac{\partial }{\partial t}u(x,y,z,T_1) = s(x,y,z),&(x,y,z)\in J,\\ u(a,y,z,t) = f_{a}(y,z,t),&(y,z)\in [c,d]\times[e,f],~t \in[T_1,T_2],\\ u(b,y,z,t) = f_{b}(y,z,t),&(y,z)\in [c,d]\times[e,f],~t \in[T_1,T_2],\\ u(x,c,z,t) = f_{c}(x,z,t),&(x,z)\in [a,b]\times[e,f],~t \in[T_1,T_2],\\ u(x,d,z,t) = f_{d}(x,z,t),&(x,z)\in [a,b]\times[e,f],~t \in[T_1,T_2],\\ u(x,y,e,t) = f_{e}(x,y,t),&(x,y)\in [a,b]\times[c,d],~t \in[T_1,T_2],\\ u(x,y,f,t) = f_{f}(x,y,t),&(x,y)\in [a,b]\times[c,d],~t \in[T_1,T_2]. \end{cases} \qquad \qquad (1) \end{align}

Consider the following notations:

A Fourth-order Scheme

A fourth-order finite difference scheme for solving the problem (1) can be considered as follows, \begin{align} \left(1+\frac{1}{12}\delta_t^2\right)^{-1}\frac{\delta_t^2}{\tau^2} u_{i,j,k}^{n}&= c^2\left(1+\frac{1}{12}\delta_x^2\right)^{-1}\frac{\delta_x^2}{h_x^2} u_{i,j,k}^{n}+ c^2\left(1+\frac{1}{12}\delta_y^2\right)^{-1}\frac{\delta_y^2}{h_y^2} u_{i,j,k}^{n}\notag \\ &+ c^2\left(1+\frac{1}{12}\delta_z^2\right)^{-1}\frac{\delta_z^2}{h_z^2} u_{i,j,k}^{n}+f_{i,j,k}^{n}. \end{align} Expanding terms, \begin{align} &\left[\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right) -\frac{r_x^2}{12} \left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_x^2 \right. \notag \\ & \left. -\frac{r_y^2}{12}\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_y^2 -\frac{r_z^2}{12}\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\delta_z^2 \right]\delta_t^2u_{i,j,k}^{n}= \notag \\ & \left[r_x^2 \left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_x^2 +r_y^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_y^2 +r_z^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\delta_z^2 \right]u_{i,j,k}^{n} \notag \\ & +\tau^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\left(1+\frac{1}{12}\delta_t^2\right)f_{i,j,k}^{n}, \end{align} where $r_x^2=c^2\frac{\tau^2}{h_x^2}$, $r_y^2=c^2\frac{\tau^2}{h_y^2}$ and $r_z^2=c^2\frac{\tau^2}{h_z^2}$.

Let \begin{align} \rho_x=\frac{1-r_x^2}{12},~\rho_y=\frac{1-r_y^2}{12},~\rho_z=\frac{1-r_z^2}{12}. \end{align}

On the other hand, note that, \begin{align} &\left| \left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)u_{i,j,k}^{n} -\frac{r_x^2}{12} \left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_x^2u_{i,j,k}^{n} -\frac{r_y^2}{12}\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_y^2u_{i,j,k}^{n} \right. \notag \\ & \left. -\frac{r_z^2}{12}\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\delta_z^2u_{i,j,k}^{n} -\left(1+\rho_x\delta_x^2\right)\left(1+\rho_y\delta_y^2\right)\left(1+\rho_z\delta_z^2\right)u_{i,j,k}^{n} \right|= \frac{1}{144}\left[r_x^2 r_y^2 \delta_x^2 \delta_y^2 +r_x^2 r_z^2 \delta_x^2 \delta_z^2+r_y^2 r_z^2 \delta_y^2 \delta_z^2\right]u_{i,j,k}^{n} \notag \\ & +\frac{1}{1728}\left[r_x^2 r_y^2 +r_x^2 r_z^2 +r_y^2 r_z^2 -r_x^2 r_y^2 r_z^2 \right] \delta_x^2 \delta_y^2 \delta_z^2 u_{i,j,k}^{n} \approx O(h^4) \end{align}

ADI Algorithm

\begin{align} &\left(1+\rho_x\delta_x^2\right)\left(1+\rho_y\delta_y^2\right)\left(1+\rho_z\delta_z^2\right)\delta_t^2 u_{i,j,k}^{n}= \left[r_x^2 \left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_x^2 \right. \notag \\ & \left. +r_y^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_y^2 +r_z^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\delta_z^2 \right]u_{i,j,k}^{n} \notag \\ & +\tau^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\left(1+\frac{1}{12}\delta_t^2\right)f_{i,j,k}^{n}. \qquad \qquad (5) \end{align}

It follows from (5) that, \begin{align} \left(1+\rho_x\delta_x^2\right)\left(1+\rho_y\delta_y^2\right)\left(1+\rho_z\delta_z^2\right) \delta_t^2 u_{i,j,k}^{n}=\mathbf{F}_{i,j,k}^{n} \end{align} where, \begin{align} \mathbf{F}_{i,j,k}^{n} &=+\left[r_x^2 \left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_x^2 +r_y^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\delta_y^2 \right. \notag \\ & \left. +r_z^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\delta_z^2 \right]u_{i,j,k}^{n} \notag \\ & +\tau^2\left(1+\frac{1}{12}\delta_x^2\right)\left(1+\frac{1}{12}\delta_y^2\right)\left(1+\frac{1}{12}\delta_z^2\right)\left(1+\frac{1}{12}\delta_t^2\right)f_{i,j,k}^{n}. \end{align}

The three-dimensional problem can be efficiently solved in three steps using ADI with the following steps,

The equation (6) also can be written as follows, \begin{align} \rho_x(\mathbf{X}_{i+1,j,k}+\mathbf{X}_{i-1,j,k})+(1-2\rho_x)\mathbf{X}_{i,j,k}=\mathbf{F}_{i,j,k}^{n}, \end{align}

Matrix form: $$ A_{j,k}=\begin{bmatrix} (1-2\rho_x) & \rho_x & 0 & \dots & \dots &0 \\ \rho_x & (1-2\rho_x) & \rho_x & 0 & \vdots &\vdots \\ 0 & \rho_x & (1-2\rho_x) & \rho_x &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & \rho_x & (1-2\rho_x) & \rho_x & 0 \\ \vdots &\vdots & 0 & \rho_x & (1-2\rho_x) & \rho_x \\ 0 &0 & \dots & 0 & \rho_x & (1-2\rho_x) \end{bmatrix}, $$ $$ \mathbf{X}_{j,k}=\begin{bmatrix} \mathbf{X}_{1,j,k} \\ \mathbf{X}_{2,j,k}\\ \vdots\\ \mathbf{X}_{N_x-2,j,k}\\ \mathbf{X}_{N_x-1,j,k} \end{bmatrix},~ \mathbf{b}_j^n= \begin{bmatrix} -\rho_x \mathbf{X}_{0,j,k}\\ 0\\ \vdots\\ 0\\ -\rho_x \mathbf{X}_{N_x,j,k}, \end{bmatrix} \text{ and } \mathbf{F}_{j}^{n}= \begin{bmatrix} \mathbf{F}_{1,j,k}^{n}\\ \mathbf{F}_{2,j,k}^{n}\\ \vdots\\ \mathbf{F}_{N_x-2,j,k}^{n}\\ \mathbf{F}_{N_x-1,j,k}^{n} \end{bmatrix}, $$ where \begin{align} \mathbf{X}_{0,j,k}&= \left(1+\rho_y\delta_y^2\right)\left(1+\rho_z\delta_z^2\right)\delta_t^2 f_{a,j,k}^{n}, \\ \mathbf{X}_{N_x,j,k}&= \left(1+\rho_y\delta_y^2\right)\left(1+\rho_z\delta_z^2\right)\delta_t^2 f_{b,j,k}^{n}. \end{align}

Therefore \begin{equation*} \mathbf{X}_{j,k}=A_j^{-1}\left(\mathbf{F}_{j,k}^{n}+\mathbf{b}_j^n\right). \end{equation*}

Next, we need to find $\mathbf{Y}_{i,k}$ for $1\leq i \leq N_x-1$ and $1\leq k\leq N_z-1$. We have \begin{align} \rho_y\left(\mathbf{Y}_{i,j-1,k}+\mathbf{Y}_{i,j+1,k}\right)+(1-2\rho_y)\mathbf{Y}_{i,j,k}=\mathbf{X}_{i,j,k}, \end{align}

In addition, $$ A_{i,k}=\begin{bmatrix} (1-2\rho_y) & \rho_y & 0 & \dots & \dots &0 \\ \rho_y & (1-2\rho_y) & \rho_y & 0 & \vdots &\vdots \\ 0 & \rho_y & (1-2\rho_y) & \rho_y &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & \rho_y & (1-2\rho_y) & \rho_y & 0 \\ \vdots &\vdots & 0 & \rho_y & (1-2\rho_y) & \rho_y \\ 0 &0 & \dots & 0 & \rho_y & (1-2\rho_y) \end{bmatrix}, $$ and \begin{equation*} \mathbf{b}_{i,k}^{n}= \begin{bmatrix} -\rho_y \mathbf{Y}_{i,0,k}\\ 0\\ \vdots\\ 0\\ -\rho_y \mathbf{Y}_{i,N_y,k} \end{bmatrix},~ \mathbf{X}_{i,k}=\begin{bmatrix} \mathbf{X}_{i,1,k} \\ \mathbf{X}_{i,2,k}\\ \vdots\\ \mathbf{X}_{i,N_y-2,k}\\ \mathbf{X}_{i,N_y-1,k} \end{bmatrix} \text{ and } \mathbf{Y}_{i,k}=\begin{bmatrix} \mathbf{Y}_{i,1,k} \\ \mathbf{Y}_{i,2,k}\\ \vdots\\ \mathbf{Y}_{i,N_y-2,k}\\ \mathbf{Y}_{i,N_y-1,k} \end{bmatrix}, \end{equation*} where \begin{align} \mathbf{Y}_{i,0,k}&=\left(1+\rho_z\delta_z^2\right)\delta_t^2f_{c,i,k}^{n},\\ \mathbf{Y}_{i,N_y,k}&=\left(1+\rho_z\delta_z^2\right)\delta_t^2 f_{d,i,k}^{n}. \end{align}

Hence, \begin{equation*} \label{eq1C6.15-1}\mathbf{Y}_{i,k}=A_{i,k}^{-1}\left(\mathbf{X}_{i,k}+\mathbf{b}_{i,k}^n\right), \end{equation*}

Finally, solving the following linear system is needed to find $\mathbf{Z}_{i,j,k}$

\begin{align} \begin{cases} \rho_z( \mathbf{Z}_{i,j,k+1}+\mathbf{Z}_{i,j,k-1})+(1-2\rho_z)\mathbf{Z}_{i,j,k}=\mathbf{Y}_{i,j,k},\\ u_{i,j,k+1}^{n+1}=\mathbf{Z}_{i,j,k}+2u_{i,j,k}^{n}-u_{i,j,k}^{n-1} \end{cases} \end{align}

As a result, \begin{equation*} \mathbf{u}_{i,j}^{n+1}=A_{i,j}^{-1}\left(\mathbf{Y}_{i,j}+\mathbf{b}_{i,j}^n\right), \end{equation*} where \begin{equation*} A_{i,j}=\begin{bmatrix} (1-2\rho_z) & \rho_z & 0 & \dots & \dots &0 \\ \rho_z & (1-2\rho_z) & \rho_z & 0 & \vdots &\vdots \\ 0 & \rho_z & (1-2\rho_z) & \rho_z &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & \rho_z & (1-2\rho_z) & \rho_z & 0 \\ \vdots &\vdots & 0 & \rho_z & (1-2\rho_z) & \rho_z \\ 0 &0 & \dots & 0 & \rho_z & (1-2\rho_z) \end{bmatrix}, \end{equation*} and \begin{equation*} \mathbf{b}_{i,j}^{n}= \begin{bmatrix} -\rho_z \delta_t^2 f_{e,i,j}^{n}\\ 0\\ \vdots\\ 0\\ -\rho_z \delta_t^2 f_{f,i,j}^{n} \end{bmatrix} \text{ and } \mathbf{u}_{i,j}^{n+1}=\begin{bmatrix} u_{i,j,1}^{n+1} \\ u_{i,j,2}^{n+1}\\ \vdots\\ u_{i,j,N_z-2}^{n+1}\\ u_{i,j,N_z-1}^{n+1} \end{bmatrix}. \end{equation*}

For stability, convergence and dispersion analyses, please see [1,2].

Example

Consider the following problem, \begin{align} \begin{cases} \dfrac{\partial^2 u}{\partial t^2} = \frac{1}{4} \Delta u +\frac{1}{4}e^{t+x+y+z},&(x,y,z)\in I,~t\in(0,1],\\ u(x,y,z,0)=g(x,y,z),&(x,y,z)\in J,\\ \dfrac{\partial }{\partial t}u(x,y,z,0) = s(x,y,z),&(x,y,z)\in J,\\ u(0,y,z,t) = e^{t+y+z},&(y,z)\in [0,1]\times[0,1],~t \in[0,1],\\ u(1,y,z,t) = e^{t+y+z+1},&(y,z)\in [0,1]\times[0,1],~t \in[0,1],\\ u(x,0,z,t) = e^{t+x+z},&(x,z)\in [0,1]\times[0,1],~t \in[0,1],\\ u(x,1,z,t) = e^{t+x+z+1},&(x,z)\in [0,1]\times[0,1],~t \in[0,1],\\ u(x,y,0,t) = e^{t+x+y},&(x,y)\in [0,1]\times[0,1],~t \in[0,1],\\ u(x,y,1,t) = e^{t+x+y+1},&(x,y)\in [0,1]\times[0,1],~t \in[0,1]. \end{cases} \end{align}

The exact solution corresponding to the problem can be found as follows, $$ u(x,t)=e^{t+x+y+z} $$

The Fourth-order Scheme

We can consider one of the following algorithms for this scheme.

The results of analyzing the numerical solution of the problem, using the explicit scheme, is available at the following tables.

What stands out from the tables is that the order of the method is 4 (in space and time). Moreover, consider the following experiment to test the stability condition.

References

  1. Das, Sambit, Wenyuan Liao, and Anirudh Gupta. "An efficient fourth-order low dispersive finite difference scheme for a 2-D acoustic wave equation." Journal of computational and Applied Mathematics 258 (2014): 151-167.

  2. Liao, Wenyuan. "On the dispersion, stability and accuracy of a compact higher-order finite difference scheme for 3D acoustic wave equation." Journal of Computational and Applied Mathematics 270 (2014): 571-583.

  3. Liao, Wenyuan. "An implicit fourth-order compact finite difference scheme for one-dimensional Burgers’ equation." Applied Mathematics and Computation 206.2 (2008): 755-764.