An Exponentially Fitted Upwind Scheme for Singularly Perturbed Differential Equations With Mixed Shift Parameters

Abstract

This paper provides numerical solutions to a class of singularly perturbed differential–difference equations characterized by mixed shift parameters. The solutions of such problems exhibit sharp boundary layers near the endpoints of the spatial domain due to the presence of a small perturbation parameter ε(0 < ε ≪ 1). Consequently, classical numerical methods fail to give accurate results on uniform meshes. To address this challenge, we propose a numerical scheme that discretizes the problem using the Crank–Nicolson method in the temporal direction and an exponentially fitted finite difference scheme in the spatial direction, both on uniform meshes. Stability and convergence analyses confirmed that the proposed scheme is uniformly convergent with respect to the perturbation parameter ε, with second‐order accuracy in the temporal and spatial directions. Three model examples were presented for simulation, and the findings indicated that the theoretical analysis aligns with the practical results. Furthermore, the numerical results demonstrated that the proposed scheme outperforms several existing methods in the literature.

1. Introduction

Singularly perturbed differential-difference equations (SPDDEs) are a class of differential equations that involve delayed or advanced arguments of the unknown function, where the highest-order derivative is multiplied by a small parameter ε, (0 < ε ≪ 1), leading to multiscale behavior. Such types of equations frequently arise in a variety of mathematical models, including physiological kinetics [1], control theory [2], evolutionary biology [3], population dynamics and epidemiology [4], and neuronal variability [5]. Particularly, spatial delay differential equations form a subclass of SPDDEs that incorporate at least one retarded and/or advanced argument in the spatial variable. These types of equations are employed to model a wide range of real-world phenomena, including physiological processes [6], population dynamics [7], disease spread [8], polymer diffusion [9], epidemiology [10], and liquid helium hydrodynamics [11].

The presence of a small parameter ε as a diffusion coefficient in the proposed problem gives rise to boundary layer regions near the left or right endpoints of the spatial domain, depending on the sign of the convection coefficient. Approximate solutions within these boundary layer regions, when computed on uniform meshes, often exhibit oscillatory behavior and/or abrupt changes. Consequently, the classical numerical methods on uniform meshes are inefficient for accurately and effectively solving such problems unless the meshes are suitably small. Nonetheless, employing very fine meshes can lead to round-off errors and significantly increase computational costs [12]. This underscores the necessity for numerical methods that remain robust regardless of parameter variations [13]. Consequently, parameter-uniformly convergent algorithms have been developed [14–21]. Examples of such parameter-independent approaches include fitted operator methods, fitted mesh techniques, and domain decomposition strategies. There are also a series of scholarly articles which are higher-order uniformly convergent numerical schemes for solving singularly perturbed problems [22–26].

For the solution of parabolic-type singularly perturbed differential–difference equations involving small and/or large delays in the time direction, many scholars have developed and examined parameter-uniformly convergent numerical techniques. A few to note are as follows: Tefera et al. [27] devised multiple fitting operator schemes on the Crank–Nicolson to discretize the continuous problem domain. Woldaregay and Duressa [28] investigated the same problem using the Crank–Nicolson scheme in the temporal direction on a uniform mesh and the hybrid finite difference scheme in the spatial direction over a Shishkin mesh. Tiruneh et al. [29] designed a higher-order uniformly convergent scheme using a nonstandard finite difference method over Crank–Nicolson for domain discretization. Hassen and Duressa [30] developed an oscillation-free first-order parameter uniform exponentially spline method to solve a singularly perturbed time-dependent delay convection–diffusion problem with Dirichlet boundary conditions.

On the other hand, several researchers have established parameter-uniformly convergent numerical schemes for singularly perturbed parabolic-type convection–diffusion equations involving small and/or large general shifts in the spatial variable. For instance: Ejere et al. [31] devised and analyzed a second-order parameter-uniformly convergent scheme to treat singularly perturbed parabolic differential equation involving large spatial delay. Ramesh and Kadalbajoo [32] discretized the problem domain using implicit Euler and midpoint upwind finite difference schemes and proved the uniform convergence of the scheme. Kumar and Kadalbajoo [33] devised a parameter-uniformly convergent scheme to solve the proposed problem using the implicit Euler and B-spline collocation method in the temporal and spatial directions, respectively. Bansal and Sharma [34] formulated a parameter-uniformly convergent scheme for time-dependent singularly perturbed convection–diffusion–reaction problems with general shift arguments in the space variable. Gupta et al. [35] developed a higher-order parameter-uniformly convergent scheme using implicit Euler on a uniform mesh and a hybrid finite difference scheme on Shishkin mesh in the temporal and spatial directions, respectively. In a series of papers [36–38], Woldaregay and Duressa evolved and analyzed different types of parameter-uniformly convergent numerical schemes to solve a class of time-dependent singularly perturbed convection–diffusion equations with general shift arguments.

The main contribution of the present work is to formulate and analyze a higher-order, reliable, computationally simple and more accurate parameter-uniformly convergent numerical scheme to solve singularly perturbed differential–difference equations having mixed shift arguments in the spatial variable. We discretized the problem using the Crank–Nicolson approach and an exponentially fitted finite difference scheme in the temporal and spatial directions on uniform meshes, respectively.

The remaining parts of the paper are organized as follows. Section 2 discusses the problem statement and some prior estimates. Section 3 describes the domain semidiscretization. Section 4 examines stability and uniform convergence analysis. In Section 5, Richardson’s extrapolation is presented .In Section 6, numerical results and discussions are presented. Section 7 contains the work’s concluding notes.

Throughout this paper, C and its subscripts denote generic positive constants independent of the perturbation parameter ε and mesh sizes h and Δt. Also, ‖.‖ denotes the standard supremum norm, defined as $‖f‖={\sup }_{\left(x,t\right)\in {\overline{\mathcal{D}}}_{xt}}\left|f\left(x,t\right)\right|$ for a function f defined on some domain ${\overline{\mathcal{D}}}_{xt}$ with the boundary $\mathrm{\partial }{\mathcal{D}}_{xt}={\overline{\mathcal{D}}}_{xt}\setminus {\mathcal{D}}_{xt}$.

2. Statement of the Problem

In this work, we consider a class of singularly perturbed differential–difference equations having small shifts subject to the initial and interval boundary conditions of the form: 1 $\begin{array}{}\left\{\begin{array}{cc}{\mathcal{L}}_{\epsilon }\mathrm{\Psi }\left(x,t\right)=f\left(x,t\right),& \forall \left(x,t\right)\in {\mathcal{D}}_{xt}={\mathcal{D}}_x\times {\mathcal{D}}_t=\left(01,\right)\times \left(0,T\right],T>0,\\ \psi \left(x,t\right)={\mathcal{S}}_0\left(x,t\right),& \forall \left(x,t\right)\in {\mu }^0=\left\{\left(x,t\right):010\le x\le ,<t\le T\right\},\\ \psi \left(x,t\right)={\mathcal{S}}_l\left(x,t\right),& \forall \left(x,t\right)\in {\mu }_{-}=\left\{\left(x,t\right):-\delta \le x\le 00,<t\le T\right\},\\ \psi \left(x,t\right)={\mathcal{S}}_r\left(x,t\right),& \forall \left(x,t\right)\in {\mu }_{+}=\left\{\left(x,t\right):110\le x\le +\eta ,<t\le T\right\},\end{array}\right.\end{array}$ where ${\mathcal{L}}_{\epsilon }\mathrm{\Psi }≔{\mathrm{\Psi }}_t\left(x,t\right)-\epsilon {\mathrm{\Psi }}_{xx}\left(x,t\right)+a\left(x\right){\mathrm{\Psi }}_x\left(x,t\right)+b\left(x\right)\mathrm{\Psi }\left(x-\delta ,t\right)+c\left(x\right)\mathrm{\Psi }\left(x,t\right)+d\left(x\right)\mathrm{\Psi }\left(x+\eta ,t\right)$ is the differential operator, ε is the perturbation parameter (0 < ε ≪ 1), ${\mathcal{D}}_x$ is the spatial domain, ${\mathcal{D}}_{xt}$ is the space-time domain, and δ and η are the step behind and step ahead parameters with the property δ, η < ε, respectively. We assume that the functions a(x), b(x), c(x), and d(x), defined on ${\mathcal{D}}_x$, and $f\left(x,t\right),{\mathcal{S}}_0\left(x,t\right)$, ${\mathcal{S}}_l\left(x,t\right)$ and ${\mathcal{S}}_r\left(x,t\right)$, defined on ${\mathcal{D}}_{xt}$, are sufficiently smooth, bounded, and independent of the perturbation parameter ε. Furthermore, the coefficients of the reaction term c(x), the delay term b(x), and the advance term d(x) are assumed to satisfy the following condition: 2 $\begin{array}{}b\left(x\right)+c\left(x\right)+d\left(x\right)\ge \beta >0,\forall \left(x\right)\in {\overline{\mathcal{D}}}_x.\end{array}$

This condition ensures the well-posedness and stability of problem (1) by enforcing a positive combined reaction-shift effect, thereby preventing unbounded growth, ensuring uniqueness, and enabling the derivation of a priori bounds [20]. In addition, when the value of the perturbation parameter ε approaches zero, the solution of the proposed problem exhibits a left or right boundary layer in the neighborhood of μ₋ or μ₊ depending on the sign a(x) < 0 or (>0) However, due to its physical consistency with natural flow behavior and its numerical advantages, such as stability and compatibility, the right boundary layer configuration is adopted in the proposed problem.

As outlined in [32] for small shift parameters δ, η < ε, we employ the Taylor series expansion to approximate the shift terms as follows: 3 $\begin{array}{}\left\{\begin{array}{c}\psi \left(x-\delta ,t\right)=\psi -\delta {\psi }_x+\frac{{\delta }^2}{2}{\psi }_{xx}+O\left({\delta }^3\right),\\ \psi \left(x+\eta ,t\right)=\psi +\eta {\psi }_x+\frac{{\eta }^2}{2}{\psi }_{xx}+O\left({\eta }^3\right).\end{array}\right.\end{array}$

Substituting (3) into (1) yields a standard singularly perturbed initial–boundary value problem of the following form: 4 $\begin{array}{}\left\{\begin{array}{cc}\mathcal{L}\mathrm{\Psi }\left(x,t\right)≔{\mathrm{\Psi }}_t\left(x,t\right)+{\mathcal{L}}_{x,C_{\epsilon }}\mathrm{\Psi }\left(x,t\right)=f\left(x,t\right),& \forall \left(x,t\right)\in {\mathcal{D}}_{xt},\\ \psi \left(x,t\right)={\mathcal{S}}_0\left(x\right),& \forall \left(x,t\right)\in {\mu }_0=\left\{\left(x,t\right):x\in {\overline{\mathcal{D}}}_x,t=0\right\},\\ \psi \left(x,t\right)={\mathcal{S}}_l\left(0,t\right),& \forall \left(x,t\right)\in {\mu }_l=\left\{\left(x,t\right):x=0,t\in {\overline{\mathcal{D}}}_t\le T\right\},\\ \psi \left(x,t\right)={\mathcal{S}}_r\left(1,t\right),& \forall \left(x,t\right)\in {\mu }_r=\left\{\left(x,t\right):x=1,t\in {\overline{\mathcal{D}}}_t\le T\right\},\end{array}\right.\end{array}$ where ${\mathcal{L}}_{x,C_{\epsilon }}\mathrm{\Psi }≔-C_{\epsilon }{\mathrm{\Psi }}_{xx}+p\left(x\right){\mathrm{\Psi }}_x+q\left(x\right)\mathrm{\Psi }$, C_ε = ε − β₁(δ²/2) − β₂(η²/2), b(x) ≥ β₁ > 0, d(x) ≥ β₂ > 0 with 0 < C_ε ≪ 1, p(x) = a(x) − δb(x) + ηd(x),  and q(x) = b(x) + c(x) + d(x). The coefficient function q(x) is assumed to be bounded such that 5 $\begin{array}{}q\left(x\right)\ge q^{\mathrm{\ast }}>0,x\in {\overline{\mathcal{D}}}_x.\end{array}$

As presented in (2), this condition guarantees the well-posedness and stability of problem (4). Moreover, decreasing the values of the parameter C_ε leads to the formation of a right boundary layer with p(x) > 0 in the neighborhood of μ_r.

The existence and uniqueness of the solution of problem (4) are established under the assumption that the given data are Hölder continuous and satisfy the following compatibility conditions as described in [39]: 6 $\begin{array}{c}{\mathcal{S}}_l\left(00,\right)={\mathcal{S}}_0\left(0\right),{\mathcal{S}}_r\left(10,\right)={\mathcal{S}}_0\left(1\right),\end{array}$ and 7 $\begin{array}{}\left\{\begin{array}{c}\frac{d{\mathcal{S}}_l\left(00,\right)}{dt}-C_{\epsilon }\frac{d^2{\mathcal{S}}_0\left(0\right)}{d{x}^2}+p\left(0\right)\frac{d{\mathcal{S}}_0\left(0\right)}{dx}+q\left(0\right){\mathcal{S}}_0\left(0\right)=f\left(00,\right),\\ \frac{d{\mathcal{S}}_r\left(10,\right)}{dt}-C_{\epsilon }\frac{d^2{\mathcal{S}}_0\left(1\right)}{d{x}^2}+p\left(1\right)\frac{d{\mathcal{S}}_0\left(1\right)}{dx}+q\left(1\right){\mathcal{S}}_0\left(1\right)=f\left(10,\right).\end{array}\right.\end{array}$

The above compatibility conditions ensure the existence of a constant C independent of the perturbation parameter C_ε such that $\forall \left(x,t\right)\in {\overline{\mathcal{D}}}_{xt}$. 8 $\begin{array}{}\left\{\begin{array}{c}\left|\mathrm{\Psi }\left(x,t\right)-{\mathrm{\Psi }}_0\left(x,0\right)\right|=\left|\mathrm{\Psi }\left(x,t\right)-{\mathcal{S}}_0\left(x\right)\right|\le Ct,\\ \left|\mathrm{\Psi }\left(x,t\right)-{\mathrm{\Psi }}_0\left(x,0\right)\right|=\left|\mathrm{\Psi }\left(x,t\right)-{\mathcal{S}}_l\left(x\right)\right|\le C\left(1-x\right).\end{array}\right.\end{array}$

$$

Statement

Proof 1.

For the detailed proofs, readers can refer to [40].

When ε = 0, equation (4) assumes the following reduced form: 9 $\begin{array}{}\left\{\begin{array}{cc}{\mathrm{\Psi }}_t^0+p\left(x\right){\mathrm{\Psi }}_x^0+q\left(x\right){\mathrm{\Psi }}^0=f\left(x,t\right),& \forall \left(x,t\right)\in {\mathcal{D}}_{xt},\\ {\mathrm{\Psi }}^0\left(x,0\right)=S_0\left(x\right),& \forall x\in {\overline{\mathcal{D}}}_x,\\ {\mathrm{\Psi }}^0\left(0,t\right)=0,& \forall \in {\overline{\mathcal{D}}}_t.\end{array}\right.\end{array}$

Consider the first-order hyperbolic equation in (9) with initial data prescribed along t = 0 and x = 0 in ${\overline{\mathcal{D}}}_{xt}$. If C_ε is sufficiently small, then the solution Ψ(x, t) satisfies 10 $\begin{array}{}\mathrm{\Psi }\left(x,t\right)={\mathrm{\Psi }}^0\left(x,t\right)+O\left(C_{\epsilon }\right),\end{array}$ where Ψ⁰(x, t) is the solution of the reduced problem (9).

For subsequent analysis, let the differential operator $\mathcal{L}$ be defined as follows: 11 $\begin{array}{}{\mathcal{L}}_{C_{\epsilon }}\mathrm{\Psi }\left(x,t\right)\equiv {\mathrm{\Psi }}_t\left(x,t\right)-C_{\epsilon }{\mathrm{\Psi }}_{xx}\left(x,t\right)+p\left(x\right){\mathrm{\Psi }}_x\left(x,t\right)+q\left(x\right)\mathrm{\Psi }\left(x,t\right).\end{array}$

$$

Statement

Lemma 1 (continuous maximum principle).

Assume that Ψ(x, t) is any sufficiently smooth function satisfying $\mathrm{\Psi }\left(x,t\right)\ge 0,\forall \left(x,t\right)\in \mathrm{\partial }{\mathcal{D}}_{xt}$ and ${\mathcal{L}}_{C_{\epsilon }}\mathrm{\Psi }\left(x,t\right)\ge 0,\forall \left(x,t\right)\in {\mathcal{D}}_{xt}$; then, $\mathrm{\Psi }\left(x,t\right)\ge 0,\forall \left(x,t\right)\in {\overline{\mathcal{D}}}_{xt}$, where the differential operator is defined in (10).

$$

Statement

Proof 2.

Suppose that $\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)\in {\mathcal{D}}_{xt}$ such that $\mathrm{\Psi }\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)={\min }_{\left(x,t\right)\in {\overline{D}}_{xt}}\mathrm{\Psi }\left(x,t\right)<0$. These conditions and elementary concepts of calculus give Ψ_t(x^∗, t^∗) = 0, Ψ_x(x^∗, t^∗) = 0 and Ψ_xx(x^∗, t^∗) ≥ 0.

Using these conditions, we defined the differential operator as follows: 12 $\begin{array}{c}{\mathcal{L}}_{C_{\epsilon }}\psi \left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)={\psi }_t\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)-C_{\epsilon }{\psi }_{xx}\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)+p\left(x^{\mathrm{\ast }}\right){\psi }_x\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)+q\left(x^{\mathrm{\ast }}\right)\psi \left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)=-C_{\epsilon }{\psi }_{xx}\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)+q\left(x^{\mathrm{\ast }}\right)\psi \left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)\\ \le 000,since,{\psi }_{xx}\ge ,\text{and\hspace{0.17em}}\psi <\text{\hspace{0.17em}at}\left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right).\end{array}$

The result ${\mathcal{L}}_{C_{\epsilon }}\psi \left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)\le 0$ contradicts the assumption ${\mathcal{L}}_{C_{\epsilon }}\psi \left(x^{\mathrm{\ast }},t^{\mathrm{\ast }}\right)\ge 0$.

Therefore, it can be concluded that the minimum of ψ(x, t) is non-negative.

$$

Statement

Lemma 2 (stability estimate).

If $\mathrm{\Psi }\left(x,t\right)\in C^{21,}\left({\overline{\mathcal{D}}}_{xt}\right)$ is the solution of the continuous problem of equation (4), then it satisfies the bound 13 $\begin{array}{}\left|\mathrm{\Psi }\left(x,t\right)\right|\le \max \left\{\left|{\mathcal{S}}_0\left(x\right)\right|,\left|{\mathcal{S}}_l\left(t\right)\right|,\left|{\mathcal{S}}_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖.\end{array}$

$$

Statement

Proof 3.

To prove this Lemma, we defined two barrier functions Ω^± as follows: 14 $\begin{array}{}{\mathrm{\Omega }}^{\pm }\left(x,t\right)=\max \left\{\left|S_0\left(x\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm \mathrm{\Psi }\left(x,t\right).\end{array}$

When t = 0, we have 15 $\begin{array}{c}{\mathrm{\Omega }}^{\pm }\left(x,0\right)=\max \left\{\left|S_0\left(x\right)\right|,\left|S_l\left(0\right)\right|,\left|S_r\left(0\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm \mathrm{\Psi }\left(x,0\right)=\max \left\{\left|S_0\left(x\right)\right|,\left|S_l\left(0\right)\right|,\left|S_r\left(0\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm S_0\left(x\right)\ge 0.\end{array}$

When x = 0, we have 16 $\begin{array}{c}{\mathrm{\Omega }}^{\pm }\left(0,t\right)=\max \left\{\left|S_0\left(0\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm \mathrm{\Psi }\left(0,t\right)=\max \left\{\left|S_0\left(0\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm S_l\left(t\right)\ge 0.\end{array}$

When x = 1, we have 17 $\begin{array}{c}{\mathrm{\Omega }}^{\pm }\left(1,t\right)=\max \left\{\left|S_0\left(1\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm \mathrm{\Psi }\left(1,t\right)=\max \left\{\left|S_0\left(0\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm S_r\left(t\right)\ge 0.\end{array}$

Having all these on the domain ${\mathcal{D}}_{xt}$, we have 18 $\begin{array}{c}{\mathcal{L}}_{C_{\epsilon }}{\mathrm{\Omega }}^{\pm }\left(x,t\right)=\frac{\mathrm{\partial }{\mathrm{\Omega }}^{\pm }\left(x,t\right)}{\mathrm{\partial }t}-C_{\epsilon }\frac{{\mathrm{\partial }}^2{\mathrm{\Omega }}^{\pm }\left(x,t\right)}{\mathrm{\partial }{x}^2}+p\left(x\right)\frac{\mathrm{\partial }{\mathrm{\Omega }}^{\pm }\left(x,t\right)}{\mathrm{\partial }x}+q\left(x\right){\mathrm{\Omega }}^{\pm }\left(x,t\right)=q\left(x\right)\frac{1}{p^{\mathrm{\ast }}}‖f‖+q\left(x\right)\times \max \left\{\left|S_0\left(x\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm {\mathcal{L}}_{C_{\epsilon }}\mathrm{\Psi }\left(x,t\right)\\ =q\left(x\right)\times \max \left\{\left|S_0\left(x\right)\right|,\left|S_l\left(t\right)\right|,\left|S_r\left(t\right)\right|\right\}+\frac{1}{p^{\mathrm{\ast }}}‖f‖\pm f\left(x,t\right)\ge 0.\end{array}$

Using Lemma 1, it follows that ${\mathrm{\Omega }}^{\pm }\left(x,t\right)\ge 0,\forall \left(x,t\right)\in {\overline{\mathcal{D}}}_{xt}$, which implies the required bounds.

$$

Statement

Lemma 3.

The bounds for the derivative of the solution Ψ(x, t) of the problem in (4) are given by 19 $\begin{array}{}\left|\frac{{\mathrm{\partial }}^{i+m}\mathrm{\Psi }\left(x,t\right)}{\mathrm{\partial }{x}^i\mathrm{\partial }{t}^m}\right|\le C\left(1+{C_{\epsilon }}^{-i}\exp \left(\frac{p^{\mathrm{\ast }}\left(1-x\right)}{C_{\epsilon }}\right)\right),\end{array}$ where 0 ≤ m ≤ 2 and 0 ≤ i + 2m ≤ 4.

$$

Statement

Proof 4.

For the proof of this lemma, we refer the reader to [41].

3. Formulation of Numerical Schemes

The main concern of this section is to discretize the proposed problem. For the sake of simplicity, we first discretized in the temporal direction and then the spatial direction using uniform meshes as follows.

3.1. Temporal Semidiscretization

Dividing the time interval [0, T] into M equal subintervals with step size Δt = T/M, we obtain the discrete time levels as follows: 20 $\begin{array}{}{\overline{\mathcal{D}}}_t^M=\left\{t^m:t^m=m\mathrm{\Delta }t,m=0123,,,,\dots M\right\}.\end{array}$

Following the procedures in [34], the temporal semidiscretization of problem (4) is done by means of the Crank–Nicolson method as follows: 21 $\begin{array}{}\left\{\begin{array}{c}\frac{{\mathrm{\Psi }}^{m+1}-{\mathrm{\Psi }}^m}{\mathrm{\Delta }t}-\frac{C_{\epsilon }}{2}\left({\mathrm{\Psi }}_{xx}^{m+1}+{\mathrm{\Psi }}_{xx}^m\right)+\frac{p\left(x\right)}{2}\left({\mathrm{\Psi }}_x^{m+1}+{\mathrm{\Psi }}_x^m\right)+\frac{q\left(x\right)}{2}\left({\mathrm{\Psi }}^{m+1}+{\mathrm{\Psi }}^m\right)=\frac{1}{2}\left(f^{m+1}+f^m\right),\\ \left\{\begin{array}{cc}\mathrm{\Psi }\left(x,0\right)=S_0\left(x\right),& \forall x\in {\mathcal{D}}_x,\\ \mathrm{\Psi }\left(0,t^{m+1}\right)=S_l\left(t^{m+1}\right),& 01\le m\le M-,\\ \mathrm{\Psi }\left(1,t^{m+1}\right)=S_r\left(1,t^{m+1}\right),& 01\le m\le M-.\end{array}\right.\end{array}\right.\end{array}$

Simplifying and separating the (m + 1)^th and (m)^th time levels give 22 $\begin{array}{}\left\{\begin{array}{c}{\mathcal{L}}_{C_{\epsilon }}\equiv -C_{\epsilon }\left({\mathrm{\Psi }}_{xx}^{m+1}\right)\left(x\right)+A\left(x\right)\left({\mathrm{\Psi }}_x^{m+1}\right)\left(x\right)+B\left(x\right){\mathrm{\Psi }}^{m+1}\left(x\right)=H_i^m\left(x\right),\\ \left\{\begin{array}{cc}\mathrm{\Psi }\left(x,0\right)=S_0\left(x\right),& \forall x\in {\mathcal{D}}_x,\\ \mathrm{\Psi }\left(0,t^{m+1}\right)=S_l\left(t^{m+1}\right),& 01\le m\le M-,\\ \mathrm{\Psi }\left(1,t^{m+1}\right)=S_r\left(1,t^{m+1}\right),& 01\le m\le M-,\end{array}\right.\end{array}\right.\end{array}$ where 23 $\begin{array}{}\left\{\begin{array}{c}{H}_i^m\left(x\right)=C_{\epsilon }\left({\mathrm{\Psi }}_{xx}^m\right)\left(x\right)-A\left(x\right)\left({\mathrm{\Psi }}_x^m\right)\left(x\right)+C\left(x\right){\mathrm{\Psi }}^m\left(x\right)+\left(f^{m+1}+f^m\right)\left(x\right),\\ A\left(x\right)=a\left(x\right)-\delta b\left(x\right)+\eta d\left(x\right),\\ B\left(x\right)=\left(b\left(x\right)+c\left(x\right)+d\left(x\right)+\frac{2}{\mathrm{\Delta }t}\right),\\ C\left(x\right)=\left(2/\mathrm{\Delta }t-b\left(x\right)-c\left(x\right)-d\left(x\right)\right).\end{array}\right.\end{array}$

To check the semidiscrete maximum principle and stability estimate of the temporal semidiscretization of the proposed scheme, we first rewrite (22) at the (m + 1)^th time level in terms of spatial variables as follows: 24 $\begin{array}{}\left\{\begin{array}{cc}{{\mathcal{L}}_{{\mathcal{C}}_{\epsilon }}}^{\mathrm{\ast }}Y:\equiv -C_{\epsilon }{Y}^{″}\left(x\right)+A\left(x\right)Y^{\prime }\left(x\right)+B\left(x\right)Y=H\left(x\right),& \forall x\in {\overline{\mathcal{D}}}_x,\\ Y\left(0\right)=S_l\left(t^{m+1}\right),Y\left(1\right)=S_r\left(t^{m+1}\right),& \forall t^{m+1}\in {\overline{\mathcal{D}}}_t^M,\end{array}\right.\end{array}$ where ${{\mathcal{L}}_{{\mathcal{C}}_{\epsilon }}}^{\mathrm{\ast }}$ is a discrete operator, Y(x) = Ψ^m+1(x).

$$

Statement

Lemma 4 (semidiscrete maximum principle).

If ${\mathrm{\Psi }}^{m+1}\left(x\right)\in C^2\left({\overline{\mathcal{D}}}_x\right)$ is the solution of the semidiscrete problem (22) such that Ψ^m+1(0) ≥ 0, Ψ^m+1(1) ≥ 0 on $\mathrm{\partial }{\overline{\mathcal{D}}}_x$ and ${\mathcal{L}}_{{\mathcal{C}}_{\epsilon }}{\mathrm{\Psi }}^{m+1}\left(x\right)\ge 0,\forall x\in {\mathcal{D}}_x$, then ${\mathrm{\Psi }}^{m+1}\left(x\right)\ge 0,\forall x\in {\overline{\mathcal{D}}}_x$.

$$

Statement

Proof 5.

Let $x^{\mathrm{\ast }}\in {\overline{\mathcal{D}}}_x$ such that ${\mathrm{\Psi }}^{m+1}\left(x^{\mathrm{\ast }}\right)={\min }_{x\in {\overline{\mathcal{D}}}_x}{\mathrm{\Psi }}^{m+1}\left(x\right)<0$, and since Ψⁿ⁺¹(0) ≥ 0, Ψ^m+1(1) ≥ 0, then x^∗ ∉ {0, 1}, that is, $x^{\mathrm{\ast }}\in {\mathcal{D}}_x$.

Moreover, since ${\psi }_x^{m+1}\left(x^{\mathrm{\ast }}\right)=0$ and ${\psi }_{xx}^{m+1}\left(x^{\mathrm{\ast }}\right)\ge 0$, we have that 25 $\begin{array}{}{\mathcal{L}}_{C_{\epsilon }}{\psi }^{m+1}\left(x^{\mathrm{\ast }}\right)=-C_{\epsilon }{\psi }_{xx}^{m+1}\left(x^{\mathrm{\ast }}\right)+A\left(x^{\mathrm{\ast }}\right){\psi }_x^{m+1}\left(x^{\mathrm{\ast }}\right)+B\left(x^{\mathrm{\ast }}\right){\psi }^{m+1}\left(x^{\mathrm{\ast }}\right)\le 0.\end{array}$

This implies ${\mathcal{L}}_{C_{\epsilon }}{\mathrm{\Psi }}^{m+1}\left(x\right)\le 0$, which contradicts our assumption ${\mathcal{L}}_{C_{\epsilon }}{\mathrm{\Psi }}^{m+1}\left(x\right)\ge 0$.

Hence, the minimum of Ψ^m+1(x) is non-negative.

The local truncation error at the (m + 1)^th time level is given by e^m+1 ≡ Ψ(x, t_m+1) − Ψ^m+1, where Ψ^m+1 is the computed solution of the semidiscrete problem (22) and Ψ(x, t_m+1) is the exact solution of the continuous problem (4). Moreover, this error measures the contribution of each time step to the global error of the temporal semidiscretization.

$$

Statement

Lemma 5 (local error estimate).

Suppose that the bound 26 $\begin{array}{}\left|\frac{{\mathrm{\partial }}^m\mathrm{\Psi }\left(x,t\right)}{\mathrm{\partial }{t}^m}\right|\le C,\forall \left(x,t\right)\in {\overline{\mathcal{D}}}_{xt}\text{\hspace{0.17em}and\hspace{0.17em}}03\le m\le ,\end{array}$ then the local error estimate in the temporal discretization is given by 27 $\begin{array}{}{‖e^{m+1}‖}_{\infty }\le C_1{\left(\mathrm{\Delta }t\right)}^3,m=1231,,,\dots ,M-.\end{array}$

$$

Statement

Proof 6.

Using the Taylor series expansion for Ψ(x, t^m) and Ψ(x, t^m+1) about (x, t^m+1/2), we get 28 $\begin{array}{}\mathrm{\Psi }\left(x,t^{m+1}\right)=\mathrm{\Psi }\left(x,t^{m+12/}\right)+\frac{\mathrm{\Delta }t}{2}{\mathrm{\Psi }}_t\left(x,t^{m+12/}\right)+\frac{{\left(\mathrm{\Delta }t\right)}^2}{8}{\mathrm{\Psi }}_{tt}\left(x,t^{m+12/}\right)+O\left({\left(\mathrm{\Delta }t\right)}^3\right),\end{array}$ 29 $\begin{array}{}\mathrm{\Psi }\left(x,t^m\right)=\mathrm{\Psi }\left(x,t^{m+12/}\right)-\frac{\mathrm{\Delta }t}{2}{\mathrm{\Psi }}_t\left(x,t^{m+12/}\right)+\frac{{\left(\mathrm{\Delta }t\right)}^2}{8}{\mathrm{\Psi }}_{tt}\left(x,t^{m+12/}\right)+O\left({\left(\mathrm{\Delta }t\right)}^3\right).\end{array}$

Now, by subtracting (29) from (28), we get 30 $\begin{array}{c}\frac{\mathrm{\Psi }\left(x,t^{m+1}\right)-\mathrm{\Psi }\left(x,t^m\right)}{\mathrm{\Delta }t}={\mathrm{\Psi }}_t\left(x,t^{m+12/}\right)+O\left({\left(\mathrm{\Delta }t\right)}^2\right)=C_{\epsilon }{\mathrm{\Psi }}_{xx}\left(x,t^{m+12/}\right)-A\left(x\right){\mathrm{\Psi }}_x\left(x,t^{m+12/}\right)\\ -B\left(x\right)\mathrm{\Psi }\left(x,t^{m+12/}\right)+H\left(x,t^{m+12/}\right)+O\left({\left(\mathrm{\Delta }t\right)}^2\right).\end{array}$

The solution to the local error is 31 $\begin{array}{}\left\{\begin{array}{c}{\mathcal{L}}^{\mathrm{\ast }}{e}^{m+1}=O\left({\left(\mathrm{\Delta }t\right)}^3\right),\\ e^{m+1}\left(0\right)=00,e^{m+1}\left(1\right)=.\end{array}\right.\end{array}$

Thus, using the maximum principle, the required result is satisfied.

Hence, ${‖e^{m+1}‖}_{\infty }\le C_1{\left(\mathrm{\Delta }t\right)}^3$ as required.

$$

Statement

Lemma 6 (global error estimate).

Under the hypothesis of Lemma 5, the global error estimate E^m+1 in the temporal direction is defined by 32 $\begin{array}{}{‖E^{m+1}‖}_{\infty }\le C_2{\left(\mathrm{\Delta }t\right)}^2.\end{array}$

$$

Statement

Proof 7.

Using the local error estimates up to the (m + 1)^th time step given by Lemma 3, we get the following global error estimate: 33 $\begin{array}{c}{‖E^{m+1}‖}_{\infty }={‖{\displaystyle \sum }_{s=1}^{m+1}{e}^s‖}_{\infty }\le {‖e^1‖}_{\infty }+{‖e^2‖}_{\infty }+\cdots +{‖e^s‖}_{\infty }\le C_2\left(m+1\right){\left(\mathrm{\Delta }t\right)}^3\\ =C_2\left(T\right){\left(\mathrm{\Delta }t\right)}^2,\text{since}\left(m+1\right)\mathrm{\Delta }t\le T\\ \le C_1{\left(\mathrm{\Delta }t\right)}^2,C_1=C_{2}T,\end{array}$ where C₁ is a constant independent of the parameter C_ε and the step length Δt.

Hence, the Crank–Nicolson method is a second-order parameter uniform convergent method.

$$

Statement

Lemma 7.

The derivative of the solution of problem (3) satisfies the bound 34 $\begin{array}{}\left|\frac{d^i\mathrm{\Psi }\left(x,t^{m+1}\right)}{d{x}^i}\right|\le C\left(1+{C_{\epsilon }}^{-i}\exp \left(\frac{p^{\mathrm{\ast }}\left(1-x\right)}{C_{\epsilon }}\right)\right),x\in {\overline{\mathcal{D}}}_x,04\le i\le .\end{array}$

$$

Statement

Proof 8.

Interested readers are referred to [42] for the proof of this lemma.

3.2. Spatial Discretization

For accurate results, an effective numerical approach is essential, especially when dealing with singular perturbation problems where small parameters cause abrupt changes in the solution. We present a finite difference scheme that is exponentially fitted and uses an adapted operator to preserve the diffusion term.

3.2.1. Exponentially Fitting Factor

In the present work, we use the exponentially fitted operator finite difference method which helps us to hinder the influence of the singular perturbation parameter (C_ε) as C_ε⟶0 in the boundary layer region.

For each m = 0, 1, 2, 3, …M − 1, we have the boundary value problem of the form 35 $\begin{array}{}{\mathcal{L}}^{\mathrm{\Delta }t}{\mathrm{\Psi }}^{m+1}\left(x\right)\equiv -C_{\epsilon }\frac{d^2}{d{x}^2}{\mathrm{\Psi }}^{m+1}\left(x\right)+A\left(x\right)\frac{d}{dx}{\mathrm{\Psi }}^{m+1}\left(x\right)+B\left(x\right){\mathrm{\Psi }}^{m+1}\left(x\right)=H^{m+1}\left(x\right),\end{array}$ with its boundary conditions 36 $\begin{array}{}{\mathrm{\Psi }}^{m+1}\left(0\right)={\varphi }_l\left(0,t^{m+1}\right),{\mathrm{\Psi }}^{m+1}\left(1,t^{m+1}\right)={\varphi }_r\left(1,t^{m+1}\right),01\le m\le M-.\end{array}$

To solve the singularly perturbed boundary value problem in (35)–(36), we use the zeroth-order asymptotic solution for the right boundary layer presented in [43]. 37 $\begin{array}{}{\mathrm{\Psi }}^{m+1}\left(x\right)={\mathrm{\Psi }}_0^{m+1}\left(0\right)+\left({\varphi }_r\left(1,t^{m+1}\right)-{\mathrm{\Psi }}_0^{m+1}\left(1\right)\right)\exp \left(-\frac{A\left(1\right)\left(1-x\right)}{C_{\epsilon }}\right)+O\left(C_{\epsilon }\right),\end{array}$ where ${\mathrm{\Psi }}_0^{m+1}$ is the solution of the reduced problem (9).

On discretizing the spatial domain ${\mathcal{D}}_x=\left[01,\right]$ into N equal number of subintervals each of length h = 1/N and x₀ = 0, x_i = ih, x_N = 1, for i = 0, 1, 2, …, N is the mesh point.

Taking h very small, the discretized form of (37) becomes 38 $\begin{array}{}\left\{\begin{array}{c}{\mathrm{\Psi }}^{m+1}\left(ih\right)={\mathrm{\Psi }}_0^{m+1}\left(ih\right)+\left({\varphi }_r\left(1,t^{m+1}\right)-{\mathrm{\Psi }}_0^{m+1}\left(1\right)\right)\exp \left(-A\left(1\right)\left(1/C_{\epsilon }-i\rho \right)\right),\\ {\mathrm{\Psi }}_{i-1}^{m+1}\left(ih\right)={\mathrm{\Psi }}_0^{m+1}\left(ih\right)+\left({\varphi }_r\left(1,t^{m+1}\right)-{\mathrm{\Psi }}_0^{m+1}\left(1\right)\right)\exp \left(-A\left(1\right)\left(1/C_{\epsilon }-\left(i-1\right)\rho \right)\right),\\ {\mathrm{\Psi }}_{i+1}^{m+1}\left(ih\right)={\mathrm{\Psi }}_0^{m+1}\left(ih\right)+\left({\varphi }_r\left(1,t^{m+1}\right)-{\mathrm{\Psi }}_0^{m+1}\left(1\right)\right)\exp \left(-A\left(1\right)\left(1/C_{\epsilon }-\left(i+1\right)\rho \right)\right),\end{array}\right.\end{array}$ where ρ = h/C_ε.

To handle the oscillations or abrupt changes to the solution of the proposed problem due to the presence of the perturbation parameter, we introduce an exponential fitting factor σ(ρ) on the diffusion term and write as follows: 39 $\begin{array}{}-\sigma \left(\rho \right)C_{\epsilon }{\delta }_x^2{\mathrm{\Psi }}^{m+1}\left(x\right)+A\left(x\right){\delta }_x^0{\mathrm{\Psi }}^{m+1}\left(x\right)+B\left(x\right){\mathrm{\Psi }}^{m+1}\left(x\right)=H^{m+1}\left(x\right).\end{array}$

We use the uniform mesh points ${\left\{x_i\right\}}_{i=0}^N$ with mesh length h = x_i+1 − x_i and any mesh function Ψ_i to denote the following finite difference operators for multiple purposes: 40 $\begin{array}{c}{\delta }^{+}{\mathrm{\Psi }}_i^m=\frac{{\mathrm{\Psi }}_{i+1}^m-{\mathrm{\Psi }}_i^m}{h},{\delta }_x^{-}{\mathrm{\Psi }}_i^m=\frac{{\mathrm{\Psi }}_i^m-{\mathrm{\Psi }}_{i-1}^m}{h},\\ {\delta }_x^0{\mathrm{\Psi }}_i^m=\frac{{\mathrm{\Psi }}_{i+1}^m-{\mathrm{\Psi }}_{i-1}^m}{2h}\\ {\delta }_x^2{\mathrm{\Psi }}_i^m=\frac{{\mathrm{\Psi }}_{i+1}-2{\mathrm{\Psi }}_i+{\mathrm{\Psi }}_{i-1}}{h^2}.\end{array}$

Using the finite difference operators in (40) in (39) gives 41 $\begin{array}{}-\sigma \left(\rho \right)C_{\epsilon }\frac{{\mathrm{\Psi }}_{i+1}^{m+1}-2{\mathrm{\Psi }}_i^{m+1}+{\mathrm{\Psi }}_{i-1}^{m+1}}{h^2}+A\left(x_i\right)\frac{{\mathrm{\Psi }}_{i+1}^m-{\mathrm{\Psi }}_{i-1}^{m+1}}{2h}+B\left(x_i\right){\mathrm{\Psi }}_i^{m+1}=H\left(x_i,t^{m+1}\right).\end{array}$

Multiplying (41) by h, where h⟶0, and truncating the term $h\left(H\left(x_i,t^{m+1}\right)-B\left(x_i\right){\mathrm{\Psi }}_i^{m+1}\right)$ give 42 $\begin{array}{}-\frac{\sigma \left(\rho \right)}{\rho }\left({\mathrm{\Psi }}_{i+1}^{m+1}-2{\mathrm{\Psi }}_i^{m+1}+{\mathrm{\Psi }}_{i-1}^{m+1}\right)+\frac{A\left(x_i\right)}{2}\left({\mathrm{\Psi }}_{i+1}^m-{\mathrm{\Psi }}_{i-1}^{m+1}\right)=0.\end{array}$

Using (38) into (43) and simplifying it, we get the required fitting factor as follows: 43 $\begin{array}{}\sigma \left(\rho \right)=\frac{\rho A\left(x_i\right)}{2}\coth \left(\frac{\rho A\left(x_i\right)}{2}\right).\end{array}$

3.2.2. The Fully Discrete Scheme

Using (39) and (43) in (41), we get the following tridiagonal system of equations: 44 $\begin{array}{}\left\{\begin{array}{c}{r}_i^0{\mathrm{\Psi }}_{i-1}+r_i^{-}{\mathrm{\Psi }}_i+r_i^{+}{\mathrm{\Psi }}_{i+1}=H_i^{m+1},\\ \left\{\begin{array}{c}{r}_i^0=\frac{\sigma \left(\rho \right)C_{\epsilon }}{h{.}^2}-\frac{A_i^{m+1}}{2h},\\ r_i^{-}=2\frac{\sigma \left(\rho \right)C_{\epsilon }}{h^2}+B_i^{m+1},\\ r_i^{+}=\frac{\sigma \left(\rho \right)C_{\epsilon }}{h^2}+\frac{A_i^{m+1}}{2h},\\ H_i^{m+1}=C_{\epsilon }\sigma \left({\mathrm{\Psi }}_{i-1}^m-2{\mathrm{\Psi }}_i^m+{\mathrm{\Psi }}_{i+1}^m\right)-\frac{A_i^m}{2}\left({\mathrm{\Psi }}_{i+1}^m-{\mathrm{\Psi }}_{i-1}^m\right)+C_i^m{\mathrm{\Psi }}_i^m+f_i^{m+1}+f_i^m.\end{array}\right.\end{array}\right.\end{array}$ where A_i = a_i − δb_i + ηd_i, B_i = b_i + c_i + d_i + 2/Δt and C_i = −b_i − c_i − d_i + 2/Δt.

4. Stability and Uniform Convergence Analysis

$$

Statement

Lemma 8 (discrete comparison principle).

There exists a comparison function $W_i^{m+1}$ such that ${\mathcal{L}}^{h,\mathrm{\Delta }t}{\mathrm{\Psi }}_i^{m+1}\le {\mathcal{L}}^{h,\mathrm{\Delta }t}{W}_i^{m+1},i=1231,,,\dots ,N-$ and if ${\mathrm{\Psi }}_0^{m+1}\le W_0^{m+1}$ and ${\mathrm{\Psi }}_N^{m+1}\le W_N^{m+1}$, then ${\mathrm{\Psi }}_i^{m+1}\le W_i^{m+1},\forall i=012,,,\dots ,N$.

$$

Statement

Proof 9.

Since the matrix associated with operator ${\mathcal{L}}^{h,\mathrm{\Delta }t}$ in (44) is of size [N + 1]², it satisfies the property of the M matrix as proved in [44]. Thus, the result of this lemma guarantees the existence and uniqueness of the discrete solution.

$$

Statement

Lemma 9.

Let $W_i^{m+1}=10+x_i,i\le \le N$, then there exists a positive constant C such that $W_i^{m+1}$ and ${\mathcal{L}}^{h,\mathrm{\Delta }t}{W}_i^{m+1}\ge C,11\le i\le N-$.

$$

Statement

Proof 10.

The proof is straightforward and provides a consistent bound in C_ε for the norm of the inverse ${\mathcal{L}}^{h,\mathrm{\Delta }t}$. 45 $\begin{array}{c}{\mathcal{L}}^{h,\mathrm{\Delta }t}{W}_i^{m+1}=-\sigma \left(\rho \right)C_{\epsilon }{\delta }_x^2{W}^{m+1}\left(x\right)+A\left(x\right){\delta }_x^0{W}^{m+1}\left(x\right)+B\left(x\right)W^{m+1}\left(x\right)=-\sigma \left(\rho \right)C_{\epsilon }{\delta }_x^2\left(1+x_i\right)+A\left(x\right){\delta }_x^0\left(1+x_i\right)+B\left(x\right)\left(1+x_i\right)\ge C.\end{array}$

$$

Statement

Lemma 10.

The solution ${\mathrm{\Psi }}_i^{m+1}$ of the discrete scheme in (44) satisfies the bound 46 $\begin{array}{}\left|{\mathrm{\Psi }}_i^{m+1}\right|\le \frac{‖H‖}{\beta }+C\max \left\{\left|S_l\left(0,t_i^{m+1}\right)\right|,\left|S_r\left(1,t_i^{m+1}\right)\right|\right\},B\left(x\right)\ge \beta >0.\end{array}$

$$

Statement

Proof 11.

Considering barrier functions ${\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}=\left(\left(‖H‖\right)/\beta \right)+C\max \left\{\left|S_l\left(0,t_i^{m+1}\right)\right|,\left|S_r\left(1,t_i^{m+1}\right)\right|\right\}\pm {\mathrm{\Psi }}_i^{m+1}$, it can be shown that ${\left({\mathrm{\Omega }}^{\pm }\right)}_0^{m+1}\ge 0$ and ${\left({\mathrm{\Omega }}^{\pm }\right)}_N^{m+1}\ge 0$. Then, for 1 ≤ i ≤ N − 1, we have 47 $\begin{array}{c}{\mathcal{L}}^{h,\mathrm{\Delta }t}{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}=-\sigma \left(\rho \right)C_{\epsilon }{\delta }_x^2{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}+A\left(x\right){\delta }_x^0{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}+B\left(x\right){\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1},=B\left(x\right)\left[\frac{‖H‖}{\beta }+C\max \left\{\left|S_l\left(0,t_i^{m+1}\right)\right|,\left|S_r\left(1,t_i^{m+1}\right)\right|\right\}\pm {\mathrm{\Psi }}_i^{m+1}\right]\ge 0.\end{array}$

Using Lemma 8, we obtain ${\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}\ge 0012,\text{for\hspace{0.17em}}i=,,,\dots ,N$.

$$

Statement

Lemma 11 (uniform stability estimate).

At any (m + 1)^th time level, if ${\varphi }_i^{m+1}$ is any mesh function such that ${\varphi }_0^{m+1}=0={\varphi }_N^{m+1}$, then 48 $\begin{array}{}\left|{\varphi }_i^{m+1}\right|\le \frac{1}{\beta }\underset{11\le j\le N-}{\max }\left|{\mathcal{L}}^{h,\mathrm{\Delta }t}{\varphi }_i^{m+1}\right|,\text{\hspace{1em}for\hspace{0.17em}}110\le i\le N-,\text{where\hspace{0.17em}}d\left(x\right)\ge \beta >.\end{array}$

$$

Statement

Proof 12.

Consider two barrier functions of the form 49 $\begin{array}{}{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}=\frac{1}{\beta }\underset{11\le i\le N-}{\max }\left|{\mathcal{L}}^{h,\mathrm{\Delta }t}{\varphi }_i^{m+1}\right|\pm {\varphi }_i^{m+1}.\end{array}$

Using the result of Lemma 10, it is clear that ${\left({\mathrm{\Omega }}^{\pm }\right)}_0^{m+1}\ge 00,{\left({\mathrm{\Omega }}^{\pm }\right)}_N^{m+1}\ge$. Having these we want to show that ${\mathcal{L}}^{h,\mathrm{\Delta }t}{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}\ge 0$ for 1 ≤ i ≤ N − 1. 50 $\begin{array}{c}{\mathcal{L}}^{h,\mathrm{\Delta }t}{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}=-\sigma \left(\rho \right)C_{\epsilon }{\delta }_x^2{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}+A\left(x\right){\delta }_x^0{\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}+B\left(x\right){\left({\mathrm{\Omega }}^{\pm }\right)}_i^{m+1}=B\left(x_i\right)\left[\frac{1}{\beta }\underset{11\le i\le N-}{\max }\left|{\mathcal{L}}^{h,\mathrm{\Delta }t}{\varphi }_i^{m+1}\right|\pm {\varphi }_i^{m+1}\right]\ge 0.\end{array}$

Consequently, using the discrete comparison principle gives the bound 51 $\begin{array}{}\left|{\varphi }_i^{m+1}\right|\le \frac{1}{\beta }\underset{11\le i\le N-}{\max }\left|{\mathcal{L}}^{h,\mathrm{\Delta }t}{\varphi }_i^{m+1}\right|.\end{array}$

Overall, this result guarantees the uniform stability of the discrete scheme.

$$

Statement

Lemma 12.

Let the coefficient and source functions A(x), B(x) and H(x, t^m+1) of the semidiscrete problem in (22) be sufficiently smooth so that Ψ^m+1(x) ∈ C⁴[0, 1]. Then, the discrete solution ${\mathrm{\Psi }}_i^{m+1}$ of (44) satisfies the following error bound: 52 $\begin{array}{}‖{\mathcal{L}}^{h,\mathrm{\Delta }t}\left({\mathrm{\Psi }}^{m+1}\left(x\right)-{\mathrm{\Psi }}_i^{m+1}\right)‖\le \frac{C{h}^2}{h+C_{\epsilon }}\left(1+{C_{\epsilon }}^{-3}\exp \left(-\frac{p^{\mathrm{\ast }}\left(1-x_i\right)}{C_{\epsilon }}\right)\right).\end{array}$

$$

Statement

Proof 13.

For a detailed proof of this lemma, we refer the reader to [37].

$$

Statement

Lemma 13.

For a fixed mesh and for all integers k, we have 53 $\begin{array}{c}\underset{C_{\epsilon }⟶0}{\lim }\underset{11\le i\le N-}{\max }\frac{\exp \left(\left(-p^{\mathrm{\ast }}{x}_i\right)/\left(C_{\epsilon }\right)\right)}{{C_{\epsilon }}^k}=0,\underset{C_{\epsilon }⟶0}{\lim }\underset{11\le i\le N-}{\max }\frac{\exp \left(\left(-p^{\mathrm{\ast }}\left(1-x_i\right)\right)/\left(C_{\epsilon }\right)\right)}{{C_{\epsilon }}^k}=0,\end{array}$ where and x_i = ih, h = 1/N, 1 ≤ i ≤ N − 1.

$$

Statement

Proof 14.

For a detailed proof of this lemma, we refer the reader to [26, 37].

Using Lemma 13 in Lemma 12 gives 54 $\begin{array}{}‖{\mathcal{L}}^{h,\mathrm{\Delta }t}\left({\mathrm{\Psi }}^{m+1}\left(x\right)-{\mathrm{\Psi }}_i^{m+1}\right)‖\le \frac{C{h}^2}{h+C_{\epsilon }}.\end{array}$

So, applying the bound in Lemma 11, we get 55 $\begin{array}{}‖{\mathrm{\Psi }}^{m+1}\left(x\right)-{\mathrm{\Psi }}_i^{m+1}‖\le \frac{C{h}^2}{h+C_{\epsilon }}.\end{array}$

Moreover, for a small value of C_ε ≪ h, the denominator of the bound in (55) is dominated by h. From this, we conclude that the proposed scheme is of first order in the spatial direction and defined by 56 $\begin{array}{}\underset{01<C_{\epsilon }\le }{\sup }‖{\mathrm{\Psi }}^{m+1}\left(x\right)-{\mathrm{\Psi }}_i^{m+1}‖\le Ch.\end{array}$

Based on the results established in Lemma 2 and (55), we can generalize the following theorem.

$$

Statement

Theorem 1.

Suppose that Ψ^m+1(x_i) and ${\mathrm{\Psi }}_i^{m+1}$ are the exact and computed solutions of (4) and (44), respectively, on the discretized domain.Then, for sufficiently small h, the following uniform error estimate holds: 57 $\begin{array}{}\underset{01<C_{\epsilon }\le }{\sup }‖{\mathrm{\Psi }}^{m+1}\left(x_i\right)-{\mathrm{\Psi }}_i^{m+1}‖\le C\left({h+\left(\mathrm{\Delta }t\right)}^2\right).\end{array}$

From this, we conclude that the proposed scheme is first-order accurate in the spatial direction and second-order accurate in the temporal direction.

5. Richardson’s Extrapolation

As stated in Theorem 1, it is seen that the proposed scheme has a low order of accuracy in the spatial direction. This affects the overall order of convergence of the scheme. To tackle this challenge, we follow a similar approach as described in [45] and use the Richardson extrapolation technique in the spatial direction.

Let ${\mathcal{D}}_x^{2N}$ be the mesh obtained from bisecting the step length h such that ${\mathcal{D}}_x^N\subseteq {\mathcal{D}}_x^{2N}$. If ${\mathrm{\Psi }}_i^m\text{\hspace{0.17em}and\hspace{0.17em}}{\overline{\mathrm{\Psi }}}_i^m$ are approximate solutions obtained from ${\mathcal{D}}_x^N\text{\hspace{0.17em}and\hspace{0.17em}}{\mathcal{D}}_x^{2N}$, respectively, then 58 $\begin{array}{}{\psi }^m\left(t\right)-{\mathrm{\Psi }}_i^m\left(t\right)=Ch+R_i^N\left(t\right),\left(x,t^m\right)\in {\mathcal{D}}_x^N\times {\mathcal{D}}_t,\end{array}$ 59 $\begin{array}{}{\psi }^m\left(t\right)-{\overline{\mathrm{\Psi }}}_i^m\left(t\right)=C\left(\frac{h}{2}\right)+R_i^{2N}\left(t\right),\left(x,t^m\right)\in {\mathcal{D}}_x^{2N}\times {\mathcal{D}}_t,\end{array}$ where ψ^m(t) is the exact solution of the proposed problem at the m^th time level, $R_m^N\left(t\right)$ and $R_m^{2N}\left(t\right)$ are the remainder terms of the error with order O(h²), and C is a constant independent of the parameters ε, h and Δt. Now, using the inequalities (58) and (59), we get the error bound of the following form: 60 $\begin{array}{}\underset{0\le i\le N}{\max }\left|{\mathrm{\Psi }}^{\text{ext}}\left(t\right)-{\mathrm{\Psi }}^m\left(t\right)\right|\le C{h}^2,\end{array}$ where Ψ^ext(t) is the approximate solution obtained from the extrapolation [15] and C is a constant independent of the parameter ε and step size h and Δt.

From these constructions, Theorem 1 can be modified in the following form.

$$

Statement

Theorem 2.

Suppose that Ψ^m+1(x_i) and ${\mathrm{\Psi }}_i^{m+1}$ are the exact and computed solutions of (4) and (44), respectively, on the discretized domain.Then, for sufficiently small h, the following uniform error estimate holds: 61 $\begin{array}{}\underset{01<C_{\epsilon }\le }{\sup }‖{\mathrm{\Psi }}^{m+1}\left(x_i\right)-{\mathrm{\Psi }}_i^{m+1}‖\le C\left({h^2+\left(\mathrm{\Delta }t\right)}^2\right).\end{array}$

This concludes that the proposed scheme is second-order accurate in both the spatial direction and temporal directions after the application of Richardson extrapolation.

6. Numerical Results and Discussion

This section deals with the validation of the proposed scheme by taking three model problems. These model problems lack analytical solutions due to different reasons such as high dimensionality or complex parameters. To address these challenges, we estimate absolute maximum errors using the double mesh principle as outlined in [46–48].

$$

Statement

Example 1.

In this example, we took the coefficient and source terms of the considered problem in (1) of the form a(x) = 2 − x², b(x) = 2, c(x) = x − 3, d(x) = 1 and f(x, t) = 10t²e^−tx(1 − x) for T = 3, Ψ₀(x) = 0, S_l(x) = S_r(x) = 0, δ = 0.6ε and η = 0.5ε [32].

$$

Statement

Example 2.

In this example, we took the coefficient and source terms of the considered problem in (1) of the form a(x) = 2 − x², b(x) = 1 + x, c(x) = 1 + x² + cos(πx), d(x) = 3 and f(x) = sin(πx) for T = 3, Ψ₀(x) = 0, S_l(x) = S_r(x) = 0, δ = 0.6ε and η = 0.5ε [35].

$$

Statement

Example 3.

In this example, we took the coefficient and source terms of the considered problem in (1) of the form a(x) = 1, b(x) = 2 − x², c(x) = 1 + x², d(x) = e^x and f(x) = 50(x(1 − x))³ for $T=200.6,S_l\left(x,t\right)=S_r\left(x,t\right)=S_0\left(x\right)=,\delta =\sqrt{\epsilon }$ and $\eta =0.5\sqrt{\epsilon }$ [36].

Let Ψ^N,M(x_i, t^m) and Ψ^2N,2M(x_i, t^m) be the approximate solutions with mesh lengths (h, Δt ≠ 0) and (h/2, Δt/2 ≠ 0), respectively, then

• •

  The absolute maximum error is given by 62 $\begin{array}{}{E}_{\epsilon }^{N,M}=\underset{0\le i,m\le N,M}{\max }\left|{\mathrm{\Psi }}^{N,M}\left(x_i,t^m\right)-{\mathrm{\Psi }}^{22N,M}\left(x_i,t^m\right)\right|.\end{array}$

• •

  The ε− uniform absolute maximum error is given by 63 $\begin{array}{}{E}^{N,M}=\underset{\epsilon }{\max }{E}_{\epsilon }^{N,M}.\end{array}$

• •

  The corresponding ε-uniform order of convergence is given by 64 $\begin{array}{}{R}^{N,M}={\log }_2\left(\frac{E^{N,M}}{E^{22N,M}}\right).\end{array}$

The tabulated results presented in Tables 1, 2, and 3, corresponding to Examples 1, 2, and 3, respectively, display the absolute maximum errors and the associated orders of convergence. For a fixed value of the perturbation parameter ε, an increase in the number of mesh points N and M results in a consistent decrease in the absolute maximum errors. Conversely, for fixed values of N and M, and the shift arguments δ and η, a reduction in ε gives nearly constant absolute maximum errors, thereby confirming that the proposed method exhibits parameter-uniform convergence. Furthermore, the two- and three-dimensional plots illustrated in Example 2 (Figures 1 and 2) demonstrate that, for a fixed mesh size and decreasing ε, the numerical solution develops a steep gradient near x = 1, confirming the existence of a right boundary layer as anticipated in the theoretical framework. The log–log plots shown in Figures 1 and 2 provide additional evidence of the uniform convergence properties of the proposed scheme, which is in complete agreement with the theoretical predictions. Upon applying the extrapolation technique, the proposed scheme achieves second-order accuracy with respect to both space and time. Furthermore, the comparison results presented in Tables 4 and 5 for Examples 1 and 2 demonstrate that the proposed method consistently outperforms several existing numerical approaches reported in the literature [32, 33, 35], and [36] (Figures 3 and 4).

Table 1 E^N,M and R^N,M for Example 1 with δ = 0.6∗ε, η = 0.5∗ε.

ε↓	N = 32	N = 64	N = 128	N = 256	N = 512
M = 60	M = 120	M = 240	M = 480	M = 960
Before Ext.
2−6	7.8834e − 03	3.1388e − 03	9.5444e − 04	2.5306e − 04	6.4253e − 05
2−8	8.1403e − 03	4.1053e − 03	1.9873e − 03	7.8599e − 04	2.3855e − 04
2−10	8.1303e − 03	4.1031e − 03	2.0591e − 03	1.0305e − 03	4.9740e − 04
2−12	8.1278e − 03	4.1018e − 03	2.0585e − 03	1.0309e − 03	5.1582e − 04
2−14	8.1271e − 03	4.1015e − 03	2.0583e − 03	1.0308e − 03	5.1578e − 04
2−16	8.1270e − 03	4.1014e − 03	2.0583e − 03	1.0308e − 03	5.1577e − 04
EN,M	8.1403e − 03	4.1053e − 03	1.9873e − 03	1.0309e − 03	5.1582e − 04
RN,M	0.9876	1.0467	0.9469	0.9990	—
After Ext.
2−6	5.1214e − 04	1.4288e − 04	5.1182e − 05	1.5237e − 05	4.2499e − 06
2−8	5.3214e − 04	1.6888e − 04	5.2882e − 05	1.6037e − 05	4.6499e − 06
2−10	5.3214e − 04	1.6888e − 04	5.2882e − 05	1.6037e − 05	4.6499e − 06
2−12	5.3214e − 04	1.6888e − 04	5.2882e − 05	1.6037e − 05	4.6499e − 06
2−14	5.3214e − 04	1.6888e − 04	5.2882e − 05	1.6037e − 05	4.6499e − 06
2−16	5.3214e − 04	1.6888e − 04	5.2882e − 05	1.6037e − 05	4.6499e − 06
EN,M	5.3214e − 04	1.6888e − 04	5.2882e − 05	1.6037e − 05	4.6499e − 06
RN,M	1.6557	1.6752	1.7862	1.7213	—

Table 2 E^N,M and R^N,M for Example 2 with δ = 0.6∗ε, η = 0.5∗ε.

ε↓	N = 32	N = 64	N = 128	N = 256	N = 512
M = 60	M = 120	M = 240	M = 480	M = 960
Before Ext.
2−6	1.0036e − 02	5.3432e − 03	2.2305e − 03	6.9293e − 04	1.8547e − 04
2−8	1.0075e − 02	5.5408e − 03	2.8976e − 03	1.4241e − 03	5.7175e − 04
2−10	1.0083e − 02	5.5458e − 03	2.9021e − 03	1.4769e − 03	7.4790e − 04
2−12	1.0083e − 02	5.5458e − 03	2.9021e − 03	1.4769e − 03	7.4790e − 04
2−14	1.0083e − 02	5.5458e − 03	2.9021e − 03	1.4769e − 03	7.4790e − 04
2−16	1.0083e − 02	5.5458e − 03	2.9021e − 03	1.4769e − 03	7.4790e − 04
EN,M	1.0083e − 02	5.5458e − 03	2.9021e − 03	1.4769e − 03	7.4790e − 04
RN,M	0.8627	0.9342	0.9745	0.9808	—
After Ext.
2−6	2.0102e − 04	6.2004e − 05	1.6925e − 05	5.0489e − 06	1.1123e − 06
2−8	2.1042e − 04	6.2022e − 05	1.7905e − 05	5.0708e − 06	1.3832e − 06
2−10	2.1042e − 04	6.2022e − 05	1.7905e − 05	5.0708e − 06	1.3832e − 06
2−12	2.1042e − 04	6.2022e − 05	1.7905e − 05	5.0708e − 06	1.3832e − 06
2−14	2.1042e − 04	6.2022e − 05	1.7905e − 05	5.0708e − 06	1.3832e − 06
2−16	2.1042e − 04	6.2022e − 05	1.7905e − 05	5.0708e − 06	1.3832e − 06
EN,M	2.1042e − 04	6.2022e − 05	1.7905e − 05	5.0708e − 06	1.3832e − 06
RN,M	1.7624	1.7924	1.8201	1.8742	—

Table 3 E^N,M and R^N,M for Example 3 with $\delta =\mathrm{0.60.5}\sqrt{\epsilon },\eta =\sqrt{\epsilon }$.

ε↓	N = 64	N = 128	N = 256	N = 512	N = 1024
M = 8	M = 16	M = 32	M = 64	M = 128
2−6	3.7701e − 03	1.9482e − 03	1.0119e − 03	5.2239e − 04	2.5485e − 04
2−8	3.7721e − 03	1.8582e − 03	1.0195e − 03	5.2629e − 04	2.5625e − 04
2−10	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−12	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−14	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−16	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−18	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−20	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−22	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
2−24	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
EN,M	3.7821e − 03	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
RN,M	0.9423	0.9512	0.9534	0.95672	—

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Table 4 Comparison of E^N,M and R^N,M for Example 1 with different schemes.

ε↓	N = 64	N = 128	N = 256	N = 512
M = 120	M = 240	M = 480	M = 960
Present method
EN,M	4.1053e − 03	1.9873e − 03	1.0309e − 03	5.1582e − 04
RN,M	1.0467	0.9469	0.9990	—
Upwind method in [33]
EN,M	1.6716e − 02	9.2022e − 03	4.9863e − 03	2.6885e − 03
RN,M	0.8612	0.8840	0.8912	—
Midpoint upwind in [33]
EN,M	1.0729e − 02	6.5942e − 03	3.8199e − 03	2.1180e − 03
RN,M	0.7022	0.7877	0.8508	—
Cubic B-spline in [33]
EN,M	7.5020e − 03	4.4966e − 03	2.4450e − 053	1.2728e − 03
RN,M	0.7384	0.8791	0.9418	—

Table 5 Comparison of E^N,M and R^N,M for Example 3 with different schemes.

ε↓	N = 64	N = 128	N = 256	N = 512
M = 120	M = 240	M = 480	M = 960
Present method
EN,M	1.9682e − 03	1.0179e − 03	5.2569e − 04	2.7085e − 04
RN,M	0.9512	0.9534	0.95672	—
The method in [32]
EN,M	4.4966e − 03	2.4450e − 03	1.2728e − 053	6.4609e − 03
RN,M	0.8791	0.9718	0.9715	—
The method in [35]
EN,M	5.4975e − 03	3.2062e − 03	1.7716e − 03	9.4594e − 04
RN,M	0.7779	0.8558	0.9053	—
The method in [36]
EN,M	4.8801e − 03	2.4414e − 03	1.2211e − 03	6.1064e − 04
RN,M	0.9492	0.9455	0.9428	—

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

7. Conclusion

In this study, a parameter-uniformly convergent numerical scheme has been developed for solving singularly perturbed differential equations with mixed shift parameters. The method employs the Crank–Nicolson scheme for time discretization and an exponentially fitted finite difference scheme for spatial discretization on uniform meshes [49]. Theoretical analysis and numerical experiments confirm that the scheme attains second-order uniform convergence with respect to the perturbation parameter ε.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Funding

No funding was received for this study.

References

1 Baker C. T., Bocharov G. A., and Rihan F. A., A Report on the Use of Delay Differential Equations, Numerical Modelling in the Biosciences, 1999, 343, Manchester Centre for Computational Mathematics, Manchester, UK.

2 Rao R. N. and Chakravarthy P. P., 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. (2014) 7, no. 1, 23–40, https://doi.org/10.4208/nmtma.2014.1316nm, 2-s2.0-84893241620.

3 Buedo-Fernández S. and Liz E., On the Stability Properties of a Delay Differential Neoclassical Model of Economic Growth, Electronic Journal of Qualitative Theory of Differential Equations. (2018) 2018, no. 43, 1–14, https://doi.org/10.14232/ejqtde.2018.1.43, 2-s2.0-85050141402.

4 KuangY., Delay Differential Equations, 1993, Academic Press, New York.

5 Stein R. B., Some Models of Neuronal Variability, Biophysical Journal. (1967) 7, no. 1, 37–68, https://doi.org/10.1016/s0006-3495(67)86574-3, 2-s2.0-84954086574.

6 Bellen A. and Zennaro M., Numerical Methods for Delay Differential Equations, 2013, Numerical Mathematics and Science.

7 Cooke K., Van den Driessche P., and Zou X., Interaction of Maturation Delay and Nonlinear Birth in Population and Epidemic Models, Journal of Mathematical Biology. (1999) 39, no. 4, 332–352, https://doi.org/10.1007/s002850050194, 2-s2.0-0033209987.

8 Wu J., Theory and Applications of Partial Functional Differential Equations, 2012, Springer Science and Business Media.

9 Liu Q., Wang X., and De Kee D., Mass Transport Through Swelling Membranes, International Journal of Engineering Science. (2005) 43, no. 19-20, 1464–1470, https://doi.org/10.1016/j.ijengsci.2005.05.010, 2-s2.0-27844439487.

10 Aziz I. and Amin R., Numerical Solution of a Class of Delay Differential and Delay Partial Differential Equations via Haar Wavelet, Applied Mathematical Modelling. (2016) 40, no. 23-24, 10286–10299, https://doi.org/10.1016/j.apm.2016.07.018, 2-s2.0-84994813247.

11 Maugin G. A., StraughanB., Book Review: Heat Waves, Springer Series in Applied Mathematical Sciences. (2014) Springer, New York.

12 Kumar D., A Parameter-Uniform Scheme for the Parabolic Singularly Perturbed Problem With a Delay in Time, Numerical Methods for Partial Differential Equations. (2021) 37, no. 1, 626–642, https://doi.org/10.1002/num.22544.

13 Kumar S., Aakansha, Singh J., and Ramos H., Parameter-Uniform Convergence Analysis of a Domain Decomposition Method for Singularly Perturbed Parabolic Problems With Robin Boundary Conditions, Journal of Applied Mathematics and Computing. (2023) 69, no. 2, 2239–2261, https://doi.org/10.1007/s12190-022-01832-w.

14 Amiraliyev G. M. and Duru H., A Note on a Parameterized Singular Perturbation Problem, Journal of computational and applied mathematics. (2005) 182, no. 1, 233–242, https://doi.org/10.1016/j.cam.2004.11.047, 2-s2.0-19644389078.

15 Demsie A. W., Tiruneh A. A., and Tsega E. G., A Fitted Linear Multistep Approach for Singularly Perturbed Parabolic-Type Reaction–Diffusion Problems Using Shishkin Meshes, International Journal of Mathematics and Mathematical Sciences. (2025) 2025, no. 1, https://doi.org/10.1155/ijmm/5553428.

16 Kadalbajoo M. K. and Gupta V., A Brief Survey on Numerical Methods for Solving Singularly Perturbed Problems, Applied Mathematics and Computation. (2010) 217, no. 8, 3641–3716, https://doi.org/10.1016/j.amc.2010.09.059, 2-s2.0-78649971331.

17 Ejere A. H., Duressa G. F., Woldaregay M. M., and Dinka T. G., An Exponentially Fitted Numerical Scheme via Domain Decomposition for Solving Singularly Perturbed Differential Equations With Large Negative Shift, Journal of Mathematics. (2022) 2022, no. 1, https://doi.org/10.1155/2022/7974134.

18 Ayele M. A., Tiruneh A. A., and Derese G. A., Uniformly Convergent Scheme for Singularly Perturbed Space Delay Parabolic Differential Equation With Discontinuous Convection Coefficient and Source Term, Journal of Mathematics. (2022) 2022, no. 1, https://doi.org/10.1155/2022/1874741.

19 Bullo T. A., Accelerated Fitted Mesh Scheme for Singularly Perturbed Turning Point Boundary Value Problems, Journal of Mathematics. (2022) 2022, no. 1, https://doi.org/10.1155/2022/3767246.

20 Ejere A. H. and Feyissa T. T., A Novel Fitted Mesh Scheme Involving Caputo–Fabrizio Approach for Singularly Perturbed Fractional-Order Differential Equations With Large Negative Shift, Chaos, Solutions and Fractals: X. (2025) 15, https://doi.org/10.1016/j.csfx.2025.100128.

21 Tsega E. G., Numerical Solution of the two-dimensional Nonlinear Schrödinger Equation Using an alternat-ing Direction Implicit Method, Computational Methods for Differential Equations. (2025) 13, no. 3, 1012–1021.

22 Ranjan R. and Prasad H. S., A Novel Approach for the Numerical Approximation to the Solution of Singularly Perturbed differential-difference Equations With Small Shifts, Journal of Applied Mathematics and Computing. (2021) 65, no. 1, 403–427, https://doi.org/10.1007/s12190-020-01397-6.

23 Daba I. T., Melesse W. G., Gelu F. W., and Kebede G. D., Numerical Investigation of Singularly Perturbed Time Lag Parabolic differential-difference Equations, Heliyon. (2025) 11, no. 1, https://doi.org/10.1016/j.heliyon.2024.e41215.

24 Vivek K. and Nageshwar Rao R., Numerical Algorithms Based on Splines for Singularly Perturbed Parabolic Partial Differential Equations With Mixed Shifts, Boundary Value Problems. (2024) 2024, no. 1, 1–21, https://doi.org/10.1186/s13661-024-01971-y.

25 Mohye M. A., Munyakazi J. B., and Dinka T. G., A Nonstandard Fitted Operator Finite Difference Method for Two-Parameter Singularly Perturbed Time-Delay Parabolic Problems, Frontiers in Applied Mathematics and Statistics. (2023) 9, https://doi.org/10.3389/fams.2023.1222162.

26 Munyakazi J. B. and Patidar K. C., Limitations of Richardson’s Extrapolation for a High Order Fitted Mesh Method for Self-Adjoint Singularly Perturbed Problems, Journal of Applied Mathematics and Computing. (2010) 32, no. 1, 219–236, https://doi.org/10.1007/s12190-009-0245-6, 2-s2.0-84867978164.

27 Tefera D. M., Tiruneh A. A., and Derese G. A., Fitted Operator on the crank-nicolson Scheme for Solving a Small Time Delayed convection-diffusion Equations, Journal of applied mathematics and informatics. (2022) 40, no. 3-4, 491–505.

28 Woldaregay M. M. and Duressa G. F., Uniformly Convergent Hybrid Numerical Method for Singularly Perturbed Delay Convection-Diffusion Problems, International Journal of Differential Equations. (2021) 2021, no. 1, 6654495–20, https://doi.org/10.1155/2021/6654495.

29 Tiruneh A. A., Derese G. A., and Tefera D. M., A Nonstandard Fitted Operator Method for Singularly Perturbed Parabolic Reaction-Diffusion Problems With a Large Time Delay, International Journal of Mathematics and Mathematical Sciences. (2022) 2022, no. 1, 5625049–11, https://doi.org/10.1155/2022/5625049.

30 Hassen Z. I. and Duressa G. F., Parameter Uniform Finite Difference Formulation With Oscillation Free for Solving Singularly Perturbed Delay Parabolic Differential Equation via Exponential Spline, BMC Research Notes. (2025) 18, no. 1, https://doi.org/10.1186/s13104-024-07005-1.

31 Ejere A. H., Duressa G. F., Woldaregay M. M., and Dinka T. G., A Uniformly Convergent Numerical Scheme for Solving Singularly Perturbed Differential Equations With Large Spatial Delay, SN Applied Sciences. (2022) 4, no. 12, https://doi.org/10.1007/s42452-022-05203-9.

32 Ramesh V. P. and Kadalbajoo M. K., Upwind and Midpoint Upwind Difference Methods for Time-dependent Differential Difference Equations With Layer Behavior, Applied Mathematics and Computation. (2008) 202, no. 2, 453–471, https://doi.org/10.1016/j.amc.2007.11.033, 2-s2.0-48049121948.

33 Kumar D. and Kadalbajoo M. K., A parameter-uniform Numerical Method for Time-Dependent Singularly Perturbed Differential-Difference Equations, Applied Mathematical Modelling. (2011) 35, no. 6, 2805–2819, https://doi.org/10.1016/j.apm.2010.11.074, 2-s2.0-79951953981.

34 Bansal K. and Sharma K. K., Parameter Uniform Numerical Scheme for Time Dependent Singularly Perturbed Convection-Diffusion-Reaction Problems With General Shift Arguments, Numerical Algorithms. (2017) 75, no. 1, 113–145, https://doi.org/10.1007/s11075-016-0199-3, 2-s2.0-84986275182.

35 Gupta V., Kumar M., and Kumar S., Higher Order Numerical Approximation for Time Dependent Singularly Perturbed Differential-Difference Convection-Diffusion Equations, Numerical Methods for Partial Differential Equations. (2018) 34, no. 1, 357–380, https://doi.org/10.1002/num.22203, 2-s2.0-85030149866.

36 Woldaregay M. M. and Duressa G. F., Parameter Uniform Numerical Method for Singularly Perturbed Parabolic Differential Difference Equations, Journal of the Nigerian Mathematical Society. (2019) 38, no. 2, 223–245.

37 Woldaregay M. M. and Duressa G. F., Robust Numerical Method for Singularly Perturbed Parabolic Differential Equations With Negative Shifts, Filomat. (2021) 35, no. 7, 2383–2401, https://doi.org/10.2298/fil2107383w.

38 Woldaregay M. M. and Duressa G. F., Uniformly Convergent Numerical Scheme for Singularly Perturbed Parabolic Pdes With Shift Parameters, Mathematical Problems in Engineering. (2021) 2021, no. 1, 6637661–15, https://doi.org/10.1155/2021/6637661.

39 Friedman A., Partial Differential Equations of Parabolic Type, 2008, Courier Dover Publications.

40 Roos H. G., Robust Numerical Methods for Singularly Perturbed Differential Equations, 2008, Springer.

41 Clavero C. and Gracia J. L., A Higher Order Uniformly Convergent Method With Richardson Extrapolation in Time for Singularly Perturbed Reaction–Diffusion Parabolic Problems, Journal of computational and applied mathematics. (2013) 252, 75–85, https://doi.org/10.1016/j.cam.2012.05.023, 2-s2.0-84899824210.

42 Demsie A. W., Tiruneh A. A., and Tsega E. G., A Uniformly Convergent Scheme for Singularly Perturbed Unsteady Reaction-Diffusion Problems, Journal of Applied Mathematics. (2025) 2025, no. 1, https://doi.org/10.1155/jama/4603501.

43 O’Malley R. E., Singular Perturbation Methods for Ordinary Differential Equations, 1991, Springer-Verlag, New York.

44 Kellogg R. B. and Tsan A., Analysis of Some Difference Approximations for a Singular Perturbation Problem Without Turning Points, Mathematics of Computation. (1978) 32, no. 144, 1025–1039, https://doi.org/10.1090/s0025-5718-1978-0483484-9, 2-s2.0-84966216257.

45 Demsie A. W., Tiruneh A. A., and Tsega E. G., An Accelerated Fourth-Order Fitted Numerical Approach for Singularly Perturbed Reaction-Diffusion Problems of Small Time Lag, Palestine Journal of Mathematics. (2025) 14, no. 3.

46 Qi X., Mori K. H., and Shin M. S., Double Mesh Method for Efficient Finite Difference Calculations, Journal of the Society of Naval Architects of Japan. (1989) 1989, no. 166, 1–7, https://doi.org/10.2534/jjasnaoe1968.1989.166_1.

47 Doolan E. P., Miller J. J., and Schilders W. H., Uniform Numerical Methods for Problems With Initial and Boundary Layers, 1980, Boole Press.

48 Farrell P., Hegarty A., Miller J. M., O’Riordan E., and Shishkin G. I., Robust Computational Techniques for Boundary Layers, 2000, Chapman and hall/CRC.

49 Clavero C. and Gracia J. L., On the Uniform Convergence of a Finite Difference Scheme for Time Dependent Singularly Perturbed Reactiondiffusion Problems, Applied Mathematics and Computation. (2010) 216, no. 5, 1478–1488, https://doi.org/10.1016/j.amc.2010.02.050, 2-s2.0-77950858865.