Convergence Analysis of Numerical Solutions of Advection-Diffusion-Reaction Equations Using a Finite Difference Method an Application to Air Pollution Problems

1. Introduction

Advection-diffusion-reaction equations (ADRE) have been widely used as mathematical models in many areas of science and engineering, for example, in areas such as water pollution, air pollution, molecular diffusion, and chemical engineering. There are now many numerical methods that have been developed to solve linear and nonlinear ADREs. Some examples are as follows. In 2012, Diego et al [1] developed a model consisting of a decoupled system of advection-diffusion-reaction equations, along with the Navier-Stokes equations of incompressible flow, and solved the model by using the finite element method. In 2012, Savovic and Djordjevich [2] proposed an explicit finite difference method to solve a one-dimensional advection-diffusion equation with variable coefficients in semi-infinite media for three dispersion problems. In 2013, Appadu and Gidey [3] proposed two time-splitting procedures that they used to solve a 2D advection-diffusion equation with constant coefficients. In 2013, Bause and Schwegler [4] proposed a higher order finite element approximation for systems of coupled convection-dominated transport equations. In 2015, Mojtabi and Deville [5] proposed a finite element method to solve analytically using separation of variables and numerical solution of a time-dependent one-dimensional linear advection-diffusion equation with Dirichlet homogeneous boundary conditions and an initial sine function. In 2016, Gharehbaghi [6] proposed an explicit and an implicit differential quadrature method to compute a numerical solution of a one-dimensional time-dependent advection-diffusion equation with variable coefficients in a semi-infinite domain. In 2017, Gyrya and Lipnikov [7] proposed an arbitrary order mimetic finite difference discretization for the diffusion equation with non-symmetric tensorial diffusion coefficient in a mixed formulation on general polygonal meshes. In 2017, Bahar and Gurarslan [8] studied the effects of Lie-Trotter and Strang splitting methods on the solution of a one-dimensional advection-diffusion equation. In 2017, Oyjindal and Pochai [9] proposed models for numerical simulations of air pollution measurements near an industrial zone. The numerical experiments consisted of different conditions including normal and controlled emissions. The models were then solved using an explicit finite difference technique and the solutions were compared with the measurements for the controlled and uncontrolled point sources at the monitoring points. In 2018, Suebyat and Pochai [10] developed a three-dimensional air pollution measurement model for a heavy traffic area under a Bangkok sky train platform and used finite difference techniques to approximate the solutions. In 2018, Al-Jawary et al. [11] proposed a semi-analytical technique for finding the exact solutions for different types of Burger’s equations and systems of equations in 1D, 2D, and 3D. Also, the method was applied to solve the diffusion and advection-diffusion equations. In 2018, Kusuma et al. [12] proposed a finite difference (FTCS) method to solve a model of pollution distribution in a street tunnel using two dimensional advection and three dimensional diffusion in a rectangular box domain. In 2018, Lou et al. [13] developed reconstructed Discontinuous Galerkin methods for solving linear advection-diffusion equations on hybrid unstructured grids based on a first-order hyperbolic system formulation. In 2018, Bhatt et al, [14] developed a Krylov subspace approximation-based locally extrapolated exponential time differencing method and studied its accuracy and efficiency for solving three-dimensional nonlinear advection-diffusion-reaction systems. In 2020, Pananu et al. [15] analyzed the convergence of the finite difference method with the implicit forward time central space (FTCS) scheme for the two-dimensional advection-diffusion-reaction equation (ADRE) and applied the scheme to a pollutant dispersion with removal mechanism model in a reservoir. In 2020, Heng and Guodong [16] improved the element-free Galerkin method and used it for solving 3D advection-diffusion problems. In 2020, Cruz-Quintero and Jurado [17] proposed a backstepping design for the boundary control of a reaction-advection-diffusion equation with constant coefficients and Neumann boundary conditions. In 2021, Para et al. [18] proposed the characteristic finite volume method for solving a convection-diffusion problem on two-dimensional triangular grids and compared the accuracy of four piecewise linear reconstruction techniques on structured triangular grids, namely, Frink, Holmes-Connell, Green-Gauss, and least squares methods. In 2021, Hong et al. [19] discussed the numerical solution of 3D unsteady advection-diffusion equations using a meshless numerical scheme with a space-time backward substitution method. In 2021, Irfan and Hidayat [20] proposed a meshless finite difference method with B-splines for the numerical solution of coupled advection-diffusion-reaction problems. In 2021, Shahid et al [21] studied an epidemic type model with advection and diffusion terms for the transmission dynamics of a computer virus model with fixed population density. In 2021, García and Jurado [22] designed an adaptive boundary control for a parabolic type reaction-advection-diffusion PDE under the assumption of unknown parameters for both advection and reaction terms and Robin and Neumann boundary conditions. In 2022, Para et al. [23] proposed an explicit characteristic-based finite volume method (FVM) for the numerical solution of advection-diffusion-reaction equations (ADRE) and applied it to solve some 1D and 2D water pollution problems which can be modeled in terms of ADREs. Then, the FVM results were compared with numerical results obtained using a finite difference method (FDM) with an implicit forward time central space (FTCS) scheme.

In this article, we consider the application of a finite difference method based on the implicit forward time central space (FTCS) scheme to the numerical solution of an advection-diffusion-reaction equation of the following form:

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial t}+\nabla ·(\mathbf{v}\varphi -\epsilon \nabla \varphi )+\kappa \varphi =q,\end{array}$$

(1)

where $\varphi$ is a scalar quantity, $\mathbf{v}=\mathbf{v}\left(\mathbf{x}\right)$ is a given advection velocity vector, $\epsilon \ge 0$ is a diffusion coefficient, $\kappa$ is a reaction coefficient, $q=q(\mathbf{x},t)$ is a prescribed source term, $\mathbf{x}$ is a position vector and time $t\in [0,T]$, where $T>0$. The present paper is organized as follows. The convergence theory is briefly provided in Section 2. The consistency and stability analysis and numerical results of the method are given in Section 3. Using the method, numerical results and their graphs for the air pollution problems written in terms of the ADRE in Eq. (1) are shown in Section 4. The conclusion of our work is presented in Section 5.

2. Convergence Theory

We consider the following initial time, space boundary value problem (IBVP) for the advection-diffusion-reaction equation

$$L\varphi (\mathbf{x},t)=\frac{\partial \varphi (\mathbf{x},t)}{\partial t}+\nabla ·\left(\mathbf{v}\left(\mathbf{x}\right)\varphi (\mathbf{x},t)-\epsilon \nabla \varphi (\mathbf{x},t)\right)+\kappa \varphi (\mathbf{x},t)=q(\mathbf{x},t),$$

(2)

where $(\mathbf{x},t)\in \mathsf{\Omega }\times (0,T]$, and with the initial time condition

$$\varphi (\mathbf{x},0)=\psi \left(\mathbf{x}\right)\mathrm{for}\mathbf{x}\in \mathsf{\Omega },$$

(3)

and either Dirichlet or Neumann boundary space conditions

$$\begin{array}{ccc}\hfill \mathrm{Dirichlet}& & \varphi (\mathbf{x},t)={\xi }_D(\mathbf{x},t)\hfill \\ \hfill \mathrm{Neumann}& & \epsilon {\displaystyle \frac{\partial \varphi (\mathbf{x},t)}{\partial \mathbf{n}}}={\xi }_N(\mathbf{x},t)\mathrm{for}(\mathbf{x},t)\in \partial \mathsf{\Omega }\times [0,T].\hfill \end{array}$$

(4)

In Eqs. (2)–(4), T is a finite time, $\partial \mathsf{\Omega }$ is the boundary of the spatial domain $\mathsf{\Omega }$, ${\xi }_D$ and ${\xi }_N$ are given functions, $\mathbf{v}\left(\mathbf{x}\right)$ is an advection vector, $\epsilon >0$ is a diffusion coefficient, $\kappa$ is a reaction constant and $q(\mathbf{x},\mathbf{t})$ is a source term. In addition, L is a linear partial differential operator from a space of continuous functions F to a space of continuous functions H, where the dependent variable $\varphi (\mathbf{x},t)\in F$, and where $L\varphi (\mathbf{x},t)\in H$ and $q(\mathbf{x},t)\in H$.

The first step in developing the finite difference method is to discretize $\mathsf{\Omega }$, $\partial \mathsf{\Omega }$ and $[0,T]$. We assume that there are $N_1$ step sizes of length $\Delta x=\frac{L_1}{N_1}$ in the x direction, $N_2$ step sizes of length $\Delta y=\frac{L_2}{N_2}$ in the y direction and $N_3$ step sizes of length $\Delta z=\frac{L_3}{N_3}$ in the z direction. For simplicity, we will assume that $\Delta x=\Delta y=\Delta z=h$ and therefore $L_1=N_{1}h$, $L_2=N_{2}h$ and $L_3=N_{3}h$. Then the discretized version of $\mathsf{\Omega }$ is the spatial grid

$${\mathsf{\Omega }}_h=\{(x_i,y_j,z_k)|x_i=ih,y_j=jh,z_k=kh,i\in J_1,j\in J_2,k\in J_3\},$$

(5)

where $J_1=\{0,1,2,...,N_1\}$, $J_2=\{0,1,2,...,N_2\}$, $J_3=\{0,1,2,...,N_3\}$. The boundary of the spatial grid is

$$\begin{array}{ccc}\hfill \partial {\mathsf{\Omega }}_h& =& \{x_0=0;y_j=jh,j\in J_2,z_k=kh,k\in J_3\}\hfill \\ & & \cup \{x_{N_1}=L_1;y_j=jh,j\in J_2,z_k=kh,k\in J_3\}\hfill \\ & & \cup \{x_i=ih,i\in J_1;y_0=0;z_k=kh,k\in J_3\}\hfill \\ & & \cup \{x_i=ih,i\in J_1;y_{N_2}=L_2;z_k=kh,k\in J_3\}\hfill \\ & & \cup \{x_i=ih,i\in J_1,y_j=jh,j\in J_2;z_0=0\}\hfill \\ & & \cup \{x_i=ih,i\in J_1,y_j=jh,j\in J_2;z_{N_3}=L_3\}.\hfill \end{array}$$

(6)

Similarly, the interval $[0,T]$ is discretized with N step sizes of size $\Delta t=\frac{T}{N}$. The set of discretized times is then

$$T_K=\{t^n|t^n=n\Delta t,n\in K\},K=\{0,1,2,\dots ,N\}.$$

(7)

Finally, we define $F_h$ and $H_h$ as discretized function spaces with domain ${\mathsf{\Omega }}_h\times T_K$. We also define a projection ${\mathcal{P}}_h\left(\varphi \right)\in F_h$ of the exact solution $\varphi (\mathbf{x},t)$ of (2) onto ${\mathsf{\Omega }}_h\times T_K$ and a projection ${\mathcal{P}}_h\left(q\right)\in H_h$ of the source term $q(\mathbf{x},t)$ onto ${\mathsf{\Omega }}_h\times T_K$.

We can then write a finite difference scheme (FDS) corresponding to the continuous partial differential equation (2) in the form

$$L_h{\varphi }^h(\mathbf{x},t)=q^h(\mathbf{x},t),(\mathbf{x},t)\in {\mathsf{\Omega }}_h\times T_K,$$

(8)

where $L_h:F_h\to H_h$ is a finite difference operator, corresponding to the partial differential operator L, which acts from the discrete function space $F_h$ to the discrete function space $H_h$. The initial condition in (3) can then be rewritten as the discretized initial condition

$${\varphi }^h(\mathbf{x},0)={\psi }^h\left(\mathbf{x}\right),\mathbf{x}\in {\mathsf{\Omega }}_h,$$

(9)

and, for the Dirichlet case, the boundary conditions in (4) can be rewritten as the discretized boundary conditions

$${\varphi }^h(\mathbf{x},t)={\xi }^h(\mathbf{x},t)\mathrm{for}(\mathbf{x},t)\in \partial {\mathsf{\Omega }}_h\times T_K.$$

(10)

Lax’s equivalence theorem [24]

An initial time, space boundary value problem (IBVP) is called a properly posed IBVP if it has a unique solution. In this paper, all advection-diffusion-reaction equations that we study will be properly posed IBVP.

Then Lax’s equivalence theorem states that necessary and sufficient conditions for the solution of a finite difference approximation to converge to the unique solution of the IBVP are the following conditions of convergence, consistency and stability.

Using the definition given above of the projection ${\mathcal{P}}_h$ of a function in the continuous space H onto the discrete space $H_h$, we define the truncation error (T.E.) between the finite difference system Eq. (8) and the IBVP Eq. (2) as follows

$$T.E.={\mathcal{P}}_h\left(L\varphi (\mathbf{x},t)\right)-L_h{\varphi }^h(\mathbf{x},t),(\mathbf{x},t)\in {\mathsf{\Omega }}_h\times T_K.$$

(11)

Introducing the norm in the set of discrete functions as the infinity norm

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \parallel {\varphi }^h{\parallel }_{F_h}}& =\underset{\begin{array}{c}i,j.k\in J\\ n\in K\end{array}}{\max }|{\varphi }^h(x_i,y_j,z_k,t^n)|,\hfill \end{array}\end{array}$$

(12)

we have the definition of convergence for numerical methods as follows

Definition 1.

Convergence

A solution ${\varphi }^h$ of the FDS in Eq. (8)-(10) converges to the solution ϕ of the IBVP in Eq. (2)-(4) if

$$\parallel {\mathcal{P}}_h\left(\varphi \right)-{\varphi }^h{\parallel }_{F_h}\to 0ash\to 0.$$

(13)

Definition 2.

Convergence with order m

The FDS in Eq. (8) converges with order m if

$$\parallel {\mathcal{P}}_h\left(\varphi \right)-{\varphi }^h{\parallel }_{F_h}\le C{h}^m,$$

(14)

where C is a positive constant that does not depend on h.

1) Consistency : A finite difference representation of a PDE is said to be consistent [24] if the difference between the PDE and its difference representation vanishes as the mesh is refined, i.e.,

Definition 3.

The FDS (8) is a consistent approximation of the IBVP (2)-(4) with order m if

$${\parallel T.E.\parallel }_{H_h}\le C{h}^m,$$

(15)

where C is a positive constant that does not depend on h.

2) Stability : The finite difference scheme defined by (8) with linear operator $L_h$ will be called stable, if there exists $h_0>0$ such that for arbitrary $h<h_0$ and for

$$\parallel f^h{\parallel }_{H_h}=\underset{\begin{array}{c}i,j,k\in J\end{array}}{\max }|\psi (x_i,y_j,z_k,t^n)|+\underset{\begin{array}{c}i,j,k\in J\\ n\in K-\{0\}\end{array}}{\max }|q(x_i,y_j,z_k,t^n)|,$$

(16)

the solution ${\varphi }^h$ of the FDS, $L_h{\varphi }^h=q^h$ in Eq. (8), exists and is unique and satisfies the inequality

$$\parallel {\varphi }^h{\parallel }_{F_h}\le C{\parallel f^h\parallel }_{H_h},$$

(17)

where the positive constant C does not depend on h.

The inequality (17) has to be true for any initial conditions $\psi (x_i,y_j,z_{,}0)$ and source terms $q(x_i,y_j,z_k,t^n)$ of Eq. (2). In particular, if $q(x_i,y_j,z_k,t^n)=0$, then the condition (17) reduces to the necessary condition for stability of the homogeneous version of (8).

For the special case of Fourier or Von Neumann Analysis [24], a solution of the FDS (8) for $q(x_i,y_j,z_k,t_n)=0$ can be written in the form

$${\varphi }^h(x_i,y_j,z_k,t^n)={\lambda }^n{e}^{I(\alpha i+\beta j+\eta k)},(i,j,k)\in J,n\in K,I=\sqrt{-1},$$

(18)

where $\alpha$, $\beta$ and $\eta$ are wave numbers and $e^{I(\alpha i+\beta j+\eta k)}$ are eigenfunctions corresponding to eigenvalues $\lambda$ of $L_h$. The necessary condition for stability of FDS (8) for $q(x_i,y_j,z_k,t^n)=0$ will then hold for all $\alpha ,\beta ,\eta \in \Re$ if the following inequality holds:

$$|\lambda |\le 1.$$

(19)

3. Consistency and Stability Analysis

Consider the following initial-boundary value problem

$$\begin{array}{cc}\hfill & \frac{\partial \varphi }{\partial t}+u\frac{\partial \varphi }{\partial x}+v\frac{\partial \varphi }{\partial y}+w\frac{\partial \varphi }{\partial z}-\epsilon \left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}\right)+\kappa \varphi =q(x,y,z,t),\hfill \end{array}$$

$$\begin{array}{cc}\hfill & (x,y,z)\in (0,L_1)\times (0,L_2)\times (0,L_3),t\in (0,T],\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \varphi (x,y,z,0)=\psi (x,y,z),\hfill \\ \hfill & \varphi (0,y,z,t)={\alpha }_1,\varphi (L_1,y,z,t)={\alpha }_2,\hfill \\ \hfill & \varphi (x,0,z,t)={\beta }_1,\varphi (x,L_2,z,t)={\beta }_2,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & \varphi (x,y,0,t)={\eta }_1,\varphi (x,y,L_3,t)={\eta }_2,\hfill \end{array}$$

(21)

in which the above partial differential equation is the component form of (1). Here $\varphi =\varphi (x,y,z,t)$ is the dependent variable and $\mathbf{v}=(u,v,w)$ is the advection velocity where $u,v$ and w are the velocities in the x-, y- and z-directions, respectively.

In this section, we will analyze the convergence of the finite difference method with the implicit FTCS scheme applied to (20). We discretize the domain of the problem: $(x_i,y_j,z_k,t^n)\in {\mathsf{\Omega }}_h\times T_K$, where the spatial grid ${\mathsf{\Omega }}_h$ is defined in (5) and the time grid $T_K=[0,T]$ has points $t^n=n\Delta t$ with $n\in K=\{0,1,2,...,N\}.$

Following the method in Section 2, we use Lax’s equivalence theorem to prove the convergence of the FTCS by proving its consistency and stability. We first show the consistency of the FTCS. The implicit FTCS scheme of Eq. (20) is as follows.

$$\begin{array}{cc}\hfill L_h{\varphi }^h=& \frac{{\varphi }_{i,j,k}^{n+1}-{\varphi }_{i,j,k}^n}{\Delta t}+u\frac{{\varphi }_{i+1,j,k}^{n+1}-{\varphi }_{i-1,j,k}^{n+1}}{2\Delta x}+v\frac{{\varphi }_{i,j+1,k}^{n+1}-{\varphi }_{i,j-1,k}^{n+1}}{2\Delta y}+w\frac{{\varphi }_{i,j,k+1}^{n+1}-{\varphi }_{i,j,k-1}^{n+1}}{2\Delta z}\hfill \\ \hfill & -\epsilon \left[\frac{{\varphi }_{i+1,j,k}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i-1,j,k}^{n+1}}{{(\Delta x)}^2}+\frac{{\varphi }_{i,j+1,k}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i,j-1,k}^{n+1}}{{(\Delta y)}^2}\right.\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & \left.+\frac{{\varphi }_{i,j,k+1}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i,j,k-1}^{n+1}}{{(\Delta z)}^2}\right]+\kappa {\varphi }_{i,j,k}^{n+1}=q_{i,j,k}^n,\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \varphi (x_i,y_j,z_k,0)=\psi (x_i,y_j,z_k),\hfill \\ \hfill & \varphi (0,y_j,z_k,t^n)={\alpha }_1,\varphi (L_1,y_j,z_k,t^n)={\alpha }_2,\hfill \\ \hfill & \varphi (x_i,0,z_k,t^n)={\beta }_1,\varphi (x_i,L_2,z_k,t^n)={\beta }_2,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & \varphi (x_i,y_j,0,t^n)={\eta }_1,\varphi (x_i,y_j,L_3,t^n)={\eta }_2,\hfill \end{array}$$

where ${\varphi }_{i,j,k}^n={\varphi }^h(x_i,y_j,z_k,t^n)$ and $q_{i,j,k}^n=q(x_i,y_j,z_k,t^n)$.

The truncation error for (22) obtained using the method is

$$\begin{array}{ccc}\hfill T.E.& =& -\frac{\Delta t}{2!}\frac{{\partial }^2\varphi }{\partial t^2}-u\frac{{(\Delta x)}^2}{3!}\frac{{\partial }^3\varphi }{\partial x^3}-v\frac{{(\Delta y)}^2}{3!}\frac{{\partial }^3\varphi }{\partial y^3}-w\frac{{(\Delta z)}^2}{3!}\frac{{\partial }^3\varphi }{\partial z^3}\hfill \\ & & +\epsilon \left(\frac{2{(\Delta x)}^2}{4!}\frac{{\partial }^4\varphi }{\partial x^4}+\frac{2{(\Delta y)}^2}{4!}\frac{{\partial }^4\varphi }{\partial y^4}+\frac{2{(\Delta z)}^2}{4!}\frac{{\partial }^4\varphi }{\partial z^4}\right)\hfill \\ & & +O({(\Delta t)}^2,{(\Delta x)}^4,{(\Delta y)}^4,{(\Delta z)}^4).\hfill \end{array}$$

(24)

Therefore, from (24), we have

$$\begin{array}{c}\hfill T.E.=O\left((\Delta t),{(\Delta x)}^2,{(\Delta y)}^2,{(\Delta z)}^2\right).\end{array}$$

Taking $(\Delta t,\Delta x,\Delta y,\Delta z)\to 0$, then we have

$$\begin{array}{c}\hfill \underset{(\Delta t,\Delta x,\Delta y,\Delta z)\to 0}{lim}\parallel T.E.\parallel =\underset{(\Delta t,\Delta x,\Delta y,\Delta z)\to 0}{lim}\parallel O\left((\Delta t),{(\Delta x)}^2,{(\Delta y)}^2,{(\Delta z)}^2\right)\parallel =0.\end{array}$$

We will now prove the stability of the finite difference method with the implicit FTCS scheme for numerical solution of (20). In order to illustrate stability of the method, we initially show that condition (19) holds for the homogeneous equation of (22) and then that condition (17) is satisfied for the nonhomogeneous equation (22). For convenience, we set $\tau =\Delta t,h=\Delta x=\Delta y=\Delta z$.

Case I: Homogeneous equation

Setting $q_{i,j,k}^n=0$ in equation (22), we obtain the discretized homogeneous version of (22). Assuming a discretized solution of the resulting equation as ${\varphi }_{i,j,k}^n={\lambda }^n{e}^{I(\alpha i+\beta j+\eta k)}$, $I=\sqrt{-1}$[24], we obtain the following equation for the eigenvalue $\lambda$

$$\begin{array}{ccc}& & \frac{{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta k)}-{\lambda }^n{e}^{I(\alpha i+\beta j+\eta k)}}{\tau }+u\frac{{\lambda }^{n+1}{e}^{I\left(\alpha \right(i+1)+\beta j+\eta k)}-{\lambda }^{n+1}{e}^{I\left(\alpha \right(i-1)+\beta j+\eta k)}}{2h}\hfill \\ & & +v\frac{{\lambda }^{n+1}{e}^{I(\alpha i+\beta (j+1)+\eta k)}-{\lambda }^{n+1}{e}^{I(\alpha i+\beta (j-1)+\eta k)}}{2h}\hfill \\ & & +w\frac{{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta (k+1\left)\right)}-{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta (k-1\left)\right)}}{2h}\hfill \end{array}$$

$$\begin{array}{ccc}& & -\epsilon \left[\frac{{\lambda }^{n+1}{e}^{I\left(\alpha \right(i+1)+\beta j+\eta k)}-2{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta k)}+{\lambda }^{n+1}{e}^{I\left(\alpha \right(i-1)+\beta j+\eta k)}}{h^2}\right.\hfill \\ & & +\frac{{\lambda }^{n+1}{e}^{I(\alpha i+\beta (j+1)+\eta k)}-2{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta k)}+{\lambda }^{n+1}{e}^{I(\alpha i+\beta (j-1)+\eta k)}}{h^2}\hfill \\ & & +\left.\frac{{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta (k+1\left)\right)}-2{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta k)}+{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta (k-1\left)\right)}}{h^2}\right]\hfill \\ & & +\kappa {\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta k)}=0.\hfill \end{array}$$

(25)

Then, we can solve (25) for the eigenvalue $\lambda$ as follows.

$$\begin{array}{ccc}& & {\lambda }^n{e}^{I(\alpha i+\beta j+\eta k)}\left[\frac{\lambda -1}{\tau }\right]-{\lambda }^{n+1}{e}^{I(\alpha i+\beta j+\eta k)}\left[-u\frac{e^{I\alpha }-e^{-I\alpha }}{2h}\right.-v\frac{e^{I\beta }-e^{-I\beta }}{2h}\hfill \\ & & -\left.w\frac{e^{I\eta }-e^{-I\eta }}{2h}+\epsilon \frac{e^{I\alpha }-2+e^{-I\alpha }}{h^2}+\epsilon \frac{e^{I\beta }-2+e^{-I\beta }}{h^2}+\epsilon \frac{e^{I\eta }-2+e^{-I\eta }}{h^2}-\kappa \right]=0,\hfill \\ & & \lambda =\frac{1}{1+\frac{\tau }{h}\left[\frac{4\epsilon }{h}\left({sin}^2\frac{\alpha }{2}+{sin}^2\frac{\beta }{2}+{sin}^2\frac{\eta }{2}\right)+I\left(usin\alpha +vsin\beta +wsin\eta \right)+\kappa h\right]}.\hfill \end{array}$$

Letting $\gamma =\frac{4\epsilon \tau }{h^2}$ and $\mu =\frac{\tau }{h}$, we then obtain,

$$\begin{array}{c}\hfill \lambda =\frac{1}{1+\gamma \left({sin}^2\frac{\alpha }{2}+{sin}^2\frac{\beta }{2}+{sin}^2\frac{\eta }{2}\right)+\mu I\left(usin\alpha +vsin\beta +wsin\eta \right)+\kappa \tau }.\end{array}$$

Consequently, the magnitude of $\lambda$ is

$$\begin{array}{ccc}\hfill |\lambda |& =& \frac{1}{\sqrt{{\left[1+\gamma \left({sin}^2\frac{\alpha }{2}+{sin}^2\frac{\beta }{2}+{sin}^2\frac{\eta }{2}\right)+\kappa \tau \right]}^2+{\left[\mu \left(usin\alpha +vsin\beta +wsin\eta \right)\right]}^2}},\hfill \\ & \le & 1.\hfill \end{array}$$

Therefore, condition (19) for Von Neumann stability, which is independent of $\tau$ and h, has been proved.

Case II: Nonhomogeneous equation

We now prove that condition (17) is satisfied for the nonhomogeneous equation (22) with $q_{i,j,k}^n\ne 0$ by finding a solution ${\varphi }_{i,j,k}^{n+1}$ for the linear system of difference equations (22) and (23) if ${\varphi }_{i,j,k}^n$ is assumed known. Substituting $\Delta x=\Delta y=\Delta z=h$ and $\Delta t=\tau$ in (22) and (23), we obtain

$$\begin{array}{ccc}\hfill L_h{\varphi }^h& =& \frac{{\varphi }_{i,j,k}^{n+1}-{\varphi }_{i,j,k}^n}{\tau }+u\frac{{\varphi }_{i+1,j,k}^{n+1}-{\varphi }_{i-1,j,k}^{n+1}}{2h}+v\frac{{\varphi }_{i,j+1,k}^{n+1}-{\varphi }_{i,j-1,k}^{n+1}}{2h}+w\frac{{\varphi }_{i,j,k+1}^{n+1}-{\varphi }_{i,j,k-1}^{n+1}}{2h}\hfill \\ & & -\epsilon \left[\frac{{\varphi }_{i+1,j,k}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i-1,j,k}^{n+1}}{h^2}+\frac{{\varphi }_{i,j+1,k}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i,j-1,k}^{n+1}}{h^2}\right.\hfill \end{array}$$

$$\begin{array}{ccc}& & +\left.\frac{{\varphi }_{i,j,k+1}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i,j,k-1}^{n+1}}{h^2}+\frac{{\varphi }_{i,j,k+1}^{n+1}-2{\varphi }_{i,j,k}^{n+1}+{\varphi }_{i,j,k-1}^{n+1}}{h^2}\right]+\kappa {\varphi }_{i,j,k}^{n+1}=q_{i,j,k}^n\hfill \\ & & \mathrm{subject}\mathrm{to}{\varphi }_{0,j,k}^n={\alpha }_1,{\varphi }_{N_1,j,k}^n={\alpha }_2;j=0,1,...,N_2;k=0,1,...,N_3\hfill \\ & & {\varphi }_{i,0,k}^n={\beta }_1,{\varphi }_{i,N_2,k}^n={\beta }_2;i=0,1,...,N_1;k=0,1,...,N_3\hfill \end{array}$$

(26)

$$\begin{array}{ccc}& & {\varphi }_{i,j,0}^n={\eta }_1,{\varphi }_{i,j,N_3}^n={\eta }_2;i=0,1,...,N_1;j=0,1,...,N_2\hfill \end{array}$$

(27)

After rearranging, we can rewrite (26) in the form

$$a{\varphi }_{i-1,j,k}^{n+1}+b{\varphi }_{i,j-1,k}^{n+1}+c{\varphi }_{i,j,k-1}^{n+1}+d{\varphi }_{i,j,k}^{n+1}+e{\varphi }_{i+1,j,k}^{n+1}+f{\varphi }_{i,j+1,k}^{n+1}+g{\varphi }_{i,j,k+1}^{n+1}={\varphi }_{i,j,k}^n,$$

(28)

$$\begin{array}{cc}\hfill \mathrm{where}& a=\left(\frac{u\tau }{2h}+\frac{\epsilon \tau }{h^2}\right),b=\left(\frac{v\tau }{2h}+\frac{\epsilon \tau }{h^2}\right),c=\left(\frac{w\tau }{2h}+\frac{\epsilon \tau }{h^2}\right),\hfill \\ \hfill & d=\left(-1-\frac{6\epsilon \tau }{h^2}-\kappa \tau \right),e=\left(-\frac{u\tau }{2h}+\frac{\epsilon \tau }{h^2}\right),f=\left(-\frac{v\tau }{2h}+\frac{\epsilon \tau }{h^2}\right),\hfill \\ \hfill & g=\left(-\frac{w\tau }{2h}+\frac{\epsilon \tau }{h^2}\right),{\varphi }_{i,j,k}^n=-{\varphi }_{i,j,k}^n-\tau q_{i,j,k}^n,\hfill \end{array}$$

and where the step sizes $\tau$ and h can be chosen so that the coefficients $a,b,c,d,e,f,g$ satisfy the conditions

$$|d|>|a|+|b|+|c|+|e|+|f|+|g|+\delta ,\delta >0.$$

(29)

Lemma 1.

If the coefficients of Eq.(28) satisfy the conditions Eq.(29) then the solution of Eq.(28) exists and is unique and satisfies the inequality

$$|{\varphi }_{i,j,k}^{n+1}|\le \max \left\{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\frac{1}{\delta }\underset{p,q,r}{\max }|{\varphi }_{p,q,r}^n|\right\}.$$

(30)

Proof 

First, we prove inequality (30). Assume that $|{\varphi }_{p,q,r}^{n+1}|=\max \{|{\varphi }_{i,j,k}^{n+1}|,i=0,1,...,N_1,j=0,1,...,N_2,k=0,1,...,N_3\}$.

Let $0<p<N_1,0<q<N_2,0<r<N_3$. Then

$$\begin{array}{ccc}\hfill |d||{\varphi }_{p,q,r}^{n+1}|& =& |-a{\varphi }_{p-1,q,r}^{n+1}-b{\varphi }_{p,q-1,r}^{n+1}-c{\varphi }_{p,q,r-1}^{n+1}-e{\varphi }_{p+1,q,r}^{n+1}-f{\varphi }_{p,q+1,r}^{n+1}-g{\varphi }_{p,q,r+1}^{n+1}+{\varphi }_{p,q,r}^n|\hfill \\ & \le & |a{\varphi }_{p-1,q,r}^{n+1}|+|b{\varphi }_{p,q-1,r}^{n+1}|+|c{\varphi }_{p,q,r-1}^{n+1}|+|e{\varphi }_{p+1,q,r}^{n+1}|+|f{\varphi }_{p,q+1,r}^{n+1}|+|g{\varphi }_{p,q,r+1}^{n+1}|\hfill \\ & & +|{\varphi }_{p,q,r}^n|\hfill \\ & =& |a||{\varphi }_{p-1,q,r}^{n+1}|+|b||{\varphi }_{p,q-1,r}^{n+1}|+|c||{\varphi }_{p,q,r-1}^{n+1}|+|e||{\varphi }_{p+1,q,r}^{n+1}|+|f||{\varphi }_{p,q+1,r}^{n+1}|\hfill \\ & & +|g||{\varphi }_{p,q,r+1}^{n+1}|+|{\varphi }_{p,q,r}^n|\hfill \\ & \le & (|a|+|b|+|c|+|e|+|f|+|g|)|{\varphi }_{p,q,r}^{n+1}|+|{\varphi }_{p,q,r}^n|,\hfill \\ \hfill |{\varphi }_{p,q,r}^{n+1}|& \le & \frac{|{\varphi }_{p,q,r}^n|}{|d|-|a|-|b|-|c|-|e|-|f|-|g|}\hfill \\ & \le & \frac{|{\varphi }_{p,q,r}^n|}{\delta }.\hfill \end{array}$$

This completes the proof.

Now we will prove the stability of the implicit FTCS scheme by showing that inequality (17) is satisfied.

Proof 

Let $i^{*}$, $j^{*}$ and $k^{*}$ be the smallest non-negative integers such that

$${\displaystyle |{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}|=\underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^{n+1}|,}$$

(31)

where $J=\{(i,j,k)|i=0,1,2,...,N_1,j=0,1,2,...,N_2,k=0,1,2,...,N_3\}.$ Now let $\overline{J}$ be the discretized boundary set of J, i.e.,

$$\begin{array}{ccc}\hfill \overline{J}& =& \{(i,j,k)\in J|i=0,N_1\}\cup \{(i,j,k)\in J|i\ne 0,i\ne N_1,j=0,N_2\}\hfill \\ & & \cup \{(i,j,k)\in J|i\ne 0,i\ne N_1,j\ne 0,j\ne N_2,k=0,N_3\}.\hfill \end{array}$$

(32)

It is obvious that if $(i^{*},j^{*},k^{*})\in \overline{J},$ then we have

$$\underset{(i,j,k)\in J}{\max }|{\varphi }_{i,j,k}^{n+1}|\le \max \{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|\}.$$

(33)

If $(i^{*},j^{*},k^{*})\in \mathcal{J}=J-\overline{J},$ then replacing $(i,j,k)$ in (26) with $(i^{*},j^{*},k^{*})$ we obtain

$$\begin{array}{ccc}\hfill & & -{\varphi }_{i^{*},j^{*},k^{*}}^n-\tau q_{i^{*},j^{*},k^{*}}^n=\left(\frac{u\tau }{2h}+\frac{\epsilon \tau }{h^2}\right){\varphi }_{i^{*}-1,j^{*},k^{*}}^{n+1}+\left(\frac{v\tau }{2h}+\frac{\epsilon \tau }{h^2}\right){\varphi }_{i^{*},j^{*}-1,k^{*}}^{n+1}+\left(\frac{w\tau }{2h}+\frac{\epsilon \tau }{h^2}\right){\varphi }_{i^{*},j^{*},k^{*}-1}^{n+1}\hfill \\ & & -\left(1+\frac{6\epsilon \tau }{h^2}+\kappa \tau \right){\varphi }_{i^{*},j^{*},k^{*}}^{n+1}-\left(\frac{u\tau }{2h}-\frac{\epsilon \tau }{h^2}\right){\varphi }_{i^{*}+1,j^{*},k^{*}}^{n+1}\hfill \\ & & -\left(\frac{v\tau }{2h}-\frac{\epsilon \tau }{h^2}\right){\varphi }_{i^{*},j^{*}+1,k^{*}}^{n+1}-\left(\frac{w\tau }{2h}-\frac{\epsilon \tau }{h^2}\right){\varphi }_{i^{*},j^{*},k^{*}+1}^{n+1}.\hfill \end{array}$$

(34)

Next we make the following hypothesis:

$$\begin{array}{c}\hfill \left(\mathbf{H}\mathbf{1}\right):\begin{array}{cc}& u\left[{\varphi }_{i^{*}-1,j^{*},k^{*}}^{n+1}-{\varphi }_{i^{*}+1,j^{*},k^{*}}^{n+1}\right]\le 0,\hfill \\ & v\left[{\varphi }_{i^{*},j^{*}-1,k^{*}}^{n+1}-{\varphi }_{i^{*},j^{*}+1,k^{*}}^{n+1}\right]\le 0,\hfill \\ & w\left[{\varphi }_{i^{*},j^{*},k^{*}-1}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}+1}^{n+1}\right]\le 0.\hfill \end{array}\end{array}$$

(35)

Without loss of generality we assume that ${\varphi }_{i^{*},j^{*},k^{*}}^{n+1}>0$ and further assume that the hypothesis (H1) holds. We can estimate the right side of (34) as

$$\begin{array}{cc}\hfill & \frac{\epsilon \tau }{h^2}\left[{\varphi }_{i^{*}-1,j^{*},k^{*}}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\right]+\frac{\epsilon \tau }{h^2}\left[{\varphi }_{i^{*},j^{*}-1,k^{*}}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\right]+\frac{\epsilon \tau }{h^2}\left[{\varphi }_{i^{*},j^{*},k^{*}-1}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\right]\hfill \\ \hfill & +\frac{\epsilon \tau }{h^2}\left[{\varphi }_{i^{*}+1,j^{*},k^{*}}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\right]+\frac{\epsilon \tau }{h^2}\left[{\varphi }_{i^{*},j^{*}+1}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\right]+\frac{\epsilon \tau }{h^2}\left[{\varphi }_{i^{*},j^{*},k^{*}+1}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\right]\hfill \\ \hfill & +\frac{u\tau }{2h}\left[{\varphi }_{i^{*}-1,j^{*},k^{*}}^{n+1}-{\varphi }_{i^{*}+1,j^{*},k^{*}}^{n+1}\right]+\frac{v\tau }{2h}\left[{\varphi }_{i^{*},j^{*}-1,k^{*}}^{n+1}-{\varphi }_{i^{*},j^{*}+1,k^{*}}^{n+1}\right]\hfill \\ \hfill & +\frac{w\tau }{2h}\left[{\varphi }_{i^{*},j^{*},k^{*}-1}^{n+1}-{\varphi }_{i^{*},j^{*},k^{*}+1}^{n+1}\right]-{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}-\kappa \tau {\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\hfill \\ \hfill & \le -{\varphi }_{i^{*},j^{*},k^{*}}^{n+1}.\hfill \end{array}$$

(36)

Then, substituting (36) into the right hand side of (34), we obtain

$${\varphi }_{i^{*},j^{*},k^{*}}^{n+1}\le {\varphi }_{i^{*},j^{*},k^{*}}^n+\tau q_{i^{*},j^{*},k^{*}}^n.$$

(37)

Therefore,

$$\begin{array}{c}\hfill \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^{n+1}|\le \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^n+\tau q_{i,j,k}^n|\le \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^n|+\tau \underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^n|,\end{array}$$

(38)

where $K=\{0,1,2,...,N\}.$ We now establish the following inequality representing the maximum principle as

$$\begin{array}{c}\hfill {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^{n+1}|\le \max \left\{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^n|+\tau \underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^n|\right\}.}\end{array}$$

(39)

We will now prove that for the given initial and boundary conditions in (27), the maximum solution

${\displaystyle \underset{(i,j,k)\in J}{\max }|{\varphi }_{i,j,k}^n|}$ is bounded for all $n\le N$ for finite $t^N=T$.

At time $t=t^0=0$, the maximum $\varphi (x_i,y_j,z_k,0)$ is ${\displaystyle \underset{(i,j,k)\in J}{\max }\psi (x_i,y_j,z_k,0)}$.

Then, from (39), the maximum at time $t=t^1=\tau$ is

$${\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^1|\le \max \left\{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\underset{\begin{array}{c}(i,j,k)\in J\\ \end{array}}{\max }|{\psi }_{i,j,k}|+\tau \underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^0|\right\}.}$$

(40)

Similarly, the maximum at time $t=t^2=2\tau$ is

$$\begin{array}{ccc}\hfill {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^2|}& \le & {\displaystyle \max \left\{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^1|+\tau \underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^1|\right\}}\hfill \\ & =& \max \{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\hfill \\ & & {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\psi }_{i,j,k}|+\tau (\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^0|+\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^1|)\},}\hfill \\ & \le & \max \{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\hfill \\ & & {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\psi }_{i,j,k}|+2\tau (\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^0|+\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^1|)\}.}\hfill \end{array}$$

(41)

Then, by iteration we find that the maximum at time $t=t^n$ is

$$\begin{array}{ccc}\hfill {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^n|}& \le & \max \{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\hfill \\ & & {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\psi }_{i,j,k}|+n\tau \max (\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^0|,\dots ,\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^{n-1}|)\},}\hfill \end{array}$$

(42)

and the maximum at time $t=t^N=T$ is

$$\begin{array}{ccc}\hfill {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\varphi }_{i,j,k}^N|}& \le & \max \{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\hfill \\ & & {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\psi }_{i,j,k}|+N\tau \max (\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^0|,\dots ,\underset{\begin{array}{c}(i,j,k)\in J\\ n\in K\end{array}}{\max }|q_{i,j,k}^{N-1}|)\},}\hfill \end{array}$$

(43)

Then, as in Eq. (16), for $n\in K$, we can define

$$\begin{array}{ccc}\hfill {\displaystyle {\parallel f\parallel }_{H_h}}& =& \underset{n\in K}{\max }\{|{\alpha }_1|,|{\alpha }_2|,|{\beta }_1|,|{\beta }_2|,|{\eta }_1|,|{\eta }_2|,\hfill \\ & & {\displaystyle \underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|{\psi }_{i,j}|+n\tau \max \{\underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|q_{i,j,k}^0|,\dots ,\underset{\begin{array}{c}(i,j,k)\in J\end{array}}{\max }|q_{i,j,k}^{n-1}|\}\}.}\hfill \end{array}$$

(44)

The above inequality is true for all $n\in K$, and hence we have that

$$\parallel {\varphi }^h{\parallel }_{{\mathsf{\Omega }}_h}\le C{\parallel f^h\parallel }_{H_h},$$

(45)

$$\mathrm{where}C=1\mathrm{and}\parallel {\varphi }^h{\parallel }_{{\mathsf{\Omega }}_h}=\underset{\begin{array}{c}i,j,k\in J\\ n\in K\end{array}}{\max }|{\varphi }_{i,j,k}^n|.$$

(46)

Hence, inequality (17) is satisfied and the solution is stable. This completes the proof.

4. Numerical Results

In this section, we use the implicit forward time central space (FTCS) finite difference method to solve two cases (homogeneous and nonhomogeneous) of three-dimensional ADRE problems of air pollution. In three dimensions, the advection-diffusion-reaction equation can be written as follows:

$$\begin{array}{cc}\hfill & \frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}+w\frac{\partial C}{\partial z}=D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}+D_z\frac{{\partial }^{2}C}{\partial z^2}-RC(x,y,z,t)+Q(x,y,z,t),\hfill \end{array}$$

(47)

where $C=C(x,y,z,t)$ [kg/m${}^3$] is the air pollutant concentration at $(x,y,z)$ [m] and time t s, u, v, w are wind velocity components [m/s] in $x-$, $y-$ and $z-$directions, respectively, $D_x$ and $D_y$ are constant diffusion coefficients in the horizontal direction [m${}^2$/s], $D_z$ is a constant diffusion coefficient in the $z-$direction (vertical) [m${}^2$/s], R is a reaction coefficient, and $Q(x,y,z,t)$ is the rate of change of concentrations of sources or sinks of air pollutants [kg/m${}^3·$s].

In this section, we use an implicit forward time central space (FTCS) scheme to solve the problems. Using the results from section 2, we can write the finite diference equation as

$$\begin{array}{ccc}\hfill & & \frac{C_{i,j,k}^{n+1}-C_{i,j,k}^n}{\Delta t}+u\frac{C_{i+1,j,k}^{n+1}-C_{i-1,j,k}^{n+1}}{2\Delta x}+v\frac{C_{i,j+1,k}^{n+1}-C_{i,j-1,k}^{n+1}}{2\Delta y}+w\frac{C_{i,j,k+1}^{n+1}-C_{i,j,k-1}^{n+1}}{2\Delta z}\hfill \\ & & -D_x\left[\frac{C_{i+1,j,k}^{n+1}-2C_{i,j,k}^{n+1}+C_{i-1,j,k}^{n+1}}{{(\Delta x)}^2}\right]-D_y\left[\frac{C_{i,j+1,k}^{n+1}-2C_{i,j,k}^{n+1}+C_{i,j-1,k}^{n+1}}{{(\Delta y)}^2}\right]\hfill \\ & & -D_z\left[\frac{C_{i,j,k+1}^{n+1}-2C_{i,j,k}^{n+1}+C_{i,j,k-1}^{n+1}}{{(\Delta z)}^2}\right]+R{C}_{i,j,k}^{n+1}=Q_{i,j,k}^n.\hfill \end{array}$$

(48)

We assume that the initial condition is

$$C(x,y,z,0)=0,0\le x\le 1;0\le y\le 1;0\le z\le 1,$$

(49)

and that the boundary conditions are as shown in Table 1.

After rearrangement, we can rewrite (48) in the form

$$\begin{array}{c}\hfill C_{i,j,k}^{n+1}+S_1{C}_{i-1,j,k}^{n+1}+S_2{C}_{i,j-1,k}^{n+1}+S_3{C}_{i,j,k-1}^{n+1}+S_4{C}_{i+1,j,k}^{n+1}+S_5{C}_{i,j+1,k}^{n+1}+S_6{C}_{i,j,k+1}^{n+1}\\ \hfill =S_7{C}_{i,j,k}^n+S_8{Q}_{i,j,k}^n,\end{array}$$

(50)

where

$$\begin{array}{cc}\hfill & S_1=\left(-\frac{u\Delta t}{2\Delta x}-\frac{D_x\Delta t}{{(\Delta x)}^2}\right)/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \\ \hfill & S_2=\left(-\frac{v\Delta t}{2\Delta y}-\frac{D_y\Delta t}{{(\Delta y)}^2}\right)/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \\ \hfill & S_3=\left(-\frac{w\Delta t}{2\Delta z}-\frac{D_z\Delta t}{{(\Delta z)}^2}\right)/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \\ \hfill & S_4=\left(\frac{u\Delta t}{2\Delta x}-\frac{D_x\Delta t}{{(\Delta x)}^2}\right)/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & S_5=\left(\frac{v\Delta t}{2\Delta y}-\frac{D_y\Delta t}{{(\Delta y)}^2}\right)/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \\ \hfill & S_6=\left(\frac{w\Delta t}{2\Delta z}-\frac{D_z\Delta t}{{(\Delta z)}^2}\right)/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \\ \hfill & S_7=1/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right)\hfill \\ \hfill & S_8=\Delta t/\left(1+\frac{2D_x\Delta t}{{(\Delta x)}^2}+\frac{2D_y\Delta t}{{(\Delta y)}^2}+\frac{2D_z\Delta t}{{(\Delta z)}^2}+R\Delta t\right).\hfill \end{array}$$

Case 1 : We consider the three-dimensional advection-diffusion equation (47). By taking $w=0,R=0,Q=0,L=1,W=1$ and $H=1,$ this reduces to a two-dimensional advection three-dimensional diffusion equation for transport of pollutants in the street tunnel problem discussed in [12,25]

$$\begin{array}{c}\hfill \frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}+D_z\frac{{\partial }^{2}C}{\partial z^2},0<t<T.\end{array}$$

(51)

This system of partial differential equation and initial condition (49) and boundary condition Table 1 comes from a model of pollution distribution in a street tunnel, where the wind flows steadily in the x and y directions and there is no flux of pollutant through the solid side-walls or the solid base and roof of the tunnel. We use the nondimensionalised parameters $\Delta x=\Delta y=\Delta z=0.1,\Delta t=0.005,D_x=D_y=D_z=0.02,0.06,0.2,0.5,u=0.6,v=0.4$ and time $T=20$. The numerical solutions and the contour plots of the problem obtained using the implicit FTCS scheme in Eq. (48) are plotted in Figure 1, Figure 2, Figure 3 and Figure 4. We have tested the accuracy of our numerical solutions by comparing them with that of Kusuma et al [12] and found that there was good agreement.

Case 2: In this problem, We consider the three-dimensional advection-diffusion equation (47). By taking $w=0,L=1,W=1$ and $H=1,$ with a decay rate (R) and source term (Q). The advection-diffusion-reaction equation that we consider is as follows.

$$\begin{array}{c}\hfill \frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_x\frac{{\partial }^{2}C}{\partial x^2}+D_y\frac{{\partial }^{2}C}{\partial y^2}+D_z\frac{{\partial }^{2}C}{\partial z^2}-RC+Q(x,y,t).\end{array}$$

(52)

This case of partial differential equation and initial condition (49) and boundary condition Table 1 comes from a model of pollution distribution in a street tunnel, where the wind flows steadily in the x and y directions and there is no flux of pollutant through the solid side-walls or the solid base and roof of the tunnel. We use the nondimensionalised parameters $\Delta x=\Delta y=\Delta z=0.1,\Delta t=0.005,D_x=D_y=D_z=0.2,u=0.6,v=0.4$ and time $T=20$.

We have computed the numerical solutions of the above system using the implicit FTCS method for the following cases.

• $R=0.05$ and $Q=0,0.007,-e^{-t},0.001sin\left(xt\right)$. The results for this case are shown in Figure 5 (a) - (d) at a height $z=0.2$ meters for $Q=0,0.007,-e^{-t}$ and $Q=0.001sin\left(xt\right)$, respectively.
• $R=0.1$ and $Q=0,0.007,-e^{-t},0.001sin\left(xt\right)$. The results for this case are shown in Figure 6 (a) - (d) at a height $z=0.2$ meters for $Q=0,0.007,-e^{-t}$ and $Q=0.001sin\left(xt\right)$, respectively.
• $R=0.5$ and $Q=0,0.007,-e^{-t},0.001sin\left(xt\right)$. The results for this case are shown in Figure 7 (a) - (d) at a height $z=0.2$ meters for $Q=0,0.007,-e^{-t}$ and $Q=0.001sin\left(xt\right)$, respectively.
• $z=0.2$ m for $R=0.05$ and $Q=0,0.007,-e^{-t},0.001sin\left(xt\right)$. The results for this case are shown in Figure 8 (a) at $y=0.5$ m and (b) at $x=0.5$ m.
• $z=0.2$ m for $Q=0.07$ and $R=0.05,0.1,0.5$. The results for this case are shown in Figure 9 (a) at $y=0.5$ m and (b) at $x=0.5$ m.

It can be seen that the numerical solutions for the different Q and R values show similar qualitative behavior but with differences in detail.

5. Conclusions

In this work, we have studied a finite difference method (FDM) for obtaining numerical solutions of advection-diffusion-reaction equations (ADRE) and applied the methods to solve air pollutant problems. For the finite difference method we have used the implicit forward time central space (FTCS) scheme. For the applications, we have solved 3-D air pollutant concentration problems for tunnels. We have investigated the stability, consistency and convergence of the finite difference method. We applied the implicit FTCS finite difference method to solve 3-D ADRE models for three problems of air pollution in tunnels for a range of different source terms Q and reaction terms R. For case 1, we checked the accuracy of our solutions by comparing them with previously published results. We assumed that $Q=R=0$, that there was advection in the x and y directions and diffusion in the x, y and z directions. We then compared our numerical solutions of $D_x=D_y=D_z=0.2$ with Kusuma et al [12] and found that there was good agreement. For case 2, we assumed that Q and R were non-zero, that there was advection in the $x-$ and $y-$ direction and diffusion in the x, y and z directions. We found that the effects of the Q and R were to rapidly reduce the pollutant concentrations in the tunnel.