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

A one-dimensional form of the wave equation presented can be found as follows, \begin{align*} \begin{cases} \dfrac{\partial^2 u}{\partial t^2} = c^2(x) \dfrac{\partial^2 u}{\partial x^2}+f(x,t),&(x,t)\in (a,b)\times(T_1,T_2],\\ u(x,T_1)=g(x),&x\in [a,b],\\ \dfrac{\partial }{\partial t}u(x,T_1) = h(x),&x\in [a,b],\\ u(a,t) = f_a(t),&t \in[T_1,T_2],\\ u(b,t) = f_b(t),&t \in[T_1,T_2]. \end{cases} \end{align*}

The following notations are introduced to generate an algorithm for the problem,

Moreover, a general form of implicit Finite-Difference scheme of the problem (1) can be written as follows, \begin{align} \frac{1}{\tau^2}\delta_t^2 u_{i}^{n}= \frac{1}{h^2}c_{i}^2\left(\alpha \delta_x^2 u_{i}^{n+1} +\beta \delta_x^2 u_{i}^{n} +\gamma \delta_x^2 u_{i}^{n-1}\right)+f_{i}^{n}, \end{align} where $\delta^2$ is the central difference operator and $\alpha,~\beta$ and $\gamma$ are constants such that $\alpha+\beta+\gamma=1$. It follows that \begin{align} u_{i}^{n+1}-2u_{i}^{n}+u_{i}^{n-1} &= \alpha c_{i}^2 \lambda^2\left( u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1}\right) +\beta c_{i}^2 \lambda^2 \left( u_{i+1}^{n} -2 u_{i}^{n}+ u_{i-1}^{n}\right) +\gamma c_{i}^2 \lambda^2 \left(u_{i+1}^{n-1}-2 u_{i}^{n-1}+u_{i-1}^{n-1}\right)+ \tau^2f_{i}^{n}. \end{align}

This also can be simplified as follows, \begin{align}\label{eq.1} -\alpha c_{i}^2 \lambda^2 \left(u_{i-1}^{n+1}+u_{i+1}^{n+1}\right)+\left(1 +2 \alpha c_{i}^2 \lambda^2 \right)u_{i}^{n+1}&=\beta c_{i}^2 \lambda^2 \left(u_{i-1}^{n}+u_{i+1}^{n} \right) +2\left(1 -\beta c_{i}^2 \lambda^2 \right)u_{i}^{n}+\gamma c_{i}^2 \lambda^2 \left(u_{i-1}^{n-1}+u_{i+1}^{n-1} \right) \notag\\ & -\left(1 +2 \gamma c_{i}^2 \lambda^2 \right)u_{i}^{n-1}+\tau^2f_{i}^{n}. \end{align}

The following matrices are introduced to write down a matrix form of the equation (2), \begin{align} A=\begin{bmatrix} 1 +2 \alpha c_{1}^2 \lambda^2 & -\alpha c_{1}^2 \lambda^2 & 0 & \dots & \dots &0 \\ -\alpha c_{2}^2 \lambda^2 &1 +2 \alpha c_{2}^2 \lambda^2 & -\alpha c_{2}^2 \lambda^2 & 0 & \vdots &\vdots \\ 0 &-\alpha c_{3}^2 \lambda^2 &1 +2 \alpha c_{3}^2 \lambda^2 & -\alpha c_{3}^2 \lambda^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & -\alpha c_{M-3}^2 \lambda^2 & 1 +2 \alpha c_{M-3}^2 \lambda^2 &-\alpha c_{M-3}^2 \lambda^2 & 0 \\ \vdots &\vdots & 0 & -\alpha c_{N_x-2}^2 \lambda^2 & 1 +2 \alpha c_{N_x-2}^2 \lambda^2 &-\alpha c_{N_x-2}^2 \lambda^2 \\ 0 &0 & \dots & 0 & -\alpha c_{N_x-1}^2 \lambda^2 & 1 +2 \alpha c_{N_x-1}^2 \lambda^2 \end{bmatrix}, \end{align} \begin{align} B=\begin{bmatrix} 2 -2 \beta c_{1}^2 \lambda^2 & \beta c_{1}^2 \lambda^2 & 0 & \dots & \dots &0 \\ \beta c_{2}^2 \lambda^2 &2 -2 \beta c_{2}^2 \lambda^2 & \beta c_{2}^2 \lambda^2 & 0 & \vdots &\vdots \\ 0 &\beta c_{3}^2 \lambda^2 &2 -2 \beta c_{3}^2 \lambda^2 & \beta c_{3}^2 \lambda^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & \beta c_{M-3}^2 \lambda^2 & 2 -2 \beta c_{M-3}^2 \lambda^2 &\beta c_{M-3}^2 \lambda^2 & 0 \\ \vdots &\vdots & 0 & \beta c_{N_x-2}^2 \lambda^2 & 2 -2 \beta c_{N_x-2}^2 \lambda^2 &\beta c_{N_x-2}^2 \lambda^2 \\ 0 &0 & \dots & 0 & \beta c_{N_x-1}^2 \lambda^2 & 2 -2 \beta c_{N_x-1}^2 \lambda^2 \end{bmatrix}, \end{align} \begin{align} C=\begin{bmatrix} -1-2 \gamma c_{1}^2 \lambda^2 & \gamma c_{1}^2 \lambda^2 & 0 & \dots & \dots &0 \\ \gamma c_{2}^2 \lambda^2 &-1-2 \gamma c_{2}^2 \lambda^2 & \gamma c_{2}^2 \lambda^2 & 0 & \vdots &\vdots \\ 0 &\gamma c_{3}^2 \lambda^2 &-1-2 \gamma c_{3}^2 \lambda^2 & \gamma c_{3}^2 \lambda^2 &0 &\vdots \\ 0 &0 & \ddots & \ddots &\ddots &\vdots \\ \vdots &\vdots & \gamma c_{M-3}^2 \lambda^2 & -1-2 \gamma c_{M-3}^2 \lambda^2 &\gamma c_{M-3}^2 \lambda^2 & 0 \\ \vdots &\vdots & 0 & \gamma c_{N_x-2}^2 \lambda^2 & -1-2 \gamma c_{N_x-2}^2 \lambda^2 &\gamma c_{N_x-2}^2 \lambda^2 \\ 0 &0 & \dots & 0 & \gamma c_{N_x-1}^2 \lambda^2 & -1-2 \gamma c_{N_x-1}^2 \lambda^2 \end{bmatrix}, \end{align} and also \begin{align} \mathbf{u}^n=\begin{bmatrix} u_{1}^{n} \\ u_{2}^{n}\\ \vdots\\ u_{N_x-2}^{n}\\ u_{N_x-1}^{n} \end{bmatrix},~ \mathbf{b}^n= \begin{bmatrix} c_{1}^2 \lambda^2 \left(\alpha u_{0}^{n+1}+\beta u_{0}^{n}+\gamma u_{0}^{n-1}\right) \\ 0\\ \vdots\\ 0\\ c_{N_x-1}^2 \lambda^2 \left(\alpha u_{N_x}^{n+1}+\beta u_{N_x}^{n}+\gamma u_{N_x}^{n-1}\right) \end{bmatrix},~ \mathbf{F}^n=\tau^2 \begin{bmatrix} f_{1}^{n}\\ f_{2}^{n}\\ \vdots\\ f_{N_x-2}^{n}\\ f_{N_x-1}^{n} \end{bmatrix}. \end{align} Therefore, the following matrix form of equation (4) can be stated based as follows, \begin{align} A \mathbf{u}^{n+1}=B \mathbf{u}^{n}+C \mathbf{u}^{n-1}+\mathbf{b}^n+\mathbf{F}^n, \end{align} which implies \begin{align} \mathbf{u}^{n+1}=A^{-1}\left(B \mathbf{u}^{n}+C \mathbf{u}^{n-1}+\mathbf{b}^n+ \mathbf{F}^n\right). \end{align}

Finite differnce method for the wave equation in time usually require computing the solution on the first two time steps. The solution on the second time step, $u(x,T_1+k)$, can be approximated as follows by means of Taylor series, \begin{align} u(x,T_1+\tau)&=u(x,T_1)+u_t(x,T_1)\tau+\frac{1}{2} u_{tt}(x,T_1)\tau^2+\frac{1}{6} u_{ttt}(x,T_1)\tau^3 +\frac{1}{24} u_{tttt}(x,T_1)\tau^4 +\mathcal{O}(\tau^5) \notag \\ & =g(x)+s(x)\tau+\frac{1}{2} u_{tt}(x,T_1)\tau^2+\frac{1}{6} u_{ttt}(x,T_1)\tau^3+\frac{1}{24} u_{tttt}(x,T_1)\tau^4+\mathcal{O}(\tau^5) \notag \\ & =g(x)+s(x)\tau+\frac{1}{2}[c^2(x) \Delta g(x)+f(x,T_1)]\tau^2 +\frac{1}{6} [c^2(x) \Delta s(x)+f_t(x,T_1)]\tau^3 \notag\\ & +\frac{1}{24}[c^2(x) \Delta [c^2(x) \Delta g(x)+f(x,T_1)] +f_{tt}(x,T_1)]\tau^4 +\mathcal{O}(\tau^5). \end{align}

Explicit Scheme

First, consider $\alpha=\gamma=0$ and $\beta=1$ in the equation (2) to get the explicit method.

Implicit Scheme

Similarly, consider $\alpha=\gamma=\dfrac{1}{4}$ and $\beta=\dfrac{1}{2}$ in the equation (2) to the obtain implicit method.

For Stability and convergence analyses, see this file.

Example

Consider the following problem, \begin{align} \begin{cases} \frac{\partial^2 u}{\partial t^2} &= \mathrm{e}^{-x} \frac{\partial^2 u}{\partial x^2}+{\mathrm{e}}^{t+x}-{\mathrm{e}}^t ,&(x,t)\in (0,1)\times(0,1],\\ u(x,0)&=\mathrm{e}^{x},&x\in [0,1],\\ \frac{\partial }{\partial t}u(x,0) &= \mathrm{e}^{x},&x\in [0,1],\\ u(0,t) &= \mathrm{e}^{t} &t \in[0,1],\\ u(1,t) &= \mathrm{e}^{t+1},&t \in[0,1]. \end{cases} \end{align} The exact solution corresponding to the problem can be found as follows, $$ u(x,t)=\mathrm{e}^{x+t} $$

The Explicit Scheme

What stands out from the tables is that the order of the method is 2 (space and time). To test the stability condition, we can try the following experiment.

The table highlights that the stability criteria were violated.

The Implicit scheme

Convergence Rate

For a fixed $N_t$, we have

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.