Front. Appl. Math. Stat. Frontiers in Applied Mathematics and Statistics Front. Appl. Math. Stat. 2297-4687 Frontiers Media S.A. 10.3389/fams.2023.1260651 Applied Mathematics and Statistics Original Research Third-degree B-spline collocation method for singularly perturbed time delay parabolic problem with two parameters Daba Imiru Takele * Melesse Wondwosen Gebeyaw Kebede Guta Demisu Department of Mathematics, Dilla University, Dilla, Ethiopia

Edited by: Justin B. Munyakazi, University of the Western Cape, South Africa

Reviewed by: Youssri Hassan Youssri, Cairo University, Egypt

Ahmed Atta, Ain Shams University, Egypt

 *Correspondence: Imiru Takele Daba imirutakele@gmail.com 08 01 2024 2023 9 1260651 18 07 2023 13 12 2023 Copyright © 2024 Daba, Melesse and Kebede. 2024 Daba, Melesse and Kebede

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.

This study deals with a fitted third-degree B-spline collocation method for two parametric singularly perturbed parabolic problems with a time lag. The proposed method comprises the Cranck-Nicolson method for time discretization and the third-degree B-spline method spatial variable discretization. Rigorous numerical experimentations were carried out on some test examples. The obtained numerical results depict that the proposed scheme is more accurate than some methods existing in the literature. Parameter convergence analysis of the scheme is carried out and shows the present scheme is (ε−μ)−uniform convergent with the order of convergence ((Δt)2 + ℓ2).

 third degree B-spline collocation method Crank-Nicholson method time delay two-parametric parameter uniform numerical method section-at-acceptance Numerical Analysis and Scientific Computation 1 Introduction

We consider the following two-parametric singularly perturbed parabolic problem with the time lag of the form:

{ε∂2u(s,t)∂s2+γa(s,t)∂u(s,t)∂s-b(s,t)u(s,t)-∂u(s,t)∂t=-c(s,t)u(s,t-τ)+f(s,t),(s,t)∈Ω,u(s,t)=θ(s,t),(s,t)∈[0,1]×[-τ,0],u(0,t)=q0(t),0≤t≤T,u(1,t)=q1(t),0≤t≤T,

where Ω=ΩsN×ΩtM=(0,1)×(0,T], T is fixed time, ε, γ(0 < ε, γ ≪ 1) are two small perturbation parameters and τ > 0 is delay parameter. For the existence and uniqueness of the solution, the functions a(s, t), b(s, t), c(s, t), f(s, t) q0(t), q1(t) and θ(s, t) are sufficiently smooth and bounded with a(s, t) ≥ α > 0, b(s, t) ≥ ϖ > 0, b(s, t)+c(s, t) ≥ β > 0. The mathematical models related to two-parameter singularly perturbed problem of type (1) arises in transport phenomena in chemistry, biology, chemical reactor theory [1], lubrication theory [2], and dc motor theory [3] and flow through unsaturated porous media [4].

Due to the dual presence of singular parameters and time delay in the problem finding oscillation-free solution of Equation (1) using contemporary computational methods are very tough task. To overcome such limitations, some scholars developed parameter uniform numerical methods. For instance, Sumit et al. [5] presented a robust numerical scheme using a hybrid monotone finite difference scheme on uniform mesh in time and a piecewise uniform Shishkin mesh in space for Equation (1). Govindarao et al. [6] developed a uniformly convergent computational method using the implicit Euler scheme for temporal discretization on a uniform mesh and the upwind difference scheme for the spatial discretization on the Shishkin type meshes for Equation (1). The authors in the study mentioned in the references [7] and [8] constructed a parameter uniform numerical method based on a uniform mesh for Equation (1). The schemes in the study mentioned in the references [5] and [6] need apriori knowledge about the location and the width of the boundary layer(s) which might be difficult to understand for beginner researchers. Exponentially fitted difference (EFD) schemes have gained popularity as a powerful technique to singularly perturbed problems. For instance, the authors in the study mentioned in the references [9–11] suggested different EFD schemes for singularly perturbed two-point boundary-value problems. The aforementioned studies above motivate us, to propose and analyze the fitted third-degree B-spline collocation scheme for Equation (1). Analytical aspects of the problem, the description of the proposed scheme, and its corresponding analysis are preserved in the subsequent sections.

1.1 Apriori estimates for the solution of the continuous problem

When τ < ε, the use of Taylor's series expansion for the term containing shift argument is valid [12]. On applying Taylor's series expansion on u(s, t − τ), we have:

u(s,t-τ)=u(s,t)-τ∂u(s,t)∂s+O(τ2),

Now taking Equation (2) into Equation (1), we have:

{Lε,γu(s,t)≡ε∂2u(s,t)∂s2+μa(s,t)∂u(s,t)∂s-η(s,t)u(s,t)-dτ(s,t)∂u(s,t)∂t=f(s,t),u(s,0)=θ0(s),s∈Ωsu(0,t)=q0(t),t∈Ωt,u(1,t)=q1(t),t∈Ωt.

where, η(s, t) = b(s, t)−c(s, t), dτ(s, t) = 1 + τc(s, t).

Lemma 1. (Continuous minimum principle). If u∈C(2,1)(Ω̄), u|∂Ω ≥ 0 and Lε,γu|Ω≤0. Then, u|Ω¯≥0.

Proof. Assume that the arbitrary function u attains its minimum value at the point (s*,t*)∈Ω̄ such that u(s*,t*)=min(s,t)∈Ω̄u(s,t) and suppose that u(s*, t*) < 0. Clearly, (s*, t*)∉∂Ω. Therefore, ∂u∂x(s*,t*)=0, ∂u∂t(s*,t*)=0 and ∂2u∂s2(s*,t*)≥0. Moreover, for (s*, t*) ∈ Ω, we have

Lε,γu(s*,t*)=ε∂2u∂s2(s*,t*)+γa(s*,t*)∂u∂s(s*,t*)-η(s*,t*)u(s*,t*)-dτ(s,t)∂u∂t(s*,t*)≥0,

which is illogicality to the supposition that Lε,μu|Ω≤0. It follows that u|Ω¯≥0.

Lemma 2. (Uniform stability estimate) Let u(s, t) be the solution of Equation (3), then we have:

∥u∥Ω̄≤max{|θ0(s,t)|,|q0(t)|,|q1(t)|}+β-1∥f∥,

where ‖.||Ω¯ is used to denote maximum norm given by ‖u‖Ω¯=max(s,t)∈Ω¯|u(s,t)|.

Proof. Let Θ± be the barrier functions given as

Θ±(s,t)=max{|θ0(s)|,|q0(t)|,|q1(t)|}+β-1∥f∥±u(s,t).

Then at the initial and boundary conditions, we have

Θ±(s,0)=max{|θ0(s)|,|q0(0)|,|q1(0)|}+β-1∥f∥±θ0(s),≥0. Θ±(0,t)=max{|θ0(0)|,|q0(t)|,|q1(t)|}+β-1∥f∥±u(0,t),=max{|θ0(0)|,|q0(t)|,|q1(t)|}+β-1∥f∥±q0(t),≥0. Θ±(1,t)=max{|θ0(1)|,|q0(t)|,|q1(t)|}+β-1∥f∥±u(1,t),=max{|θ0(1)|,|q0(t)|,|q1(t)|}+β-1∥f∥±q1(t),≥0.

Applying the differential operator in (3) on Θ±(s, t), we have

ℒε,γΘ±(s,t)             =ε∂2Θ±(s,t)∂s2+γ(s,t)∂Θ±(s,t)∂s                                       −η(s,t)Θ±(s,t)−dτ(s,t)∂Θ±(s,t)∂t,              =ε∂2∂s2[max{|θ0(s)|,|q0(t)|,|q1(t)|}                                          +β−1‖f‖±u(s,t)]              +γa(s,t)∂∂s[max{|θ0(s)|,|q0(t)|,|q1(t)|}                                          +β−1‖f‖±u(s,t)]              −η(s,t)[max{|θ0(s)|,|q0(t)|,|q1(t)|}                                          +β−1‖f‖±u(s,t)]              −dτ(s,t)∂∂t[max{|θ0(s)|,|q0(t)|,|q1(t)|}                                          +β−1‖f‖±u(s,t)],              =−η(s,t)[max{|θ0(s)|,|q0(t)|,|q1(t)|}                                          +β−1‖f‖±u(s,t)]±ℒε,μu(s,t),              =−η(s,t)[max{|θ0(s,t)|,|q0(t)|,|q1(t)|}                                          +β−1‖f‖±u(s,t)]±f(s,t)

Since η(s, t) ≥ β > 0 we have η(s, t)β−1 > 1, and the fact ||f|| ≥ f(s, t) implies that Lε,μΘ±(s,t)≤0. Hence, Lemma (1) confirms Θ±(s,t)≥0,∀(s,t)∈Ω̄, which yields the desired result.

The above Lemma (1) and Lemma (2) are guarantees for the existence and uniqueness solution of Equation (3), respectively.

 2 Formulation of the numerical method 2.1 Temporal semi-discretization

We describe the uniform mesh for the domain Ωt as

ΩtM={tj=jΔt,j=0,⋯,M,Δt=T/M}.

Applying the Crank-Nicholson method on t−direction of Equation (3) yields

{εUss(s,tj+1)+γa(s,tj+1)Us(s,tj+1)-p(s,tj+1)U(s,tj+1)=g(s,tj+1),U(s,tj+1)=θ0(s),s∈Ω¯s,U(0,tj+1)=q0(tj+1),0≤tj+1≤T,U(1,tj+1)=q1(tj+1),0≤tj+1≤T,

where

g(s,tj+1)=-εUss(s,tj)-γa(s,tj)Us(s,tj)+Q(s,tj)U(s,tj)+f(s,tj+1)+f(s,tj),

p(s,tj+1)=η(s,tj+1)+2dτ(s,tj+1)Δt and Q(s,tj)=η(s,tj)+2dτ(s,tj)Δt.

The local truncation error (LTE) of the temporal semi-discretization of Equation (6) is defined as ej+1 = u(s, tj+1) − U(s, tj+1), where U(s, tj+1) is the solution of the following BVP

{εUss(s,tj+1)+γa(s,tj+1)Us(s,tj+1)-p(s,tj+1)U(s,tj+1)=g(s,tj+1),U(0,tj+1)=q0(tj+1), 0≤tj+1≤T,U(1,tj+1)=q1(tj+1), 0≤tj+1≤T.

Now, we state the bounds for the errors in the local and global as follows.

Lemma 3. (LTE) If

|∂kU(s,t)∂sk|≤C,(s,t)∈Ω̄,  0≤k≤2,

the LTE of Equation (6) is given by

∥ej+1∥∞≤C(Δt)3,1≤j≤M.

Proof. See [6].

Lemma 4. Underneath Lemma (3), the global error estimate (GEE)at tj is given by

∥Ej∥∞≤C(Δt)2,   j≤T/Δt.

Proof. From Lemma 3, we have

‖Ej||∞=‖∑i=1jei||∞,j≤T/Δt≤‖e1||∞+‖e2||∞+‖e3||∞+⋯+‖ej||∞                      ≤Cj(Δt)3              (by  Lemma 3)                          ≤C(jΔt)((Δt)2)                                                       ≤CT((Δt)2)         (jΔt≤T)                                        ≤C((Δt)2),                                                              

where C(C > 0) is constant independent of ε, μ, and Δt.

 2.2 Spatial discretization

To approximate Equation (6), we employ the third-degree B-spline collocation method. We divide the space domain Ω¯s as Ω¯s={0=s0<s1<⋯<sN=1} with si = iℓ, where ℓ = 1/N. The third-degree B-Spline (Bi(s)) can be defined as [13]

Bi(s)=1ℓ3{(s-si-2)3,                                        s∈[si-2,si-1],ℓ3+3ℓ2(s-si-1)+3ℓ(s-si-1)2-3(s-si-1)3,s∈[si-1,si],ℓ3+3ℓ2(si+1-s)+3ℓ(si+1-s)2-3(si+1-s)3,s∈[si,si+1],(si+2-s)3,                                            s∈[si+1,si+2],0,                                            otherwise,

where {Bi}i=-1N+1 form a basis over Ω¯s. Let S3(Ω¯s) be the set of all third-degree B-spline functions over the subintervals Ω¯s. Let λ = {B−1, B0, B1, ⋯   , BN+1} and ξ3(Ω¯s) be the set of all linear components of Bi′s. The function Bi′s are linearly independent, thus ξ3(Ω¯s) is dimensional space of S3(Ω¯s) [14].

In the third-degree B-spline collocation method, we approximate the exact solution U(s, tj+1) by S(s) in the form:

S(s)≈∑i=-1N+1αi(t)Bi(s),

where αi(t)'s are time-dependent parameters to be determined from the collocation method together with using the boundary and initial conditions. The values of Bi(s) and its derivatives are presented in Table 1.

Coefficients of third-degree B-splines and their derivatives at knots.

s si−2 si−1 si si+1 si+2
Bi(s) 0 1 4 1 0
Bi′(s) 0 −3/ℓ 0 3/ℓ 0
Bi ″(s) 0 6/ℓ2 −12/ℓ2 6/ℓ2 0

From Equation (10) and Table 1, we have

{(S)ij=αi-1j+4αij+αi+1j,(Ss)ij=3ℓ(αi+1j-αi-1j),(Sss)ij=6ℓ2(αi-1j-2αij+αi+1j).

At s = si undertaking the notation U(si, tj+1) = Û(si) Equation (6) and introducing a fitting factor σ(si) in the resulting equation to handle the influence of perturbation parameter on solution profile, we have

{Lε,γN,MÛ(si)=σ(si)εÛss(si)+γa(s)Ûs(si)-p(si)Û(si)=g(si),Û0(si)=θ0(si),0≤i≤N,Û(0)=qj+1(0),  0≤j<M,Û(1)=qj+1(1),  0≤j<M.

In the corresponding time level, plugging Equation (11) at the knots into Equation (12) and simplifying, we obtain

Yi-αi-1j+1+Yicαij+1+Yi+αi+1j+1=Di-αi-1j+Dicαi+Di+αi+1j+ℓ2(fij+fij+1),i=0,1,2,...,N,

where

Yi-=6σiε-3ℓγai-ℓ2pi,Yic=-12σiε-4ℓ2pi,Yi+=6σiε+3ℓγai-ℓ2pi. Di-=-6σiε+3ℓγai+ℓ2pi,Dic=12σiε+4ℓ2pi,Di+=-6σiε-3ℓγai+ℓ2pi.

The fitting factor σi is given as [15, 16]

σi=γaiρ2coth(γaiρ2),

where ρ = ℓ/ε.

At each (j + 1)th time level, this gives (N + 1) equations with (N + 3) unknowns [α-1j+1,α0j+1,α1j+1,⋯,αN+1j+1]. Using boundary conditions in Equation (11) and the first and last equation, of Equation (13), we get

Yαij+1=Dαij+H,

where

Y=((Y0c-4Y0-)(Y0+-Y0-)Y1-Y1cY1+⋯⋯⋯⋯⋯⋯YN-1-YN-1cYN-1+(YN--YN+)(YNc-4YN+)), D=((D0c-4D0-)(D0+-D0-)D1-D1cD1+⋯⋯⋯⋯⋯⋯DN-1-DN-1cDN-1+(DN--DN+)(DNc-4DN+)), H=(ℓ2(f0j+f0f+1)-Y0-q0(tj+1)+D0-q0(tj)ℓ2(f1j+f1f+1)⋮ℓ2(fN-1j+fN-1f+1)ℓ2(fNj+fNf+1)-YN+q1(tj+1)+DN+q1(tj)), αij+1=(α0j+1,αij+1,α1j+1,⋯,αNj+1)t,  andαij=(α0j,αij,α1j,⋯,αNj)t,

which is a diagonally dominant and non-singular matrix. Consequently, we get a exclusively solvable system of equations. So, we can solve the system in the study mentioned in the reference [15] for α′s, substituting these values into Equation (10), and we obtain the required approximate solution.

 3 Convergence analysis

Lemma 5. [15] The B-Splines B−1, B0, B1, ..., BN+1, satisfy the inequality

∑i=-1N+1|Bi(s)|≤10,s∈Ωs.

Proof. See [15]

Lemma 6. Let S(s) be the collocation approximation from the space of third-degree splines ξ3(Ω¯s) to the solution Û(s) of Equation (15) at the (j + 1)th time step. If g ∈ C2[0, 1], the parameter uniform error estimate is given by:

sup0<ε≤1max0≤i≤N|Û(si)-S(si)|≤Cℓ2,

where C is a positive constant independent of ε, μ, and N.

Proof. Let R(s) be the unique spline interpolate from ξ3(Ω¯N) to the solution Û(si) of the Equation (12) given by

R(s)=∑i=-1N+1α^iBi(s).

If g(s)∈C2(Ω^s), Û(si)∈C4(Ω¯s), and adopting the Hall error estimates [17] for this case, we get

‖(Û(s)-R(s))||∞≤c0‖Û(4)||∞ℓ4,‖(Û′(s)-R′(s))||∞≤c1‖Û(4)||∞ℓ3,‖(Û″(s)-R″(s))||∞≤c2‖Û(4)||∞ℓ2,

where the c0, c1, and c2 are the constants independent of ℓ and N. Therefore, from the estimates of the above equations and the collocating conditions Lε,γN,MS(si)=Lε,γMÛ(si)=g(si), it follows

            |ℒε,γN,MS(si)−ℒε,γN,MR(si)|=|ℒε,γMU^(si)−ℒε,γN,MR(si)|       =|ε(Û″(si)−σ(si)R″(si))+γa(si)(Û′(si)−R′(si))|                  +|−p(si)(Û(si)−R(si))|,=|ε(Û″(si)−σ(si)R″(si))+εσ(si)Û″(si)−εσ(si)Û″(si)|          +|γa(si)(Û′(si)−R′(si))|+|p(si)(Û(si)−R(si))|,              =|ε(σ(si)−1)Û″(si)+εσ(si)(Û″(si)−R″(si))|          +|γa(si)(Û′(si)−R′(si))|+|p(si)(Û(si)−R(si))|,          ≤|ε| |σ−1|‖Û″(s)‖∞+|ε| |σ|‖Û″(s)−R″(s)‖∞        +‖γa(s)‖∞‖Û′(s)−R′(s)‖∞+‖p(s)‖∞‖Û(s)−R(s)‖∞,    ≤|ε| |σ−1|‖Û(2)(s)‖∞+(|ε| |σ|c2ℓ2+‖μa(s)‖∞c1ℓ3                                                     +‖p(s)‖∞c0ℓ4)‖ ∞Û(4)‖.

Since ε ≪ 1 and as ℓ → 0, we obtain the following estimate

|Lε,γN,MS(si)-Lε,γN,MR(si)|≤|ε‖σ-1‖|Û(2)||∞.

From the power series expansion of hyperbolic cotangent function, we have scoths≈1+s23-s445 and then |scoths − 1| ≤ Cs2. This and Equation (14) indicate

|σ-1|≤Cℓ2.

Substituting Equation (19) in Equation (18), we get

‖Lε,γN,MS(s)-Lε,γN,MR(s)‖≤Cℓ2.

Using relationship in Equation (15), Lε,γN,M(Û(si)-R(si)) yields,

Y(α-α^)=(F-F^),

where α=αij+1,   α¯=α¯ij+1,

α-α^=(α0-α^0,α1-α^1,⋯,αN-α^N),t F=Dαij+H, F-F^=(ℓ2(g(s0)-ĝ(s0)),ℓ2(g(s1)-ĝ(s1)),⋯,ℓ2(g(sN)-ĝ(sN)))t.

Using Equation (20), we get

‖F-F^||∞≤Cℓ4.

It is observed that the coefficient matrix A is diagonally dominant as they gratify the subsequent relations::

|yi,i|-(|yi,i-1|+|yi,i+1|)=2piℓ2≥2(ν)ℓ2>0,as  y(s)≥ν>0.

In addition, from the first and last rows of Y, we have

|y0,0|-|y0,1|=36σ0+6p0ℓ, |yN,N|-|yN,N-1|=36σN-6pNℓ.

The above diagonal dominant property for smaller values of ℓ [18] gives

‖Y-1||∞≤Cℓ2.

From (20)-(22), we get

|αi-α^i|≤Cℓ2,0≤i≤N.

Now to estimate |α-1-α^-1| and |αN+1-α^N+1|, using the boundary conditions, we get

|α-1-α^-1|≤Cℓ2. |αN+1-α^N+1|≤Cℓ2.

Therefore

max-1≤i≤N+1|αi-α^i|≤Cℓ2.

The above inequality allows us to estimate: ||S(s)−R(s)||∞, so ||Û(s)−S(s)||∞. In specific,

S(s)-R(s)=∑i=-1N+1(αi-α^i)Bi(s).

Thus, substituting Equation (24) and Lemma 5 in (25), we get

|S(s)-R(s)|≤max-1≤i≤N+1|αi-α^i|∑i=-1N+1|Bi(s)|≤Cℓ2.

This gives

‖S(s)-R(s)||∞≤Cℓ2.

Using the triangle inequality rule, we obtain

‖Û(s)-S(s)||∞≤‖Û(s)-R(s)||∞+‖R(s)-S(s)||∞.

Using first term of Equation (17) and Equation (26) into Equation (27), we have

sup0<ε≤1max0≤i≤N|Û(si)-S(si)|≤Cℓ2.

Theorem 1. Let u(si, tj) and S(si, tj) be the solution of the problem (3) and the collocation approximation from the space ξ3(Ωs) to the solution u(si, tj), respectively. If g(s,tj)∈C2(Ω¯s),

‖u(si,tj)-S(si,tj)‖≤C((Δt)2+ℓ2). 4 Test examples, numerical computations, and discussions

Some numerical experiments have been presented using the following examples.

Example 1. Consider the following problem of the form (1) [5, 6]

{ε∂2u∂s2+γ(1+s)∂u∂s-u(s,t)-∂u∂t=-u(s,t-1)+16s2(1-s2),(s,t)∈Ωs×(0,2],u(s,t)=0,(s,t)∈Ωs×(-1,0],u(0,t)=0,u(1,t)=0,t∈(0,2].

Example 2. Consider the following problem of the form (1) [5, 6]

{ε∂2u∂s2+γ(1+s(1-s)+t2)∂u∂s-(1+5st)u(s,t)-∂u∂t=-u(s,t-1)+s(1-s)(et-1),(s,t)∈Ωs×(0,2],u(s,t)=0,(s,t)∈Ωs×(-1,0],u(0,t)=0,u(1,t)=0,t∈(0,2].

Since the exact solution for Examples (1) and (2) is not exist, we use the double mesh principle to compute the maximum absolute errors, for each ε and μ, using the following formula [15, 19]:

Eε,γN,M=max0≤i≤N;t∈[0,T]|UN,M(si,tj)-U2N,2M(si,tj)|,

where UN,M(si,tj) is the numerical solution with (N, M) mesh points, and U2N,2M(si,tj) is the numerical solution at the finer mesh with (2N, 2M) mesh points. The (ε, γ)-uniform absolute errors are calculated using the following formula:

EN,M=maxε,γEε,γN,M.

The numerical rate of convergence and the (ε, γ)-uniform rate of convergence are computed by using the following formulas, respectively.

ρε,γN,M=log2(Eε,γN,MEε,γ2N,2M)  and  ρN,M=maxε,γρε,γN,M.

The computed results by the proposed scheme and comparison with previous schemes for the considered problems are presented in Tables 2–5. From these tables, we observed that the presented method provides more accurate results than results in the study mentioned in the references [5] and [6]. 3-D view numerical simulation for Examples (1) and (2) is presented in Figures 1, 2, respectively. These figures show that as (ε, γ) goes very small, twin boundary layers formed at s = 0 and s = 1. We have plotted a log-log plot of the Eε,γN,M verses N in Figures 3A, B for (1)-(2), respectively. This figure shows that the scheme is stable and (ε, γ)−uniformly convergent.

Comparison of Eε,γN,M and ρε,γN,M for fixed γ = 10−4 and varying ε for Example (1) with the study mentioned in the reference [6] using Shishkin (S-) mesh and Bakhvalov-Shishkin (BS-) mesh.

ε↓ M = N = 32 64 128 256 512 1024
Proposed Scheme
10−4 9.1896e-005 4.0758e-005 1.1623e-005 2.5647e-006 6.3395e-007 1.5761e-007
1.1729 1.8101e 2.1801 2.0163 2.0080
10−6 4.4930e-05 1.1228e-005 2.8067e-006 1.0981e-006 5.5981e-007 1.8088e-007
2.0006 2.0002 1.3539 9.7200e-01 1.6299
10−8 4.4931e-05 1.1228e-005 2.8068e-006 7.0168e-007 1.7542e-007 4.3855e-008
2.0006 2.0001 2.0000 2.0000 2.0000
10−10 4.4931e-05 1.1228e-005 2.8068e-006 7.0168e-007 1.7542e-007 4.3855e-008
2.0006 2.0001 2.0000 2.0000 2.0000
Results using S-mesh [6]
10−4 2.8183e-003 1.4253e-003 7.1638e-004 3.5904e-004 1.7973e-004 -
0.9835 0.9925 0.9965 0.9983 0.9992
10−6 2.8274e-003 1.4272e-003 7.1717e-004 3.5952e-004 1.7999e-004 -
0.9863 0.9927 0.9962 0.9981 0.9991
10−8 2.8287e-003 1.4280e-003 7.1748e-004 3.5961e-004 1.8002e-004 -
0.9860 0.9930 0.9965 0.9982 0.9991
10−10 2.8286e-003 1.4280e-003 7.1747e-004 3.5960e-004 1.8002e-004 -
0.9861 0.9930 0.9965 0.9982 0.9991
Results using BS-mesh [6]
10−4 1.9776e-003 9.9679e-004 5.0025e-004 2.5059e-004 1.2544e-004 -
0.9892 0.9946 0.9973 0.9986 0.9993
10−6 1.9789e-003 9.9691e-004 5.0033e-004 2.5064e-004 1.2544e-004 -
0.9891 0.9945 0.9973 0.9986 0.9993
10−8 1.9789e-003 9.9691e-004 5.0033e-004 2.5064e-004 1.2544e-004 -
0.9891 0.9945 0.9973 0.9981 0.9993
10−10 1.9789e-003 9.9691e-004 5.0033e-004 2.5064e-004 1.2544e-004 -
0.9891 0.9945 0.9973 0.9986 0.9993

Comparison of Eε,γN,M, EN,M, and ρN,M for ϵ = 10−4 and varying γ for Example (1) with the study mentioned in the reference [6] using S- mesh and BS- mesh.

γ↓ N = 32 64 128 256 512 1024
M = 32 64 128 256 512 1024
Proposed Scheme [6]
10−4 9.1896e-005 4.0758e-005 1.1623e-005 2.5647e-006 6.3395e-007 1.5761e-007
10−6 4.4389e-005 1.1075e-005 2.7674e-006 1.1918e-006 6.0607e-007 1.9911e-007
10−8 4.4103e-005 1.0828e-005 2.6141e-006 7.2964e-007 3.6784e-007 1.8261e-007
10−10 4.4103e-005 1.0828e-005 2.6141e-006 7.2964e-007 3.6786e-007 1.8430e-007
10−12 4.4103e-005 1.0828e-005 2.6141e-006 7.2964e-007 3.6786e-007 1.8430e-007
EN,M 9.1896e-005 4.0758e-005 1.1623e-005 2.5647e-006 6.3395e-007 1.9911e-007
ρN,M 1.1729 1.8101 2.1801 2.0163 1.6708
Results using S-mesh [6]
EN,M 2.8183e-003 1.4253e-003 7.1638e-004 3.5904e-004 1.7973e-004 -
ρN,M 0.9835 0.9925 0.9965 0.9983 0.9992 -
Results using BS-mesh [6]
EN,M 1.9653e-003 9.9001e-004 4.9686e-004 2.4891e-004 1.2531e-4 -
ρN,M 0.9892 0.9946 0.9973 0.9986 0.9901 -

Comparison of Eε,γN,M and ρε,γN,M for γ = 10−3 and varying ε for Example (2) with the results mentioned in the reference [5].

ε↓ N = 32 64 128 256 512
M = 8 32 128 512 2048
10−0 8.8039e-005 3.0501e-005 8.1957e-006 2.0851e-006 4.5917e-007
1.5293 1.8959 1.9748 2.183
10−2 6.1690e-004 2.1201e-004 5.6926e-005 1.4481e-005 3.6251e-006
1.5409 1.8970 1.9749 1.9981
10−4 6.3732e-004 2.1895e-004 5.8789e-005 1.4955e-005 3.6850e-006
1.5414 1.8970 1.9749 2.0209
10−6 6.3813e-004 2.1945e-004 5.9048e-005 1.5001e-005 3.6970e-006
1.5400 1.8939 1.9768 2.0206
10−8 6.3813e-004 2.1945e-004 5.9049e-005 1.5002e-005 3.6970e-006
1.5400 1.8939 1.9768 2.0207
10−10 6.3813e-004 2.1945e-004 5.9049e-005 1.5002e-005 3.6970e-006
1.5400 1.8939 1.9768 2.0207
EN,M 6.3813e-004 2.1945e-004 5.9049e-005 1.5002e-005 3.6970e-006
ρN,M 1.5400 1.8939 1.9768 2.0207
Results [5]
EN,M 1.1161e-002 2.4750e-003 8.9349e-004 3.3001e-004 9.5237e-005
ρN,M 2.1729 1.4699 1.4369 1.7929

Comparison of Eε,γN,M, EN,M, ρε,γN,M, and ρN,M, for ϵ = 10−4 and varying γ for Example (2).

γ↓ N = 32 64 128 256 512
M = 8 32 128 512 2048
Proposed Scheme
10−4 6.3734e-004 2.1896e-004 5.8792e-005 1.4955e-005 3.7550e-006
1.5414 1.8970 1.9750 1.9937
10−6 6.3735e-004 2.1897e-004 5.8793e-005 1.4955e-005 3.7550e-006
1.5414 1.8970 1.9750 1.9937
10−8 6.3735e-004 2.1897e-004 5.8793e-005 1.4955e-005 3.7550e-006
1.5414 1.8970 1.9750 1.9937
10−10 6.3735e-004 2.1897e-004 5.8793e-005 1.4955e-005 3.7550e-006
1.5414 1.8970 1.9750 1.9937
EN,M 6.3735e-004 2.1897e-004 5.8793e-005 1.4955e-005 3.7550e-006
ρN,M 1.5414 1.8970 1.9750 1.9937

3-D view numerical solution profile for Example 1 at N = 512 = M. (A) γ = 10−6, ε = 10−1. (B) γ = 10−1, ε = 10−6.

3-D view numerical solution profile for Example 2 at N = 512 = M. (A) γ = 10−6, ε = 10−1. (B) γ = 10−1, ε = 10−6.

Loglog plot of the maximum point-wise errors for Examples 1 and 2 at different values of ε. (A) Example 1 with γ = 10−4. (B) Example 2 with γ = 10−3.

 5 Conclusion

In this study, second order numerical method for the two-parametric singularly perturbed time delay parabolic problem on a uniform mesh is presented. The problem is discretized by the Crank-Nicolson method in time and the third-degree B-spline method in space variable. The presented scheme is proven to be an (ε, γ)−uniformly convergent accuracy of order O((Δt)2 + ℓ2). To validate the theoretical findings, some test examples are considered. The computed results are in agreement with the theoretical investigations. Furthermore, the numerical results show that the presented scheme gives better results than available schemes in the literature.

 Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

 Author contributions

ITD: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review and editing. WGM: Conceptualization, Investigation, Methodology, Validation, Writing - review and editing. GDK: Investigation, Methodology, Validation, Writing - review and editing.

 Funding

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

The authors would like to thank the editor and referees 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 Chen J O'Malley RE Jr. On the asymptotic solution of a two-parameter boundary value problem of chemical reactor theory. SIAM J Appl Math. (1974) 26:717–29. 10.1137/0126064 DiPrima RC. Asymptotic methods for an infinitely long slider squeeze-film bearing. J Lubricat Technol. (1968) 90:173–83. 10.1115/1.3601534 Vasil'eva AB. Asymptotic methods in the theory of ordinary differential equations containing small parameters in front of the highest derivatives. USSR Comput Math Phys. (1963) 3:823–63. 10.1016/0041-5553(63)90381-1 Bhathawala PH Verma AP. A two-parameter singular perturbation solution of one dimension flow through unsaturated porous media. Proc Indian Natl Sci Acad. (1975) 43:380–4. Saini S Kumar S Kuldeep K Kumar M. A robust numerical method for a two-parameter singularly perturbed time delay parabolic problem. Comput Appl Math. (2020) 39. 10.1007/s40314-020-01236-1 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. 10.1007/s40995-019-00697-2 Negero NT. A uniformly convergent numerical scheme for two parameters singularly perturbed parabolic convection-diffusion problems with a large temporal lag. Results Appl Math. (2022) 16:100338. 10.1016/j.rinam.2022.100338 Negero NT. A parameter-uniform efficient numerical scheme for singularly perturbed time-delay parabolic problems with two small parameters. Partial Diff Equat Appl Math. (2023) 7:100518. 10.1016/j.padiff.2023.100518 Vigo-Aguiara J Natesan S. An efficient numerical method for singular perturbation problems. J Comput Appl Math. (2006) 192:132–41. 10.1016/j.cam.2005.04.042 Natesan S Vico-Acuiar J Ramanujam N. A numerical algorithm for singular perturbation problems exhibiting weak boundary layers. Comput Math Appl. (2003) 45:469–79. 10.1016/S0898-1221(03)80031-7 Vico-Acuiar J Natesan S. A parallel boundary value technique for singularly perturbed two-point boundary value problems. J Supercomput. (2004) 27:195–206. 10.1023/B:SUPE.0000009322.23950.53 Tian H. The exponential asymptotic stability of singularly perturbed delay differential equations with a bounded lag. J Math Anal Appl. (2002) 270:143–9. 10.1016/S0022-247X(02)00056-2 Prenter P. Spline and Variational Methods. New York, NY: Wiley (1975). Schoenberg I. On spline functions. J Appl Comput Mech. (1966). Daba IT Duressa GF. Collocation method using artificial viscosity for time dependent singularly perturbed differential–difference equations. Math Comput Simulat. (2022) 192:201–20. 10.1016/j.matcom.2021.09.005 O'Malley RE. Singular Perturbation Methods for Ordinary Differential Equations. Berlin: Springer (1991). 10.1007/978-1-4612-0977-5 Hall C. On error bounds for spline interpolation. J Approx Theory. (1968) 1:209–18. 10.1016/0021-9045(68)90025-7 Varah JM. A lower bound for the smallest singular value of a matrix. Linear Algebra Appl. (1975) 11:3–5. 10.1016/0024-3795(75)90112-3 Kadalbajoo MK Yadaw AS. Parameter-uniform finite element method for two-parameter singularly perturbed parabolic reaction-diffusion problems. Int J Comput Methods. (2012) 9:1250047. 10.1142/S0219876212500478