Numerical treatment of singularly perturbed parabolic partial differential equations with nonlocal boundary condition

This paper presents numerical treatments for a class of singularly perturbed parabolic partial differential equations with nonlocal boundary conditions. The problem has strong boundary layers at x = 0 and x = 1. The nonstandard finite difference method was developed to solve the considered problem in the spatial direction, and the implicit Euler method was proposed to solve the resulting system of IVPs in the temporal direction. The nonlocal boundary condition is approximated by Simpsons $\frac{1}{3}$ rule. The stability and uniform convergence analysis of the scheme are studied. The developed scheme is second-order uniformly convergent in the spatial direction and first-order in the temporal direction. Two test examples are carried out to validate the applicability of the developed numerical scheme. The obtained numerical results reflect the theoretical estimate.

1. Introduction

Differential equations that involve a small parameter in their higher order derivative term are said to be singularly perturbed problems (SPPs) or singularly perturbed differential equations (SPDEs). Many mathematical models, starting from fluid dynamics to mathematical biology, are modeled using (SPPs). For example, high Reynold's number flow in fluid dynamics, heat transport problems with large Péclet numbers, elastic vibration, etc. [1] and the references therein. Such mathematical problems are extremely difficult to solve exactly. Thus, for treating such problems numerical methods are preferable. Various scientific and engineering processes can be modeled as integral terms over the spatial domain that appear inside or outside of the boundary conditions [2, 3]. Such problems are said to be nonlocal problems. Many physical phenomena are formulated as nonlocal mathematical models. For instance, problems in thermodynamics [4], plasma physics [5], heat conduction [6], underground water flow, and populace dynamics [7] can be reduced to nonlocal problems with integral conditions. SPPs having nonlocal boundary conditions in which the highest order derivative term is multiplied by way of a small parameter are referred to as SPPs with integral boundary conditions. Such problems exhibit boundary layer phenomena wherein the solution changes. However, the numerical treatments of SPPs attract the attention of researchers due to the boundary layer behavior of the solution. Since the small parameter multiplies the highest derivative, the small regions adjoin the domain of interest's boundaries or any interior stage at which the variable quantity undergoes a very unexpected change. As a result, these problems have strong boundary layers, which ensures that there are small areas where the solution rapidly changes within very small layers near the boundary or within the problem domain [8]. Numerically treating such SPPs with nonlocal boundary conditions is a challenging task due to a very small perturbation parameter (ε).

In recent times, a class of SPPs involving nonlocal boundary conditions have been obtained great attention from scholars. To mention few of them, the authors in Bahuguna and Dabas [9], Feng et al. [10], and Li and Sun [11] studied the existence and uniqueness of a class of SPPs with nonlocal boundary conditions. The authors in Raja and Tamilselvan [12] developed a finite difference scheme for solving a class of a system of singularly perturbed reaction diffusion equations with integral boundary conditions. Debala and Duressa [13] built a uniformly convergent numerical scheme for solving SPPs with nonlocal boundary conditions. Numerical methods for solving singularly perturbed delay differential equations (SPDDEs) are considered in Sekar and Tamilselvan [14–17]. The authors developed finite difference schemes with suitable piecewise uniform Shiskin meshes. The authors in Debela and Duressa [18] used an exponentially fitted numerical scheme to solve SPDDEs of the convection-diffusion kind with nonlocal boundary conditions. Debela and Duressa [19] improved the order of accuracy for the method proposed in Debela and Duressa [18]. Kumar and Kumari [20] developed the method based on the idea of B-spline functions and an efficient numerical method on a piecewise-uniform mesh was recommended to approximate the solutions of SPPs having a delay of unit magnitude with an integral boundary condition. In the literature, only few authors considered a class of singularly perturbed parabolic partial differential equations (SPPPDEs) with integral boundary conditions. Sekar and Tamislevan [21] investigate a numerical solution for singularly perturbed delay partial differential equations (SPDPDEs) of the reaction-diffusion type with integral boundary conditions. They developed the standard finite difference on a rectangular piecewise uniform mesh for spatial discretization and a backward difference scheme in time derivative. Gobena and Duressa [22] constructed and analyzed an accurate numerical method for solving SPDPDEs with integral boundary conditions.

In general, the classical numerical methods used for solving SPDEs are not well-posed and fail to provide an exact solution when a perturbation parameter (ε) goes to zero. Therefore, it is essential to develop a numerical method that offers suitable results for small values of the perturbation parameter. As far as the researchers' knowledge, singularly perturbed parabolic partial differential equations with nonlocal boundary conditions are first being considered and have not yet been treated numerically. In this study, we investigate a uniformly convergent numerical method for solving the problem under consideration. We used the nonstandard finite difference method for space direction and the implicit Euler method for time direction.

The contents of the paper are arranged in the following manner: A brief introduction of the given problem is discussed in Section 1. In Section 2, the properties of continuous solutions are given. In Section 3, a numerical method is formulated by using the method of lines for the given problem. Stability and convergence analysis for developed numerical methods are also studied. Numerical results and discussions are given in Section 4. In Section 5, the conclusion of the paper is given.

Notation: In this paper, N and M denote the number of mesh intervals in spatial and temporal discretization, respectively. C is a generic positive constant independent of ε, N, and M. The norm used for studying the convergence of numerical solutions is the maximum norm defined as ||z(s, t)||: = sup|z(s, t)|, (s, t) ∈ D.

2. Properties of continuous problem

In this paper, we consider a class of singularly perturbed 1D parabolic partial differential equations of the reaction-diffusion type with non-local boundary conditions,

$\begin{array}{ll}\{\begin{array}{l}ℒz(s,t)=\left(-\epsilon \frac{{\partial }^2}{\partial s^2}+\frac{\partial }{\partial t}+a(s,t)\right)z(s,t)=f(s,t)(s,t)\in D,\hfill \\ z(s,0)={\varphi }_b(s),{\varphi }_b(s,t)\in {\Gamma }_b=\left\{(s,0)\right\},\hfill \\ z(0,t)={\varphi }_l(t),{\varphi }_l(s,t)\in {\Gamma }_l=\left\{(0,t);0\le t\le T\right\},\hfill \\ Kz(s,t)=z(1,t)-\epsilon {\displaystyle {\int }_0^{1}g(s)z(s,t)ds={\varphi }_r(s,t),{\varphi }_r(s,t)\in }\hfill \\ {\Gamma }_r=\left\{(1,t);0\le t\le T\right\}.\hfill \end{array}& (1)\end{array}$

where $(s,t)\in D={\Omega }_x\times {\Omega }_t=(0,1)\times (0,T],\bar{D}=\left[0,1\right]\times \left[0,T\right],$ and ε is a small parameter (0 < ε ≪ 1). Suppose that a(s, t) ≥ α > 0, f(s, t), ϕ_l, ϕ_r, ϕ_b are sufficiently smooth functions and g(s) is nonnegative monotone function and satisfies ${\displaystyle {\int }_0^{1}g(s)ds<1}$. The existence and uniqueness of the problem (1) can be established under the assumption that the data are Hölder continuous and imposing proper compatibility conditions at the corners [23]. Note that ϕ_l and ϕ_r are only functions of t, while ϕ_b is a function of x only. The problems necessarily satisfies the following sufficient compatibility conditions ϕ_b(0, 0) = ϕ_l(0, 0), ϕ_b(1, 0) = ϕ_r(1, 0), and

$\begin{array}{l}-\epsilon \frac{{\partial }^2{\varphi }_b(0,0)}{\partial s^2}+a(0,0){\varphi }_b(0,0)+\frac{\partial {\varphi }_l(0,0)}{\partial t}=f(0,0),\\ -\epsilon \frac{{\partial }^2{\varphi }_b(1,0)}{\partial s^2}+a(1,0){\varphi }_b(1,0)+\frac{\partial {\varphi }_r(1,0)}{\partial t}=f(1,0).\end{array}$

Note that ϕ_l, ϕ_r, and ϕ_b are assumed to be sufficiently smooth for Equation (1) to make sense, namely ${\varphi }_l,{\varphi }_r\in C^1(\left[0,T\right])$, and ${\varphi }_b\in C^{(2,1)}({\Gamma }_b)$.

Next, we analyze some properties of the continuous solution (Equation 1) which guarantee the existence and uniqueness of the analytical solution. A replication of this belonging in semi-discrete form can be used to present the approximate solution, which we provide in the following section.

Lemma 1. (Continuous Maximum Principle) Let $\text{Ψ}(s,t)\in C^{(0,0)}(\bar{D})\cap C^{(1,0)}(D)\cap C^{(2,1)}(D)$ be a sufficiently smooth function such that $\text{Ψ}(0,t)\ge 0,\text{Ψ}(s,0)\ge 0,{K}\text{Ψ}(1,t)\ge 0,{L}\text{Ψ}(s,t)\ge 0,\forall (s,t)\in D$. Then Ψ(s, t) ≥ 0, $\forall (s,t)\in \bar{D}$, where ${L}\text{Ψ}(s,t)={\text{Ψ}}_t(s,t)-\epsilon {\text{Ψ}}_{ss}(s,t)+a\text{Ψ}(s,t)$.

Proof. Assume (s^*, t^*) be defined as $\text{Ψ}(s^{*},t^{*})=\underset{(s,t)\in \bar{D}}{min}\text{Ψ}(s,t)$ and suppose that Ψ(s^*, t^*) ≤ 0. It is known that (s^*, t^*) ∉ ∂D. Thus, ${L}\text{Ψ}(s^{*},t^{*})={\text{Ψ}}_t(s^{*},t^{*})-\epsilon {\text{Ψ}}_{ss}(s^{*},t^{*})+a(s,t)\text{Ψ}(s^{*},t^{*})$. Since $\text{Ψ}(s^{*},t^{*})=\underset{(s,t)\in \bar{D}}{min}\text{Ψ}(s,t)$, which indicates that Ψ(s^*, t^*) = 0, ${\text{Ψ}}_t(s^{*},t^{*})=0$, ${\text{Ψ}}_{ss}(s^{*},t^{*})\ge 0$ and implies that ${L}\text{Ψ}(s^{*},t^{*})<0$, which is contradicts with the above assumption. ${L}\text{Ψ}(s^{*},t^{*})>0,\forall s\in D$. So that, Ψ(s, t) ≥ 0, ∀(s, t) ∈ D.                                                                                                                □

Lemma 2. (Stability Result) Assume z(s, t) is the solution to the continuous problem in Equation (1). Then we have the bound

$z(s,t)\le {\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\},$

where ||f|| = max{f(s, t)}.

Proof. We prove this by using the maximum principle Lemma (1) and by constructing the barrier functions ${\theta }^{\pm }(s,t)=CM\pm z(s,t),(s,t)\in \bar{D}$, where $M={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\}$. At initial, we have

$\begin{array}{l}{\theta }^{\pm }(s,0)={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,0),{\varphi }_r(s,0)\right\}\right\}\\ \pm z(s,0)\\ ={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s)\right\}\pm {\varphi }_b(s)\ge 0.\end{array}$

For x = 0, we have

$\begin{array}{l}{\theta }^{\pm }(0,t)={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(0),max\left\{{\varphi }_l(0,t),{\varphi }_r(0,t)\right\}\right\}\\ \pm z(0,t)\\ ={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_l(t)\right\}\pm {\varphi }_l(t)\ge 0.\end{array}$

For x = 1, we have

$\begin{array}{l}{K}{\theta }^{\pm }(1,t)={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(1),max\left\{{\varphi }_l(1,t),{K}{\varphi }_r(1,t)\right\}\right\}\\ \pm {K}z(1,t)\\ ={\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_r(1,t)\right\}\pm {\varphi }_r(1,t)\ge 0.\end{array}$

For 0 < s < 1, we have

$\begin{array}{l}{L}{\theta }^{\pm }(s,t)\\ ={\theta }_t^{\pm }(s,t)-\epsilon {\theta }_{ss}^{\pm }(s,t)+a(s,t){\theta }^{\pm }(s,t),\\ ={\left[{\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\}\pm z(s,t)\right]}_t\\ -\epsilon {\left[{\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\}\pm z(s,t)\right]}_{ss}\\ +a(s,t)({\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\}\pm z(s,t))\\ =max\left\{{{\varphi }_l}_t(s,t),{{\varphi }_r}_t(s,t)\right\}\pm z_t(s,t)-\epsilon max\left\{{{\varphi }_b}_{ss}(s),{{\varphi }_l}_{ss}(s,t),{{\varphi }_r}_{ss}(s,t)\right\}\\ \pm -\epsilon u_{ss}(s,t)+\alpha {\alpha }^{-1}\left|\right|f\left|\right|+\alpha max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\}\\ \pm \alpha z(s,t)\\ \ge 0,\end{array}$

where ε > 0, a(s, t) ≥ α > 0. This implies that ${L}{\theta }^{\pm }(s,t)\ge 0$. Hence, by Lemma 1, we have, ${\theta }^{\pm }(s,t)\ge 0,\forall (s,t)\in \bar{D}$, which indicates

$z(s,t)\le {\alpha }^{-1}\left|\right|f\left|\right|+max\left\{{\varphi }_b(s),max\left\{{\varphi }_l(s,t),{\varphi }_r(s,t)\right\}\right\}.$                                     □

The sufficient conditions for the existence of a unique solution is given in Lemma 3 and Theorem 1.

Lemma 3. If the coefficient satisfies $a(s,t),f(s,t)\in C^0(\bar{D})$ and boundary conditions satisfies ${\varphi }_l\in C^1(\left[0,T\right]),{\varphi }_b\in C^{(2,1)}({\Gamma }_b),{\varphi }_r\in C^1(\left[0,T\right])$ and suppose that the compatibility conditions are satisfied. Then, the problem (Equation 1) has a unique solution z(s, t) which is satisfy $z(s,t)\in C^{(2,1)}(\bar{D})$.

Proof. Refer to Ladyženskaja et al. [23]                                                                □

To estimate the error for the fitted numerical technique below, the idea that the solution of Equation (1) is more regular than the one guaranteed by using the result in Theorem 1. To attain this greater regularity, stronger compatibility conditions are imposed at the corners.

Theorem 1. If the coefficient satisfies $a(s,t),f(s,t)\in C^{(2,1)}(\bar{D})$ and boundary conditions satisfies ${\varphi }_l\in C^2(\left[0,T\right]),{\varphi }_b\in C^{(4,2)}({\Gamma }_b),{\varphi }_r\in C^2(\left[0,T\right])$, Then the problem (Equation 1) having a unique solution z which satisfies $z\in C^{(4,2)}(\bar{D})$. And also the derivatives of solution z are bounded, ∀i, j ∈ Z ≥ 0 such that 0 ≤ i + 2j ≤ 4,

$\left|\right|\frac{{\partial }^{i+j}z}{\partial s^i\partial t^j}\left|\right|\le C{\epsilon }^{\frac{-i}{2}}.$

Proof. The boundedness of the solutions and its derivative is given as follows. Under the stretched transformation $\tilde{s}=\frac{s}{\sqrt{\epsilon }}$ problem (Equation 1) can be rewritten as

$\begin{array}{ll}\{\begin{array}{l}ℒ\tilde{z}(\tilde{s},t)=\left(-\epsilon \frac{{\partial }^2}{\partial {\tilde{s}}^2}+\frac{\partial }{\partial t}+\tilde{a}(\tilde{s},t)\right)\tilde{z}(\tilde{s},t)=f(\tilde{s},t),(\tilde{s},t)\in {\tilde{D}}_{\epsilon }\hfill \\ \tilde{z}(\tilde{s},t)={\varphi }_l(\tilde{s},t),(\tilde{s},t)\in {\tilde{\Gamma }}_l\hfill \\ K\tilde{z}(\tilde{s},t)=\tilde{z}(1,t)-\epsilon {\displaystyle {\int }_0^{1}g(s)\tilde{z}(\tilde{s},t)ds={\varphi }_r(\tilde{s},t),}(\tilde{s},t)\in {\tilde{\Gamma }}_r\hfill \\ \tilde{z}(\tilde{s},t)={\varphi }_b(\tilde{s},t),(\tilde{s},t)\in {\tilde{\Gamma }}_b\hfill \end{array}& (2)\end{array}$

where ${\tilde{D}}_{\epsilon }=(0,\frac{1}{\sqrt{\epsilon }})\times (0,T)$, and the boundary condition $\tilde{\Gamma }$ to Γ, where Equation (2) is independent of ε. Then, by taking the idea of estimation (10.6) of Ladyženskaja et al. [23] (p. 352), we will obtain

$\left|\right|\frac{{\partial }^{i+j}\tilde{z}}{\partial {\tilde{s}}^i\partial t^j}|{|}_{{Ñ}_{\delta }}\le C\left(1+\left|\right|\tilde{z}|{|}_{{Ñ}_{2\delta }}\right),$

∀ Ñ_δ in ${\tilde{D}}_{\epsilon }$. Here, Ñ_δ, δ > 0 is a neighborhood with diameter δ in ${\tilde{D}}_{\epsilon }$. Returning to the original variable

$\left|\right|\frac{{\partial }^{i+j}z}{\partial {\tilde{s}}^i\partial t^j}|{|}_{\bar{D}}\le C{\epsilon }^{\frac{-i}{2}}(1+\left|\right|z|{|}_{\bar{D}}).$

Hence, the proof is complete by using the bound on z of Lemma 2.                □

The bounds of the derivatives of the solution given in Theorem 1 were derived from classical results. They are not adequate for the proof of the ε -uniform error estimate. Stronger bounds on these derivatives are now obtained by a method originally devised in Shishkin [24]. The main idea is to decompose the solution z into smooth and singular components.

Lemma 4. If the coefficient satisfies $a(s,t),f(s,t)\in C^{(4,2)}(\bar{D})$, and the boundary conditions satisfies ${\varphi }_l\in C^{(3)}(\left[0,T\right]),{\varphi }_b\in C^{(6,3)}({\Gamma }_b),{\varphi }_r\in C^{(3)}(\left[0,T\right])$. Then we have

$\begin{array}{l}\left|\right|\frac{{\partial }^{i+j}v}{\partial s^i\partial t^j}|{|}_{\bar{D}}\le C\left(1+{\epsilon }^{1-i/2}\right),(s,t)\in D\\ \left|\frac{{\partial }^{i+j}{w}_l}{\partial s^i\partial t^j}\right|\le C{\epsilon }^{\frac{-i}{2}}{e}^{\frac{s}{\sqrt{\epsilon }}},(s,t)\in D\\ \left|\frac{{\partial }^{i+j}{w}_r}{\partial s^i\partial t^j}\right|\le C{\epsilon }^{\frac{-i}{2}}{e}^{\frac{-(1-s)}{\sqrt{\epsilon }}},(s,t)\in D\end{array}$

where C is a constant independent of parameter ε, $(s,t)\in \bar{D},i,j\ge 0,0\le i+2j\le 4.$

Proof. For proof, the interested reader can refer to Elango et al. [21].                □

3. Numerical scheme

3.1. Spatial semi-discretization

The fundamental idea of non-standard discrete modeling techniques is the development of the exact finite difference technique. Micken presented methods and rules for developing nonstandard FDMs for various types of problems [25]. To develop a discrete scheme in keeping with Mickens' guidelines, the denominator characteristic for the discrete derivatives should be described in terms of more complicated functions with larger step sizes than those used in classical methods. These complicated functions are a general property of the method that may be useful when constructing dependable methods for such problems.

To construct a genuine finite difference scheme for the problem of the form in Equation (1), we use the methods described in Woldaregay and Duressa [26]. The constant coefficient given in Equation (3) without the time variable is considered as follows.

$\begin{array}{ll}-\epsilon \frac{d^{2}z(s)}{d{s}^2}+az(s)=0.& (3)\end{array}$

By solving Equation (3), we obtain two independent solutions $e^{{\mu }_{1}s}$ and $e^{{\mu }_{2}s}$, where

${\mu }_{1,2}=\pm \sqrt{\alpha /\epsilon }.$

The spatial domain [0, 1] is discretized on uniform mesh length Δs = h as follows. $D^N=\left\{s_i=s_0+ih,i=1(1)N,s_0=0,s_N=1,h=1/N\right\}$, N is taken as number of mesh points in the spatial discretization. The approximate solution of z(s_i) will be denoted by Z_i. Here, the main aim is to compute difference equations that have similar results with the problem (Equation 1) at the mesh point s_i which is given by $Z_i=B_1{e}^{{\mu }_1{s}_i}+B_2{e}^{{\mu }_2{s}_i}$. Applying the procedures given in Mickens [25], we get

$\begin{array}{ll}det\left[\begin{array}{ccc}{Z}_{i-1}& exp({\mu }_1{s}_{i-1})& exp({\mu }_2{s}_{i-1})\\ Z_i& exp({\mu }_1{s}_i)& exp({\mu }_2{s}_i)\\ Z_{i+1}& exp({\mu }_1{s}_{i+1})& exp({\mu }_2{s}_{i+1})\end{array}\right]=0.& (4)\end{array}$

After simplification, Equation (4) becomes

$\begin{array}{ll}{Z}_{i-1}-2cosh\left(\sqrt{\frac{\alpha }{\epsilon }}h\right)Z_i+Z_{i+1}=0.& (5)\end{array}$

which is an exact difference scheme for Equation (3). By performing some arithmetic manipulation and making rearrangement on Equation (5) for the variable coefficient problem, we obtain

$\begin{array}{ll}-\epsilon \frac{Z_{i+1}-2Z_i+Z_{i-1}}{{\lambda }_i^2}+a_i{Z}_i=0.& (6)\end{array}$

The denominator function ${\lambda }_i^2$ becomes

$\begin{array}{ll}{\lambda }_i^2=\frac{4}{{\beta }_i^2}{sinh}^2(\frac{{\beta }_i}{2}h),& (7)\end{array}$

where λ² is a function of ε, β_i, h, and ${\beta }_i=\sqrt{\frac{a_i}{\epsilon }}$.

For more information about nonstandard finite difference methods for reaction diffusion problems, an interested reader can refer to the study of Munyakazi and Patidar [27].

By using Equation (7), and applying the nonstandard finite difference method to a semi-discrete problem, we have

$\begin{array}{ll}\frac{d{Z}_i(t)}{dt}-\epsilon \frac{Z_{i+1}(t)-2Z_i(t)+Z_{i-1}(t)}{{\lambda }_i^2(\epsilon ,h,t)}+a_i{Z}_i(t)=f(s_i,t).& (8)\end{array}$

with boundary conditions

$\begin{array}{ll}\{\begin{array}{l}{Z}_i={\varphi }_{li}(t),i=0,\hfill \\ Z_i={\varphi }_b,i=1(1)N-1,\hfill \\ K^N{Z}_N=Z_{N}7-\epsilon {\displaystyle {\sum }_{i=1}^N\frac{g_{i-1}{Z}_{i-1}^{j+1}+4g_i{Z}_i^{j+1}+g_{i+1}{Z}_{i+1}^{j+1}}{3}}\hfill \\ h={\varphi }_{rN},i=N.\hfill \end{array}& (9)\end{array}$

Here, for i = N, the integral boundary condition ${\displaystyle {\int }_0^{1}g(s)z(s)ds}$ approximated by composite Simpson's integration rule.

$\begin{array}{l}{\displaystyle {\int }_0^{1}g(s)z(s)ds=}\\ \frac{h}{3}(g(0)z(0)+g(N)z(N)+2{\displaystyle \sum _{i=1}^{N-1}g}(s_{2i})z(s_{2i})\\ +4{\displaystyle \sum _{i=1}^{N}g}(s_{2i-1})z(s_{2i-1}))\\ ={\varphi }_r.& (10)\end{array}$

Substituting Equation (10) in to Equation (9), we obtain

$\begin{array}{l}z(N)-\frac{h}{3}(g(0)z(0)+g(N)z(N)+2{\displaystyle \sum _{i=1}^{N-1}g}(s_{2i})z(s_{2i})\\ +4{\displaystyle \sum _{i=1}^{N}g}(s_{2i-1})z(s_{2i-1}))={\varphi }_r.& (11)\end{array}$

Equation (11) can be rewritten as

$\begin{array}{l}-\frac{4\epsilon h}{3}{\displaystyle \sum _{i=1}^N}g(s_{2i-1})z(s_{2i-1})-\frac{2\epsilon h}{3}{\displaystyle \sum _{i=1}^{N-1}}g(s_{2i})z(s_{2i})\\ +(1-\frac{\epsilon h}{3}g(N))z(N)={\varphi }_r+\frac{\epsilon h}{3}g(0)z(0).\end{array}$

Assume that the approximation of z(s_i, t) is denoted as Z_i(t), by using the non-standard finite difference approximation. At this level, Equation (1) is reduced in the form of semi-discrete as follows.

$\begin{array}{ll}\{\begin{array}{l}{ℒ}^h{Z}_i(t)=\frac{d{Z}_i(t)}{dt}\hfill \\ -\epsilon \frac{Z_{i+1}(t)-2Z_i(t)+Z_{i-1}(t)}{{\lambda }_i^2(\epsilon ,h,t)}+a_i{Z}_i(t)=f(s_i,t),\hfill \\ Z_i(0)={\varphi }_b(s_i),\hfill \\ Z_0(t)={\varphi }_l(0,t),\hfill \\ K{Z}_N(t)={\varphi }_r(N,t).\hfill \end{array}& (12)\end{array}$

Equation (12) is the system of IVPs and its compact form is written as

$\begin{array}{ll}\frac{d{Z}_i(t)}{dt}+B{Z}_i(t)=F_i(t),& (13)\end{array}$

where B is (N − 1) × (N − 1) tridiagonal matrix, Z_i(t) and F_i(t) are (N − 1) entries of the column vector. The entries of B and F respectively given as

$\{\begin{array}{l}{b}_{i,i}=\frac{2\epsilon }{{\lambda }_i^2(\epsilon ,h,t)}+a(s_i),i=1(1)N-1\hfill \\ b_{i,i-1}=-\frac{2\epsilon }{{\lambda }_i^2(\epsilon ,h,t)},i=2(1)N-1\hfill \\ b_{i,i+1}=-\frac{2\epsilon }{{\lambda }_i^2(\epsilon ,h,t)},i=1(1)N-1,\hfill \end{array}$

and

$\{\begin{array}{l}{F}_1(t)=f_1(t)-(a(s_1)+\frac{2\epsilon }{{\lambda }_1^2(\epsilon ,h,t)}){\varphi }_l(0,t),\hfill \\ F_i(t)=f_i(t),i=2(1)N-1\hfill \\ F_{N-1}(t)=f_{N-1}(t)-\frac{2\epsilon }{{\lambda }_{N-1}^2(\epsilon ,h,t)}{\varphi }_{rN}(t)\hfill \end{array}$

3.2. Stability and convergence analysis

Here, we present the maximum principle and uniform stability estimate of the semi-discrete operator ${{L}}^h$ and its convergence analysis.

Lemma 5. (Semi-discrete Maximum Principle): Assume that Z₀(t) ≥ 0, ${K}{Z}_N(t)\ge 0$. Then ${{L}}^h{Z}_i(t)\ge 0$ ∀ i = 1(1)N − 1, implies that Z_i(t) ≥ 0 ∀ i = 0(1)N.

Proof. Assume there exists q ∈ {0, ⋯ , N} such that $Z_q(t)=\underset{0\le i\le N}{min}{Z}_i(t)$. Suppose Z_q(t) ≤ 0, which implies q ≠ 0, N. Also, we have Z_q+1 − Z_q > 0 and Z_q − Z_q−1 < 0. Here, we have

${{L}}^h{Z}_q(t)=\frac{d{Z}_q(t)}{dt}-\epsilon \frac{Z_{q+1}(t)-2Z_q(t)+Z_{q-1}(t)}{{\lambda }_q^2}+a_q{Z}_q(t).$

By using the above assumption, we get that ${{L}}^h{Z}_i(t)<0$, for i = 1(1)N − 1. Thus, the assumption Z_i(t) < 0, i = 0(1)N is not correct. Hence, Z_i(t) ≥ 0 ∀ i = 0(1)N.                □

This Lemma 5 is used to obtain the bounds of the discrete solution given in Lemma 6. In general, the discrete maximum principle is widely used to show the boundedness and positivity of a discrete solution.

Lemma 6. The solution Z_i(t) of the semidiscrete problem in Equations (12) or (13) satisfies the following bound.

$\begin{array}{l}\left|Z_i(t)\right|=\frac{1}{\alpha }{\displaystyle \underset{i}{max}}\left|{{L}}^h{Z}_i(t)\right|+{\displaystyle \underset{i}{max}}\left\{{\varphi }_b(s_i),{\displaystyle \underset{i}{max}}\left\{{\varphi }_l(s_i,t),{\varphi }_r(s_i,t)\right\}\right\}.\end{array}$

Proof. Suppose $q=\frac{1}{\alpha }\underset{i}{max}\left|{{L}}^h{Z}_i(t)\right|+\underset{i}{max}\left\{{\varphi }_b(s_i),\underset{i}{max}\left\{{\varphi }_l(s_i,t),{\varphi }_r(s_i,t)\right\}\right\}$ and define a comparison function ${\theta }_i^{\pm }(t)$ as

${\theta }_i^{\pm }(t)=q\pm Z_i(t).$

For the points on the boundary, we have

$\begin{array}{l}{\theta }_0^{\pm }(t)=q\pm Z_0(t)=q\pm {\varphi }_l(0,t)\ge 0,\\ {{K}}^N{\theta }_N^{\pm }(t)=q\pm {{K}}^N{Z}_N(t)=q\pm {\varphi }_r(1,t)\ge 0.\end{array}$

For 1 ≤ i ≤ N − 1, we have

$\begin{array}{l}{ℒ}^h{\theta }_i^{\pm }(t)=\frac{d(q\pm Z_i(t))}{dt}\\ -\epsilon \frac{\left(q\pm Z_{i-1}(t)-2(q\pm Z_i(t))+q\pm Z_{i+1}(t)\right)}{{\lambda }^2}\\ +a_i(q\pm Z_i(t))\\ =a_{i}q\pm {ℒ}^h{Z}_i(t)\\ =a_i({\alpha }^{-1}\underset{i}{\max }\left|{ℒ}^h{Z}_i(t)\right|\\ +\underset{i}{\max }\left\{{\varphi }_b(s_i),\underset{i}{\max }\left\{{\varphi }_l(s_i,t),{\varphi }_r(s_i,t)\right\}\right\})\pm f_{i,j}\\ \ge 0,\text{ since }{a}_i\ge \alpha .\end{array}$

From Lemma 5, we get, ${\theta }_i^{\pm }(t)\ge 0$, ∀ $(s_i,t)\in {\bar{\Omega }}_x^N\times (0,T)$.                □

Next, we present the convergence analysis of spatial discretization. We denoted Z_i(t) as approximate solution for the spatial semidiscretization to the exact solution z(s, t) at s = s_i, i = 0(1)N. Let us define the backward and forward finite differences in space as:

$D^{-}z(s_i,t)=\frac{z(s_i,t)-z(s_{i-1})}{h},D^{+}z(s_i,t)=\frac{z(s_{i+1},t)-z(s_i,t)}{h},$

respectively, and the second order central finite difference operator as

${\delta }^{2}z(s_i,t)=D^{+}{D}^{-}z(s_i,t)=\frac{D^{+}z(s_i,t)-D^{-}z(s_i,t)}{h}.$

Lemma 7. Let N be a fixed mesh. Then, for ε → 0, we have

$\begin{array}{l}{\displaystyle \underset{\epsilon \to 0}{lim}}{\displaystyle \underset{1\le i\le N-1}{max}}\frac{exp(-p{s}_i/\sqrt{\epsilon })}{{\epsilon }^{m/2}}=0\text{and}\\ {\displaystyle \underset{\epsilon \to 0}{lim}}{\displaystyle \underset{1\le i\le N-1}{max}}\frac{exp(-p(1-s_i)/\sqrt{\epsilon })}{{\epsilon }^{m/2}}=0.\end{array}$

where m = 1, 2, 3, ⋯ .

Proof. Refer to Munyakazi and Patidar [27]                □

Theorem 2. Let the coefficient function a(s) and the function f(s, t) in Equation (12) be sufficiently smooth and $z(s,t)\in C^4(\bar{D})$. Then the semidiscrete solution Z_i(t) of Equation (12) satisfies

$\begin{array}{l}\left|{{L}}^h(z(s_i,t)-Z_i(t))\right|\le C{h}^2.\end{array}$

Proof. The truncation error in spatial direction is considered as

$\begin{array}{l}{{L}}^h(z(s_i,t)-Z_i(t))={{L}}^{h}z(s_i,t)-{{L}}^h{Z}_i(t)\\ =-\epsilon \frac{d^2}{d{s}^2}z(s_i,t)+\frac{D_s^{+}{D}_s^{+}{h}^2}{{\lambda }^2}z(s_i,t)\\ =-\epsilon \frac{d^2}{d{s}^2}z(s_i,t)+\frac{\epsilon }{{\lambda }^2}(h^2\frac{d^2}{d{s}^2}z(s_i,t)+\frac{h^4}{12}\frac{d^4}{d{s}^4}z(s_i,t)).& (14)\end{array}$

Note that we have used Taylor expansions of z_i−1(t) and z_i+1(t). A truncated Taylor expansion of $\frac{1}{{\lambda }^2}$ of order five becomes

$\begin{array}{ll}\frac{1}{{\lambda }^2}=\frac{{\beta }^2}{4}\left(\frac{4}{{\beta }^2{h}^2}-\frac{1}{3}+\frac{{\beta }^2{h}^2}{60}\right).& (15)\end{array}$

Using Equation (15) in Equation (14), we obtain

$\begin{array}{l}{{L}}^h(z(s_i,t)-Z_i(t))=\frac{\epsilon }{12}\left(\frac{d^4}{d{s}^4}z(s_i,t)-{\beta }^2\frac{d^2}{d{s}^2}z(s_i,t)\right)h^2\\ +\epsilon {\beta }^2{h}^4\left(\frac{{\beta }^2}{240}\frac{d^{2}z(s_i,t)}{d{s}^2}-\frac{1}{144}\frac{d^{4}z(s_i,t)}{d{s}^4}\right)+h^6\frac{\epsilon {\beta }^4}{2880}\frac{d^{4}z(s_i,t)}{d{s}^4}.& (16)\end{array}$

We use Lemma (7), to obtain the boundedness of Equation (16). Using Lemma (7) and Theorem (1), we obtain

$\left|{{L}}^h(z(s_i,t)-Z_i(t))\right|\le C{N}^{-2}.$

The truncation error at s = s_N, become

$\begin{array}{l}{K}^N\left(Z(s_i)-z(s_i)\right)=K^{N}Z(s_N)-K^{N}s(s_i),\\ ={\varphi }_r-K^{N}Z(s_N),\\ =Kz(s_i)-K^{N}Z(s_N),\\ =z(s_N)-\epsilon {\displaystyle {\int }_0^{1}g(s)z(s)ds}-\left(Z(s_N)-\epsilon {\displaystyle {\int }_{s_0}^{s_N}g(s)z(s)ds}\right),\\ =\epsilon {\displaystyle {\int }_{s_0}^{s_N}g(s)z(s)ds}-\epsilon {\displaystyle \sum _{i=1}^N\frac{g_{i-1}{z}_{i-1}+4g_i{z}_i+g_{i+1}{z}_{i+1}}{3}}h,\\ =\epsilon [{\displaystyle {\int }_{s_0}^{s_1}g(s)z(s)ds}+{\displaystyle {\int }_{s_1}^{s_2}g(s)z(s)ds}+\cdots \\ +{\displaystyle {\int }_{s_N}^{s_{N+1}}g(s)z(s)ds}]\\ -\epsilon [\frac{g_0{z}_0+4g_1{z}_1+g_2{z}_2}{3}h+\cdots \\ +\frac{g_{N-1}{z}_{N-1}+4g_N{z}_N+g_{N+1}{z}_{N+1}}{3}h],\\ \left|K^N\left(Z(s_i)-z(s_i)\right)\right|\\ =\left|C\epsilon \left(h^4{z}^{(4)}({\xi }_1)+h^4{z}^{(4)}({\xi }_2)+\cdots +h^4{z}^{(4)}({\xi }_N)\right)\right|,\\ \left|K^N\left(Z(s_i)-z(s_i)\right)\right|\\ \le C\epsilon h^4\left(z^{(4)}({\xi }_1)+z^{(4)}({\xi }_2)+\cdots +z^{(4)}({\xi }_N)\right),\\ \le C\epsilon h^4\left|\left|\frac{d^{4}z({\xi }_i)}{d{x}^4}\right|\right|\le C{h}^2=C{N}^{-2}.\square & (17)\end{array}$

Theorem 3. The semidiscrete solutions satisfy the uniform error bound

$\begin{array}{ll}{\displaystyle \underset{0<\epsilon \ll 1}{sup}}{\displaystyle \underset{i}{max}}|z(s_i,t)-Z_i(t){|}_{\bar{D}}\le C{N}^{-2}.& (18)\end{array}$

Proof. The proof follows from Theorem (1) and Lemma (7) under the properties of boundedness of a semi-discrete solution and the required bound is satisfied.

3.3. Temporal discretization

A mesh with length Δt = t_j+1 − t_j,j = 0(1)M − 1 is constructed on the time domain D_t = [0, T], where M is a positive integer. The IVPs Equation (13) are discretized using the implicit Euler method on a uniform mesh. By denoting the approximation of z_i(t_j) as $Z_i^j$, we construct the time discretization as follows.

$\begin{array}{ll}\frac{Z_i^j-Z_i^{j-1}}{\Delta t}=B{Z}_i^j+F_i^j& (19)\end{array}$

with the initial condition Z₀(t) = ϕ_l(t_j), and by rearranging Equation (19), we obtain

$\begin{array}{ll}{Z}_i^j={\left[I+\Delta tB\right]}^{-1}\left[\Delta t{F}_i^j+Z_i^{i-1}\right].& (20)\end{array}$

Lemma 8. Suppose $\left|\frac{{\partial }^{i}z(s,t_j)}{\partial t^i}\right|\le C,\forall (s,t)\in \bar{D},i=0,1,2.$ Then the local truncation error associated with the time direction satisfies $\left|e_j\right|\le C{(\Delta t)}^2.$

Proof. The local truncation error is defined as

$\begin{array}{l}{e}_j=z(t_j)-Z_i^j\\ =z(t_j)-{\left[I+\Delta tB\right]}^{-1}\left[\Delta t{F}_i^j+Z_i^{j-1}\right].\end{array}$

Using Taylor expansion, we obtain z(t_j−1) as

$z(t_{j-1})=z(t_j)-\Delta t{z}_t(t_j)+\frac{{(\Delta t)}^2}{2}{z}_{tt}(t_j)+\frac{{(\Delta t)}^3}{3!}{z}_{ttt}(t_j)+{O}({(\Delta t)}^4).$

However, z_t(t_j) = F(t_j) − B(t_j)z(t_j). Thus,

$\begin{array}{l}z(t_{j-1})=z(t_j)-\Delta t\left[F(t_j)-B(t_j)z(t_j)\right]+\frac{{(\Delta t)}^2}{2}{z}_{tt}(t_j)\\ +\frac{{(\Delta t)}^3}{3!}{z}_{ttt}(t_j)+{O}({(\Delta t)}^4).\end{array}$

Now, the local truncation error e_j becomes

$\begin{array}{l}{e}_j=z(t_j)-{\left[I+\Delta tB\right]}^{-1}\left[\Delta t{F}_i^j+Z_i^{j-1}\right]\\ =z(t_j)-{\left[I+\Delta tB\right]}^{-1}\left[\left[I+\Delta tB(t_j)\right]z(t_j)+\frac{{(\Delta t)}^2}{2}{z}_{tt}(t_j)+\cdots \right]\\ ={\left[I+\Delta tB\right]}^{-1}\left[\frac{{(\Delta t)}^2}{2}{z}_{tt}(t_j)-\frac{{(\Delta t)}^3}{3!}{z}_{ttt}(t_j)+{O}({(\Delta t)}^4)\right].\end{array}$

Since the matrix B is invertible, using the relation (Δt)² > (Δt)³ for small Δt and z(t_j) ≤ C, we obtain

$\begin{array}{l}\left|\right|e_j\left|\right|\le \left|\right|{\left[I+\Delta tB\right]}^{-1}\left|\right|\left|\right|\frac{{(\Delta t)}^2}{2}{z}_{tt}(t_j)-\frac{{(v)}^3}{3!}{z}_{ttt}(t_j)+{O}({(\Delta t)}^4)\left|\right|\\ \le \left|\right|{\left[I+\Delta tB\right]}^{-1}\left|\right|{(\Delta t)}^2\le C{(\Delta t)}^2.\end{array}$

                                                                                                                □

Lemma 9. The global error estimate in the time direction is given by ||E_j+1|| ≤ CΔt, ∀j ≤ T/Δt, where $E_{j+1}=\underset{i}{max}|Z_i(t_{j+1})-Z_{i,j+1}{|}_D$.

Proof. The global error estimate at (j+1)^th time step is obtained by using the local error estimate up to j^th time step as follows.

$\begin{array}{l}\left|\right|E_{j+1}\left|\right|=\left|\right|{\displaystyle \sum _{i=1}}{e}_j\left|\right|j\le T/\Delta t\\ \le \left|\right|e_1\left|\right|+\left|\right|e_2\left|\right|+\left|\right|e_3\left|\right|+\left|\right|e_4\left|\right|+\cdots +\left|\right|e_j\left|\right|\\ \le C_1(j\Delta t)\Delta t\\ \le C_{1}T\Delta t\text{since }j\Delta t\le T\\ \le C\Delta t.\end{array}$

Hence,

$\begin{array}{ll}\left|\right|E^{j+1}\left|\right|={\displaystyle \underset{i}{max}}|Z_i(t_{j+1})-Z_i^{j+1}{|}_D\le C\Delta t.& (21)\end{array}$

where C is a positive constant independent of ε and Δt. By taking the supremum ∀ ε ∈ (0, 1], we obtained

$\begin{array}{ll}{\displaystyle \underset{0<\epsilon \ll 1}{sup}}{\displaystyle \underset{i}{max}}|Z_i(t_{j+1})-Z_i^{j+1}{|}_D\le C\Delta t.& (22)\end{array}$

                                                                                                                □

We summarizes the results of this work by considering the error estimate obtained in Equations (18) and (22) and we conclude by the following theorem.

Theorem 4. The error estimate for the solution of the continuous and fully discrete problems is given by

${\displaystyle \underset{0<\epsilon <<1}{sup}}{\displaystyle \underset{0\le i\le N}{max}}{\displaystyle \underset{0\le i\le M}{max}}\left|\right|z(s,t)-Z_i^{j+1}\left|\right|\le C\left(N^{-2}+\Delta t\right),$

where z(s, t) and $Z_i^{j+1}$ are the solutions to problems Equations (1) and (12), respectively.

Proof. The error estimation of the fully discrete scheme is given as follows.

$\begin{array}{l}{\displaystyle \underset{\epsilon }{sup}}{\displaystyle \underset{i,j}{max}}|z(s_i,t_j)-Z_i^j|={\displaystyle \underset{\epsilon }{sup}}{\displaystyle \underset{i,j}{max}}|z(s_i,t_j)-Z_i(t_j)+Z_i(t_j)-Z_i^j|\\ \le {\displaystyle \underset{\epsilon }{sup}}{\displaystyle \underset{i,j}{max}}|z(s_i,t_j)-Z_i(t_j)|+{\displaystyle \underset{\epsilon }{sup}}{\displaystyle \underset{i,j}{max}}|Z_i(t_j)-Z_i^j|.\end{array}$

Then, by combining the bound given in Theorem 3 and Lemma 9, the theorem gets proved.

4. Numerical examples, results, and discussions

Here, we developed an algorithm for the proposed method for the problem and perform experiments to validate the theoretical justifications and results. Since the exact solutions of the given examples are not known, we use double mesh techniques to obtain the maximum pointwise error of the developed scheme. Now, let U^N,Δt be a conducted solution of a problem using mesh points N and time step size Δt. Again, $U_{i,j}^{2N,\Delta t/2}$ be a conducted solution on double mesh points of 2N and half of time step size Δt/2.

We calculate the maximum absolute error as $E_{\epsilon }^{N,\Delta t}=\underset{i,j}{max}|Z_{i,j}^{N,\Delta t}-Z_{i,j}^{2N,\Delta t/2}|,$ and the parameter uniform error estimate by using $E^{N,\Delta t}=\underset{\epsilon }{max}(E_{\epsilon }^{N,\Delta t}).$ We calculate the rate of convergence of the developed scheme by using $P_{\epsilon }^{N,\Delta t}={log}_2(E_{\epsilon }^{N,\Delta t})-{log}_2(E_{\epsilon }^{2N,\Delta t/2}).$ The parameter rate of convergence is calculated as $P^{N,\Delta t}={log}_2(E^{N,\Delta t})-{log}_2(E^{2N,\Delta t/2}).$

Example 1.

$\{\begin{array}{l}\frac{\partial z(s,t)}{\partial t}-\epsilon \frac{{\partial }^{2}z(s,t)}{\partial s^2}+\frac{1+s^2}{2}z(s,t)=e^{-t}-1\hfill \\ +\sin (\pi s),(s,t)\in (0,1)\times (0,1]\hfill \\ z(s,0)=0,s\in (0,1),\hfill \\ z(0,t)=0,t\in (0,1],\hfill \\ Kz(1,t)=z(1,t)-\epsilon {\displaystyle {\int }_0^1\frac{s}{6}z(s,t)ds=0,t\in (0,1].}\hfill \end{array}$

Example 2.

$\{\begin{array}{l}\frac{\partial z(s,t)}{\partial t}-\epsilon \frac{{\partial }^{2}z(s,t)}{\partial s^2}+\frac{1+s^2}{2}z(s,t)=t^3,(s,t)\in (0,1)\times (0,1]\hfill \\ z(s,0)=0,s\in (0,1),\hfill \\ z(0,t)=0,t\in (0,1],\hfill \\ Kz(1,t)=z(1,t)-\epsilon {\displaystyle {\int }_0^1\cos (s)z(s,t)ds=0,}t\in (0,1].\hfill \end{array}$

The solutions of the above two examples exhibit strong boundary layers near x = 0 and x = 1. We presented the surface plots for numerical solutions of Examples 1 and 2 in Figures 1, 2 respectively, which display the presence of boundary layers formation on the left and right side of the spatial domain for different values of ε. The maximum pointwise error and rate of convergence of the proposed schemes of Examples 1 and 2 are given in Tables 1, 2 respectively for various values of the perturbation parameter ε, mesh number N and time step size Δt. From these tables, one can observe that the developed scheme is parameter uniform convergent, which supports the theoretical results. Figure 3 indicates the Log-Log plots for the maximum absolute error vs. mesh number N for the singular perturbation parameter ε. One can observe that as ε goes very small, the developed method converges uniformly independent of the perturbation parameter ε.

FIGURE 1

Figure 1. 3-D graph of numerical solution for Example (1) which displays the existing layer. (A) ε = 10⁻². (B) ε = 10⁻¹².

FIGURE 2

Figure 2. 3-D graph of numerical solution for Example (2) that displays the existing layer. (A) ε = 10⁻². (B) ε = 10⁻¹².

TABLE 1

Table 1. Maximum absolute error and rate of convergence of the scheme for Example (1).

TABLE 2

Table 2. Maximum absolute error and rate of convergence of the scheme for Example (2).

FIGURE 3

Figure 3. The Log-Log plot of the maximum absolute error for different values of ε for Examples 1 and 2, respectively. (A) Log-Log plot for Example (1). (B) Log-Log plot for Example (2).

5. Conclusion

This paper investigates a numerical treatment for a class of singularly perturbed parabolic partial differential equations of the reaction-diffusion type with nonlocal boundary conditions. To solve the problem at hand, we employed the method of lines. A nonstandard finite difference scheme is used to semi-discretize the spatial direction, and the implicit Euler method is used to discretize the results of initial value problems. To deal with the integral boundary condition, we used a composite Simpson's $\frac{1}{3}$ rule. The stability of the evolved numerical scheme is established, and the scheme's uniform convergence is demonstrated. To validate the problem's applicability, two test examples are carried out for numerical computation for different values of the perturbation parameter ε and mesh points. The entire procedure has been demonstrated to be second-order uniformly convergent in the spatial direction and first-order in the temporal direction. The theoretical estimation is reflected in our numerical results.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Acknowledgments

Authors are grateful to their reviewers for their contributions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Kumar M, Paul. Methods for solving singular perturbation problems arising in science and engineering. Math Comput Model. (2011) 54:556–75. doi: 10.1016/j.mcm.2011.02.045

CrossRef Full Text | Google Scholar

2. Dehghan M, Tatari M. Use of radial basis functions for solving the second-order parabolic equation with nonlocal boundary conditions. Math Comput Simulat. (2008) 24:924–38. doi: 10.1002/num.20297

CrossRef Full Text | Google Scholar

3. Amiraliyev G, Amiraliyeva I, Kudu M. A numerical treatment for singularly perturbed differential equations with integral boundary condition. Appl Math Comput. (2007) 185:574–82. doi: 10.1016/j.amc.2006.07.060

CrossRef Full Text | Google Scholar

4. Day W. Extensions of a property of the heat equation to linear thermoelasticity and other theories. Q Appl Math. (1982) 40:319–30. doi: 10.1090/qam/678203

CrossRef Full Text | Google Scholar

5. Bouziani A. Mixed problem with boundary integral conditions for a certain parabolic equation. J Appl Math Stochastic Anal. (1996) 9:323–30. doi: 10.1155/S1048953396000305

CrossRef Full Text | Google Scholar

6. Cannon J. The solution of the heat equation subject to the specification of energy. Q Appl Math. (1963) 21:155–60. doi: 10.1090/qam/160437

CrossRef Full Text | Google Scholar

7. Kudu M, Amiraliyev GM. Finite difference method for a singularly perturbed differential equations with integral boundary condition. Int J Math Comput. (2015) 26:71–9.

Google Scholar

8. OMalley RE. Ludwig prandtls boundary layer theory. In: Historical Developments in Singular Perturbations. Springer (2014). p. 1–26.

Google Scholar

9. Bahuguna D, Dabas J. Existence and uniqueness of a solution to a semilinear partial delay differential equation with an integral condition. Nonlinear Dyn Syst Theory. (2008) 8:7–19.

Google Scholar

10. Feng M, Ji D, Ge W. Positive solutions for a class of boundary-value problem with integral boundary conditions in Banach spaces. J Comput Appl Math. (2008) 222:351–63. doi: 10.1016/j.cam.2007.11.003

CrossRef Full Text | Google Scholar

11. Li H, Sun F. Existence of solutions for integral boundary value problems of second-order ordinary differential equations. Boundary Value Problems. (2012) 2012:1–7. doi: 10.1186/1687-2770-2012-147

CrossRef Full Text | Google Scholar

12. Raja V, Tamilselvan A. Fitted finite difference method for third order singularly perturbed convection diffusion equations with integral boundary condition. Arab J Math Sci. (2019) 25:231–42. doi: 10.1016/j.ajmsc.2018.10.002

CrossRef Full Text | Google Scholar

13. Debela, Duressa GF. Uniformly convergent numerical method for singularly perturbed convection-diffusion type problems with nonlocal boundary condition. Int J Num Methods Fluids. (2020) 92:1914–26. doi: 10.1002/fld.4854

CrossRef Full Text | Google Scholar

14. Sekar E, Tamilselvan A. Finite difference scheme for singularly perturbed system of delay differential equations with integral boundary conditions. J Korean Soc Ind Appl Math. (2018) 22:201–15. doi: 10.12941/jksiam.2018.22.201

CrossRef Full Text | Google Scholar

15. Sekar E, Tamilselvan A. Finite difference scheme for third order singularly perturbed delay differential equation of convection diffusion type with integral boundary condition. J Appl Math Comput. (2019) 61:73–86. doi: 10.1007/s12190-019-01239-0

CrossRef Full Text | Google Scholar

16. Sekar E, Tamilselvan A. Singularly perturbed delay differential equations of convection-diffusion type with integral boundary condition. J Appl Math Comput. (2019) 59:701–22. doi: 10.1007/s12190-018-1198-4

CrossRef Full Text | Google Scholar

17. Sekar E, Tamilselvan A. Third order singularly perturbed delay differential equation of reaction diffusion type with integral boundary condition. J Appl Math Comput Mech. (2019) 18:99–110. doi: 10.17512/jamcm.2019.2.09

CrossRef Full Text | Google Scholar

18. Debela, Duressa GF. Exponentially fitted finite difference method for singularly perturbed delay differential equations with integral boundary condition. Int J Eng Appl Sci. (2019) 11:476–93. doi: 10.24107/ijeas.647640

CrossRef Full Text | Google Scholar

19. Debela, Duressa GF. Accelerated fitted operator finite difference method for singularly perturbed delay differential equations with non-local boundary condition. J Egypt Math Soc. (2020) 28:1–16. doi: 10.1186/s42787-020-00076-6

CrossRef Full Text | Google Scholar

20. Kumar D, Kumari P. A parameter-uniform collocation scheme for singularly perturbed delay problems with integral boundary condition. J Appl Math Comput. (2020) 63:813–28. doi: 10.1007/s12190-020-01340-9

CrossRef Full Text | Google Scholar

21. Elango S, Tamilselvan A, Vadivel R, Gunasekaran N, Zhu H, Cao J, et al. Finite difference scheme for singularly perturbed reaction diffusion problem of partial delay differential equation with nonlocal boundary condition. Adv Diff Equat. (2021) 2021:1–20. doi: 10.1186/s13662-021-03296-x

CrossRef Full Text | Google Scholar

22. Gobena WT, Duressa GF. Parameter-uniform numerical scheme for singularly perturbed delay parabolic reaction diffusion equations with integral boundary condition. Int J Diff Equat. (2021) 2021:9993644. doi: 10.1155/2021/9993644

CrossRef Full Text | Google Scholar

23. Ladyženskaja OA, Solonnikov VA, Ural'ceva NN. Linear and quasi-linear equations of parabolic type. In: American Mathematical Society, Vol. 23. Petersburg, VA (1988).

Google Scholar

24. Shishkin GI. Approximation of the solutions of singularly perturbed boundary-value problems with a parabolic boundary layer. USSR Comput Math Math Phys. (1989) 29:1–10. doi: 10.1016/0041-5553(89)90109-2

CrossRef Full Text | Google Scholar

25. Mickens RE. Advances in the applications of nonstandard finite difference schemes. In: World Scientific. (2005).

Google Scholar

26. Woldaregay MM, Duressa GF. Uniformly convergent numerical method for singularly perturbed delay parabolic differential equations arising in computational neuroscience. Kragujevac J Math. (2022) 46:65–84. doi: 10.46793/KgJMat2201.065W

CrossRef Full Text | Google Scholar

27. Munyakazi JB, Patidar KC. A fitted numerical method for singularly perturbed parabolic reaction-diffusion problems. Comput Appl Math. (2013) 32:509–19. doi: 10.1007/s40314-013-0033-7

CrossRef Full Text | Google Scholar