- Research Note
- Open access
- Published:

# An efficient numerical approach for singularly perturbed time delayed parabolic problems with two-parameters

*BMC Research Notes*
**volumeÂ 17**, ArticleÂ number:Â 158 (2024)

## Abstract

### Objectives

The main objective of this work is to design an efficient numerical scheme is proposed for solving singularly perturbed time delayed parabolic problems with two parameters.

### Results

The scheme is constructed via the implicit Euler and non-standard finite difference method to approximate the time and space derivatives, respectively. Besides, to enhance the accuracy and order of convergence of the method Richardson extrapolation technique is employed. Intensive numerical experimentation has been done on some model examples. Further, the layer behavior of the solutions is presented using graphs and observed to agree with the existing theories. Finally, error analysis of the scheme is done and observed that the proposed method is parameter uniform convergent with the order of convergence \(\left( (\Delta t )^2+h^2\right) \)

## Introduction

Many real-life problems in the fields of engineering and applied mathematics involve small positive parameter(s). Many researchers, scholars and engineers considered differential equations with one-parameter \(\varepsilon (0<\varepsilon \ll 1)\), multiplying the highest derivative of the equation. This type of equation is one-parameter singularly perturbed differential equation. A number of research papers which deal with asymptotic, analytical and numerical solutions of such type of ordinary and partial differential equations can be found in the literature. Asymptotic and numerical solutions of two-parameter singularly perturbed differential equations are not studied extensively like their one-parameter counterpart.

It is a well-known fact that the presence of small parameters in the differential equations exhibit boundary and/or interior layers in the solution. Oâ€™Malley introduced singularly perturbed two-parameter problems and examined asymptotic expansion of their solutions [1,2,3,4]. In the subsequent years, several numerical methods were developed to improve the accuracy of the asymptotic methods proposed by Oâ€™Malley and his co-researchers [5]- [8]. From many researchers in the literature, author in [9] has proposed a robust fitted operator finite difference method for a two-parameter singular perturbation problem. Recently, quadratic B-spline collocation method for two-parameter singularly perturbed problem on exponentially graded mesh was proposed in [10].

The numerical solution of two-parametric singularly perturbed parabolic partial differential equations with non-smooth data were studied in [11,12,13]. Parameter-uniform numerical methods have been proposed by different researchers to solve two-parametric singularly perturbed parabolic partial differential equations with smooth data. To mention those methods, implicit Euler method in time on a uniform mesh and upwind finite difference method on Shishkin mesh [14], Rotheâ€™s method for temporal discretization on a uniform mesh and finite element method for spatial discretization on a Shishkin mesh [15], an implicit Euler method in time and upwind finite difference method in space has been introduced using layer adapted equidistributed/moving meshes in [16], the classical implicit Euler method for time discretization and upwind scheme on the Shishkin-Bakhvalov mesh for spatial discretization was proposed in [17], an implicit Euler method for time discretization and a non-standard finite difference method on uniform mesh for spatial discretization was proposed in [18], the implicit Euler method for time stepping on a uniform mesh and a special hybrid monotone difference operator for spatial discretization on a specially designed piecewise uniform Shishkin mesh was considered in [19], an implicit Euler scheme on a uniform mesh in the temporal direction and the quadratic B-spline collocation scheme on an exponentially graded mesh in the spatial direction was constructed in [20], Crank-Nicolson scheme for the time derivative and cubic spline in tension for the spatial derivatives on a layer resolving non-uniform Bakhvalov-type mesh for a singularly perturbed two small parameters was presented in [22] to solve singularly perturbed two-parameter parabolic convection-diffusionâ€“reaction problems. The authors in [40] have developed robust weak Galerkin finite element method for two parameter singularly perturbed parabolic problems on nonuniform meshes. The scholars in [43] have designed two hybrid computational algorithms to find approximate solutions for singularly perturbed parabolic convection-diffusionâ€“reaction problems with two small parameters. All the aforementioned works for two-parameter singularly perturbed parabolic problems deals with initial-Dirichlet boundary conditions having no delay terms. However, the robust numerical solution of two-parameter singularly perturbed parabolic problems with delay terms are still limited.

Due to the wide application area of singularly perturbed problems with delay(s) have gained remarkable attention from researchers. For example, in population ecology, time delay represents the hatching period or duration of gestation; in genetic repression modeling, time delays play an important role in processes of transcription and translation as well as spatial diffusion of reactants and in control systems, delay terms account for the time delay in actuation and in information transmission and processing. Many other examples can be found in [34]. Also, the numerical treatment of such problems premeditated by many scholars. For instance, one-parameter singularly perturbed parabolic problems with time delay were studied in [23, 24, 36,37,38] and references therein. High-order finite difference technique for delay pseudo-parabolic equations is discussed in [41]. The authors [17] have established the finite difference scheme on Bakhvalov-type mesh for the singularly perturbed pseudo-parabolic problems with time-delay. Cimen and Amiraliyev [25] have suggested a uniform convergence method for singularly perturbed delay ordinary differential equation. Sumit et al. [26] have presented a robust numerical scheme using a hybrid monotone finite difference scheme on a rectangular mesh which is a product of uniform mesh in time and a layer-adapted Shishkin mesh in space for two-parametric singularly perturbed parabolic problems with time delay. Govindarao et al. [27] have 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 (standard Shishkin mesh, Bakhvalov-Shishkin mesh) for singularly perturbed two parameter time delay parabolic problem. Kumar and Kumar [45] Have constructed numerical method using a hybrid monotone finite difference scheme on a rectangular mesh which is a product of uniform mesh in time and a layer-adapted Shishkin mesh in space for singularly perturbed two-parameter parabolic partial differential equations with time delay. The authors in [39] have suggested the high order difference approximation with Identity Expansions technique to construct a second order finite difference scheme and combine this with standard backward Euler difference scheme in a special way on a piecewise-uniform Shishkin mesh to solve coupled system of singularly perturbed first order ordinary differential equations. Multigrid techniques to solve the linear systems arising from finite difference discretization on Shishkin meshes of 2D singularly perturbed problems have presented in [44]

The schemes in [26] and [27] 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 solve boundary value problems. For instance, the authors in [31,32,33] have suggested different EFD schemes for singularly perturbed two-point boundary-value problems.

Motivated by the proceeding works, in this work, we constructed the non-standard difference method together with the Richardson extrapolation for two-parameter singularly perturbed parabolic problem with time delay. Let \(\Omega _{s}^{N}=(0,1)\) and \(\Omega _{t}^{M}=(0,T]\) be the space and time domain respectively, where *T* is fixed time. Then, the rectangular domain is the tensor product defined by \(\Omega =\Omega _{s}^{N} \times \Omega _{t}^{M}\). Consider the following two-parameter singularly perturbed parabolic problem

subject to initial condition

and boundary conditions

where \(\varepsilon (0<\varepsilon \ll 1)\) and \(\mu (0\le \mu \ll 1)\) are two small perturbation parameters and \(\tau >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*) \(q_0(t), q_1(t)\) and \(\theta (s,t)\) are sufficiently smooth and bounded with \(a(s,t)\ge \alpha>0, b(s,t)\ge \beta >0\). The mathematical models related to two-parameter singularly perturbed problem of type (1) arises in transport phenomena in chemistry, biology, chemical reactor theory [4], lubrication theory [6] and dc motor theory [5] and flow through unsaturated porous media [8]. In two-parameter SPPs, the diffusion and convection terms are multiplied by the perturbation parameters.

Two-parameter singularly perturbed boundary value problems was first introduced by Oâ€™Malley, see [1,2,3,4]. He has shown that the layer behaviour for these problems depends not only on the parameters \(\varepsilon \) and \(\mu \), but significantly depends on the ratios of \(\varepsilon \) to different powers of \(\mu \). However, the ratio of \(\mu ^2\) to \(\varepsilon \) plays a significant role in the study of these problems. The particular cases, \(\mu = 0\) in which the parabolic type layers each of width O\((\sqrt{\varepsilon })\) appear at both lateral boundary of the domain, and \(\mu = 1\) in which an exponential layer of width O\((\varepsilon )\) appears near the left lateral boundary have been extensively considered analytically as well as numerically. For small \(\varepsilon \) and \(\mu \), the solution to the problem (1)â€“(3) exhibits boundary layers in the neighborhoods of \(s=0\) and \(s=1\).

**Notations:** Throughout this paper *N* and *M* denote the number of mesh points in *s* and *t* direction respectively. *C* denotes a generic positive constant independent of the singular perturbation parameters \(\varepsilon , \mu \) and the mesh sizes.

## Analytical aspects of the problem

### Lemma 1

(Continuous minimum principle). Assume that \(z(s,t)\in C^{(2,1)}({\bar{\Omega }})\) be sufficiently smooth function such that \(z(s,t)\ge 0\) on \((s,t)\in [0,1]\times [-\tau ,0]\), \(z(0,t)\ge 0\), \(z(1,t)\ge 0\) on \(0\le t\le T\), and \({\mathcal {L}}_{\varepsilon ,\mu }z(s,t)\le 0, (s,t)\in \Omega \). Then, \(z(s,t)\ge 0 \), \((s,t)\in \bar{\Omega }\).

### Proof

See [35] \(\square \)

The following Lemma proves the stability estimate to obtain unique solution.

### Lemma 2

(Uniform stability estimate) Let *z*(*s*,Â *t*) be the solution to the continuous problem (1)â€“(3), then it satisfies the bound

where \(\left\| . \right\| _{\overline{\Omega }} \) is used to denote maximum norm given by \(\left\| z\right\| _{\overline{\Omega }} =\max _{(s,t)\in \overline{\Omega }} |z(s,t)|. \)

### Proof

See [35]. \(\square \)

## Description of the numerical method

We first discretized the time direction using an implicit Euler method with uniform step size \(\Delta t\) which leads to a system of boundary value problem. Then, the discretization of space direction is made using the non-standard finite difference method.

### Time semi-discretization

We have two intervals \([-\tau , 0]\) and [0,Â *T*] on the time direction and we use the uniform mesh with time step \(\Delta t\)

where *M* is number of mesh points in t-direction in the interval [0,Â *T*] and *m* is the number of mesh points in \([-\tau , 0]\). The step length \(\Delta t\) satisfies \(m\Delta t = \tau \), where *m* is a positive integer and \(t_{j} = j\Delta t, j \ge - m\). To discretize the time variable for Eq. (1), we use the implicit Euler method, which is given by

subject to semi-discrete initial and boundary conditions

where \(d(s,t_{j+1})=b(s,t_{j+1})+\frac{1}{\Delta t}\) and \(g(s,t_{j+1})=-c(s,t_{j+1})Z(s,t_{j+1-m})+f(s,t_{j+1})-\frac{Z(s,t_{j})}{\Delta t}\). By using the initial condition, we can evaluate the right-hand side as

The local truncation error of the time semi-discretization is given by \(e_{j+1}=z(s,t_{j+1})-Z(s,t_{j+1})\), where \(Z(s,t_{j+1})\) is the solution of the following boundary value problem

with the boundary conditions

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

### Lemma 3

(Local error estimate) If

the local error estimate of the time discretization is given by

### Proof

One can find the proof of lemma in [27].\(\square \)

The global error is the measure of the contribution of the local error estimate at each time step and is given by \( E_{j}=z(s,t_j)-Z(s,t_j).\)

### Lemma 4

Under Lemma (3), the global error estimate at \(t_{j}\) is given by

We conclude that time semi-discretization is first-order uniformly convergent. The \(\varepsilon -\)convergence analysis of the numerical scheme which we propose requires that we use bounds on the solution and its derivatives. The solutions of the characteristic equation for time semi-discrete problem

are \(r_0(s)<0\) and \(r_1(s)>0\) which are used to describe the boundary layers at \(s=0\) and \(s=1\), respectively. The quantities \(\mu _0\) and \(\mu _1\) are defined as \(\mu _0=-\max \limits _{[0,1]}r_0(s)\) and \(\mu _1=\max \limits _{[0,1]}r_1(s)\).

### Remark 1

The situations of two external layers are characterized by \(\mu ^2\ll 1\) or \(\mu ^2/\varepsilon \rightarrow 0\) as \(\varepsilon \rightarrow 0\), which implies that \(\mu _0 \approx \mu _1 \approx \min \limits _{[0,1]} \sqrt{\frac{(b(s,t_{j+1})+1/\Delta t)}{\varepsilon }}\) and we have boundary layers at \(s=0\) and \(s=1\). The boundary layer at \(s=0\) is encountered in the case when \(\varepsilon \ll \mu ^2\) as \(\mu \rightarrow 0\). In this case, \(\mu _1\approx 0\) and \(\mu _0 \approx \min \limits _{[0,1]} \frac{\mu a(s,t_j)}{\varepsilon }\).

The next theorem, gives the \(\varepsilon \)-uniform bounds for the derivatives of the solution *u* with respect to *x*, needs to study the uniform convergence at spatial discretization.

### Theorem 1

Up to a certain order *k* that depends on the smoothness of of the functions *a*(*s*,Â *t*),Â *b*(*s*,Â *t*),Â *f*(*s*,Â *t*) and for any real constant \(p\in (0,1)\), we have the following bound

### Proof

For the details of the proof, see [15]. \(\square \)

### Spatial discretization

In this section, we discretize the spatial variable of (4) using the nonstandard finite difference method of Mickens [29] as discussed below. On the spatial domain [0,Â 1], we introduce the equidistant meshes with uniform mesh length *h* such that \(s_{0}=0, s_{i}=ih,i=1(1)N-1,s_{N}=1,h=\frac{1}{N}\), where *h* is the step size and *N* is the number of mesh points in the spatial direction. The spatial Discretization of (4) yields

We use the notation \(Z(s_i,t_{j+1})\equiv Z_i\) as the approximation of \(z(s_i,t_{j+1})\equiv z_i\) for the sake of simplicity. From the theory of non-standard finite difference method, we can discretize (8) in space to obtain the discrete problem in the form

where \(d_i=b_i+\frac{1}{\Delta t}\) and \(g_i=-c_iZ_i^{-m}+f_i-\frac{Z_i}{\Delta t}\) for \(i=1,...,N-1, \;\; j=0,...,M.\) According to Mickens [28, 29], the concept of sub-equations is the major tool to derive the denominator function for the differential equation. From (9), we take the homogeneous form of the constant coefficient sub-equation

where \(\delta =\mu a_i\) and \(\eta =d_i=b_i+\frac{1}{\Delta t}\). Equation (10) has two linearly independent analytical solutions, namely, \(\exp (\lambda _{1} x)\) and \(\exp (\lambda _{2} x)\), where

Following Mickenâ€™s, we construct the second-order difference equation for (10) as follows

Simplifying the determinant in (12), we get

which is the exact scheme for (10) in the sense that (13) has the same general solution

as (10). It is noted that to construct the non-standard finite difference method, we need to find a suitable denominator function which replaces \(h^2\). To this end, the extraction of the denominator function from (13) is not straightforward. However, the fact that the layer behaviors of the solution of problem (1) and that of the problem (10) in the case when \(\eta \equiv 0\) are similar [9]. Thus, for the latter case, that is, \(\eta \equiv 0\), we have the exact scheme of the form

from hyperbolic identity. Multiplying both sides of equation (15) by \(\exp (\frac{\delta h}{2\varepsilon })\), we have

Adding the term \((Z_{i+1}+Z_{i})\) and subtracting it and after some manipulations, we have

From (16), we have

where the denominator function is given by

Adopting the denominator function for the variable coefficient, we can write as

Thus, we get the following fully discrete problem

with the following discrete initial and boundary conditions

where the denominator function is given as

The discrete scheme in (20) and the discrete conditions in (21) can be written in matrix form as

where *Z* and *G* are column vectors of \(N-1\) and the matrix *W* is a tri-diagonal matrix of \((N-1)\times (N-1)\). The entries of the coefficient matrix *W* are given by

The entries of column vectors *G* and *Z* are given as follows

Next, we prove some useful attributes the discrete scheme in (22). The discrete operator \({\mathcal {L}}_{\varepsilon ,\mu }^{N,M}\) defined in (9) satisfies the following discrete minimum principle.

### Theorem 2

Assume that \({\mathcal {L}}_{\varepsilon ,\mu }^{N,M}\) be discrete operator given in (9) and \(\Pi _i^j\) be any mesh function that satisfies the initial condition \(\Pi _i^{-j}\ge 0, 1\le i \le N-1, \;\; 0\le j \le m\) and boundary conditions \(\Pi _0^j\ge 0, \; \Pi _N^j\ge 0\), \(0\le j \le M\). If \({\mathcal {L}}_{\varepsilon ,\mu }^{N,M}\Pi _i^j\le 0\) for all \((i,j)\in \Omega ^{N,M}\), then \(\Pi _i^j\ge 0\) in \(\bar{\Omega }^{N,M}.\)

### Proof

Let *s* and *p* be indices such that \(\Pi _{s}^{p}= \min \limits _{{\forall (i,j)}}\Pi _i^j\) for \(\Pi _i^j\in \bar{\Omega }^{N,M}\). Assume that \(\Pi _{s}^{p}<0\). Then, it is easy to see that \((s,p) \in \{1,\cdots ,N-1\}\times \{1,\cdots ,M\}\), because otherwise \(\Pi _{s}^{p}\ge 0\). It follows that \(\Pi _{s+1}^p-\Pi _s^p\ge 0\) and \( \Pi _s^p-\Pi _{s-1}^p\le 0\). Therefore, now

which is a contradiction and thus, the assumption \(\Pi _s^{p+1}<0, ~\forall (s,p)\) is wrong. Thus, \(\Pi _s^{p+1}>0\) implies that \(\Pi _i^j\ge 0\) in \(\bar{\Omega }^{N,M}.\) \(\square \)

Using this discrete minimum principle, we now show that the present method also satisfies the uniform stability estimate given in the following lemma.

### Lemma 5

The discrete operator \({\mathcal {L}}_{\varepsilon ,\mu }^{N,M}\) is uniformly stable, in the sense that if \(P_i^{j+1}\) is any mesh function such that \(P_0^{j+1}=P_N^{j+1}=0\), then

### Proof

Let define the two barrier functions \((\Psi ^{\pm })_i^j\) by \((\Psi ^{\pm })_i^j=Y \pm P_i^{j+1},\) where

We have \((\Psi ^{\pm })_0^{j+1}=Y \pm P_0^{j+1}=Y \pm q_0(t_{j+1})\ge 0\), and \((\Psi ^{\pm })_N^{j+1}=Y \pm P_N^{j+1}=Y \pm q_1(t_{j+1})\ge 0.\)

On the discretized domain \(1\le i \le N-1\), we have

Using the fact that where \(0<\beta \le b_i^{j+1}<(b_i^{j+1}+1/\Delta t)\), we have \({\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(\Psi ^{\pm })_i^{j+1}\le 0\). By the discrete minimum principle in Theorem (2), we obtain \((\Psi ^{\pm })_i^{j+1}\ge 0,\; 0\le i \le N.\) \(\square \)

### Convergence analysis

We use the following lemma to prove the uniform convergence analysis of the discrete problem (9).

### Lemma 6

For all positive integers *k* on a fixed mesh, we have

where \(x_{i}=ih, h=\frac{1}{N}, \hspace{0.2cm} \forall i=1,...,N-1.\)

### Proof

For the proof, see [30]. \(\square \)

The following analysis concerns the space variable *x*. The local truncation error in the spatial variable of the proposed method is given by

Simplifying the above expression, we obtain

Taylor series expansions of the terms \(z_{i}^{j+1}, z_{i+1}^{j+1}\) and \(z_{i-1}^{j+1}\) on space direction are given as following

Adding the first two equations in the above Taylor series expansion, we deduce the following

From the first Taylor series expansion, we have

Substituting equations (26)-(27) in their respective positions into (25) and rearranging, we obtain

Simplifying (28), we obtain

Using a truncated Taylor series expansion of the denominator function [18]

in (30) gives

Applying the relation \(h>h^2\) in (31) and using the bounds on derivatives in Theorem (1) together with Lemma (6) and using the fact that as \(\varepsilon \rightarrow 0\), both the terms \(\mu _0^k\exp (-p\mu _0s_i)\) and \(\mu _1^k\exp (-p\mu _1(1-s_i))\rightarrow 0\) for all \(k\in \{0,1,2,\cdots ,\}\). Thus, the discrete problem satisfies the following bound

where *C* is constant independent of the perturbation parameters and mesh sizes. Invoking the uniform stability in Lemma (5), we obtain the result

Therefore, main convergence of the fully discretized scheme is given in the following theorem.

### Theorem 3

Let \(z_i^{j+1}\in C^{4,2}(\bar{\Omega })\) be the solution to problem in (1)-(3)) and \(Z_i^{j+1}\) be the solution to discrete problem in (20) with its discrete boundary conditions in (21). Then, the overall error bound satisfies

From the above theorem, we deduce that the developed method is first-order convergent, independent of the parameters \(\varepsilon \) and \(\mu \) both in space and time directions.

Next, we use Richardson extrapolation to boost the accuracy and rate of convergence for the proposed method.

### Richardson extrapolation

Richardson extrapolation is a technique used to accelerate the accuracy and rate of convergence of the proposed method. From (32), we have

where \(z(s_i,t_{j+1})\) and \(Z_i^{j+1}\) are exact and approximate solutions, respectively. Assume \(\Omega ^{N,M}\subset \Omega ^{2N,2M}\) where \(\Omega ^{N,M}\) is the mesh obtained from the mesh intervals *h* and \(\Delta t\) and \(\Omega ^{2N,2M}\) is the mesh obtained by bisecting the mesh intervals *h* and \(\Delta t\). Denoting the numerical solution obtained with the mesh points \(\Omega ^{2N,2M}\) by \(\tilde{Z}_i^{j+1}\). From (32), it is clear that for the mesh \((s_i,t_{j+1})\in \Omega ^{N,M}\)

As \({\tilde{s}}_{i+1}-{\tilde{s}}_{i}={\tilde{h}}=\frac{h}{2}\) for \({\tilde{s}}_i\in \Omega ^{2N}\) and \({\tilde{t}}_{j+1}-{\tilde{t}}_{j}=\tilde{\Delta }t=\frac{\Delta t}{2}\) for \({\tilde{t}}_{j+1}\in \Omega ^{2\,M}.\) For the mesh \(({\tilde{s}}_i,{\tilde{t}}_{j+1})\in \Omega ^{2N,2M}\), we have

where \(R^{N,M}\) and \(R^{2N,2\,M}\) are the remainder terms with the truncation error of \(O(h^2+\Delta t^2)\). Combining the inequalities in (34) and (35) gives us with \(z(s_i,t_{j+1})-(2\tilde{Z}_i^{j+1}-z(s_i,t_{j+1}))\le C(h^2+\Delta t^2)\), which yields that

is also an extrapolated numerical solution. Therefore, we have the error bound for extrapolated solution summarized in the theorem as follows.

### Theorem 4

Let \(z_i^{j+1}\) be the solution to the continuous problem (1) and (2) and \((Z_i^{j+1})^{extr}\) be the extrapolated solution. Then, the new error bound takes the form

### Proof

The proof is given in [21]. \(\square \)

Therefore, using Richardson extrapolation, first-order uniformly convergent method is changed into second-order uniformly convergent. Thus, the proposed method is second-order convergent.

## Test examples, numerical computations and discussions

In this section, we carry out numerical experiments to corroborate the performance of the proposed method with the theoretical results discussed in the previous sections. Since the exact solution for the Examples (1) and (2) are not available, we use the double mesh principle to calculate maximum absolute errors, for each \((\varepsilon ,\mu )\), using the following formula

before extrapolation (BE) and after extrapolation (AE), we use the formula

where \(Z^{N,M}(s_i,t_j)\) is the numerical solution with (*N*,Â *M*) mesh points and \(Z^{2N,2M}(s_i,t_j)\) is the numerical solution at the finer mesh with (2*N*,Â 2*M*) mesh points before extrapolation (BE). The numerical solutions after extrapolation (AE) are \((Z^{N,M})^{extr}(s_i,t_j)\) using the mesh points (*N*,Â *M*) with mesh sizes *h* and \(\Delta t\) and \((Z^{2N,2M})^{extr}(s_i,t_j)\) using the mesh points (2*N*,Â 2*M*) with mesh sizes \(\frac{h}{2}\) and \(\frac{\Delta t}{2}\). The \((\varepsilon ,\mu )\)-uniform errors before and after extrapolation were calculated using the following formulas, respectively

Furthermore, we compute the numerical rate of convergence before and after extrapolation with the following formulas, respectively

The \((\varepsilon ,\mu )\)-uniform rate of convergence before and after extrapolation were calculated using the following formulas, respectively

### Example 1

Consider singularly perturbed two-parameter parabolic problem [26, 27]

### Example 2

Consider singularly perturbed two-parameter parabolic problem [26, 27]

The computed \((e_{\varepsilon , \mu }^{N,M}),\; (e_{\varepsilon , \mu }^{N,M})^{extr}, \;(e^{N,M}),\; (e^{N,M})^{extr}, \; \rho _{\varepsilon , \mu }^{N,M} \), \((\rho _{\varepsilon , \mu }^{N,M})^{extr}\; \rho ^{N,M} \), and \((\rho ^{N,M})^{extr}\) for Examples (1) and (2) are tabulated in Tables 1, 2, 3, 4, 5, 6, 7 and 8 for different values of \(\varepsilon , \mu \) and mesh points. From these tables, one can observe that the results obtained after extrapolation provides more accurate results obtained before extrapolation and results in [26, 27]. From the table of values, we deduce that when the mesh points increases the maximum absolute errors decreases. Numerical simulation for Examples (1 and 2) are displayed in Figs. (1 and 2), respectively. From these figures, we observe that as \((\varepsilon ,\mu )\) goes very small a twin boundary layers are created at \(s=0\) and \(s=1\). For a visual understanding of the theoretical order of convergence graphically, the maximum absolute errors for Examples (1) and (2) are plotted using log-log scale in Figs. (3 and 4), respectively.

## Conclusion

In this study, a robust numerical method for the two-parametric singularly perturbed time-delayed parabolic problem on a uniform mesh is presented. The problem is discretized by an implicit Euler method in the time variable and the non-standard finite difference method in the space variable. The method is analyzed for parameter uniform convergence. To boost the accuracy, the Richardson extrapolation technique has been applied. The numerical solutions displayed in the Tables show that the present method is parameter uniform convergence of second-order and it agrees with the theoretical order of convergence. The performance of the proposed method is examined by comparing the results with those of previous studies. It has been found that the proposed scheme provides more accurate and stable results. To substantiate the suitability of the proposed method, graphs have been plotted for the two examples by taking different values of the parameters \(\varepsilon \) and \(\mu \). The drawback of the proposed method is that difficult to apply to higher-order singular perturbation problems. The proposed method is easy to implement and, with a little modification, can easily be extended to nonlinear, discontinuous data, and other families of the problem under consideration.

## Data availability

Not applicable.

## References

Oâ€™Malley RE Jr. Singular perturbations of boundary value problems for linear ordinary differential equations involving two-parameters. J Math Anal Appl. 1967;19(2):291â€“308.

Oâ€™Malley RE Jr. Two-parameter singular perturbation problems for second-order equations. J Math Mech. 1967;16(10):1143â€“64.

Oâ€™Malley RE Jr. Boundary value problems for linear systems of ordinary differential equations involving many small parameters. J Math Mech. 1969;18:835â€“56.

Chen J, Oâ€™Malley RE Jr. On the asymptotic solution of a two-parameter boundary value problem of chemical reactor theory. J SIAM Appl Math. 1974;26(4):717â€“29.

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(4):823â€“63.

DiPrima RC. Asymptotic methods for an infinitely long slider squeeze-film bearing. J Lubr Technol. 1968;90(1):173â€“83.

Bigge J, Bohl E. Deformations of the bifurcation diagram due to discretization. Math Comput. 1985;45(172):393â€“403.

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(5):380â€“4.

Patidar KC. A robust fitted operator finite difference method for a two-parameter singular perturbation problem. J Differ Equ Appl. 2008;14(12):1197â€“214.

Shivhare M, Podila PC, Kumar D. Quadratic B-spline collocation method for two-parameter singularly perturbed problem on exponentially graded mesh. Int J Comput Math. 2021;98(12):2461â€“81.

Chandru M, Prabha T, Das P, Shanthi V. A numerical method for solving boundary and interior layers dominated parabolic problems with discontinuous convection coefficient and source terms. Differ Equ Dynam Syst. 2019;27(1):91â€“112.

Chandru M, Das P, Ramos H. Numerical treatment of two-parameter singularly perturbed parabolic convection diffusion problems with non-smooth data. Math Methods Appl Sci. 2018;41(14):5359â€“87.

Kumar D, Kumari P. Uniformly convergent scheme for two-parameter singularly perturbed problems with non-smooth data. Numer Methods Partial Differ Equ. 2021;37(1):796â€“817.

Oâ€™Riordan E, Pickett ML, Shishkin GI. Parameter-uniform finite difference schemes for singularly perturbed parabolic diffusion-convection-reaction problems. Math Comput. 2006;75(255):1135â€“54.

Kadalbajoo MK, Yadaw AS. Parameter-uniform finite element method for two-parameter singularly perturbed parabolic reaction-diffusion problems. Int J Comput Methods. 2012;09(04):1250047.

Das P, Mehrmann V. Numerical solution of singularly perturbed convection-diffusion-reaction problems with two small parameters. BIT Numer Math. 2016;56(1):51â€“76.

Jha A, Kadalbajoo MK. A robust layer adapted difference method for singularly perturbed two-parameter parabolic problems. Int J Comput Math. 2015;92(6):1204â€“21.

Munyakazi JB. A robust finite difference method for two-parameter parabolic convection-diffusion problems. Appl Math Inf Sci. 2015;9(6):2877â€“83.

Gupta V, Kadalbajoo MK, Dubey RK. A parameter-uniform higher order finite difference scheme for singularly perturbed time-dependent parabolic problem with two small parameters. Int J Comput Math. 2019;96(3):474â€“99.

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(1):186â€“215.

Bullo TA, Duressa GF, Delga GA. Robust finite difference method for singularly perturbed two-parameter parabolic convection-diffusion problems. Int J Comput Methods. 2021;18(02):2050034.

Mekonnen TB, Duressa GF. A fitted mesh cubic spline in tension method for singularly perturbed problems with two parameters. Int J Math Math Sci. 2022;2022:11.

Gelu FW, Duressa GF. A uniformly convergent collocation method for singularly perturbed delay parabolic reaction-diffusion problems. Abstract Appl Anal. 2021;2021:11.

Govindarao L, Mohapatra J, Das A. A fourth-order numerical scheme for singularly perturbed delay parabolic problem arising in population dynamics. J Appl Math Comput. 2020;63:171â€“95.

Cimen E, Amiraliyev GM. Uniform convergence method for a delay differential problem with layer behaviour. Mediterr J Math. 2019;16(3):57,1-15.

Sumit Kumar S, Kuldeep Kumar M. A robust numerical method for a two-parameter singularly perturbed time delay parabolic problem. Comput Appl Math. 2020;39(209).

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(5):2373â€“83.

Mickens RE. Nonstandard finite difference models of differential equations. Singapore: World Scientific; 1994.

Mickens RE. Advances in the applications of nonstandard finite difference schemes. Singapore: World Scientific; 2005.

Patidar KC, Sharma KK. Uniformly convergent nonstandard finite difference methods for singularly perturbed differential-difference equations with delay and advance. Int J Numer Methods Eng. 2006;66(2):272â€“96.

Vigo-Aguiara J. An efficient numerical method for singular perturbation problems. J Comput Appl Math. 2006;192:132â€“41.

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.

Vico-Acuiar J, Natesan S. A parallel boundary value technique for singularly perturbed two-point boundary value problems. J Supercomput. 2004;27:195â€“206.

Wu J. Theory and applications of partial functional differential equations, vol. 119. Berlin: Springer; 2012.

Daba IT, Melesse WG, Kebede GD. A fitted numerical approach for singularly perturbed two-parameter parabolic problem with time delay. Comput Math Methods. 2023.

Kumar S, Kumar M. High order parameter-uniform discretization for singularly perturbed parabolic partial differential equations with time delay. Comput Math Appl. 2014;68(10):1355â€“67.

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(5):1849â€“66.

Sahoo SK, Gupta V. Parameter robust higher-order finite difference method for convection-diffusion problem with time delay. Numer Methods Partial Differ Equ. 2023;39(6):4145â€“73.

Chandra S, Sekhara Rao, Kumar S. Second order global uniformly convergent numerical method for a coupled system of singularly perturbed initial value problems. Appl Math Comput. 2012.

Singh J, Kumar N, Jiwari R. A robust weak Galerkin finite element method for two parameter singularly perturbed parabolic problems on nonuniform meshes. J Comput Sci. 2024:102241.

Amiraliyev GM, Cimen E, Amirali I, Cakir M. High-order finite difference technique for delay pseudo-parabolic equations. J Comput Appl Math. 2016;321:1â€“7.

Gunes B, Duru H. A computational method for the singularly perturbed delay pseudo-parabolic differential equations on adaptive mesh. Int J Comput Math. 2023;100(8):1667â€“82.

Ansari KJ, Izadi M, Noeiaghdam S. Enhancing the accuracy and efficiency of two uniformly convergent numerical solvers for singularly perturbed parabolic convection- diffusion-reaction problems with two small parameters. Demonstratio Math. 2024;57:20230144.

Gaspar FJ, Clavero C, Lisbona F. Some numerical experiments with multigrid methods on Shishkin meshes. J Comput Appl Math. 2002;138(1):21â€“35.

Sumit SK, Kuldeep MK. A robust numerical method for a two-parameter singularly perturbed time delay parabolic problem. Comput Appl Math. 2020;39(3):209.

## Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

## Author information

### Authors and Affiliations

### Corresponding author

## Ethics declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declared no potential Conflict of interest concerning the research, authorship and publication of this article.

## Additional information

### Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the articleâ€™s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the articleâ€™s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

## About this article

### Cite this article

Daba, I.T., Melesse, W.G., Gelu, F.W. *et al.* An efficient numerical approach for singularly perturbed time delayed parabolic problems with two-parameters.
*BMC Res Notes* **17**, 158 (2024). https://doi.org/10.1186/s13104-024-06813-9

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s13104-024-06813-9