Novel Numerical Scheme for Singularly Perturbed Time Delay Convection-Diffusion Equation

Abstract

This paper deals with numerical treatment of singularly perturbed parabolic differential equations having large time delay. The highest order derivative term in the equation is multiplied by a perturbation parameter $\epsilon$, taking arbitrary value in the interval $0,1$. For small values of $\epsilon$, solution of the problem exhibits an exponential boundary layer on the right side of the spatial domain. The properties and bounds of the solution and its derivatives are discussed. The considered singularly perturbed time delay problem is solved using the Crank-Nicolson method in temporal discretization and exponentially fitted operator finite difference method in spatial discretization. The stability of the scheme is investigated and analysed using comparison principle and solution bound. The uniform convergence of the scheme is discussed and proven. The formulated scheme converges uniformly with linear order of convergence. The theoretical analysis of the scheme is validated by considering numerical test examples for different values of $\epsilon$.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

A large number of mathematical models appear in different areas of science and engineering such as in control theory, epidemiology, and laser optics that take into account not only the present state of a physical system but also its past history [1]. These models are described by certain classes of functional differential equations often called delay differential equations. Examples of delays include the time taken for a signal to travel to the controlled object, driver reaction time, the time for the body to produce red blood cells, and cell division time in the dynamics of viral exhaustion or persistence. In the life sciences, delays are often introduced to account for hidden variables and processes which, although not well understood, are known to cause a time lag [1]. Time delays are natural components of the dynamic processes of biology, ecology, physiology, economics, epidemiology, and mechanics [2], and to ignore them is to ignore reality [1].

Singularly perturbed delay differential equations are differential equations in which its highest order derivative term is multiplied by small perturbation parameter $\epsilon$ and having at least one delay term. A simplified mathematical model for the control system of a furnace used to process metal sheets is represented by singularly perturbed parabolic delay differential equation of the form [3] $$\begin{array}{}\text{(1)}& \frac{\partial u}{\partial t}-\epsilon \frac{{\partial }^{2}u}{\partial x^2}=vgux,t-\tau \frac{\partial u}{\partial x}+cfux,t-\tau -ux,t,\end{array}$$where $u$ is the temperature distribution in a metal sheet, moving at an instantaneous material strip velocity $v$ and heated by distributed temperature source specified by the function $f$; both $v$ and $f$ are dynamically adapted by a controlling device monitoring the current temperature distribution. The finite speed of the controller, however, introduces a fixed delay of length $\tau$. When $\tau =0$, this problem becomes a thermal problem without time delay.

The presence of the singular perturbation parameter $\epsilon$ on the highest order derivative term leads to occurrences of oscillations in the computed solution while using central FDM, Galerkin FEM, and collocation methods on a uniform mesh [4]. While using these methods to avoid the oscillations, an unacceptably large number of mesh points are required when ε is very small. This is not practical due to the limited size of computer memory and rounding off error. Therefore, to overcome this drawback associated with classical numerical methods, authors use the fitted mesh technique which have a fine mesh in the boundary layer region or fitted operator technique. Recently, numerical treatment of singularly perturbed parabolic PDEs has gotten great attention, and different authors in [5–8] have developed uniformly convergent numerical schemes (scheme that converges independent of the influence of the singular perturbation parameter). Currently, authors in [9–13] have developed numerical schemes for solving singularly perturbed parabolic time delay reaction-diffusion equations. In papers [14, 15], the authors considered a singularly perturbed parabolic convection diffusion equation with time delay and degenerate coefficient. They treat the problems using fitted mesh techniques. Different authors in [16–23] developed numerical schemes using fitted mesh techniques for treating singularly perturbed parabolic time delay convection-diffusion equations. In [24], Podila and Kumar used nonstandard FDM for treating singularly perturbed parabolic time delay convection-diffusion equations. Their scheme gives linear order of convergence.

The main motive of this paper is to formulate an accurate and uniformly convergent numerical scheme for the singularly perturbed parabolic time delay convection-diffusion-reaction equation using exponentially fitted operator FDM and to establish the stability and uniform convergence of the scheme. The proposed scheme uses the procedure of Roth (i.e., first discretizing in temporal direction followed by discretization in spatial direction) using the Crank-Nicolson method in temporal direction and exponentially fitted operator FDM on spatial direction. In this method, it is not required to have any restriction on the mesh generation.

Notation 1.

Throughout this paper, $N$ and $M$ are denoted for the number of mesh intervals in space and time direction discretization, respectively; $C$ is denoted for positive constant independent of $\epsilon ,N,\text{and}M$. The norm $.$ denotes the maximum norm defined as $g=\underset{x,t}{\max }gx,t$.

2. Continuous Problem

Consider a class of second order singularly perturbed parabolic convection-diffusion-reaction equations having a term with large time delay, on the domain $D={\Omega }_x\times {\Omega }_t=0,1\times 0,T$ and on smooth boundary $\partial D=D_l\cup D_b\cup D_r$ is given by $$\begin{array}{}\text{(2)}& \frac{\partial u}{\partial t}-\epsilon \frac{{\partial }^{2}u}{\partial x^2}+ax\frac{\partial u}{\partial x}+bx,tux,t=-cx,tux,t-\tau +fx,t,\end{array}$$with the boundary conditions $$\begin{array}{}\text{(3)}& u0,t={\varphi }_{l}t,t\in {\Omega }_t,u1,t={\varphi }_{r}t,t\in {\Omega }_t,\end{array}$$and the interval condition $$\begin{array}{}\text{(4)}& ux,t={\varphi }_{b}x,t,x,t\in 0,1\times -\tau ,0,\end{array}$$where $\epsilon \in 0,1$ is a singular perturbation parameter and $\tau >0$ is the delay parameter. The terminal time $T$ is assumed to satisfy the condition $T=k\tau$ for some positive integer $k$. The functions $ax,bx,t,cx,t,\text{and}fx,t\text{and}{\varphi }_{l}t,{\varphi }_{b}x,t,\text{and}{\varphi }_{r}t$ are assumed to be sufficiently smooth and bounded and satisfy $$\begin{array}{}\text{(5)}& ax\ge \alpha >0,bx,t\ge \varsigma >0,cx,t\ge \gamma >0,x,t\in \overline{D}.\end{array}$$

This condition ensures that the solution of the problem in (2)–(4) exhibits a boundary layer of thickness $O\epsilon$ on the right side of the spatial domain.

Our objective in this paper is to formulate an accurate and parameter uniformly convergent numerical scheme and to discuss the uniform stability and the parameter uniform convergence of the scheme for the considered problem in (2)–(4).

2.1. Bounds on the Solution and Its Derivatives

The existence and uniqueness of the solution of (2)–(4) can be established by assuming that the data is Holder continuous and imposing appropriate compatibility conditions at the corner points, using the assumptions of sufficiently smoothness of${\varphi }_{l}t,{\varphi }_{b}x,t$, and${\varphi }_{r}t$.

The required compatibility condition at the corner points and the delay term are $$\begin{array}{}\text{(6)}& {\varphi }_{l}0={\varphi }_{b}0,0,{\varphi }_{r}0={\varphi }_{b}1,0,\text{(7)}& \begin{array}{c}\frac{d{\varphi }_{l}0}{dt}-\epsilon \frac{{\partial }^2{\varphi }_{b}0,0}{\partial x^2}=a0\frac{\partial {\varphi }_{b}0,0}{\partial x}+b0,0{\varphi }_{b}0,0=-c0,0u0,-\tau +f0,0,\\ \frac{d{\varphi }_{r}0}{dt}-\epsilon \frac{{\partial }^2{\varphi }_{b}1,0}{\partial x^2}=a1\frac{\partial {\varphi }_{b}1,0}{\partial x}+b1,0{\varphi }_{b}1,0=-c0,0u1,-\tau +f1,0,\end{array}\end{array}$$so that the data matches at the corner points. For $ax,bx,t,cx,t$, and $fx,t$ to be continuous on domain $D$, then (2)–(4) has the unique solution $ux,t\in C^{2}D\cap C^0\overline{D}$. In particular, when the compatibility conditions are not satisfied, a unique classical solution still exists but is not differentiable on all of $\partial D$.

The reduced problem corresponding to the singularly perturbed parabolic delay PDE (2)–(4) is given as $$\begin{array}{}\text{(8)}& \frac{\partial u^0}{\partial t}+ax\frac{\partial u^0}{\partial x}+bx,t{u}^{0}x,t=-cx,t{u}^{0}x,t-\tau +fx,t,x,t\in D,\end{array}$$with the boundary conditions $$\begin{array}{}\text{(9)}& u^{0}0,t={\varphi }_{l}t,t\in {\Omega }_t,u^{0}1,t\ne {\varphi }_{r}t,t\in {\Omega }_t,\end{array}$$and the interval condition $$\begin{array}{}\text{(10)}& u^{0}x,t={\varphi }_{b}x,t,x,t\in 0,1\times -\tau ,0.\end{array}$$

It is in the form of hyperbolic delay PDEs. The solution $ux,t$ of (2)–(4) becomes very close to the solution $u^{0}x,t$ of (8)–(10) as $\epsilon ⟶0$.

Lemma 2.

The solution $ux,t$ of (2)–(4) satisfies the estimate $$\begin{array}{}\text{(11)}& ux,t-{\varphi }_{b}x,0\le Ct,\hspace{1em}x,t\in \overline{D},\end{array}$$where $C$ is a constant independent of $\epsilon$.

Proof.

The result follows from the compatibility condition. See the detailed proof in [16].

Let $L$ be a differential operator that denotes the differential equation in (2)–(4) as $L=\partial /\partial t-\epsilon {\partial }^2/\partial x^2+ax\partial /\partial x+bx,t$.

Lemma 3.

Suppose the function $zx,t\in C^{2}D\cap C^0\overline{D}$, which satisfies $zx,t\ge 0,$$x,t\in \partial D$, and $Lzx,t\ge 0,x,t\in D$ implies that $zx,t\ge 0,x,t\in \overline{D}$.

Proof.

Assume there exists $x^{\ast },t^{\ast }\in \overline{D}$, such that $$\begin{array}{}\text{(12)}& z{x}^{\ast },t^{\ast }=\underset{x,t\in \overline{D}}{\min }zx,t<0.\end{array}$$

It is clear that the point $x^{\ast },t^{\ast }\notin \partial D,$ which implies that $x^{\ast },t^{\ast }\in D$. Since $z{x}^{\ast },t^{\ast }={\min }_{x,t\in \overline{D}}zx,t$ from extreme values in calculus which implies $\partial /\partial xz{x}^{\ast },t^{\ast }=0,$$\partial /\partial tz{x}^{\ast },t^{\ast }=0,\text{and}{\partial }^2/\partial x^{2}z{x}^{\ast },t^{\ast }\ge 0$ giving that $Lz{x}^{\ast },t^{\ast }<0,$ which is contradiction to the assumption that made $Lz{x}^{\ast },t^{\ast }\ge 0,\forall x,t\in D$. Therefore, $zx,t\ge 0,\forall x,t\in \overline{D}$.

Lemma 4.

Let $ux,t$ be the solution of the continuous problem in (2)–(4). Then, we obtain the bound $$\begin{array}{}\text{(13)}& ux,t\le {\zeta }^{-1}f+\max {\varphi }_{l}t,{\varphi }_{b}x,t,{\varphi }_{r}t,\end{array}$$where $\zeta$ is the lower bound of $bx,t$.

Proof.

By defining the barrier functions ${\vartheta }^{\pm }x,t$ as ${\vartheta }^{\pm }x,t={\zeta }^{-1}f+\max {\varphi }_{l}t,{\varphi }_{b}x,t,{\varphi }_{r}t\pm ux,t$ and applying the maximum principle, we obtain the required bound.

Lemma 5.

The bound on the derivative of the solution $ux,t$ of the problem in (2)–(4) with respect to $x$ and $t$ is given by $$\begin{array}{c}\text{(14)}& \frac{{\partial }^{k}ux,t}{\partial x^k}\le C1+{\epsilon }^{-k}\exp -\frac{\alpha 1-x}{\epsilon },\hspace{1em}x,t\in \overline{D},0\le k\le 4,\frac{{\partial }^{l}ux,t}{\partial t^l}\le C,\hspace{1em}x,t\in \overline{D},l=0,1,2,,\end{array}$$where $\alpha$ is lower bound of $ax$.

Proof.

See on [5, 16].

3. Numerical Scheme

In general, for singularly perturbed problems, there are two strategies for designing numerical methods which have a small error in the boundary layer region [25]. The first approach is the class of fitted mesh methods which uses fine mesh in the boundary layer region and coarse mesh in outer layer region. The stability and convergence analyses of this approach are well developed. The second approach is the fitted operator methods in which it uses uniform mesh and an exponentially fitting factor for stabilizing the term containing the singular perturbation parameter. In this approach, the difference schemes reflect the qualitative behaviour of the solution inside the boundary layer region. In this article, we formulate an exponentially fitted operator finite difference scheme to solve the problem in (2)–(4).

3.1. Temporal Semidiscretization

The time domain $0,T$ is discretized using a uniform mesh with time step $\Delta t$ as ${\Omega }_t^M=t_j=j\Delta t,j=0,1,\cdots ,M,t_M=T,\Delta t=T/M$ and${\Omega }_t^m=t_j=j\Delta t,j=0,1,\cdots ,m,t_m=\tau ,\Delta t=\tau /m$ where $M$ is the number of mesh points in time direction in the interval $0,T$ and $m$ is the number of mesh points in $-\tau ,0$. Note that $T=k\tau$ for some positive integer $k$.

For approximating the temporal derivative term of (2)–(4), we use the averaged Crank-Nicolson method, which gives a system of boundary value problems $$\begin{array}{}\text{(15)}& 1+\frac{\Delta t}{2}{L}^{\Delta t}{U}_{j+1}x=\begin{array}{c}1-\frac{\Delta t}{2}{L}^{\Delta t}{U}_{j}x-cx,t_{j+1/2}{\varphi }_{b}x+fx,t_{j+1/2},\hfill \\ \text{for}j=0,1,\cdots ,m,x\in {\Omega }_x\hfill \\ 1-\frac{\Delta t}{2}{L}^{\Delta t}{U}_{j}x-cx,t_{j+1/2}{U}_{j+1/2-m}x+fx,t_{j+1/2},\hfill \\ \text{for}j=m+1,\cdots ,M-1,x\in {\Omega }_x,\hfill \end{array}\end{array}$$with the discretized boundary conditions $$\begin{array}{}\text{(16)}& U_{j+1}0={\varphi }_l{t}_{j+1},U_{j+1}1={\varphi }_r{t}_{j+1},\hspace{1em}j=0,1,2,\cdots ,M,\end{array}$$where $$\begin{array}{}\text{(17)}& L^{\Delta t}{U}_{j+1}x=-\epsilon \frac{d^2}{{dx}^2}{U}_{j+1}x+ax\frac{d}{dx}{U}_{j+1}x+bx,t_{j+1/2}{U}_{j+1}x.\end{array}$$

Here, $U_{j+1}x$ is denoted for the approximation of $ux,t_{j+1}$ at the $j+1\text{th}$ time level and$bx,t_{j+1/2}=1/2bx,t_j+bx,t_{j+1}$, similarly for $cx,t_{j+1/2}\text{and}fx,t_{j+1/2}$.

The semidiscrete difference operator $1+\Delta t/2L^{\Delta t}{U}_{j+1}x$ in (15)–(16) satisfies the maximum principle as follows.

Lemma 6.

Semidiscrete maximum principle. Let $Z_{j+1}x$ be a sufficiently smooth function on ${\overline{\Omega }}_x$. If $Z_{j+1}0\ge 0,Z_{j+1}1\ge 0\text{and}1+\Delta t/2L^{\Delta t}{Z}_{j+1}x\ge 0,\forall x\in {\Omega }_x$, then $Z_{j+1}x\ge 0,\forall x\in {\overline{\Omega }}_x$.

Proof.

Assume that there exist $x^{\ast }\in {\overline{\Omega }}_x=0,1$ such that $Z_{j+1}{x}^{\ast }={\min }_{x\in {\overline{\Omega }}_x}{Z}_{j+1}x<0$. From the above assumption, it is clear that $x^{\ast }\notin 0,1\text{implies}\text{that}{x}^{\ast }\in 0,1$. Since $Z_{j+1}{x}^{\ast }={\min }_{x\in {\overline{\Omega }}_x}{Z}_{j+1}x$, using the property in calculus, we have $d/dx{Z}_{j+1}{x}^{\ast }=0$ and $d^2/{dx}^2{Z}_{j+1}{x}^{\ast }\ge 0$, then we obtain $1+\Delta t/2L^{\Delta t}{Z}_{j+1}{x}^{\ast }<0$ which is in contradiction to $1+\Delta t/2L^{\Delta t}{Z}_{j+1}{x}^{\ast }\ge 0,\forall x\in {\Omega }_x$. Therefore, we conclude $Z_{j+1}x\ge 0,\forall x\in {\overline{\Omega }}_x.$.

Next, let us analyse the truncation error for the temporal discretization made above.

Let the local error at each time step be denoted by $e_{j+1}x≔ux,t_{j+1}-U_{j+1}x,j=0,1,2,\cdots ,M$.

Lemma 7.

Suppose that $$\begin{array}{}\text{(18)}& \frac{{\partial }^{k}ux,t}{\partial t^k}\le C,x,t\in \overline{D},\hspace{1em}k=0,1,2.\end{array}$$

The local truncation error in the temporal direction is given by $$\begin{array}{}\text{(19)}& e_{j+1}\le C_1{\Delta t}^3.\end{array}$$

Proof.

Using Taylor’s series approximation for $ux,t_j$ and $ux,t_{j+1}$ centering at ${\text{t}}_{j+1/2},$ we obtain $$\begin{array}{}\text{(20)}& \begin{array}{c}ux,t_j=ux,t_{j+1/2}-\frac{\Delta t}{2}{u}_{t}x,t_{j+1/2}+\frac{{\Delta t}^2}{8}{u}_{tt}x,t_{j+1/2}+O{\Delta t}^3,\\ ux,t_{j+1}=ux,t_{j+1/2}+\frac{\Delta t}{2}{u}_{t}x,t_{j+1/2}+\frac{{\Delta t}^2}{8}{u}_{tt}x,t_{j+1/2}+O{\Delta t}^3.\end{array}\end{array}$$

From the approximation in (20), we obtain $$\begin{array}{}\text{(21)}& \frac{ux,t_{j+1}-ux,t_j}{\Delta t}=u_{t}x,t_{j+1/2}+O{\Delta t}^2.\end{array}$$

Using the approximation into (2), we obtain $$\begin{array}{}\text{(22)}& \frac{ux,t_{j+1}-ux,t_j}{\Delta t}=\epsilon u_{xx}x,t_{j+1/2}-ax{u}_{x}x,t_{j+1/2}-bx,t_{j+1/2}ux,t_{j+1/2}+fx,t_{j+1/2}+O{\Delta t}^2+\begin{array}{c}-cx,t_{j+1/2}{\varphi }_{b}x,\hfill \\ \text{for}j=0,1,\cdots ,m,\hfill \\ -cx,t_{j+1/2}ux,t_{j+1/2-m},\hfill \\ \text{for}j=m+1,\cdots ,M-1,\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(23)}& ux,t_{j+1/2}=\frac{ux,t_{j+1}+ux,t_j}{2}+O{\Delta t}^2,fx,t_{j+1/2}=\frac{fx,t_{j+1}+fx,t_j}{2}+O{\Delta t}^2.\end{array}$$

Since the error $e_{j+1}x≕ux,t_{j+1}-U_{j+1}x$ satisfies the semidiscrete differential operator $$\begin{array}{}\text{(24)}& \begin{array}{c}1+\frac{\Delta t}{2}{L}^{\Delta t}{e}_{j+1}x=O{\Delta t}^3,\\ e_{j+1}0=0=e_{j+1}1.\end{array}\end{array}$$

Hence, by applying the maximum principle, we obtain $$\begin{array}{}\text{(25)}& e_{j+1}\le C_1{\Delta t}^3.\end{array}$$

Next, we need to show the bound for the global error of the temporal discretization. Let us denote ${\text{TE}}_{j+1}$ as the global error up to the $j+1\text{th}$ time step.

Lemma 8.

The global error at $t_{j+1}$ is given by $$\begin{array}{}\text{(26)}& \text{TE}\le C{\Delta t}^2,j=1,2,\cdots ,M-1.\end{array}$$

Proof.

Using the local error up to the $j+1\text{th}$ time step given in the above lemma, we obtain the global error at the $j+1\text{th}$ time step as $$\begin{array}{}\text{(27)}& \text{TE}=\sum _{l=1}^{j+1}{e}_l\le e_1+e_2+\cdots +e_{j+1}\le C_{1}T{\Delta t}^2,\text{since}j+1\Delta t\le T=C{\Delta t}^2,C_{1}T=C,\end{array}$$where $C$ is a constant independent of $\epsilon \text{and}\Delta t$. See the detailed proof in [26].

Next, we set a bound for the derivatives of solution of (15)–(16).

Lemma 9.

For each $j=0,1,\cdots ,M-1$, the solution $U_{j+1}x$ of the boundary value problems in (15)–(16) satisfies the bound $$\begin{array}{}\text{(28)}& \frac{d^k{U}_{j+1}x}{{dx}^k}\le C1+{\epsilon }^{-k}\exp -\frac{\alpha 1-x}{\epsilon },\hspace{1em}x\in {\overline{\Omega }}_x,0\le k\le 4.\end{array}$$

Proof.

See the proof in [5].

3.2. Spatial Discretization

The spatial domain $0,1$ is discretized into $N$ equal number of subintervals, each of length $h$. Let $0=x_0,x_N=1$, and $x_i=ih,i=0,1,2,\cdots ,N$, be the mesh points. For spatial discretization, we apply an exponentially fitted operator finite difference method which helps us to hinder the influence of the singular perturbation parameter.

First, let us find the exponential fitting factor for anonymous BVPs and then apply discretization in the spatial direction.

3.2.1. Computing the Exponential Fitting Factor

To develop the numerical method for (15)–(16), we use the technique designed in the theory of asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer occurs on the right side of the domain. From the theory of singular perturbation in [27], the zero order asymptotic solution of the singularly perturbed boundary value problems of the form $$\begin{array}{}\text{(29)}& \begin{array}{c}-\epsilon u^{″}x+ax{u}^{\prime }x+bxux=fx,x\in {\Omega }_x=0,1,\hfill \\ u0=a_l,u1=a_r,\hfill \end{array}\end{array}$$is given by $$\begin{array}{}\text{(30)}& ux=u_{0}x+\frac{a1}{ax}{a}_r-u_{0}1\exp -{\int }_x^1\frac{ax}{\epsilon }-\frac{bx}{ax}dx+O\epsilon .\end{array}$$

Using Taylor’s series expansion for $ax$ and $bx$ about $x=1$ and restriction to their first terms and simplifying gives $$\begin{array}{}\text{(31)}& ux=u_{0}x+a_r-u_{0}1\exp -\frac{a11-x}{\epsilon }+O\epsilon ,\end{array}$$where $u_0$ is the solution of the reduced problems (obtained by setting $\epsilon =0$) in (29). Considering $h$ is reasonably small and evaluating the result in (31) at $x_i$ gives $$\begin{array}{}\text{(32)}& uih=u_{0}0+a_r-u_{0}1\exp -a1\frac{1}{\epsilon }-i\rho .\end{array}$$

Similarly, we have $$\begin{array}{}\text{(33)}& \begin{array}{c}{u}_{i-1}=u_{0}0+a_r-u_{0}1\exp -a1\frac{1}{\epsilon }-i-1\rho ,\\ u_{i+1}=u_{0}0+a_r-u_{0}1\exp -a1\frac{1}{\epsilon }-i+1\rho ,\end{array}\end{array}$$where $\rho =h/\epsilon$.

Consider a uniform grid ${\overline{\Omega }}_x^N={x_i}_{i=0}^N$ and denote $h=x_{i+1}-x_i$. For any mesh function $v_i$, define the following difference operators: $$\begin{array}{}\text{(34)}& D^{+}{v}_i=\frac{v_{i+1}-v_i}{h},D^{-}{v}_i=\frac{v_i-v_{i-1}}{h},D^0{v}_i=\frac{v_{i+1}-v_{i-1}}{2h},D^{+}{D}^{-}{v}_i=\frac{v_{i+1}-2v_i+v_{i-1}}{h^2}.\end{array}$$

To handle the effect of the perturbation parameter, we multiply artificial viscosity (exponentially fitting factor $\sigma \rho$) on the diffusive part of the problem. Introducing an exponentially fitting factor $\sigma \rho$ in (29) and applying the central finite difference scheme gives $$\begin{array}{}\text{(35)}& -\epsilon \sigma \rho D^{+}{D}^{-}u{x}_i+a{x}_i{D}^{0}u{x}_i+b{x}_{i}u{x}_i=f{x}_i.\end{array}$$

Multiplying (35) by $h$ and considering $h$ is small and truncating the term $hf{x}_i-b{x}_{i}u{x}_i$ gives $$\begin{array}{}\text{(36)}& \frac{\sigma \rho }{\rho }{u}_{i-1}-2u_i+u_{i+1}+\frac{a{x}_i}{2}{u}_{i+1}-u_{i-1}=0.\end{array}$$

From (32), we have $$\begin{array}{}\text{(37)}& u_{i\pm 1}=u_{0}0+a_r-u_{0}1\exp -a1\frac{1}{\epsilon }-i\pm 1\rho .\end{array}$$

Substituting (37) and (33) into (36) and simplifying, the exponential fitting factor is obtained as $$\begin{array}{}\text{(38)}& \sigma \rho =\frac{\rho a{x}_i}{2}\coth \frac{\rho a1}{2}.\end{array}$$

3.2.2. The Discrete Scheme

Using the central finite difference method for the spatial discretization of (15)–(16) and applying the exponential fitting factor in (38), for $i=1,2,\cdots ,N-1$, the fully discrete scheme becomes $$\begin{array}{}\text{(39)}& 1+\frac{\Delta t}{2}{L}^{\Delta t,h}{U}_{i,j+1}=\begin{array}{c}1-\frac{\Delta t}{2}{L}^{\Delta t,h}{U}_{i,j}-\Delta tc{x}_i,t_{j+1/2}{\varphi }_b{x}_i+\Delta tf{x}_i,t_{j+1/2},\hfill \\ \text{for}j=0,1,\cdots ,m,\hfill \\ 1-\frac{\Delta t}{2}{L}^{\Delta t,h}{U}_{i,j}-\Delta tc{x}_i,t_{j+1/2}{U}_{i,j+1/2-m}+\Delta tf{x}_i,t_{j+1/2},\hfill \\ \text{for}j=m+1,\cdots ,M-1,\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(40)}& L^{\Delta t,h}{U}_{i,j+1}=-\sigma \rho \epsilon D^{+}{D}^{-}{U}_{i,j+1}+a{x}_i{D}^0{U}_{i,j+1}+b{x}_i,t_{j+1/2}{U}_{i,j+1}.\end{array}$$

In explicit form, the scheme is rewritten as $$\begin{array}{}\text{(41)}& r_i^{-}{U}_{i-1,j+1}+r_i^c{U}_{i,j+1}+r_i^{+}{U}_{i+1,j+1}=s_i^{-}{U}_{i-1,j}+s_i^c{U}_{i,j}+s_i^{+}{U}_{i+1,j}+\begin{array}{c}-\Delta tc{x}_i,t_{j+1/2}{\varphi }_b{x}_i+\Delta tf{x}_i,t_{j+1/2},\hfill \\ \text{for}j=0,1,\cdots ,m,\hfill \\ -\Delta tc{x}_i,t_{j+1/2}{U}_{i,j+1/2-m}+\Delta tf{x}_i,t_{j+1/2},\hfill \\ \text{for}j=m+1,\cdots ,M-1,\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(42)}& r_i^{-}=-\frac{\Delta t}{2}\frac{\epsilon \sigma \rho }{h^2}+\frac{a{x}_i}{2h},r_i^c=\Delta t\frac{\epsilon \sigma \rho }{h^2}+\frac{1}{2}b{x}_i,t_{j+1/2}+1,r_i^{+}=-\frac{\Delta t}{2}\frac{\epsilon \sigma \rho }{h^2}-\frac{a{x}_i}{2h},\text{(43)}& s_i^{-}=\frac{\Delta t}{2}\frac{\epsilon \sigma \rho }{h^2}+\frac{a{x}_i}{2h},s_i^c=1-\Delta t\frac{\epsilon \sigma \rho }{h^2}+\frac{1}{2}b{x}_i,t_{j+1/2},s_i^{+}=\frac{\Delta t}{2}\frac{\epsilon \sigma \rho }{h^2}-\frac{a{x}_i}{2h}.\end{array}$$

3.3. Stability and Uniform Convergence Analysis

First, we need to prove the discrete comparison principle for the discrete scheme in (39).

Lemma 10.

Discrete comparison principle. There exist a comparison function $Z_{i,j+1}$ such that $1+\Delta t/2L^{\Delta t,h}{U}_{i,j+1}\le 1+\Delta t/2L^{\Delta t,h}{Z}_{i,j+1},i=1,2,\cdots ,N-1$ and if $U_{0,j+1}\le Z_{0,j+1}$ and $U_{N,j+1}\le Z_{N,j+1},$ then $U_{i,j+1}\le Z_{i,j+1},i=0,1,2,\cdots ,N$.

Proof.

The matrix associated with operator $1+\Delta t/2L^{\Delta t,h}{U}_{i,j+1}$ is of size $N+1\times N+1$ with its entries for $i=1,2,\cdots ,N-1$ are $$\begin{array}{}\text{(44)}& \begin{array}{c}{a}_{i,i-1}=-\frac{\Delta t}{4h}a{x}_i\coth \frac{ha1}{2\epsilon }+1<0,\\ a_{i,i+1}=-\frac{\Delta t}{4h}a{x}_i\coth \frac{ha1}{2\epsilon }-1\le 0,\text{since}\coth \frac{ha1}{2\epsilon }\ge 1,\\ a_{i,i}=\frac{\Delta t}{2h}a{x}_i\coth \frac{ha1}{2\epsilon }+bx,t_{j+1/2}+1>0.\end{array}\end{array}$$

So, the coefficient matrix satisfies the property of $M$ matrix. So, the inverse matrix exists and it is nonnegative. This guarantees the existence and uniqueness of the discrete solution. See the detailed proof in [28].

Lemma 11.

Let $Z_{i,j+1}=1+x_i$, for, $0\le i\le N$. Then, there exist a positive constant $C$ such that $1+\Delta t/2L^{\Delta t,h}{Z}_{i,j+1}\ge C,\text{for}1\le i\le N-1$.

Proof.

The proof is a simple computation, enables one to give a bound, that is uniform in $\epsilon$ for the norm of the inverse of $1+\Delta t/2L^{\Delta t,h}$. $$\begin{array}{}\text{(45)}& 1+\frac{\Delta t}{2}{L}^{\Delta t,h}{Z}_{i,j+1}=1-\frac{\Delta t}{2}\sigma \rho \epsilon D^{+}{D}^{-}{Z}_{i,j+1}+\frac{\Delta t}{2}a{x}_i{D}^0{Z}_{i,j+1}+\frac{\Delta t}{2}b{x}_i,t_{j+1/2}{Z}_{i,j+1}=1-\frac{\Delta t}{2}\sigma \rho \epsilon D^{+}{D}^{-}1+x_i+\frac{\Delta t}{2}a{x}_i{D}^{0}1+x_i+\frac{\Delta t}{2}b{x}_i,t_{j+1/2}1+x_i=1+x_i-\frac{\Delta t}{2h^2}\sigma \rho \epsilon 1+x_{i-1}-21+x_i+1+x_{i+1}+\frac{\Delta t}{4h}a{x}_{i}1+x_{i+1}-1+x_{i-1}+\frac{\Delta t}{2}b{x}_i,t_{j+1/2}1+x_i=1+x_i+\frac{\Delta t}{2}a{x}_i+\frac{\Delta t}{2}b{x}_i,t_{j+1/2}1+x_i\ge C,\text{since}ax\ge \alpha >0,b{x}_i,t_{j+1/2}\ge \zeta >0.\end{array}$$

Lemma 12.

(uniform stability estimate). The solution $U_{i,j+1}$ of the discrete scheme in (39) satisfies the bound $$\begin{array}{}\text{(46)}& U_{i,j+1}\le \frac{1+\Delta t/2L^{\Delta t,h}{U}_{i,j+1}}{1+\Delta t\zeta /2}+C\max {\varphi }_l{t}_{j+1},{\varphi }_r{t}_{j+1},\end{array}$$

where $b{x}_i,t_{j+1/2}\ge \zeta >0$.

Proof.

Let us construct a barrier function ${\pi }_{i,j+1}^{\pm }$ as ${\pi }_{i,j+1}^{\pm }=1+\Delta t/2L^{\Delta t,h}{U}_{i,j+1}/1+\Delta t\zeta /2+C\max {\varphi }_l{t}_{j+1},{\varphi }_r{t}_{j+1}\pm U_{i,j+1}$. We can easily show that${\pi }_{0,j+1}^{\pm }\ge 0\text{and}{\pi }_{N,j+1}^{\pm }\ge 0$. Then, $$\begin{array}{}\text{(47)}& 1+\frac{\Delta t}{2}{L}^{\Delta t,h}{\pi }_{i,j+1}^{\pm }={\pi }_{i,j+1}^{\pm }+\frac{\Delta t}{2}-\epsilon \sigma \rho D^{+}{D}^{-}{\pi }_{i,j+1}^{\pm }+a{x}_i{D}^0{\pi }_{i,j+1}^{\pm }+b{x}_i,t_{j+1/2}{\pi }_{i,j+1}^{\pm }=1+\frac{\Delta t}{2}b{x}_i,t_{j+1/2}\frac{1+\Delta t/2L^{\Delta t,h}{U}_{i,j+1}}{1+\Delta t/2\zeta }+C\max {\varphi }_l{t}_{j+1},{\varphi }_r{t}_{j+1}\pm U_{i,j+1}\pm \frac{\Delta t}{2}-\epsilon \sigma \rho D^{+}{D}^{-}{U}_{i,j+1}+a{x}_i{D}^0{U}_{i,j+1}+b{x}_i,t_{j+1/2}{U}_{i,j+1}=1+\frac{\Delta t}{2}b{x}_i,t_{j+1/2}\frac{1+\Delta t/2L^{\Delta t,h}{U}_{i,j+1}}{1+\Delta t/2\zeta }+C\max {\varphi }_l{t}_{j+1},{\varphi }_r{t}_{j+1}\pm 1+\frac{\Delta t}{2}{L}^{\Delta t,h}{U}_{i,j+1}\ge 0.\end{array}$$

Using the discrete comparison principle, we obtain ${\pi }_{i,j+1}^{\pm }\ge 0,i=0,1,2,\cdots ,N.$

Lemma 13.

If $Z_{i,j+1}$ be any mesh function such that $Z_{0,j+1}=Z_{N,j+1}=0.$Then, $$\begin{array}{}\text{(48)}& Z_{i,j+1}\le \frac{1}{\zeta }{\max }_{1\le k\le N-1}{L}^{\Delta t,h}{Z}_{k,j+1}.\end{array}$$

Proof.

Consider two barrier functions of the form ${\pi }_{i,j+1}^{\pm }=P\pm Z_{i,j+1},$ where $P=1/\zeta {\max }_{1\le k\le N-1}{L}^{\Delta t,h}{Z}_{k,j+1}$. It is easily shown that ${\pi }_{0,j+1}^{\pm }\ge 0,{\pi }_{N,j+1}^{\pm }\ge 0$. Next, we show that $L^{\Delta t,h}{\pi }_{i,j+1}^{\pm }\ge 0$. $$\begin{array}{}\text{(49)}& L^{\Delta t,h}{\pi }_{i,j+1}^{\pm }=L^{h,\Delta t}P\pm Z_{i,j+1}=\pm L^{\Delta t,h}{Z}_{i,j+1}+\frac{b{x}_i,t_{j+1/2}}{\zeta }{\max }_{1\le k\le N-1}{L}^{\Delta t,h}{Z}_{k,j+1}\ge 0.\end{array}$$

Hence, using the discrete comparison principle gives $$\begin{array}{}\text{(50)}& Z_{i,j+1}\le \frac{1}{\zeta }{\max }_{1\le k\le N-1}{L}^{\Delta t,h}{Z}_{k,j+1}.\end{array}$$

Next, we consider the semidiscrete problem in (15)–(16) and the fully discrete scheme in (39) to find the truncation error of the spatial direction discretization.

Theorem 14.

Let the coefficient functions $ax$ and $bx,t_{j+1/2}$ in (15)–(16) be sufficiently smooth functions so that $U_{j+1}x\in C^{4}0,1$. Then, the computed solution $U_{i,j+1}$ of the problem in (39) satisfies the bound $$\begin{array}{}\text{(51)}& L^{\Delta t,h}{U}_{j+1}{x}_i-U_{i,j+1}\le \frac{{Ch}^2}{h+\epsilon }1+{\epsilon }^{-3}\exp -\frac{\alpha 1-\eta }{\epsilon },\end{array}$$where $\eta \in x_0,x_N$ such that $\exp -\alpha 1-\eta /\epsilon ={\max }_{x_i}\exp -\alpha 1-x_i/\epsilon$.

Proof.

The local truncation error in space discretization is given as $$\begin{array}{}\text{(52)}& L^{\Delta t,h}{U}_{j+1}{x}_i-U_{i,j+1}=-\epsilon \sigma \rho \frac{d^2}{{dx}^2}-D^{+}{D}^{-}{U}_{j+1}{x}_i+a{x}_i\frac{d}{dx}-D^0{U}_{j+1}{x}_i\le \epsilon a{x}_i\frac{\rho }{2}\coth a1\frac{\rho }{2}-1D^{+}{D}^{-}{U}_{j+1}{x}_i+\epsilon \frac{d^2}{{dx}^2}-D^{+}{D}^{-}{U}_{j+1}{x}_i+a{x}_i\frac{d}{dx}-D^0{U}_{j+1}{x}_i,\end{array}$$where $\sigma \rho =a{x}_i\rho /2\coth a1\rho /2\text{and}\rho =h/\epsilon$.

Now for $\rho >0$, $C_1$ and $C_2$ are constants and we have $\rho \coth \rho -1\le C_1{\rho }^2,\text{for}\rho \le 1$. For $\rho ⟶\infty$, since $\underset{\rho \to \infty }{\lim }\coth \rho =1$ gives $\rho \coth \rho -1\le C_1\rho$. In general, for all $\rho >0,$ we write $$\begin{array}{}\text{(53)}& C_1\frac{{\rho }^2}{\rho +1}\le \rho \coth \rho -1\le C_2\frac{{\rho }^2}{\rho +1},\end{array}$$giving that $$\begin{array}{}\text{(54)}& \epsilon a{x}_i\frac{\rho }{2}\coth a1\frac{\rho }{2}-1\le \epsilon \frac{{h/\epsilon }^2}{h/\epsilon +1}=\frac{h^2}{h+\epsilon }.\end{array}$$

Using Taylor series expansion, we obtain the bound for $$\begin{array}{}\text{(55)}& D^{+}{D}^{-}{U}_{j+1}{x}_i\le C{U}_{j+1}^{\prime \prime }\eta ,\frac{d}{dx}-D^0{U}_{j+1}{x}_i\le {Ch}^2{U}_{j+1}^{\prime \prime \prime }\eta ,\frac{d^2}{{dx}^2}-D^{+}{D}^{-}{U}_{j+1}{x}_i\le {Ch}^2{U}_{j+1}^4\eta ,\end{array}$$where $U_{j+1}^k\eta ={\max }_{\eta \in x_0,x_N}{U}_{j+1}^k{x}_i$.

Using the bounds for the differences of the derivatives in (54) and (55), we obtain

$L^{\Delta t,h}{U}_{j+1}{x}_i-U_{i,j+1}\le C{h}^2/h+\epsilon U_{j+1}^{\prime \prime }\eta +C\epsilon h^2{U}_{j+1}^4\eta +{Ch}^2{U}_{j+1}^{\prime \prime \prime }\eta \le {Ch}^2/h+\epsilon U_{j+1}^{\prime \prime }\eta +{Ch}^2\epsilon U_{j+1}^4\eta +U_{j+1}^{\prime \prime \prime }\eta .$Using the bounds for the derivatives of the solution in Lemma 9 gives

$L^{\Delta t,h}{U}_{j+1}{x}_i-U_{i,j+1}\le {Ch}^2/h+\epsilon 1+{\epsilon }^{-2}\exp -\alpha 1-\eta /\epsilon +{Ch}^2\epsilon +{\epsilon }^{-3}\exp -\alpha 1-\eta /\epsilon +1+{\epsilon }^{-3}\exp -\alpha 1-\eta /\epsilon .$ Since ${\epsilon }^{-3}\ge {\epsilon }^{-2}$, we obtain $$\begin{array}{}\text{(56)}& L^{\Delta t,h}{U}_{j+1}{x}_i-U_{i,j+1}\le \frac{{Ch}^2}{h+\epsilon }1+{\epsilon }^{-3}\exp -\frac{\alpha 1-\eta }{\epsilon },\end{array}$$which gives the required bound.

Lemma 15.

For a fixed number of mesh numbers $N$, and for $\epsilon ⟶0,$ it holds $$\begin{array}{}\text{(57)}& \underset{\epsilon \to 0}{\lim }{\max }_{1\le i\le N-1}\frac{\exp -\alpha x_i/\epsilon }{{\epsilon }^m}=0,\underset{\epsilon \to 0}{\lim }{\max }_{1\le i\le N-1}\frac{\exp -\alpha 1-x_i/\epsilon }{{\epsilon }^m}=0,\hspace{1em}m=1,2,3,\cdots ,\end{array}$$where $x_i=ih,\forall i=1,2,\cdots ,N-1.$.

Proof.

On the discrete domain ${x_i}_0^N$, for the interior grid points, we have $$\begin{array}{}\text{(58)}& \begin{array}{c}{\max }_{1\le i\le N-1}\frac{\exp -\alpha x_i/\epsilon }{{\epsilon }^m}\le \frac{\exp -\alpha x_1/\epsilon }{{\epsilon }^m}=\frac{\exp -\alpha h/\epsilon }{{\epsilon }^m},\\ {\max }_{1\le i\le N-1}\frac{\exp -\alpha 1-x_i/\epsilon }{{\epsilon }^m}\le \frac{\exp -\alpha 1-x_{N-1}/\epsilon }{{\epsilon }^m}=\frac{\exp -\alpha h/\epsilon }{{\epsilon }^m},\end{array}\end{array}$$

Since $x_1=h,1-x_{N-1}=h.$ Repeatedly using L’Hospital’s rule gives $$\begin{array}{}\text{(59)}& \underset{\epsilon \to 0}{\lim }\frac{\exp -\alpha h/\epsilon }{{\epsilon }^m}=\underset{\zeta =1/\epsilon \to \infty }{\lim }\frac{{\zeta }^m}{\exp \alpha h\zeta }=\underset{\zeta =1/\epsilon \to \infty }{\lim }\frac{m!}{{\alpha h}^m\exp \alpha h\zeta }=0.\end{array}$$

This completes the proof.

Theorem 16.

The numerical solution $U_{i,j+1}$ of the problem in (39) satisfies the following uniform error bound $$\begin{array}{}\text{(60)}& \underset{0<\epsilon \le 1}{\sup }{\max }_{0\le i\le N}{U}_{j+1}{x}_i-U_{i,j+1}\le Ch.\end{array}$$

Proof.

Using Lemma 15 into (56) gives $$\begin{array}{}\text{(61)}& L^{\Delta t,h}{U}_{j+1}{x}_i-U_{i,j+1}\le \frac{{Ch}^2}{h+\epsilon }.\end{array}$$

Hence, using the result in Lemma 13, we obtain $$\begin{array}{}\text{(62)}& U_{j+1}{x}_i-U_{i,j+1}\le \frac{{Ch}^2}{h+\epsilon }.\end{array}$$

Using the supremum over all $\epsilon \in 0,1$, we obtain $$\begin{array}{}\text{(63)}& \underset{0<\epsilon \le 1}{\sup }{\max }_{0\le i\le N}{U}_{j+1}{x}_i-U_{i,j+1}\le Ch.\end{array}$$

Remark 17.

For the case $\epsilon >h$, the developed scheme gives second order convergence. For small values of $\epsilon \ll h$, the scheme is first order uniformly convergent in spatial direction.

Theorem 18.

Let $u\text{and}U$ be the solution of (2)–(4) and (39), respectively. Then, the following uniform error bound holds: $$\begin{array}{}\text{(64)}& \underset{0<\epsilon \le 1}{\sup }u-U\le Ch+{\Delta t}^2.\end{array}$$

Proof.

The combination of temporal and spatial error bounds gives the required result.

4. Numerical Examples, Results, and Discussions

Here, we illustrate the proposed scheme using model examples. The exact solutions of the considered examples are not known. We investigate the theoretical results by performing experiments using the proposed scheme.

Example 19.

Consider singularly perturbed time delay parabolic PDEs: $$\begin{array}{}\text{(65)}& \frac{\partial u}{\partial t}-\epsilon \frac{{\partial }^{2}u}{\partial x^2}+2-x^2\frac{\partial u}{\partial x}+xux,t+ux,t-\tau =10t^2\exp -tx1-x,\end{array}$$with interval condition $ux,t=0,\text{on}x,t\in 0,1\times -1,0$ and the boundary conditions $$\begin{array}{}\text{(66)}& u0,t={\varphi }_{l}t=0,u1,t={\varphi }_{r}t=0,\hspace{1em}\text{for}t\in 0,2.\end{array}$$

Example 20.

Consider singularly perturbed time delay parabolic PDEs: $$\begin{array}{}\text{(67)}& \frac{\partial u}{\partial t}-\epsilon \frac{{\partial }^{2}u}{\partial x^2}+2-x^2\frac{\partial u}{\partial x}+x+1t+1ux,t+ux,t-\tau =10t^2\exp -tx1-x,\end{array}$$with interval condition $ux,t=0,\text{on}x,t\in 0,1\times -1,0$ and the boundary conditions $$\begin{array}{}\text{(68)}& u0,t={\varphi }_{l}t=0,u1,t={\varphi }_{r}t=0,\hspace{1em}\text{for}t\in 0,2.\end{array}$$

Since the exact solutions of the considered examples are not known, with the help of double mesh techniques, we compute the maximum absolute error $E_{\epsilon }^{N,M}$ of the scheme. Let $U_{i,j}^{N,M}$ be the computed solution of the problem using $N,M$ mesh points and let $U_{i,j}^{2N,2M}$ be the computed solution on double number of mesh points $2N,2M$ by including the midpoints $x_{i+1/2}=x_{i+1}+x_i/2$ and $t_{j+1/2}=t_{j+1}+t_j/2$ into the mesh points. We calculate the maximum absolute error as $$\begin{array}{}\text{(69)}& E_{\epsilon }^{N,M}={\max }_{i,j}{U}_{i,j}^{N,M}-U_{i,j}^{2N,2M},\end{array}$$and the $\epsilon$-uniform (parameter uniform) error estimate by using $$\begin{array}{}\text{(70)}& E^{N,M}={\max }_{\epsilon }{E}_{\epsilon }^{N,M}.\end{array}$$

We calculate the rate of convergence of the proposed scheme using $$\begin{array}{}\text{(71)}& r_{\epsilon }^{N,M}={\log }_2{E}_{\epsilon }^{N,M}-{\log }_2{E}_{\epsilon }^{2N,2M}.\end{array}$$and the $\epsilon$-uniform rate of convergence using $$\begin{array}{}\text{(72)}& r^{N,M}={\log }_2{E}^{N,M}-{\log }_2{E}^{2N,2M}.\end{array}$$

In Figures 1 and 2, the numerical solution of Examples 19 and 20 is given, respectively, for different values of perturbation parameter$\epsilon$. As we observe on these figures, a strong boundary layer is maintained on the right side of the spatial domain as $\epsilon$ goes small starting from $\epsilon =2^{-5}\text{to}{2}^{-20}$. In these figures, we observe that the proposed scheme approximates the solution without creating oscillations or divergence. Note that time delay does not have an effect on position and size of the boundary layer since the layer occurs on the spatial domain direction. If the delay was on the spatial variable, we may have an interior layer which is created because of the delay parameter. In Tables 1 and 2, the maximum absolute error of Examples 19 and 20 is given, respectively, for different values of perturbation parameter $\epsilon$ and mesh numbers $N,M$. As we observed the results in these tables for each number of mesh $N,M$ as $\epsilon$ becomes small, the maximum absolute error after showing growth becomes stable and identical. This indicates the stability and uniform convergence of the scheme. In the last two rows of these tables, the $\epsilon$-uniform error and $\epsilon$-uniform rate of convergence of the scheme are given. As one observes, the scheme gives first order uniform convergence. The second order convergence of the temporal discretization is depicted in Table 3. In Tables 4 and 5, comparison of the $\epsilon$-uniform error and $\epsilon$-uniform rate of convergence of the proposed scheme with results of some recently published papers is given. As we observed, the proposed scheme gives a more accurate result than the results in [16, 18–20] and [24]. The proposed scheme has a limitation for solving nonlinear singularly perturbed problems.

[figure omitted; refer to PDF]

Table 1

Example 19, maximum absolute errors of the scheme.

$\epsilon \downarrow$	$N=M=16$	32	64	128	256	512
1	1.0742$e$-04	6.6552$e$-05	4.4167$e$-05	2.7392$e$-05	1.5041$e$-05	7.8477$e$-06
$$2^{-2}$$	3.8959$e$-04	3.3429$e$-04	2.5726$e$-04	1.5142$e$-04	8.1424$e$-05	4.2142$e$-05
$$2^{-4}$$	1.6240$e$-03	5.6521$e$-04	4.3969$e$-04	2.6289$e$-04	1.4270$e$-04	7.4226$e$-05
$$2^{-6}$$	8.3743$e$-03	2.0823$e$-03	8.0261$e$-04	3.8986$e$-04	1.9220$e$-04	9.5980$e$-05
$$2^{-8}$$	1.0404$e$-02	4.8511$e$-03	1.9751$e$-03	7.6437$e$-04	3.0100$e$-04	1.2812$e$-04
$$2^{-10}$$	1.0409$e$-02	4.9874$e$-03	2.6147$e$-03	1.3838$e$-03	5.8980$e$-04	2.2292$e$-04
$$2^{-12}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3950$e$-04	3.6693$e$-04
$$2^{-14}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3985$e$-04	3.7720$e$-04
$$2^{-16}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3985$e$-04	3.7720$e$-04
$$2^{-18}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3985$e$-04	3.7720$e$-04
$$2^{-20}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3985$e$-04	3.7720$e$-04
$$E^{N,M}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3985$e$-04	3.7720$e$-04
$$r^{N,M}$$	1.0615	0.9310	0.8795	0.9425	0.9719	—

Table 2

Example 20, maximum absolute errors of the scheme.

$\epsilon \downarrow$	$N=16$	32	64	128	256	512
	$M=10$	20	40	80	160	320
1	9.0612$e$-05	7.8840$e$-05	5.0913$e$-05	2.8361$e$-05	1.4905$e$-05	7.6151$e$-06
$$2^{-2}$$	5.0971$e$-04	3.4843$e$-04	1.9886$e$-04	1.0569$e$-04	5.4425$e$-05	2.7604$e$-05
$$2^{-4}$$	1.3108$e$-03	6.8037$e$-04	3.4323$e$-04	1.7164$e$-04	8.5788$e$-05	4.2882$e$-05
$$2^{-6}$$	4.2800$e$-03	1.6597$e$-03	6.2601$e$-04	2.5497$e$-04	1.1234$e$-04	5.2343$e$-05
$$2^{-8}$$	5.4803$e$-03	3.2949$e$-03	1.4175$e$-03	5.1091$e$-04	1.8211$e$-04	7.0533$e$-05
$$2^{-10}$$	5.4832$e$-03	3.3921$e$-03	1.8283$e$-03	9.1599$e$-04	3.8114$e$-04	1.3511$e$-04
$$2^{-12}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7807$e$-04	2.3354$e$-04
$$2^{-14}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7832$e$-04	2.4067$e$-04
$$2^{-16}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7832$e$-04	2.4067$e$-04
$$2^{-18}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7832$e$-04	2.4067$e$-04
$$2^{-20}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7832$e$-04	2.4067$e$-04
$$E^{N,M}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7832$e$-04	2.4067$e$-04
$$r^{N,M}$$	0.6928	0.8910	0.9549	0.9803	0.9909	—

Table 3

Maximum absolute error and rate of convergence of the scheme.

Example 19				Example 20
$\epsilon \downarrow$	$N=16$	64	256	$N=16$	64	256
	$M=16$	32	64	$M=10$	20	40
$$2^{-10}$$	1.0409$e$-02	2.6147$e$-03	5.8980$e$-04	5.4832$e$-03	1.8283$e$-03	3.8004$e$-04
	1.9931	2.1483	—	1.5845	2.0621	—
$$2^{-12}$$	1.0409$e$-02	2.6159$e$-03	7.3950$e$-04	5.4832$e$-03	1.8292$e$-03	4.7807$e$-04
	1.9925	1.8227	—	1.5838	1.9359	—
$$2^{-14}$$	1.0409$e$-02	2.6159$e$-03	7.3985$e$-04	5.4832$e$-03	1.8292$e$-03	4.7832$e$-04
	1.9925	1.8220	—	1.5838	1.9352	—
$$2^{-16}$$	1.0409$e$-02	2.6159$e$-03	7.3985$e$-04	5.4832$e$-03	1.8292$e$-03	4.7832$e$-04
	1.9925	1.8220	—	1.5838	1.9352	—
$$2^{-18}$$	1.0409$e$-02	2.6159$e$-03	7.3985$e$-04	5.4832$e$-03	1.8292$e$-03	4.7832$e$-04
	1.9925	1.8220	—	1.5838	1.9352	—
$$2^{-20}$$	1.0409$e$-02	2.6159$e$-03	7.3985$e$-04	5.4832$e$-03	1.8292$e$-03	4.7832$e$-04
	1.9925	1.8220	—	1.5838	1.9352	—

Table 4

Comparison of uniform error and uniform rate of convergence of Example 19.

Schemes $\downarrow$		$N=M=16$	32	64	128	256
Proposed scheme	$$E^{N,M}$$	1.0409$e$-02	4.9874$e$-03	2.6159$e$-03	1.4219$e$-03	7.3985$e$-04
	$$r^{N,M}$$	1.0615	0.9310	0.8795	0.9425	0.9719
Scheme in [19]	$$E^{N,M}$$	3.41$e$-02	1.84$e$-02	9.38$e$-03	4.67$e$-03	2.31$e$-03
	$$r^{N,M}$$	0.8901	0.9720	1.0062	1.0155	1.0063
Scheme in [16]	$$E^{N,M}$$	4.9485$e$-02	3.3203$e$-02	2.1165$e$-02	1.3320$e$-02	7.9345$e$-03
	$$r^{N,M}$$	0.5757	0.6496	0.6681	0.7474	0.7908

Table 5

Comparison of uniform error and uniform rate of convergence of Example 20.

Schemes $\downarrow$		$N=16$	32	64	128	256
		$M=10$	20	40	80	160
Proposed scheme	$$E^{N,M}$$	5.4832$e$-03	3.3921$e$-03	1.8292$e$-03	9.4367$e$-04	4.7832$e$-04
	$$r^{N,M}$$	0.6928	0.8910	0.9549	0.9803	0.9909
Scheme in [20]	$$E^{N,M}$$	1.86$e$-02	1.00$e$-02	5.48$e$-03	2.86$e$-03	1.46$e$-03
	$$r^{N,M}$$	0.89	0.87	0.94	0.97	1.11
Scheme in [18]	$$E^{N,M}$$	1.6119$e$-02	9.9504$e$-03	5.8541$e$-03	3.3439$e$-03	1.8650$e$-03
	$$r^{N,M}$$	0.6960	0.7653	0.8079	0.8424	0.8660
Scheme in [24]	$$E^{N,M}$$	7.4252$e$-03	4.0993$e$-03	2.1528$e$-03	1.1033$e$-03	5.5845$e$-04
	$$r^{N,M}$$	0.8570	0.9291	0.9644	0.9822	—

5. Conclusion

In this paper, singularly perturbed parabolic convection-diffusion-reaction equation with large time delay is considered. The solution of the considered problem exhibits boundary layer of thickness $\epsilon$ on the right side of the spatial domain. The bounds and properties of the analytical solution are discussed. To stabilize the influence of the singular perturbation on the discrete solution, and exponential fitting parameter is developed using the derivation from the asymptotic solution. The numerical scheme is developed using the Crank-Nicolson method in temporal discretization and an exponentially fitted central finite difference method on the spatial discretization. The existence of the discrete solution is discussed using the comparison principle. The stability and uniform convergence of the scheme are investigated well theoretically. Numerical results are depicted using maximum absolute error, $\epsilon$-uniform error, and $\epsilon$-uniform rate of convergence in tables which are in good agreement with the theoretical analysis. The developed scheme gives stable and uniformly convergent result with linear order of convergence.

References

[1] C. T. Baker, G. A. Bocharov, F. A. Rihan, A Report on the Use of Delay Differential Equations in Numerical Modelling in the Biosciences, 1999.

[2] D. A. Goodwin, C. F. Meares, L. H. DeRiemer, C. I. Diamanti, R. L. Goode, J. J. Baumert, D. J. Sartoris, R. L. Lantieri, H. D. Fawcett, "Clinical studies with in-111 BLEDTA, a tumor-imaging conjugate of bleomycin with a bifunctional chelating agent," Journal of nuclear medicine, vol. 22 no. 9, pp. 787-792, 1981.

[3] P. Wang, C. Keng, "Asymptotic stability of a time-delayed diffusion system," Journal of Applied Mechanics, vol. 30 no. 4, pp. 500-504, 1963.

[4] M. M. Woldaregay, G. F. Duressa, "Higher-order uniformly convergent nu equations with mixed small shifts," International Journal of Differential Equations, vol. 2020,DOI: 10.1155/2020/6661592, 2020.

[5] C. Clavero, J. C. Jorge, F. Lisbona, "A uniformly convergent scheme on a non-uniform mesh for convection-diffusion parabolic problems," Journal of Computational and Applied Mathematics, vol. 154 no. 2, pp. 415-429, DOI: 10.1016/S0377-0427(02)00861-0, 2003.

[6] S. Gowrisankar, S. Natesan, "Robust numerical scheme for singularly perturbed convection-diffusion parabolic initial-boundary-value problems on equidistributed grids," Computer Physics Communications, vol. 185 no. 7, pp. 2008-2019, DOI: 10.1016/j.cpc.2014.04.004, 2014.

[7] M. K. Kadalbajoo, V. Gupta, A. Awasthi, "A uniformly convergent B-spline collocation method on a nonuniform mesh for singularly perturbed one-dimensional time-dependent linear convection-diffusion problem," Journal of Computational and Applied Mathematics, vol. 220 no. 1-2, pp. 271-289, DOI: 10.1016/j.cam.2007.08.016, 2008.

[8] D. A. Turuna, M. M. Woldaregay, G. F. Duressa, "Uniformly convergent numerical method for singularly perturbed convection-diffusion problems," Kyungpook National University, vol. 60 no. 3, pp. 629-645, 2020.

[9] A. R. Ansari, S. A. Bakr, G. I. Shishkin, "A parameter-robust finite difference method for singularly perturbed delay parabolic partial differential equations," Journal of Computational and Applied Mathematics, vol. 205 no. 1, pp. 552-566, DOI: 10.1016/j.cam.2006.05.032, 2007.

[10] E. B. Bashier, K. C. Patidar, "A novel fitted operator finite difference method for a singularly perturbed delay parabolic partial differential equation," Applied Mathematics and Computation, vol. 217 no. 9, pp. 4728-4739, DOI: 10.1016/j.amc.2010.11.028, 2011.

[11] E. B. Bashier, K. C. Patidar, "A second-order fitted operator finite difference method for a singularly perturbed delay parabolic partial differential equation," Journal of Difference Equations and Applications, vol. 17 no. 5, pp. 779-794, DOI: 10.1080/10236190903305450, 2011.

[12] G. Lolugu, J. Mohapatra, A. Das, "A fourth-order numerical scheme for singularly perturbed delay parabolic problem arising in population dynamics," Journal of Applied Mathematics and Computing, vol. 63, 2020.

[13] N. R. Rao, P. P. Chakravarthy, "A fitted Numerov method for singularly perturbed parabolic partial differential equation with a small negative shift arising in control theory," Numerical Mathematics: Theory, Methods and Applications, vol. 7 no. 1, pp. 23-40, 2014.

[14] S. Kumar, J. Vigo-Aguiar, "A parameter-uniform grid equidistribution method for singularly perturbed degenerate parabolic convection-diffusion problems," Journal of Computational and Applied Mathematics, vol. no. article 113273,DOI: 10.1016/j.cam.2020.113273, 2020.

[15] P. Rai, Y. Swati, "Robust numerical schemes for singularly perturbed delay parabolic convection-diffusion problems with degenerate coefficient," International Journal of Computer Mathematics, vol. 98, pp. 195-221, 2020.

[16] A. Das, S. Natesan, "Uniformly convergent hybrid numerical scheme for singularly perturbed delay parabolic convection-diffusion problems on Shishkin mesh," Applied Mathematics and Computation, vol. 271, pp. 168-186, DOI: 10.1016/j.amc.2015.08.137, 2015.

[17] L. Govindarao, J. Mohapatra, S. R. Sahu, "Uniformly convergent numerical method for singularly perturbed two parameter time delay parabolic problem," International Journal of Applied and Computational Mathematics, vol. 5 no. 3,DOI: 10.1007/s40819-019-0672-5, 2019.

[18] S. Gowrisankar, S. Natesan, "ε -Uniformly convergent numerical scheme for singularly perturbed delay parabolic partial differential equations," International Journal of Computer Mathematics, vol. 94 no. 5, pp. 902-921, 2017.

[19] D. Kumar, P. Kumari, "A parameter-uniform numerical scheme for the parabolic singularly perturbed initial boundary value problems with large time delay," Journal of Applied Mathematics and Computing, vol. 59 no. 1-2, pp. 179-206, DOI: 10.1007/s12190-018-1174-z, 2019.

[20] D. Kumar, P. Kumari, "A parameter-uniform scheme for singularly perturbed partial differential equations with a time lag," Numerical Methods for Partial Differential Equations, vol. 36 no. 4, pp. 868-886, DOI: 10.1002/num.22455, 2020.

[21] G. Lolugu, J. Mohapatra, "Numerical analysis and simulation of delay parabolic partial differential equation involving a small parameter," Engineering Computations, vol. 37 no. 1, pp. 289-312, DOI: 10.1108/EC-03-2019-0115, 2019.

[22] G. Lolugu, R. Sahu, J. Mohapatra, "Uniformly convergent numerical method for singularly perturbed time delay parabolic problem with two small parameters," Iranian Journal of Science and Technology, Transaction A, Science, vol. 43 no. 5, pp. 2373-2383, 2019.

[23] G. Lolugu, J. Mohapatra, "A second order numerical method for singularly perturbed delay parabolic partial differential equation," Engineering Computations, vol. 36 no. 2, pp. 420-444, DOI: 10.1108/EC-08-2018-0337, 2019.

[24] P. C. Podila, K. Kumar, "A new stable finite difference scheme and its convergence for time-delayed singularly perturbed parabolic PDEs," Computational & Applied Mathematics, vol. 39 no. 3,DOI: 10.1007/s40314-020-01170-2, 2020.

[25] J. J. Miller, E. O'riordan, G. I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions, 2012.

[26] M. K. Kadalbajoo, A. Awasthi, "Crank-Nicolson finite difference method based on a midpoint upwind scheme on a non-uniform mesh for time-dependent singularly perturbed convection-diffusion equations," International Journal of Computer Mathematics, vol. 85 no. 5, pp. 771-790, DOI: 10.1080/00207160701459672, 2008.

[27] R. E. OMalley, Singular perturbation methods for ordinary differential equations,DOI: 10.1007/978-1-4612-0977-5, 1991.

[28] T. Kellogg, A. Tsan, "Analysis of some difference approximations for a singular perturbation problem without turning points," Mathematics of Computation, vol. 32 no. 144, pp. 1025-1039, DOI: 10.1090/S0025-5718-1978-0483484-9, 1978.