Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 15 July 2024
Sec. Numerical Analysis and Scientific Computation
This article is part of the Research Topic Advanced Numerical Methods for Solving Singularly Perturbed Differential Equations View all articles

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

  • 1Departement of Mathematics, Wollega University, Nekemte, Ethiopia
  • 2Departement of Mathematics, Jimma University, Jimma, Ethiopia

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 Δ¯=ΔΔ, 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

{𝕃ε,μwε2wx2(x,t)+μa(x,t)wx(x,t)-b(x,t)w(x,t)-wt(x,t)=                                -c(x,t)w(x,t-γ)+ϖ(x,t),(x,t)Δx,                    w|b=Υb(x,t),                    w|l=Υl(t),    w|r=Υr(t),    (1)

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)Δ¯ satisfying:

0<αa(x,t),0<βb(x,t),0<ςc(x,t),    (x,t)Δ¯.

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

w(0,0)=Υl(0),w(1,0)=Υr(0),w(0,-γ)=Υl(γ),w(1,-γ)=Υr(-γ),ε(Υb)xx(0,0)+μa(0,0)(Υb)x(0,0)-b(0,0)(Υb)(0,0)-(Υb)t(0,0)=-c(0,0)(Υb)(0,-γ)+ϖ(0,0),ε(Υb)xx(1,0)+μa(1,0)(Υb)x(1,0)-b(1,0)(Υb)(1,0)-(Υb)t(1,0)=-c(1,0)(Υb)(1,-γ)+ϖ(1,0),

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)C1+d/2,2+d(Δ¯) 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 [−γ, ] 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. [39]).

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. [1017]). 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, 1821]).

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(ε), 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, μ2/ε << 1 and μ2/ε >> 1, are contributed by the ratio of μ2 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, 2429].

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 [3032]. 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 Ξ(x,t)C2(Δ)C0(Δ¯), such that Ξ(x, t) ≥ 0, ∀(x, t) ∈ ∂Δ = ∧b ∪ ∧l ∪ ∧r & 𝕃ε,μ Ξ (x, t) ≤ 0 then, Ξ(x,t)0,(x,t)Δ¯.

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

||w(x,t)-Υb(x,0)||C1t,    (x,t)Δ¯,    (2)

where C1 is a constant which does not depend on ε.

Proof. As Υb(x, 0) is smooth on Δ¯, 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:

||w(x,t)||C*,for any    (x,t)Δ¯.    (3)

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

||w(x,t)||=||w(x,t)-Υb(x,0)+Υb(x,0)||,                    ||w(x,t)-Υb(x,0)||+||Υb(x,0)||,                    C1T+supx[0,1]||Υb(x,0)||C*.

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

||w(x,t)||Δ¯|ζ-1|||ϖ(x,t)||+max{||Υb(x,t)||,||Υl(t)||,||Υr(t)||},(x,t)Δ¯,    (4)

where 0 < ζ ≤ (b(x, t)+c(x, t)) and ||·||Δ¯ are used to denote maximum norm given in Equation (4) by ||w(x,t)||Δ¯=max(x,t)Δ¯|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:

||k+mwxktm||{Cε-k2 for μ2ε𝔑α   ,C(με)k(μ2ε)m for μ2ε𝔑α   ,    (5)

where 𝔑min(x,t)Δ¯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 ΛtM=[0,T], we use uniform mesh with time step k as follows:

ΛtM={tj=kj,j=0(1)M,tM=T,k=TM},
Λγm={tj=jk,j=0(1)m,tm=γ,k=γm},

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, tj+1/2) to derive the Crank–Nicolson in the temporal direction as follows:

wj+1(x)=wj+1/2(x)+k2dwj+1/2(x)dt+k28d2wj+1/2(x)dt2+k348d3wj+1/2(x)dt3+...    (6)
wj(x)=wj+1/2(x)-k2dwj+1/2(x)dt+k28d2wj+1/2(x)dt2-k348d3wj+1/2(x)dt3+...    (7)

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

dwj+1/2(x)dt=wj+1(x)-wj(x)k+TEj+1/2(x),    (8)

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

TEj+1/2(x)=k324d3wj+1/2(x)dt3+Higher order terms    (9)

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

ε2[wxxj+1+wxxj]+μ2[aj+1(x)wxj+1+aj(x)wxj]12[bj+1(x)wj+1(x)+bj(x)wj(x)]=wj+1(x)wj(x)k+TEj+1/2(x)+12[ϖj+1(x)+ϖj(x)]12[cj+1(x)wj+1m(x)+cj(x)wjm(x)],ε2[wxxj+1+wxxj]+μ2[aj+1(x)wxj+1+(a(x))j+1wj+1(x)+(a(x))jwj(x)+aj(x)wxj]=12[bj+1(x)wj+1(x)+bj(x)wj(x)]+wj+1(x)wj(x)k+12[ϖj+1(x)+ϖj(x)]12[cj+1(x)wj+1m(x)+cj(x)wjm(x)]+μ2[(a(x))j+1wj+1(x)+(a(x))jwj(x)],ε2[wxxj+1+wxxj]+μ2[aj+1(x)wj+1(x)+aj(x)wj(x)]=    12[bj+1(x)wj+1(x)+bj(x)wj(x)]+wj+1(x)wj(x)k+12[ϖj+1(x)+ϖj(x)]    12[cj+1(x)wj+1m(x)+cj(x)wjm(x)]+μ2[(a(x))j+1wj+1(x)+(a(x))jwj(x)],

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

{𝕃ε,μMwj+1(x):=ε2[wxxj+1+wxxj]+μ2[aj+1(x)wj+1(x)+aj(x)wj(x)]-    12[bj+1(x)wj+1(x)+bj(x)wj(x)]=wj+1(x)-wj(x)k+12[ϖj+1(x)+ϖj(x)]-    12[cj+1(x)wj+1-m(x)+cj(x)wj-m(x)]+μ2[(a(x))j+1wj+1(x)+(a(x))jwj(x)],wj+1(0)=Υl(tj+1),    wj+1(1)=Υr(tj+1),    0jM,wj+1(x)=Υb(x,tj+1),    (x,tj+1)[0,1]×[-γ,0],    (10)

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

{𝕃ε,μMwj+1(x):=ε2(wxx(x))j+1+μ2[(aj+1(x)wj+1(x))-(a(x))j+1wj+1(x)]                                                -(12bj+1(x)+1k)wj+1(x)=𝔊j+1(x),wj+1(0)=Υl(tj+1),    wj+1(1)=Υr(tj+1),    0jM,wj+1(x)=Υb(x,tj+1),    (x,tj+1)[0,1]×[-γ,0],    (11)

where 𝔊j+1(x)=-ε2(wxx(x))j+μ2[(a(x))jwj(x)-(aj(x)wj(x))]+(12bj(x)-1k)wj(x)+12[ϖj+1(x)+ϖj(x)]-12[cj+1(x)wj+1-m(x)+cj(x)wj-m(x)]. By using the initial condition, we can evaluate the right hand-side of Equation (11) as follows:

𝔊j+1(x)={-ε2(wxx(x))j+μ2[(a(x))jwj(x)-(aj(x)wj(x))]+(12bj(x)-1k)wj(x)-12[cj+1(x)Υb(x,tj+1-m)+cj(x)Υb(x,tj-m)]+12[ϖj+1(x)+ϖj(x)],j=0(1)m,-ε2(wxx(x))j+μ2[(a(x))jwj(x)-(aj(x)wj(x))]+(12bj(x)-1k)wj(x)-12[cj+1(x)w(x,tj+1-m)+cj(x)w(x,tj-m)]+12[ϖj+1(x)+ϖj(x)],j=m+1(1)M.

Lemma 3.1. (Semi-discrete Minimum Principle). At (j + 1)th time level, let Ψ(x,tj+1)C2(Δ¯x) be smooth function such that Ψ(0, tj+1) ≥ 0 and Ψ(1, tj+1) ≥ 0, then 𝕃ε,μMΨ(x,tj+1)0,xΔ implies Ψ(x,tj+1)0xΔ¯x.

Proof. Let (χ*,tj+1){(x,tj+1):xΔ¯x} be such that Ψ(χ*,tj+1)=minxΔ¯xΨ(x,tj+1) and let Ψ(χ*,tj+1)<0. Then, (χ*,tj+1){(0,tj+1),(1,tj+1)}, also Ψx(χ*,tj+1)=0 and Ψxx(χ*,tj+1)0. By using this assumption and Equation (11), we arrive on 𝕃ε,μMΨ(x,tj+1)>0, which is a contradiction to the hypothesis, 𝕃ε,μMΨ(x,tj+1)0. Hence, our asumption is wrong, so Ψ(χ*,tj+1)0, and this implies that Ψ(x,tj+1)0,xΔ¯x.

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

Ej+1=maxj+1T/k||wj+1(x)-Wj+1(x)||,    ej+1=|wj+1(x)-Wj+1(x)|,    (12)

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 wj+1(x) as the initial data rather than Wj+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)Δ¯=[0,1]×[0,T], independent of ε, μ, N & M, then the LTE associated with Equation (10) satisfies ||ej+1|| ≤ C(k3).

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

Lemma 3.3. The global error estimate Ej+1 defined in Equation (12) at the (j + 1)th time level is bounded by Ej+1C(k2), ∀jT/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:

ελ2+μaj+1(x)λ-βj+1(x)=0,    (13)

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

λ1=-μaj+1(x)2ε+(μaj+1(x)2ε)2+βj+1(x)ε0,λ2=-μaj+1(x)2ε-(μaj+1(x)2ε)2+βj+1(x)ε0.

These roots represent the boundary layer behavior of the solution in the neighborhood of x = 0 and x = 1 [4, 6, 14]. Let ϱ1=-maxx[0,1]λ2(x) and ϱ2=minx[0,1]λ1(x), we have three cases: (i) when μ2ε0, as ε0,ϱ1=ϱ2=βj+1ε, where 0 < β* < βj+1 and (ii) when εμ20, as μ0,ϱ1=ϱ2=μa*ε and ϱ2 = 0, where, 0 < a* < aj+1, (iii) when ϱ2 << ϱ1, the solutions exhibit a stronger boundary layer along x = 0 and x = 1 [4].

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

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

|dqwj+1(x)dxq|C(1+ϱ1qe-τϱ1x+ϱ2qe-τϱ2(1-x)),q=0(1)r,r>1    and    0<τ<1.    (14)

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: ΛxN={xi=ih,i=0(1)N,x0=0,xN=1,h=1/N}. Integrating Equation (10), in the interval [xi−1, xi+1], i = 1(1)N − 1, we obtain as follows:

ε2xi1xi+1[wxxj+1+wxxj](x)dx+μ2xi1xi+1[aj+1(x)wj+1(x)+aj(x)wj(x)]dx     12xi1xi+1[bj+1(x)wj+1(x)+bj(x)wj(x)]dx=1kxi1xi+1[wj+1(x)wj(x)]dx+     12xi1xi+1[ϖj+1(x)+ϖj(x)]dx12xi1xi+1[cj+1(x)wj+1m(x)+cj(x)wjm(x)]dx+μ2xi1xi+1[(a(x))j+1wj+1(x)+(a(x))jwj(x)]dx,ε2[(wx(xi+1))j+1(wx(xi1))j+1+(wx(xi+1))j(wx(xi1))j]+μ2[aj+1(xi+1)wj+1(xi+1)aj+1(xi1)wj+1(xi1)+aj(xi+1)wj(xi+1)aj(xi1)wj(xi1)]=xi1xi+1Θj(x)dx,    (15)

where

Θj(x)dx=12[bj+1(x)wj+1(x)+bj(x)wj(x)]+1k[wj+1(x)-wj(x)]+12[ϖj+1(x)+ϖj(x)]-12[cj+1(x)wj+1-m(x)+cj(x)wj-m(x)]+μ2[(a(x))j+1wj+1(x)+(a(x))jwj(x)].

In Equation (15), there is definite integral, and to approximate this integral, we use Simpson's 1/3rd formula, and by using the notation wj(xi)=Wij, Equation (15) is reduced to:

ε2[(Wx)i+1j+1-(Wx)i-1j+1+(Wx)i+1j-(Wx)i-1j]+μ2[ai+1j+1Wi+1j+1-ai-1j+1Wi-1j+1+ai+1jWi+1j-ai-1jWi-1j]                              =h3[Θi+1j+4Θij+Θi-1j].    (16)

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

(wx)ij+1Wi+1j+1-Wi-1j+12h,(wx)i+1j+13Wi+1j-4Wij+1+Wi-1j+12h,(wx)i-1j+1-Wi+1j+1+4Wij+1-3Wi-1j+12h,    (17)

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:

εh[Wi+1j+1-2Wij+1+Wi-1j+1+Wi+1j-2Wij+Wi-1j]+μ2[ai+1j+1Wi+1j+1-ai-1j+1Wi-1j+1+ai+1jWi+1j-ai-1jWi-1j]                             =h3[Θi+1j+4Θij+Θi-1j].    (18)

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

εσ(ρ)h[Wi+1j+1-2Wij+1+Wi-1j+1+Wi+1j-2Wij+Wi-1j]+μ2[ai+1j+1Wi+1j+1-ai-1j+1Wi-1j+1+ai+1jWi+1j-ai-1jWi-1j]                           =h3[Θi+1j+4Θij+Θi-1j].    (19)

Let ρ=hε and taking the limit of Equation (19) as h → 0

limh0σ(ρ)ρ[Wi+1j+1-2Wij+1+Wi-1j+1+Wi+1j-2Wij+Wi-1j]+limh0μ2[ai+1j+1Wi+1j+1-ai-1j+1Wi-1j+1+ai+1jWi+1j-ai-1jWi-1j]=0.    (20)

3.3 Determination of fitting factor

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

Wj+1(x)={W0j+1(x)+[W0j+1(0)-Υlj+1(0)]exp(-μaj+1(0)ε(x)),at Left Layer,W0j+1(x)+[W0j+1(1)-Υrj+1(1)]exp(μaj+1(1)ε(1-x)),at Right Layer,

where W0j+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 = xi = h gives

Wij+1={(W0)ij+1+[W0j+1(0)-Υlj+1(0)]exp(-μaj+1(0)ε(ih)),left layer,(W0)ij+1+[W0j+1(1)-Υrj+1(1)]exp(μaj+1(1)ε(1-ih)),right layer.    (21)

For left layer, substituting Equation (21) into Equation (20) and taking ρ=hε

σρ[e(-μaj+1(0)ρ)-2+e(μaj+1(0)ρ)]=-μ2aj+1(0)[e(-μaj+1(0)ρ)+e(μaj+1(0)ρ)],-σρ[e(μaj+1(0)ρ2)-e(-μaj+1(0)ρ2)]2=-μ2aj+1(0)[e(μaj+1(0)ρ2)-e(-μaj+1(0)ρ2)],

After some simplification of the above expression, we have

σ(ρ)=ρμaj+1(0)2coth(μaj+1(0)ρ2).    (22)

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:

σ(ρ)={ρμaj+1(0)2coth(μaj+1(0)ρ2), for x=0   ,ρμaj+1(1)2coth(μaj+1(1)ρ2), for x=1   .    (23)

where ρ = h/ε. Substituting Equation (22) into Equation (19) and solving for Θij, we evaluate at:

{[σ(ρ)ρ-μ2ai-1j+1-h6bi-1j+1-h3k-h6μ(a)i-1j+1]Wi-1j+1+[-2σ(ρ)ρ-2h3bij+1+4h3k-2h3μ(a)ij+1]Wij+1+[σ(ρ)ρ+μ2ai+1j+1-h6bi+1j+1-h3k-h6μ(a)i+1j+1]Wi+1j+1=[-σ(ρ)ρ+μ2ai-1j+h6bi-1j-h3k+h6μ(a)i-1j]Wi-1j+[2σ(ρ)ρ+2h3bij-4h3k+2h3μ(a)ij]Wij+[-σ(ρ)ρ-μ2ai+1j+h6bi+1j-h3k+h6μ(a)i+1j]Wi+1j+h6[ϖi+1j+1+ϖi+1j]+2h3[ϖij+1+ϖij]+h6[ϖi-1j+1+ϖi-1j]-h6[ci-1j+1Wi-1j+1-m+ci-1jWi-1j-m]-2h3[cij+1Wij+1-m+cijWij-m]                 -h6[ci+1j+1Wi+1j+1-m+ci-1jWi+1j-m].    (24)

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

{𝕃ε,μN,MWij+1=Ai-Wi-1j+1+AicWij+1+Ai+Wi+1j+1=Bi-Wi-1j+BicWij+Bi+Wi+1j+Fjj,                     Wij+1=Υb(xi,tj+1),i=1,...,N-1,-(m+1)j-1,                     W0j+1=Υl(tj+1),    WN+1j+1=Υr(tj+1),0jM-1.    (25)

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

{Ai-=σ(ρ)ρ-μ2ai-1j+1-h6bi-1j+1-h3k-h6μ(a)i-1j+1,Aic=-2σ(ρ)ρ-2h3bij+1-4h3k-2h3μ(a)ij+1,Ai+=σ(ρ)ρ+μ2ai+1j+1-h6bi+1j+1-h3k-h6μ(a)i+1j+1,Bi-=-σ(ρ)ρ+μ2ai-1j+h6bi-1j-h3k+h6μ(a)i-1j,Bic=2σ(ρ)ρ+2h3bij-4h3k+2h3μ(a)ij,Bi+=-σ(ρ)ρ-μ2ai+1j+h6bi+1j-h3k+h6μ(a)i+1j,    (26)

and the right hand side is given as follows:

Fij={h6[ci1j+1Υb(xi1,tj+1)+2cij+1Υb(xi,tj+1)+ci+1j+1Υb(xi+1,tj+1)]h6[ci1jΥb(xi1,tj)+2cijΥb(xi,tj)+ci+1jΥb(xi+1,tj)]+h6[ϖi1j+1+2ϖij+1+ϖi+1j+1]+h6[ϖi1j+2ϖij+ϖi+1j],for tj<m,j=0(1)M1,h6[ci1j+1W(xi1,tj+1m)+2cij+1W(xi,tj+1m)+ci+1j+1W(xi+1,tj+1m)]h6[ci1jW(xi1,tjm)+2cijW(xi,tjm)+ci+1jW(xi+1,tjm)]+h6[ϖi1j+1+2ϖij+1+ϖi+1j+1]+h6[ϖi1j+2ϖij+ϖi+1j],for tj>m,j=0(1)M1.

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

|Aic||Ai-|+|Ai+|.    (27)

This shows the coefficient matrix associated with 𝕃ε,μ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, tj+1) ≥ 0, Ξ(1, tj+1) ≥ 0, and 𝕃,μN,MΞ(xi,tj+1)0 on ΔxN×ΔtM. Then, Ξ(xi, tj) ≥ 0 at each point Δ¯xN×Δ¯tM, where ΔxN and ΔtM are discretized domain.

Proof. Let(xϑ, tj+1) ∈ {(xi, tj+1):i = 1(1)N} such that Ξj+1(xϑ)=min1iNΞj+1(xi) and let Ξj+1(xϑ)<0. From this (xϑ, tj+1)∉{(0, tj+1), (1, tj+1)}, which implies that (xϑ, tj+1) ∈ {(xi, tj+1):i = 1(1)N − 1}. Then, from Equation (25), we have

𝕃ε,μN,MΞϑj+1=Ai-Ξϑ-1j+1+AicΞϑj+1+Ai+Ξϑ+1j+10,

since |Ai-|>0,|Aic|>0 and |Ai+|>0, which is a contradiction. As a result, Ξj+1(xϑ)0. Therefore, Ξj+1(xi)0,(xi,tj+1){(xi,tj+1):i=1(1)N-1}.

According to Lemma (4.1), the discrete operator 𝕃ε,μ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 Wij+1 of the discrete scheme in Equation (25) on 𝕃ε,μN,M gratifies the following estimate:

||Wij+1||ζ-1||𝕃ε,μN,MWij+1||+max{|Υl(tj+1)|,|Υr(tj+1)|,|Υb(xi,tj+1)|},    (28)

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

Proof. Let Γ=ζ-1||𝕃ε,μN,MWij+1||+max{|Υl(tj+1)|,|Υr(tj+1)|,|Υb(xi,tj+1)|} and describe the barrier functions as follows:

(±)ij+1=Γ±Wij+1.    (29)

The discrete function ±(xi,tj+1) given in Equation (29) on the boundary and interval functions is as follows:

±(0,tj+1)=Γ+max{|Υl(tj+1)|,|Υr(tj+1)|,|Υb(0,tj+1)|}±Υl(tj+1)0,±(1,tj+1)=Γ+max{|Υl(tj+1)|,|Υr(tj+1)|,|Υb(1,tj+1)|}±Υr(tj+1)0,±(xi,0)=Γ+max{|Υl(0)|,|Υr(0)|,|Υb(xi,0)|}±Υb(xi,0)0.

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

𝕃ε,μN,M(±)ij+1=Γ±Wij+1=Ai-(Γ±Wi-1j+1)+                                                                Aic(Γ±Wij+1)+Ai+(Γ±Wi+1j+1),                                                                =(Ai-+Aic+Ai+)Γ±Hij0,

where Hij 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:

Ti(h)=Ai-Wi-1j+1+AicWij+1+Ai+Wi+1j+1-(Bi-Wi-1j+BicWij+Bi+Wi+1j+F)    (30)

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

Ti(h)={[Ai-Wi-1j+1-h(εWxxj+1+μaj+1(x)Wxj+1-βj+1(x)Wj+1)i-1]+[AicWij+1-h(εWxxj+1+μaj+1(x)Wxj+1-βj+1(x)Wj+1)i]+[Ai+Wi+1j+1-h(εWxxj+1+μaj+1(x)Wxj+1-βj+1(x)Wj+1)i+1]+[Bi-Wi-1j-h(εWxxj+μaj(x)Wxj-βj(x)Wj)i-1]+[BicWij-h(εWxxj+μaj(x)Wxj-βj(x)Wj)i]+[Bi+Wi+1j-h(εWxxj+μaj(x)Wxj-βj(x)Wj)i+1]    (31)

Replacing Taylor's series expansion of

{Wi-1j+1Wij+1-h1!dWij+1dx+h22!d2Wij+1dx2-h33!d3Wij+1dx3+h44!d4Wij+1dx4+O(h5),Wi+1j+1Wij+1+h1!dWij+1dx+h22!d2Wij+1dx2+h33!d3Wij+1dx3+h44!d4Wij+1dx4+O(h5),Wi-1jWij-h1!dWijdx+h22!d2Wijdx2-h33!d3Wijdx3+h44!d4Wijdx4+O(h5),Wi+1jWij+h1!dWijdx+h22!d2Wijdx2+h33!d3Wijdx3+h44!d4Wijdx4+O(h5),(Wx)i-1j+1(Wx)ij+1-h1!(Wxx)ij+1+h22!(Wxxx)ij+1-h33!(Wxxxx)ij+1+O(h4),(Wx)i+1j+1(Wx)ij+1+h1!(Wxx)ij+1+h22!(Wxxx)ij+1+h33!(Wxxxx)ij+1+O(h4),(Wx)i-1j(Wx)ij-h1!(Wxx)ij+h22!(Wxxx)ij-h33!(Wxxxx)ij+O(h4),(Wx)i+1j(Wx)ij+h1!(Wxx)ij+h22!(Wxxx)ij+h33!(Wxxxx)ij+O(h4),(Wxx)i-1j+1(Wxx)ij+1-h1!(Wxxx)ij+1+h22!(Wxxxx)ij+1+O(h3),(Wxx)i+1j+1(Wxx)ij+1+h1!(Wxxx)ij+1+h22!(Wxxxx)ij+1+O(h3),(Wxx)i+1j(Wxx)ij+h1!(Wxxx)ij+h22!(Wxxxx)ij+O(h3),(Wxx)i-1j(Wxx)ij+h1!(Wxxx)ij+h22!(Wxxxx)ij+O(h3),    (32)

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

Ti(h)=[ξ0j+1Wij+1+ξ0jWij]+            [ξ1j+1Wx,ij+1+ξ1jWx,ij]+            [ξ2j+1Wxx,ij+1+ξ2jWxx,ij]+            [ξ3j+1Wxxx,ij+1+ξ3jWxxx,ij]+higher order terms,    (33)

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

{ξ0j+1=(Ai-+Aic+Ai+)+h(βi-1j+1+βij+1+βi+1j+1),ξ0j=-(Bi-+Bic+Bi+)+h(βi-1j+βij+βi+1j),ξ1j+1=h(Ai+-Ai-)-hμ(ai-1j+1+aij+1+ai+1j+1)-h2(βi-1j+1-βi+1j+1),ξ1j=h(Bi--Bi+)-hμ(ai-1j+aij+ai+1j)+h2(βi+1j-βi-1j),ξ2j+1=h22(Ai-+Ai+)-3hε+h2μ(ai-1j+1-ai+1j+1)+h32(βi-1j+1+βi+1j+1),ξ2j=-h22(Ai++Ai-)-3hε+h2μ(ai-1j-ai+1j)+h32(βi-1j+βi+1j),ξ3j+1=h36(Ai+-Ai-)-h32μ(ai-1j+1+ai+1j+1)-h46(βi-1j+1-βi+1j+1),ξ3j=h36(Bi--Bi+)-h32μ(ai-1j+ai+1j)-h46(βi+1j-βi-1j).

Substituting and simplifying the values of Ai-,Aic,Ai+,Bi-,Bic,Bi+ and restricting the expansion of the coefficients to the first term gives ξ0j+1=ξ0j=ξ1j+1=ξ1j=0 and

{ξ2j+1=h2σ(ρ)ρ-3hε-h36(11(a)ij+1+5βij+1)+O(h3),ξ2j=h2σ(ρ)ρ-3hε-h36(11(a)ij+5βji)+O(h3).    (34)

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

|Ti(h)|ξ2j+1(Wxx)ij+1+ξ2j(Wxx)ijhε(σ-3)(Wxx)ij+1,    since    ρ=hε.    (35)

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

coth(ψ)=1ψ+ψ3-ψ345+O(ψ5)ψcoth(ψ)=1+ψ23+O(ψ4),

where ψ=ρμaj+1(xi)2. Now |σ-3|=|1+(ρμaj+1(xi)2)23-3|Cρ2μ2. From this, Equations (34, 35) are reduced to

|Ti(h)|hε(Cρ2μ2)(Wxx)ij+1=Ch3μ2ε(Wxx)ij+1.    (36)

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

|Ti(h)|{    Ch3με when μ2ε   ,Ch3μ2ε when εμ2   .    (37)

Lemma 4.3. Let the solutions to Equations (1, 25) be w(xi, tj+1) and 𝕃ε,μN,MWij+1, respectively. Next, in the spatial direction, the proposed scheme satisfies the estimation [6]:

|𝕃ε,μN,MWij+1-w(xi,tj+1)|{Ch2 when μ2ε   ,  Ch when εμ2   .    (38)

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

Θ±(xi,tj+1)=Ti(h)±(𝕃ε,μN,MWij+1-w(xi,tj+1)).

Now from boundary conditions, Θ±(x0,tj+1)0,Θ±(xN,tj+1)0 and 𝕃ε,μN,MΘ±(xi,tj+1)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(xi, tj+1) and Wij+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:

max0iN,0jM||w(xi,tj+1)-Wij+1||{C(h2+k2) when μ2ε   ,    C(h+k2) when εμ2   ,    (39)

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 Wj+1(xi) 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 μ2/ε ≤ 1 and shows convergence, which is independent of the perturbation parameters. Additionally, the method shows first-order convergence for ε/μ2 << 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ε,μN,M and rate of convergence R¯ε,μN,M of the scheme using the double mesh principle given in the study by Doolan et al. [39] as follows:

Eε,μN,M=max0iN,0jM|Wi,jN,M-Wi,j2N,2M|,R¯ε,μN,M=log(Eε,μN,M)-log(Eε,μ2N,2M)log(2),

where Wi,jN,M and Wi,j2N,2M are numerical solutions computed on the mesh N × M and 2N × 2M, respectively. Additionally, the uniform rate of convergence R¯N,M and uniform maximum point-wise error EN, M are given in the study by Duressa and Mekonnen [17] by the following formula:

EN,M=maxε,μEε,μN,M,    R¯N,M=log(EN,M)-log(E2N,2M)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]

{ε2wx2+μ(1+x)wx-w(x,t)-wt=    -w(x,t-γ)+16x2(1-x)2,(x,t)(0,1)×(0,2],w(0,t)=0,    w(1,t)=0,    t(0,2],w(x,t)=0,    (x,t)[0,1]×[-γ,0].

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

{ε2wx2+μ(1+x-x2+t2)wx-(1+5xt)w(x,t)-wt=    -w(x,t-γ)+(x-x2)(et-1),    (x,t)(0,1)×(0,2],w(0,t)=0,    w(1,t)=0,    t(0,2],w(x,t)=0,    (x,t)[0,1]×[-γ,0].

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

{ε2wx2+μ(1+x2)wx-(1+5xt)w(x,t)-wt=-w(x,t-γ)+16x2(1-x)2,w(0,t)=0,    w(1,t)=0,    t(0,2],w(x,t)=0,    (x,t)[0,1]×[-γ,0].

The computed maximum point-wise error Eε,μN,M and rate of convergence R¯ε,μN,M for examples (5.1–5.3) are presented in Tables 18 with different values of ε, μ, N, and M. Tables 15, 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 ε/μ2 << 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 13, 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 46 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
www.frontiersin.org

Table 1. Computed Eε,μN,M,EN,M,andR¯ε,μN,M for example (5.1) at μ = 10−9 and comparison with Negero [4] and Sumit et al. [18].

Table 2
www.frontiersin.org

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

Table 3
www.frontiersin.org

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

Table 4
www.frontiersin.org

Table 4. Computed Eε,μN,M,EN,M, and R¯N,M for example (5.1) at ε = 2−10 and comparison with the study by Singh et al. [7] for varying values of μ.

Table 5
www.frontiersin.org

Table 5. Computed Eε,μN,M,R¯ε,μN,M,EN,M,R¯N,M and CPU time (in seconds) for example (5.2) at μ = 10−9 and comparison with the study by Negero [4] and Sumit et al. [18].

Table 6
www.frontiersin.org

Table 6. Computed Eε,μN,M,EN,M, and R¯N,M for example (5.2) at μ = 10−9 after extrapolation.

Table 7
www.frontiersin.org

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

Table 8
www.frontiersin.org

Table 8. Computed Eε,μN,M,EN,M, and R¯N,M for example (5.3) at ε = 10−10 after extrapolation.

Figure 1
www.frontiersin.org

Figure 1. The numerical solution profile: (A) for N = M = 64, ε = 10−8, μ = 1, and (B) for N = M = 128, ε = 10−12, and μ = 10−10 for example (5.1).

Figure 2
www.frontiersin.org

Figure 2. The numerical solution profile: (A) for N = M = 64, ε = 10−8, and μ = 1 and (B) for N = 256, M = 128, ε = 10−6, and μ = 10−10 for example (5.2).

Figure 3
www.frontiersin.org

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

Figure 4
www.frontiersin.org

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
www.frontiersin.org

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
www.frontiersin.org

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

Keywords: singularly perturbed problems, delay parabolic differential equation, numerical integration, Simpson's 1/3rd rule, parameter-uniform

Citation: Cheru SL, Duressa GF and Mekonnen TB (2024) Numerical integration method for two-parameter singularly perturbed time delay parabolic problem. Front. Appl. Math. Stat. 10:1414899. doi: 10.3389/fams.2024.1414899

Received: 09 April 2024; Accepted: 28 June 2024;
Published: 15 July 2024.

Edited by:

Eun-Jae Park, Yonsei University, Republic of Korea

Reviewed by:

Sunil Kumar, Indian Institute of Technology (BHU), India
Appanah Rao Appadu, Nelson Mandela University, South Africa

Copyright © 2024 Cheru, Duressa and Mekonnen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Tariku Birabasa Mekonnen, seenaa29@gmail.com

Disclaimer: 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.