Numerical integration method for two-parameter singularly perturbed time delay parabolic problem

This study presents an (ε, μ)−uniform numerical method for a two-parameter singularly perturbed time-delayed parabolic problems. The proposed approach is based on a fitted operator finite difference method. The Crank–Nicolson method is used on a uniform mesh to discretize the time variables initially. Subsequently, the resulting semi-discrete scheme is further discretized in space using Simpson's 1/3rd rule. Finally, the finite difference approximation of the first derivatives is applied. The method is unique in that it is not dependent on delay terms, asymptotic expansions, or fitted meshes. The fitting factor's value, which is used to account for abrupt changes in the solution, is calculated using the theory of singular perturbations. The developed scheme is demonstrated to be second-order accurate and uniformly convergent. The proposed method's applicability is validated by three model examples, which yielded more accurate results than some other methods found in the literature.

1 Introduction

This study examines a 1−D, singularly perturbed delay initial boundary value problem (IBVP) with two small parameters. Defining $\overline{\Delta }=\Delta \cup \partial \Delta$, where Δ = Δ_x = (0, 1) × Δ_t = (0, T] and ∂Δ = ∧_b ∪ ∧_l ∪ ∧_r with ∧_b = [0, 1] × [−γ, 0], ∧_l = {0} × (0, T], and ∧_r = {1} × (0, T], we consider

$\begin{array}{ll}\{\begin{array}{l}{𝕃}_{\epsilon ,\mu }w\equiv \epsilon \frac{{\partial }^{2}w}{\partial x^2}(x,t)+\mu a(x,t)\frac{\partial w}{\partial x}(x,t)-b(x,t)w(x,t)-\frac{\partial w}{\partial t}(x,t)=\\ -c(x,t)w(x,t-\gamma )+\varpi (x,t),(x,t)\in {\Delta }_x,\\ w{|}_{{\displaystyle \underset{b}{\wedge }}}={{\rm Y}}_b(x,t),\\ w{|}_{{\displaystyle \underset{l}{\wedge }}}={{\rm Y}}_l(t),w{|}_{{\displaystyle \underset{r}{\wedge }}}={{\rm Y}}_r(t),\end{array}& (1)\end{array}$

where 0 < ε ≪ 1 and 0 ≤ μ ≤ 1 are the perturbation parameters. The coefficients a(x, t), b(x, t), c(x, t), ϖ(x, t), Υ_b(x, t), Υ_l(t), and Υ_r(t) are presumed to be sufficiently smooth functions of x and t for $(x,t)\in \overline{\Delta }$ satisfying:

$0<\alpha \le a(x,t),0<\beta \le b(x,t),0<\varsigma \le c(x,t),(x,t)\in \overline{\Delta }.$

At the corner points, (0, 0), (0, 1), (−γ, 0), and (−γ, 1), the regularity and compatibility conditions are

$\begin{array}{l}\begin{array}{l}w(0,0)={{\rm Y}}_l(0),w(1,0)={{\rm Y}}_r(0),w(0,-\gamma )=\\ {{\rm Y}}_l(\gamma ),w(1,-\gamma )={{\rm Y}}_r(-\gamma ),\\ \epsilon {({{\rm Y}}_b)}_{xx}(0,0)+\mu a(0,0){({{\rm Y}}_b)}_x(0,0)-b(0,0)({{\rm Y}}_b)(0,0)\\ -{({{\rm Y}}_b)}_t(0,0)=-c(0,0)({{\rm Y}}_b)(0,-\gamma )+\varpi (0,0),\\ \epsilon {({{\rm Y}}_b)}_{xx}(1,0)+\mu a(1,0){({{\rm Y}}_b)}_x(1,0)-b(1,0)({{\rm Y}}_b)(1,0)\\ -{({{\rm Y}}_b)}_t(1,0)=-c(1,0)({{\rm Y}}_b)(1,-\gamma )+\varpi (1,0),\end{array}\end{array}$

for Δ = (0, 1) × (0, T], and Υ_b(x, t) (initial-boundary data) satisfies appropriate compatibility criteria at the two corners, (0, 0) and (1, 0). The problem shown in Equation (1) has a unique solution in the domain under consideration based on the assumptions mentioned above. To establish the existence of continuous solution $w(x,t)\in C^{1+d/2,2+d}(\overline{\Delta })$ with d ∈ (0, 1), the assumptions of Hölder continuity are made on all the coefficients and IBVPs involved in Equation (1).

Having a time delay γ > 0, in Equation (1), the entire time domain [−γ, kγ] for k ∈ ℕ is divided into several subdomains. For the analysis purpose, k = 2 is used. Because, in the region, [0, 1] × [0, γ] the solution depends on the solution in [0, 1] × [−γ, 0], which is known, and hence, the analysis is similar to any singularly perturbed partial differential equations (SPPDEs) having no time delay [1, 2]. Moreover, our analysis is conducted for the interval [0, 1] × [γ, 2γ], where w(x, t − γ) is the solution in [0, 1] × [0, γ], which is known.

The modeling of many physical, biological, and chemical systems, including population dynamics, problems in control theory, epidemiology, circadian rhythms, respiratory system, chemostat models, and tumor growth, frequently involves SPPDEs with delay. We can incorporate some historical behavior into these models due to the delay terms, which help us create more useful models for the phenomena. Time delays, for instance, are crucial to the processes of transcription and translation; in population ecology, they represent the hatching period or duration of gestation; in control systems, delay terms account for the time delay in actuation and information transmission and processing (see refs. [3–9]).

SPPDEs connect an unknown function with its derivatives evaluated at the same time. Compared with convectional instantaneous SPPDEs, SPPDEs with delay term offer more realistic models for phenomena in many scientific domains that exhibit time-lag or after-effect. Several authors have studied SPPDEs without delay in detail (see refs. [10–17]). However, numerical methods for SPPDEs with delay are still in their early stages of development. In recent years, several numerical methods have been developed to solve SPPDEs with delay term (see refs. [1, 4, 6, 7, 18–21]).

The models for two-parameter singularly perturbed problems (TP-SPPs) can be observed in the chemical reactor theory [22], transport phenomena appearing in chemistry and biology [23], and lubrication theory [12].

When the parameter μ = 1, the considered problem belongs to the class of one-parameter time-delayed convection-diffusion problem, and a boundary layer with width O(ε) arises near x = 0. When μ = 0, Equation (1) comes under the category of reaction-diffusion time-delayed problem, and two boundary layers, each with width $O(\sqrt{\epsilon })$, are observed at both ends x = 0 and x = 1. When the parameter μ ≠ 0, 1, we have two-parameter time-delayed SPPDEs for which two distinct cases, μ²/ε << 1 and μ²/ε >> 1, are contributed by the ratio of μ² to ε and two boundary layers are observed at x = 0 and x = 1 (see refs. [17, 24, 25]).

First, we refer to some earlier studies on TP-SPPs. Jha and Kadalbajoo [26] used the implicit Euler method for time discretization and upwind scheme on Shishkin–Bakhvalov mesh for spatial discretization to develop a finite difference method (FDM) for a time-dependent singularly perturbed convection-diffusion problem involving two small parameters. Shivhare et al. [10] developed second-order accurate uniformly convergent schemes using quadratic B-splines on graded mesh. To find more literature studies on this, readers may refer to refs. [11, 14, 17, 22, 24–29].

Recently, certain scholars studied two-parameter singularly perturbed time-delay parabolic problems. Ayele et al. [6] examined a numerical solution of a time-delayed TP-SPPs. The authors developed the scheme by discretizing the temporal variable on a uniform mesh using the θ− method, while a cubic spline scheme is applied by introducing a fitting factor for the spatial variable discretization. First-order parameter-uniform and nearly second-order accurate results were developed by Sumit et al. [18] in the time and space directions, respectively. Negero [4] studied TP-SPPs with time delay using the Crank–Nicolson method for time variables and the cubic splines method for spatial variables, with a fitted operator FDM. Singh et al. [7] used the Crank–Nicolson scheme in the temporal direction and B-spline collocation in the spatial direction to generate parameter-uniform numerical solutions for Equation (1).

The two non-classical FDMs used for solving SPPs are fitted mesh and fitted operator methods. The non-standard fitted operator FDMs are designed to address numerical instabilities and chaotic behavior that often affect many numerical techniques. These methods are based on the principle of dynamical consistency, which can vary depending on the specific system being analyzed [30–32]. Fitted mesh FDMs utilize non-uniform meshes, and due to the design of the mesh grid, the order of convergence is typically affected by a logarithmic factor [33]. The schemes in the studies by Singh et al. [7] and Sumit et al. [18] need a previous knowledge about the location and width of the boundary layer(s), which might be difficult to understand for beginner researchers. The exponentially fitted operator method has gained popularity as a powerful technique to solve singularly perturbed time delay parabolic problem [33, 34].

The aforementioned studies serve as our inspiration for proposing and analyzing the parameter-uniform numerical solution for TP-SPPs with time delay. As far as the author is aware, no numerical method described in the literature incorporates numerical integration (specifically, Simpson's 1/3rd rule) for solving time-delayed two-parameter singularly perturbed parabolic problems. Exponential fitted numerical integration for a class of ordinary differential equations-based quasilinear SPPs is presented in the study by Alam et al. [35], Ranjan and Prasad [36], and Reddy and Reddy [37]. In this study, we employ the Crank–Nicolson method on a uniform mesh to discretize the time variables at the first step. Subsequently, we further discretize these semi-discrete problems in space using Simpson's 1/3rd rule, and finally, we utilize the finite difference approximation of the first derivatives. This yielded a scheme with three recurrence relation.

The objective of this study is to improve accuracy by using the Crank–Nicolson scheme for temporal direction and 1/3rd Simpson's rule for space direction, with the creation of fitted operator FDM. The developed scheme is straightforward yet novel, as it offers a more accurate solution using a uniform mesh, in contrast to previous research that employed a piecewise uniform mesh (Shishkin mesh). The main benefit of the proposed scheme is that it does not rely on an extremely fine mesh or any other parameters, such as transition or delay arguments. Furthermore, the original problem does not need to be reduced to a first-order IBVP through asymptotic expansion. Compared with previous methods, the current approach is a computationally advanced integration scheme and produces better results in terms of accuracy.

The rest of the study is organized as follows: We examine properties of the continuous problem in Section 2. The semi-discretized and fully discretized scheme of the continuous problem are presented in Section 3. In Section 4, the parameter-uniform convergence analysis is presented. Section 5 presents the numerical results, and Section 6 summarizes the main conclusions of the study.

2 Solution bounds for continuous problem

To conduct a convergence analysis of the proposed approach, we will establish derivative bounds for the solution using the minimum principle for 𝕃_ε,μw.

Lemma 2.1. (Minimum principle). Let $\Xi (x,t)\in C^2(\Delta )\cap C^0(\overline{\Delta })$, such that Ξ(x, t) ≥ 0, ∀(x, t) ∈ ∂Δ = ∧_b ∪ ∧_l ∪ ∧_r & 𝕃_ε,μ Ξ (x, t) ≤ 0 then, $\Xi (x,t)\ge 0,\forall (x,t)\in \overline{\Delta }$.

Proof. One can get the proof in the study by Ayele et al. [6] and Mohye et al. [33].

The following lemma provides the bounds for solution of Equation (1) and its derivatives.

Lemma 2.2. The solution w(x, t) of the problem (1) is such that

$\begin{array}{ll}\left|\right|w(x,t)-{{\rm Y}}_b(x,0)\left|\right|\le C_{1}t,(x,t)\in \overline{\Delta },& (2)\end{array}$

where C₁ is a constant which does not depend on ε.

Proof. As Υ_b(x, 0) is smooth on $\overline{\Delta }$, the proof is similar to the study by Singh et al. [7] and Sumit et al. [18].

Lemma 2.3. Let w(x, t) be the solution of Equation (1), it satisfies:

$\begin{array}{ll}\left|\right|w(x,t)\left|\right|\le C_{*},\text{for any }(x,t)\in \overline{\Delta }.& (3)\end{array}$

Proof. From Lemma (2.2), we have ||w(x, t)−Υ_b(x, 0)|| ≤ C₁t ≤ C₁T. It follows from Equations (2) and (3) that

$\begin{array}{l}\begin{array}{l}\left|\right|w(x,t)\left|\right|=\left|\right|w(x,t)-{{\rm Y}}_b(x,0)+{{\rm Y}}_b(x,0)\left|\right|,\\ \le \left|\right|w(x,t)-{{\rm Y}}_b(x,0)\left|\right|+\left|\right|{{\rm Y}}_b(x,0)\left|\right|,\\ \le C_{1}T+{\displaystyle \underset{x\in \left[0,1\right]}{sup}}\left|\right|{{\rm Y}}_b(x,0)\left|\right|\le C_{*}.\end{array}\end{array}$

Lemma 2.4. (Uniform Stability Estimate) Let w(x, t) be the solution to the continuous problem (1), it satisfies the bound

$\begin{array}{c}\left|\right|w(x,t)|{|}_{\overline{\Delta }}\le |{\zeta }^{-1}\left|\right|\left|\varpi (x,t)\right||\\ +max\left\{\left|\right|{{\rm Y}}_b(x,t)\left|\right|,\left|\right|{{\rm Y}}_l(t)\left|\right|,\left|\right|{{\rm Y}}_r(t)\left|\right|\right\},\forall (x,t)\in \overline{\Delta },& (4)\end{array}$

where 0 < ζ ≤ (b(x, t)+c(x, t)) and $||·|{|}_{\overline{\Delta }}$ are used to denote maximum norm given in Equation (4) by $||w(x,t)|{|}_{\overline{\Delta }}=\underset{(x,t)\in \overline{\Delta }}{max}|w(x,t)|$.

Proof. See Negero [4] and Ayele et al. [6].

Lemma 2.5. For k, m ∈ ℤ⁺ satisfying 0 ≤ k + 2m ≤ 4, the solution w(x, t) and its derivatives of Equation (1) satisfy the bound:

$\begin{array}{ll}\left|\right|\frac{{\partial }^{k+m}w}{\partial x^k\partial t^m}\left|\right|\le \{\begin{array}{ll}\hfill C{\epsilon }^{-\frac{k}{2}}& \text{ for }{\mu }^2\le \frac{\epsilon 𝔑}{\alpha },\\ C{\left(\frac{\mu }{\epsilon }\right)}^k{\left(\frac{{\mu }^2}{\epsilon }\right)}^m& \text{ for }{\mu }^2\ge \frac{\epsilon 𝔑}{\alpha },\end{array}& (5)\end{array}$

where $𝔑\approx \underset{(x,t)\in \overline{\Delta }}{min}\frac{b(x,t)}{a(x,t)}$ and C are independent of ε and μ.

Proof. For the proof of Lemma (2.5), refer to Singh et al. [7] and Sumit et al. [18].

3 Construction of numerical methods

We first discretized the time direction using the Crank–Nicolson method with uniform step size k, which leads to a system of boundary value problem. Then, the discretization of space direction is made using numerical integration (Simpson's 1/3rd rule).

3.1 Temporal semi-discretization

For discretizing the domain ${\Lambda }_t^M=\left[0,T\right]$, we use uniform mesh with time step k as follows:

${\Lambda }_t^M=\left\{t_j=kj,j=0(1)M,t_M=T,k=\frac{T}{M}\right\},$

${\Lambda }_{\gamma }^m=\left\{t_j=jk,j=0(1)m,t_m=\gamma ,k=\frac{\gamma }{m}\right\},$

where M and m are the number of mesh points in the intervals [0, T] and [−γ, 0], respectively.

The Taylor series expansion was applied to w(x, t) at the point (x, t_j+1/2) to derive the Crank–Nicolson in the temporal direction as follows:

$\begin{array}{c}{w}^{j+1}(x)=w^{j+1/2}(x)+\frac{k}{2}\frac{d{w}^{j+1/2}(x)}{dt}+\frac{k^2}{8}\frac{d^2{w}^{j+1/2}(x)}{d{t}^2}+\\ \frac{k^3}{48}\frac{d^3{w}^{j+1/2}(x)}{d{t}^3}+...& (6)\end{array}$

$\begin{array}{c}{w}^j(x)=w^{j+1/2}(x)-\frac{k}{2}\frac{d{w}^{j+1/2}(x)}{dt}+\frac{k^2}{8}\frac{d^2{w}^{j+1/2}(x)}{d{t}^2}-\\ \frac{k^3}{48}\frac{d^3{w}^{j+1/2}(x)}{d{t}^3}+...& (7)\end{array}$

By substracting Equation (6) from Equation (7), we eliminate the term w^j+1/2(x) and we have

$\begin{array}{ll}\frac{d{w}^{j+1/2}(x)}{dt}=\frac{w^{j+1}(x)-w^j(x)}{k}+T{E}^{j+1/2}(x),& (8)\end{array}$

where the term TE^j+1/2(x) is local truncation error (LTE) of the scheme and is given by

$\begin{array}{ll}T{E}^{j+1/2}(x)=\frac{k^3}{24}\frac{d^3{w}^{j+1/2}(x)}{d{t}^3}+\text{Higher order terms}& (9)\end{array}$

which is order three. The semi-discretization of Equation (1) is given below:

$\begin{array}{l}\frac{\epsilon }{2}\left[w_{xx}^{j+1}+w_{xx}^j\right]+\frac{\mu }{2}\left[a^{j+1}(x)w_x^{j+1}+a^j(x)w_x^j\right]-\\ \frac{1}{2}\left[b^{j+1}(x)w^{j+1}(x)+b^j(x)w^j(x)\right]=\frac{w^{j+1}(x)-w^j(x)}{k}+\\ T{E}^{j+1/2}(x)+\\ \frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]-\frac{1}{2}[c^{j+1}(x)w^{j+1-m}(x)+\\ c^j(x)w^{j-m}(x)],\\ \Rightarrow \frac{\epsilon }{2}\left[w_{xx}^{j+1}+w_{xx}^j\right]\\ +\frac{\mu }{2}[a^{j+1}(x)w_x^{j+1}+{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)+{(a^{\prime }(x))}^j{w}^j(x)+\\ a^j(x)w_x^j]\\ =\frac{1}{2}\left[b^{j+1}(x)w^{j+1}(x)+b^j(x)w^j(x)\right]+\\ \frac{w^{j+1}(x)-w^j(x)}{k}+\frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]-\\ \frac{1}{2}\left[c^{j+1}(x)w^{j+1-m}(x)+c^j(x)w^{j-m}(x)\right]+\\ \frac{\mu }{2}\left[{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)+{(a^{\prime }(x))}^j{w}^j(x)\right],\\ \Rightarrow \frac{\epsilon }{2}\left[w_{xx}^{j+1}+w_{xx}^j\right]+\frac{\mu }{2}{\left[a^{j+1}(x)w^{j+1}(x)+a^j(x)w^j(x)\right]}^{\prime }=\\ \frac{1}{2}\left[b^{j+1}(x)w^{j+1}(x)+b^j(x)w^j(x)\right]+\frac{w^{j+1}(x)-w^j(x)}{k}+\\ \frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]-\\ \frac{1}{2}\left[c^{j+1}(x)w^{j+1-m}(x)+c^j(x)w^{j-m}(x)\right]+\\ \frac{\mu }{2}\left[{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)+{(a^{\prime }(x))}^j{w}^j(x)\right],\end{array}$

where (′) denotes derivative of the function. The foregoing expression is rewritten as follows:

$\begin{array}{ll}\{\begin{array}{l}{𝕃}_{\epsilon ,\mu }^M{w}^{j+1}(x):=\\ \frac{\epsilon }{2}\left[w_{xx}^{j+1}+w_{xx}^j\right]+\frac{\mu }{2}{\left[a^{j+1}(x)w^{j+1}(x)+a^j(x)w^j(x)\right]}^{\prime }-\\ \frac{1}{2}\left[b^{j+1}(x)w^{j+1}(x)+b^j(x)w^j(x)\right]=\frac{w^{j+1}(x)-w^j(x)}{k}+\\ \frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]-\\ \frac{1}{2}\left[c^{j+1}(x)w^{j+1-m}(x)+c^j(x)w^{j-m}(x)\right]+\\ \frac{\mu }{2}\left[{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)+{(a^{\prime }(x))}^j{w}^j(x)\right],\\ w^{j+1}(0)={{\rm Y}}_l(t_{j+1}),w^{j+1}(1)={{\rm Y}}_r(t_{j+1}),0\le j\le M,\\ w^{j+1}(x)={{\rm Y}}_b(x,t_{j+1}),(x,t_{j+1})\in \left[0,1\right]\times \left[-\gamma ,0\right],\end{array}& (10)\end{array}$

It is possible to further simplify Equation (10) as follows:

$\begin{array}{ll}\{\begin{array}{l}{𝕃}_{\epsilon ,\mu }^M{w}^{j+1}(x):=\frac{\epsilon }{2}{(w_{xx}(x))}^{j+1}+\\ \frac{\mu }{2}\left[{(a^{j+1}(x)w^{j+1}(x))}^{\prime }-{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)\right]\\ -\left(\frac{1}{2}{b}^{j+1}(x)+\frac{1}{k}\right)w^{j+1}(x)={𝔊}^{j+1}(x),\\ w^{j+1}(0)={{\rm Y}}_l(t_{j+1}),w^{j+1}(1)={{\rm Y}}_r(t_{j+1}),0\le j\le M,\\ w^{j+1}(x)={{\rm Y}}_b(x,t_{j+1}),(x,t_{j+1})\in \left[0,1\right]\times \left[-\gamma ,0\right],\end{array}& (11)\end{array}$

where ${𝔊}^{j+1}(x)=-\frac{\epsilon }{2}{(w_{xx}(x))}^j+$$\frac{\mu }{2}\left[{(a^{\prime }(x))}^j{w}^j(x)-{(a^j(x)w^j(x))}^{\prime }\right]+$$\left(\frac{1}{2}{b}^j(x)-\frac{1}{k}\right)w^j(x)+$$\frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]-$$\frac{1}{2}\left[c^{j+1}(x)w^{j+1-m}(x)+c^j(x)w^{j-m}(x)\right]$. By using the initial condition, we can evaluate the right hand-side of Equation (11) as follows:

$\begin{array}{l}{𝔊}^{j+1}(x)=\{\begin{array}{l}-\frac{\epsilon }{2}{(w_{xx}(x))}^j+\\ \frac{\mu }{2}\left[{(a^{\prime }(x))}^j{w}^j(x)-{(a^j(x)w^j(x))}^{\prime }\right]+(\frac{1}{2}{b}^j(x)-\frac{1}{k})\\ w^j(x)-\frac{1}{2}\left[c^{j+1}(x){{\rm Y}}_b(x,t_{j+1-m})+c^j(x){{\rm Y}}_b(x,t_{j-m})\right]+\\ \frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right],j=0(1)m,\\ -\frac{\epsilon }{2}{(w_{xx}(x))}^j+\\ \frac{\mu }{2}\left[{(a^{\prime }(x))}^j{w}^j(x)-{(a^j(x)w^j(x))}^{\prime }\right]+(\frac{1}{2}{b}^j(x)-\frac{1}{k})\\ w^j(x)-\frac{1}{2}\left[c^{j+1}(x)w(x,t_{j+1-m})+c^j(x)w(x,t_{j-m})\right]+\\ \frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right],j=m+1(1)M.\end{array}\end{array}$

Lemma 3.1. (Semi-discrete Minimum Principle). At (j + 1)th time level, let $\text{Ψ}(x,t_{j+1})\in C^2({\overline{\Delta }}_x)$ be smooth function such that Ψ(0, t_j+1) ≥ 0 and Ψ(1, t_j+1) ≥ 0, then ${𝕃}_{\epsilon ,\mu }^M\text{Ψ}(x,t_{j+1})\le 0,\forall x\in \Delta$ implies $\text{Ψ}(x,t_{j+1})\ge 0\forall x\in {\overline{\Delta }}_x$.

Proof. Let $({\chi }^{*},t_{j+1})\in \left\{(x,t_{j+1}):x\in {\overline{\Delta }}_x\right\}$ be such that $\text{Ψ}({\chi }^{*},t_{j+1})={min}_{x\in {\overline{\Delta }}_x}\text{Ψ}(x,t_{j+1})$ and let $\text{Ψ}({\chi }^{*},t_{j+1})<0$. Then, $({\chi }^{*},t_{j+1})\notin \left\{(0,t_{j+1}),(1,t_{j+1})\right\}$, also ${\text{Ψ}}_x({\chi }^{*},t_{j+1})=0$ and ${\text{Ψ}}_{xx}({\chi }^{*},t_{j+1})\ge 0$. By using this assumption and Equation (11), we arrive on ${𝕃}_{\epsilon ,\mu }^M\text{Ψ}(x,t_{j+1})>0$, which is a contradiction to the hypothesis, ${𝕃}_{\epsilon ,\mu }^M\text{Ψ}(x,t_{j+1})\le 0$. Hence, our asumption is wrong, so $\text{Ψ}({\chi }^{*},t_{j+1})\ge 0$, and this implies that $\text{Ψ}(x,t_{j+1})\ge 0,\forall x\in {\overline{\Delta }}_x$.

Definition 3.1. [19] The global and local errors of the time semi-discretization method in Equation (11) are given as follows:

$\begin{array}{c}{E}^{j+1}={\displaystyle {max}_{j+1\le T/k}}\left|\right|w^{j+1}(x)-W^{j+1}(x)\left|\right|,e^{j+1}=\\ |w^{j+1}(x)-W^{j+1}(x)|,& (12)\end{array}$

respectively. The contribution of each time step to the global error of the time semi-discretization is measured by the local error. After every time semi-discretization step, the LTE is calculated using the precise values w^j+1(x) as the initial data rather than W^j+1(x). The consistency and stability of the Crank–Nicolson method, which are required to assess the order of convergence in the time direction, are obtained using the following lemmas.

Lemma 3.2. If the solution w(x, t) of Equation (1) and its time derivatives are bounded in $(x,t)\in \overline{\Delta }=\left[0,1\right]\times \left[0,T\right]$, independent of ε, μ, N & M, then the LTE associated with Equation (10) satisfies ||e^j+1|| ≤ C(k³).

Proof. The proof is a straightforward consequence of Equations (8, 9).

Lemma 3.3. The global error estimate E^j+1 defined in Equation (12) at the (j + 1)th time level is bounded by E^j+1 ≤ C(k²), ∀j ≤ T/k, with a constant C that is independent of ε, μ, x, and t.

Proof. The proof is on Clavero and Jorge [2], Das and Mehrmann [27], and Clavero et al. [38].

To evaluate the proposed approach, one needs to know how the derivatives of the exact solution behave in the semi-discretized form of Equation (10). The characteristic equation of the homogenous part of Equation (10) is obtained after some rearrangements and is given as follows:

$\begin{array}{ll}\epsilon {\lambda }^2+\mu a^{j+1}(x)\lambda -{\beta }^{j+1}(x)=0,& (13)\end{array}$

where β^j+1(x) = b^j+1(x)+2/k, and Equation (13) has two roots λ₁ and λ₂, where

$\begin{array}{l}{\lambda }_1=-\frac{\mu a^{j+1}\left(x\right)}{2\epsilon }+\sqrt{{\left(\frac{\mu a^{j+1}\left(x\right)}{2\epsilon }\right)}^2+\frac{{\beta }^{j+1}\left(x\right)}{\epsilon }}\ge 0,\\ {\lambda }_2=-\frac{\mu a^{j+1}\left(x\right)}{2\epsilon }-\sqrt{{\left(\frac{\mu a^{j+1}\left(x\right)}{2\epsilon }\right)}^2+\frac{{\beta }^{j+1}\left(x\right)}{\epsilon }}\le 0.\end{array}$

These roots represent the boundary layer behavior of the solution in the neighborhood of x = 0 and x = 1 [4, 6, 14]. Let ${\varrho }_1=-{max}_{x\in \left[0,1\right]}{\lambda }_2(x)$ and ${\varrho }_2={min}_{x\in \left[0,1\right]}{\lambda }_1(x)$, we have three cases: (i) when $\frac{{\mu }^2}{\epsilon }\to 0$, as $\epsilon \to 0,{\varrho }_1={\varrho }_2=\sqrt{\frac{{\beta }^{j+1}}{\epsilon }}$, where 0 < β^* < β^j+1 and (ii) when $\frac{\epsilon }{{\mu }^2}\to 0$, as $\mu \to 0,{\varrho }_1={\varrho }_2=\frac{\mu a^{*}}{\epsilon }$ and ϱ₂ = 0, where, 0 < a^* < a^j+1, (iii) when ϱ₂ << ϱ₁, the solutions exhibit a stronger boundary layer along x = 0 and x = 1 [4].

Next, we give the semi-discrete bound of the solution w^j+1(x) of Equation (10).

Lemma 3.4. The solution w^j+1(x) of Equation (10) satisfies the derivative bound as follows:

$\begin{array}{c}\left|\frac{d^q{w}^{j+1}(x)}{d{x}^q}\right|\le C\left(1+{\varrho }_1^q{e}^{-\tau {\varrho }_{1}x}+{\varrho }_2^q{e}^{-\tau {\varrho }_2(1-x)}\right),q=0(1)r,r>1\\ \text{and }0<\tau <1.& (14)\end{array}$

Proof. For the proof, see [27, 28].

3.2 Spatial discretization

We discretized the spatial domain [0, 1] into N equal number of mesh elements as follows: ${\Lambda }_x^N=\left\{x_i=ih,i=0(1)N,x_0=0,x_N=1,h=1/N\right\}$. Integrating Equation (10), in the interval [x_i−1, x_i+1], i = 1(1)N − 1, we obtain as follows:

$\begin{array}{l}\frac{\epsilon }{2}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}\left[w_{xx}^{j+1}+w_{xx}^j\right]}(x)dx+\\ \frac{\mu }{2}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}{\left[a^{j+1}(x)w^{j+1}(x)+a^j(x)w^j(x)\right]}^{\prime }}dx-\\ \frac{1}{2}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}\left[b^{j+1}(x)w^{j+1}(x)+b^j(x)w^j(x)\right]}dx=\\ \frac{1}{k}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}\left[w^{j+1}(x)-w^j(x)\right]}dx+\\ \frac{1}{2}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]}dx-\\ \frac{1}{2}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}\left[c^{j+1}(x)w^{j+1-m}(x)+c^j(x)w^{j-m}(x)\right]}dx\\ +\frac{\mu }{2}{\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}\left[{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)+{(a^{\prime }(x))}^j{w}^j(x)\right]}dx,\\ \Rightarrow \frac{\epsilon }{2}\left[{\left(w_x(x_{i+1})\right)}^{j+1}-{\left(w_x(x_{i-1})\right)}^{j+1}+{\left(w_x(x_{i+1})\right)}^j-{\left(w_x(x_{i-1})\right)}^j\right]+\\ \frac{\mu }{2}[a^{j+1}(x_{i+1})w^{j+1}(x_{i+1})-a^{j+1}(x_{i-1})w^{j+1}(x_{i-1})+a^j(x_{i+1})w^j(x_{i+1})-\\ a^j(x_{i-1})w^j(x_{i-1})]={\displaystyle {\int }_{x_{i-1}}^{x_{i+1}}{\Theta }^j}(x)dx,& (15)\end{array}$

where

$\begin{array}{l}\begin{array}{l}{\Theta }^j(x)dx=\frac{1}{2}\left[b^{j+1}(x)w^{j+1}(x)+b^j(x)w^j(x)\right]+\\ \frac{1}{k}\left[w^{j+1}(x)-w^j(x)\right]+\\ \frac{1}{2}\left[{\varpi }^{j+1}(x)+{\varpi }^j(x)\right]-\\ \frac{1}{2}\left[c^{j+1}(x)w^{j+1-m}(x)+c^j(x)w^{j-m}(x)\right]+\\ \frac{\mu }{2}\left[{(a^{\prime }(x))}^{j+1}{w}^{j+1}(x)+{(a^{\prime }(x))}^j{w}^j(x)\right].\end{array}\end{array}$

In Equation (15), there is definite integral, and to approximate this integral, we use Simpson's 1/3rd formula, and by using the notation $w^j(x_i)=W_i^j$, Equation (15) is reduced to:

$\begin{array}{ll}\begin{array}{l}\frac{\epsilon }{2}\left[{(W_x)}_{i+1}^{j+1}-{(W_x)}_{i-1}^{j+1}+{(W_x)}_{i+1}^j-{(W_x)}_{i-1}^j\right]+\\ \frac{\mu }{2}\left[a_{i+1}^{j+1}{W}_{i+1}^{j+1}-a_{i-1}^{j+1}{W}_{i-1}^{j+1}+a_{i+1}^j{W}_{i+1}^j-a_{i-1}^j{W}_{i-1}^j\right]\\ =\frac{h}{3}\left[{\Theta }_{i+1}^j+4{\Theta }_i^j+{\Theta }_{i-1}^j\right].\end{array}& (16)\end{array}$

The first derivative of w^j+1(x) and w^j(x) with respect to x at the grid point x_i, i = 1(1)N − 1 is approximated by the following formula:

$\begin{array}{ll}\begin{array}{l}{(w_x)}_i^{j+1}\approx \frac{W_{i+1}^{j+1}-W_{i-1}^{j+1}}{2h},\\ {(w_x)}_{i+1}^{j+1}\approx \frac{3W_{i+1}^j-4W_i^{j+1}+W_{i-1}^{j+1}}{2h},\\ {(w_x)}_{i-1}^{j+1}\approx \frac{-W_{i+1}^{j+1}+4W_i^{j+1}-3W_{i-1}^{j+1}}{2h},\end{array}& (17)\end{array}$

The formula in Equation (17) is also used for the jth time level. Substituting Equation (17) into Equation (16) and simplifying, we arrive at the following relation:

$\begin{array}{ll}\begin{array}{l}\frac{\epsilon }{h}\left[W_{i+1}^{j+1}-2W_i^{j+1}+W_{i-1}^{j+1}+W_{i+1}^j-2W_i^j+W_{i-1}^j\right]+\\ \frac{\mu }{2}\left[a_{i+1}^{j+1}{W}_{i+1}^{j+1}-a_{i-1}^{j+1}{W}_{i-1}^{j+1}+a_{i+1}^j{W}_{i+1}^j-a_{i-1}^j{W}_{i-1}^j\right]\\ =\frac{h}{3}\left[{\Theta }_{i+1}^j+4{\Theta }_i^j+{\Theta }_{i-1}^j\right].\end{array}& (18)\end{array}$

To treat the effect of perturbation parameter, exponential fitting factor (σ(ρ)) is multiplied Equation (18) on the term containing ε as follows:

$\begin{array}{ll}\begin{array}{l}\frac{\epsilon \sigma (\rho )}{h}\left[W_{i+1}^{j+1}-2W_i^{j+1}+W_{i-1}^{j+1}+W_{i+1}^j-2W_i^j+W_{i-1}^j\right]+\\ \frac{\mu }{2}\left[a_{i+1}^{j+1}{W}_{i+1}^{j+1}-a_{i-1}^{j+1}{W}_{i-1}^{j+1}+a_{i+1}^j{W}_{i+1}^j-a_{i-1}^j{W}_{i-1}^j\right]\\ =\frac{h}{3}\left[{\Theta }_{i+1}^j+4{\Theta }_i^j+{\Theta }_{i-1}^j\right].\end{array}& (19)\end{array}$

Let $\rho =\frac{h}{\epsilon }$ and taking the limit of Equation (19) as h → 0

$\begin{array}{ll}\begin{array}{l}{\displaystyle \underset{h\to 0}{lim}}\frac{\sigma (\rho )}{\rho }\left[W_{i+1}^{j+1}-2W_i^{j+1}+W_{i-1}^{j+1}+W_{i+1}^j-2W_i^j+W_{i-1}^j\right]+\\ {\displaystyle \underset{h\to 0}{lim}}\frac{\mu }{2}\left[a_{i+1}^{j+1}{W}_{i+1}^{j+1}-a_{i-1}^{j+1}{W}_{i-1}^{j+1}+a_{i+1}^j{W}_{i+1}^j-a_{i-1}^j{W}_{i-1}^j\right]=0.\end{array}& (20)\end{array}$

3.3 Determination of fitting factor

From the theory of singular perturbations [4, 5], the solution of Equation (11) is written as follows:

$\begin{array}{c}{W}^{j+1}(x)=\\ \{\begin{array}{ll}{W}_0^{j+1}(x)+\left[W_0^{j+1}(0)-{{\rm Y}}_l^{j+1}(0)\right]exp& \left(-\frac{\mu a^{j+1}(0)}{\epsilon }(x)\right),\\ \hfill \text{at Left Layer},\\ W_0^{j+1}(x)+\left[W_0^{j+1}(1)-{{\rm Y}}_r^{j+1}(1)\right]exp& \left(\frac{\mu a^{j+1}(1)}{\epsilon }(1-x)\right),\\ \hfill \text{at Right Layer},\end{array}\end{array}$

where $W_0^{j+1}(x)$ is the solution of reduced problem (see refs. [4, 6]). Using Taylor series expansion for a(x) near the point x = 0 and x = 1 and considering h is reasonably small and evaluating the result in Equation (20) at x = x_i = h gives

$\begin{array}{c}{W}_i^{j+1}=\\ \{\begin{array}{ll}{(W_0)}_i^{j+1}+\left[W_0^{j+1}(0)-{{\rm Y}}_l^{j+1}(0)\right]exp& \left(-\frac{\mu a^{j+1}(0)}{\epsilon }(ih)\right),\\ \hfill \text{left layer},\\ {(W_0)}_i^{j+1}+\left[W_0^{j+1}(1)-{{\rm Y}}_r^{j+1}(1)\right]exp& \left(\frac{\mu a^{j+1}(1)}{\epsilon }(1-ih)\right),\\ \hfill \text{right layer}.\end{array}& (21)\end{array}$

For left layer, substituting Equation (21) into Equation (20) and taking $\rho =\frac{h}{\epsilon }$

$\begin{array}{l}\begin{array}{l}\frac{\sigma }{\rho }\left[e^{(-\mu a^{j+1}(0)\rho )}-2+e^{(\mu a^{j+1}(0)\rho )}\right]=-\frac{\mu }{2}{a}^{j+1}(0)\\ \left[e^{(-\mu a^{j+1}(0)\rho )}+e^{(\mu a^{j+1}(0)\rho )}\right],\\ \Rightarrow -\frac{\sigma }{\rho }{\left[e^{\left(\frac{\mu a^{j+1}(0)\rho }{2}\right)}-e^{\left(-\frac{\mu a^{j+1}(0)\rho }{2}\right)}\right]}^2=\\ \frac{-\mu }{2}{a}^{j+1}(0)\left[e^{\left(\frac{\mu a^{j+1}(0)\rho }{2}\right)}-e^{\left(-\frac{\mu a^{j+1}(0)\rho }{2}\right)}\right],\end{array}\end{array}$

After some simplification of the above expression, we have

$\begin{array}{ll}\sigma (\rho )=\frac{\rho \mu a^{j+1}(0)}{2}coth\left(\frac{\mu a^{j+1}(0)\rho }{2}\right).& (22)\end{array}$

We use the same process for the left layers to create the right layer as well. After computing σ(ρ) and taking limits on both sides as h → 0, we have:

$\begin{array}{ll}\sigma (\rho )=\{\begin{array}{ll}\frac{\rho \mu a^{j+1}(0)}{2}coth\left(\frac{\mu a^{j+1}(0)\rho }{2}\right),& \text{ for }x=0,\\ \frac{\rho \mu a^{j+1}(1)}{2}coth\left(\frac{\mu a^{j+1}(1)\rho }{2}\right),& \text{ for }x=1.\end{array}& (23)\end{array}$

where ρ = h/ε. Substituting Equation (22) into Equation (19) and solving for ${\Theta }_i^j$, we evaluate at:

$\begin{array}{ll}\{\begin{array}{l}\left[\frac{\sigma (\rho )}{\rho }-\frac{\mu }{2}{a}_{i-1}^{j+1}-\frac{h}{6}{b}_{i-1}^{j+1}-\frac{h}{3k}-\frac{h}{6}\mu {(a^{\prime })}_{i-1}^{j+1}\right]W_{i-1}^{j+1}+\\ \left[-\frac{2\sigma (\rho )}{\rho }-\frac{2h}{3}{b}_i^{j+1}+\frac{4h}{3k}-\frac{2h}{3}\mu {(a^{\prime })}_i^{j+1}\right]W_i^{j+1}+\\ \left[\frac{\sigma (\rho )}{\rho }+\frac{\mu }{2}{a}_{i+1}^{j+1}-\frac{h}{6}{b}_{i+1}^{j+1}-\frac{h}{3k}-\frac{h}{6}\mu {(a^{\prime })}_{i+1}^{j+1}\right]W_{i+1}^{j+1}=\\ \left[-\frac{\sigma (\rho )}{\rho }+\frac{\mu }{2}{a}_{i-1}^j+\frac{h}{6}{b}_{i-1}^j-\frac{h}{3k}+\frac{h}{6}\mu {(a^{\prime })}_{i-1}^j\right]W_{i-1}^j+\\ \left[\frac{2\sigma (\rho )}{\rho }+\frac{2h}{3}{b}_i^j-\frac{4h}{3k}+\frac{2h}{3}\mu {(a^{\prime })}_i^j\right]W_i^j+\\ \left[-\frac{\sigma (\rho )}{\rho }-\frac{\mu }{2}{a}_{i+1}^j+\frac{h}{6}{b}_{i+1}^j-\frac{h}{3k}+\frac{h}{6}\mu {(a^{\prime })}_{i+1}^j\right]W_{i+1}^j+\\ \frac{h}{6}\left[{\varpi }_{i+1}^{j+1}+{\varpi }_{i+1}^j\right]+\frac{2h}{3}\left[{\varpi }_i^{j+1}+{\varpi }_i^j\right]+\frac{h}{6}\left[{\varpi }_{i-1}^{j+1}+{\varpi }_{i-1}^j\right]-\\ \frac{h}{6}\left[c_{i-1}^{j+1}{W}_{i-1}^{j+1-m}+c_{i-1}^j{W}_{i-1}^{j-m}\right]-\frac{2h}{3}\left[c_i^{j+1}{W}_i^{j+1-m}+c_i^j{W}_i^{j-m}\right]\\ -\frac{h}{6}\left[c_{i+1}^{j+1}{W}_{i+1}^{j+1-m}+c_{i-1}^j{W}_{i+1}^{j-m}\right].\end{array}& (24)\end{array}$

Now, Equation (24) can be expressed as a three-term recurrence relation as follows:

$\begin{array}{ll}\{\begin{array}{l}{𝕃}_{\epsilon ,\mu }^{N,M}{W}_i^{j+1}=A_i^{-}{W}_{i-1}^{j+1}+A_i^c{W}_i^{j+1}+A_i^{+}{W}_{i+1}^{j+1}=B_i^{-}{W}_{i-1}^j+\\ B_i^c{W}_i^j+B_i^{+}{W}_{i+1}^j+F_j^j,\\ W_i^{j+1}={{\rm Y}}_b(x_i,t_{j+1}),i=1,...,N-1,-(m+1)\le j\le -1,\\ W_0^{j+1}={{\rm Y}}_l(t_{j+1}),W_{N+1}^{j+1}={{\rm Y}}_r(t_{j+1}),0\le j\le M-1.\end{array}& (25)\end{array}$

For i = 0, 1, ..., N and j = 0, 1, ..., M where

$\begin{array}{ll}\{\begin{array}{l}{A}_i^{-}=\frac{\sigma (\rho )}{\rho }-\frac{\mu }{2}{a}_{i-1}^{j+1}-\frac{h}{6}{b}_{i-1}^{j+1}-\frac{h}{3k}-\frac{h}{6}\mu {(a^{\prime })}_{i-1}^{j+1},\\ A_i^c=-\frac{2\sigma (\rho )}{\rho }-\frac{2h}{3}{b}_i^{j+1}-\frac{4h}{3k}-\frac{2h}{3}\mu {(a^{\prime })}_i^{j+1},\\ A_i^{+}=\frac{\sigma (\rho )}{\rho }+\frac{\mu }{2}{a}_{i+1}^{j+1}-\frac{h}{6}{b}_{i+1}^{j+1}-\frac{h}{3k}-\frac{h}{6}\mu {(a^{\prime })}_{i+1}^{j+1},\\ B_i^{-}=-\frac{\sigma (\rho )}{\rho }+\frac{\mu }{2}{a}_{i-1}^j+\frac{h}{6}{b}_{i-1}^j-\frac{h}{3k}+\frac{h}{6}\mu {(a^{\prime })}_{i-1}^j,\\ B_i^c=\frac{2\sigma (\rho )}{\rho }+\frac{2h}{3}{b}_i^j-\frac{4h}{3k}+\frac{2h}{3}\mu {(a^{\prime })}_i^j,\\ B_i^{+}=-\frac{\sigma (\rho )}{\rho }-\frac{\mu }{2}{a}_{i+1}^j+\frac{h}{6}{b}_{i+1}^j-\frac{h}{3k}+\frac{h}{6}\mu {(a^{\prime })}_{i+1}^j,\end{array}& (26)\end{array}$

and the right hand side is given as follows:

$\begin{array}{l}{F}_i^j=\{\begin{array}{l}-\frac{h}{6}[c_{i-1}^{j+1}{{\rm Y}}_b(x_{i-1},t_{j+1})+2c_i^{j+1}{{\rm Y}}_b(x_i,t_{j+1})+\\ c_{i+1}^{j+1}{{\rm Y}}_b(x_{i+1},t_{j+1})]\\ -\frac{h}{6}\left[c_{i-1}^j{{\rm Y}}_b(x_{i-1},t_j)+2c_i^j{{\rm Y}}_b(x_i,t_j)+c_{i+1}^j{{\rm Y}}_b(x_{i+1},t_j)\right]+\\ \frac{h}{6}\left[{\varpi }_{i-1}^{j+1}+2{\varpi }_i^{j+1}+{\varpi }_{i+1}^{j+1}\right]+\frac{h}{6}\left[{\varpi }_{i-1}^j+2{\varpi }_i^j+{\varpi }_{i+1}^j\right],\\ \text{for }{t}_j<m,j=0(1)M-1,\\ -\frac{h}{6}[c_{i-1}^{j+1}W(x_{i-1},t_{j+1-m})+2c_i^{j+1}W(x_i,t_{j+1-m})+\\ c_{i+1}^{j+1}W(x_{i+1},t_{j+1-m})]\\ -\frac{h}{6}[c_{i-1}^{j}W(x_{i-1},t_{j-m})+2c_i^{j}W(x_i,t_{j-m})+\\ c_{i+1}^{j}W(x_{i+1},t_{j-m})]+\\ \frac{h}{6}\left[{\varpi }_{i-1}^{j+1}+2{\varpi }_i^{j+1}+{\varpi }_{i+1}^{j+1}\right]+\\ \frac{h}{6}\left[{\varpi }_{i-1}^j+2{\varpi }_i^j+{\varpi }_{i+1}^j\right],\text{for }{t}_j>m,j=0(1)M-1.\end{array}\end{array}$

To solve the problem in Equation (1), the required scheme developed in Equation (25) is referred to as a fitted operator FDM obtained through Simpson's 1/3rd rule. By multiplying both sides of Equation (24) by negative, ensuring that the row sums are non-negative and the off-diagonal entries are non-positive, one can establish an M-matrix criterion using the provided coefficients. For small mesh size, it can be observed from the tridiagonal system of Equations (25, 26) that

$\begin{array}{ll}\left|A_i^c\right|\ge \left|A_i^{-}\right|+\left|A_i^{+}\right|.& (27)\end{array}$

This shows the coefficient matrix associated with ${𝕃}_{\epsilon ,\mu }^{N,M}$ is irreducibly diagonally dominant as shown in Equation (27) and non-singular. As a result, the matrix is an M−matrix, and with the specified boundary conditions, matrix inverse can be used to solve Equation (25).

4 Parameter-uniform convergence analysis

We have shown that the continuous solution and its derivatives are bounded. We can also estimate and control the errors caused by the discrete approximation in the time variable. To ensure the stability and consistency of the developed scheme, we investigate the error estimate in the spatial variable and the total discrete scheme.

Lemma 4.1. (Discrete Minimum Principle-) Assume Ξ^j+1(x) be a mesh function satisfying, Ξ(0, t_j+1) ≥ 0, Ξ(1, t_j+1) ≥ 0, and ${⋎}_{𝕃,\mu }^{N,M}\Xi (x_i,t_{j+1})\le 0$ on ${\Delta }_x^N\times {\Delta }_t^M$. Then, Ξ(x_i, t_j) ≥ 0 at each point ${\overline{\Delta }}_x^N\times {\overline{\Delta }}_t^M$, where ${\Delta }_x^N$ and ${\Delta }_t^M$ are discretized domain.

Proof. Let(x_ϑ, t_j+1) ∈ {(x_i, t_j+1):i = 1(1)N} such that ${\Xi }^{j+1}(x_{\vartheta })={min}_{1\le i\le N}{\Xi }^{j+1}(x_i)$ and let ${\Xi }^{j+1}(x_{\vartheta })<0$. From this (x_ϑ, t_j+1)∉{(0, t_j+1), (1, t_j+1)}, which implies that (x_ϑ, t_j+1) ∈ {(x_i, t_j+1):i = 1(1)N − 1}. Then, from Equation (25), we have

${𝕃}_{\epsilon ,\mu }^{N,M}{\Xi }_{\vartheta }^{j+1}=A_i^{-}{\Xi }_{\vartheta -1}^{j+1}+A_i^c{\Xi }_{\vartheta }^{j+1}+A_i^{+}{\Xi }_{\vartheta +1}^{j+1}\ge 0,$

since $|A_i^{-}|>0,|A_i^c|>0$ and $|A_i^{+}|>0$, which is a contradiction. As a result, ${\Xi }^{j+1}(x_{\vartheta })\ge 0$. Therefore, ${\Xi }^{j+1}(x_i)\ge 0,\forall (x_i,t_{j+1})\in \left\{(x_i,t_{j+1}):i=1(1)N-1\right\}$.

According to Lemma (4.1), the discrete operator ${𝕃}_{\epsilon ,\mu }^{N.M}$ follows discrete minimum principle, which results in a monotone coefficient matrix. Furthermore, because it is irreducibly diagonally dominant, it is an M− matrix. This guarantees the existence of a distinct discrete solution.

Lemma 4.2. (Estimating Discrete Uniform Stability-) The solution $W_i^{j+1}$ of the discrete scheme in Equation (25) on ${𝕃}_{\epsilon ,\mu }^{N,M}$ gratifies the following estimate:

$\begin{array}{c}\left|\right|W_i^{j+1}\left|\right|\le {\zeta }^{-1}\left|\right|{𝕃}_{\epsilon ,\mu }^{N,M}{W}_i^{j+1}\left|\right|+\\ max\left\{\left|{{\rm Y}}_l(t_{j+1})\right|,\left|{{\rm Y}}_r(t_{j+1})\right|,\left|{{\rm Y}}_b(x_i,t_{j+1})\right|\right\},& (28)\end{array}$

for i = 1(1)N − 1 and 0 < ζ ≤ b(x, t) + c(x, t).

Proof. Let $\Gamma ={\zeta }^{-1}||{𝕃}_{\epsilon ,\mu }^{N,M}{W}_i^{j+1}||+max\left\{|{{\rm Y}}_l(t_{j+1})|,|{{\rm Y}}_r(t_{j+1})|,|{{\rm Y}}_b(x_i,t_{j+1})|\right\}$ and describe the barrier functions as follows:

$\begin{array}{ll}{({▽}^{\pm })}_i^{j+1}=\Gamma \pm W_i^{j+1}.& (29)\end{array}$

The discrete function ${▽}^{\pm }(x_i,t_{j+1})$ given in Equation (29) on the boundary and interval functions is as follows:

$\begin{array}{l}\begin{array}{l}{▽}^{\pm }(0,t_{j+1})=\Gamma +max\left\{\left|{{\rm Y}}_l(t_{j+1})\right|,\left|{{\rm Y}}_r(t_{j+1})\right|,\left|{{\rm Y}}_b(0,t_{j+1})\right|\right\}\\ \pm {{\rm Y}}_l(t_{j+1})\ge 0,\\ {▽}^{\pm }(1,t_{j+1})=\Gamma +max\left\{\left|{{\rm Y}}_l(t_{j+1})\right|,\left|{{\rm Y}}_r(t_{j+1})\right|,\left|{{\rm Y}}_b(1,t_{j+1})\right|\right\}\pm \\ {{\rm Y}}_r(t_{j+1})\ge 0,\\ {▽}^{\pm }(x_i,0)=\Gamma +max\left\{\left|{{\rm Y}}_l(0)\right|,\left|{{\rm Y}}_r(0)\right|,\left|{{\rm Y}}_b(x_i,0)\right|\right\}\\ \pm {{\rm Y}}_b(x_i,0)\ge 0.\end{array}\end{array}$

Additionally, on the discretized domain, it is given as follows:

$\begin{array}{l}\begin{array}{l}{𝕃}_{\epsilon ,\mu }^{N,M}{\left({▽}^{\pm }\right)}_i^{j+1}=\Gamma \pm W_i^{j+1}=A_i^{-}\left(\Gamma \pm W_{i-1}^{j+1}\right)+\\ A_i^c\left(\Gamma \pm W_i^{j+1}\right)+A_i^{+}\left(\Gamma \pm W_{i+1}^{j+1}\right),\\ =\left(A_i^{-}+A_i^c+A_i^{+}\right)\Gamma \pm H_i^j\ge 0,\end{array}\end{array}$

where $H_i^j$ is the right hand side of Equation (25). Using Lemma (4.1), Equation (28) holds true.

Truncation error for Discrete scheme

To establish a parametric uniform convergence of the discrete scheme of Equation (25), let Ti(h) represent an LTE of the proposed discrete scheme. The LTE is given as follows:

$\begin{array}{c}{T}_i(h)=A_i^{-}{W}_{i-1}^{j+1}+A_i^c{W}_i^{j+1}+A_i^{+}{W}_{i+1}^{j+1}-\\ \left(B_i^{-}{W}_{i-1}^j+B_i^c{W}_i^j+B_i^{+}{W}_{i+1}^j+F\right)& (30)\end{array}$

To calculate the spatial truncation error, use Equation (25) at x_ι, ι = i, i±1 in Equation (30).

$\begin{array}{l}{T}_i(h)=\\ \{\begin{array}{l}\left[A_i^{-}{W}_{i-1}^{j+1}-h{(\epsilon W_{xx}^{j+1}+\mu a^{j+1}(x)W_x^{j+1}-{\beta }^{j+1}(x)W^{j+1})}_{i-1}\right]+\\ \left[A_i^c{W}_i^{j+1}-h{(\epsilon W_{xx}^{j+1}+\mu a^{j+1}(x)W_x^{j+1}-{\beta }^{j+1}(x)W^{j+1})}_i\right]+\\ \left[A_i^{+}{W}_{i+1}^{j+1}-h{(\epsilon W_{xx}^{j+1}+\mu a^{j+1}(x)W_x^{j+1}-{\beta }^{j+1}(x)W^{j+1})}_{i+1}\right]+\\ \left[B_i^{-}{W}_{i-1}^j-h{(\epsilon W_{xx}^j+\mu a^j(x)W_x^j-{\beta }^j(x)W^j)}_{i-1}\right]+\\ \left[B_i^c{W}_i^j-h{(\epsilon W_{xx}^j+\mu a^j(x)W_x^j-{\beta }^j(x)W^j)}_i\right]+\\ \left[B_i^{+}{W}_{i+1}^j-h{(\epsilon W_{xx}^j+\mu a^j(x)W_x^j-{\beta }^j(x)W^j)}_{i+1}\right]\end{array}& (31)\end{array}$

Replacing Taylor's series expansion of

$\begin{array}{ll}\{\begin{array}{l}{W}_{i-1}^{j+1}\cong W_i^{j+1}-\frac{h}{1!}\frac{d{W}_i^{j+1}}{dx}+\frac{h^2}{2!}\frac{d^2{W}_i^{j+1}}{d{x}^2}-\frac{h^3}{3!}\frac{d^3{W}_i^{j+1}}{d{x}^3}+\frac{h^4}{4!}\frac{d^4{W}_i^{j+1}}{d{x}^4}\\ +O(h^5),\\ W_{i+1}^{j+1}\cong W_i^{j+1}+\frac{h}{1!}\frac{d{W}_i^{j+1}}{dx}+\frac{h^2}{2!}\frac{d^2{W}_i^{j+1}}{d{x}^2}+\frac{h^3}{3!}\frac{d^3{W}_i^{j+1}}{d{x}^3}+\frac{h^4}{4!}\frac{d^4{W}_i^{j+1}}{d{x}^4}\\ +O(h^5),\\ W_{i-1}^j\cong W_i^j-\frac{h}{1!}\frac{d{W}_i^j}{dx}+\frac{h^2}{2!}\frac{d^2{W}_i^j}{d{x}^2}-\frac{h^3}{3!}\frac{d^3{W}_i^j}{d{x}^3}+\frac{h^4}{4!}\frac{d^4{W}_i^j}{d{x}^4}+O(h^5),\\ W_{i+1}^j\cong W_i^j+\frac{h}{1!}\frac{d{W}_i^j}{dx}+\frac{h^2}{2!}\frac{d^2{W}_i^j}{d{x}^2}+\frac{h^3}{3!}\frac{d^3{W}_i^j}{d{x}^3}+\frac{h^4}{4!}\frac{d^4{W}_i^j}{d{x}^4}+O(h^5),\\ {(W_x)}_{i-1}^{j+1}\cong {(W_x)}_i^{j+1}-\frac{h}{1!}{(W_{xx})}_i^{j+1}+\frac{h^2}{2!}{(W_{xxx})}_i^{j+1}-\frac{h^3}{3!}{(W_{xxxx})}_i^{j+1}+O(h^4),\\ {(W_x)}_{i+1}^{j+1}\cong {(W_x)}_i^{j+1}+\frac{h}{1!}{(W_{xx})}_i^{j+1}+\frac{h^2}{2!}{(W_{xxx})}_i^{j+1}+\frac{h^3}{3!}{(W_{xxxx})}_i^{j+1}+O(h^4),\\ {(W_x)}_{i-1}^j\cong {(W_x)}_i^j-\frac{h}{1!}{(W_{xx})}_i^j+\frac{h^2}{2!}{(W_{xxx})}_i^j-\frac{h^3}{3!}{(W_{xxxx})}_i^j+\\ O(h^4),\\ {(W_x)}_{i+1}^j\cong {(W_x)}_i^j+\frac{h}{1!}{(W_{xx})}_i^j+\frac{h^2}{2!}{(W_{xxx})}_i^j+\frac{h^3}{3!}{(W_{xxxx})}_i^j+\\ O(h^4),\\ {(W_{xx})}_{i-1}^{j+1}\cong {(W_{xx})}_i^{j+1}-\frac{h}{1!}{(W_{xxx})}_i^{j+1}+\frac{h^2}{2!}{(W_{xxxx})}_i^{j+1}+O(h^3),\\ {(W_{xx})}_{i+1}^{j+1}\cong {(W_{xx})}_i^{j+1}+\frac{h}{1!}{(W_{xxx})}_i^{j+1}+\frac{h^2}{2!}{(W_{xxxx})}_i^{j+1}+O(h^3),\\ {(W_{xx})}_{i+1}^j\cong {(W_{xx})}_i^j+\frac{h}{1!}{(W_{xxx})}_i^j+\frac{h^2}{2!}{(W_{xxxx})}_i^j+O(h^3),\\ {(W_{xx})}_{i-1}^j\cong {(W_{xx})}_i^j+\frac{h}{1!}{(W_{xxx})}_i^j+\frac{h^2}{2!}{(W_{xxxx})}_i^j+O(h^3),\end{array}& (32)\end{array}$

Equation (32) with Equation (31), we have

$\begin{array}{ll}\begin{array}{l}{T}_i(h)=\left[{\xi }_0^{j+1}{W}_i^{j+1}+{\xi }_0^j{W}_i^j\right]+\\ \left[{\xi }_1^{j+1}{W}_{x,i}^{j+1}+{\xi }_1^j{W}_{x,i}^j\right]+\\ \left[{\xi }_2^{j+1}{W}_{xx,i}^{j+1}+{\xi }_2^j{W}_{xx,i}^j\right]+\\ \left[{\xi }_3^{j+1}{W}_{xxx,i}^{j+1}+{\xi }_3^j{W}_{xxx,i}^j\right]+\text{higher order terms},\end{array}& (33)\end{array}$

where the coefficients in Equation (33) are given as follows:

$\begin{array}{l}\{\begin{array}{l}{\xi }_0^{j+1}=\left(A_i^{-}+A_i^c+A_i^{+}\right)+h\left({\beta }_{i-1}^{j+1}+{\beta }_i^{j+1}+{\beta }_{i+1}^{j+1}\right),\\ {\xi }_0^j=-\left(B_i^{-}+B_i^c+B_i^{+}\right)+h\left({\beta }_{i-1}^j+{\beta }_i^j+{\beta }_{i+1}^j\right),\\ {\xi }_1^{j+1}=h\left(A_i^{+}-A_i^{-}\right)-h\mu \left(a_{i-1}^{j+1}+a_i^{j+1}+a_{i+1}^{j+1}\right)-\\ h^2\left({\beta }_{i-1}^{j+1}-{\beta }_{i+1}^{j+1}\right),\\ {\xi }_1^j=h\left(B_i^{-}-B_i^{+}\right)-h\mu \left(a_{i-1}^j+a_i^j+a_{i+1}^j\right)+\\ h^2\left({\beta }_{i+1}^j-{\beta }_{i-1}^j\right),\\ {\xi }_2^{j+1}=\frac{h^2}{2}\left(A_i^{-}+A_i^{+}\right)-3h\epsilon +\\ h^2\mu \left(a_{i-1}^{j+1}-a_{i+1}^{j+1}\right)+\frac{h^3}{2}\left({\beta }_{i-1}^{j+1}+{\beta }_{i+1}^{j+1}\right),\\ {\xi }_2^j=-\frac{h^2}{2}\left(A_i^{+}+A_i^{-}\right)-3h\epsilon +\\ h^2\mu \left(a_{i-1}^j-a_{i+1}^j\right)+\frac{h^3}{2}\left({\beta }_{i-1}^j+{\beta }_{i+1}^j\right),\\ {\xi }_3^{j+1}=\frac{h^3}{6}\left(A_i^{+}-A_i^{-}\right)-\\ \frac{h^3}{2}\mu \left(a_{i-1}^{j+1}+a_{i+1}^{j+1}\right)-\\ \frac{h^4}{6}\left({\beta }_{i-1}^{j+1}-{\beta }_{i+1}^{j+1}\right),\\ {\xi }_3^j=\frac{h^3}{6}\left(B_i^{-}-B_i^{+}\right)-\\ \frac{h^3}{2}\mu \left(a_{i-1}^j+a_{i+1}^j\right)-\frac{h^4}{6}\left({\beta }_{i+1}^j-{\beta }_{i-1}^j\right).\end{array}\end{array}$

Substituting and simplifying the values of $A_i^{-},A_i^c,A_i^{+},B_i^{-},B_i^c,B_i^{+}$ and restricting the expansion of the coefficients to the first term gives ${\xi }_0^{j+1}={\xi }_0^j={\xi }_1^{j+1}={\xi }_1^j=0$ and

$\begin{array}{ll}\{\begin{array}{l}{\xi }_2^{j+1}=\frac{h^2\sigma \left(\rho \right)}{\rho }-3h\epsilon -\frac{h^3}{6}\left(11{\left(a^{\prime }\right)}_i^{j+1}+5{\beta }_i^{j+1}\right)+O\left(h^3\right),\\ {\xi }_2^j=\frac{h^2\sigma \left(\rho \right)}{\rho }-3h\epsilon -\frac{h^3}{6}\left(11{\left(a^{\prime }\right)}_i^j+5{\beta }_j^i\right)+O\left(h^3\right).\end{array}& (34)\end{array}$

Thus, by utilizing boundedness of the coefficient functions, we obtain

$\begin{array}{c}\left|T_i(h)\right|\le {\xi }_2^{j+1}{(W_{xx})}_i^{j+1}+{\xi }_2^j{(W_{xx})}_i^j\le h\epsilon (\sigma -3){(W_{xx})}_i^{j+1},\\ \text{since }\rho =\frac{h}{\epsilon }.& (35)\end{array}$

From the power series expansion of coth(ψ), we derive the following inequality:

$coth(\psi )=\frac{1}{\psi }+\frac{\psi }{3}-\frac{{\psi }^3}{45}+O({\psi }^5)\Rightarrow \psi coth(\psi )=1+\frac{{\psi }^2}{3}+O({\psi }^4),$

where $\psi =\frac{\rho \mu a^{j+1}(x_i)}{2}$. Now $|\sigma -3|=|1+\frac{{\left(\frac{\rho \mu a^{j+1}(x_i)}{2}\right)}^2}{3}-3|\le C{\rho }^2{\mu }^2$. From this, Equations (34, 35) are reduced to

$\begin{array}{ll}\left|T_i(h)\right|\le h\epsilon (C{\rho }^2{\mu }^2){(W_{xx})}_i^{j+1}=C{h}^3\frac{{\mu }^2}{\epsilon }{(W_{xx})}_i^{j+1}.& (36)\end{array}$

Thus, by applying Equation (14) given in Lemma (3.4) and considering the μ²/ε, Equation (36) reduced to:

$\begin{array}{ll}\left|T_i(h)\right|\le \{\begin{array}{ll}C{h}^3\frac{\mu }{\epsilon }& \text{ when }{\mu }^2\le \epsilon ,\\ C{h}^3\frac{{\mu }^2}{\epsilon }& \text{ when }\epsilon \le {\mu }^2.\end{array}& (37)\end{array}$

Lemma 4.3. Let the solutions to Equations (1, 25) be w(x_i, t_j+1) and ${𝕃}_{\epsilon ,\mu }^{N,M}{W}_i^{j+1}$, respectively. Next, in the spatial direction, the proposed scheme satisfies the estimation [6]:

$\begin{array}{ll}|{𝕃}_{\epsilon ,\mu }^{N,M}{W}_i^{j+1}-w(x_i,t_{j+1})|\le \{\begin{array}{ll}C{h}^2& \text{ when }{\mu }^2\le \epsilon ,\\ Ch& \text{ when }\epsilon \le {\mu }^2.\end{array}& (38)\end{array}$

Proof. For i = 1(1)N, let us construct a barrier function as follows:

${\Theta }^{\pm }(x_i,t_{j+1})=T_i(h)\pm \left({𝕃}_{\epsilon ,\mu }^{N,M}{W}_i^{j+1}-w(x_i,t_{j+1})\right).$

Now from boundary conditions, ${\Theta }^{\pm }(x_0,t_{j+1})\ge 0,{\Theta }^{\pm }(x_N,t_{j+1})\ge 0$ and ${𝕃}_{\epsilon ,\mu }^{N,M}{\Theta }^{\pm }(x_i,t_{j+1})\le 0$ for i = 1(1)N − 1. Hence, by applying Lemma (4.1) and Equations (5, 37), the required estimate can be achieved.

As of now, the above-mentioned bound yields the main result stated in the following theorem.

Theorem 4.4. Let w(x_i, t_j+1) and $W_i^{j+1}$ be solutions of the continuous problem Equation (1) and the discrete problem Equation (25), respectively. Then, the error estimate for the fully discrete scheme is as follows:

$\begin{array}{ll}{\displaystyle \underset{0\le i\le N,0\le j\le M}{max}}\left|\right|w(x_i,t_{j+1})-W_i^{j+1}\left|\right|\le \{\begin{array}{ll}C(h^2+k^2)& \text{ when }{\mu }^2\le \epsilon ,\\ C(h+k^2)& \text{ when }\epsilon \le {\mu }^2,\end{array}& (39)\end{array}$

where C is a constant unaffected by any perturbation parameter.

Proof. The proof of this theorem begins by considering the left side of Equation (39). By applying the triangular inequality using semi-discrete solution $W^{j+1}(x_i)$ and utilizing the error bounds from Lemma (4.3) and Lemma (3.3), we arrive at the result [33].

Therefore, the suggested approach reaches a second-order rate of convergence for μ²/ε ≤ 1 and shows convergence, which is independent of the perturbation parameters. Additionally, the method shows first-order convergence for ε/μ² << 1, which occurs when the spatial convergence rate takes precedence over the temporal direction as shown in Equation (38).

5 Numerical examples and results

In this section, the proposed method is implemented on three test problems, and we have demonstrated its effectiveness by comparing the outcomes with the previous finding reported in the studies by Sumit et al. [18], Negero [4], Govindarao et al. [20], and Singh et al. [7]. To support the theoretical analysis, we have provided the tabular results for errors and order of convergence. As the considered examples have no exact solution, we calculate the maximum point-wise error $E_{\epsilon ,\mu }^{N,M}$ and rate of convergence ${\overline{R}}_{\epsilon ,\mu }^{N,M}$ of the scheme using the double mesh principle given in the study by Doolan et al. [39] as follows:

$\begin{array}{c}{E}_{\epsilon ,\mu }^{N,M}={\displaystyle \underset{0\le i\le N,0\le j\le M}{max}}|W_{i,j}^{N,M}-W_{i,j}^{2N,2M}|,{\overline{R}}_{\epsilon ,\mu }^{N,M}=\\ \frac{log\left(E_{\epsilon ,\mu }^{N,M}\right)-log\left(E_{\epsilon ,\mu }^{2N,2M}\right)}{log(2)},\end{array}$

where $W_{i,j}^{N,M}$ and $W_{i,j}^{2N,2M}$ are numerical solutions computed on the mesh N × M and 2N × 2M, respectively. Additionally, the uniform rate of convergence ${\overline{R}}^{N,M}$ and uniform maximum point-wise error E^{N, M} are given in the study by Duressa and Mekonnen [17] by the following formula:

$E^{N,M}={\displaystyle \underset{\epsilon ,\mu }{max}}{E}_{\epsilon ,\mu }^{N,M},{\overline{R}}^{N,M}=\frac{log\left(E^{N,M}\right)-log\left(E^{2N,2M}\right)}{log(2)}.$

Despite the fact that the method was intended to handle both small and large time delay(γ), due to source and time constraints, we only offer examples for large delay. However, you could use either one [40]. The time delay in the examples under consideration is γ = 1.

Example 5.1. Consider the problems in the study by Negero [4] and Singh et al. [7]

$\begin{array}{l}\{\begin{array}{l}\epsilon \frac{{\partial }^{2}w}{\partial x^2}+\mu (1+x)\frac{\partial w}{\partial x}-w(x,t)-\frac{\partial w}{\partial t}=\\ -w(x,t-\gamma )+16x^2{(1-x)}^2,(x,t)\in (0,1)\times (0,2],\\ w(0,t)=0,w(1,t)=0,t\in (0,2],\\ w(x,t)=0,(x,t)\in \left[0,1\right]\times \left[-\gamma ,0\right].\end{array}\end{array}$

Example 5.2. Consider the problems in the study by Negero [4] and Singh et al. [7]

$\begin{array}{l}\{\begin{array}{l}\epsilon \frac{{\partial }^{2}w}{\partial x^2}+\mu (1+x-x^2+t^2)\frac{\partial w}{\partial x}-(1+5xt)w(x,t)-\frac{\partial w}{\partial t}=\\ -w(x,t-\gamma )+(x-x^2)(e^t-1),(x,t)\in (0,1)\times (0,2],\\ w(0,t)=0,w(1,t)=0,t\in (0,2],\\ w(x,t)=0,(x,t)\in \left[0,1\right]\times \left[-\gamma ,0\right].\end{array}\end{array}$

Example 5.3. Consider two-parameter SPPDE with delay (γ = 1) on (0, 2] × [0, 1] [7]

$\begin{array}{l}\{\begin{array}{l}\epsilon \frac{{\partial }^{2}w}{\partial x^2}+\mu (1+x^2)\frac{\partial w}{\partial x}-(1+5xt)w(x,t)-\\ \frac{\partial w}{\partial t}=-w(x,t-\gamma )+16x^2{(1-x)}^2,\\ w(0,t)=0,w(1,t)=0,t\in (0,2],\\ w(x,t)=0,(x,t)\in \left[0,1\right]\times \left[-\gamma ,0\right].\end{array}\end{array}$

The computed maximum point-wise error $E_{\epsilon ,\mu }^{N,M}$ and rate of convergence ${\overline{R}}_{\epsilon ,\mu }^{N,M}$ for examples (5.1–5.3) are presented in Tables 1–8 with different values of ε, μ, N, and M. Tables 1–5, 7 compare numerical results with those given in the studies by Negero [4], Singh et al. [7], Sumit et al. [18], and Govindarao et al. [20] by fixing or varying ε or μ. This comparison indicates that, once compared with studies found in the literature, the suggested approach provides a more accurate solution. Consequently, tabulated numerical results confirm that our theoretical estimates in Lemma (3.3) and Theorem (4.4) can be achieved by the proposed scheme, that is 1/3rd Simpson's scheme. Tables 5, 7 display the error for examples (5.2) and (5.3), respectively, for the case ε/μ² << 1 by taking a fixed value of μ and varying the values of ϵ. According to the theoretical analysis, our method is expected to be first-order convergent in space under these conditions and provides a robust validation of the method's performance. We have conducted additional numerical experiments to verify this behavior, ensuring that the spatial convergence rate dominates over the temporal direction. As we double the number of mesh points, presented in Tables 6, 8, errors become reduced, and as we approach downward in each column, errors become constant. These results confirm parameter uniformity and second-order accuracy in the relevant cases. All tables show that convergence is independent of perturbation parameters, and that the maximum absolute error decreases as the number of mesh points increases. The time taken by the CPU (in seconds) to deliver the outputs (maximum point-wise error and convergence rate) is also provided, which shows the efficiency of our proposed scheme. The numerical solution surface plots for the discussed test problems (40–42) are presented in Figures 1–3, respectively, with different values of N, M, ε, and μ. The problem displays left and right boundary layers based on the size of ε or μ, as a result of the small parameters. In addition, Figures 4–6 show the log–log plot of the maximum point-wise errors for examples (5.1), (5.2), and (5.3), respectively. Figure 4A depict the log–log plot of maximum point-wise errors for example (5.1), taking the values in Table 1. A comparison of error between the study by Negero [4] and Sumit et al. [18] and the current method using a log–log plot is shown in Figure 4B, and it can be noticed that the proposed method yields better accuracy.

Table 1

Table 1. Computed $E_{\epsilon ,\mu }^{N,M},E^{N,M},and{\overline{R}}_{\epsilon ,\mu }^{N,M}$ for example (5.1) at μ = 10⁻⁹ and comparison with Negero [4] and Sumit et al. [18].

Table 2

Table 2. Computed maximum point-wise error and rate of convergence for example (5.1) at μ = 10⁻⁴ and comparison with the study by Govindarao et al. [20] for N = M.

Table 3

Table 3. Maximum point-wise error, convergence rate, and CPU time (in seconds) for example (5.1) with ε = 10⁻⁴ and N = M.

Table 4

Table 4. Computed $E_{\epsilon ,\mu }^{N,M},E^{N,M},$ and ${\overline{R}}^{N,M}$ for example (5.1) at ε = 2⁻¹⁰ and comparison with the study by Singh et al. [7] for varying values of μ.

Table 5

Table 5. Computed $E_{\epsilon ,\mu }^{N,M},{\overline{R}}_{\epsilon ,\mu }^{N,M},E^{N,M},{\overline{R}}^{N,M}$ and CPU time (in seconds) for example (5.2) at μ = 10⁻⁹ and comparison with the study by Negero [4] and Sumit et al. [18].

Table 6

Table 6. Computed $E_{\epsilon ,\mu }^{N,M},E^{N,M},$ and ${\overline{R}}^{N,M}$ for example (5.2) at μ = 10⁻⁹ after extrapolation.

Table 7

Table 7. Maximum point-wise error, convergence rate, and CPU time (in seconds) for example (5.3) with ε = 2⁻¹⁰ and varying μ values.

Table 8

Table 8. Computed $E_{\epsilon ,\mu }^{N,M},E^{N,M},$ and ${\overline{R}}^{N,M}$ for example (5.3) at ε = 10⁻¹⁰ after extrapolation.

Figure 1

Figure 1. The numerical solution profile: (A) for N = M = 64, ε = 10⁻⁸, μ = 1, and (B) for N = M = 128, ε = 10⁻¹², and μ = 10⁻¹⁰ for example (5.1).

Figure 2

Figure 2. The numerical solution profile: (A) for N = M = 64, ε = 10⁻⁸, and μ = 1 and (B) for N = 256, M = 128, ε = 10⁻⁶, and μ = 10⁻¹⁰ for example (5.2).

Figure 3

Figure 3. Surface plots of the numerical solutions: (A) for N = M = 64, ε = 2⁻¹⁸, and μ = 1 and (B) for N = M = 128, ε = 2⁻¹⁰, and μ = 2⁻¹⁰ for example (5.3).

Figure 4

Figure 4. Log–log plot for example (5.1): (A) for Table 1 and (B) comparison of log–log plot with the existing literature.

Figure 5

Figure 5. Log–log plot for example (5.2): (A) for Tables 5, 6 and (B) comparison of log–log plot with the existing literature.

Figure 6

Figure 6. Log–log plot for example (5.3): (A) for Tables 7, 8 and (B) comparison of log–log plot with the existing literature.

For μ = 1, the considered problem belongs to the class of one-parameter time-delayed convection–diffusion parabolic problems, and a boundary layer of width O(ε) arises in the vicinity of x = 0, as shown in Figures 1A, 2A, 3A.

6 Conclusion

A two-parameter singularly perturbed time-delayed parabolic convection–diffusion problem is considered by the Crank–Nicolson method for the discretization of the time derivative, whereas in the discretization of the spatial variable, a numerical integration (1/3 rd Simpson's) formula is used by introducing an appropriate fitting factor. The method's novelty resides in its independence from delay terms, asymptotic expansions, and fitted meshes. The presence of two small parameters causes the problem to exhibit left and right boundary layers depending on the size of ε or μ. The method is demonstrated to have parameter uniform second-order convergence in both time and space. The performance of the proposed scheme is investigated by comparing the results, and it is discovered that the accuracy of the numerical results is comparable to or better than that of the existing difference schemes [4, 7, 18, 20], as verified both theoretically and numerically. To further demonstrate the versatility and applicability of the proposed approach, we highlight that the analysis method used in this study can be extended to high-order singularly perturbed time-delayed parabolic problems with Robin boundary conditions.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

SC: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. GD: Conceptualization, Investigation, Methodology, Validation, Writing – review & editing. TM: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Acknowledgments

The authors would like to thank the editor and reviewers for careful reading and giving prolific comments.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Kumar S, Kumar M. High order parameter-uniform discretization for singularly perturbed parabolic PDEs with time delay. Comput Math Appl. (2014) 68:1355–67. doi: 10.1016/j.camwa.2014.09.004

Crossref Full Text | Google Scholar

2. Clavero C, Jorge JC. An efficient and uniformly convergent scheme for one-dimensional parabolic singularly perturbed semilinear systems of reaction-diffusion type. Numer Algorithms. (2020) 85:1005–27. doi: 10.1007/s11075-019-00850-3

Crossref Full Text | Google Scholar

3. Wu J. Theory and Applications of Partial Functional Differential Equations, Vol. 119. Berlin: Springer Science & Business Media (2012).

Google Scholar

4. Negero NT. Fitted cubic spline in tension difference scheme for two-parameter singularly perturbed delay parabolic partial differential equations. Partial Differ Equ Appl Math. (2023) 8:100530. doi: 10.1016/j.padiff.2023.100530

Crossref Full Text | Google Scholar

5. Schmitt J, DiPrima R. Asymptotic methods for an infinite slider bearing with a discontinuity in film slope. J Lubr Technol. (1976) 98:446–52. doi: 10.1115/1.3452885

Crossref Full Text | Google Scholar

6. Ayele MA, Tiruneh AA, Derese GA. Fitted cubic spline scheme for two-parameter singularly perturbed time-delay parabolic problems. Results Appl Math. (2023) 18:100361. doi: 10.1016/j.rinam.2023.100361

PubMed Abstract | Crossref Full Text | Google Scholar

7. Singh S, Kumari P, Kumar D. An effective numerical approach for two parameter time-delayed singularly perturbed problems. Comput Appl Math. (2022) 41:337. doi: 10.1007/s40314-022-02046-3

Crossref Full Text | Google Scholar

8. Singh J, Kumar S, Kumar M. A domain decomposition method for solving singularly perturbed parabolic reaction-diffusion problems with time delay. Numer Methods Partial Differ Equ. (2018) 34:1849–66. doi: 10.1002/num.22256

Crossref Full Text | Google Scholar

9. Rao SCS, Kumar S. Second order global uniformly convergent numerical method for a coupled system of singularly perturbed initial value problems. Appl Math Comput. (2012) 219:3740–53. doi: 10.1016/j.amc.2012.09.075

Crossref Full Text | Google Scholar

10. Shivhare M, Podila PC, Kumar D. A uniformly convergent quadratic B-spline collocation method for singularly perturbed parabolic partial differential equations with two small parameters. J Math Chem. (2021) 59:186–215. doi: 10.1007/s10910-020-01190-7

Crossref Full Text | Google Scholar

11. Bhathawala P, Verma A. A two-parameter singular perturbation solution of one dimension flow through unsaturated porous media. Appl Math. (1975) 43:380–4.

Google Scholar

12. DiPrima RC. Asymptotic methods for an infinitely long slider squeeze-film bearing. J Lubr Technol. (1968) 90:173–83. doi: 10.1115/1.3601534

Crossref Full Text | Google Scholar

13. Jazar RN. Perturbation Methods in Science and Engineering. Berlin: Springer (2021). doi: 10.1007/978-3-030-73462-6

Crossref Full Text | Google Scholar

14. Mekonnen TB, Duressa GF. Nonpolynomial spline method for singularly perturbed parabolic problem with two small parameters. Math Probl Eng. (2023) 2023:4798517. doi: 10.1155/2023/4798517

Crossref Full Text | Google Scholar

15. Hemker P, Shishkin G, Shishkina L. ε-uniform schemes with high-order time-accuracy for parabolic singular perturbation problems. IMA J Numer Anal. (2000) 20:99–121. doi: 10.1093/imanum/20.1.99

PubMed Abstract | Crossref Full Text | Google Scholar

16. Shishkin GI. On finite difference fitted schemes for singularly perturbed boundary value problems with a parabolic boundary layer. J Math Anal Appl. (1997) 208:181–204. doi: 10.1006/jmaa.1997.5314

Crossref Full Text | Google Scholar

17. Duressa GF, Mekonnen TB. An exponentially fitted method for two parameter singularly perturbed parabolic boundary value problems. Commun Korean Math Soc. (2023) 38:299–318. doi: 10.4134/CKMS.c220020

PubMed Abstract | Crossref Full Text | Google Scholar

18. Sumit, Kumar S, Kuldeep, Kumar M. A robust numerical method for a two-parameter singularly perturbed time delay parabolic problem. Comput Appli Math. (2020) 39:209. doi: 10.1007/s40314-020-01236-1

Crossref Full Text | Google Scholar

19. Priyadarshana S, Mohapatra J, Pattanaik S. Parameter uniform optimal order numerical approximations for time-delayed parabolic convection diffusion problems involving two small parameters. Comput Appl Math. (2022) 41:233. doi: 10.1007/s40314-022-01928-w

Crossref Full Text | Google Scholar

20. Govindarao L, Sahu SR, Mohapatra J. Uniformly convergent numerical method for singularly perturbed time delay parabolic problem with two small parameters. Iran J Sci Technol Trans A Sci. (2019) 43:2373–83. doi: 10.1007/s40995-019-00697-2

PubMed Abstract | Crossref Full Text | Google Scholar

21. Gaspar F, Clavero C, Lisbona F. Some numerical experiments with multigrid methods on Shishkin meshes. J Comput Appl Math. (2002) 138:21–35. doi: 10.1016/S0377-0427(01)00365-X

Crossref Full Text | Google Scholar

22. Chen J, O'Malley R Jr. On the asymptotic solution of a two-parameter boundary value problem of chemical reactor theory. SIAM J Appl Math. (1974) 26:717–29. doi: 10.1137/0126064

Crossref Full Text | Google Scholar

23. Bigge J, Bohl E. Deformations of the bifurcation diagram due to discretization. Math Comput. (1985) 45:393–403. doi: 10.1090/S0025-5718-1985-0804931-X

PubMed Abstract | Crossref Full Text | Google Scholar

24. O'Malley RE Jr. Two-Parameter Singular Perturbation Problems. Stanford, CA: Stanford University (1966).

PubMed Abstract | Google Scholar

25. O'malley R. Two-parameter singular perturbation problems for second-order equations. J Math Mech. (1967) 16:1143–1164.

Google Scholar

26. Jha A, Kadalbajoo MK. A robust layer adapted difference method for singularly perturbed two-parameter parabolic problems. Int J Comput Math. (2015) 92:1204–21. doi: 10.1080/00207160.2014.928701

Crossref Full Text | Google Scholar

27. Das P, Mehrmann V. Numerical solution of singularly perturbed convection-diffusion-reaction problems with two small parameters. BIT Numer Math. (2016) 56:51–76. doi: 10.1007/s10543-015-0559-8

Crossref Full Text | Google Scholar

28. LinßT, Roos HG. Analysis of a finite-difference scheme for a singularly perturbed problem with two small parameters. J Math Anal Appl. (2004) 289:355–66. doi: 10.1016/j.jmaa.2003.08.017

PubMed Abstract | Crossref Full Text | Google Scholar

29. O'Malley RE. Introduction to Singular Perturbations. Academic Press (1974).

Google Scholar

30. Agbavon KM, Appadu AR, Khumalo M. On the numerical solution of Fisher's equation with coefficient of diffusion term much smaller than coefficient of reaction term. Adv Differ Equ. (2019) 2019:1–33. doi: 10.1186/s13662-019-2080-x

Crossref Full Text | Google Scholar

31. Appadu AR, Inan B, Tijani YO. Comparative study of some numerical methods for the Burgers-Huxley equation. Symmetry. (2019) 11:1333. doi: 10.3390/sym11111333

Crossref Full Text | Google Scholar

32. Agbavon KM, Appadu AR, Inan B, Tenkam HM. Convergence analysis and approximate optimal temporal step sizes for some finite difference methods discretising Fisher's equation. Front Appl Math Stat. (2022) 8:921170. doi: 10.3389/fams.2022.921170

Crossref Full Text | Google Scholar

33. Mohye MA, Munyakazi JB, Dinka TG. A nonstandard fitted operator finite difference method for two-parameter singularly perturbed time-delay parabolic problems. Front Appl Math Stat. (2023) 9:1222162. doi: 10.3389/fams.2023.1222162

Crossref Full Text | Google Scholar

34. Kehinde OO, Munyakazi JB, Appadu AR. A NSFD discretization of two-dimensional singularly perturbed semilinear convection-diffusion problems. Front Appl Math Stat. (2022) 8:861276. doi: 10.3389/fams.2022.861276

Crossref Full Text | Google Scholar

35. Alam MJ, Prasad HS, Ranjan R. An exponentially fitted integration scheme for a class of quasilinear singular perturbation problems. J Math Comput Sci. (2021) 11:3052–66. doi: 10.2891/jmcs/5589

Crossref Full Text | Google Scholar

36. Ranjan R, Prasad H. An efficient method of numerical integration for a class of singularly perturbed two point boundary value problems. WSEAS Trans Math. (2018) 17:265–73.

Google Scholar

37. Reddy Y, Reddy KA. Numerical integration method for general singularly perturbed two point boundary value problems. Appl Math Comput. (2002) 133:351–73. doi: 10.1016/S0096-3003(01)00246-6

Crossref Full Text | Google Scholar

38. Clavero C, Jorge J, Lisbona F. A uniformly convergent scheme for convection-diffusion parabolic problems. J Comput Appl Math. (2003) 154:415–29. doi: 10.1016/S0377-0427(02)00861-0

PubMed Abstract | Crossref Full Text | Google Scholar

39. Doolan EP, Miller JJ, Schilders WH. Uniform Numerical Methods for Problems with Initial and Boundary Layers. Dublin: Boole Press (1980).

PubMed Abstract | Google Scholar

40. Daba IT, Melesse WG, Kebede GD. Third-degree B-spline collocation method for singularly perturbed time delay parabolic problem with two parameters. Front Appl Math Stat. (2024) 9:1260651. doi: 10.3389/fams.2023.1260651

Crossref Full Text | Google Scholar