# Accelerated Exponentially Fitted Operator Method for Singularly Perturbed Problems with Integral Boundary Condition

Authors: Habtamu Garoma Debela; Gemechis File Duressa
Journal: International Journal of Differential Equations
 (2020)

---

## Abstract

In this paper, we consider a class of singularly perturbed differential equations of convection diffusion type with integral boundary condition. An accelerated uniformly convergent numerical method is constructed via exponentially fitted operator method using Richardson extrapolation techniques and numerical integration methods to solve the problem. The integral boundary condition is treated using numerical integration techniques. Maximum absolute errors and rates of convergence for different values of perturbation parameter and mesh sizes are tabulated for the numerical example considered. The method is shown to beε-uniformly convergent.


---

## Body

## 1. Introduction

Boundary value problems involving integral boundary conditions have received considerable attention in recent years [1, 2]. For a discussion of existence and uniqueness results and for applications of problems with integral boundary conditions, one can refer [3–5] and the references therein. In [2, 6, 7], some approximating or numerical treatment aspects of this kind of problems have been considered. However, the methods or algorithms developed so far mainly concerned with the regular cases (i.e., when the boundary layers are absent). Boundary value problems with integral boundary conditions in which the leading derivative term is multiplied by a small parameter are called singularly perturbed problems with integral boundary conditions. The solutions of such types of problems manifest boundary layer phenomena where the solution changed abruptly. As a result, numerical analysis of singular perturbation cases has been far from trivial because of the boundary layer behavior of the solution. The solutions of the problems with boundary layer undergo rapid changes within very thin layers near the boundary or inside the problem domain [2, 8–10], and hence classical numerical methods for solving such problems are unstable and fail to give good results when the perturbation parameter is small (i.e., for h≥ε) [10]. Therefore, it is important to develop a numerical method that gives good results for small values of the perturbation parameter where others fails to give good result and convergent independent of the values of the perturbation parameter. In recent years, the authors [11–15] have developed various numerical schemes on uniform meshes for singularly perturbed first and second order differential equations with integral boundary conditions.As far as the researchers’ knowledge is concerned numerical solution of the singularly perturbed boundary value problem containing integral boundary condition via the accelerated exponential fitted operator method is first being considered. Hence, this paper proposed a uniformly convergent numerical method based on exponential fitted operator and numerical integration methods to solve the problem under consideration.


## 2. Statement of the Problem

Consider the following singularly perturbed problem with integral boundary condition:(1)Ly≡ε2y″x+εaxy′x−bxyx=fx,0<x<l,(2)y0=μ0,(3)L1y=yl−∫l0l1gxyxdx=μ1,0≤l0<l1≤l,where ε,0<ε≪1 is a perturbation parameter, μ0 and μ1 are given constants, and the functions ax≥α>0,bx≥β>0,gx, and fx are sufficiently smooth functions in the 0,l. The solution yx of problems (1)–(3) has in general boundary layers at x=0 and x=l for ε near 0.In this paper, we analyze the exponential fitted operator scheme with numerical integration techniques on a uniform mesh for the numerical solution of equations (1)–(3). Uniform convergence is proved in the discrete maximum norm. Finally, we formulate the algorithm for solving the discrete problem and give the illustrative numerical results.


## 3. Properties of Continuous Solution

The differential operator for equation (1) is given by(4)Lε≡ε2d2dx2+εaddx−b,and it satisfies the following minimum principle for boundary value problems (BVPs). The following lemmas [9] 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 is any sufficiently smooth function satisfying v0≥0 and vl≥0. Then, Lvx≤0 for all x∈Ω=0,l implies that vx>0 for all 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∗=ε2v″x∗+εax∗v′x∗−bx∗vx∗≥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).
The following lemma shows the bound derivatives of the solution fork=0,1,2,3,4.Lemma 2.
Letyε be the solution of Pε. Then, for k=0,1,2,3,4,(5)yεkx≤C1+ε−kexp−c0xε+exp−c1l−xε,for allx∈Ω¯,where c0=1/2a20+4b0+a0 and c1=1/2a2l+4bl−al.Proof.
For the proof refer [2].


## 4. Formulation of the Method

For small values ofε, the solution yx of problems (1)–(3) has in general boundary layers at x=0 and x=l (see [2]).The linear ordinary differential equation (1) cannot, in general, be solved analytically because of the dependence of axandbx on the spatial coordinate x. We divide the interval 0,1 into N equal parts with constant mesh length h. Let 0=x0,x1,…,xN=1 be the mesh points. Then, we have xi=ih,i=0,1,2,…,N. If we consider the interval x∈0,1 and the coefficients of equation (1) are evaluated at the midpoint of each interval, then we will obtain the differential equation:(6)ε2y″x+εaxy′x−bxyx=fx,x∈Ω=0,1,y0=y0=μ0.In order to evaluate the fitting factor, divide both side of equation (6) by ε, and we obtain(7)εy″x+axy′x+pxyx=rx,where px=−bx/εandrx=fx/ε.To find the numerical solution of equation (7), we use the theory applied in the asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer is in the left side of the domain, i.e., near x=0. From the theory of singular perturbations given by O’Malley [16] and using Taylor’s series expansion for ax about x=0 and restriction to their first terms, we get the asymptotic solution as(8)yx=y0x+μ0−y00exp−a0εx+Oε,where y0x is the solution of the reduced problem (obtained by setting ε=0) of equation (7) which is given by(9)axy0′x+pxy0x=rx,withy01=θ1,where θ1=μ1+∫l0l1gxyxdx.Consideringh small enough, the discretized form of equation (8) becomes(10)yih=y00+μ0−y00exp−a0εih,which simplifies to(11)yih=y00+μ0−y00exp−iρa0,where ρ=h/εandh=1/N.To handle the effect of the perturbation parameter artificial viscosity, exponentially fitting factorσρ is multiplied on the term containing the perturbation parameter as(12)εσρy″x+axy′x+pxyx=rx,with boundary conditions y00=μ0andyN=θ1.Next, we consider the difference approximation of equation (7) on a uniform grid Ω¯N=xii=0N and denote h=xi+1−xi.For any mesh functionzi, define the following difference operators:(13)D+zi=zi+1−zih,D−zi=zi−zi−1h,D0zi=zi+1−zi−12h,D+D−zi=zi+1−2zi+zi−1h2.By applying the central finite difference scheme on equation (12) takes the form(14)εσρD+D−yxi+axiD0yxi+pxiyxi=rxi,with the boundary conditions y00=μ0andyN=θ1.Using operator, equation (7) is rewritten as(15)Lhyi=ri,with the boundary conditions y0=μ0andyN=θ1, where(16)Lhyi=εσρyi+1−2yi+yi−1h2+axiyi+1−yi−12h+pxiyi=ri.Multiplying equation (16) by h and considering h small and truncating the term hri−pxiyxi, results(17)εσρρyi+1−2yi+yi−1+axi2yi+1−yi−1=0.Now using Taylor’s series foryi and yi+1 up to first term and substituting the results in equation (17) into equation (14) and simplifying, the exponential fitting factor is obtained as(18)σρ=ρa02cothρa02.Similarly, by following the above mentioned procedure, the fitting factor for the right layer becomes(19)σρ=ρa12cothρa12.Assume thatΩ¯N denote partition of 0,1 into N subintervals such that 0=x0<x1<⋯<xN=1 with xi=ih,h=1/N,i=0,1,2,…,N.Case 1.
For the left layer, consider equation (6) on the domain Ω=0,1 which is given by(20)ε2y″x+εaxy′x−bxyx=fx.
Hence, the required finite difference scheme becomes(21)ε2σρh2yi−1+−2ε2σρh2−εaxih−bxiyi+ε2σρh2+εaxi2hyi+1=fi,for i=0,1,2,…,N.
The numerical scheme in equation (21) can be written in three-term recurrence relation as(22)Eiyi−1+Fiyi+Giyi+1=Hi,i=1,2,…,N,where Ei=ε2σ/h2,Fi=−2ε2σ/h2−εai/h−bi,Gi=ε2σ/h2+εai/2h,andHi=fi.Case 2.
For the right layer, consider equation (6) on the domain Ω=0,1 which is given by(23)ε2y″x+εaxy′x−bxyx=fx.
Hence, the required finite difference scheme becomes(24)ε2σρh2−εaxihyi−1+−2ε2σρh2+εaxih−bxiyi+ε2σρh2yi+1=fi,for i=0,1,2,…,N.
The numerical scheme in equation (24) can be written in three-term recurrence relation as(25)Eiyi−1+Fiyi+Giyi+1=Hi,i=1,2,…,N,where Ei=ε2σ/h2−εai/h,Fi=−2ε2σ/h2+εai/h−bi,Gi=ε2σ/h2,andHi=fi,
Fori=N, we considered Simpson’s rule, i.e., treating the integral boundary condition using Simpson’s rules. Suppose gxyx is a function defined on the interval 0,l, and let xii=02N be a uniform partition of 0,l with step length h.
The composite Simpson’s rule approximates the integral ofgxyx is given by(26)∫l0l1gxyxdx=h3gl0yl0+gl1yl1+2∑i=12N−1gx2iyx2i+4∑i=12Ngx2i−1yx2i−1.
Substituting equation (26) into equation (3) gives(27)yl−h3gl0yl0+gl1yl1+2∑i=12N−1gx2iyx2i+4∑i=12Ngx2i−1yx2i−1=μ1.
Sincey0=μ0, from equations (2) and (27), we obtain:(28)−4h3∑i=12Ngx2i−1yx2i−1−2h3∑i=12N−1gx2iyx2i+yl−h3gl1yl1=μ1+h3gl0yl0.
Therefore, the problem in equation (1) with given boundary condition in equations (2) and (3) can be solved using the scheme equations (22), (25), and (28) which forms N×N system of algebraic equations.


## 5. Uniform Convergence Analysis

In this section, we need to show the discrete scheme in equations (22), (25), and (28) satisfy the discrete minimum principle and uniform convergence.Lemma 3 (discrete minimum principle).
Letvi be any mesh function that satisfies v0≥0,vN≥0, and Lhvi≤0,i=1,2,…,N−1, then vi≥0,i=0,1,2,…,N.Proof.
The proof is obtained by contradiction. Letj be such that vj=minivi and suppose that vj<0. Clearly, j∉0,N,vj+1−vj≥0 and vj−vj−1≤0.
Therefore,(29)Lhvj=εh2vj+1−2vj+vj−1+ajhvj+1−vj+pjvj=εh2vj+1−vj−vj−vj−1+ajhvj+1−vj+pjvj≥0,where the strict inequality holds if vj+1−vj>0. This is a contradiction, and therefore vj≥0. Since j is arbitrary, we have vi≥0,i=0,1,2,…,N.
We proved above the discrete operatorLh satisfy the minimum principle. Next, we analyze the uniform convergence analysis. Let us define the forward, backward, and second order finite difference operators as follows:(30)D+vj=vj+1−vjh,D−vj=vj−vj−1h,δ2vj=D+D−vj=D+vj−D−vjh.Theorem 1.
Letyxi and yi be, respectively, the exact solution of (1)–(3) and numerical solutions of equation (15). Then, for sufficiently large N,the following parameter uniform error estimate holds:(31)sup0<ε≤1yxi−yi≤CN−2.Proof.
Let us consider the local truncation error defined as(32)Lhyxi−yi=εσρd2dx2−D+D−yxi+axiddx−D+yxi,where εσρ=a0N−1/2cotha0N−1/2ε since ρ=N−1/ε. In our assumption ε≤h=N−1.
By consideringN is fixed and taking the limit for ε⟶0, we obtain the following:(33)limε⟶0εσρ=limε⟶0a0N−12cotha0N−12ε=CN−1.
From Taylor series expansion, the bound for the difference becomes(34)d2dx2−D+D−yxi≤CN−3d4yxidx4,ddx−D+yxi≤CN−2d3yxidx3,where dkyxi/dxk=supxi∈x0,xNdkyxi/dxk,k=3,4.
Now using the bounds and the assumptionε≤N−1, equation (32) reduces to(35)Lhyxi−yi=εσρd2dx2−D+D−yxi+axiddx−D+yxi≤εσρd2dx2−D+D−yxi+axiddx−D+yxi≤CN−3d4yxidx4+CN−2d3yxidx3.
Here, the target is to show the scheme convergence independent on the number of mesh points.
By using the bounds for the derivatives of the solution in Lemma2, we obtain(36)Lhyxi−yi≤CN−3d4yxidx4+CN−2d3yxidx3≤CN−31+ε−4exp−c0xε+exp−c1l−xε+CN−21+ε−3exp−c0xε+exp−c1l−xε≤CN−21+ε−4exp−c0xε+exp−c1l−xε,sinceε−4≥ε−3.
Most of the time during analysis, one encounters with exponential terms involvingε divided by the power function in ε which are always the main cause of worry. For their careful consideration while proving the ε-uniform convergence, we prove for the right layer case as follows.Lemma 4.
For a fixed mesh and forε⟶0, it holds:(37)limε⟶0max1≤j≤N−1exp−c0x/εε3+exp−c1l−x/εε3=0,m=1,2,3,…,where xj=jh,h=1/N,∀j=1,2,…,N−1.Proof.
Consider the partition0,1:=0=x0<x1<⋯<xN−1<xN=1 for the interior grid points, we have(38)max1≤j≤N−1exp−c0xj/εεm≤exp−c0x1/εεm=exp−c0h/εεm,max1≤j≤N−1exp−c11−xj/εεm≤exp−c11−xN−1/εεm=exp−c1h/εεm,asx1=h,1−xN−1=1−N−1h=h.
The repeated application of L’Hospital’s rule gives(39)limε⟶0exp−ch/εεm=limσ=1/ε⟶∞σmexpchσ=limσ=1/ε⟶∞m!chmexpchσ=0.
This complete the proof.Theorem 2.
Under the hypothesis of boundness of discrete solution (i.e., it satisfies the discrete minimum principle), Lemma4 and Theorem 1, the discrete solution satisfy the following bound:(40)sup0≤ε≤1maxiyi−Yi≤CN−1.Proof.
Results from boundness of solution, Lemma4 and Theorem 1, give the required estimates.


### 5.1. Richardson Extrapolation

This technique is acceleration technique which involves combination of two computed approximations of a solution. The combination goes out to be an improved approximation. From the local truncation term, we have(41)yxi−yi≤Ch,where yxiandyi are exact and approximate solutions, respectively, and C is constant free from mesh size h.LetΩ4N be the mesh found by dividing each mesh interval in Ω2N and symbolize the calculation of the solution on Ω4N by y¯i. Consider equation (41) works for any h≠0, which implies(42)yxi−yi≤Ch+R2N,xi∈Ω2N.So that it works for anyh/2≠0 yielding(43)yxi−y¯i≤Ch2+R4N,xi∈Ω4N,where the remainders R2N and R4N are Oh2. Combination of inequalities in equations (42) and (43) leads to yxi−2y¯i−yi≈Oh2 which proposes that(44)yiext=2y¯i−yi,is also a rough calculation of yxi. By means of this approximation to estimate the truncation error, we obtain(45)yxi−yiext≤Ch2,where C is free of mesh size h. Thus, using Richardson extrapolation first-order convergent method is accelerated into second-order convergent as provided in (45). Thus, we can say that the proposed method is second-order convergent.


## 5.1. Richardson Extrapolation

This technique is acceleration technique which involves combination of two computed approximations of a solution. The combination goes out to be an improved approximation. From the local truncation term, we have(41)yxi−yi≤Ch,where yxiandyi are exact and approximate solutions, respectively, and C is constant free from mesh size h.LetΩ4N be the mesh found by dividing each mesh interval in Ω2N and symbolize the calculation of the solution on Ω4N by y¯i. Consider equation (41) works for any h≠0, which implies(42)yxi−yi≤Ch+R2N,xi∈Ω2N.So that it works for anyh/2≠0 yielding(43)yxi−y¯i≤Ch2+R4N,xi∈Ω4N,where the remainders R2N and R4N are Oh2. Combination of inequalities in equations (42) and (43) leads to yxi−2y¯i−yi≈Oh2 which proposes that(44)yiext=2y¯i−yi,is also a rough calculation of yxi. By means of this approximation to estimate the truncation error, we obtain(45)yxi−yiext≤Ch2,where C is free of mesh size h. Thus, using Richardson extrapolation first-order convergent method is accelerated into second-order convergent as provided in (45). Thus, we can say that the proposed method is second-order convergent.


## 6. Numerical Example and Results

To validate the established theoretical results, we perform numerical experiments using the model problem of the form in equations (1)–(3).Example 1.
Consider the model singularly perturbed boundary value problem:(46)ε2y″+2εx+1y′−−x2−2x+4y=−ex+ex2,which subject to the conditions(47)y0=1,y1−∫0.51cosxyxdx=0.
We haveyj≡yjh (the approx. solution obtained via the present method) for different values of h and ε, since the exact solution is not available; the maximum errors (denoted by Eεh) are evaluated using the formula given by the double mesh principle [17]:(48)Eεh:=max0≤j≤Nyjh−y2j2h.
Furthermore, we will tabulate the errors:(49)Eh=max0<ε≤1Eεh.
The numerical rates of convergence are computed using the following formula [17]:(50)rεh:=log2Eεh−log2Eεh/2,and the numerical rate of “ε-uniform convergence” is computed using(51)Rh=log2Eh−log2Eh/2.


## 7. Discussion and Conclusion

This study introduces the uniformly convergent numerical method based on the exponential fitted operator method for solving singularly perturbed boundary value problems with integral boundary conditions. The behavior of the continuous solution of the problem is studied and the derivatives of the solution are also bounded. The numerical scheme is developed on a uniform mesh. The integral boundary condition is treated using numerical integration techniques, namely, Simpson’s rule; 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, a model problem/example is 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 and 2) and compared with the results of the previously developed numerical methods existing in the literature (Table 2). Furthermore, the ε-uniform convergence of the method is shown by the log-log plot of the ε-uniform error (Figure 1) and the numerical solution for various values of N and ε are given (see Figures 2–4). Unlike other fitted operator finite difference methods constructed in standard ways, the method that we presented in this paper is fairly simple to construct.Table 1
Maximum absolute errors and rates of convergence for our example.
ε
N
=
32
N
=
64
N
=
128
N
=
256
N
=
512
10
−
4
6.2800e − 04
1.5350e − 04
3.7899e − 05
9.4122e − 06
2.3501e − 06
2.0325
2.0180
2.0096
2.0049
10
−
8
6.3295e − 04
1.5487e − 04
3.8258e − 05
9.5046e − 06
2.3685e − 06
2.0310
2.0172
2.0091
2.0047
10
−
12
6.3295e − 04
1.5487e − 04
3.8258e − 05
9.5046e − 06
2.3685e − 06
2.0310
2.0172
2.0091
2.0047
10
−
16
6.3295e − 04
1.5487e − 04
3.8258e − 05
9.5046e − 06
2.3685e − 06
2.0310
2.0172
2.0091
2.0047
10
−
20
6.3295e − 04
1.5487e − 04
3.8258e − 05
9.5046e − 06
2.3685e − 06
2.0310
2.0172
2.0091
2.0047
E
h
6.3295e − 04
1.5487e − 04
3.8258e − 05
9.5046e − 06
2.3685e − 06
R
h
2.0310
2.0172
2.0091
2.0047Table 2
Comparison of maximum absolute errors and rates of convergence of our example.
ε
N
=
16
N
=
32
N
=
16
N
=
32
Present method
Method in [2]
10
−
3
2.4362e − 03
5.7378e − 04
1.5304e − 02
7.8100e − 03
2.0861
0.97
10
−
5
2.6259e − 03
6.3247e − 04
1.5304e − 02
7.8810e − 03
2.0537
0.96
10
−
7
2.6276e − 03
6.3295e − 04
1.5304e − 02
7.8810e − 03
2.0536
0.96
E
h
2.6276e − 03
6.3295e − 04
1.5304e − 02
7.8810e − 03
R
h
2.0536
0.96Figure 1
ε-uniform convergence with fitted operator FDM in log-log scale for our example.Figure 2
The numerical solution for various values ofN and ε=10−3.Figure 3
The numerical solution for various values ofN and ε=10−5.Figure 4
The numerical solution for various values ofN and ε=10−7.



---

*Source: 9268181-2020-03-16.xml*