# Solution of Singularly Perturbed Differential-Difference Equations with Mixed Shifts Using Galerkin Method with Exponential Fitting

Authors: D. Kumara Swamy; K. Phaneendra; Y. N. Reddy
Journal: Chinese Journal of Mathematics
 (2016)

---

## Abstract

Galerkin method is presented to solve singularly perturbed differential-difference equations with delay and advanced shifts using fitting factor. In the numerical treatment of such type of problems, Taylor’s approximation is used to tackle the terms containing small shifts. A fitting factor in the Galerkin scheme is introduced which takes care of the rapid changes that occur in the boundary layer. This fitting factor is obtained from the asymptotic solution of singular perturbations. Thomas algorithm is used to solve the tridiagonal system of the fitted Galerkin method. The method is analysed for convergence. Several numerical examples are solved and compared to demonstrate the applicability of the method. Graphs are plotted for the solutions of these problems to illustrate the effect of small shifts on the boundary layer solution.


---

## Body

## 1. Introduction

Singularly perturbed differential-difference equations (SPDDEs) arise very frequently in the mathematical modelling of real life situations in science and engineering [1–3]. In the mathematical modelling of a physical system as in control theory, the presence of small time parasitic parameters like moments of inertia, resistances, inductances, and capacitances increases the order and stiffness of these systems. The suppression of these small constants results in the reduction of the order of the system. Such systems are termed as singular perturbation systems and when these systems take into account the past history as well as the present state of the physical system then they are called singularly perturbed delay differential equations. Delay differential equations arise in first-exit time problems in neurobiology and in mathematical formulation of various practical phenomena in biosciences. A differential-difference equation with the presence of shift terms induces large amplitudes and exhibits oscillations, resonance, turning point behaviour, and boundary and interior layers. Hence, to control such behaviour, we need some simple and efficient numerical techniques.Lange and Miura [3–7] published a series of papers extending the method of matched asymptotic expansions initially developed for ordinary differential equations to obtain approximate solution of singularly perturbed differential-difference equations.Numerical analysis of singularly perturbed differential-difference turning point problems was initiated by Kadalbajoo and Sharma. In a series of papers, [8–10], they gave many robust numerical techniques for the solution of such type of problems. Kadalbajoo and Sharma [8] elucidate a numerical method to solve boundary value problems for singularly perturbed differential-difference equation with mixed shifts. Kadalbajoo and Sharma [9] proposed a numerical method to solve boundary value problems for a singularly perturbed differential-difference equation of a mixed type, that is, which contains both types of terms having negative shifts as well as positive shifts, and considered the case in which the solution of the problem exhibits rapid oscillations. Kadalbajoo and Sharma [10] described a numerical approach based on finite difference method to solve a mathematical model arising from a model of neuronal variability. Kadalbajoo and Kumar [11] used B-spline collocation method with fitted mesh for the solution of singularly perturbed differential-difference equations with small delay.Patidar and Sharma [12] combined fitted-operator methods with Micken’s nonstandard finite difference techniques for the numerical approximations of singularly perturbed linear delay differential equations. Kadalbajoo et al. [13] derived ε-uniformly convergent fitted methods for the solution of singularly perturbed differential-difference equation (SPDDE). Kumar and Sharma [14] presented a numerical scheme based on B-spline collocation to approximate the solution of boundary value problems for singularly perturbed differential-difference equations with delay and advance.With this motivation, an exponentially fitting factor is introduced in Galerkin method for the solution of singularly perturbed differential-difference equation with delay and advanced parameters. In Section2, description of the problem is given. In Section 3, numerical scheme for the solution of the problem is presented and Section 4 deals with convergence analysis of the proposed scheme. To demonstrate the efficiency of the proposed method, numerical experiments are carried out for several test problems and the results are given in Section 5. Finally the conclusions are given in the last section.


## 2. Description of the Problem

Consider a linear singularly perturbed differential-difference equation of the following form:(1)εy′′x+pxy′x+qxyx-δ+rxyx+sxyx+η=fxon (0, 1), under the boundary conditions(2)yx=φxon-δ≤x≤0,y1=γxon1≤x≤1+η.Here ε is a small parameter such that 0<ε≪1,px,qx,rx,sx,fx,ϕx, and γx are smooth functions, and 0<δ=oε,0<η=oε are, respectively, the delay (negative shift) and the advance (positive shift) parameters. If px-δqx+ηsx>0, the solution of (1) with (2) exhibits layer at the left end of the interval and if px-δqx+ηsx<0, the layer exists at the right end of the interval.Since the solutionyx of boundary value problem equations (1) and (2) is sufficiently differentiable, we expand the terms yx-δ and yx+η using Taylor series; we get(3)yx-δ≈yx-δy′x+Oδ2,(4)yx+η≈yx+ηy′x+Oη2.Using (3) and (4) in (1), we get(5)εy′′x+axy′x+bxyx=fx+Oδ2+η2.Equation (5) is an asymptotically equivalent second-order singular perturbation problem of (1) with boundary conditions(6a)y0=φ0,(6b)y1=γ1.Since 0<δ≪1 and 0<η≪1, the transition from (1) to (5) is admitted. This replacement is significant from the computational point of view. For more details on the validity of this transition, one can refer El’sgol’ts and Norkin [15]. Thus, the solution of (5) provides a good approximation to the solution of (1).Here,(7)ax=px-δqx+ηsx,bx=qx+rx+sx.


## 3. Numerical Scheme

### 3.1. Left-End Boundary Layer Problem

Let0=x0,x1,…,xN=1 be a decomposition of the considered interval 0,1 into N equal intervals with constant mesh length h. Then we have the nodes xi=ih, for 0,1,…,N. Assume that ax,bx, and fx are sufficiently continuously differentiable functions in 0,1. If ax≥M>0 in 0,1 where M is a positive constant, (5) has a unique solution yx which, in general, displays a boundary layer of width Oε at x=0.Lemma 1 (Doolan et al. [16] and O’Malley [17]).
Lety~x=y0+z0 be the zeroth-order asymptotic approximation to the solution of (5), where y0 represents the zeroth-order approximate outer solution (i.e., the solution of the reduced problem of (5)) and z0 represents the zeroth-order approximate solution in the boundary layer region of (5).
Then for a fixed positive integeri,(8)limh→0⁡yih≈y00+φ0-y00e-a0iρwhereρ=hε.Proof.
Lety0x be the solution of the reduced problem of (5)(9)pxy0′x+qxy0x=fx,y01=γ1and z0t is the solution of the boundary value problem (cf. O’Malley [17])(10)z0′′t+p0z0′t=0,z00=φ0-y00,z0∞=0wheret=xε.
From the theory of singular perturbation, the zeroth-order asymptotic approximation to the solution of (4) is (cf. O’Malley [17])(11)y~x≈y0x+a0axφ0-y00e-∫0xax/εdx.
As we are considering the differential equations on sufficiently small subintervals, the coefficients could be assumed to be locally constant. Hence,(12)y~x≈y0x+φ0-y00e-a0/εx.So, at the nodal points, we have(13)y~xi≈y0xi+φ0-y00ea0/εxi,i=0,1,…,N;that is,(14)y~ih≈y0ih+φ0-y00e-a0/εih.Therefore(15)limh→0⁡y~ih≈y00+φ0-y00exp-a0iρfori=0,1,2,3,…,N,where ρ=h/ε.
Now, we consider the difference scheme [18] by Galerkin method as follows:
Select a set of basis functionsφjx,j=0,1,…,N, which will define an interpolation scheme for the approximate solution over a grid of points a=x0<x1<⋯<xN+1=b. For simplicity, we use piecewise Lagrange polynomials lix of first degree as the basis functions. These interpolating polynomials are(16)l0ξ=ξ-ξ1ξ0-ξ1=1-ξ2,(17)l1ξ=ξ-ξ0ξ1-ξ0=1+ξ2in local element coordinates -1≤ξ≤1.TheN nodal values of the approximate solution y at the interior nodes x1,x2,…,xN are determined using this basis. The given boundary conditions determine the value of yx at the end nodes x0 and xN+1. The Galerkin method is now employed to obtain the integral equations; we have(18)εy′′x+axy′x+bxyx,φj=fx,φjforj=1,2,…,Nwhich is an integral equation ∫x0xNεy′′x+axy′x+bxyxϕjdx=∫x0xNfxφjdx.Sincey is sum of piecewise linear Lagrange polynomials, the second-order derivatives appearing in (17) vanish except at the element boundaries xi, where they become infinite.By integration by parts, (18) becomes(19)-εdydx,dϕjdx+axdydx+bxy,ϕj+εdydxϕjab=fx,ϕj.Using the substitution of trial function y(x)=ϕ0x+∑i=1Nyiϕjx into the integral equation (19), we have(20)∑i=1Nyiεdϕidx,dϕjdx-∑i=1Nyiaxdϕidx+bxϕi,ϕj=-αdl0dx,dϕjdx-βdlN+1dx,dϕjdx+αaxdl0dx+bxl0,ϕj+βaxdlN+1dx+bxlN+1,ϕj+εdydxϕjab-fx,ϕjfor j=1,2,…,N.It can be observed that all quantities on the right side of (20) can be computed from known boundary data to obtain N equations in the N unknown values yi at the interior nodes.The integrals in (20) can be solved by taking advantage of local coordinate ξ system.Since(21)x=b-a2ξ+b+a2=h2ξ+b+a2,ξ=2hx-b+a2⟹dξdx=2h,we have, by simple integration,(22)∫abdϕidxdϕjdxdx=∫-11dl0dξdξdxdl1dξdξdxh2dξ=-1h,fori=j-1,∫abdϕidxdϕjdxdx=∫-11dl1dξdξdxdl0dξdξdxh2dξ=-1h,fori=j+1,∫abdϕidxdϕjdxdx=∫-11dl1dξdξdxdl1dξdξdxh2dξ=1h,fori=j,∫abdϕidxdϕjdxdx=0fori-j>1.By assuming ax,bx, and fx as constants, the integral equation (20) gives, for a typical internal node j,(23)yj-1-εh+a2-bh6+yj2εh-2bh3+yj+1-εh-a2-bh6=-fjh.Equation (23), when rearranged, gives the following system of difference equations:(24)εyj+1-2yj+yj-1h2+ayj+1-yj-12h+byj-1+4yj+yj+16=fj;for1≤j≤N-1.Now, introduce a fitting factor σ in the Galerkin scheme as follows:(25)σεyj+1-2yj+yj-1h2+ajyj+1-yj-12h+bjyj-1+4yj+yj+16=fjfor1≤j≤N-1 with y0=ϕ0,yN=γ1. Here σ is a fitting factor which is to be determined in such a way that the solution of (25) converges uniformly to the solution of (5). Multiplying (25) by h and taking the limit as h→0 (in [16]), we get (26)limh→0⁡σρyj+1-2yj+yj-1+12ajyj+1-yj-1=0iffj-byjis  bounded(27)∴limh→0⁡σρyjh+h-2yjh+yjh-h+12ajyjh+h-yjh-h=0.Now, approximating the solutiony(x) by zeroth-order asymptotic approximation y~(x) and using Lemma 1, we have(28)limh→0⁡y~ih=y00+φ0-y00e-a0iρ,limh→0⁡y~jh+h=y00+φ0-y00e-a0jρe-a0ρ,limh→0⁡y~jh-h=y00+φ0-y00e-a0jρea0ρ.Using the above equations in (27), we get(29)σ=ρ2a0coth⁡a20-εb0a0ρ2.From (25), we have(30)yj-1εσh2-aj2h+bj6-yj2εσh2-2bj3+yj+1εσh2+aj2h+bj6=fjfor j=1,2,…,N-1.Equation (30) can be written as a three-term recurrence relation as follows:(31)Ejyj-1-Fjyi+Gjyj+1=Hj;j=1,2,…,N-1,where(32)Ej=εσh2-aj2h+bj6,Fj=2εσh2-2bj3,Gj=εσh2+aj2h+bj6,Hj=fj.The tridiagonal system equation (31) is solved using Thomas algorithm.


### 3.2. Right-End Layer Problems

We now discuss the method for singularly perturbed two-point boundary value problems with right-end boundary layer of the underlying interval. Assume thatax,bx, and fx are sufficiently continuously differentiable functions in 0,1. Furthermore, assume that ax≤M<0 in 0,1, where M is a negative constant. Under these assumptions, (5) has a unique solution yx which, in general, displays a boundary layer of width Oε at x=1.Lemma 2.
Letyx=y0+z0 be the zeroth-order asymptotic approximation to the solution of (5), where y0 represents the zeroth-order approximate outer solution and z0 represents the zeroth-order approximate solution in the boundary layer region.
Then for a fixed positive integeri,(33)limh→0⁡yih≈y00+β-y01expa11ε-iρwhereρ=hε.Proof.
The proof is based on asymptotic analysis (Doolan et al. [16] and O’Malley [17]) and is similar to the proof of Lemma 1.Applying the same procedure as in Section3 and using Lemma 2, we get the tridiagonal system equation (20) with fitting factor as(34)σ=ρ2a0coth⁡a21-εb1a1ρ2.


## 3.1. Left-End Boundary Layer Problem

Let0=x0,x1,…,xN=1 be a decomposition of the considered interval 0,1 into N equal intervals with constant mesh length h. Then we have the nodes xi=ih, for 0,1,…,N. Assume that ax,bx, and fx are sufficiently continuously differentiable functions in 0,1. If ax≥M>0 in 0,1 where M is a positive constant, (5) has a unique solution yx which, in general, displays a boundary layer of width Oε at x=0.Lemma 1 (Doolan et al. [16] and O’Malley [17]).
Lety~x=y0+z0 be the zeroth-order asymptotic approximation to the solution of (5), where y0 represents the zeroth-order approximate outer solution (i.e., the solution of the reduced problem of (5)) and z0 represents the zeroth-order approximate solution in the boundary layer region of (5).
Then for a fixed positive integeri,(8)limh→0⁡yih≈y00+φ0-y00e-a0iρwhereρ=hε.Proof.
Lety0x be the solution of the reduced problem of (5)(9)pxy0′x+qxy0x=fx,y01=γ1and z0t is the solution of the boundary value problem (cf. O’Malley [17])(10)z0′′t+p0z0′t=0,z00=φ0-y00,z0∞=0wheret=xε.
From the theory of singular perturbation, the zeroth-order asymptotic approximation to the solution of (4) is (cf. O’Malley [17])(11)y~x≈y0x+a0axφ0-y00e-∫0xax/εdx.
As we are considering the differential equations on sufficiently small subintervals, the coefficients could be assumed to be locally constant. Hence,(12)y~x≈y0x+φ0-y00e-a0/εx.So, at the nodal points, we have(13)y~xi≈y0xi+φ0-y00ea0/εxi,i=0,1,…,N;that is,(14)y~ih≈y0ih+φ0-y00e-a0/εih.Therefore(15)limh→0⁡y~ih≈y00+φ0-y00exp-a0iρfori=0,1,2,3,…,N,where ρ=h/ε.
Now, we consider the difference scheme [18] by Galerkin method as follows:
Select a set of basis functionsφjx,j=0,1,…,N, which will define an interpolation scheme for the approximate solution over a grid of points a=x0<x1<⋯<xN+1=b. For simplicity, we use piecewise Lagrange polynomials lix of first degree as the basis functions. These interpolating polynomials are(16)l0ξ=ξ-ξ1ξ0-ξ1=1-ξ2,(17)l1ξ=ξ-ξ0ξ1-ξ0=1+ξ2in local element coordinates -1≤ξ≤1.TheN nodal values of the approximate solution y at the interior nodes x1,x2,…,xN are determined using this basis. The given boundary conditions determine the value of yx at the end nodes x0 and xN+1. The Galerkin method is now employed to obtain the integral equations; we have(18)εy′′x+axy′x+bxyx,φj=fx,φjforj=1,2,…,Nwhich is an integral equation ∫x0xNεy′′x+axy′x+bxyxϕjdx=∫x0xNfxφjdx.Sincey is sum of piecewise linear Lagrange polynomials, the second-order derivatives appearing in (17) vanish except at the element boundaries xi, where they become infinite.By integration by parts, (18) becomes(19)-εdydx,dϕjdx+axdydx+bxy,ϕj+εdydxϕjab=fx,ϕj.Using the substitution of trial function y(x)=ϕ0x+∑i=1Nyiϕjx into the integral equation (19), we have(20)∑i=1Nyiεdϕidx,dϕjdx-∑i=1Nyiaxdϕidx+bxϕi,ϕj=-αdl0dx,dϕjdx-βdlN+1dx,dϕjdx+αaxdl0dx+bxl0,ϕj+βaxdlN+1dx+bxlN+1,ϕj+εdydxϕjab-fx,ϕjfor j=1,2,…,N.It can be observed that all quantities on the right side of (20) can be computed from known boundary data to obtain N equations in the N unknown values yi at the interior nodes.The integrals in (20) can be solved by taking advantage of local coordinate ξ system.Since(21)x=b-a2ξ+b+a2=h2ξ+b+a2,ξ=2hx-b+a2⟹dξdx=2h,we have, by simple integration,(22)∫abdϕidxdϕjdxdx=∫-11dl0dξdξdxdl1dξdξdxh2dξ=-1h,fori=j-1,∫abdϕidxdϕjdxdx=∫-11dl1dξdξdxdl0dξdξdxh2dξ=-1h,fori=j+1,∫abdϕidxdϕjdxdx=∫-11dl1dξdξdxdl1dξdξdxh2dξ=1h,fori=j,∫abdϕidxdϕjdxdx=0fori-j>1.By assuming ax,bx, and fx as constants, the integral equation (20) gives, for a typical internal node j,(23)yj-1-εh+a2-bh6+yj2εh-2bh3+yj+1-εh-a2-bh6=-fjh.Equation (23), when rearranged, gives the following system of difference equations:(24)εyj+1-2yj+yj-1h2+ayj+1-yj-12h+byj-1+4yj+yj+16=fj;for1≤j≤N-1.Now, introduce a fitting factor σ in the Galerkin scheme as follows:(25)σεyj+1-2yj+yj-1h2+ajyj+1-yj-12h+bjyj-1+4yj+yj+16=fjfor1≤j≤N-1 with y0=ϕ0,yN=γ1. Here σ is a fitting factor which is to be determined in such a way that the solution of (25) converges uniformly to the solution of (5). Multiplying (25) by h and taking the limit as h→0 (in [16]), we get (26)limh→0⁡σρyj+1-2yj+yj-1+12ajyj+1-yj-1=0iffj-byjis  bounded(27)∴limh→0⁡σρyjh+h-2yjh+yjh-h+12ajyjh+h-yjh-h=0.Now, approximating the solutiony(x) by zeroth-order asymptotic approximation y~(x) and using Lemma 1, we have(28)limh→0⁡y~ih=y00+φ0-y00e-a0iρ,limh→0⁡y~jh+h=y00+φ0-y00e-a0jρe-a0ρ,limh→0⁡y~jh-h=y00+φ0-y00e-a0jρea0ρ.Using the above equations in (27), we get(29)σ=ρ2a0coth⁡a20-εb0a0ρ2.From (25), we have(30)yj-1εσh2-aj2h+bj6-yj2εσh2-2bj3+yj+1εσh2+aj2h+bj6=fjfor j=1,2,…,N-1.Equation (30) can be written as a three-term recurrence relation as follows:(31)Ejyj-1-Fjyi+Gjyj+1=Hj;j=1,2,…,N-1,where(32)Ej=εσh2-aj2h+bj6,Fj=2εσh2-2bj3,Gj=εσh2+aj2h+bj6,Hj=fj.The tridiagonal system equation (31) is solved using Thomas algorithm.


## 3.2. Right-End Layer Problems

We now discuss the method for singularly perturbed two-point boundary value problems with right-end boundary layer of the underlying interval. Assume thatax,bx, and fx are sufficiently continuously differentiable functions in 0,1. Furthermore, assume that ax≤M<0 in 0,1, where M is a negative constant. Under these assumptions, (5) has a unique solution yx which, in general, displays a boundary layer of width Oε at x=1.Lemma 2.
Letyx=y0+z0 be the zeroth-order asymptotic approximation to the solution of (5), where y0 represents the zeroth-order approximate outer solution and z0 represents the zeroth-order approximate solution in the boundary layer region.
Then for a fixed positive integeri,(33)limh→0⁡yih≈y00+β-y01expa11ε-iρwhereρ=hε.Proof.
The proof is based on asymptotic analysis (Doolan et al. [16] and O’Malley [17]) and is similar to the proof of Lemma 1.Applying the same procedure as in Section3 and using Lemma 2, we get the tridiagonal system equation (20) with fitting factor as(34)σ=ρ2a0coth⁡a21-εb1a1ρ2.


## 4. Convergence Analysis

Writing the tridiagonal system equation (31) in matrix-vector form, we get(35)AY=C,in which A=mij,1≤i,j≤N-1, is a tridiagonal matrix of order N-1, with(36)mii+1=σεh-ai2+hbi6,mii=-2σεh-2hbi3,mii-1=σεh-ai2+hbi6,and C=di is a column vector with di=hHi, where i=1,2,…,N-1 with local truncation error(37)τi=h2bi6yi′′+Oh3and Y=y0,y1,…,yNt.We also have(38)AY--Th=C,where Y-=y0-,y1-,…,yN-t denotes the actual solution and Th=T0(h0),T1(h1),…,TN(hN)t is the local truncation error.From (35) and (38), we get(39)AY--Y=Th.Thus, the error equation is(40)AE=Th,where E=Y--Y=e0,e1,…,eNt.Let theith row elements sum of matrix A be Si; then we have(41)Si=∑j=1N-1mij=-εσh+5hbi6+ai2,fori=1,Si=∑j=1N-1mij=∑j=1N-1mij=hbi=Bi0,fori=2,3,…,N-2,Si=∑j=1N-1mij=-σεh+5hbi6-ai2,fori=N-1.We can choose h sufficiently small so that the matrix A is irreducible and monotone. It follows that A-1 exists and its elements are nonnegative.Hence, from (40), we get(42)E=A-1Th,(43)E≤A-1·Th.Also from the theory of matrices we have(44)∑i=1N-1m-k,iSi=1,k=11N-1,where m-k,i is k,i element of the matrix A-1 for some i0 between 1 and N-1.Therefore,(45)∑i=1N-1m-k,i≤1min1≤i≤N-1⁡Si=1Bio≤1Bio,where Bio=hbi. We define A-1=max1≤k≤N-1⁡∑i=1N-1m-k,i and Th=max1≤i≤N-1⁡Tih.From (37), (40), (43), and (45), we get(46)ej=∑i=1N-1m-k,iTih,j=11N-1which implies(47)ej≤Oh2Bio,j=11N-1.Therefore,(48)E=Oh;that is, our method reduces to a first-order convergent for uniform mesh.


## 5. Numerical Examples

To demonstrate the applicability of the method, we have applied the method on four boundary value problems. These examples have been chosen because they have been widely discussed in literature and exact solutions are available for comparison.The exact solution of the boundary value problem(49)εy′′x+pxy′x+qxyx-δ+rxyx+sxyx+η=fxunder the boundary conditions(50)yx=φx,-δ≤x≤0,yx=γx,1≤x≤1+η,is(51)yx=c1em1x+c2em2x+fc,where c=q+r+s, c1=-f+γc+em2f-ϕc/cem1-em2, c2=f-γc+em1-f+ϕc/cem1-em2,(52)m1=-p-qδ+sη+p-qδ+sη2-4cε2ε,m2=-p-qδ+sη-p-qδ+sη2-4cε2ε.Example 1.
Consider the model boundary value problem of the type given by (1)-(2) having the boundary layer at the left end(53)εy′′+y′-2yx-δ-5y+yx+η=0with boundary conditions y(x)=1,-δ≤x≤0, and y(x)=1, 1≤x≤1+η.
The maximum absolute errors are given in Tables1 and 2 for different values of the delay and advanced parameters with perturbation parameter. The effect of the small parameters on the boundary layer solutions is shown in Figures 1 and 2.Table 1
The maximum absolute errors in solution of Example1.
ε↓ N = 8 N = 16 N = 32 N = 64 N = 128 N = 256 Present methodδ=0.5ε, η=0.5ε 1 0 - 1 0.03707062 0.01012487 0.00237565 0.0005848 0.00014565 3.6379e − 005 1 0 - 2 0.06876102 0.04407092 0.02342131 0.00865832 0.002129327 0.00050922 1 0 - 3 0.08556830 0.05231609 0.02808039 0.01392938 0.00720360 0.00356728 1 0 - 4 0.08740211 0.05431975 0.03006119 0.01590962 0.00810919 0.00401592 1 0 - 5 0.08758721 0.05452232 0.03026147 0.01610995 0.00830856 0.00421487 1 0 - 6 0.08760574 0.05454260 0.03028152 0.01613000 0.00832852 0.00423478  Results in Kadalbajoo and Sharma [10] 1 0 - 1 0.12011566 0.07181396 0.04482982 0.02694612 0.01516093 0.00775036 1 0 - 2 0.18727108 0.10697821 0.05904116 0.30796891 0.01567964 0.00799076 1 0 - 3 0.20429729 0.11915028 0.06879232 0.03655236 0.01893849 0.00963304 1 0 - 4 0.20614146 0.12048418 0.06989944 0.03721375 0.01932774 0.00984236 1 0 - 5 0.20632746 0.12061888 0.07001167 0.03728089 0.01936732 0.00986365 1 0 - 6 0.20634608 0.12063236 0.07002291 0.03728761 0.01937129 0.00986578Table 2
The maximum errors in solution of Example1 with ε=0.1.
N = 8 N = 32 N = 128 N = 512 Present method δ↓ η = 0.5 ε 0.00 0.03724038 0.00231155 0.00014202 8.8820e − 006 0.05 0.03707062 0.00237565 0.00014565 9.0928e − 006 0.09 0.03669100 0.00241210 0.00014764 9.2248e − 006 η↓ δ = 0.5 ε 0.00 0.03720017 0.00234619 0.00014399 8.9938e − 006 0.05 0.03707062 0.00237565 0.00014565 9.0928e − 006 0.09 0.03690625 0.00239551 0.00014675 9.1626e − 006  Results in Kadalbajoo and Sharma [10] δ↓ η = 0.5 ε 0.00 0.09190267 0.03453494 0.01164358 0.00300463 0.05 0.10233615 0.03823132 0.01295871 0.00335137 0.09 0.11018870 0.04110846 0.01400144 0.00362925 η↓ δ = 0.5 ε 0.00 0.09720079 0.03640446 0.01229476 0.00317786 0.05 0.10233615 0.03823132 0.01295871 0.00335137 0.09 0.10632014 0.03965833 0.01348348 0.00349050Figure 1
Numerical solution of Example1 for different values of δ with ε=0.1,η=0.05.Figure 2
Numerical solution of Example1 for different values of η with ε=0.1,δ=0.05.Example 2.
Consider the boundary value problem having the boundary layer at the left end(54)εy′′+y′-2yx-δ+y-yx+η=-1with boundary conditions yx=1,-δ≤x≤0, and yx=1,1≤x≤1+η.
The maximum absolute errors are given in Table3 for different values of ε with the delay and advance parameter values δ=0.5ε,η=0.5ε. The effect of the small parameters on the boundary layer solutions is shown in Figures 3 and 4.Table 3
The maximum errors in solution of Example2.
ε↓ N = 8 N = 16 N = 32 N = 64 N = 128 N = 256 Present methodδ=0.5ε, η=0.5ε 1 0 - 1 0.00470784 0.00123359 0.00029344 7.2455e − 005 1.8125e − 005 4.5281e − 006 1 0 - 2 0.01585466 0.00697845 0.00240001 0.00093287 0.00024806 5.8741e − 005 1 0 - 3 0.01910533 0.01015879 0.00512445 0.00244082 0.00105444 0.00037664 1 0 - 4 0.01944288 0.01049291 0.00545649 0.00277170 0.00138358 0.00067752 1 0 - 5 0.01947677 0.01052645 0.00548983 0.00280492 0.00141674 0.00071065 1 0 - 6 0.01948016 0.01052981 0.00549316 0.00280824 0.00142006 0.000713967  Results in Kadalbajoo and Sharma [10] 1 0 - 1 0.08579690 0.05129568 0.03202130 0.01924723 0.01098354 0.00553597 1 0 - 2 0.13376506 0.07641301 0.04217226 0.02199778 0.01119974 0.00570769 1 0 - 3 0.14592663 0.08510734 0.04913737 0.02610883 0.01352749 0.00688074 1 0 - 4 0.14724390 0.08606013 0.04992817 0.02658125 0.01380553 0.00703026 1 0 - 5 0.14737676 0.08615634 0.05000834 0.02662921 0.01383380 0.00704546 1 0 - 6 0.14739006 0.08616597 0.05001637 0.02663401 0.01383663 0.00704699Figure 3
Numerical solution of Example2 for different values of δ with ε=0.1,η=0.05.Figure 4
Numerical solution of Example2 for different values of η with ε=0.1,δ=0.05.Example 3.
Consider the boundary value problem of the type given by (1)-(2) having the boundary layer at the right end (55)εy′′-y′-2yx-δ+y-2yx+η=0with boundary conditions yx=1,-δ≤x≤0, and yx=-1,1≤x≤1+η.
The maximum absolute errors are given in Table4 with ε=0.1 for different values of the delay and advance parameters. The effect of the small parameters on the boundary layer solutions is shown in Figures 5 and 6.Table 4
The maximum errors in solution of Example3 with ε=0.1.
N = 8 N = 32 N = 128 N = 512 Present method δ↓ η = 0.5 ε 0.00 0.031377538 0.001800241 0.000112071 7.0036e − 006 0.05 0.029748010 0.001700026 0.000105418 6.5860e − 006 0.09 0.028294285 0.001611053 9.9793e − 005 6.2344e − 006 η↓ δ = 0.5 ε 0.00 0.027910529 0.001587651 9.8361e − 005 6.1442e − 006 0.05 0.029748010 0.001700026 0.000105418 6.5860e − 006 0.09 0.031068500 0.001781207 0.000110800 6.9223e − 006  Results in Kadalbajoo and Sharma [10] δ↓ η = 0.5 ε 0.00 0.09930002 0.03685072 0.01331683 0.00342882 0.05 0.09997296 0.03218424 0.01167102 0.00299572 0.09 0.10044578 0.02850398 0.01038902 0.00266379 η↓ δ = 0.5 ε 0.00 0.10055269 0.02759534 0.01007834 0.00258299 0.05 0.09997296 0.03218424 0.01167102 0.00299572 0.09 0.09944067 0.03591410 0.01297367 0.00334044Figure 5
Numerical solution of Example3 for different values of η with ε=0.1,η=0.05.Figure 6
Numerical solution of Example3 for different values of η with ε=0.1,δ=0.05.Example 4.
Consider the boundary value problem having the boundary layer at the left end(56)εy′′+0.5y′-3yx-δ-2y+2yx+η=1with boundary conditions yx=1,-δ≤x≤0, and yx=0,1≤x≤1+η.
The maximum absolute errors are given in Table5 and Table 6 for different values of delay and advanced parameters with perturbation parameter. The effect of the small parameters on the boundary layer solutions is shown in Figures 7 and 8.Table 5
The maximum errors in solution of Example4.
ε↓ N = 8 N = 16 N = 32 N = 64 N = 128 N = 256 δ = 0.5 ε, η=0.5ε 1 0 - 1 0.02753382 0.00616170 0.00156771 0.00038927 9.7155e − 005 2.4278e − 005 1 0 - 2 0.09909124 0.05582849 0.02171656 0.00554645 0.00131747 0.00032522 1 0 - 3 0.10181648 0.06264442 0.03547474 0.01902508 0.00950250 0.00372083 1 0 - 4 0.10207612 0.06281638 0.03556178 0.01906740 0.00989895 0.00504914 1 0 - 5 0.10210236 0.06283382 0.03557070 0.01907016 0.00989794 0.00504604 1 0 - 6 0.10210499 0.06283557 0.03557159 0.01907044 0.00989798 0.00504673Table 6
The maximum errors in solution of Example4 with ε=0.1.
N = 8 N = 32 N = 128 N = 512 δ↓ η = 0.5 ε 0.00 0.025347510 0.001425327 8.9204e − 005 5.5742e − 006 0.05 0.027533826 0.001567710 9.7155e − 005 6.0690e − 006 0.09 0.028669770 0.001645550 0.000102186 6.3826e − 006  η↓ δ = 0.5 ε 0.00 0.026174618 0.001478341 9.2083e − 005 5.7527e − 006 0.05 0.027533826 0.001567710 9.7155e − 005 6.0690e − 006 0.09 0.028348272 0.001623113 0.00010057 6.2854e − 006Figure 7
Numerical solution of Example4 for different values of η with ε=0.1,η=0.05.Figure 8
Numerical solution of Example4 for different values of η with ε=0.1,δ=0.05.


## 6. Discussion and Conclusion

An exponentially fitted Galerkin method has been presented for solving singularly perturbed differential-difference equations with delay as well as advance parameters. To demonstrate the applicability of the method, three examples with left-end and one with right-end boundary layer have been solved for different values of the delay, advance, and perturbation parameters. The numerical results are taken by using MATLAB coding and solutions have been compared with the exact solutions and maximum absolute errors are presented in tables. To show the efficiency of the method, we have compared results of the proposed scheme with the results of Kadalbajoo and Sharma [10]. The rate of convergence in the examples is given in Table 7. It is observed that the present method approximates the exact solution very well for which other classical finite difference methods fail to give good results. The effect of the delay and advance parameters on the solutions has also been investigated and presented by using graphs. When the solution of the boundary value problem exhibits layer behaviour on the left side, the effect of delay or advance on the solution in the boundary layer region is negligible while, in the outer region, it is considerable and the change in the advance affects the solution in similar fashion; that is, the width of the layer in the outer region increases as the advance parameter increases and decreases as the delay increases, but reversely (Figures 1–4). When the solution of the boundary value problem exhibits layer behaviour on the right side, the changes in delay or advance affect the solution in boundary layer region as well as outer region. The thickness of the layer increases as the size of the delay increases while it decreases as the size of the advance increases (Figures 5 and 6).Table 7
The rate of convergence of the examples withδ=0.5ε and η=0.5ε.
ε↓ N = 8 N = 16 N = 32 N = 64 N = 128 Example1 1 0 - 2 0.6418 0.9120 1.4357 2.0237 2.0640 1 0 - 3 0.7098 0.8977 1.0114 0.9513 1.0139 1 0 - 4 0.6862 0.8536 0.9180 0.9723 1.0138 1 0 - 5 0.6839 0.8494 0.9095 0.9553 0.9791 1 0 - 6 0.6836 0.8489 0.9087 0.9536 0.9758  Example2 1 0 - 2 1.1839 1.5399 1.3633 1.9110 2.0782 1 0 - 3 0.9112 0.9873 1.0700 1.2109 1.4852 1 0 - 4 0.8898 0.9434 0.9772 1.0024 1.0301 1 0 - 5 0.8877 0.9392 0.9688 0.9854 0.9954 1 0 - 6 0.8875 0.9388 0.9680 0.9837 0.9920  Example3 1 0 - 2 1.1839 1.5398 1.3633 1.9112 2.0779 1 0 - 3 0.9112 0.9872 1.0700 1.2109 1.4853 1 0 - 4 0.8898 0.9433 0.9771 1.0024 1.0300 1 0 - 5 0.8877 0.9391 0.9688 0.9854 0.9954 1 0 - 6 0.8875 0.9387 0.9679 0.9837 0.9920  Example4 1 0 - 2 0.8278 1.3622 1.9692 2.0738 2.0183 1 0 - 3 0.7007 0.8204 0.8989 1.0015 1.3527 1 0 - 4 0.7004 0.8208 0.8992 0.9458 0.9712 1 0 - 5 0.7004 0.8209 0.8994 0.9461 0.9720 1 0 - 6 0.7004 0.8209 0.8994 0.9461 0.9718



---

*Source: 1935853-2016-10-16.xml*