Skip to main content


Front. Appl. Math. Stat., 25 October 2022
Sec. Numerical Analysis and Scientific Computation
Volume 8 - 2022 |

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

  • 1Department of Applied Mathematics, Adama Science and Technology University, Adama, Ethiopia
  • 2Department of Mathematics, Jimma University, Jimma, Ethiopia

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 13 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 [1417]. 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,

{z(s,t)=(ε2s2+t+a(s,t))z(s,t)=f(s,t)  (s,t)D,z(s,0)=ϕb(s),  ϕb(s,t)Γb={(s,0)},z(0,t)=ϕl(t),  ϕl(s,t)Γl={(0,t);0tT},Kz(s,t)=z(1,t)ε01g(s)z(s,t)ds=ϕr(s,t),  ϕr(s,t)  Γr={(1,t);0tT}.    (1)

where (s,t)D=Ωx×Ωt=(0,1)×(0,T], D̄=[0,1]×[0,T], 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 01g(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


Note that ϕl, ϕr, and ϕb are assumed to be sufficiently smooth for Equation (1) to make sense, namely ϕl,ϕrC1([0,T]), and ϕbC(2,1)(Γ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 Ψ(s,t)C(0,0)(D̄)C(1,0)(D)C(2,1)(D) be a sufficiently smooth function such that Ψ(0,t)0,Ψ(s,0)0,KΨ(1,t)0,LΨ(s,t)0,(s,t)D. Then Ψ(s, t) ≥ 0, (s,t)D̄, where LΨ(s,t)=Ψt(s,t)-εΨss(s,t)+aΨ(s,t).

Proof. Assume (s*, t*) be defined as Ψ(s*,t*)=min(s,t)D̄Ψ(s,t) and suppose that Ψ(s*, t*) ≤ 0. It is known that (s*, t*) ∉ ∂D. Thus, LΨ(s*,t*)=Ψt(s*,t*)-εΨss(s*,t*)+a(s,t)Ψ(s*,t*). Since Ψ(s*,t*)=min(s,t)D̄Ψ(s,t), which indicates that Ψ(s*, t*) = 0, Ψt(s*,t*)=0, Ψss(s*,t*)0 and implies that LΨ(s*,t*)<0, which is contradicts with the above assumption. LΨ(s*,t*)>0,sD. 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


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

Proof. We prove this by using the maximum principle Lemma (1) and by constructing the barrier functions θ±(s,t)=CM±z(s,t),(s,t)D̄, where M=α-1||f||+max{ϕb(s),max{ϕl(s,t),ϕr(s,t)}}. At initial, we have

θ±(s,0)=α-1||f||+max{ϕb(s),max{ϕl(s,0),ϕr(s,0)}}                   ±z(s,0)               =α-1||f||+max{ϕb(s)}±ϕb(s)0.

For x = 0, we have

θ±(0,t)=α-1||f||+max{ϕb(0),max{ϕl(0,t),ϕr(0,t)}}                   ±z(0,t)                =α-1||f||+max{ϕl(t)}±ϕl(t)0.

For x = 1, we have

Kθ±(1,t)=α-1||f||+max{ϕb(1),max{ϕl(1,t),Kϕr(1,t)}}                   ±Kz(1,t)                   =α-1||f||+max{ϕr(1,t)}±ϕr(1,t)0.

For 0 < s < 1, we have


where ε > 0, a(s, t) ≥ α > 0. This implies that Lθ±(s,t)0. Hence, by Lemma 1, we have, θ±(s,t)0,(s,t)D̄, which indicates

z(s,t)α-1||f||+max{ϕb(s),max{ϕl(s,t),ϕr(s,t)}}.                                     □

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)C0(D̄) and boundary conditions satisfies ϕlC1([0,T]),ϕbC(2,1)(Γb),ϕrC1([0,T]) 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)C(2,1)(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)C(2,1)(D̄) and boundary conditions satisfies ϕlC2([0,T]),ϕbC(4,2)(Γb),ϕrC2([0,T]), Then the problem (Equation 1) having a unique solution z which satisfies zC(4,2)(D̄). And also the derivatives of solution z are bounded, ∀i, jZ ≥ 0 such that 0 ≤ i + 2j ≤ 4,


Proof. The boundedness of the solutions and its derivative is given as follows. Under the stretched transformation s~=sε problem (Equation 1) can be rewritten as

{z˜(s˜,t)=(ε2s˜2+t+a˜(s˜,t))z˜(s˜,t)=f(s˜,t),   (s˜,t)D˜εz˜(s˜,t)=ϕl(s˜,t),                                               (s˜,t)Γ˜lKz˜(s˜,t)=z˜(1,t)ε01g(s)z˜(s˜,t)ds=ϕr(s˜,t),        (s˜,t)Γ˜rz˜(s˜,t)=ϕb(s˜,t),                                                (s˜,t)Γ˜b    (2)

where D~ε=(0,1ε)×(0,T), and the boundary condition Γ~ 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


∀ Ñδ in D~ε. Here, Ñδ, δ > 0 is a neighborhood with diameter δ in D~ε. Returning to the original variable


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)C(4,2)(D̄), and the boundary conditions satisfies ϕlC(3)([0,T]),ϕbC(6,3)(Γb),ϕrC(3)([0,T]). Then we have

||i+jvsitj||D̄C(1+ε1-i/2),  (s,t)D|i+jwlsitj|Cε-i2esε,           (s,t)D|i+jwrsitj|Cε-i2e-(1-s)ε,       (s,t)D

where C is a constant independent of parameter ε, (s,t)D̄, i,j0,0i+2j4.

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.

-εd2z(s)ds2+az(s)=0.    (3)

By solving Equation (3), we obtain two independent solutions eμ1s and eμ2s, where


The spatial domain [0, 1] is discretized on uniform mesh length Δs = h as follows. DN={si=s0+ih,  i=1(1)N,s0=0,sN=1,h=1/N}, N is taken as number of mesh points in the spatial discretization. The approximate solution of z(si) will be denoted by Zi. Here, the main aim is to compute difference equations that have similar results with the problem (Equation 1) at the mesh point si which is given by Zi=B1eμ1si+B2eμ2si. Applying the procedures given in Mickens [25], we get

det[Zi-1exp(μ1si-1)exp(μ2si-1)Ziexp(μ1si)exp(μ2si)Zi+1exp(μ1si+1)exp(μ2si+1)]=0.    (4)

After simplification, Equation (4) becomes

Zi-1-2cosh(αεh)Zi+Zi+1=0.    (5)

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

-εZi+1-2Zi+Zi-1λi2+aiZi=0.    (6)

The denominator function λi2 becomes

λi2=4βi2sinh2(βi2h),    (7)

where λ2 is a function of ε, βi, h, and βi=aiε.

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

dZi(t)dt-εZi+1(t)-2Zi(t)+Zi-1(t)λi2(ε,h,t)+aiZi(t)=f(si,t).    (8)

with boundary conditions

{Zi=ϕli(t),  i=0,Zi=ϕb,                   i=1(1)N1,KNZN=ZN7εi=1Ngi1Zi1j+1+4giZij+1+gi+1Zi+1j+13  h=ϕrN, i=N.    (9)

Here, for i = N, the integral boundary condition 01g(s)z(s)ds approximated by composite Simpson's integration rule.

01g(s)z(s)ds=h3(g(0)z(0)+g(N)z(N)+2i=1N1g(s2i)z(s2i)+4i=1Ng(s2i1)z(s2i1))=ϕr.    (10)

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

z(N)h3(g(0)z(0)+g(N)z(N)+2i=1N1g(s2i)z(s2i)+4i=1Ng(s2i1)z(s2i1))=ϕr.    (11)

Equation (11) can be rewritten as


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

{hZi(t)=dZi(t)dtεZi+1(t)2Zi(t)+Zi1(t)λi2(ε,h,t)+aiZi(t)=f(si,t),Zi(0)=ϕb(si),Z0(t)=ϕl(0,t),KZN(t)=ϕr(N,t).    (12)

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

dZi(t)dt+BZi(t)=Fi(t),    (13)

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

{bi,i=2ελi2(ε,h,t)+a(si),    i=1(1)N1bi,i1=2ελi2(ε,h,t),    i=2(1)N1bi,i+1=2ελi2(ε,h,t),    i=1(1)N1,


{F1(t)=f1(t)(a(s1)+2ελ12(ε,h,t))ϕl(0,t),Fi(t)=fi(t),    i=2(1)N1FN1(t)=fN1(t)2ελN12(ε,h,t)ϕrN(t)

3.2. Stability and convergence analysis

Here, we present the maximum principle and uniform stability estimate of the semi-discrete operator Lh and its convergence analysis.

Lemma 5. (Semi-discrete Maximum Principle): Assume that Z0(t) ≥ 0, KZN(t)0. Then LhZi(t)0i = 1(1)N − 1, implies that Zi(t) ≥ 0 ∀ i = 0(1)N.

Proof. Assume there exists q ∈ {0, ⋯ , N} such that Zq(t)=min0iNZi(t). Suppose Zq(t) ≤ 0, which implies q ≠ 0, N. Also, we have Zq+1Zq > 0 and ZqZq−1 < 0. Here, we have


By using the above assumption, we get that LhZi(t)<0, for i = 1(1)N − 1. Thus, the assumption Zi(t) < 0, i = 0(1)N is not correct. Hence, Zi(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 Zi(t) of the semidiscrete problem in Equations (12) or (13) satisfies the following bound.


Proof. Suppose q=1αmaxi|LhZi(t)|+maxi{ϕb(si),maxi{ϕl(si,t),ϕr(si,t)}} and define a comparison function θi±(t) as


For the points on the boundary, we have


For 1 ≤ iN − 1, we have

hθi±(t)=d(q±Zi(t))dt                 ε(q±Zi1(t)2(q±Zi(t))+q±Zi+1(t))λ2                 +ai(q±Zi(t))                 =aiq±hZi(t)                 =ai(α1maxi|hZi(t)|                 +maxi{ϕb(si),maxi{ϕl(si,t),ϕr(si,t)}})±fi,j                  0,      since  aiα.

From Lemma 5, we get, θi±(t)0, ∀ (si,t)Ω̄xN×(0,T).                □

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

D-z(si,t)=z(si,t)-z(si-1)h, D+z(si,t)=z(si+1,t)-z(si,t)h,

respectively, and the second order central finite difference operator as


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

limε0max1iN-1exp(-psi/ε)εm/2=0     andlimε0max1iN-1exp(-p(1-si)/ε)εm/2=0.

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)C4(D̄). Then the semidiscrete solution Zi(t) of Equation (12) satisfies


Proof. The truncation error in spatial direction is considered as

Lh(z(si,t)-Zi(t))=Lhz(si,t)-LhZi(t)=-εd2ds2z(si,t)+Ds+Ds+h2λ2z(si,t)=-εd2ds2z(si,t)+ελ2(h2d2ds2z(si,t)+h412d4ds4z(si,t)).    (14)

Note that we have used Taylor expansions of zi−1(t) and zi+1(t). A truncated Taylor expansion of 1λ2 of order five becomes

1λ2=β24(4β2h2-13+β2h260).    (15)

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

Lh(z(si,t)-Zi(t))=ε12(d4ds4z(si,t)-β2d2ds2z(si,t))h2+ εβ2h4(β2240d2z(si,t)ds2-1144d4z(si,t)ds4)+h6εβ42880d4z(si,t)ds4.    (16)

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


The truncation error at s = sN, become

KN(Z(si)z(si))=KNZ(sN)KNs(si),=ϕrKNZ(sN),=Kz(si)KNZ(sN),=z(sN)ε01g(s)z(s)ds(Z(sN)εs0sNg(s)z(s)ds),=εs0sNg(s)z(s)dsεi=1Ngi1zi1+4gizi+gi+1zi+13h,=ε[s0s1g(s)z(s)ds+s1s2g(s)z(s)ds++sNsN+1g(s)z(s)ds]ε[g0z0+4g1z1+g2z23h++gN1zN1+4gNzN+gN+1zN+13h],|KN(Z(si)z(si))|=|Cε(h4z(4)(ξ1)+h4z(4)(ξ2)++h4z(4)(ξN))|,|KN(Z(si)z(si))|Cεh4(z(4)(ξ1)+z(4)(ξ2)++z(4)(ξN)),Cεh4||d4z(ξi)dx4||Ch2=CN2.                             (17)

Theorem 3. The semidiscrete solutions satisfy the uniform error bound

sup0<ε1maxi|z(si,t)-Zi(t)|D̄CN-2.    (18)

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 = tj+1tj,j = 0(1)M − 1 is constructed on the time domain Dt = [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 zi(tj) as Zij, we construct the time discretization as follows.

Zij-Zij-1Δt=BZij+Fij    (19)

with the initial condition Z0(t) = ϕl(tj), and by rearranging Equation (19), we obtain

Zij=[I+ΔtB]-1[ΔtFij+Zii-1].    (20)

Lemma 8. Suppose |iz(s,tj)ti|C, (s,t)D̄,i=0,1,2. Then the local truncation error associated with the time direction satisfies |ej|C(Δt)2.

Proof. The local truncation error is defined as

ej=z(tj)-Zij    =z(tj)-[I+ΔtB]-1[ΔtFij+Zij-1].

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


However, zt(tj) = F(tj) − B(tj)z(tj). Thus,


Now, the local truncation error ej becomes


Since the matrix B is invertible, using the relation (Δt)2 > (Δt)3 for small Δt and z(tj) ≤ C, we obtain



Lemma 9. The global error estimate in the time direction is given by ||Ej+1|| ≤ CΔt, ∀jTt, where Ej+1=maxi|Zi(tj+1)-Zi,j+1|D.

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

||Ej+1||=||i=1ej||    jT/Δt                  ||e1||+||e2||+||e3||+||e4||++||ej||                  C1(jΔt)Δt                  C1TΔt         since  jΔtT                  CΔt.


||Ej+1||=maxi|Zi(tj+1)-Zij+1|DCΔt.    (21)

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

sup0<ε1maxi|Zi(tj+1)-Zij+1|DCΔt.    (22)


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


where z(s, t) and Zij+1 are the solutions to problems Equations (1) and (12), respectively.

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


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 UNt be a conducted solution of a problem using mesh points N and time step size Δt. Again, Ui,j2N,Δ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εN,Δt=maxi,j|Zi,jN,Δt-Zi,j2N,Δt/2|, and the parameter uniform error estimate by using EN,Δt=maxε(EεN,Δt). We calculate the rate of convergence of the developed scheme by using PεN,Δt=log2(EεN,Δt)-log2(Eε2N,Δt/2). The parameter rate of convergence is calculated as PN,Δt=log2(EN,Δt)-log2(E2N,Δt/2).

Example 1.

{z(s,t)tε2z(s,t)s2+1+s22z(s,t)=et1+sin(π s),  (s,t)(0,1)×(0,1]z(s,0)=0,                                    s(0,1),z(0,t)=0,                                     t(0,1],Kz(1,t)=z(1,t)ε01s6z(s,t)ds=0,   t(0,1].

Example 2.

{z(s,t)tε2z(s,t)s2+1+s22z(s,t)=t3,     (s,t)(0,1)×(0,1]z(s,0)=0,                                          s(0,1),z(0,t)=0,                                           t(0,1],Kz(1,t)=z(1,t)ε01cos(s)z(s,t)ds=0,    t(0,1].

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. 3-D graph of numerical solution for Example (1) which displays the existing layer. (A) ε = 10−2. (B) ε = 10−12.


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


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


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


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


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.


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/

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

Keywords: singularly perturbed problems, partial differential equations, reaction-diffusion, method of lines, uniform convergence, nonlocal boundary condition

Citation: Wondimu GM, Woldaregay MM, Dinka TG and Duressa GF (2022) Numerical treatment of singularly perturbed parabolic partial differential equations with nonlocal boundary condition. Front. Appl. Math. Stat. 8:1005330. doi: 10.3389/fams.2022.1005330

Received: 28 July 2022; Accepted: 03 October 2022;
Published: 25 October 2022.

Edited by:

Lijian Jiang, Tongji University, China

Reviewed by:

Majid Darehmiraki, Behbahan Khatam Alanbia University of Technology, Iran
Khalid Hattaf, Centre Régional des Métiers de l'Education et de la Formation (CRMEF), Morocco

Copyright © 2022 Wondimu, Woldaregay, Dinka and Duressa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Mesfin Mekuria Woldaregay,