Skip to main content

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

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) \)

Peer Review reports

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

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu } z(s,t)&{}\equiv \varepsilon \frac{\partial ^{2} z(s,t)}{\partial s^{2}}+\mu a(s,t)\frac{\partial z(s,t)}{\partial s}-b(s,t)z(s,t)-\frac{\partial z(s,t)}{\partial t}\\ {} &{}=-c(s,t)z(s,t-\tau )+f(s,t), \;\; (s,t)\in \Omega , \end{aligned} \end{array}\right. } \end{aligned}$$
(1)

subject to initial condition

$$\begin{aligned} z(s,t)=\theta (s,t), \;\;\; (s,t)\in [0,1]\times [-\tau ,0], \end{aligned}$$
(2)

and boundary conditions

$$\begin{aligned} \begin{aligned} z(0,t)=q_0(t), \;\;\;\; 0\le t\le T,\\ z(1,t)=q_1(t), \;\;\;\; 0\le t\le T, \end{aligned} \end{aligned}$$
(3)

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(st), b(st), c(st), f(st) \(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(st) be the solution to the continuous problem (1)–(3), then it satisfies the bound

$$\begin{aligned} \Vert z \Vert _{{\bar{\Omega }}} \le \max \big \{ \vert \theta (s,t) \vert , \vert q_0(t) \vert , \vert q_1(t) \vert \big \}+\beta ^{-1}\Vert f\Vert , \end{aligned}$$

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\)

$$\begin{aligned} \begin{array}{ll} \Omega ^{M}_{t}=\big \{t_{j}:t_{j}=j\Delta t, \; j=0,\cdots ,M, t_{M}=T, \Delta t=T/M\big \}, \\ \Omega ^{m}_{t}=\big \{t_{n}:t_{n}=n\Delta t, \; n=0,\cdots ,m, t_{m}=\tau , \Delta t=\tau /m \big \}, \end{array} \end{aligned}$$

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

$$\begin{aligned} {\mathcal {L}}^M_{\varepsilon ,\mu }Z^{j+1}\equiv \varepsilon Z_{ss}(s,t_{j+1})+\mu a(s,t_{j+1})Z_s(s,t_{j+1})-d(s,t_{j+1})Z(s,t_{j+1})=g(s,t_{j+1}), \end{aligned}$$
(4)

subject to semi-discrete initial and boundary conditions

$$\begin{aligned} {\left\{ \begin{array}{ll} Z(s,t_{j+1})=\theta (s,t_{j+1}), &{} (s,t_{j+1})\in [0,1]\times [-\tau ,0],\\ Z(0,t_{j+1})=q_0(t_{j+1}), &{} 0\le t_{j+1}\le T,\\ Z(1,t_{j+1})=q_1(t_{j+1}), &{} 0\le t_{j+1}\le T, \end{array}\right. } \end{aligned}$$
(5)

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

$$\begin{aligned} g(s,t_{j+1})= {\left\{ \begin{array}{ll} \frac{Z(s,t_j)}{\Delta t}-c(s,t_{j+1})\theta (s,t_{j+1-m})+f(s,t_{j+1}), &{} j=0,1,...,m,\\ \frac{Z(s,t_j)}{\Delta t}-c(s,t_{j+1})Z(s,t_{j+1-m})+f(s,t_{j+1}), &{} j=m+1,...,M. \end{array}\right. } \end{aligned}$$

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

$$\varepsilon {Z_{ss}}\left( {s,{t_{j + 1}}} \right) + \mu a\left( {s,{t_{j + 1}}} \right){Z_s}\left( {s,{t_{j + 1}}} \right) - d\left( {s,{t_{j + 1}}} \right)Z\left( {s,{t_{j + 1}}} \right) = g\left( {s,{t_{j + 1}}} \right),$$
(6)

with the boundary conditions

$$\begin{aligned} {\left\{ \begin{array}{ll} Z(0,t_{j+1})=q_0(t_{j+1}), &{} 0\le t_{j+1}\le T,\\ Z(1,t_{j+1})=q_1(t_{j+1}), &{} 0\le t_{j+1}\le T. \end{array}\right. } \end{aligned}$$
(7)

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

Lemma 3

(Local error estimate) If

$$\begin{aligned} \bigg |\frac{\partial ^k Z(s,t)}{\partial s^k} \bigg |\le C, \hspace{1cm} (s,t)\in \bar{\Omega }, ~~ 0\le k\le 2, \end{aligned}$$

the local error estimate of the time discretization is given by

$$\begin{aligned} \Vert e_{j+1} \Vert _\infty \le C\Delta t^2, \;\;\; 1\le j \le M. \end{aligned}$$

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

$$\begin{aligned} \Vert E_{j}\Vert _{\infty } \le C\Delta t, ~~~ j\le T/\Delta t. \end{aligned}$$

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

$$\begin{aligned} \varepsilon r^2(s,t_j)+\mu a(s,t_{j+1})r(s,t_{j+1})-(b(s,t_{j+1})+1/\Delta t)=0 \end{aligned}$$

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(st), b(st), f(st) and for any real constant \(p\in (0,1)\), we have the following bound

$$\begin{aligned} \bigg \Vert \frac{\partial ^{k} z(s)}{\partial s^k} \bigg \Vert _{\bar{\Omega }} \le C\big (1+\mu _0^ke^{-p\mu _0s}+\mu _1^ke^{-p\mu _1(1-s)}\big ), \hspace{0.5cm} 0\le k \le q. \end{aligned}$$

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

$$\begin{aligned} \begin{array}{c} \varepsilon \frac{\partial ^{2} Z(s_i,t_{j+1})}{\partial s^{2}}+\mu a(s_i,t_{j+1})\frac{\partial Z(s_i,t_{j+1})}{\partial s}-b(s_i,t_{j+1})Z(s_i,t_{j+1})-\frac{Z(s_i,t_{j+1}-Z(s_i,t_{j}}{\Delta t}\\ =-c(s_i,t_{j+1})Z(s_i,t_{j+1-m})+f(s_i,t_{j+1}). \end{array} \end{aligned}$$
(8)

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

$$\begin{aligned} \begin{array}{ll} {\mathcal {L}}^{N,M}_{\varepsilon ,\mu }Z_i\equiv \varepsilon \bigg (\frac{ Z_{i-1}-2Z_i+Z_{i+1}}{\gamma ^2_i}\bigg )+\mu a_i\bigg (\frac{Z_{i+1}-Z_i}{h}\bigg ) -d_iZ_i=g_i. \end{array} \end{aligned}$$
(9)

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

$$\begin{aligned} \varepsilon \bigg (\frac{ Z_{i-1}-2Z_i+Z_{i+1}}{\gamma ^2_i}\bigg )+\delta \bigg (\frac{Z_{i+1}-Z_i}{h}\bigg ) -\eta Z_i=0, \end{aligned}$$
(10)

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

$$\begin{aligned} \lambda _{1,2}=\frac{-\delta \pm \sqrt{\delta ^{2}-4\varepsilon \eta }}{2\varepsilon }. \end{aligned}$$
(11)

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

$$\begin{aligned} {\left| \begin{array}{ccc} Z_{i-1}&{}Z_{1},_{i-1}&{}Z_{2},_{i-1}\\ Z_{i}&{}Z_{1},_{i}&{}Z_{2},_{i}\\ Z_{i+1}&{}Z_{1},_{i+1}&{}Z_{2},_{i+1} \end{array} \right| } =\left| \begin{array}{ccc} Z_{i-1}&{}\exp (\lambda _{1} x_{i-1})&{}\exp (\lambda _{2} x_{i-1})\\ Z_{i}&{}\exp (\lambda _{1} x_{i})&{}\exp (\lambda _{2} x_{i})\\ Z_{i+1}&{}\exp (\lambda _{1} x_{i+1})&{}\exp (\lambda _{2} x_{i+1}) \end{array}\right| =0. \end{aligned}$$
(12)

Simplifying the determinant in (12), we get

$$\begin{aligned} -\exp (\frac{\delta h}{2\varepsilon })Z_{i+1}+2\cosh \bigg (\frac{h\sqrt{\delta ^{2}+4\varepsilon \eta }}{2\varepsilon }\bigg )Z_{i}-\exp (\frac{-\delta h}{2\varepsilon })Z_{i-1}=0, \end{aligned}$$
(13)

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

$$\begin{aligned} Z_{i}=C_{1}\exp (\lambda _{1} x_{i})+C_{2}\exp (\lambda _{2} x_{i}) \end{aligned}$$
(14)

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

$$\begin{aligned} -\exp (\frac{\delta h}{2\varepsilon })Z_{i+1}+\bigg (\exp (\frac{\delta h}{2\varepsilon })+\exp (\frac{-\delta h}{2\varepsilon })\bigg )Z_{i}-\exp (\frac{-\delta h}{2\varepsilon })Z_{i-1}=0, \end{aligned}$$
(15)

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

$$\begin{aligned} \exp (\frac{\delta h}{\varepsilon })Z_{i+1}-\bigg (\exp (\frac{\delta h}{\varepsilon })+1\bigg )Z_{i}+Z_{i-1}=0, \end{aligned}$$

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

$${Z_{i + 1}} - 2{Z_i} + {Z_{i - 1}} + \left( {\exp (\frac{{\delta h}}{\varepsilon }) - 1} \right)\left( {{Z_{i + 1}} - {Z_i}} \right) = 0.$$
(16)

From (16), we have

$$\begin{aligned} \varepsilon \frac{ Z_{i-1}-2Z_i+Z_{i+1}}{\gamma ^{2}}+\delta \frac{Z_{i+1}-Z_i}{h}=0, \end{aligned}$$
(17)

where the denominator function is given by

$$\begin{aligned} \gamma ^{2}(h,\varepsilon ,\mu )\equiv \gamma ^{2}=\frac{h\varepsilon }{\mu a}\bigg (\exp (\frac{\mu a h}{\varepsilon })-1\bigg ). \end{aligned}$$
(18)

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

$$\begin{aligned} \gamma _i^{2}(h,\varepsilon ,\mu )\equiv \gamma _i^{2}=\frac{h\varepsilon }{\mu a_i}\bigg (\exp (\frac{\mu a_i h}{\varepsilon })-1\bigg ). \end{aligned}$$
(19)

Thus, we get the following fully discrete problem

$$\begin{aligned} \begin{array}{l} \varepsilon \frac{ Z_{i-1}^{j+1}-2Z_i^{j+1}+Z_{i+1}^{j+1}}{\gamma _i^2}+\mu a_i^{j+1}\bigg (\frac{Z_{i+1}^{j+1}-Z_i^{j+1}}{h}\bigg )-b_i^{j+1}Z_i^{j+1}-\frac{ Z_i^{j+1}-Z_i^{j}}{\Delta t}\\ =-c_i^{j+1}Z_i^{j+1-m}+f_i^{j+1}, \end{array} \end{aligned}$$
(20)

with the following discrete initial and boundary conditions

$$\begin{aligned} {\left\{ \begin{array}{ll} Z_i^{-j}=\theta (s_i,-t^{j}), &{} i=1,\cdots ,N-1, \;\; j=0,\cdots , m,\\ Z_0^{j+1}=q_0(t_{j+1}), &{} t_{j+1}\in [0,T],\\ Z_N^{j+1}=q_1(t_{j+1}), &{} t_{j+1}\in [0,T], \end{array}\right. } \end{aligned}$$
(21)

where the denominator function is given as

$$\begin{aligned} \gamma _i^{2}(h,\varepsilon ,\mu )\equiv \gamma _i^{2}=\frac{h\varepsilon }{\mu a_i^{j+1}}\bigg (\exp (\frac{\mu a_i^{j+1} h}{\varepsilon })-1\bigg ). \end{aligned}$$

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

$$\begin{aligned} WZ=G, \;\;\; i=1,2,...,N-1, \;\;\; j=0,...,M, \end{aligned}$$
(22)

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

$$\begin{aligned} {\left\{ \begin{array}{ll} W_{i},_{i-1}=\frac{\varepsilon }{\gamma ^2_{i}}, &{} i=1,...,N-2,\\ W_{i},_{i}=-(\frac{2\varepsilon }{\gamma ^2_i}+\frac{\mu a_i^{j+1}}{h}+\frac{1}{\Delta t}+b_i^{j+1}), &{} i=1,...,N-1,\\ W_{i},_{i+1}=\frac{\varepsilon }{\gamma ^2_{i}}+\frac{\mu a_i^{j+1}}{h}, &{} i=1,...,N-1. \end{array}\right. } \end{aligned}$$
(23)

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

$$\begin{aligned} {\left\{ \begin{array}{ll} G_0^{j+1}=q_0(t_{j+1}),\\ G_i^{j+1}=-c_i^{j+1}Z_i^{j+1-m}+f_i^{j+1}-\frac{Z_i^{j}}{\Delta t}, &{} i=1(1)N-1,\\ G_N^{j+1}=q_1(t_{j+1}),\\ Z=[Z_0,Z_1,\cdots ,Z_N]^{T}. \end{array}\right. } \end{aligned}$$
(24)

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

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}\Pi _{s}^p&=\frac{\varepsilon }{\gamma ^2_{s}}(\Pi _{s+1}^{p+1}-2\Pi _{s}^{p+1}+\Pi _{s-1}^{p+1})+\frac{\mu a_s^{p+1}}{h}(\Pi _{s+1}^{p+1}-\Pi _s^{p+1})-d_s^{p+1}\Pi _s^{p+1},\\&=\frac{\varepsilon }{\gamma ^2_s}[(\Pi _{s+1}^{p+1}-\Pi _s^{p+1})-(\Pi _s^{p+1}-\Pi _{s-1}^{p+1})]+\frac{\mu a_s^{p+1}}{h}(\Pi _{s+1}^{p+1}-\Pi _s^{p+1})-d_s^{p+1}\Pi _s^{p+1},\\&\ge 0, \end{aligned} \end{aligned}$$

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

$$\begin{aligned} \vert P_i^{j+1} \vert \le \beta ^{-1} \max \limits _{1\le i \le N-1}\vert {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}P_i^{j+1}\vert , \;\;\; \text {for}\;\; 0<i<N. \end{aligned}$$

Proof

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

$$\begin{aligned} =\beta ^{-1} \max \limits _{1\le i \le N-1}\vert {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}P_i^{j+1}\vert . \end{aligned}$$

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

$$\begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(\Psi ^{\pm })_i^{j+1}\equiv -\frac{(b_i^{j+1}+1/\Delta t)}{\beta }\max \limits _{1\le i \le N-1}\vert {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}P_i^{j+1}\vert \pm {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}P_i^{j+1}. \end{aligned}$$

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

$$\begin{aligned} \lim \limits _{\varepsilon \rightarrow \ 0} \max \limits _{1<i<N-1}\frac{\exp (\frac{-Cx_{i}}{\sqrt{\varepsilon }})}{\varepsilon ^{\frac{k}{2}}}=0 \;\;\; \text{ and } \;\;\; \lim \limits _{\varepsilon \rightarrow \ 0} \max \limits _{1<i<N-1}\frac{\exp (\frac{-C(1-x_{i})}{\sqrt{\varepsilon }})}{\varepsilon ^{\frac{k}{2}}}=0, \end{aligned}$$

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

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(Z_i^{j+1}-z_i^{j+1})&=({\mathcal {L}}_{\varepsilon ,\mu } -{\mathcal {L}}_{\varepsilon ,\mu }^{N,M})z_i^{j+1}, \\&={\mathcal {L}}_{\varepsilon ,\mu }z_i^{j+1}-{\mathcal {L}}_{\varepsilon ,\mu }^{N,M}z_i^{j+1},\\&=\varepsilon \left( z^{\prime \prime }\right) _i^{j+1}+\mu a_i^{j+1} \left( z^{\prime }\right) _i^{j+1}-b_i^{j+1}z_i^{j+1}\\&-\bigg [\varepsilon \bigg (\frac{z_{i+1}^{j+1}-2z_i^{j+1}+z_{i-1}^{j+1}}{\gamma _i^2}\bigg )+\mu a_i^{j+1}\bigg (\frac{z_{i+1}^{j+1}-z_i^{j+1}}{h}\bigg )-b_i^{j+1}z_i^{j+1}\bigg ]. \end{aligned} \end{aligned}$$

Simplifying the above expression, we obtain

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(Z_i^{j+1}-z_i^{j+1})&= \varepsilon \left( z^{\prime \prime }\right) _i^{j+1}+\mu a_i^{j+1}\left( z^{\prime }\right) _i^{j+1}\\ {}&-\varepsilon \bigg (\frac{z_{i+1}^{j+1}-2z_i^{j+1}+z_{i-1}^{j+1}}{\gamma _i^2}\bigg )-\mu a_i^{j+1}\bigg (\frac{z_{i+1}^{j+1}-z_i^{j+1}}{h}\bigg ). \end{aligned} \end{aligned}$$
(25)

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

$$\begin{aligned} \begin{aligned} z_{i+1}^{j+1}&=z_i(t)+h\left( z^{\prime }\right) _i^{j+1}+\frac{h^2}{2}\left( z^{\prime \prime }\right) _i^{j+1}+\frac{h^3}{6}\left( z^{\prime \prime \prime }\right) _i^{j+1}+\frac{h^4}{24}\left( z^{\prime \prime \prime \prime }\right) _i^{j+1}+\cdots \\ z_{i-1}^{j+1}&=z_i^{j+1}-h\left( z^{\prime }\right) _i^{j+1}+\frac{h^2}{2}\left( z^{\prime \prime }\right) _i^{j+1}-\frac{h^3}{6}\left( z^{\prime \prime \prime }\right) _i^{j+1}+\frac{h^4}{24}\left( z^{\prime \prime \prime \prime }\right) _i^{j+1}+\cdots . \end{aligned} \end{aligned}$$

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

$$\begin{aligned} z_{i+1}^{j+1}-2z_i^{j+1}+z_{i-1}^{j+1}=h^2\left( z^{\prime \prime }\right) _i^{j+1}+\frac{h^4}{12}\left( z^{\prime \prime \prime \prime }\right) _i^{j+1}+\cdots . \end{aligned}$$
(26)

From the first Taylor series expansion, we have

$$\begin{aligned} z_{i+1}^{j+1}-z_i^{j+1}=h\left( z^{\prime }\right) _i^{j+1}+\frac{h^2}{2}\left( z^{\prime \prime }\right) _i^{j+1}+\frac{h^3}{6}\left( z^{\prime \prime \prime }\right) _i^{j+1}+\cdots . \end{aligned}$$
(27)

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

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(Z_i^{j+1}-z_i^{j+1})&= \varepsilon \left( z^{\prime \prime }\right) _i^{j+1}+\mu a_i^{j+1}\left( z^{\prime }\right) _i^{j+1}- \frac{\varepsilon }{\gamma _i^2} \bigg (h^2\left( z^{\prime \prime }\right) _i^{j+1}+\frac{h^4}{12}\left( z^{\prime \prime \prime \prime }\right) _i^{j+1}\bigg )\\ {}&-\frac{\mu a_i^{j+1}}{h}\bigg (h\left( z^{\prime }\right) _i^{j+1}+\frac{h^2}{2}\left( z^{\prime \prime }\right) _i^{j+1}+\frac{h^3}{6}\left( z^{\prime \prime \prime }\right) _i^{j+1}\bigg ). \end{aligned} \end{aligned}$$
(28)

Simplifying (28), we obtain

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(Z_i^{j+1}-z_i^{j+1})&= \varepsilon \left( z^{\prime \prime }\right) _i^{j+1} -\frac{\varepsilon }{\gamma _i^2} \bigg (h^2\left( z^{\prime \prime }\right) _i^{j+1}+\frac{h^4}{12}\left( z^{\prime \prime \prime \prime }\right) _i^{j+1}\bigg )-\bigg (\frac{\mu a_i(t)}{2}\left( z^{\prime \prime }\right) _i^{j+1}\bigg ) h\\bigg(\frac{\mu a_i^{j+1}}{6}\left( z^{\prime \prime \prime }\right) _i^{j+1}\bigg )h^2. \end{aligned} \end{aligned}$$
(29)

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

$$\begin{aligned} \frac{1}{\gamma _i^2}=\frac{1}{h^2}-\frac{\mu a_i^{j+1}}{2\varepsilon h}+\frac{\mu ^{2} \left( a_i^{j+1}\right) ^{2}}{12\varepsilon ^{2} } \end{aligned}$$
(30)

in (30) gives

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(Z_i^{j+1}-z_i^{j+1})&=\bigg (\mu a_i^{j+1}\left( z^{\prime \prime }\right) _i^{j+1}\bigg )h-\bigg (\frac{(\mu a_i^{j+1})^2}{12\varepsilon }\left( z^{\prime \prime }\right) _i^{j+1}+\frac{\mu a_i(^{j+1}}{6}\left( z^{\prime \prime \prime \prime }\right) _i^{j+1}\bigg )h^2. \end{aligned} \end{aligned}$$
(31)

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

$$\begin{aligned} \vert {\mathcal {L}}_{\varepsilon ,\mu }^{N,M}(Z_i^{j+1}-z_i^{j+1})\vert \le Ch. \end{aligned}$$

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

$$\begin{aligned} \vert Z_i^{j+1}-z_i^{j+1}\vert \le Ch. \end{aligned}$$

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

$$\begin{aligned} \max \limits _{0\le i\le N; 0\le j\le M} \vert Z_i^{j+1}-z(s_i,t_{j+1})\vert \le C(h+\Delta t). \end{aligned}$$
(32)

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

$$\begin{aligned} \vert Z_i^{j+1}-z(s_i,t_{j+1})\vert \le C(h+\Delta t), \end{aligned}$$
(33)

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}\)

$$\begin{aligned} Z_i^{j+1}-z(s_i,t_{j+1}) \le C(h+\Delta t)+R^{N,M}, \end{aligned}$$
(34)

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

$$\begin{aligned} \tilde{Z}_i^{j+1}-z(s_i,t_{j+1}) \le C(\frac{h}{2}+\frac{\Delta t}{2})+R^{2N,2M}, \end{aligned}$$
(35)

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

$$\begin{aligned} (Z_i^{j+1})^{extr}=2\tilde{Z}_i^{j+1}-Z_i^{j+1}, \end{aligned}$$
(36)

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

$$\begin{aligned} \sup \limits _{0<\varepsilon \le 1} \ \max \limits _{0\le i\le N; 0\le j\le M} \vert (Z_i^{j+1})^{extr}-z(s_i,t_{j+1})\vert \le C((\Delta t)^2+h^2). \end{aligned}$$

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

$$\begin{aligned} e_{\varepsilon ,\mu }^{N,M}=\max \limits _{0\le i\le N;t\in [0,T]}\big \vert Z^{N,M}(s_i,t_j)-Z^{2N,2M}(s_i,t_j)\big \vert , \end{aligned}$$

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

$$\begin{aligned} (e_{\varepsilon ,\mu }^{N,M})^{extr}=\max \limits _{0\le i\le N;t\in [0,T]}\big \vert (Z^{N,M})^{extr}(s_i,t_j)-(Z^{2N,2M})^{extr}(s_i,t_j)\big \vert , \end{aligned}$$

where \(Z^{N,M}(s_i,t_j)\) is the numerical solution with (NM) mesh points and \(Z^{2N,2M}(s_i,t_j)\) is the numerical solution at the finer mesh with (2N, 2M) mesh points before extrapolation (BE). The numerical solutions after extrapolation (AE) are \((Z^{N,M})^{extr}(s_i,t_j)\) using the mesh points (NM) with mesh sizes h and \(\Delta t\) and \((Z^{2N,2M})^{extr}(s_i,t_j)\) using the mesh points (2N, 2M) 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

$$\begin{aligned} e^{N,M}=\max \limits _{\varepsilon ,\mu }e_{\varepsilon ,\mu }^{N,M}\;\;\; \text{ and } \;\;\; (e^{N,M})^{extr}=\max \limits _{\varepsilon ,\mu }(e_{\varepsilon ,\mu }^{N,M})^{extr}. \end{aligned}$$

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

$$\begin{aligned} \rho _{\varepsilon ,\mu }^{N,M}=\log _{2}\bigg (\frac{e_{\varepsilon ,\mu }^{N,M}}{e_{\varepsilon ,\mu }^{2N, 2M}}\bigg ) \;\;\; \text{ and } \;\;\; (\rho _{\varepsilon ,\mu }^{N,M})^{extr}=\log _{2}\bigg (\frac{(e_{\varepsilon ,\mu }^{N,M})^{extr}}{(e_{\varepsilon ,\mu }^{2N, 2M})^{extr}}\bigg ). \end{aligned}$$

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

$$\begin{aligned} \rho ^{N,M}=\max \limits _{\varepsilon ,\mu }\rho _{\varepsilon ,\mu }^{N,M}\;\;\; \text{ and } \;\;\; \rho _{extr}^{N,M}=\max \limits _{\varepsilon ,\mu }(\rho _{\varepsilon ,\mu }^{N,M})^{extr}. \end{aligned}$$

Example 1

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

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} &{}\varepsilon \frac{\partial ^{2} z}{\partial s^{2}}+\mu (1+s)\frac{\partial z}{\partial s}-z(s,t)-\frac{\partial z}{\partial t} =-z(s,t-1)+16s^2(1-s^2), \hspace{0.5cm} (s,t)\in (0,1)\times (0,2],\\ {} &{}z(s,t)=0, \hspace{2cm} (s,t)\in (0,1)\times (-1,0],\\ &{}z(0,t)=0, \hspace{2cm} z(1,t)=0, \hspace{2cm} t\in (0,2]. \end{aligned} \end{array}\right. } \end{aligned}$$

Example 2

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

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} &{}\varepsilon \frac{\partial ^{2} z}{\partial s^{2}}+\mu (1+s(1-s)+t^2)\frac{\partial z}{\partial s}-(1+5st)z(s,t)-\frac{\partial z}{\partial t}=-z(s,t-1)+s(1-s)(e^t-1),\\ {} &{} (s,t)\in (0,1)\times (0,2],\;\;\;\\ &{}z(s,t)=0, \hspace{2cm} (s,t)\in (0,1)\times (-1,0],\\ {} &{}z(0,t)=0, \hspace{2cm} z(1,t)=0, \hspace{2cm} t\in (0,2]. \end{aligned} \end{array}\right. } \end{aligned}$$
Table 1 Comparison of \(e^{N,M}_{\varepsilon ,\mu }, \; (e_{\varepsilon ,\mu }^{N,M})^{extr}\) and \(\rho ^{N,M}_{\varepsilon ,\mu }\) for fixed \(\mu =10^{-4}\) and varying \(\varepsilon \) for Example (1) with \(e^{N,M}_{\varepsilon ,\mu }\) in [27] using Shishkin (S-) mesh and Bakhvalov-Shishkin (BS-) mesh
Table 2 Comparison of \((e^{N,M}),\; (e^{N,M})^{extr}, \; \rho ^{N,M} \), and \((\rho ^{N,M})^{extr}\) for Example (1) with [26] using \(\mu =10^{-3}\)
Table 3 Comparison of \((e^{N,M}),\; (e^{N,M})^{extr}, \; \rho ^{N,M} \), and \((\rho ^{N,M})^{extr}\) for Example (1) with [26] using \(\mu =10^{-9}\)
Table 4 \(e^{N,M}_{\varepsilon ,\mu }, \; (e_{\varepsilon ,\mu }^{N,M})^{extr}\), \(\rho ^{N,M}_{\varepsilon ,\mu }\) and \(\left( \rho ^{N,M}_{\varepsilon ,\mu }\right) ^{extr} \) for \(\mu =10^{-6}\) and varying \(\varepsilon \) for Example (2)
Table 5 Comparison of \((e^{N,M}),\; (e^{N,M})^{extr}, \; \rho ^{N,M} \), and \((\rho ^{N,M})^{extr}\) for Example (2) with [26] using using \(\mu =10^{-3}\)
Table 6 Comparison of \((e^{N,M}),\; (e^{N,M})^{extr}, \; \rho ^{N,M} \), and \((\rho ^{N,M})^{extr}\) for Example (2) with [26] using \(\mu =10^{-9}\)
Table 7 Comparison of \(e^{N,M}_{\varepsilon ,\mu }, \; (e_{\varepsilon ,\mu }^{N,M})^{extr}\), \(\rho ^{N,M}_{\varepsilon ,\mu }\), \((e^{N,M})^{extr}\) and \((\rho ^{N,M})^{extr}\) for \(\epsilon =10^{-4}\) and varying \(\mu \) for Example (1) with [27]
Table 8 Comparison of \(e^{N,M}_{\varepsilon ,\mu }, \; (e_{\varepsilon ,\mu }^{N,M})^{extr}\) and \(\rho ^{N,M}_{\varepsilon ,\mu }\) for \(\epsilon =10^{-4}\) and varying \(\mu \) for Example (2)
Fig. 1
figure 1

3-D view of the numerical solution profiles for Example 1 at \(N=64=M\)

Fig. 2
figure 2

3-D view of the numerical solution profiles for Example (2) at \(N=64=M\)

Fig. 3
figure 3

Loglog plot of the maximum point-wise errors at \(\mu =10^{-4}\) for Example (1)

Fig. 4
figure 4

Loglog plot of the maximum point-wise errors at \(\mu =10^{-6}\) for Example (2)

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

  1. 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.

    Article  Google Scholar 

  2. O’Malley RE Jr. Two-parameter singular perturbation problems for second-order equations. J Math Mech. 1967;16(10):1143–64.

    Google Scholar 

  3. 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.

    Google Scholar 

  4. 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.

    Article  Google Scholar 

  5. 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.

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    Article  Google Scholar 

  8. 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.

    Google Scholar 

  9. 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.

    Article  Google Scholar 

  10. 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.

    Article  Google Scholar 

  11. 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.

    Article  Google Scholar 

  12. 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.

    Article  Google Scholar 

  13. 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.

    Article  Google Scholar 

  14. 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.

    Article  Google Scholar 

  15. 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.

    Article  Google Scholar 

  16. 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.

    Article  Google Scholar 

  17. 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.

    Article  Google Scholar 

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

    Google Scholar 

  19. 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.

    Article  Google Scholar 

  20. 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.

    Article  CAS  Google Scholar 

  21. 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.

    Article  Google Scholar 

  22. 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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  24. 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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  26. 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).

  27. 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.

    Article  Google Scholar 

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

    Google Scholar 

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

    Book  Google Scholar 

  30. 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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  32. 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.

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    Google Scholar 

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

  36. 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.

    Article  Google Scholar 

  37. 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.

    Article  Google Scholar 

  38. 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.

    Article  Google Scholar 

  39. 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.

  40. 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.

  41. 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.

    Article  Google Scholar 

  42. 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.

    Article  Google Scholar 

  43. 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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  45. 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.

    Article  Google Scholar 

Download references

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

Authors

Corresponding author

Correspondence to Imiru Takele Daba.

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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

Download citation

  • Received:

  • Accepted:

  • Published:

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

Keywords

Mathematics Subject Classification