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.

Second-Order Finite-Difference Methods 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} \end{align}

Consider the following notations:

Explicit Scheme

An explicit finite difference scheme for solving the problem (1) can be considered as follows, \begin{align} \frac{1}{\tau^2} \delta_t^2 u_{i,j,k}^{n}&= c_{i,j,k}^2\left(\frac{1}{h_x^2} \delta_x^2+\frac{1}{h_y^2} \delta_y^2+\frac{1}{h_z^2} \delta_z^2\right)u_{i,j,k}^{n}+f_{i,j,k}^n. \end{align} To develop a numerical algorithm, we can rewrite this as follows, \begin{align} u_{i,j,k}^{n+1}&=c_{i,j,k}^2\left[\lambda_x^2\left(u_{i+1,j,k}^{n}+u_{i-1,j,k}^{n}\right)+\lambda_y^2\left(u_{i,j+1,k}^{n}+u_{i,j-1,k}^{n}\right)+\lambda_z^2\left(u_{i,j,k+1}^{n}+u_{i,j,k-1}^{n}\right)\right] \notag \\ & +2\left(1-c_{i,j,k}^2\left(\lambda_x^2+\lambda_y^2+\lambda_z^2\right)\right)u_{i,j,k}^{n}-u_{i,j,k}^{n-1}+\tau^2f_{i,j,k}^n, \end{align}

Implicit Scheme

A general implicit finite difference form of the problem can be expressed as follows, \begin{align} \frac{1}{\tau^2} \delta_t^2 u_{i,j,k}^{n}&= c_{i,j,k}^2 \left(\frac{1}{h_x^2} \delta_x^2+\frac{1}{h_y^2} \delta_y^2 +\frac{1}{h_z^2} \delta_z^2 \right) \left[\alpha u_{i,j,k}^{n+1}+\beta u_{i,j,k}^{n} +\gamma u_{i,j,k}^{n-1}\right]+f_{i,j,k}^n, \end{align} where $\alpha,~\beta$ and $\gamma$ are constants such that $\alpha+\beta+\gamma=1$.

The equation (2) can be expanded as follows, \begin{align} \left(1-\alpha c^2_{i,j,k}\left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 +\lambda_z^2 \delta_z^2 \right)\right)u_{i,j,k}^{n+1}&= \left[c_{i,j,k}^2 \beta \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 +\lambda_z^2 \delta_z^2 \right)+2\right]u_{i,j,k}^{n} \notag \\ & +\left[c_{i,j,k}^2\gamma \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 +\lambda_z^2 \delta_z^2 \right)-1\right]u_{i,j,k}^{n-1} +\tau^2f_{i,j,k}^n. \end{align} Note that, \begin{align} &\left(1-\alpha c^2_{i,j,k}\left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 +\lambda_z^2 \delta_z^2 \right)\right)u_{i,j,k}^{n} -\left(1-\alpha c^2_{i,j,k} \lambda_x^2 \delta_x^2\right)\left(1-\alpha c^2_{i,j,k} \lambda_y^2 \delta_y^2\right)\left(1-\alpha c^2_{i,j,k} \lambda_z^2 \delta_z^2\right)u_{i,j,k}^{n}= \notag \\ & \alpha^2 c^4_{i,j,k} \left(\lambda_x^2 \lambda_y^2 \delta_x^2\delta_y^2+\lambda_x^2 \lambda_z^2 \delta_x^2\delta_z^2 +\lambda_y^2 \lambda_z^2 \delta_y^2\delta_z^2\right)u_{i,j,k}^{n} -\alpha^3 c^6_{i,j,k} \lambda_x^2 \lambda_y^2 \lambda_z^2 \delta_x^2\delta_y^2\delta_z^2 u_{i,j,k}^{n}\approx O(h^4), \qquad \qquad (3). \end{align}

Furthermore, for simplicity, consider \begin{align} \mathbf{F}_{i,j,k}^{n}&= \left[c_{i,j,k}^2 \beta \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 +\lambda_z^2 \delta_z^2 \right)+2\right]u_{i,j,k}^{n} +\left[c_{i,j,k}^2\gamma \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 +\lambda_z^2 \delta_z^2 \right)-1\right]u_{i,j,k}^{n-1} +\tau^2f_{i,j,k}^n. \end{align}

As a result, \begin{align} \left(1-\alpha c^2_{i,j,k} \lambda_x^2 \delta_x^2\right)\left(1-\alpha c^2_{i,j,k} \lambda_y^2 \delta_y^2\right)\left(1-\alpha c^2_{i,j,k} \lambda_z^2 \delta_z^2\right)u_{i,j,k}^{n+1}=\mathbf{F}_{i,j,k}^{n}. \end{align}

ADI Algorithm

The first time step after the initial time level is needed for proceeding the introduced algorithm and $u(x,y,z,T_1+\tau)$ is approximated by means of Taylor series, \begin{align} u(x,y,z,T_1+\tau)&= g(x,y,z)+s(x,y,z)\tau+\frac{1}{2}[c^2(x,y,z) \Delta g(x,y,z)+f(x,y,z,T_1)]\tau^2 \notag\\ & +\frac{1}{6} [c^2(x,y,z) \Delta s(x,y,z)+f_t(x,y,z,T_1)]\tau^3 \notag\\ & +\frac{1}{24}[c^2(x,y,z) \Delta [c^2(x,y,z) \Delta g(x,y,z)+f(x,y,z,T_1)] +f_{tt}(x,y,z,T_1)]\tau^4 +\mathcal{O}(\tau^5) \end{align}

To numerical algorithm based,, first, for current time step, fix $n$ and solve the following system for $i=1,2,\ldots,N_x-1$, $j=1,2,\ldots,N_y-1$ and $k=1,2,\ldots,N_z-1$, \begin{align} \left(1-\alpha c^2_{i,j,k} \lambda_x^2 \delta_x^2\right)\mathbf{X}_{i,j,k}=\mathbf{F}_{i,j,k}^{n}, \qquad \qquad (4) \end{align} In other words, \begin{align}\label{eq1C1.11} -\alpha c_{i,j,k}^2 \lambda_x^2 \left(\mathbf{X}_{i-1,j,k}+\mathbf{X}_{i+1,j,k}\right)+\left(1 +2 \alpha c_{i,j,k}^2 \lambda_x^2\right) \mathbf{X}_{i,j,k} =\mathbf{F}_{i,j,k}^{n}. \end{align} In matrix form: $$ A_{j,k}=\begin{bmatrix} 1 +2 \alpha c_{1,j,k}^2\lambda_x^2 & -\alpha c_{1,j,k}^2\lambda_x^2 & 0 & \dots & \dots &0 \\ -\alpha c_{2,j,k}^2\lambda_x^2 &1 +2 \alpha c_{2,j,k}^2\lambda_x^2 & -\alpha c_{2,j,k}^2\lambda_x^2 & 0 & \vdots &\vdots \\ 0 &-\alpha c_{3,j,k}^2\lambda_x^2 &1 +2 \alpha c_{3,j,k}^2\lambda_x^2 & -\alpha c_{3,j,k}^2\lambda_x^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & -\alpha c_{N_x-3,j,k}^2\lambda_x^2 & 1 +2 \alpha c_{N_x-3,j,k}^2\lambda_x^2 &-\alpha c_{N_x-3,j,k}^2\lambda_x^2 & 0 \\ \vdots &\vdots & 0 & -\alpha c_{N_x-2,j,k}^2\lambda_x^2 & 1 +2 \alpha c_{N_x-2,j,k}^2\lambda_x^2 &-\alpha c^2_{N_x-2,j,k}\lambda_x^2 \\ 0 &0 & \dots & 0 & -\alpha c_{N_x-1,j,k}^2\lambda_x^2 & 1 +2 \alpha c_{N_x-1,j,k}^2\lambda_x^2 \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,k}^{n}= \begin{bmatrix} \alpha c^2_{1,j,k} \lambda_x^2 \mathbf{X}_{0,j,k}\\ 0\\ \vdots\\ 0\\ \alpha c^2_{N_x-1,j,k}\lambda_x^2 \mathbf{X}_{N_x,j,k}, \end{bmatrix}, \mathbf{F}_{j,k}^{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 $\mathbf{X}_{0,j,k}$ and $\mathbf{X}_{N_x,j,k}$ for $1\leq j\leq N_x-1$ and $1\leq j\leq N_y-1$ can be stated as follows, \begin{align} \label{eq1C1.12} \mathbf{X}_{0,j,k}&= \left(1-\alpha c^2_{i,j,k} \lambda_y^2 \delta_y^2\right)\left(1-\alpha c^2_{i,j,k} \lambda_z^2 \delta_z^2\right)f_{a,j,k}^{n+1}, \\ \label{eq1C1.13} \mathbf{X}_{N_x,j,k}&=\left(1-\alpha c^2_{i,j,k} \lambda_y^2 \delta_y^2\right)\left(1-\alpha c^2_{i,j,k} \lambda_z^2 \delta_z^2\right)f_{a,j,k}^{n+1}. \end{align} and, $f_{a,j,k}^{n}\approx f_{a}(jh_y,kh_z,n\tau)$ and $f_{b,j,k}^{n}\approx f_{b}(jh_y,kh_z,n\tau)$.

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

Next, we need to solve the following system, \begin{align} \left(1-\alpha c^2_{i,j,k} \lambda_y^2 \delta_y^2\right)\mathbf{Y}_{i,j,k}=\mathbf{X}_{i,j,k}. \end{align} This can be expanded as follows, \begin{align} -\alpha c_{i,j,k}^2 \lambda_y^2 \left(\mathbf{Y}_{i,j-1,k}+\mathbf{Y}_{i,j+1,k}\right)+\left(1 +2 \alpha c_{i,j,k}^2\lambda_y^2\right)\mathbf{Y}_{i,j,k}=\mathbf{X}_{i,j,k}. \end{align} Hence, \begin{equation*} \label{eq1C2.13}\mathbf{Y}_{i,k}=A_{i,k}^{-1}\left(\mathbf{X}_{i,k}+\mathbf{b}_{i,k}^n\right), \end{equation*} where \begin{equation*}\label{eq1C2.14} A_{i,k}=\begin{bmatrix} 1 +2 \alpha c_{i,1,k}^2\lambda_y^2 & -\alpha c_{i,1,k}^2\lambda_y^2 & 0 & \dots & \dots &0 \\ -\alpha c_{i,2,k}^2\lambda_y^2 &1 +2 \alpha c_{i,2,k}^2\lambda_y^2 & -\alpha c_{i,2,k}^2\lambda_y^2 & 0 & \vdots &\vdots \\ 0 &-\alpha c_{i,3,k}^2\lambda_y^2 &1 +2 \alpha c_{i,3,k}^2\lambda_y^2 & -\alpha c_{i,3,k}^2\lambda_y^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & -\alpha c_{i,N_y-3,k}^2\lambda_y^2 & 1 +2 \alpha c_{i,N_y-3,k}^2\lambda_y^2 &-\alpha c_{i,N_y-3,k}^2\lambda_y^2 & 0 \\ \vdots &\vdots & 0 & -\alpha c_{i,N_y-2,k}^2\lambda_y^2 & 1 +2 \alpha c_{i,N_y-2,k}^2\lambda_y^2 &-\alpha c^2_{i,N_y-2,k}\lambda_y^2 \\ 0 &0 & \dots & 0 & -\alpha c_{i,N_y-1,k}^2\lambda_y^2 & 1 +2 \alpha c_{i,N_y-1,k}^2\lambda_y^2 \end{bmatrix}, \end{equation*} and \begin{equation*}\label{eq1C2.15} \mathbf{b}_{i,k}^{n}= \begin{bmatrix} \alpha c^2_{i,1,k}\lambda_y^2 \mathbf{Y}_{i,0,k}\\ 0\\ \vdots\\ 0\\ \alpha c^2_{i,N_y-1,k}\lambda_y^2 \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-1,k}\\ \mathbf{X}_{i,N_y,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*} Moreover, $\mathbf{Y}_{i,0,k}$ and $\mathbf{Y}_{i,N_y,k}$ can be approximated as follows, \begin{align} \label{eq1C1.17}\mathbf{Y}_{i,0,k}&=-\alpha c_{i,0,k}^2 \lambda_z^2 \left(f_{c,i,k-1}^{n+1}+f_{c,i,k+1}^{n+1}\right)+\left(1 +2 \alpha c_{i,0,k}^2\lambda_z^2\right)f_{c,i,k}^{n+1},\\ \label{eq1C1.18}\mathbf{Y}_{i,N_y,k}&=-\alpha c_{i,N_y,k}^2 \lambda_z^2 \left(f_{d,i,k-1}^{n+1}+f_{d,i,k+1}^{n+1}\right)+\left(1 +2 \alpha c_{i,N_y,k}^2\lambda_z^2\right)f_{d,i,k}^{n+1}, \end{align} where $f_{c,i,k}^{n}\approx f_{c}(ih_x,kh_z,n\tau)$ and $f_{d,i,k}^{n}\approx f_{d}(ih_x,kh_z,n\tau)$.

Finally, we need to solve the following system for $u_{i,j,k}^{n}$ for each $i=1,2,\ldots,N_x-1$, $j=1,2,\ldots,N_y-1$ and $k=1,2,\ldots,N_z-1$, \begin{align} \label{eq1C1.19}\left(1-\alpha c^2_{i,j,k} \lambda_z^2 \delta_z^2\right)u_{i,j,k}^{n+1}=\mathbf{Y}_{i,j,k}. \end{align} In other words, \begin{align} \label{eq1C1.20}-\alpha c_{i,j,k}^2 \lambda_z^2 \left(u_{i,j,k-1}^{n+1}+u_{i,j,k+1}^{n+1}\right)+\left(1 +2 \alpha c_{i,j,k}^2\lambda_z^2\right)u_{i,j,k}^{n+1}=\mathbf{Y}_{i,j,k}. \end{align} Hence, \begin{equation*} \label{eq1C1.21}u_{i,j}^{n+1}=A_{i,j}^{-1}\left(\mathbf{Y}_{i,j,k}+\mathbf{b}_{i,j}^n\right), \end{equation*} where \begin{equation*}\label{eq1C1.22} A_{i,j}=\begin{bmatrix} 1 +2 \alpha c_{i,j,1}^2\lambda_z^2 & -\alpha c_{i,j,1}^2\lambda_z^2 & 0 & \dots & \dots &0 \\ -\alpha c_{i,j,2}^2\lambda_z^2 &1 +2 \alpha c_{i,j,2}^2\lambda_z^2 & -\alpha c_{i,j,2}^2\lambda_z^2 & 0 & \vdots &\vdots \\ 0 &-\alpha c_{i,j,3}^2\lambda_z^2 &1 +2 \alpha c_{i,j,3}^2\lambda_z^2 & -\alpha c_{i,j,3}^2\lambda_z^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & -\alpha c_{i,j,N_z-3}^2\lambda_z^2 & 1 +2 \alpha c_{i,j,N_z-3}^2\lambda_z^2 &-\alpha c_{i,j,N_z-3}^2\lambda_z^2 & 0 \\ \vdots &\vdots & 0 & -\alpha c_{i,j,N_z-2}^2\lambda_z^2 & 1 +2 \alpha c_{i,j,N_z-2}^2\lambda_z^2 &-\alpha c^2_{i,j,N_z-2}\lambda_y^2 \\ 0 &0 & \dots & 0 & -\alpha c_{i,j,N_z-1}^2\lambda_z^2 & 1 +2 \alpha c_{i,j,N_z-1}^2\lambda_y^2 \end{bmatrix}, \end{equation*} and \begin{equation*} \mathbf{b}_{i,j}^{n}= \begin{bmatrix} \alpha c^2_{i,j,1}\lambda_z^2 f_{e,i,j}^{n+1}\\ 0\\ \vdots\\ 0\\ \alpha c^2_{i,j,N_z-1}\lambda_z^2 f_{f,i,j}^{n+1} \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*}

Note that $f_{e,i,j}^{n}\approx f_{e}(ih_x,jh_y,n\tau)$ and $f_{f,i,j}^{n}\approx f_{e}(ih_x,jh_y,n\tau)$.

For Stability and convergence analyses, see this file.

Example

Consider the following problem, \begin{align} \begin{cases} \dfrac{\partial^2 u}{\partial t^2} = \frac{1}{3}e^{x+y+z} \Delta u +e^{-t-x-y-z}-e^{-t},&(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,1) = 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 Explicit Scheme

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

The Implicit Scheme

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.