Numerical solution of 2D singularly perturbed reaction–diffusion system with multiple scales

Abstract

In this article, a robust numerical method is studied to approximate singularly perturbed system of reaction–diffusion problems with multiple scales. The analytical properties of the exact solution have been studied. The numerical method consists of the classical central difference scheme on a Shishkin mesh for spatial semidiscretization processes and the implicit-Euler scheme on a uniform time stepping for temporal derivative. The error estimate is deduced, which exhibits that the numerical approximation is uniformly convergent of almost second-order in spatial variable and first-order in temporal variable. Numerical experiments are given which reveals the effectiveness of the proposed scheme.

Introduction

Singular perturbation problems (SPP) emerge regularly in a few branches of engineering and applied mathematics. The interesting examples of these system are the reaction–diffusion enzyme model, neutron transport problem and model for saturated flow in fractured porous media. For instance, consider the following model for saturated flow in fractured porous media from [1]:

$$\left\{\begin{array}{c}({\beta }_{c_1}+m_1\beta )\frac{\partial p_1}{\partial t}-\frac{k_1}{\mu }\Delta p_1+\frac{\beta }{\mu }(p_1-p_2)=f_1(x,t),\hfill \\ ({\beta }_{c_2}+m_2\beta )\frac{\partial p_2}{\partial t}-\frac{k_2}{\mu }\Delta p_2+\frac{\beta }{\mu }(p_2-p_1)=f_2(x,t),\hfill \end{array}\right.$$

where $p_1,p_2$ are the pressure of liquid in the pores of the first and second-order respectively, $\beta$ is the coefficient of compressibility of the liquid and $\mu$ is the viscosity of the liquid. Here ${\beta }_{c_1}$ and ${\beta }_{c_2}$ are positive constants, whereas $k_1,k_2$ are the porosity of the system of pores of first and second-order respectively, and $m_1,m_2$ are the values of the first and second-order porosity at standard pressure. More details on application of these system can be seen in [2, Chapter 12].

In this article, the system of 2D reaction–diffusion initial–boundary-value problem (IBVP) on the domain $\mathcal{G}:=\mathcal{D}\times (0,T],\mathcal{D}={\Omega }_x\times {\Omega }_y={(0,1)}^2$ is considered:

$$\left\{\begin{array}{c}L\overrightarrow{u}\equiv \frac{\partial \overrightarrow{u}}{\partial t}(x,y,t)+L_{\mathcal{E}}\overrightarrow{u}(x,y,t)=\overrightarrow{f}(x,y,t),(x,y,t)\in \mathcal{G},\hfill \\ \overrightarrow{u}(x,y,0)={\overrightarrow{u}}_0(x,y),(x,y)\in \mathcal{D},\hfill \\ \overrightarrow{u}(x,y,t)=\overrightarrow{0},(x,y,t)\in \{0,1\}\times [0,1]\times [0,T],\hfill \\ \overrightarrow{u}(x,y,t)=\overrightarrow{0},(x,y,t)\in [0,1]\times \{0,1\}\times [0,T],\hfill \end{array}\right.$$

where the spatial differential operator $L_{\mathcal{E}}$ is given by

$$L_{\mathcal{E}}\equiv -{\mathcal{E}}^2\Delta +B(x,y,t).$$

The diffusion and reaction coefficients are given by $\mathcal{E}=\text{diag}({\epsilon }_1,{\epsilon }_2),0<{\epsilon }_1\le {\epsilon }_2\ll 1$ and $B={\left\{b_{ij}\right\}}_{i,j=1}^2$, respectively. Here, we assume that $B$ is strongly diagonally dominant with

$$\left\{b_{ii}>0,{\Vert \frac{b_{ij}}{b_{ii}}\Vert }_{\infty }<1,\text{for}i\ne j.\right.$$

In particular, assume that there exist $\xi$ and $\beta =\beta \left(\xi \right)$ such that

$${\Vert \frac{b_{ij}}{b_{ii}}\Vert }_{\infty }<\xi <1,{\beta }^2≔(1-\xi )min\{b_{11},b_{22}\}.$$

In the past, SPPs with boundary and interior layers have been studied and analyzed with a wide range of numerical schemes. For instances; Natesan et al. [3], [4] developed numerical approximation for turning point problem and weak boundary layer problems. A parallel boundary-value technique has been developed for SPP by Vigo-Aguiar and Natesan in [5]. Simos and Vigo-Aguiar studied the modified Runge Kutta method for the numerical solution of the Schrödinger equation in [6]. Vigo-Aguiar and Natesan [7] analyzed an efficient numerical method for SPP by splitting the model problem in the system of two first-order differential equations. Ramos et al. [8] solved nonlinear singularly perturbed problems by using a non-standard algorithm.

A brief literature about the computational methods for the singularly perturbed system of reaction–diffusion problems (SPSRDP) is presented here. The fitted difference scheme for SPSRDP involving small diffusion parameters of different values is considered in [9], [10]. To improve the accuracy of the numerical approximation of the previous article, a special approach for the error estimate of the numerical method is followed in [11], [12]. In [13], optimal error estimate for SPSRDP by using mesh equidistribution technique is derived. For time-dependent one-dimensional SPSRDP, a uniform convergent numerical scheme is proposed in [14]. A numerical scheme for time-dependent 1D SPSRDP is studied in [15]. Recently, the authors [16] studied upwind based finite difference scheme for time-dependent singularly perturbed system of convection–diffusion problem (SPSCDP) which produces first-order accurate results. To improve the order of convergence of the approximate solution, we applied extrapolation idea in [17].

There are very few articles available in the literature, which deals with the numerical solution of singularly perturbed 2D parabolic IBVPs. In [18], a monotone finite difference scheme is considered to approximate the solution of time-dependent SPP of convection–diffusion type with implicit-Euler method as time integrators, whereas a higher-order numerical method is proposed in [19]. An ADI numerical scheme for time-dependent SPP of degenerate type is studied in [20]. Again, for system of PDE in two-dimension, very few articles are available in the literature: In [21], [22], the classical central difference scheme is considered to solve singularly perturbed system of elliptic BVPs. Recently, the authors proposed numerical methods for time dependent SPSCDP involving single diffusion parameters in [23]. Numerical simulations for time dependent SPSCDP and SPSRDP involving single parameters have been analyzed in [24], [25]. Therefore, it is always interesting to study singularly perturbed systems in more generalized sense, i.e., model problems with diffusion parameters of different magnitudes.

The primary intent of this study is to develop a robust finite difference method and analyze the error estimates for 2D parabolic SPSRDP involving distinct diffusion parameters. Prior knowledge about the behavior of the exact solution is required to obtain the error estimates. Because of the multi-scale nature of the model problem, it is difficult to extract the analytical properties of the exact solution. For this reason, here, we use the Shishkin-type decomposition of the exact solution and obtain the bounds for the derivatives of the solution for regular and layer components. Initially, we discretize the spatial derivatives by employing the central difference scheme on a layer-adapted Shishkin mesh. Later on, we approximate the resulting stiff initial-value problems by the implicit-Euler scheme. We derive the error estimate, which is almost second-order accurate in the spatial variables and first-order in the temporal variable.

This article is presented as follows: In Section 2, we discuss some results concerning the analytic behavior of the exact solution and its derivatives and later on a Shishkin-type decomposition has been established. In Section 3, the spatial semidiscretization has been introduced and we provide its truncation error analysis. For the time discretization, the implicit-Euler scheme has been used and the convergence result of the fully discrete scheme has been established in Section 4. In Section 5, we carried out some numerical experiments to show the accuracy and efficiency of the proposed method, and validate the theoretical error estimates. Finally, we concluded our work in Section 6.

Throughout this article, $\overrightarrow{C}={(C,C)}^T,C$ denotes a generic positive constant, which is independent of ${\epsilon }_1,{\epsilon }_2$ and $N,\Delta t$. We denote the partial derivative of $u$ with respect to $x$ by $u_x$, whereas $u^{\prime }$ denotes the time derivative.

Next, we introduce some notations for barrier functions. Let us define the following notations

$$\left\{\begin{array}{c}{\mathcal{B}}_{{\epsilon }_1}\left(x\right)=e^{-\beta x∕{\epsilon }_1}+e^{-\beta (1-x)∕{\epsilon }_1},{\mathcal{B}}_{{\epsilon }_1}\left(y\right)=e^{-\beta y∕{\epsilon }_1}+e^{-\beta (1-y)∕{\epsilon }_1},\hfill \\ {\mathcal{B}}_{{\epsilon }_2}\left(x\right)=e^{-\beta x∕{\epsilon }_2}+e^{-\beta (1-x)∕{\epsilon }_2},{\mathcal{B}}_{{\epsilon }_2}\left(y\right)=e^{-\beta y∕{\epsilon }_2}+e^{-\beta (1-y)∕{\epsilon }_2},\hfill \end{array}\right.$$

where $\beta$ is defined in (1.3) also assume that ${\epsilon }_2\le \beta ∕4$, from a practical point of view, this is always true.

For a vector-valued function $\overrightarrow{u}≔\left(u_1,u_2\right)$, set $\left|\overrightarrow{u}\right|={\left(\left|u_1\right|,\left|u_2\right|\right)}^T$ and define the maximum norm as ${\Vert \overrightarrow{u}\Vert }_{\infty }=max\left\{{\Vert u_1\Vert }_{\infty },{\Vert u_2\Vert }_{\infty }\right\}$. And define the maximum norm on ${\overline{\mathcal{G}}}^{N,M}≔{\{x_i,y_j,t_n\}}_{i,j,n=0}^{N,N,M}$, for a mesh function $U≔{\left\{U(x_i,y_j,t_n)\right\}}_{i,j,n=0}^{N,N,M}$ by ${\Vert U\Vert }_{\infty }={max}_{{\overline{\mathcal{G}}}^{N,M}}\left|U(x_i,y_j,t_n)\right|.$ For the mesh function $\overrightarrow{U}≔\left(U_1,U_2\right)$, define the maximum norm as ${\Vert \overrightarrow{U}\Vert }_{\infty }=max\left\{{\Vert U_1\Vert }_{\infty },{\Vert U_2\Vert }_{\infty }\right\}$.