# Exponential Fitted Operator Method for Singularly Perturbed Convection-Diffusion Type Problems with Nonlocal Boundary Condition

Authors: Habtamu Garoma Debela
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

This paper presents the study of singularly perturbed differential equations of convection-diffusion type with nonlocal boundary condition. The proposed numerical scheme is a combination of the classical finite difference method for the boundary conditions and exponential fitted operator method for the differential equations at the interior points. Maximum absolute errors and rates of convergence for different values of perturbation parameter and mesh sizes are tabulated for the numerical examples considered. The method is shown to be first-order accuracy independent of the perturbation parameterε.


---

## Body

## 1. Introduction

Singularly perturbed differential equations are typically characterized by the presence of a small positive parameterε multiplying some or all of the highest order terms in differential equations. Such types of problems arise frequently in mathematical models of different areas of physics, chemistry, biology, engineering science, economics, and even sociology. The well-known examples are the heat transfer problem with large Peclet numbers, semiconductor theory, chemical reactor theory, reaction-diffusion process, theory of plates, optimal control, aerodynamics, seismology, oceanography, meteorology, and geophysics. Solutions of such equations usually possess thin boundary or interior layers where the solutions change very rapidly, while away from the layers, the solutions behave regularly and change slowly. More details about these problems can be found in [1–4] and also the literature cited there. Due to the presence of these boundary layers, the usual numerical treatment of singularly perturbed problems gives rise to computational difficulties. Standard numerical methods are not appropriate for practical applications when the perturbation parameter ε is sufficiently small. Therefore, it is necessary to develop suitable numerical methods that are uniformly convergent with respect to ε: To solve these problems, there are generally two types of approaches, such as fitted operator methods that reflect the nature of the solution in the boundary layers and fitted mesh methods which use layer-adapted meshes. In recent years, many authors have worked for solving singularly perturbed problems with one or two boundary layers using uniformly convergent numerical methods [5–11]. Boundary value problems including nonlocal conditions which connect the values of the unknown solution at the boundary with values in the interior are known as nonlocal boundary value problems. The study of this kind of problems was initiated by Il’in and Moiseev in [12, 13], motivated by the work of Bitsadze and Samarskii on nonlocal linear elliptic boundary value problems [14]. These problems have been used to represent mathematical models of a large number of phenomena, such as problems of semiconductors in electronics, the vibrations of a guy wire of a uniform cross-section, heat transfer problems, problems of hydromechanics, catalytic processes in chemistry and biology, the diffusion-drift model of semiconducting devices, and some other physical phenomena [15–17]. The existence and uniqueness of the solutions of nonlocal boundary value problems have been studied by many authors [18, 19]. Some approaches for the numerical solution of singularly perturbed nonlocal boundary value problems have been proposed in [20–28]. Uniformly convergent numerical methods of order second and high for solving different singularly perturbed problems have been studied in [29–34]. To the best of our knowledge, the problem under consideration has not been done using the fitted operator method. Motivated by papers [35, 36], we develop a uniformly convergent numerical method for solving singularly perturbed problem under consideration.


## 2. Statement of the Problem

Consider the following singularly perturbed problem with nonlocal condition of the form(1)Lyx≡εy′′x+axy′x=fx,0<x<l,with the given conditions(2)y′0=Aε,(3)y0+γyl1=Byl+d,where 0<ε≪1 is a small positive parameter, A,B,γ, and d are given constants, l1 and l are given real numbers, and Ω=0,1 and Ω¯=0,l. We assume that ax≥a>0 and fx are sufficiently smooth functions on Ω¯. Under these assumptions, singularly perturbed nonlocal Equations (1)–(3) possess a unique solution indicating a boundary layer of exponential type at x=0.


## 3. Properties of Continuous Solution

The following lemmas [13] are necessary for the existence and uniqueness of the solution and for the problem to be well posed.Lemma 1 (continuous minimum principle).
Assume thatvx∈C2Ω¯ is any function satisfying v0≥0, vl≥0, and Lvx≤0, ∀x∈Ω=0,l. Then, vx≥0,∀x∈Ω=0,l.Proof.
Letx∗ be such that vx∗=minx∈0,lvx and assume that vx∗<0. Clearly, x∗∉0,l; therefore, v′x∗=0 and v′′x∗≥0. Moreover, Lvx∗=εv′′x∗+ax∗v′x∗≥0, which is a contradiction. It follows that vx∗>0, and thus, vx≥0,∀x∈0,l.
The uniqueness of the solution is implied by this minimum principle. Its existence follows trivially (as for linear problems, the uniqueness of the solution implies its existence). This principle is now applied to prove that the solution of Equations (1)–(3) is bounded. The following lemma shows the bound for the derivatives of the solution.Lemma 2.
Letyε be the solution of Pε. Then, for k=0,1,2,3,(4)∣yεkx∣≤C1+ε−kexp−aεx,∀x∈0,lProof.
The homogeneous differential equation of Equation (1) is(5)εy′′x+axy′x=0.
The characteristic equation of Equation (5) is(6)εm2+am=0⇒m=0
or(7)m=−aε.
The asymptotic solution of Equation (5) is given by(8)ux=A+Bexp−aεx,where A and B are arbitrary constant.
To get thekth derivative of the asymptotic solution of the homogeneous part of Equation (5),(9)u′x=Cε−1exp−aεx,u′′x=Cε−2exp−aεx,u′′x=Cε−3exp−aεx.
In general, fork=1,2,3,(10)ukx=Cε−kexp−aεx.
The reduced problem obtained from Equation (1) takes axv0′x=fx, where v00=y0 and has the solution(11)v0x=y0+∫0xftatdtv0x≤y0+∫0xftatdt≤C+fζaζ∫0xdt≤C+fζaζx,x∈0,l,≤C;
from the assumptions ona and f, it is clear that for k=0,1,2,3,(12)v0kx≤C,∀x∈0,1.
So, from the relationyε=v0+u, we have yεk=v0k+uk, and from the relation of triangular inequality,(13)yεk≤v0k+uk≤C+Cε−kexp−aεx≤C1+ε−kexp−aεx.
Therefore, it is well accepted that the solution of Equation (1) has a boundary layer near x=0 and its derivatives satisfy(14)∣yεkx∣≤C1+ε−kexp−aεx,∀x∈0,l.


## 4. Formulation of the Method

Consider the homogeneous differential equation with constant coefficientεy′′x+ay′x=0 whose solution is given by(15)yx=A+Bexp−aεx,whereA and B are constants which will be determined depending on the given conditions. Now, dividing the interval 0,l into N equal parts with constant mesh length h=l/N, we obtain xi=x0+ih, for i=1,2,⋯,N, where x0=0,xN=l.To demonstrate the procedure, we consider Equation (1), at discrete nodes xi(16)εyi′′x+axiyi′x=fxi.Approximating Equation (16) by central difference approximations, we obtain(17)εyi−1−2yi+yi+1h2+aiyi+1−yi−12h=fi.Under the assumption thatfi is bounded, introducing the fitting parameter σ onto the higher order difference approximation of Equation (17), multiply both sides by h, and evaluating its limit gives(18)σ=ρalimh→0yi+1−yi−12limh→0yi+1−2yi+yi−1,where ρ=h/ε.Evaluating Equation (15) at each nodal points xi, we obtain(19)limh→0yi=A1−exp−1εexp−iρ−Bexp−1ε1−exp−iρ,limh→0yi+1=A1−exp−1εexp−iρ−Bexp−1ε1−exp−iρ,limh→0yi−1=A1−exp−1εexp−iρ−Bexp−1ε1−exp−iρ,(20)σ=ρalimh→0yi+1−yi−12limh→0yi+1−2yi+yi−1=aρ2cothaρ2.Hence, from Equations (17) and (20), we get(21)εσh2yi−1−2yi+yi+1+ai2hyi+1−yi−1=fi.This can be rewritten as three-term recurrence relation:(22)Eiyi−1+Fiyi+Giyi+1=Hi,i=1,2,⋯,N−1,where(23)Ei=εσh2−ai2h,Fi=−2εσh2,Gi=εσh2+ai2h,Hi=fi.Since the problem involves nonlocal boundary conditions, we considered the following cases, to obtain two equations at each end conditions.Fori=0, Equation (22) becomes(24)E0Y−1+F0Y0+G0Y1=H0.Here, in Equation (24), the term Y−1 is out of the domain, so that using Equation (2), we have(25)Y′0=μ0ε=Y1−Y−12h⇒Y−1=Y1−2hY′0.Putting Equation (25) into Equation (24) gives(26)E0Y0+E0+G0Y1=H0+2hE0Y′0.Fori=N, Equation (22) becomes(27)ENYN−1+FNYN+GNYN+1=HN.Here, in Equation (27), the term YN+1 is out of the domain, so that using Equation (3), we have(28)yN+1=y0B+γYl1B−dB.Putting Equation (28) into Equation (27) gives(29)GNBY0+GNγBYl1=HN+GNdB.Therefore, Equation (1) with the given boundary conditions in Equations (2) and (3) can be solved using the schemes in Equations (22), (26), and (29) which gives N×N system of algebraic equations.


## 5. Uniform Convergence Analysis

In this section, we need to show the discrete scheme in Equation (22) and satisfy the discrete minimum principle, uniform stability estimates, and uniform convergence.Lemma 3 (discrete minimum principle).
LetYi be any mesh function that satisfies Y0≥0, YN≥0, and LεNYi≤0, i=1,2,3,⋯,N−1, then Yi≥0, for i=0,1,2,⋯,N.Proof.
The proof is by contradiction. Letj be such that Yj=miniYi and suppose that Yj≤0. Clearly, j∉0,N. Yj+1−Yj≥0 and Yj−Yj−1≤0.
Therefore,(30)LεNYj=εYj+1−2Yj+Yj−1Ψi2+ajYj+1−Yjh=εΨi2Yj+1−2Yj+Yj−1+ajhYj+1−Yj=εΨi2Yj+1−Yj−Yj−Yj−1+ajhYj+1−Yj≥0,where the strict inequality holds if Yj+1−Yj>0. This is a contradiction, and therefore, Yj≥0. Since j is arbitrary, we have Yi≥0,fori=0,1,2,⋯,N.
From the discrete minimum principle, we obtain anε− uniform stability property for the operator LεN.Lemma 4 (uniform stability estimate).
Ifϕj is any mesh function such that(31)ϕ0=ϕN=0,thenϕj≤1amax1≤i≤N−1LεNϕi,j=0,1,2,⋯,N.Proof.
We introduce two mesh functionsψj+,ψj− defined by(32)ψj±=1amax1≤i≤N−1∣LεNϕi∣±ϕj.It follows that(33)ψ±0=1amax1≤i≤N−1∣LεNϕi∣±ϕ0=1amax1≤i≤N−1∣εδ2ϕi+aiD0ϕi∣±ϕ0=1amax1≤i≤N−1∣εδ2ϕi+aiD0ϕi∣≥0,ψ±N=1amax1≤i≤N−1LεNϕi±ϕN=1amax1≤i≤N−1∣εδ2ϕi+aiD0ϕi∣±ϕN=1amax1≤i≤N−1∣εδ2ϕi+aiD0ϕi∣≥0,for all j=1,2,⋯,N−1,(34)LεNψj±=1amax1≤i≤N−1∣LεNϕi∣±LεNϕj≤0.
From a discrete minimum principle, ifψ0≥0,ψN≥0, and LεNψj≤0,∀0<j<N, then ψj±≥0,0≤j≤N.We provide above that the discrete operatorLεN satisfy the minimum principle. Next, we analyze the uniform convergence analysis.Theorem 5.
Letyxi and Yi be, respectively, the exact solution of Equations (1)–(3) and numerical solutions of Equation (22). Then, for sufficiently large N, the following parameter uniform error estimate holds:(35)LNyi−Yi≤Ch1+supx∈0,1exp−axi/εε3.Proof.
Let us consider the local truncation error defined as(36)LNyxi−Yi=εσρy′′xi−D+D−yxi+axiy′xi−D0yxi=ερa02cothρa02−1D+D−yxi+εy′′xi−D+D−yxi+axiy′xi−D0yxi,where σρ=a0ρ/2cotha0ρ/2 and ρ=N−1/ε.
Now, forz>0, C1 and C2 are constants, and we have ∣zcothz−1∣≤C1z2,z≤1. Similarly, for z→∞, since limz→∞cothz=1, ∣zcothz−1∣≤C1z is given.
In general, for allz>0, as [11], we write(37)C1z2z+1≤zcothz−1≤C2z2z+1,implying that(38)εa0ρ2cotha0ρ2−1≤εN−1/ε2N−1/ε+1=N−2N−1+ε.
Using Taylor series expansion, the bound foryxi−1 and yxi+1 at xi as(39)yxi−1=yxi−hy′xi+h22!y′′xi−h33!y3xi+h44!y4xi+Oh5,yxi+1=yxi+hy′xi+h22!y′′xi+h33!y3xi+h44!y4xi+Oh5.
We obtain the bound for(40)D+D−yxi≤Cy′′xi,∣y′′xi−D+D−yxi∣≤CN−2y4xi.
Similarly, for the first derivative term,(41)y′xi−D0yxi≤CN−2y3xi,where ∣ykxi∣=supxi∈x0,xN∣ykxi∣, k=2,3,4.
Using the bounds in Equations (40) and (41), we obtain(42)∣LNyxi−Yi∣≤CN−2N−1+ε∣y′′xi∣+εCN−2∣y4xi∣+CN−2∣y3xi∣≤CN−2N−1+ε∣y′′xi∣+CN−2ε∣y4xi∣+∣y3xi∣.
Now, using the bounds for the derivatives of the solution in Lemma (2), we have(43)∣LNyxi−Yi∣≤CN−2N−1+ε1+ε−2exp−cxiε+CN−2ε1+ε−4exp−cxiε+1+ε−3exp−cxiε≤CN−2N−1+ε1+ε−2exp−cxiε+CN−2ε+ε−3exp−cxiε+1+ε−3exp−cxiε,which simplifies to(44)∣LNyxi−Yi∣≤CN−2N−1+ε1+ε−3exp−c1xiε,since ε−3≥ε−2.Lemma 6.
For a fixed mesh and forε→0, it holds(45)limε→0max1≤i≤N−1exp−axi/εεm=0,m=1,2,3,⋯,limε→0max1≤i≤N−1exp−a1−xi/εεm=0,m=1,2,3,⋯.Proof.
Consider the partition0,1≔0=x0<x1<⋯<xN−1<xN=1 for the interior grid points, we have(46)max1≤i≤N−1exp−axi/εεm≤exp−ax1/εεm=exp−ah/εεm,max1≤i≤N−1exp−a1−xi/εεm≤exp−a1−xN−1/εεm=exp−ah/εεm,asx1=1−xN−1=h.
The repeated application of L’Hospital’s rule gives(47)limε→0exp−ah/εεm=limr=1/ε→∞rmexpahr=limr=1/ε→∞m!ahmexpahr=0.
This completes the proof.Theorem 7.
Letyxi and Yi be the exact solution of Equations (1)–(3) and numerical solutions of Equation (22), respectively. Then, the following error bound holds(48)sup0<ε<<1∣yxi−Yi∣≤CN−2N−1+ε≤CN−1.Proof.
By substituting the results in Lemma6 into Theorem 5 and applying the discrete maximum principle, we obtain the required bound.
For the caseε>N−1, the scheme secures second-order convergence, and we expect to lose an order of convergence for ε≤N−1, and in fact, it turns out that the scheme is first-order uniformly convergent.


## 6. Numerical Example and Results

To validate the established theoretical results, we perform numerical experiments using the model problems of the form in Equations (1)–(3).Example 8.
Consider the model singularly perturbed boundary value problem:(49)εy′′x+2y′x=ε−2e−x,0<x<1,subject to the boundary conditions(50)y′0=1ε,y0+13y14+y1=1.
Its exact solution is(51)yx=d1+d2e−2x/ε+e−x,where d1=−3/7e−1+1/3e−1/4+1+e−2x/ε+1/3e−1/2εd2, d2=−1+ε/2.Example 9.
Consider the model singularly perturbed boundary value problem:(52)εy′′x+2y′x=ε−2e−x,0<x<1,subject to the boundary conditions(53)y′0=1ε,y0+23y34+y1=1.
Its exact solution is(54)yx=d1+d2e−2x/ε+e−x,where d1=−3/8e−1+2/3e−3/4+1+e−2x/ε+2/3e−3/2εd2 and d2=−1+ε/2.
We define the pointwise absolute errorsEεN and the computed ε-uniform maximum pointwise error EN as follows:(55)EεN=Y−y,EN=maxεEεN,where Y is the numerical approximation to y for various values of N and ε. We also define the computed ε-uniform convergence rate:(56)RN=log2ENE2N.


## 7. Discussion and Conclusion

This study introduces a uniformly convergent numerical method based on the exponential fitted operator method for solving singularly perturbed boundary value problems with nonlocal boundary condition. The behavior of the continuous solution of the problem is studied and shows that it satisfies the continuous stability estimate, and the derivatives of the solution are also bounded. The numerical scheme is developed on uniform mesh. The nonlocal boundary condition is treated using the finite difference formula, and the results are compared accordingly. The stability of the developed scheme is established, and its uniform convergence is proved. To validate the applicability of the method, two model problems are considered for numerical experimentation for different values of the perturbation parameter and mesh points. The numerical results are tabulated in terms of maximum absolute errors, numerical rate of convergence, and uniform errors (see Tables1–4) and compared with the results of the previously developed numerical methods existing in the literature (Tables 2 and 4). Further, the behavior of the numerical solutions is plotted in Figure 1; Figure 2 shows that as the mesh size decreases or as the number of nodal points increases, point wise absolute error decreases; and Figure 3 shows that our method is ε-uniformly convergent independent of the mesh size for h≥ε, where the classical numerical methods fail in the log-log scale. Unlike other fitted operator finite difference methods constructed in standard ways, the method that we presented in this paper is fairly simple to construct. Moreover, the method is more accurate and gives a good result where existing numerical methods fail (that is for the values where the perturbation parameter, ε, is much less than the mesh size, h).Table 1
Maximum pointwise absolute errors, computedε-uniform errors, and convergence rates for Example 8.
εN=16N=32N=64N=128N=25610−44.5340e−032.2114e−031.0919e−045.4256e−042.7046e−0410−84.5342e−032.2115e−031.0920e−035.4258e−042.7044e−0410−124.5342e−032.2115e−031.0920e−035.4258e−042.7044e−0410−164.5342e−032.2115e−031.0920e−035.4258e−042.7044e−0410−204.5342e−032.2115e−031.0920e−035.4258e−042.7044e−04EN4.5342e−032.2115e−031.0920e−035.4258e−042.7044e−04RN1.03581.01811.00911.0045Table 2
Comparison ofε-uniform maximum absolute errors and ε-uniform rate of convergence for Example 8.
εN=16N=32N=64N=128N=256Present methodEN4.5342e−032.2115e−031.0920e−035.4258e−032.7044e−04RN1.03581.01811.00911.0045Method in [36]EN0.01040760.00517700.00258180.00128920.000644RN1.011.001.001.00Table 3
Maximum pointwise absolute errors, computedε-uniform errors, and convergence rates for Example 9.
εN=16N=32N=64N=128N=25610−44.1080e−032.0018e−039.8802e−044.9081e−042.4426e−0410−84.1082e−032.0019e−039.8807e−044.9083e−042.4463e−0410−124.1082e−032.0019e−039.8807e−044.9083e−042.4462e−0410−164.1082e−032.0019e−039.8807e−044.9083e−042.4462e−0410−204.1082e−032.0019e−039.8807e−044.9083e−042.4462e−04EN4.1082e−032.0019e−039.8807e−044.9083e−042.4462e−04RN1.03711.01871.00941.0047Table 4
Comparison ofε-uniform maximum absolute errors and ε-uniform rate of convergence for Example 9.
εN=16N=32N=64N=128N=256Present methodEN4.1082e−032.0019e−039.8807e−044.9083e−042.4462e−04RN1.03711.01871.00941.0047Method in [36]EN0.01164990.00579490.00289000.00144310.0007211RN1.011.001.001.00Figure 1
The behavior of the numerical solution atε=2−4 and N=32 of Examples 8 and 9, respectively.
(a)(b)Figure 2
Pointwise absolute error of Examples8 and 9 at ε=10−20 with different mesh points N, respectively.
(a)(b)Figure 3
ε-uniform convergence with exponential fitted operator method in log-log scale for Examples 8 and 9, respectively.
(a)(b)



---

*Source: 5559486-2021-03-25.xml*