 Research note
 Open Access
 Published:
Accurate numerical scheme for singularly perturbed parabolic delay differential equation
BMC Research Notes volume 14, Article number: 358 (2021)
Abstract
Objectives
Numerical treatment of singularly perturbed parabolic delay differential equation is considered. Solution of the equation exhibits a boundary layer, which makes it difficult for numerical computation. Accurate numerical scheme is proposed using \(\theta\)method in time discretization and nonstandard finite difference method in space discretization.
Result
Stability and uniform convergence of the proposed scheme is investigated. The scheme is uniformly convergent with linear order of convergence before Richardson extrapolation and second order convergent after Richardson extrapolation. Numerical examples are considered to validate the theoretical findings.
Introduction
Singularly perturbed delay differential equation (SPDDE) is a differential equation in which its highest order derivative term is multiplied by small perturbation parameter and involving at least one delay term. Such type of equations have variety of applications in modelling of neuronal variability [15], in control theory [3], in description of human pupillight reflex [14] and so on. Recently, a number of papers have been published on numerical treatment of time dependent singularly perturbed differential difference equations. In papers [1, 2, 4,5,6,7,8, 10,11,12,13] and [17] different authors have developed numerical scheme for treating SPDDE. The schemes in above listed papers have only linear order of convergence. In this paper, we construct second order uniformly convergent numerical scheme using nonstandard FDM with Richardson extrapolation.
Notation: In this paper, the symbols \(C, C_{1}\) and \(C_{2}\) denotes a positive constant independent of the perturbation parameter and number of mesh points. The norm \(\left\ .\right\\) denotes the maximum norm.
Considered equation
We consider a singularly perturbed parabolic delay differential equation of the form
on the domain \(D=\Omega \times \Lambda = (0,1) \times (0,T]\) for some fixed number \(T>0\) with initial and intervalboundary conditions
where \(0< \varepsilon \ll 1\) is singular perturbation parameter and \(\delta\) is delay satisfying \(\delta < \varepsilon\). The functions \(a(x), b(x),f(x,t), u_{0}(x), \phi (x,t)\) and \(\psi (1,t)\) are assumed to be sufficiently smooth and bounded with \(b(x) \ge b^{*} > 0,\) for some constant \(b^{*}\).
Solution of (1)–(2) exhibits boundary layer [4] and position of the layer depends on the conditions: If \(a(x)< 0\) left layer exist. If \(a(x)>0\) right layer exist.
Some preliminary of the analytical solution
For the case \(\delta <\varepsilon\), using Taylor’s series approximation for the terms containing delay \(u(x\delta ,t)\) and \(u_{x}(x\delta ,t)\) is valid [16]. Since, we assumed \(\delta < \varepsilon\), we approximate (1)–(2) by
with initial and boundary conditions
where \(c_{\varepsilon }(x)=\varepsilon \frac{\delta ^2}{2} b(x)+\delta a(x)\) and \(p(x)=a(x)\delta b(x)\).
For small values of \(\delta\), (1)–(2) and (3)–(4) are asymptotically equivalent. We assume \(0<c_{\varepsilon }(x) \le \varepsilon ^2\frac{\delta ^2}{2} b^{*}+\delta a^{*} =c_\varepsilon\), where \(b^{*}\) and \(a^{*}\) are the lower bound for b(x) and a(x) respectively. We assume also \(p(x) \ge p^{*} >0,\) implies occurrence of boundary layer near \(x=1\).
Lemma 2.1
([6], Theorem 2.1) The derivatives of the solution u(x, t) of (3)–(4) is bounded as
where C a positive constant independent of the parameter \(c_{\varepsilon }.\)
Numerical scheme
Temporal semidiscretization
We subdivide the time domain [0, T] into M intervals as \(t_{0}=0,\; t_{j}=j\Delta t, j=0,1,2,\ldots,M1\), where \(\Delta t= T/(M1)\). We use \(\theta\)—method for semidiscretizing (3)–(4). In general, stable numerical scheme is obtained for \(\frac{1}{2}\le \theta \le 1\). In case \(\theta =\frac{1}{2}\) it becomes Crank Nicolson method which is second order convergent. In this discretization, for each \(j=0,1,2,\ldots,M1\) we obtain a system of BVPs
where
with the boundary conditions
Lemma 3.1
(Global error estimate.) The global error estimate up to \(t_{j+1}\) time step is given by
Spatial discretization
Exact finite difference
To construct exact finite difference scheme, we follow the procedure in [9]. Consider the constant coefficient homogeneous differential equations of the form
Equation (9) has two independent solutions namely \(\exp (\lambda _{1}x)\) and \(\exp (\lambda _{2}x)\) where \(\lambda _{1,2}=\frac{p^{*}\pm \sqrt{(p^{*})^2+4c_{\varepsilon }b}}{2c_{\varepsilon }}.\)
Let \(x_{i}=x_{0}+ih,\; i=1,2,\ldots,N, \;\; x_{0}=0,\; x_{N}=1,\; h=\frac{1}{N}\) where N is the number of mesh intervals. We denote the approximate solution of u(x) at mesh point \(x_{i}\) by \(U_{i}\). Our main objective is to calculate a difference equation which has the same general solution as the differential equation in (9) has at the mesh point \(x_{i}\) given by \(U_{i}=A_{1}\exp (\lambda _{1}x_{i})+A_{2}\exp (\lambda _{2}x_{i})\). Using the theory of difference equations for second order linear difference equations in [9], we obtain
is an exact difference scheme for (9). For \(\varepsilon \rightarrow 0\), after arithmetic adjustment, we obtain
From (12) the denominator function for second derivative discretization is \(\gamma ^2 =\frac{hc_{\varepsilon }}{p^{*}}\left (\exp \left (\frac{hp^{*}}{c_\varepsilon }\right )1\right )\). For variable coefficient equation, \(\gamma ^2\) can be written as
Discrete scheme
The spatial domain \(\bar{\Omega }=[0,1]\), discretized with uniform mesh length \(\Delta x=h\) such that \(\Omega ^{N}=\{x_{i}=x_{0}+ih,\; i=1,2,\ldots,N1, \; x_{0}=0,\; x_{N}=1,\; h=\frac{1}{N}\}\) where N is the number of mesh intervals. Using the discretization in (12) into the scheme in (6), we obtain
where
and \(g_{i,j+1}=(1\theta )\Delta tL^{h,\Delta t}U_{i,j}+\Delta t[\theta f(x_{i},t_{j+1})+(1\theta )f(x_{i},t_{j})]\).
Stability analysis and uniform convergence
We need to show that the scheme in (14) satisfies the discrete maximum principle, uniform stability estimates and uniform convergence.
Lemma 3.2
(Discrete maximum principle.) Let \(U_{i,j+1}\) be any mesh function satisfying \(U_{0,j+1}\ge 0,\; U_{N,j+1}\ge 0.\) Then, \(( 1+\Delta t \theta L^{h,\Delta t})U_{i,j+1}\ge 0, i=1,2,\ldots,N1\) implies that \(U_{i,j+1}\ge 0, \forall i=0,1,\ldots,N.\)
Proof
Suppose there exist \(k\in \{0,1,\ldots,N\}\) such that \(U_{k,j+1}=\min _{0\le i\le N}U_{i,j+1}<0\), which implies \(k\ne 0,N\). Also we assume that \(U_{k+1,j+1}U_{k,j+1}>0\) and \(U_{k,j+1}U_{k1,j+1}<0\). Using the assumptions made above, we obtain \(( 1+\Delta t \theta L^{h,\Delta t})U_{k,j+1}< 0\), for \(k=1,2,3,\ldots,N1\). Thus the supposition \(U_{i,j+1}<0\), for \(i=0,1,\ldots,N\) is wrong. Hence, we obtain \(U_{i,j+1}\ge 0 ,\; \forall i=0,1,\ldots,N.\) \(\square\)
Lemma 3.3
(Uniform stability estimate.) Solution \(U_{i,j+1}\) of the discrete scheme in ( 14) satisfies the bound
Proof
Let us construct a barrier functions as \(\pi ^{\pm }_{i,j+1}=\left\ g_{i,j+1}\right\ (1+\Delta t \theta b^{*})^{1}+\max \left \{\left \phi (0,t_{j+1})\right , \left \psi (1,t_{j+1})\right \right \}\pm U_{i,j+1}\). We can easily show that \(\pi ^{\pm }_{0,j+1}\ge 0\), \(\pi ^{\pm }_{N,j+1}\ge 0\) and \((1+\Delta t \theta L^{h,\Delta t})\pi ^{\pm }_{i,j+1}\ge 0.\) By the discrete maximum principle, we obtain \(\pi ^{\pm }_{i,j+1}\ge 0, \;\forall i=0,1,2,\ldots,N\). \(\square\)
Let us define the differences operators in space as \(D^{+}V_{j+1}(x_{i})=\frac{V_{j+1}(x_{i+1})V_{j+1}(x_{i})}{h}\), \(D^{}V_{j+1}(x_{i})=\frac{V_{j+1}(x_{i})V_{j+1}(x_{i1})}{h}\) and \(D^{+}D^{}V_{j+1}(x_{i})=\frac{(D^{+}D^{})V_{j+1}(x_{i})}{h}\).
Theorem 3.1
The solution \(U_{i,j+1}\) of ( 6 ) satisfies the truncation error bound
Proof
We consider the truncation error
The estimate \(c_{\varepsilon }\left \frac{h^2}{\gamma _{i}^2}1\right \le Ch\) used in the above expression is proved in [1]. Using bound of the derivatives of the solution in Lemma 2.1 and since \(c^3_{\varepsilon }\le c^2_{\varepsilon }\), we obtain
Lemma 3.4
For a fixed mesh N, and for \(c_{\varepsilon } \rightarrow 0,\) we obtain
Using Lemma 3.4 into Theorem 3.1 gives \(\left ( 1+\Delta t \theta L^{h,\Delta t})\left ( U_{j+1}(x_{i}) U_{i,j+1}\right ) \right \le Ch\). Applying the discrete maximum principle in Lemma 3.2, we obtain the error bound as \(\left U_{j+1}(x_{i}) U_{i,j+1} \right \le Ch\).
Theorem 3.2
Solution of the scheme in ( 14 ) satisfies the uniform error bound
Proof
The uniform error bound of the scheme follows from the results of Theorem 3.1, Lemma 3.4 and the bound from temporal discretization. \(\square\)
Richardson extrapolation
We apply the Richardson extrapolation technique in spatial direction to accelerate the rate of convergence of the scheme. Let \(U_{i,j+1}^{2N,M}\) denoted for an approximate solution on 2N and M number of mesh points by including the mid points \(x_{i+1/2}\) into the mesh points, which gives that \(U^{ext}_{i,j+1}=2U_{i,j+1}^{2N,M}U_{i,j+1}\) is an extrapolated solution. The uniform error bound becomes
Numerical results and discussion
We considered two numerical examples of the form in (1)–(2) from [6, 18] to illustrate the theoretical findings of the proposed scheme.
Example 4.1
\(a(x)=2x^2,\; b(x)=x^2+1+\cos (\pi x)\) and \(f(x)=10t^{2}\exp (t)(1x)\) for \(u_{0}(x)=0,\; 0\le x\le 1\) and \(\phi (x,t)=0,\; x\in [\delta , 0],\; \psi (1,t)=0\) for final time \(T=1\).
Example 4.2
\(a(x)=2x^2,\; b(x)=3x\) and \(f(x)=\exp (t)\sin (\pi x(1x))\) for \(u_{0}(x)=0, 0\le x\le 1\) and \(\phi (x,t)=0,\;x\in [\delta , 0],\;\; \psi (1,t)=0\) for final time \(T=1\).
The exact solution of the examples are not known. We use the double mesh procedure to calculate maximum absolute error as \(E^{N,M}_{\varepsilon ,\delta }=\max _{ i,j} U^{N,M}_{i,j}U^{2N,2M}_{i,j}.\) The uniform error estimate is calculated using \(E^{N,M}=\max _{\varepsilon ,\delta }\left E^{N,M}_{\varepsilon ,\delta }\right\). The uniform rate of convergence is calculated using \(r^{N,M}=\log _{2}\left (E^{N,M}/E^{2N,2M}\right ).\)
Solution of Examples 4.1 and 4.2 exhibits a right boundary layer. As one observes in Fig. 1, as the perturbation parameter, \(\varepsilon\) goes small; the boundary layer formation becomes more visible. In Tables 1 and 2, the maximum absolute error, the uniform error and the uniform rate of convergence of the scheme before and after Richardson extrapolation is given for different values of \(\varepsilon\) and mesh numbers. As one observes the results in the tables, the maximum absolute error before Richardson extrapolation are independent of \(\varepsilon\) as, the parameter \(\varepsilon\) goes small. The scheme before Richardson extrapolation have linear order of convergence and the scheme after Richardson extrapolation have second order of convergence.
Conclusion
In this paper, second order uniformly convergent numerical scheme is developed for solving singularly perturbed parabolic delay differential equation. The developed scheme is based on non standard FDM. Stability of the scheme is investigated using construction of barrier function for the solution bound. Uniform convergence of the scheme is proved. Applicability of the scheme is investigated by considering two test examples. Effects of the perturbation parameter on the solution is shown using figures and tables. The scheme is accurate, stable and uniformly convergent.
Limitations

The proposed scheme is not layer resolving method (i.e. there is no sufficient number of mesh points in the boundary layer region).
Availability of data and materials
No additional data is used for this research work.
Abbreviations
 SPDDE:

Singularly perturbed delay differential equation
 FDM:

Finite difference method
References
Bansal K, Sharma KK. Parameterrobust numerical scheme for timedependent singularly perturbed reaction–diffusion problem with large delay. Numer Funct Anal Optim. 2018;39(2):127–54.
Chakravarthy PP, Kumar K. An adaptive mesh method for time dependent singularly perturbed differential–difference equations. Nonlinear Eng. 2019;8(1):328–39.
Glizer VY. Asymptotic analysis and solution of a finitehorizon \(H^{\infty }\) control problem for singularlyperturbed linear systems with small state delay. J Optim Theory Appl. 2003;117(2):295–325.
Gupta V, Kumar M, Kumar S. Higher order numerical approximation for time dependent singularly perturbed differential–difference convection–diffusion equations. Numer Methods Partial Differ Equ. 2018;34:357–80.
Kumar D, Kadalbajoo MK. A parameteruniform numerical method for timedependent singularly perturbed differential–difference equations. Appl Math Model. 2011;35(6):2805–19.
Kumar S, Kumar BR. A domain decomposition Taylor Galerkin finite element approximation of a parabolic singularly perturbed differential equation. Appl Math Comput. 2017;293:508–22.
Kumar S, Kumar BR. A finite element domain decomposition approximation for a semilinear parabolic singularly perturbed differential equation. Int J Nonlinear Sci Numer Simul. 2017;18(1):41–55.
Kumar R, Kumar S. Convergence of threestep Taylor Galerkin finite element scheme based monotone Schwarz iterative method for singularly perturbed differential–difference equation. Numer Funct Anal Optim. 2015;36(8):1029–45.
Mickens RE. Advances in the applications of nonstandard finite difference schemes. Singapore: World Scientific; 2005.
Parthiban S, Valarmathi S, Franklin V: A numerical method to solve singularly perturbed linear parabolic second order delay differential equation of reaction–diffusion type. Malaya J Matematik. 2015:412–20.
Rai P, Sharma KK. Singularly perturbed parabolic differential equations with turning point and retarded arguments Int J Appl Math. 2015;45(4):404–9.
Ramesh VP, Priyanga B. Higher order uniformly convergent numerical algorithm for timedependent singularly perturbed differential–difference equations. Differ Equ Dyn Syst. 2019;27:1–25.
Ramesh VP, Kadalbajoo MK. Upwind and midpoint upwind difference methods for timedependent differential difference equations with layer behavior. Appl Math Comput. 2008;202(2):453–71.
Senthilkumar LS, Mahendran R, Subburayan V. A second order convergent initial value method for singularly perturbed system of differential–difference equations of convection diffusion type. J Math Comput Sci. 2022;25(1):73–83.
Stein RB. Some models of neuronal variability. Biophys J. 1967;7(1):37–68.
Tian H. The exponential asymptotic stability of singularly perturbed delay differential equations with a bounded lag. J Math Anal Appl. 2002;270(1):143–9.
Woldaregay MM, Duressa GF. Uniformly convergent numerical method for singularly perturbed delay parabolic differential equations arising in computational neuroscience. Kragujevac J Math. 2022;46(1):65–84.
Woldaregay MM, Duressa GF. Almost secondorder uniformly convergent numerical method for singularly perturbed convection–diffusion–reaction equations with delay. Appl Anal. 2021. https://doi.org/10.1080/00036811.2021.1961756.
Acknowledgements
The authors thanks the editor and reviewers for their constructive comments.
Funding
No funding organization for this research work.
Author information
Authors and Affiliations
Contributions
Both authors contributed equally. Both authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
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
Woldaregay, M.M., Duressa, G.F. Accurate numerical scheme for singularly perturbed parabolic delay differential equation. BMC Res Notes 14, 358 (2021). https://doi.org/10.1186/s13104021057694
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13104021057694
Keywords
 Boundary layer
 Nonstandard finite difference
 Singularly perturbed
Mathematics Subject Classification
 Primary 65M06
 65M12
 Secondary 65M15