# Uniform Hybrid Difference Scheme for Singularly Perturbed Differential-Difference Turning Point Problems Exhibiting Boundary Layers

Authors: Wondwosen Gebeyaw Melesse; Awoke Andargie Tiruneh; Getachew Adamu Derese
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

In this paper, a class of linear second-order singularly perturbed differential-difference turning point problems with mixed shifts exhibiting two exponential boundary layers is considered. For the numerical treatment of these problems, first we employ a second-order Taylor’s series approximation on the terms containing shift parameters and obtain a modified singularly perturbed problem which approximates the original problem. Then a hybrid finite difference scheme on an appropriate piecewise-uniform Shishkin mesh is constructed to discretize the modified problem. Further, we proved that the method is almost second-orderɛ-uniformly convergent in the maximum norm. Numerical experiments are considered to illustrate the theoretical results. In addition, the effect of the shift parameters on the layer behavior of the solution is also examined.


---

## Body

## 1. Introduction

Many real-life phenomena in different fields of science are modeled mathematically by delay differential or differential-difference equations (DDEs). Equations of this type arise widely in scientific fields such as physics, biosciences, ecology, control theory, economics, material science, medicine, and robotics, in which time evolution depends not only on present states but also on the states at or near a given time in the past. DDEs are also prominent in describing several aspects of infectious disease dynamics such as primary infection, drug therapy, and immune response. In addition, statistical analysis of ecological data has shown that there is evidence of delay effects in the population dynamics of many species; for the detail theory of DDEs, one can refer the books [1, 2].If we restrict the class of DDEs in which the highest derivative is multiplied by a small parameter, then it is said to be a singularly perturbed differential-difference equations (SPDDEs) [3]. In the past, less attention had been given for the solutions of SPDDEs. However, in recent years, there has been a growing interest in the treatment of such problems. This is due to their importance in the modeling of processes in various fields such as optical bistable devices [4], variational problems in control theory [5, 6], the hydrodynamics of liquid helium [7], the first exit-time problem [8], describing the human pupil-light reflex [9], microscale heat transfer [10], and a variety of models for physiological processes or diseases [11, 12].The study of different classes of SPDDEs was initiated by Lange and Miura [8, 13, 14], where they used extension of the method of matched asymptotic expansions for approximating the solution. But in all the cases, they excluded the occurrence of turning points and left it for future study. On the other hand, Kadalbajoo and Sharma [3, 15, 16] initiated the numerical study of SPDDEs with mixed shifts by constructing a variety of numerical schemes. In recent years, different scholars further developed numerical schemes for SPDDEs with mixed shifts, to mention few [17–21]. Most of the works developed so far focuses only on SPDDEs without turning points. In contrast, there are few works on singularly perturbed delay turning point problems. The papers by Ria and Sharma [22–25] are the first and also the only notable works in the treatment of such problems when the solutions exhibit both interior and boundary layers, where the authors used a fitted mesh and fitted operator methods and obtained an almost first-order uniform convergence. Therefore, it is natural to develop a robust numerical method for turning point problems having a better accuracy and efficiency.In this paper, we consider the following second-order linear singularly perturbed differential-difference problem containing mixed shifts and with a turning point atx = 0:(1)−εy″x−axy′x+bxyx−cxyx−δ−dxyx+η=fx,∀x∈Ω=−1,1,(2)yx=ϕx,−1−δ≤x≤−1,yx=ψx,1≤x≤1+η,where a(x), b(x), c(x), d(x), f(x), ϕ(x), and ψ(x) are sufficiently smooth functions on Ω = (−1, 1), 0 < ɛ ≪ 1 is the singular perturbation parameter, and 0 < δ ≪ 1 and 0 < η ≪ 1 are the delay and advance parameters, respectively, together with the following assumptions:(3)a0=0,a′0<0,(4)bx≥b0>0,∀x∈Ω¯=−1,1,(5)bx−cx−dx≥β>0,cx≥K1>0,dx≥K2>0,∀x∈Ω¯,(6)a′x≥a′02,∀x∈Ω¯,under these assumptions, problems (1) and (2) possesses a unique solution having boundary layers of exponential type at x = ± 1, i.e., at both end points [25].Here for the numerical treatment of (1) and (2), we propose an hybrid finite difference scheme on an appropriate piecewise-uniform Shishkin mesh and analyze the stability and uniform convergence the proposed method. Further, we investigate the effect of the shift parameters on the behavior of the solution.Throughout this paper,M (sometimes subscripted) denotes a generic positive constant independent of the singular perturbation parameter ɛ and in the case of discrete problems, also independent of the mesh parameter N. The maximum norm (i.e., f=max−1≤x≤1fx) is used for studying the convergence of the approximate solution to the exact solution of the problem.


## 2. The Continuous Problem

Using Taylor’s series expansion to approximate the terms containing shift arguments gives us(7)yx−δ≈yx−δy′x+δ22y″x,yx+η≈yx+ηy′x+η22y″x.Substituting (7) into (1) and (2) and simplifying gives the following asymptotically equivalent two-point boundary value problem(8)−ε+δ22cx+η22dxy″−ax−δcx+ηdxy′+bx−cx−dxy≈fx,y−1≈ϕ−1,y1≈ψ1.Since (8) is an approximation version of (1) and (2), it is better to use different notation (say u(x)) for the solution of this approximate equation. Thus, (8) can be written as(9)Lu≡−Cεxu″x−Axu′x+Bxux=fx,u−1=ϕ−1=ϕ,u1=ψ1=ψ,where Cεx=ε+δ2/2cx+η2/2dx, A(x) = a(x) − δc(x) + ηd(x), and B(x) = b(x) − c(x) − d(x). Moreover, the terms a(x), c(x), d(x), δ, and η are such that |A(x) |≥ 2α > 0, for τ < |x| ≤ 1, for some τ > 0, and later on, we will use the term Cɛ to denote the constant part of Cɛ(x) (since c(x), d(x) are bounded and δ, η are small parameters, we have Cɛ = O(ɛ)).The solution of problem (9) is an approximation to the solution of the original problem (1) and (2).We establish some a priori results about the solutions and their derivatives for the modified problem (9). Hereinafter, we divide the interval Ω¯ in to three subintervals as Ω1 = [−1, −τ], Ω2 = [−τ, τ], and Ω3 = [τ, 1] such that Ω¯=Ω1∪Ω2∪Ω3, where 0 < τ ≤ 1/2.First, we consider the following property of the operatorL of (9).Lemma 1 (Maximum principle).
Letπ(x) be any sufficiently smooth function satisfying π(−1) ≥ 0 and π(1) ≥ 0, such that Lπ(x) ≥ 0 for all x ∈Ω. Then, π(x) ≥ 0 for all x∈Ω¯.Proof.
Letx∗ be an arbitrary point in Ω = (−1, 1) such that πx∗=minx∈Ω¯πx and assume that π(x∗) < 0. Clearly x∗ ∉ {−1, 1}, and from the definition of x∗, we have π′(x∗) = 0 and π″(x∗) ≥ 0. But then,(10)Lπx∗=−Cεxπ″x∗−Ax∗π′x∗+Bx∗πx∗≤0,which is a contradiction. It follows that our assumption π(x∗) < 0 is wrong, so π(x∗) ≥ 0. Since x∗ is an arbitrary point, π(x) ≥ 0, ∀x∈Ω¯=−1,1.Using the maximum principle, it is easy to prove that:Lemma 2 (Stability Result).
Letu(x) be the solution of the TPP (9). Then ∀Cɛ > 0 we have(11)u≤β−1f+maxϕ,ψ,∀x∈Ω¯=−1,1.Proof.
First we consider the barrier functionsφ± defined by(12)φ±x=β−1f+maxϕ,ψ±ux.
Then it is easy to show thatφ±(−1) ≥ 0 and φ±(1) ≥ 0, and(13)Lφ±x=−Cεxφ±x″−Axφ±x′+Bxφ±x=Bxβ−1f+maxϕ,ψ±Lux=Bxβ−1f+maxϕ,ψ±fx≥f±fx+βmaxϕ,ψ≥0.
Therefore, from Lemma1, we obtain φ±(x) ≥ 0 for all x ∈ [−1, 1], which gives the desired estimate.The following theorem gives estimates foru and its derivatives in the interval Ω1 and Ω3, which excludes the turning point x = 0.Theorem 1.
LetA(x), B(x), and fx∈CmΩ¯, m > 0, |A(x)| ≥ 2α > 0, ∀x ∈ Ω1 ∪ Ω1. Then there exists a positive constant M, such that for A(x) > 0 on Ω1, the solution u(x) of problem (9) satisfies(14)uix≤M1+Cε−iexp−α1−xCε,fori=1,…,m+1,and for A(x) < 0 on Ω3, we have(15)uix≤M1+Cε−iexp−α1−xCε,fori=1,…,m+1.Proof.
For the proof of this theorem, the reader can refer [25].Ifλ = B(0)/a′(0) < 0, then the solution u(x) is smooth near the turning point x = 0 [26]. Using this, the following theorem gives the bound for the derivatives of the solution in the interval Ω2 which contains the turning point x = 0.Theorem 2.
Letλ < 0. If u(x) is the solution of (9) and satisfies all conditions from (3) to (6), let A, B, and f∈CmΩ¯, for m > 0. Then there exists a positive constant M, such that(16)uix≤M,fori=1,…,mand∀x∈Ω2.Proof.
For the proof one can refer [25, 26].Finally, to prove the uniform convergence of the proposed numerical method, we need to consider the following theorem which provides bounds for the smooth and singular components of the exact solutionu of problem (9).Theorem 3.
LetA, B and f∈C4Ω¯, and assume that the solution u(x) of problem (9) is decomposed in to the smooth and singular components as(17)ux≔vx+wx,∀x∈Ω¯.
Then for alli, 0 ≤ i ≤ 4 the smooth component satisfies(18)vix≤M1+Cε−i−3ex,α,∀x∈Ω¯,and the singular component satisfies(19)wix≤MCε−iex,α,∀x∈Ω¯,where e(x, α) = (exp(−α(1 + x)/Cɛ) + exp(−α(1 − x)/Cɛ)).Proof.
For the proof of this theorem, the reader can refer [25, 27].


## 3. Discrete Problem

In this section, we describe the piecewise-uniform Shishkin mesh for the discretization of the domain and study the behavior of the hybrid difference scheme used for the modified problem (9).


### 3.1. Piecewise-Uniform Shishkin Mesh

Consider the domainΩ¯=−1,1 and let N = 8k and k > 0 is a positive integer. Since the TPP (9) has two boundary layers at x = ±1, we construct a piecewise-uniform Shishkin mesh by subdividing the domain Ω¯ into three subintervals ΩL = [−1, −1 + τ], ΩC = [−1 + τ, 1 − τ], and ΩR = [1 − τ, 1] such that Ω¯=ΩL∪ΩC∪ΩR, where the transition parameter τ satisfies 0 < τ ≤ 1/2 and defined by(20)τ=min12,2CεαlnN,whereCε=ε+K1δ22+K2η22.Then the discrete meshΩ¯N is obtained by putting a uniform mesh with N/4 mesh elements in both ΩL and ΩR, and a uniform mesh with N/2 mesh elements in ΩC, such that Ω¯N=xi∈−1,1:i=0,1,…,N.Lethi = xi − xi−1, for i = 1, …, N, denotes the variable step size. Since the mesh is piecewise-uniform, then the mesh elements are given by(21)xi=−1+τ+i−N4h,fori=1,…,N4,i−N2H,fori=N4,…,3N4,1−τ+i−3N4h,fori=3N4,…,N,where h = 4τ/N and H = 4(1 − τ)/N are uniform mesh lengths for the piecewise-uniform meshes.IfCɛ > MN−1, then the mesh becomes equally spaced and the transition parameter becomes τ = 1/2 such that(22)hi=H=h=2N−1,Cε−1<4lnNα.On the other hand, forCɛ ≤ MN−1, the mesh is piecewise-uniform and τ=2Cε/αlnN. Here, we have 2N−1 ≤ H ≤ 4N−1 and(23)hCε=8αN−1lnN,e−α1+xN/4/Cε=e−α1−x3N/4/Cε=N−2.


### 3.2. Hybrid Difference Scheme

Before describing the scheme, for a given mesh functiony(xi) = yi, we define the forward, backward, and central difference operators D+, D−, and D0 by(24)D+yi=yi+1−yihi+1,D−yi=yi−yi−1hi,D0yi=yi+1−yi−1hi¯,respectively, and the second-order central difference operator δ2 by(25)δ2yi=2D+yi−D−yihi¯=2hihi+1hi¯hi+1yi−1−hi¯yi+hiyi+1,where hi¯=hi+hi+1, for i = 1, …, N − 1.Further, we define the midpoint upwind schemesLM±N and the classical central difference scheme LCN used to approximate the continuous operator L as(26)LM+Nyi=−Cε,i+1/2δ2yi−Ai+1/2D+yi+Bi+1/2yi+1/2=fi+1/2,LM−Nyi=−Cε,i−1/2δ2yi−Ai−1/2D−yi+Bi−1/2yi−1/2=fi−1/2,LCNyi=−Cε,iδ2yi−AiD0yi+Biyi=fi,where Ai±1/2 = (Ai + Ai±1)/2 and similarly for Cɛ,i±1/2, Bi±1/2, and fi±1/2.Now, we propose the hybrid difference scheme to discretize (9), which consists of the classical central difference scheme when Cɛ > MN−1, and a proper combination of the midpoint upwind schemes in the outer region ΩC and the central difference scheme in the layer regions ΩL and ΩR, whenever Cɛ ≤ MN−1. Hence, the proposed hybrid scheme on Ω¯N takes the following form:(27)LHNUi=fi¯,fori=1,…,N−1,U0=u0,UN=uN,where(28)LHNUi=LCNUi,i=1,…,N−1,andCε>MN−1,LM+NUi,i=N4,…,N2,andCε≤MN−1,LM−NUi,i=N2+1,…,3N4,andCε≤MN−1,LCNUi,i=1,…,N4−1,i=3N4+1,…,N−1,andCε≤MN−1,and the right hand side vector fi¯ as(29)fi¯=fi,i=1,…,N−1,andCε>MN−1,fi+1/2,i=N4,…,N2,andCε≤MN−1,fi−1/2,i=N2+1,…,3N4,andCε≤MN−1,fi,i=1,…,N4−1,i=3N4+1,…,N−1,andCε≤MN−1.After rearranging the terms in (27), we obtain the following system of equations:(30)LHNUi=ri−Ui−1+ricUi+ri+Ui+1=fi¯,where the coefficients are given by(31)ri−=−2Cε,ihihi¯+Aihi¯,ri+=−2Cε,ihi+1hi¯−Aihi¯,ric=2Cε,ihihi+1+Bi,ifLHN≡LCN,(32)ri−=−2Cε,i+1/2hihi¯,ri+=−2Cε,i+1/2hi+1hi¯−Ai+1/2hi+1+Bi+1/22,ric=2Cε,i+1/2hihi+1+Ai+1/2hi+1+Bi+1/22,ifLHN≡LM+N,(33)ri−=−2Cε,i−1/2hihi¯+Ai−1/2hi+Bi−1/22,ri+=−2Cε,i−1/2hi+1hi¯,ric=2Cε,i−1/2hihi+1−Ai−1/2hi+Bi−1/22,ifLHN≡LM−N.In general, central difference schemes can be unstable on coarser meshes, but we use this scheme only on the fine part of the Shishikn mesh and thus attain stability under the mild assumption on the minimum number of mesh pointsN, considered in the following lemma.Lemma 3.
Assume that there exists a positive integerN0 and for all N ≥ N0 such that(34)NlnN≥4Aα,αN2≥B,holds true. Then the discrete operator defined by (27) satisfies a discrete maximum principle, i.e., if Ui and Di are mesh functions that satisfie U0 ≤ D0, UN ≤ DN, and LHNUi≤LHNDi, for i = 1, …, N − 1, then Ui ≤ Di, for i = 0, …, N. Hence, (27) has a unique solution.Proof.
To prove the results, it is enough to show that the operator given by (30) is an M-matrix. For this, we need to show that (30) satisfies(35)ri−<0,ri+<0,ri−+ric+ri+>0,for all the operators defined in (31)–(33). Here, we separately consider the following two cases based on the relation between Cɛ and N.Case 1: whenCɛ > MN−1, the mesh is uniform, and we used central difference scheme on the entire domain. Thus, for M≥A and using (22) into (31), we get(36)ri−=−Cε,ih2+Ai2h=1h2−Cε,i+12Aih≤1h2−Cε+AN−1<0,ri+=−Cε,ih2−Ai2h=1h2−Cε,i−12Aih≤1h2−Cε+AN−1<0,and some simple calculations gives ri−+ric+ri+=Bi>0, for all 1 ≤ i ≤ N − 1.Case 2: whenCɛ ≤ MN−1, different operators are used in the layer regions and the outer regions.
In the layer regions, it is apparent thatri+<0 and ri−<0, for 1 ≤ i < N/4 and 3N/4 < i ≤ N − 1, respectively. Further, using the first assumption of (34) and (23) in (31) we get(37)ri−=−Cε,ih2+Ai2h≤1h−Cεh+Ai2≤1h−αN8lnN+A2<0,ri+=−Cε,ih2−Ai2h≤1h−Cεh−Ai2≤1h−αN8lnN+A2<0,for 1≤i≤N/4−1 and 3N/4+1≤i≤N−1, respectively.
In both the layer regions, we simply obtainri−+ric+ri+=Bi>0.
Finally, in the outer regions, it is straightforward thatri−<0 for N/4 ≤ i ≤ N/2 and ri+<0 for N/2 + 1 ≤ i ≤ 3N/4. In addition, using H ≤ 4N−1 and the second assumption of (34) in (32) and (33) gives us(38)ri+=−2Cε,i+1/2hi+1hi¯−Ai+1/2hi+1+Bi+1/22≤−2Cε,i+1/2hi+1hi¯−αN4+B2<0,ri−=−2Cε,i−1/2hihi¯+Ai−1/2hi+Bi−1/22≤−2Cε,i−1/2hihi¯−αN4+B2<0,for N/4≤i≤N/2 and N/2+1≤i≤3N/4, respectively.
Moreover, for allN/4 ≤ i ≤ 3N/4, it is easy to verify that ri−+ric+ri+>0.
For all the cases, it is verified that the operator (30) satisfies the conditions in (35). Hence, the matrix is an M-matrix. Therefore, the solution of (27) exists and the maximum principle easily follows. For more details the reader can refer [28, 29].Whenever, the conditions of the maximum principle are satisfied, we can take {Di} as a barrier function for {Ui}.


## 3.1. Piecewise-Uniform Shishkin Mesh

Consider the domainΩ¯=−1,1 and let N = 8k and k > 0 is a positive integer. Since the TPP (9) has two boundary layers at x = ±1, we construct a piecewise-uniform Shishkin mesh by subdividing the domain Ω¯ into three subintervals ΩL = [−1, −1 + τ], ΩC = [−1 + τ, 1 − τ], and ΩR = [1 − τ, 1] such that Ω¯=ΩL∪ΩC∪ΩR, where the transition parameter τ satisfies 0 < τ ≤ 1/2 and defined by(20)τ=min12,2CεαlnN,whereCε=ε+K1δ22+K2η22.Then the discrete meshΩ¯N is obtained by putting a uniform mesh with N/4 mesh elements in both ΩL and ΩR, and a uniform mesh with N/2 mesh elements in ΩC, such that Ω¯N=xi∈−1,1:i=0,1,…,N.Lethi = xi − xi−1, for i = 1, …, N, denotes the variable step size. Since the mesh is piecewise-uniform, then the mesh elements are given by(21)xi=−1+τ+i−N4h,fori=1,…,N4,i−N2H,fori=N4,…,3N4,1−τ+i−3N4h,fori=3N4,…,N,where h = 4τ/N and H = 4(1 − τ)/N are uniform mesh lengths for the piecewise-uniform meshes.IfCɛ > MN−1, then the mesh becomes equally spaced and the transition parameter becomes τ = 1/2 such that(22)hi=H=h=2N−1,Cε−1<4lnNα.On the other hand, forCɛ ≤ MN−1, the mesh is piecewise-uniform and τ=2Cε/αlnN. Here, we have 2N−1 ≤ H ≤ 4N−1 and(23)hCε=8αN−1lnN,e−α1+xN/4/Cε=e−α1−x3N/4/Cε=N−2.


## 3.2. Hybrid Difference Scheme

Before describing the scheme, for a given mesh functiony(xi) = yi, we define the forward, backward, and central difference operators D+, D−, and D0 by(24)D+yi=yi+1−yihi+1,D−yi=yi−yi−1hi,D0yi=yi+1−yi−1hi¯,respectively, and the second-order central difference operator δ2 by(25)δ2yi=2D+yi−D−yihi¯=2hihi+1hi¯hi+1yi−1−hi¯yi+hiyi+1,where hi¯=hi+hi+1, for i = 1, …, N − 1.Further, we define the midpoint upwind schemesLM±N and the classical central difference scheme LCN used to approximate the continuous operator L as(26)LM+Nyi=−Cε,i+1/2δ2yi−Ai+1/2D+yi+Bi+1/2yi+1/2=fi+1/2,LM−Nyi=−Cε,i−1/2δ2yi−Ai−1/2D−yi+Bi−1/2yi−1/2=fi−1/2,LCNyi=−Cε,iδ2yi−AiD0yi+Biyi=fi,where Ai±1/2 = (Ai + Ai±1)/2 and similarly for Cɛ,i±1/2, Bi±1/2, and fi±1/2.Now, we propose the hybrid difference scheme to discretize (9), which consists of the classical central difference scheme when Cɛ > MN−1, and a proper combination of the midpoint upwind schemes in the outer region ΩC and the central difference scheme in the layer regions ΩL and ΩR, whenever Cɛ ≤ MN−1. Hence, the proposed hybrid scheme on Ω¯N takes the following form:(27)LHNUi=fi¯,fori=1,…,N−1,U0=u0,UN=uN,where(28)LHNUi=LCNUi,i=1,…,N−1,andCε>MN−1,LM+NUi,i=N4,…,N2,andCε≤MN−1,LM−NUi,i=N2+1,…,3N4,andCε≤MN−1,LCNUi,i=1,…,N4−1,i=3N4+1,…,N−1,andCε≤MN−1,and the right hand side vector fi¯ as(29)fi¯=fi,i=1,…,N−1,andCε>MN−1,fi+1/2,i=N4,…,N2,andCε≤MN−1,fi−1/2,i=N2+1,…,3N4,andCε≤MN−1,fi,i=1,…,N4−1,i=3N4+1,…,N−1,andCε≤MN−1.After rearranging the terms in (27), we obtain the following system of equations:(30)LHNUi=ri−Ui−1+ricUi+ri+Ui+1=fi¯,where the coefficients are given by(31)ri−=−2Cε,ihihi¯+Aihi¯,ri+=−2Cε,ihi+1hi¯−Aihi¯,ric=2Cε,ihihi+1+Bi,ifLHN≡LCN,(32)ri−=−2Cε,i+1/2hihi¯,ri+=−2Cε,i+1/2hi+1hi¯−Ai+1/2hi+1+Bi+1/22,ric=2Cε,i+1/2hihi+1+Ai+1/2hi+1+Bi+1/22,ifLHN≡LM+N,(33)ri−=−2Cε,i−1/2hihi¯+Ai−1/2hi+Bi−1/22,ri+=−2Cε,i−1/2hi+1hi¯,ric=2Cε,i−1/2hihi+1−Ai−1/2hi+Bi−1/22,ifLHN≡LM−N.In general, central difference schemes can be unstable on coarser meshes, but we use this scheme only on the fine part of the Shishikn mesh and thus attain stability under the mild assumption on the minimum number of mesh pointsN, considered in the following lemma.Lemma 3.
Assume that there exists a positive integerN0 and for all N ≥ N0 such that(34)NlnN≥4Aα,αN2≥B,holds true. Then the discrete operator defined by (27) satisfies a discrete maximum principle, i.e., if Ui and Di are mesh functions that satisfie U0 ≤ D0, UN ≤ DN, and LHNUi≤LHNDi, for i = 1, …, N − 1, then Ui ≤ Di, for i = 0, …, N. Hence, (27) has a unique solution.Proof.
To prove the results, it is enough to show that the operator given by (30) is an M-matrix. For this, we need to show that (30) satisfies(35)ri−<0,ri+<0,ri−+ric+ri+>0,for all the operators defined in (31)–(33). Here, we separately consider the following two cases based on the relation between Cɛ and N.Case 1: whenCɛ > MN−1, the mesh is uniform, and we used central difference scheme on the entire domain. Thus, for M≥A and using (22) into (31), we get(36)ri−=−Cε,ih2+Ai2h=1h2−Cε,i+12Aih≤1h2−Cε+AN−1<0,ri+=−Cε,ih2−Ai2h=1h2−Cε,i−12Aih≤1h2−Cε+AN−1<0,and some simple calculations gives ri−+ric+ri+=Bi>0, for all 1 ≤ i ≤ N − 1.Case 2: whenCɛ ≤ MN−1, different operators are used in the layer regions and the outer regions.
In the layer regions, it is apparent thatri+<0 and ri−<0, for 1 ≤ i < N/4 and 3N/4 < i ≤ N − 1, respectively. Further, using the first assumption of (34) and (23) in (31) we get(37)ri−=−Cε,ih2+Ai2h≤1h−Cεh+Ai2≤1h−αN8lnN+A2<0,ri+=−Cε,ih2−Ai2h≤1h−Cεh−Ai2≤1h−αN8lnN+A2<0,for 1≤i≤N/4−1 and 3N/4+1≤i≤N−1, respectively.
In both the layer regions, we simply obtainri−+ric+ri+=Bi>0.
Finally, in the outer regions, it is straightforward thatri−<0 for N/4 ≤ i ≤ N/2 and ri+<0 for N/2 + 1 ≤ i ≤ 3N/4. In addition, using H ≤ 4N−1 and the second assumption of (34) in (32) and (33) gives us(38)ri+=−2Cε,i+1/2hi+1hi¯−Ai+1/2hi+1+Bi+1/22≤−2Cε,i+1/2hi+1hi¯−αN4+B2<0,ri−=−2Cε,i−1/2hihi¯+Ai−1/2hi+Bi−1/22≤−2Cε,i−1/2hihi¯−αN4+B2<0,for N/4≤i≤N/2 and N/2+1≤i≤3N/4, respectively.
Moreover, for allN/4 ≤ i ≤ 3N/4, it is easy to verify that ri−+ric+ri+>0.
For all the cases, it is verified that the operator (30) satisfies the conditions in (35). Hence, the matrix is an M-matrix. Therefore, the solution of (27) exists and the maximum principle easily follows. For more details the reader can refer [28, 29].Whenever, the conditions of the maximum principle are satisfied, we can take {Di} as a barrier function for {Ui}.


## 4. Convergence Analysis of the Proposed Method

In this section, we establish theɛ-uniform error estimate of the hybrid scheme (27). For this, we consider the two cases Cɛ > MN−1 and Cɛ ≤ MN−1 separately.For both the cases, analogous to the continuous solutionu, we decompose the discrete solution U into a smooth component V and a singular component W, such that U≔V + W, where V is the solution of the nonhomogeneous problem given by(39)LHNVi=f¯,fori=1,…,N−1,V0=v−1,VN=v1,and W the solution of the homogeneous problem(40)LHNWi=0,fori=1,…,N−1,W0=w−1,WN=w1.Then the error at each mesh point is(41)Ui−uxi=Vi−vxi+Wi−wxi,which implies(42)Ui−uxi≤Vi−vxi+Wi−wxi,and so the error in the smooth and singular components of the solution can be estimated separately.First, to bound the errors we need to consider the truncation error of associated with the discrete operators in (27). For any smooth function y(x), the truncation errors LM±N applied to y at xi±1/2 and LCN applied to y at yi, becomes T1±:=LM±Nyi−Lyxi±1/2 and T2:=LCNyi−Lyxi respectively, where yi : = y(xi). Thus, the bounds are given in the following Lemma.Lemma 4.
Lety(x) be a smooth function defined on [−1, 1]. Then there exists a positive constant M such that(43)T1+≤MCε∫xi−1xi+1y‴tdt+Mhi+1∫xixi+1y‴tdt,T1−≤MCε∫xi−1xi+1y‴tdt+Mhi∫xi−1xiy‴tdt,T2≤Mh∫xi−1xi+1Cεy4t+y‴tdt,forhi=hi+1=h.Proof.
By repeated use of the fundamental theorem of calculus, one can obtain the proof as in Lemma 3.3 of [28].
To bound the truncation error of the scheme the comparison principle of Lemma3 alone is not enough, so we will consider the following lemma which enables us to bound the error.Lemma 5.
LetZi = 2 + xi for 0 ≤ i ≤ N be the mesh function for (27). Then there exists a positive constant M such that(44)LHNZi≥M,for1≤i≤N−1.Proof.
The proof is an easy computation.
Sometimes, the truncation error contains a term of magnitude greater than the desired order of convergence, when this happens we shall combine Lemma3 with the following results.
WheneverCɛ ≤ MN−1, we define the auxiliary discrete function on the mesh Ω¯N by using the mesh elements given in (21) as(45)Si≔21+γhε−N/4,fori=0,…,34N,1+γhε−N/4+1+γhε−N−i,fori=34N,…,N,where γ is a positive constant.Lemma 6.
For anyγ > 0 the discrete function {Si} from (45) is such that(46)e−γ1−xi/ε≤Si≤CN−2γ/α,fori≤34N,C,fori>34N.(47)LHNSi=BiSi,i<N4,Bi+1/2Si,N4≤i≤N2,Bi−1/2Si,N2<i<3N4,Bi−1/2−γh+HSi,i=3N4,1+ηhCε−N−i−1γCε−Ai−γ−Aiγ2h2Cε2+BiSi,i>3N4.Proof.
The lower bound forSi follows from the inequality e−t ≤ (1 + t)−1 which holds true for t ≥ 0. The upper bound for Si is obvious for i > 3N/4. For i ≤ 3N/4, it follows from the inequality 1+t−1≤e−t+t2, which holds true for t ≥ 0. Setting t ≔ γh/ɛ and using (23), we get(48)Si=21+γhε−N/4≤2e−γh/ε+γh/ε2N/4=2eγh/εγh/ε−1N/4=2eγ8N−1lnN/αγ8N−1lnN/α−1N/4=2e2γ/αlnN8γ/αN−1lnN−1=2elnN16γ2/α2N−1lnN−2γ/α=2N−2γ/α.N16γ2/α2N−1lnN≤CN−2γ/α,because the sequence N16γ2/α2N−1lnN is bounded for N ≥ 8. This proves the upper bound for Si. The property (47) can be obtained by direct calculation.Lemma 7.
LetSi¯:=2α/βHS3N/4+Si be a discrete function, where {Si} is from (45) with γ = α, then(49)0<Si¯≤CN−1,fori≤34N,C,fori>34N.LHNSi¯≥CNe−α1−xi/ε,fori≤34N,Cε−1e−α1−xi/ε,fori>34N.Proof.
The proof is similarly to Lemma 3.3 of [30].Lemma 8.
LetSiˇ:=2α/βh+HS3N/4+Si be a discrete function, where {Si} is from (45) with γ = 2α, then(50)0<Siˇ≤CN−3,fori≤34N,C,fori>34N,Siˇ≥e−α1−xi/εfori=0,N,andLHNSiˇ≥0.Proof.
For the proof, one can follow similarly like Lemma 3.4 of [30].Remark 1.
Because of the symmetry of the mesh and the adaptive nature of the hybrid scheme, it is easy to derive a similar result like Lemmas7 and 8 using the mesh function {SN−i} related to the layer function e−α1+x/Cε.
Now we have assembled the tools for the proof of theɛ-uniform convergence.Theorem 4.
Assume that the conditions of (34) holds true. Then the hybrid scheme (27) satisfies the following error estimates.Case 1: forCɛ > MN−1, we have(51)Ui−uxi≤MN−2ln3N,fori=0,…,N.Case 2: forCɛ ≤ MN−1, we have(52)Ui−uxi≤MN−2,forN4≤i≤3N4,MN−2ln2N,for0≤i<N4and3N4<i≤N,where Ui the solution of the discrete problem (27) and u(xi) is the solution of the continuous problem (9) at the mesh points in Ω¯N.Proof.
Here we estimate the error bounds by considering the two cases separately as follows:Case 1: whenCɛ > MN−1, we employed the central difference scheme on the entire domain. First let us compute the nodal error for the smooth part Vi. Using Lemma 4 and the bound of v in Theorem 3, the truncation error of the scheme becomes(53)LHNVi−vxi≤Mh∫xi−1xi+1Cεv4t+v′′′tdt≤Mh2Cε+1+Mh∫xi−1xi+1e−α1+t/Cε+e−α1−t/Cεdt=Mh2Cε+1+MhCεe−α1+xi/Cε+e−α1−xi/CεsinhαhCε≤Mh2Cε+1+Mh2e−α1+xi/Cε+e−α1−xi/Cε,sinht≤Mtfor0≤t≤1≤Mh2Cε+1+Mh2≤Mh2.
Here, usingh = 2N−1 on the above inequality, we obtain the following estimate:(54)LHNVi−vxi≤MN−2,fori=1,…,N−1.
Now, let us takeDi ≔ MN−2(2 + xi) as a barrier function for Vi−vxi, then from (39), it is easy to see that V0−vx0=0≤D0, VN−vxN=0≤DN, and from (54) together with Lemma 5, we observe that LHNVi−vxi≤LHNDi, for i = 1, …, N − 1. Thus, invoking Lemma 3 we get(55)Vi−vxi≤MN−2,fori=0,…,N.
Next, we analyze the error bounds for the singular componentWi. The local truncation error is bounded in standard way as done above. More precisely,(56)LHNWi−wxi=LCNWi−wxi=L−LCNwxi≤Mh∫xi−1xi+1Cεw4t+w‴tdt.
Application of Theorem3 and using (22) on the above inequality and simplifying gives(57)LHNWi−wxi≤Mh2Cε−3≤MN−2ln3N,fori=1,…,N−1.
Now, arguing similarly like the smooth part, we obtain(58)Wi−wxi≤MN−2ln3N,fori=0,…,N.
Using (55) and (58) in to (42) gives the required result of the first case (51).Case 2: forCɛ ≤ MN−1, the mesh becomes piecewise-uniform and we employed combinations of midpoint upwind and central difference schemes. First, let us compute the nodal error for the smooth part Vi. Similarly like the smooth part of Case 1, the truncation error becomes(59)LHNVi−vxi=Lvxi+1/2−LM+Nvxi,fori=N4,…,N2,Lvxi−1/2−LM−Nvxi,fori=N2+1,…,3N4,L−LCNvxi,otherwise,then using the bound for the truncation error of Lemma 4 and the bound for vx from Theorem 3, we get(60)LHNVi−vxi≤MCε+hi+1hi+hi+1,fori=N4,…,N2,MCε+hihi+hi+1,fori=N2+1,…,3N4,Mh2,otherwise.
Since,Cɛ ≤ MN−1 and hi ≤ 4N−1, then using these in the above inequality gives us(61)LHNVi−vxi≤MN−2,fori=1,…,N−1.
Now, arguing similarly like the smooth part of the previous case, we obtain(62)Vi−vxi≤MN−2,fori=0,…,N.
Next, we analyze the error bounds for the singular componentWi. A different argument is used to bound W−w in the outer region and layer regions.
In the subinterval with no boundary layerΩ¯C, both W and w are small, and by the triangle inequality, we have(63)W−wxi≤Wxi+wxi,so it suffices to bound W(xi) and wxi separately. Theorem 3 for i = N/4, …, 3N/4 gives(64)wxi≤Me−α1+xi/Cε+e−α1−xi/Cε≤Me−α1+xN/4/Cε+e−α1−x3N/4/Cε.
Then using (23) in the above inequality, we get(65)wxi≤MN−2,fori=N4,…,3N4.
To obtain a similar bound forW(xi), we set Di≔M1N−1Si¯+N−1SN¯−i+Siˇ+SN˘−i for i = 0, …, N, where Si¯ and Siˇ are from Lemmas 7 and 8, respectively. Now for sufficiently large M1, using Theorem 3 in (40) and Lemma 8, we get(66)W0=w−1≤Me−α1+x0/Cε+e−α1−x0/Cε≤M1S0ˇ+SNˇ≤D0,(67)WN=w−1≤Me−α1+xN/Cε+e−α1−xN/Cε≤M1SNˇ+S0ˇ≤DN.
Further, fori = 1, …, N − 1, the property of the discrete operator from Lemmas 7 and 8 implies(68)LHNDi≥M1MCε−1N−1e−α1+xi/Cε+e−α1+xi/Cε,fori=1,…,N4−1,M1e−α1+xi/Cε+e−α1−xi/Cε,fori=N4,…,3N4,M1e−α1+xi/Cε+MCε−1N−1e−α1+xi/Cε,fori=3N4+1,…,N−1,since Cɛ ≤ MN−1 implies MCε−1N−1≥1, using this in the above inequality, we get(69)LHNDi≥M1e−α1+xi/Cε+e−α1+xi/Cε≥0=LHNWi,∀1≤i≤N−1.
From (66)–(69), we can easily observe that Di can be a barrier function for Wi for M1 sufficiently large. Therefore, by the discrete maximum principle of Lemma 3 we get(70)Wi≤Di,fori=0,…,N.
In particular, in the coarser region, (70) and Lemmas 7 and 8 together imply that(71)Wi≤Di≤MN−2,fori=N4,…,3N4.
Therefore, combining (63), (65), and (71) we get(72)Wi−wxi≤MN−2,fori=N4,…,3N4.
It remains to prove the bound for the singular component in the layer regionsΩ¯L and Ω¯R. First we estimate the bounds for the singular component in Ω¯R. Since we employ central difference scheme in Ω¯R, so as we did for the smooth component, we use the truncation error to bound the error. Thus, from Lemma 4 and Theorem 3 we get(73)LHNWi−wxi≤Mh∫xi−1xi+1Cεw4t+w‴tdt≤MhCε−3∫xi−1xi+1e−α1−t/Cεdt=MhCε−2e−α1−xi+1/Cε−e−α1−xi−1/Cε=MhCε−2e−α1−xi/CεsinhαhCε.
Clearly, the first assumption of (34) implies αh/Cɛ ≤ 1 and since sinh t ≤ Mt for 0 ≤ t ≤ 1, so sinh(αh/Cɛ) ≤ Mαh/Cɛ. Thus, for i = 3N/4 + 1, …, N − 1 the above inequality is reduced to(74)LHNWi−wxi≤MhCε2Cε−1e−α1−xi/Cε.
Further, takingi = 3N/4 in (72), we get W3N/4−w3N/4≤MN−2, and for i = N, the boundary condition in (40) gives WN−wxN=0. Now let(75)Di≔M2N−2+hCε2Si¯,fori=3N4,…,N,be the mesh function, where Si¯ is from Lemma 7. If M2 is chosen large enough, our estimates shows that Di is a barrier function for Wi−wxi. So by using the discrete maximum principle of Lemma 3 and Lemma 7, together with (23), we get(76)Wi−wxi≤MN−2ln2N,fori=3N4,…,N.
Similarly the proof follows for the left boundary layer region,Ω¯L, i.e.,(77)Wi−wxi≤MN−2ln2N,fori=0,…,N4.
Finally, properly using (62), (72), (76), and (77) in to (42) gives the required bound of the second case (52), which completes the proof.


## 5. Test Problems and Numerical Results

To demonstrate the applicability of the proposed method, we have implemented it on two boundary value problems of the form (1) and (2). Since the exact solution for the problems are not available, the pointwise error e^iN and maximum absolute error (ÊN) are calculated by using the double mesh principle given by(78)e^iN=UNxi−U2Nxi,E^N=maxie^iN,where UN and U2N denotes the numerical solutions obtained using N and 2N meshes points, respectively. Further, we determine the corresponding rate of convergence by(79)R^N=log2E^N/E^2N.Example 1.
Consider the following homogeneous SPDDE with a turning point atx = 0:(80)−εy″x+xy′x+3yx−yx−δ−0.5yx+η=0,x∈−1,1yx=1,−1−δ≤x≤−1,yx=1,1≤x≤1+η.
In Tables1 and 2, the maximum pointwise errors and the rate of convergences for Example 1 are displayed for different values of N and ɛ, respectively. The plots of the approximate solutions for N = 64, ɛ = 0.1, with different values of δ and η are shown in Figures 1 and 2.Table 1
Result of Example1 (maximum pointwise error) for δ = η = 0.5 ∗ ɛ.
ɛ
↓/N ⟶
32
64
128
256
512
1024
1.00
1.90E − 04
4.74E − 05
1.19E − 05
2.96E − 06
7.41E − 07
1.85E − 07
0.50
3.67E − 04
9.16E − 05
2.29E − 05
5.73E − 06
1.43E − 06
3.58E − 07
10–1
5.40E − 03
1.35E − 03
3.36E − 04
8.41E − 05
2.10E − 05
5.25E − 06
10–2
1.83E − 02
6.16E − 03
2.01E − 03
6.72E − 04
2.11E − 04
6.51E − 05
10–3
1.90E − 02
6.39E − 03
2.12E − 03
6.89E − 04
2.17E − 04
6.62E − 05
10–4
1.91E − 02
6.41E − 03
2.13E − 03
6.92E − 04
2.18E − 04
6.74E − 05
10–5
1.91E − 02
6.42E − 03
2.13E − 03
6.92E − 04
2.19E − 04
6.74E − 05
10–6
1.91E − 02
6.42E − 03
2.13E − 03
6.92E − 04
2.19E − 04
6.74E − 05
10–7
1.91E − 02
6.42E − 03
2.13E − 03
6.92E − 04
2.19E − 04
6.74E − 05
E
N
1.91E − 02
6.42E − 03
2.13E − 03
6.92E − 04
2.19E − 04
6.74E − 05Table 2
Result of Example1 (rate of convergence) for δ = η = 0.5 ∗ ɛ.
ɛ
↓/N ⟶
32
64
128
256
512
1024
1.00
1.9997
1.9999
1.9999
2.0000
2.0000
2.0001
0.50
2.0006
2.0000
2.0000
2.0000
2.0000
2.0000
10–1
1.9974
2.0094
1.9977
2.0007
2.0002
2.0000
10–2
1.5723
1.6139
1.5822
1.6686
1.6986
1.7256
10–3
1.5711
1.5918
1.6220
1.6627
1.7161
1.6956
10–4
1.5707
1.5920
1.6213
1.6627
1.6971
1.7252
10–5
1.5707
1.5920
1.6212
1.6626
1.6970
1.7249
10–6
1.5707
1.5920
1.6212
1.6626
1.6969
1.7249
10–7
1.5707
1.5920
1.6212
1.6626
1.6969
1.7249
R
N
1.5707
1.5920
1.6212
1.6626
1.6969
1.7249Figure 1
Plot of the solutions of Example1 for ɛ = 0.1 and N = 64.Figure 2
Plot of the solutions of Example1 for ɛ = 0.1 and N = 64.Example 2.
Consider the following nonhomogeneous SPDDE with a turning point atx = 1/2:(81)−εy″x−21−2xy′x+6yx−yx−δ−yx+η=41−4x,x∈0,1yx=1,−δ≤x≤0,yx=1,1≤x≤1+η.
In Tables3 and 4 the maximum pointwise errors and the rate of convergences for Example 2 are displayed for different values of N and ɛ, respectively. The plots of the approximate solutions for N = 256, ɛ = 0.1, with different values of δ and η are shown in Figures 3 and 4.Table 3
Result of Example2 (maximum pointwise error) for δ = η = 0.5 ∗ ɛ.
ɛ
↓/N ⟶
32
64
128
256
512
1024
1.00
2.78E − 04
6.95E − 05
1.74E − 05
4.34E − 06
1.09E − 06
2.71E − 07
0.50
6.92E − 04
1.73E − 04
4.33E − 05
1.08E − 05
2.70E − 06
6.76E − 07
10–1
1.12E − 02
3.86E − 03
9.60E − 04
2.39E − 04
5.98E − 05
1.50E − 05
10–2
5.44E − 02
1.83E − 02
6.11E − 03
1.97E − 03
6.31E − 04
1.94E − 04
10–3
5.68E − 02
1.91E − 02
6.35E − 03
2.06E − 03
6.52E − 04
2.01E − 04
10–4
5.71E − 02
1.92E − 02
6.38E − 03
2.07E − 03
6.55E − 04
2.02E − 04
10–5
5.71E − 02
1.92E − 02
6.38E − 03
2.08E − 03
6.56E − 04
2.02E − 04
10–6
5.71E − 02
1.92E − 02
6.39E − 03
2.08E − 03
6.56E − 04
2.02E − 04
10–7
5.71E − 02
1.92E − 02
6.39E − 03
2.08E − 03
6.56E − 04
2.02E − 04
E
N
5.71E − 02
1.92E − 02
6.39E − 03
2.08E − 03
6.56E − 04
2.02E − 04Table 4
Result of Example2 (rate of convergence) for δ = η = 0.5 ∗ ɛ.
ɛ
↓/N ⟶
32
64
128
256
512
1024
1.00
1.9986
2.0000
2.0000
2.0000
2.0000
2.0000
0.50
1.9999
2.0002
2.0001
2.0000
2.0000
2.0000
10–1
1.5395
2.0098
2.0024
2.0006
1.9999
2.0000
10–2
1.5677
1.5876
1.6326
1.6424
1.6989
1.7256
10–3
1.5700
1.5915
1.6217
1.6628
1.6969
1.7310
10–4
1.5699
1.5916
1.6212
1.6627
1.6971
1.7251
10–5
1.5699
1.5916
1.6211
1.6626
1.6969
1.7249
10–6
1.5699
1.5916
1.6211
1.6626
1.6969
1.7249
10–7
1.5699
1.5916
1.6211
1.6626
1.6969
1.7249
R
N
1.5699
1.5916
1.6211
1.6626
1.6969
1.7249Figure 3
Plot of the solutions of Example2 for ɛ = 0.1 and N = 256.Figure 4
Plot of the solutions of Example2 for ɛ = 0.1 and N = 256.


## 6. Discussion

Singularly perturbed differential-difference turning point problems exhibiting twin boundary layers which contain small mixed shifts on the reaction coefficients are considered. For the numerical treatment of these problems, first we employ a second-order Taylor’s series approximation on the terms containing shift parameters and obtain a modified singularly perturbed problem which approximates the original problem. And then an efficient hybrid difference scheme on an appropriate piecewise-uniform Shishkin mesh is developed for the modified problem. The proposed method is analyzed for stability and convergence, and it has been shown that the method isɛ-uniformly convergent with an almost second-order rate of convergence. Further, two numerical experiments are examined to support the theoretical results and to illustrate the effect of the small shifts on the layer behavior of the solutions.Tables1–4 present the computed maximum pointwise error and the rate of convergence for the considered examples. The results demonstrate that the method is robust, i.e., converges for all ɛ. We also observed that the maximum pointwise error and the rate of convergence stabilizes as ɛ ⟶ 0 for each appropriate N. Further, the numerical results clearly support the theoretical error bounds and order of convergence.In addition, to demonstrate the effect of the small shifts on the behavior of the solution graphs of the considered problems are plotted in Figures1–4 for different values of δ and η. It is observed that the boundary layers are maintained but layer gets shifted as delay/advance argument changes. Shifts in the layer depend upon the type of shift as well as on the value of the coefficients of the term containing delay/advance.



---

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