Parameter uniform finite difference formulation with oscillation free for solving singularly perturbed delay parabolic differential equation via exponential spline

Abstract

Objective

In this work, singularly perturbed time dependent delay parabolic convection-diffusion problem with Dirichlet boundary conditions is considered. The solution of this problem exhibits boundary layer at the right of special domain. In this layer the solution experiences steep gradients or oscillation so that traditional numerical methods may fail to provide smooth solutions. We developed oscillation free parameter uniform exponentially spline numerical method to solve the considered problem.

Results

In the temporal direction, the implicit Euler method is applied, and in the spatial direction, an exponential spline method with uniform mesh is applied. To handle the effect of perturbation parameter, an exponential fitting factor is introduced. For the developed numerical scheme, stability and uniform error estimates are examined. It is shown that the scheme is uniformly convergent of linear order in the maximum norm. Numerical examples are provided to illustrate the theoretical findings.

Introduction

Delay differential equations (DDEs) are a class of differential equation where the unknown function or its derivative at a certain time depends on the solution and possibly its derivatives at earlier times. The delay in these equations represents the transport delay, incubation period, gestation time, etc. If a small positive parameter $\varepsilon$ multiply the highest derivative term of DDEs and involves at least a delay term, the DDEs are said to be singularly perturbed delay differential equations (SPDDEs). When the delay parameter magnitude is larger than the perturbation parameter, the equations are said to SPDDEs with large delay, otherwise they are said to be SPDDEs with small delay. These problems can be found various applications in science and engineering such as control systems [1], chemical reactions [2], epidemiology [3], optics and physiology [4], and neural networks [5].

Singularly perturbed delay parabolic convection-diffusion problems (SPDPCDPs) are a type of SPDDE. SPDPCDPs are significant in modeling various physical phenomena, particularly in systems where the current state influences future behavior, such as in fluid dynamics and heat transfer. These problems are crucial in modeling systems where there is a significant disparity between the rates of convection, diffusion, and delay effects. For instance, in chemical engineering [2], they can describe the behavior of reactive transport in porous media where the reaction rates and transport processes differ vastly, and delays in the system (due to transport lags or reaction time) affect the overall dynamics. Analyzing these problems helps in understanding how the interplay between fast and slow processes, coupled with delay effects, influences the stability and evolution of the system.

The solution for SPDDEs has a boundary layer because of the perturbation parameter $\varepsilon$. Thus, the solution shows large variation, oscillation, in small region of the domain. It is long familiar that most classical numerical methods are unable to give accurate result on uniform mesh for such problems especially as $\varepsilon$ goes to zero unless a very fine mesh is considered, which is computationally expensive. Hence there is the need for methods which are stable and uniformly convergent irrespective of the values of $\varepsilon$ and mesh size [6]. Finite difference numerical methods that exhibit uniform convergence and stable are mainly developed using fitted meshes and fitted operators. While fitted operator methods maintain a uniform mesh, fitted-mesh methods concentrate on selecting a fine mesh in the layer region(s). There are also other nonclassical finite difference numerical methods to solve singularly perturbed delay differential equations such as adaptive mesh refinement and domain decomposition methods.

SPDPCDPs are studied by various authors. Kaushik and Sharma [7] approximated SPDPCDPs using a weighted difference time discretization and central difference space discretization on a piecewise Shishkin mesh. They have shown that the method is stable and uniformly convergent with respect to $\varepsilon$. SPDPCDPs are estimated by Das and Natesan [8] using for time derivative an implicit-Euler scheme and for spatial derivatives a hybrid scheme which made up of midpoint upwind scheme and the central difference scheme. To solve SPDPCDPs, Gowrisankar and Natesan [9] used the upwind finite difference scheme for spatial derivatives and the backward-Euler scheme for time derivatives. They proved the proposed method is parameter uniform convergent of first order. Using the exponentially fitted operator finite difference method for spatial discretization and the Crank-Nicolson method for temporal discretization, Woldaregay et al. [10] solved SPDPCDPs through both approaches. They have shown that the proposed scheme converges uniformly with first order of convergence. Negero and Duressa [11] studied second order convergent scheme to approximate SPDPCDPs. After a year the authors [12] constructed a second order accurate scheme for solving SPDPCDPs. Negero and Duressa [13] estimated the solution of SPDPCDPs using Crank-Nicolson’s time discretization scheme and exponentially fitted cubic spline scheme for spatial discretization. Recently the following authors have developed parameter uniform convergent numerical scheme to solve SPDPCDPs. Hassen and Duressa [14] approximated SPDPCDPs using Crank-Nicolson time discretization and upwind finite difference for spatial derivative using Peano kernel theorem convergent analysis. Fitted computational method is developed by Tesfaye et al. [15] to solve SPDPCDPs. After a year the authors [16] solved SPDPCDPs by employing backward Euler scheme for time derivatives. They used a higher-order finite difference method to approximate the second-order derivative and non-symmetric finite difference schemes to approximate the first-order derivative terms. Hassen and Duressa [17] developed a parameter uniform convergent numerical scheme to solve SPDPCDPs by employing implicit Euler approach in the time direction and extended cubic B-spline collocation in the space direction. Kumar and Gowrisankar [18] have suggested an efficient numerical method for SPDPCDPs. The authors proved that the proposed numerical method converges uniformly with first-order up to logarithm in the spatial variable and also first-order in the temporal variable. Readers can refer different numerical scheme for solving SPDDEs in [19,20,21,22,23,24,25,26,27,28,29,30,31,32].

In this paper, our aim is to develop a parameter uniform numerical scheme for SPDPCDPs large delay version with Dirichlet boundary conditions. The proposed scheme comprises of implicit Euler in temporal direction and exponential spline scheme in spatial direction. We provided an exponentially fitting factor to manage the perturbation parameter’s effects. The novelty of the presented scheme is that, unlike Shishkin and Bakhvalov mesh types, it does not depend on a specially designed mesh and needs no prior knowledge regarding the boundary layer’s width and position. Results from the suggested scheme are more precise, consistent, and uniformly convergent.

Notation

The symbols $N_t$ and $N_z$ are denoted for the number of mesh elements, mesh parameters, in time and space direction, respectively; the symbol C is denoted for a generic positive constant which is independent of perturbation parameter and mesh parameters. The norm $\left\| .\right\|$ denotes supremum or maximum norm, i.e., $\left\| \pi (z,t)\right\| =\max _{(z,t)\in \Omega }|\pi (z,t)|$.

Continuous problem

Let $\Omega _z=(0,1), \Omega _t=(0,\mathcal {T}]$ be are spatial and temporal domain respectively, and $\Omega =\Omega _z\times \Omega _t$ for $\mathcal {T}>0$. We consider the following SPDPCDP of the form:

$$\begin{aligned} {\left\{ \begin{array}{ll} f_t(z,t)+\mathscr {L}_\varepsilon f(z,t)=-c(z,t)f(z,t-\kappa ) +g(z,t),\quad (z,t)\in \Omega ,\\ f(z,t) = B_b(z,t), \qquad (z,t)\in \varphi _b=[0,1]\times [-\kappa ,0],\\ f(0,t) = B_l(t),\qquad t\in \varphi _l=\{(0,t):0\le t\le \mathcal {T}\},\\ f(1,t) = B_r(t),\qquad t\in \varphi _r=\{(1,t):0\le t\le \mathcal {T}\}, \end{array}\right. } \end{aligned}$$

(1)

where $\mathscr {L}_\varepsilon f(z,t)=-\varepsilon f_{zz}(z,t)+a(z)f_z(z,t)+b(z,t)f(z,t).$ The delay parameter in the given problem is denoted by $\kappa > 0$, and the perturbation parameter by $\varepsilon \in (0,1].$ The functions a(z), b(z, t), c(z, t), and g(z, t) on $\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}=[0,1]\times [0,\mathcal {T}]$ and $B_b(z,t)$, $B_l(t)$, and $B_r(t)$ on $\varphi = \varphi _l \cup \varphi _b \cup \varphi _r$ are assumed sufficiently smooth and bounded which satisfy $a(z)\ge \gamma > 0$, $b(z,t)\ge \eta$ and $c(z,t)\ge \beta$. Assume $\mathcal {T}$ satisfy $\mathcal {T}=\varsigma \kappa$, $\varsigma$ is positive integer. Under these circumstances the problem exhibits a boundary layer at the right side of the spatial domain.

The Hölder continuous of the data, together with the compatibility condition at the corner points delay term [33] as stated below, can ensure the existence and uniqueness of the solution to the problem (1).

$$\begin{aligned} B_l(0)=B_b(0,0),\qquad B_r(0) =B_b(1,0), \end{aligned}$$

(2)

$$\begin{aligned} \frac{dB_l(0)}{dt}-\varepsilon \frac{\partial ^2B_b(0,0)}{\partial z^2}+a(0)\frac{\partial B_b(0,0)}{\partial z}+b(0,0)B_b(0,0)&=-c(0,0)B_b(0,-\kappa ) + g(0,0),\nonumber \\ \frac{dB_l(0)}{dt}-\varepsilon \frac{\partial ^2B_b(1,0)}{\partial z^2}+a(1)\frac{\partial B_b(1,0)}{\partial z}+b(1,0)B_b(1,0)&=-c(1,0)B_b(1,-\kappa ) + g(1,0). \end{aligned}$$

(3)

Let $\mathscr {L}f(z,t)=f_t(z,t)+\mathscr {L}_\varepsilon f(z,t)$, then the differential operator $\mathscr {L}$ satisfies the next Lemma.

Lemma 1

(Continuous Maximum Principle) Let $y(z,t)\in C^2(\Omega )\cap C^0(\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt})$,satisfies $\mathscr {L}y(z,t)\ge 0$ for all $(z,t)\in \Omega$ and $y(z,t)\ge 0$ for all $(z,t)\in \varphi$. Then $y(z,t)\ge 0,$ for all $(z,t)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}$.

Proof

Let $(z^*,t^*)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}$, such that $y(z^*,t^*)=\min \limits _{(z,t)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}}y(z,t)$ and suppose that $y(z^*,t^*)<0$. Obviously $(z^*,t^*)\notin \varphi$ and $(z^*,t^*)\in \Omega$. From calculus property, we have $y_z(z^*,t^*)=0$, $y_t(z^*,t^*)=0$, and $y_{zz}(z^*,t^*)>0$. Hence from (1), we have

$$\begin{aligned} \mathscr {L}y(z^*,t^*)=y_t(z^*,t^*)-\varepsilon y_{zz}(z^*,t^*)+a(z^*)y_z(z^*,t^*)+b(z^*,t^*)y(z^*,t^*)< 0. \end{aligned}$$

This contradicts the hypothesis $\mathscr {L}y(z,t)\ge 0$. Therefore $y(z,t)\ge 0$ for all $(z,t)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}$. $\square$

Lemma 2

The solution y(z, t) of the problem (1) satisfies

$$\begin{aligned} \left| f(z,t)-B_b(z,0)\right| \le Ct,\quad (z,t)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}, \end{aligned}$$

where a constant $C>0$, does not depend on $\varepsilon$.

Proof

Refer [8]. $\square$

Lemma 3

The solution f(z, t) of the problem (1) satisfies $|f(z,t)|\le C,(z,t)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}$, where a constant $C>0$, does not depend on $\varepsilon$.

Proof

From Lemma 2, we have

$$\begin{aligned} |f(z,t)|&\le |f(z,t)-B_b(z,0)|+|B_b(z,0)|\\&\le Ct + |B_b(z,0)|\\&\le C,\quad \text {since}\; t\in (0,\mathcal {T}]\;\text {and}\; B_b(z,0)\in C^2(\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}). \end{aligned}$$

This completes the proof. $\square$

The stability of the continuous differential operator $\mathscr {L}$ and an $\varepsilon$-uniform bound for the problem (1) in the maximum norm are provided by the following Lemma. The Lemma result follows from the maximum principle.

Lemma 4

(Stability result for Continuous Problem) The solution f(z, t) of (1) satisfies

$$\begin{aligned} |f(z,t)|\le \eta ^{-1}\left\| \mathscr {L}f(z,t)\right\| +\max \{|B_l(t),B_r(t),B_b(z,t)|\}. \end{aligned}$$

Proof

Let $K=\max \{|B_l(t),B_r(t),B_b(z,t)|\}$. For the barrier function $\Lambda ^\pm (z,t)=\eta ^{-1}\left\| \mathscr {L}f(z,t)\right\| +K\pm f(z,t)$, we have

$$\begin{aligned} \Lambda ^\pm (0,t)&=\eta ^{-1}\left\| \mathscr {L}f(z,t)\right\| +\max \{|B_l(t),B_r(t),B_b(0,t)|\}\pm f(0,t)>0,\\ \Lambda ^\pm (1,t)&=\eta ^{-1}\left\| \mathscr {L}f(z,t)\right\| +\max \{|B_l(t),B_r(t),B_b(1,t)|\}\pm f(1,t) >0,\\ \mathscr {L}\Lambda ^\pm (z,t)&=b(z,t)[ \eta ^{-1}\left\| \mathscr {L}f(z,t)\right\| +K] \pm \mathscr {L}f(z,t)\\&\ge \left\| \mathscr {L}f(z,t)\right\| +b(z,t)K \pm \mathscr {L}f(z,t)\\&\ge \left\| \mathscr {L}f(z,t)\right\| \pm \mathscr {L}f(z,t)\\&\ge 0. \end{aligned}$$

Hence by applying the Lemma 1, we get the required result. $\square$

Furthermore, bounds on the solution and its derivatives are provided by the following theorem.

Theorem 1

The solution f(z, t) to the problem (1) and its derivatives satisfies

$$\begin{aligned} \left| \frac{\partial ^{j+i}f(z,t)}{\partial z^j\partial t^i}\right| \le C\left( 1+\varepsilon ^{-j}\exp (-\gamma (1-z)/\varepsilon )\right) ,\qquad (z,t)\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega \hspace{-0.83328pt}}\hspace{0.83328pt}, \end{aligned}$$

where i and j are non-negative integers such that $0\le i+j\le 5$.

Proof

Refer on [34] $\square$

Numerical scheme

Time discretization

We engage a uniform mesh on the time domain $[0,\mathcal {T}]$ with time step size $\Delta t$ as $\Omega _t^{N_t}=\{t_m=m\Delta t: m=0(1)N_t,\; t_{N_t}=\mathcal {T},\; \Delta t=\mathcal {T}/N_t\}$ and $\Omega _t^P=\{t_s=-s\Delta t: s=0(1)P,\; t_P=-\kappa \}$, where $N_t$ and P are th number of mesh elements in $[0,\mathcal {T}]$ and $[-\kappa ,0]$ respectively. In order to handle the term with delay, a special mesh is selected so that it coincides with a mesh point in $\Omega _t^P$. We use the implicit Euler scheme for time derivatives, so we obtain

$$\begin{aligned} \frac{F^{m+1}(z)-F^{m}(z)}{\Delta t}-\varepsilon \frac{d^2F^{m+1}(z)}{dz^2}+a(z)\frac{dF^{m+1}(z)}{dz}+b^{m+1}(z)F^{m+1}(z)\\ =-c^{m+1}(z)F^{m+1-P}(z)+g^{m+1}(z) \end{aligned} $$

(4)

Consistently, we write (4) which gives semi-discrete scheme as

$$\begin{aligned} {\left\{ \begin{array}{ll} \left( I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}\right) F^{m+1}(z) =H^m(z),\\ F^{m+1}(0)=B_l(t_{m+1}),\quad F^{m+1}(1)=B_r(t_{m+1}),\\ F^{-s}(z,t)=B_b(z,-t_s),s=0(1)P,z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}, \end{array}\right. } \end{aligned}$$

(5)

where $H^m(z)=F^m(z)+\Delta t\left( -c^{m+1}(z)F^{m+1-P}(z)+g^{m+1}(z)\right)$, $\mathscr {L}_\varepsilon ^{*\Delta t}=-\varepsilon \frac{d^2}{dz^2}+a(z)\frac{d}{dz}+b^{m+1}(z)$, and $F^{m}(z)$ is the approximation of f(z, t) at $t=t_m=m\Delta t$, i.e., $F^m(z)\approx f(z,t_m)$.

In (5) let $\mathscr {L}_\varepsilon ^{\Delta t}=I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}$. Then $\mathscr {L}_\varepsilon ^{\Delta t}$ satisfies the next maximum principle.

Lemma 5

(Semi-discrete Maximum Principle) Let $\mu (z,t_{m+1})\in C^2(\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt})$. Assume $\mu (0,t_{m+1})\ge 0$, $\mu (1,t_{m+1})\ge 0$, and $\mathscr {L}_\varepsilon ^{\Delta t}\mu (z,t_{m+1})\ge 0$ for all $z\in \Omega _z$, then $\mu (z,t_{m+1})\ge 0$ for all $\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}$.

Proof

Let $(z^*,t_{m+1}) \in \{(z,t_{m+1}): z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}\}$, and $\min \limits _{z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}}\mu (z,t_{m+1}) = \mu (z^*,t_{m+1})<0$. Clearly, $(z^*,t_{m+1}) \notin \{(0,t_{m+1}),(1,t_{m+1})\}$. Also we have, $\frac{d\mu (z^*,t_{m+1})}{dz} =0,\; \frac{d\mu (z^*,t_{m+1})}{dt}=0$, and $\frac{d^2\mu (z^*,t_{m+1})}{dz^2} \ge 0$. Then

$$\begin{aligned} \mathscr {L}_\varepsilon ^{\Delta t}\mu (z^*,t_{m+1})&=\frac{d\mu (z^*,t_{m+1})}{dt} +\Delta t\left( -\varepsilon \frac{d^2\mu (z^*,t_{m+1})}{dz^2}+a(z^*)\frac{d\mu (z^*,t_{m+1})}{dz} \right. \\&\qquad +\left. b^{m+1}(z^*)\mu (z^*,t_{m+1})\right) >0. \end{aligned} $$

This contradicts to the assumption made. Thus $\mu (z^*,t_{m+1})\ge 0$,and hence $\mu (z,t_{m+1})\ge 0$ for all $z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}$. $\square$

The local truncation error $e_{m+1}$ for the scheme (5) is given by $e_{m+1}=f(z,t_{m+1})-\hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(z)$, where $\hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(z)$ is the solution of

$$\begin{aligned} {\left\{ \begin{array}{ll} \left( I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}\right) \hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(z) = f(z,t_m)+\Delta t\left( -c(z,t_{m+1})f(z,t_{m+1})+g(z,t_{m+1})\right) ,\\ \hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(0)=B_l(t_{m+1}),\quad \hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(1)=B_r(t_{m+1}),\\ \hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{-s}=B_b(z,-t_s),\quad s=0(1)P, z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}. \end{array}\right. } \end{aligned}$$

(6)

The local error in the time direction is estimated in the following lemma.

Lemma 6

(Local error) The local error $e_{m+1}$ at $t_{m+1}$ associated to the scheme (5) satisfies the bound $\left\| e_{m+1}\right\| \le C(\Delta t)^2$.

Proof

The function $\hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(z)$ satisfies

$$\begin{aligned} \left( I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}\right) \hspace{0.83328pt}\overline{\hspace{-0.83328pt}F\hspace{-0.83328pt}}\hspace{0.83328pt}^{m+1}(z) -\Delta t\left( -c(z,t_{m+1})f(z,t_{m+1})+g(z,t_{m+1})\right) =f(z,t_m). \end{aligned}$$

(7)

From Taylor series expansion, we get

$$\begin{aligned} f(z,t_{m})&=f(z,t_{m+1})-\Delta t f_t(z,t_{m+1})+\mathcal {O}((\Delta t)^2),\\&=f(z,t_{m+1})-\Delta t[-\mathscr {L}_\varepsilon ^{*\Delta t}f(z,t_{m+1})\\&\quad -c(z,t_{m+1})f(z,t_{m+1-P}) + g(z,t_{m+1})] + \mathcal {O}((\Delta t)^2)\\&=\left( I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}\right) f(z,t_{m+1})\\&\quad + c(z,t_{m+1})f(z,t_{m+1-P}) + g(z,t_{m+1}) + \mathcal {O}((\Delta t)^2). \end{aligned} $$

(8)

From (7) and (8) one can observe that $e_{m+1}$ is the solution of

$$\begin{aligned} {\left\{ \begin{array}{ll} \left( I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}\right) e_{m+1}=\mathcal {O}((\Delta t)^2),\\ e_{m+1}(0)=0=e_{m+1}(1). \end{array}\right. } \end{aligned}$$

(9)

Clearly the operator $I+\Delta t\mathscr {L}_\varepsilon ^{*\Delta t}$ satisfies semi-discrete maximum principle. Thus we obtain $\left\| e_{m+1}\right\| \le C(\Delta t)^2$. $\square$

$E_{m+1}=f(z,t_{m+1})-F^{m+1}(z)$ defines the global error in time direction by providing the error contribution at each time step.

Lemma 7

(Global Error) The global error $E_{m+1}$ at $t_{m+1}$ associated to (5) satisfies $\le \left\| E_{m+1} \right\| \le C\Delta t,\quad m=0(1)N_{t}-1$.

Proof

Using Lemma 6, we get

$$\begin{aligned} \left\| E_{m+1}\right\|&=\left\| \sum _{\iota =0}^{m+1}e_{\iota +1} \right\| \\&\le \left\| e_1\right\| + \left\| e_2\right\| + \cdots + \left\| e_{m+1}\right\| \\&\le (m+1)C(\Delta t)^2\\&=C((m+1)\Delta t)\Delta t,\quad (m+1)\Delta t\le \mathcal {T}\\&\le C\mathcal {T}\Delta t\\&\le C\Delta t. \end{aligned}$$

As a result, the temporal discretization process is first order uniform convergent. $\square$

The Lemmas 5,6,and 7 show the stability and consistency of the scheme (5). The derivative bound of the solution utilized to demonstrate the convergence of the method is provided by the following theorem.

Theorem 2

The solution $F^{m+1}(z)$ of (5) satisfies the estimate

$$\begin{aligned} \left| \frac{d^jF^{m+1}(z)}{dz^j}\right| \le C\left( 1+\varepsilon ^{-j}\exp (-\gamma (1-z)/\varepsilon )\right) ,\quad j=0(1)4,\quad z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}\end{aligned}$$

Proof

Refer [16, 35] $\square$

Space discretization

We divide $\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}=[0,1]$ in to $N_z$ equal number of subdomain with length of $h=1/N_z$ as $0=z_0,z_1,\dots ,z_{N_z}=1$ and $z_n=nh,n=0(1)N_z$. Define $\Omega _z^{N_z}=\{z_n=nh,n=0(1)N_z\}$ and $\Omega ^N=\Omega _z^{N_z}\times \Omega _t^{N_t}$. At $z_n=nh$ for $(m+1)^{th}$ time level, let $F^{m+1}_n$ be an approximation to $F^{m+1} (z_n)$ found by exponential spline function $Q_n(z)$ passing through the points $(z_n,F^{m+1}_n)$ and $(z_{n+1},F^{m+1}_{n+1})$. Omit the superscript $m+1$ for convenience, i.e.,$F^{m+1}_n=F_n$. For each $n^{th}$ segment, the exponential spline function has the following form [36]:

$$\begin{aligned} Q_n(z)=\phi _n\exp (\lambda (z-z_n)) + \chi _n\exp (-\lambda (z-z_n)) + \psi _n(z-z_n) + \vartheta _n,\quad n=1(1)N_z-1, \end{aligned}$$

(10)

where $\phi _n,\chi _n,\psi _n$,and $\vartheta _n$ are constants to be determined and $\lambda$ is a free parameter used to advance the accuracy of the scheme. Here $Q_n(z)\in C^2[\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}]$ interpolate $F_n$ at $z_n,n=0(1)N_z$ depends on $\lambda$ and reduce to cubic spline $Q_n(z)$ in $\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}$ as $\lambda \rightarrow 0$. To obtain the coefficients introduced in (10), $Q_n(z)$ should satisfy the condition of first derivative continuity at the common nodes. Define

$$\begin{aligned} Q_n(z_n)=F_n,\quad Q_n(z_{n+1})=F_{n+1},\quad Q_n^{''}(z_n)=M_n,\quad Q_n^{''}(z_{n+1})=M_{n+1} \end{aligned}$$

(11)

From (9) and (10) after some manipulation, we get

$$\begin{aligned} \phi _n&=h^2\frac{M_{n+1}-\exp (-\xi )M_n}{2\xi ^2 \sinh (\xi )},\quad \chi _n=h^2\frac{\exp (\xi )M_n-M_{n+1}}{2\xi ^2 \sinh (\xi )}\\ \psi _n&=\frac{F_{n+1}-F_n}{h}-h\frac{M_{n+1}-M_n}{\xi ^2},\quad \vartheta _n=F_n-h^2\frac{M_n}{\xi ^2}, \end{aligned} $$

(12)

where $\xi =\lambda h$ and $n=1(1)N_z-1.$ From the first derivative continuity $Q_{n-1}^{'}(z_n)=Q_n^{'}(z_n)$ for $n=1(1)N_z-1$, we obtain the following relations:

$$\begin{aligned} \frac{F_{n-1}-2F_n+F_{n+1}}{h^2}=L_1M_{n-1} + 2L_2M_n + L_1M_{n+1}, \end{aligned}$$

(13)

where $L_1=\frac{\sinh (\xi )-\xi }{\xi ^2\sinh (\xi )}$, $L_2=\frac{\xi \cosh (\xi )-\sinh (\xi )}{\xi ^2\sinh (\xi )}$, and $M_\nu =F^{''}(z_\nu ),\nu =n,n\pm 1$. Equation (4) can be written as

$$\begin{aligned} -\varepsilon F^{''}(z)=a(z)F^{'}(z) + p(z)F(z) - G(z), \end{aligned}$$

(14)

where $F(z)=F^{m+1}(z),p(z)=p^{m+1}(z)=\frac{1}{\Delta t} + b^{m+1}(z)$, $G(z)=G^m(z) = \frac{1}{\Delta t}F^m(z)-c^{m+1}(z)F^{m+1-P}(z)+g^{m+1}(z)$. Equation (14) discretized by the exponential spline by introducing a fitting factor $\sigma (\rho )$, we get

$$\begin{aligned} \varepsilon \sigma (\rho ) M_\nu =a_\nu F_\nu ^{'} + p_\nu F_\nu -G_\nu , \quad \text {for}\;\nu =n,n\pm 1, \end{aligned}$$

(15)

where $\rho =h/\varepsilon$. Placing exact solution in (15) and substitute the resulting in to (13), we obtain

$$ \begin{aligned}&-\varepsilon \sigma (\rho )\frac{F(z_{n-1})-2F(z_n) + F(z_{n+1})}{h^2} + L_1\left[ a_{n-1}F^{'}(z_{n-1}) + p_{n-1}F(z_{n-1})\right] \\&\qquad + 2L_2\left[ a_{n}F^{'}(z_{n}) + p_{n}F(z_{n})\right] + L_1\left[ a_{n+1}F^{'}(z_{n+1}) + p_{n+1}F(z_{n+1})\right] \\&=L_1G_{n-1} + 2L_2G_{n} + L_1G_{n+1} + TE(h), \quad n=1(1)N_z-1, \end{aligned}$$

(16)

where TE(h) is local truncation error [37], given as

$$\begin{aligned} TE(h)=\frac{h^4}{3}(L_2-2L_1)a_nF^{(3)}(\varrho _n) + \varepsilon \sigma (\rho )\frac{h^4}{12}(-12L_1+1)a_nF^{(4)}(\varrho _n) + \mathcal {O}(h^6). \end{aligned}$$

(17)

Clearly $TE(h)=\mathcal {O}(h^4)$ for $L_1+L_2=1/2$. If $L_2=5/12, L_1=1/12$, we have $TE(h)=\varepsilon \sigma (\rho )\frac{h^6}{240}F^{(6)}(\varrho _n),\varrho _n \in [z_{n-1},z_{n+1}]$. Using non-symmetric finite difference approximation for first derivative [38] $F(z_n)$, we have

$$\begin{aligned} F^{'}(z_n)&=\frac{F(z_{n+1})-F(z_{n-1})}{2h} + \mathcal {O}(h^2),\\ F^{'}(z_{n-1})&= \frac{-F(z_{n+1})+4F(z_n)-3F(z_{n-1})}{2h} + \mathcal {O}(h^2),\\ F^{'}(z_{n+1})&= \frac{3F(z_{n+1})-4F(z_n)+F(z_{n-1})}{2h} + \mathcal {O}(h^2). \end{aligned} $$

(18)

Finally by substituting the approximation of (18) in to the approximation of (16), we get the following full discretized scheme:

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n\equiv A_n^-F^{m+1}_{n-1} + A_n^0F^{m+1}_{n} + A_n^+F^{m+1}_{n+1}= G^{*m}_n,\quad n=1(1)N_z-1,\\ F^{m+1}_0=B_l(t_{m+1}),\quad F^{m+1}_{N_z}=B_r(t_{m+1}), \quad m=0(1)N_t-1,\\ F^{-s}_n=B_b(z_n,-t_s),\quad n=0(1)N_z,s=0(1)P, \end{array}\right. } \end{aligned}$$

(19)

where

$$\begin{aligned} A_n^-&=-\frac{\varepsilon \sigma (\rho )}{h^2} + L_1\left( -\frac{3a_{n-1}}{2h} + \frac{a_{n+1}}{2h} + p^{m+1}_{n-1}\right) - L_2\frac{a_n}{h},\\ A_n^0&= \frac{2\varepsilon \sigma (\rho )}{h^2} + L_1\left( \frac{2a_{n-1}}{h} - \frac{2a_{n+1}}{h}\right) + 2L_2p^{m+1}_n,\\ A_n^+&=-\frac{\varepsilon \sigma (\rho )}{h^2} + L_1\left( \frac{3a_{n+1}}{2h} - \frac{a_{n-1}}{2h} + p^{m+1}_{n+1}\right) + L_2\frac{a_n}{h},\\ G^{*m}_n&=L_1G^m_{n-1}+2L_2G^{m}_n+L_1G^m_{n+1},\quad G^m_{n-1}= \frac{1}{\Delta t}F^m_{n-1}-c^{m+1}_{n-1}F^{m+1-P}_{n-1}+g^{m+1}_{n-1}, \\ G^m_{n}&= \frac{1}{\Delta t}F^m_{n}-c^{m+1}_{n}F^{m+1-P}_{n}+g^{m+1}_{n},\quad G^m_{n+1}= \frac{1}{\Delta t}F^m_{n+1}-c^{m+1}_{n+1}F^{m+1-P}_{n+1}+g^{m+1}_{n+1}. \end{aligned}$$

Calculating fitting factor

The fitting factor $\sigma (\rho )$ is determined in such a way that the solution of (19) converges uniformly to the solution of (1). Multiplying (19) by h and assuming the limit as $h \rightarrow 0$, we obtain

$$\begin{aligned} \lim \limits _{h\rightarrow 0}\frac{\sigma (\rho )}{\rho }\left[ F^{m+1}_{n-1}-2F^{m+1}_{n}+F^{m+1}_{n+1}\right] =a_0(L_1+L_2)\lim \limits _{h\rightarrow 0}\left[ F^{m+1}_{n+1}-F^{m+1}_{n-1}\right] . \end{aligned}$$

(20)

From singular perturbation theory [39] concerning to the right boundary layer, we have

$$\begin{aligned} F^{m+1}(z)=F^{m+1}_0(z)+(B_r(t_{m+1})-F^{m+1}_0(1))\exp (-a(1)(1-z)/\varepsilon ) + \mathcal {O}(\varepsilon ), \end{aligned}$$

(21)

where $F^{m+1}_0(z)$ is the solution of reduced problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{F^{m+1}_0(z)-F^{m}_0(z)}{\Delta t}+a(z)\frac{dF^{m+1}_0(z)}{dz}+b^{m+1}(z)F^{m+1}_0(z) =-c^{m+1}(z)F^{m+1-P}_0(z)+g^{m+1}(z),\\ F_0^{-s}=B_b(z,-t_s),\quad s=0(1)P, z\in \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}. \end{array}\right. } \end{aligned}$$

Let $0=z_0<z_1<z_2<\cdots <z_{N_z}=1$, such that $z_n=nh,n=0(1)N_z$. Assume h is sufficiently small. Then discretization of (21) gives

$$\begin{aligned} F^{m+1}_n=F^{m+1}(nh)=F^{m+1}_0(nh)+(B_r(t_{m+1})-F^{m+1}_0(1))\exp (-a(1)(1/\varepsilon -n\rho )). \end{aligned}$$

(22)

Similarly, we have

$$\begin{aligned} F^{m+1}_{n+1}=F^{m+1}_0((n+1)h)+(B_r(t_{m+1})-F^{m+1}_0(1))\exp (-a(1)(1/\varepsilon -(n+1)\rho )),\\ F^{m+1}_{n-1}=F^{m+1}_0((n-1)h)+(B_r(t_{m+1})-F^{m+1}_0(1))\exp (-a(1)(1/\varepsilon -(n-1)\rho )). \end{aligned} $$

(23)

Substituting (22)- (23) into (20) and then simplifying, we obtain

$$\begin{aligned} \sigma (\rho )=a(1)\rho (L_1+L_2)\coth \left( \frac{a(1)\rho }{2}\right) . \end{aligned}$$

(24)

Stability and convergence analysis

The uniform stability and convergence analysis for (19) are covered in this section. Firstly, we establish the existence of the unique discrete solution for the scheme (19) by proving the discrete comparison principle.

Lemma 8

(Discrete Comparison Principle) Let $U^{m+1}_n$ and $F^{m+1}_n$ be two mesh functions, satisfying $\mathscr {L}_\varepsilon ^{\Delta t,h}U^{m+1}_n\le \mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n$, for $n=1(1)N_z-1$, $U^{m+1}_0 \le F^{m+1}_0$, and $U^{m+1}_{N_z} \le F^{m+1}_{N_z}$, then $U^{m+1}_{n} \le F^{m+1}_{n}$ for $n=0(1)N_z$.

Proof

The matrix associate with $\mathscr {L}_\varepsilon ^{\Delta t,h}$ is size of $(N_z+1)\times (N_z+1)$, where for $n=1$,and $n=N_z-1$, the terms demanding $F^{m+1}_0$ and $F^{m+1}_{N_z}$ shifted to the right side. The coefficient matrix satisfies the property of an irreducible M matrix. Hence it has positive inverse. So the existence of unique solution for (19) ensured. The reader can refer to [35] for further details. $\square$

Lemma 9

(Stability result for discrete problem) The solution $F^{m+1}_n$ of (19), satisfies the estimate

$$\begin{aligned} \left| F^{m+1}_n\right| \le \pi ^{-1}\left\| \mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n \right\| + \max \left\{ \left| F^{m+1}_0\right| ,\left| F^{m+1}_{N_z}\right| \right\} , \end{aligned}$$

where $\pi$ is the lower bound of $p^{m+1}_n$.

Proof

Let $\theta =\pi ^{-1}\left\| \mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n \right\| + \max \left\{ \left| F^{m+1}_0\right| ,\left| F^{m+1}_{N_z}\right| \right\}$. For the barrier function $\Theta ^{\pm }_n=\theta \pm F^{m+1}_n$, we have $\Theta ^{\pm }_0=\theta \pm F^{m+1}_0 \ge 0$, $\Theta ^{\pm }_{N_z}=\theta \pm F^{m+1}_{N_z} \ge 0$. Moreover for $0<n<N_z$ (19) can be written as

$$\begin{aligned} \mathscr {L}_\varepsilon ^{\Delta t,h}\Theta _n^{\pm }&=\varepsilon \sigma (\rho )\left[ (\theta \pm F^{m+1}_{n-1})-2(\theta \pm F^{m+1}_{n})+(\theta \pm F^{m+1}_{n+1})\right] \\&\qquad + L_1\left[ a_{n-1}(\theta \pm F^{m+1}_{n-1})^{'}+p^{m+1}_{n-1}(\theta \pm F^{m+1}_{n-1})\right] \\&\qquad + 2L_2\left[ a_{n}(\theta \pm F^{m+1}_{n})^{'}+p^{m+1}_{n}(\theta \pm F^{m+1}_{n})\right] \\&\qquad + L_1\left[ a_{n+1}(\theta \pm F^{m+1}_{n+1})^{'}+p^{m+1}_{n+1}(\theta \pm F^{m+1}_{n+1})\right] \\&=\left( L_1p^{m+1}_{n-1} + 2L_2p^{m+1}_{n} + L_1p^{m+1}_{n+1}\right) \theta \pm \mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n\\&\ge 0,\qquad \text {as}\;p^{m+1}_n >\pi . \end{aligned}$$

By applying discrete comparison principle, we obtained the required result. $\square$

From the power expansion of $\coth (q)$ and its property one can deduce that $|q\coth (q)-1|\le C\frac{q^2}{q^2+1}$. Thus, for $C_1$ and $C_2$ constants

$$\begin{aligned} C_1\frac{q^2}{q^2+1}\le q\coth (q)-1\le C_2\frac{q^2}{q^2+1},\; \text {and}\;\varepsilon \frac{(h/\varepsilon )^2}{(h/\varepsilon )^2+1}=\frac{h^2}{h+\varepsilon }. \end{aligned}$$

(25)

Using (25), we can have

$$\begin{aligned} \left| \varepsilon \left( a(1)\rho (L_1+L_2)\coth \left( \frac{a(1)\rho }{2}\right) -1\right) D^+D^-F^{m+1}(z_n)\right| \le C\frac{h^2}{h+\varepsilon }\left\| (F^{m+1})^{''}(z_n) \right\| . \end{aligned}$$

(26)

Then for $n=1(1)N_z-1$, using (26) we have

$$\begin{aligned} |-\varepsilon (\sigma (\rho )D^+D^-F^{m+1}_n-(F^{m+1})^{''}(z_n))|&= |-\varepsilon (\sigma (\rho )-1)D^+D^-F^{m+1}_n \\ &\quad -\varepsilon (D^+D^-F^{m+1}_n-(F^{m+1})^{''}(z_n))|\\&\le \varepsilon \left| \left( a(1)\rho (L_1+L_2)\coth \left( \frac{a(1)\rho }{2}\right) \right. \right. \\&\quad \left. \left. -1\right) D^+D^-F^{m+1}(z_n)\right| \\&\quad +\varepsilon |(F^{m+1})^{''}(z_n)-D^+D^-F^{m+1}_n|\\&\le C\frac{h^2}{h+\varepsilon }\left\| (F^{m+1})^{''}(z_n) \right\| \\&\quad + C\varepsilon h^2\left\| (F^{m+1})^{(4)}(z_n) \right\| . \end{aligned} $$

(27)

From Taylor series expansion, we have

$$\begin{aligned} |e^{'}_{n-1}|&=|(F^{m+1})^{'}(z_{n-1})-(F^{m+1})^{'}_{n-1}| \le Ch\left\| (F^{m+1})^{''}(x_j) \right\| , \end{aligned}$$

(28)

$$\begin{aligned} |e^{'}_{n}|&=|(F^{m+1})^{'}(z_{n})-(F^{m+1})^{'}_{n}| \le Ch^2\left\| (F^{m+1})^{(3)}(x_j) \right\| ,\end{aligned}$$

(29)

$$\begin{aligned} |e^{'}_{n+1}|&=|(F^{m+1})^{'}(z_{n+1})-(F^{m+1})^{'}_{n+1}| \le Ch\left\| (F^{m+1})^{''}(x_j) \right\| , \end{aligned}$$

(30)

where $\left\| (F^{m+1})^{''}(x_j) \right\| =\max _{0\le n\le N_z}|(F^{m+1})^{''}(z_n)|$,

$\left\| (F^{m+1})^{(3)}(x_j) \right\| =\max _{0\le n\le N_z}|(F^{m+1})^{(3)}(z_n)|$. The following theorem provides the spatial direction truncation error bound for the proposed scheme.

Theorem 3

(Error in the spatial direction) Let $a(z),b(z,t_{m+1})$, and$c(z,t_{m+1})$ are sufficiently smooth functions so that $F^{m+1}(z)\in C^4[\hspace{0.83328pt}\overline{\hspace{-0.83328pt}\Omega _z\hspace{-0.83328pt}}\hspace{0.83328pt}]$. Then the solution $F^{m+1}_n$ of (19) satisfies the estimate

$$\begin{aligned} \left| \mathscr {L}_\varepsilon ^{\Delta t}F^{m+1}(z_n)-\mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n\right| \le Ch\left( 1+\varepsilon ^{-3}\exp \left( -\frac{\gamma (1-z_n)}{\varepsilon }\right) \right) . \end{aligned}$$

Proof

The local truncation error bound for (19) at node $z_n$ is

$$\begin{aligned} \left| \mathscr {L}_\varepsilon ^{\Delta t}F^{m+1}(z_n)-\mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n\right|&\le \left\| TE(h) \right\| + |-\varepsilon (\sigma (\rho )D^+D^-F^{m+1}_n \\&\quad - (F^{m+1})^{''}(z_n))|+ |e^{'}_{n-1}|+ |e^{'}_{n}|+ |e^{'}_{n+1}|. \end{aligned} $$

By bounding (17) and using the bounds (27)–(30) we obtain

$$\begin{aligned} \left| \mathscr {L}_\varepsilon ^{\Delta t}F^{m+1}(z_n)-\mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n\right|&\le Ch^4(L_2-2L_1)\left\| (F^{m+1})^{(3)}(\varrho _n)\right\| \\&+ Ch^5(-12L_1+1)\left\| (F^{m+1})^{(4)} (\varrho _n)\right\| \\&+ C\frac{h^2}{h+\varepsilon } \left\| (F^{m+1})^{''}(z_n) \right\| + C\varepsilon h^2\left\| (F^{m+1})^{(4)}(z_n) \right\| \\&+ L_1Ch\left\| (F^{m+1})^{''}(x_j) \right\| +L_2Ch^2\left\| (F^{m+1})^{''}(x_j) \right\| . \end{aligned} $$

(31)

Now apply Theorem 2 on (30),

$$\begin{aligned} \left| \mathscr {L}_\varepsilon ^{\Delta t}F^{m+1}(z_n)-\mathscr {L}_\varepsilon ^{\Delta t,h}F^{m+1}_n\right|&\le Ch^4(L_2-2L_1)\left( 1+\varepsilon ^{-3}\exp (-\gamma (1-z_n)/\varepsilon )\right) \\&\quad + Ch^5(-12L_1+1)\left( 1+\varepsilon ^{-4} \exp (-\gamma (1-z_n)/\varepsilon )\right) \\&\quad + C\frac{h^2}{h+\varepsilon }\left( 1+\varepsilon ^{-2} \exp (-\gamma (1-z_n)/\varepsilon )\right) \\&\quad + C\varepsilon h^2 \left( 1+\varepsilon ^{-4} \exp (-\gamma (1-z_n)/\varepsilon )\right) \\&\quad +CL_1h\left( 1+\varepsilon ^{-2} \exp (-\gamma (1-z_n)/\varepsilon )\right) \\&\quad +CL_2h^2\left( 1+\varepsilon ^{-3} \exp (-\gamma (1-z_n)/\varepsilon )\right) \\&\le Ch \left( 1+\varepsilon ^{-3} \exp (-\gamma (1-z_n)/\varepsilon )\right) ,\\ \text {since}\;\varepsilon ^{-2}<\varepsilon ^{-3},\text {and}\;h^5<h^4<h^2<h. \end{aligned} $$

Now that the proof is completed. $\square$

Lemma 10

For a fixed mesh number $N_z$ and for $\varepsilon \rightarrow 0$, we have

$$\begin{aligned} \lim \limits _{\varepsilon \rightarrow 0}\max \limits _{1\le n\le N_z-1}\varepsilon ^{-\lambda }\exp \left( -\frac{\gamma (1-z_n)}{\varepsilon }\right) =0,\lambda \; \text {is positive integer}. \end{aligned}$$

Proof

Refer [14]. $\square$

Theorem 4

Let $F^{m+1}(z_n)$ and $F^{m+1}_n$ be the solutions of (5) and (19) respectively, then we have uniform error estimate

$$\begin{aligned} \left\| F^{m+1}(z_n)-F^{m+1}_n \right\| =\sup \limits _{0<\varepsilon \ll 1}|F^{m+1}(z_n)-F^{m+1}_n|\le Ch. \end{aligned}$$

Proof

Use the result in Lemma 10 to Theorem 3, then applying Lemma 8 gives the required result. $\square$

The uniform error bound of the scheme in the maximum norm is provided by the following theorem.

Theorem 5

Let $f(z_n,t_m)$ and $F^m_n$ are to be the solution of the continuous problem (1) and the discrete problem (19) respectively. Then we have

$$\begin{aligned} \left\| f(z_n,t_m)-F^m_n \right\| \le C(\Delta t + h),\quad n=0(1)N_z,m=0(1)N_t. \end{aligned}$$

Proof

The proof follows from the combination of Lemma 7 and Theorem 4. $\square$

Numerical examples, results, and discussion

We consider two problems of singularly perturbed parabolic differential equations with large delays in order to illustrate the applicability of the method.

Example 1

Consider the following problem [40]:

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial f}{\partial t}-\varepsilon \frac{\partial ^2f}{\partial z^2} + (2-z^2)\frac{\partial f}{\partial z} + (z+1)(t+1)f(z,t) \\ \qquad =-f(z,t-1) + 10t^2\exp (-t)z(1-z),(z,t)\in (0,1)\times (0,2],\\ f(0,1) = 0, f(1,t)=0, t\in (0,2],\\ f(z,t) = 0, (z,t)\in [0,1]\times [-1,0], \end{array}\right. } \end{aligned}$$

whose exact solution is not known.

Example 2

Now consider the following problem [41]:

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial f}{\partial t}-\varepsilon \frac{\partial ^2f}{\partial z^2} + \frac{(5-z^2)}{3}\frac{\partial f}{\partial z} + tf(z,t) \\ \qquad =-f(z,t-1) + t^3z(1-z)\sin (\pi z),(z,t)\in (0,1)\times (0,2],\\ f(0,1) = 0, f(1,t)=0, t\in (0,2],\\ f(z,t) = 0, (z,t)\in [0,1]\times [-1,0], \end{array}\right. } \end{aligned}$$

whose exact solution is not known.

Fig. 1

Numerical solution of Example 1 with boundary layer formation, $N=64$

Fig. 2

Numerical solution of Example 2 with boundary layer formation, $N=64$

Figures 1 and 2 show the physical behavior surface graph of the numerical solutions to Examples 1 and 2, respectively. From these graphs, we observe that the solution exhibits a boundary layer near $x=1$ for different values of $\varepsilon$. The graphs also display the effect of the perturbation parameter $\varepsilon$. That is, as the values of $\varepsilon$ decrease, the width of the boundary layer decreases. The log-log plot in Figs. 3 also approves the $\varepsilon$-uniform convergence of the proposed method. It is evident from these figures that as mesh numbers increases, the maximum point wise error decreases monotonically.

As the exact solutions to the problems are not known, we use the double mesh technique [42] to compute the maximum pointwise error $\tilde{E}_\varepsilon ^{N_z,N_t}$ and order of convergence $\tilde{R}_\varepsilon ^{N_z,N_t}$ as follows

$$\begin{aligned} \tilde{E}_\varepsilon ^{N_z,N_t}&=\max _{\begin{array}{c} 0\le j\le N_z \\ 0 \le i\le N_t \end{array}} \left| F^{N_z,N_t}(z_j,t_i)-F^{2N_z,2N_t}(z_j,t_i)\right| ,\\ \tilde{R}_\varepsilon ^{N_z,N_t}&= \log _2\tilde{E}_\varepsilon ^{N_z,N_t}-\log _2\tilde{E}_\varepsilon ^{2N_z,2N_t}, \end{aligned} $$

(32)

where $F^{N_z,N_t}(z_j,t_i)$ and $F^{2N_z,2N_t}(z_j,t_i)$ are the numerical approximation of the exact solution $f(z_j,t_i)$ on the mesh $\Omega ^N$ and $\Omega ^{2N}$ respectively. $\Omega ^{2N}$ is obtained by doubling the mesh $\Omega ^N$ such that the mid points $z_{j+1/2}=(z_{j+1}+z_j)/2$ and $t_{i+1/2}=(t_{i+1}+t_i)/2$ are included in to the mesh points. From (32), we compute $\varepsilon$-uniform maximum error $\tilde{E}^{N_z,N_t}$ and the corresponding $\varepsilon$-uniform convergence rate $\tilde{R}^{N_z,N_t}$ as

$$\begin{aligned} \tilde{E}^{N_z,N_t}&= \max _{\varepsilon }\tilde{E}_\varepsilon ^{N_z,N_t},\\ \tilde{R}^{N_z,N_t}&= \log _2\tilde{E}^{N_z,N_t}-\log _2\tilde{E}^{2N_z,2N_t}. \end{aligned} $$

(33)

Table 1 Maximum point wise error $\tilde{E}_\varepsilon ^{N_z,N_t}$, uniform error $\tilde{E}^{N_z,N_t}$,convergence rate $\tilde{R}^{N_z,N_t}$, and CPU time in second for Example 1, $N=M$

Table 2 Maximum point wise error $\tilde{E}_\varepsilon ^{N_z,N_t}$, uniform error $\tilde{E}^{N_z,N_t}$, convergence rate $\tilde{R}^{N_z,N_t}$, and CPU time in second for Example 2, $N=M$

The numerical results presented in Tables 1 and 2 show the maximum point-wise error, $\varepsilon$-uniform maximum error, rate of convergence, and CPU time in seconds for the proposed method for Examples 1 and 2, respectively. From these tables, we confirm that the suggested method is linear-order $\varepsilon$-uniformly convergent, according to the error analyses carried out in this work. Furthermore, we see that as the mesh numbers increase, the maximum point-wise error decreases, and as the values of $\varepsilon$ decrease, a stable and bounded maximum error is established. Thus, the proposed scheme is $\varepsilon$-uniformly convergent.

Table 3 Comparison of maximum point wise errors for Example 1

Table 4 Comparison of maximum point wise errors for Example 2

Tables 3 and 4 show the comparison of the maximum point wise error of the methods that existed in the literature [40, 41] for Examples 1 and 2, respectively. In Table 3, as the $\varepsilon$ gets smaller, we observe that the proposed method holds a more accurate $\varepsilon$ uniform convergence than the method in [40]. Similarly, in Table 4, the comparison shows that the results of the proposed scheme are more accurate $\varepsilon$ uniform convergence than the method in [41].

Fig. 3

Log-log plot of maximum point wise error 3 a for Example  1 and 3 B for Example  2

Conclusion

In this work, a non-classical numerical method is developed to solve a class of singularly perturbed delay parabolic convection-diffusion problems with Dirichlet boundary conditions. The solutions of these problems display boundary layer at the right side of spatial domain as $\varepsilon \rightarrow 0$. A delay term is handled by constructing a mesh in such a way that the delay argument coincides with a mesh point. The method is based on exponential spline on uniform mesh with exponential fitting factor. It is shown that the method is uniformly convergent independent of mesh parameters and perturbation parameter and provides uniform first-order convergence. The proposed method has the advantage of being applicable to personal computers with very low CPU computing time for the required number of mesh points. For the proposed method the evaluation of exponential functions is computationally intensive, particularly if high precision is needed. This difficulty can add to the overall CPU time required for solving the differential equations. The theoretical result is validated numerically by two test examples that are presented. The performance of the method was compared with some existing literatures and gave more accurate result.