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 Two Dimensional Wave Equation in Time Domain

A two-dimensional form of the wave equation presented can be found as follows,

\begin{align} \begin{cases} \frac{\partial^2 u}{\partial t^2} = c^2(x,y) \Delta u+f(x,y,t),&(x,y)\in I,~t\in(T_1,T_2],\\ u(x,y,T_1)=g(x,y),&(x,y)\in J,\\ \frac{\partial }{\partial t}u(x,y,T_1) = s(x,y),&(x,y)\in J,\\ u(a,y,t) = f_{a}(y,t),&y\in [c,d],~t \in[T_1,T_2],\\ u(b,y,t) = f_{b}(y,t),&y\in [c,d],~t \in[T_1,T_2],\\ u(x,c,t) = f_{c}(x,t),&x\in [a,b],~t \in[T_1,T_2],\\ u(x,d,t) = f_{d}(x,t),&x\in [a,b],~t \in[T_1,T_2], \end{cases} \end{align}

where $I=(a,b)\times(c,d)$ and $J=[a,b]\times[c,d]$.

Consider the following notations:

Implicit Scheme

An explicit finite difference scheme for solving the two-dimensional problem can be considered as follows, \begin{align} \frac{1}{\tau^2} \delta_t^2 u_{i,j}^{n}= c_{i,j}^2\left(\frac{1}{h_x^2} \delta_x^2+\frac{1}{h_y^2} \delta_y^2\right)u_{i,j}^{n} +f_{i,j}^n, \end{align} which also can be expressed as the following form, \begin{align} u_{i,j}^{n+1}=c_{i,j}^2\left[\lambda_x^2\left(u_{i+1,j}^{n}+u_{i-1,j}^{n}\right)+\lambda_y^2\left(u_{i,j+1}^{n}+u_{i,j-1}^{n}\right) \right] +2\left(1-c_{i,j}^2\left(\lambda_x^2+\lambda_y^2\right)\right)u_{i,j}^{n}-u_{i,j}^{n-1}+\tau^2f_{i,j}^n, \end{align}

Implicit Scheme

Moreover, 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}^{n}= c_{i,j}^2 \left(\frac{1}{h_x^2} \delta_x^2+\frac{1}{h_y^2} \delta_y^2 \right) \left[\alpha u_{i,j}^{n+1}+\beta u_{i,j}^{n}+\gamma u_{i,j}^{n-1}\right]+f_{i,j}^n, \end{align} where $\alpha,~\beta$ and $\gamma$ are constants such that $\alpha+\beta+\gamma=1$. The equation (3) also can be written in the following form, \begin{align} \left(1-\alpha c^2_{i,j}\left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 \right)\right)u_{i,j}^{n+1}= \left[c_{i,j}^2 \beta \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 \right)+2\right]u_{i,j}^{n} +\left[c_{i,j}^2\gamma \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 \right)-1\right]u_{i,j}^{n-1} +\tau^2f_{i,j}^n. \end{align}

Note that, \begin{align} |\left(1-\alpha c^2_{i,j}\left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 \right)\right)u_{i,j}^{n} -\left(1-\alpha c^2_{i,j} \lambda_x^2 \delta_x^2\right)\left(1-\alpha c^2_{i,j} \lambda_y^2 \delta_y^2\right)u_{i,j}^{n}| =\alpha^2 c^4_{i,j} \lambda_x^2 \lambda_y^2 \delta_x^2\delta_y^2u_{i,j}^{n}, \end{align} is fourth-order and can be neglected.

Furthermore, for simplicity, consider the right hand side of (4) as $\mathbf{F}_{i,j}^{n}$, i.e. \begin{align} \mathbf{F}_{i,j}^{n}= \left[c_{i,j}^2 \beta \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 \right)+2\right]u_{i,j}^{n} +\left[c_{i,j}^2\gamma \left(\lambda_x^2 \delta_x^2 +\lambda_y^2 \delta_y^2 \right)-1\right]u_{i,j}^{n-1}+\tau^2f_{i,j}^n. \end{align} It follows from (2) and (5), \begin{align} \left(1-\alpha c^2_{i,j} \lambda_x^2 \delta_x^2\right)\left(1-\alpha c^2_{i,j} \lambda_y^2 \delta_y^2\right)u_{i,j}^{n+1}=\mathbf{F}_{i,j}^{n}. \end{align}

An algorithm based on ADI [1, 2] is going be developed to approximate the solution of equation (6). The two-dimensional problem can be efficiently solved in two steps using ADI with the following steps,

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

A numerical algorithm based on a second-order alternating direction implicit (ADI) method is going to be presented. For current time step, fix $n$ and solve the following system for $i=1,2,\ldots,N_x-1$ and $j=1,2,\ldots,N_y-1$, \begin{align} \left(1-\alpha c^2_{i,j} \lambda_x^2 \delta_x^2\right)\mathbf{X}_{i,j}=\mathbf{F}_{i,j}^{n}, \end{align}

The equation (9) can also be written in the following form, \begin{align} -\alpha c_{i,j}^2 \lambda_x^2 \left(\mathbf{X}_{i-1,j}+\mathbf{X}_{i+1,j}\right)+\left(1 +2 \alpha c_{i,j}^2 \lambda_x^2\right) \mathbf{X}_{i,j} =\mathbf{F}_{i,j}^{n}. \end{align}

Therefore, the following matrices are needed to generate an algorithm, \begin{align} A_{j}=\begin{bmatrix} 1 +2 \alpha c_{1,j}^2\lambda_x^2 & -\alpha c_{1,j}^2\lambda_x^2 & 0 & \dots & \dots &0 \\ -\alpha c_{2,j}^2\lambda_x^2 &1 +2 \alpha c_{2,j}^2\lambda_x^2 & -\alpha c_{2,j}^2\lambda_x^2 & 0 & \vdots &\vdots \\ 0 &-\alpha c_{3,j}^2\lambda_x^2 &1 +2 \alpha c_{3,j}^2\lambda_x^2 & -\alpha c_{3,j}^2\lambda_x^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & -\alpha c_{N_x-3,j}^2\lambda_x^2 & 1 +2 \alpha c_{N_x-3,j}^2\lambda_x^2 &-\alpha c_{N_x-3,j}^2\lambda_x^2 & 0 \\ \vdots &\vdots & 0 & -\alpha c_{N_x-2,j}^2\lambda_x^2 & 1 +2 \alpha c_{N_x-2,j}^2\lambda_x^2 &-\alpha c^2_{N_x-2,j}\lambda_x^2 \\ 0 &0 & \dots & 0 & -\alpha c_{N_x-1,j}^2\lambda_x^2 & 1 +2 \alpha c_{N_x-1,j}^2\lambda_x^2 \end{bmatrix}, \end{align} \begin{align} \mathbf{X}_{j}=\begin{bmatrix} \mathbf{X}_{1,j} \\ \mathbf{X}_{2,j}\\ \vdots\\ \mathbf{X}_{N_x-2,j}\\ \mathbf{X}_{N_x-1,j} \end{bmatrix},~ \mathbf{b}_{j}^{n}= \begin{bmatrix} \alpha c^2_{1,j} \lambda_x^2 \mathbf{X}_{0,j}\\ 0\\ \vdots\\ 0\\ \alpha c^2_{N_x-1,j}\lambda_x^2 \mathbf{X}_{N_x,j}, \end{bmatrix}, \mathbf{F}_{j}^{n}= \begin{bmatrix} \mathbf{F}_{1,j}^{n}\\ \mathbf{F}_{2,j}^{n}\\ \vdots\\ \mathbf{F}_{N_x-2,j}^{n}\\ \mathbf{F}_{N_x-1,j}^{n} \end{bmatrix}. \end{align} where $\mathbf{X}_{0,j}$ and $\mathbf{X}_{N_x,j}$ for $1\leq j\leq N_x-1$ can be approximated as follows, \begin{align}\label{eq1B1.16} \mathbf{X}_{0,j}&=\left(1-\alpha c^2_{i,j} \lambda_y^2 \delta_y^2\right) f_{a,j}^{n+1}, \\ \label{eq1B1.17} \mathbf{X}_{N_x,j}&=\left(1-\alpha c^2_{i,j} \lambda_y^2 \delta_y^2\right) f_{b,j}^{n+1}, \end{align} and $f_{a,j}^{n}\approx f_{a}(jh_y,kh_z)$ and $f_{b,j}^{n}\approx f_{b}(jh_y,kh_z)$.

Therefore, the following matrix form can be written as follows, \begin{align} \mathbf{X}_{j}=A_j^{-1}\left(\mathbf{F}_{j}^{n}+\mathbf{b}_j^n\right). \end{align} Having found $\mathbf{X}_{j}$ for each $j=1,2,\ldots,N_y-1$, it's time to solve the following system, \begin{align} \left(1-\alpha c^2_{i,j} \lambda_y^2 \delta_y^2\right)u_{i,j}^{n+1}=\mathbf{X}_{i,j}. \end{align} This also can be written in the following form, \begin{align} -\alpha c_{i,j}^2 \lambda_y^2 \left(u_{i,j-1}^{n+1}+u_{i,j+1}^{n+1}\right)+\left(1 +2 \alpha c_{i,j}^2\lambda_y^2\right)u_{i,j}^{n+1}=\mathbf{X}_{i,j}. \end{align}

The following matrix form can be written \begin{equation} \label{eq1B1.23}u_{i}^{n+1}=A_{i}^{-1}\left(\mathbf{X}_{i}+\mathbf{b}_{i}^n\right), \end{equation} where \begin{equation}\label{eq1B1.24} A_{i}=\begin{bmatrix} 1 +2 \alpha c_{i,1}^2\lambda_y^2 & -\alpha c_{i,1}^2\lambda_y^2 & 0 & \dots & \dots &0 \\ -\alpha c_{i,j,2}^2\lambda_y^2 &1 +2 \alpha c_{i,j,2}^2\lambda_y^2 & -\alpha c_{i,j,2}^2\lambda_y^2 & 0 & \vdots &\vdots \\ 0 &-\alpha c_{i,j,3}^2\lambda_y^2 &1 +2 \alpha c_{i,j,3}^2\lambda_y^2 & -\alpha c_{i,j,3}^2\lambda_y^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & -\alpha c_{i,j,N_y-3}^2\lambda_y^2 & 1 +2 \alpha c_{i,j,N_y-3}^2\lambda_y^2 &-\alpha c_{i,j,N_y-3}^2\lambda_y^2 & 0 \\ \vdots &\vdots & 0 & -\alpha c_{i,j,N_y-2}^2\lambda_y^2 & 1 +2 \alpha c_{i,j,N_y-2}^2\lambda_y^2 &-\alpha c^2_{i,j,N_y-2}\lambda_y^2 \\ 0 &0 & \dots & 0 & -\alpha c_{i,N_y-1}^2\lambda_y^2 & 1 +2 \alpha c_{i,N_y-1}^2\lambda_y^2 \end{bmatrix}, \end{equation} and also \begin{equation}\label{eq1B1.25} \mathbf{b}_{i}^{n}= \begin{bmatrix} \alpha c^2_{i,1}\lambda_y^2 f_{e,i,j}^{n+1}\\ 0\\ \vdots\\ 0\\ \alpha c^2_{i,N_y-1}\lambda_y^2 f_{f,i,j}^{n+1} \end{bmatrix} \text{ and } \mathbf{u}_{i}^{n+1}=\begin{bmatrix} u_{i,1}^{n+1} \\ u_{i,j,2}^{n+1}\\ \vdots\\ u_{i,j,N_y-2}^{n+1}\\ u_{i,N_y-1}^{n+1} \end{bmatrix}, \end{equation} where $f_{c,i}^{n}\approx f_{e}(ih_x,n\tau)$ and $f_{d,i}^{n}\approx f_{d}(ih_x,n\tau)$.

For Stability and convergence analyses, see this file.

Example

Consider the following problem, \begin{align*} \begin{cases} \frac{\partial^2 u}{\partial t^2} = \frac{1}{2}e^{x+y} \Delta u+e^{-t}\left(e^{-x-y}-1\right) ,&(x,y)\in I,~t\in(0,0],\\ u(x,y,0)=e^{-x-y},&(x,y)\in J,\\ \frac{\partial }{\partial t}u(x,y,T_1) = -e^{-x-y},&(x,y)\in J,\\ u(0,y,t) = e^{-t-y},&y\in [0,1],~t \in[0,1],\\ u(1,y,t) = e^{-t-y-1},&y\in [0,1],~t \in[0,1],\\ u(x,0,t) = e^{-t-x},&x\in [0,1],~t \in[0,1],\\ u(x,1,t) = e^{-t-x-1},&x\in [0,1],~t \in[0,1], \end{cases} \end{align*} where $I=(0,1)\times(0,1)$ and $J=[0,1]\times[0,1]$.

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

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.