# Exponentially Fitted Numerical Method for Singularly Perturbed Differential-Difference Equations

Authors: Habtamu Garoma Debela; Solomon Bati Kejela; Ayana Deressa Negassa
Journal: International Journal of Differential Equations
 (2020)

---

## Abstract

This paper presents a numerical method to solve singularly perturbed differential-difference equations. The solution of this problem exhibits layer or oscillatory behavior depending on the sign of the sum of the coefficients in reaction terms. A fourth-order exponentially fitted numerical scheme on uniform mesh is developed. The stability and convergence of the proposed method have been established. The effect of delay parameter (small shift) on the boundary layer(s) has also been analyzed and depicted in graphs. The applicability of the proposed scheme is validated by implementing it on four model examples. Maximum absolute errors in comparison with the other numerical experiments are tabulated to illustrate the proposed method.


---

## Body

## 1. Introduction

Differential equations with a small parameterε0<ε≤1 multiplying the highest order derivative are called singularly perturbed differential equations. Mathematically, any ordinary differential equation in which the highest derivative is multiplied by a small positive parameter and containing at least one delay/advance parameter is known as a singularly perturbed differential-difference equation [1]. Such type of equations arises frequently from the mathematical modeling of various practical phenomena, for example, in the modeling of the human pupil-light reflex [2], the study of bistable devices [3], and vibrational problems in control theory [4]. When perturbation parameter ε is very small, most numerical methods for solving such problems may be unstable and give inaccurate results. So, it is important to develop suitable numerical methods to solve singularly perturbed delay differential equations.Hence, in the recent times, many researchers have been trying to develop numerical methods for solving singularly perturbed delay differential equations. For example, the authors in [5] proposed a computational method of first order for singularly perturbed delay reaction-diffusion equations with layer or oscillatory behavior. The authors in [6] presented a fourth-order finite difference scheme for second-order singularly perturbed differential-difference equations with negative shift. The authors in [7] presented exponentially fitted a second-order finite difference scheme for a class of singularly perturbed delay differential equations with large delay. In [8], the numerical solution of singularly perturbed differential-difference equations with dual layer was presented. Recently, the authors in [9] presented a computational method for solving a singularly perturbed delay differential equation with twin layers or oscillatory behavior. But, still there is a lack of accuracy because the treatment of singularly perturbed problems is not trivial and the solution depends on perturbation parameter ε and mesh size h [10–12]. Due to this, numerical treatment of singularly perturbed delay differential equations needs improvement. Therefore, it is important to develop a more accurate and convergent numerical method for solving singularly perturbed delay differential equations.Consider the following singularly perturbed delay differential equation of the form:(1)εy″+axyx−δ+bxyx=fx,0<x<1,subject to the interval and boundary conditions(2)yx=ϕx,for−δ≤x≤0andy1=φ,where ε is the perturbation parameter, 0<ε≪1, and δ is a small delay parameter of Oε,0<δ≪1. Also, ax,bx,fx,andϕx are bounded smooth functions and φ is a given constant. The layer or oscillatory behavior of the problem under consideration is maintained for δ≠0 but sufficiently small, depending on the sign of ax+bx, for all x∈0,1. If ax+bx<0, the solution of the problem in equations (1) and (2) exhibits layer behavior, and if ax+bx<0, it exhibits oscillatory behavior. Therefore, if the solution exhibits layer behavior, there will be two boundary layers which will occur at both end points x=0 and x=1[12].Thus, the purpose of this study is to develop stable, convergent, and more accurate numerical method for solving singularly perturbed delay differential equations.


## 2. Description of the Method

By using Taylor series expansion in the neighborhood ofx, we have(3)yx−δ≈yx−δy′x+Oδ2.Substituting equation (3) into equation (1), we obtain an asymptotically equivalent singularly perturbed two-point boundary value problem:(4)εy″x+pxy′x+qxyx=fx,where px=−δaxandqx=ax+bx under boundary conditions(5)y0=ϕ0,y1=φ.Discretizing the given interval [0, 1] in toN equal parts with constant mesh size h, we have xi=x0+ih,i=0,1,2,…,N.Using Taylor’s expansions ofyi+1 and yi−1 up to Oh5, we get the finite difference approximations for yi′ and yi″ as(6)yi′=yi+1−yi−12h−h26yi″′+τ1,where τ1=−h4/20y5ξ1, for ξ1∈xi−1,xi, and(7)yi″=yi+1−2yi+yi−1h2−h212yi4+τ2,where τ2=−h4/360y6ξ2, for ξ2∈xi−1,xi.Substituting equations (6) and (7) into equation (4) and simplifying, we obtain(8)εh2yi+1−2yi+yi−1+pi2hyi+1−yi−1−h26piyi−εh212yi4+qiyi=fi+τi,where τi=h4pi/120y5ξ1+ε/360y6ξ2+Oh5 is the local truncation error and pxi=pi,qxi=qi,fxi=fi, and yi≈yxi.By successively differentiating both sides of equation (4) and evaluating at xi and using into equation (8), we obtain(9)εh2yi+1−2yi+yi−1+pi2hyi+1−yi−1+εh2pi26ε+h262pi′+qi−pi2εyi″+Aiyi′+Bi+biyi=fi+Ci,where(10)Ai=h26εpipi′+qi−h212piεpi′+qi−pi″−2qi′,Bi=h26εpiqi′−h212piqi′ε−qi″,Ci=fi+aih212εfi′+h212fi″+τi.Now introducing a fitting parameterσ and using central difference approximation for yi″ and yi′ in equation (9), we have(11)σεh21+h2pi212ε+h2122pi′+qiyi+1−2yi+yi−1+12hpi+Aiyi+1−yi−1+Bi+qiyi=fi+Ci.Multiplying both sides of equation (11) by h and taking the limit as h⟶0, we obtain(12)σ12ρ12+ρ2pi2limh⟶0yi+1−2yi+yi−1+pi2limh⟶0yi+1−yi−1=0,where ρ=h/ε.From the theory of singular perturbations and [13], we have two cases for px>0 and px<0:Case 1.
Forpx<0 (right-end boundary layer), we have(13)limh⟶0yi+1−2yi+yi−1=ϕ0−y00e−p0−1/ε+iρe−p0ρ+ep0ρ−2,limh⟶0yi+1−yi−1=ϕ0−y00e−p0−1/ε+iρe−p0ρ−ep0ρ.
Thus, from equation (12), we get(14)σ0=6ρp012+ρ2p20cothp0ρ2.Case 2.
Forpx>0 (left-end boundary layer), we have(15)limh⟶0yi+1−2yi+yi−1=φ0−y01e−p1−1/ε+iρe−p1ρ+ep1ρ−2,limh⟶0yi+1−yi−1=φ0−y01e−p1−1/ε+iρe−p1ρ−ep1ρ.
Thus, from equation (12), we get(16)σ1=6ρp112+ρ2p21cothp1ρ2.
In general, for discretization, we take a variable fitting parameter as(17)σi=6ρipi12+ρipi2cothpiρi2,where ρi=h/ε.
Simplifying equation (11), we get the tridiagonal system of the equation as follows:(18)LN≡Eiyi−1−Fiyi+Giyi+1=Hi,fori=1,2,...,N−1,where(19)Ei=σiεh2+σipi212+σiε122pi′+qi−12hpi+Ai,Fi=2σiεh2+σipi26+σiε62pi′+qi−pi−Bi,Gi=σiεh2+σipi212+σiε122pi′+qi+12hpi+Ai,Hi=fi+Ci.
The tridiagonal system in equation (18) can be easily solved by the Thomas algorithm with help of Matlab 2013.


## 3. Stability and Convergence Analysis

Case1: Layer behavior (i.e., ax+bx<0, for x∈0,1; thus, qx<0, since ε>0).Lemma 1.
Ify0≥0 and Lyx≤0, for all x∈0,1, then the solution yx≥0 for all x∈0,1 for equations (4) and (5).Proof.
Supposet∈0,1, such that yt=minx∈0,1yx and yt<0. Since t∉0,1 and is a point of minima, then yt=0 and y″t≥0.
Therefore, we haveLyt=εy″t+pty′t+qtyt>0, since yt<0 (by assumption) and qt<0. But this is a contradiction. Then, it follows that yx≥0 for x∈0,1.Theorem 1.
If the solution of the problem in equations (4) and (5) satisfies(20)yx≤Kmaxy0,maxx∈0,1Lyx,for some constant K≥1, then the solution is stable.Proof.
We define two functions:ψ±x=Kmaxy0,maxx∈0,1Lyx±yx. Then, we get ψ±0≥0 and(21)Lψ±x=Kqxmaxy0,maxx∈0,1Lyx±Lyx≤0,since qx<0 and for suitable choice of K. Therefore, by Lemma 1, we get ψ±x≥0, for all x∈0,1. So, yx≤Kmaxy0,maxx∈0,1Lyx.
Hence, the stability of the solutions of the problem in equations (4) and (5) is proved for the case of the layer behavior.Lemma 2.
The finite difference operatorLN in equation (13) satisfies the discrete minimum principle; i.e., if wi is any mesh function such that w0≥0 and LNwi≤0, for all xi∈0,1, then wi≥0 for all xi∈0,1.Proof.
Suppose there exists a positive integerk such that wk<0andwi=min0≤i≤Nwi. Then, from equation (13), we have(22)LNwk=Ekwk−1−Fkwk+Gkwk+1=σkεh2+σkpk2112+σkε122pk′+qkwk+1−2wk+wk−1+12hpk+Akwk+1−wk−1+pk+Bkwk=σkεh2+σkpk2112+σkε122pk′+qkwi+1−wk+σkεh2+σkpk2112+σkε122pk′+qkwk−1−wk+12hpk+h26εpkpk′+qk−h212pkεpk′+qk−pk″−2qk′wi+1−wi−1+pk+h26εpiqi′−h212piqi′ε−qk″wk.
For sufficiently smallh and for suitable value of pk, we obtain LNwk>0. However, wk<0 (by the assumption) and Bk⟶pk<0. But, this is a contradiction. Hence, wi≥0for allxi∈0,1.Theorem 2.
The finite difference operatorLN in equation (13) is stable for ax+bx<0, if wi is any mesh function, then wi≤Kmaxw0,maxxi∈0,1Lwi, for some constant K≥1.Proof.
We define two functions:(23)ψ±=Kmaxw0,maxx∈0,1Lwi±wi.
Then, similar to Theorem1, we get ψ±0≥0 and(24)Lψ±=Kqimaxw0,maxx∈0,1Lwi±Lwi≤0,since ai+bi<0⇒qx<0 and K≥1, therefore, by Lemma 2, we get ψi±≥0, for all xi∈0,1. ⇒ψi±≡Cmaxw0,maxxi∈0,1Lwi±wi≥0.
Thus,wi≤Cmaxw0,maxxi∈0,1Lwi.
This proves the stability of the scheme for the case of layer behavior.
Case2: oscillatory behavior (i.e., ax+bx>0, for x∈0,1; thus, qx>0, since ε>0).
The continuous maximum principle and stability of the solution of equations (4) and (5) are presented as follows for the case of oscillatory behavior.Lemma 3.
Ify0≥0 and Lyx≥0, for all x∈0,1, then the solution yx≥0 for all x∈0,1 for equations (4) and (5).Proof.
Supposet∈0,1, such that yt=maxx∈0,1yx and yt<0. Since t∉0,1 and is a point of maxima, therefore, y′t=0 and y″t≤0. Therefore, we have Lyt=εy″t+pty′t+qtyt<0, since yt<0 (by assumption) and qt<0. But this is a contradiction. Hence, yx≥0 for x∈0,1.Theorem 3.
If the solution of the problem in equations (4) and (5) satisfies(25)yx≤Kmaxy0,maxx∈0,1Lyx,for some constant, k≥1.Proof.
The proof is analogous to Theorem1. Hence, the stability of the solutions of the problem in equations (4) and (5) is proved for the case of oscillatory behavior. Now, we present the stability of the discrete problem in equation (13) for the case of oscillatory behavior.Lemma 4.
The finite difference operatorLN in equation (13) satisfies the discrete maximum principle, if wi is any mesh function such that w0≥0 and LNwi≥0, for all xi∈0,1, then wi≥0 for all xi∈0,1.Proof.
Suppose there exists a positive integerk such that wk<0andwi=max0≤i≤Nwi. Then, from equation (13), we have(26)LNwk=Ekwk−1−Fkwk+Gkwk+1=σkεh2+σkpk2112+σkε122pk′+qkwk+1−2wk+wk−1+12hpk+Akwk+1−wk−1+pk+Bkwk=σkεh2+σkpk2112+σkε122pk′+qkwi+1−wk+σkεh2+σkpk2112+σkε122pk′+qkwk−1−wk+12hpk+h26εpkpk′+qk−h212pkεpk′+qk−pk″−2qk′wi+1−wi−1+pk+h26εpiqk′−h212piqi′ε−qk″wk.
For sufficiently smallh and for suitable value of pk, we obtain LNwk<0. However, wk<0(by the assumption) and Bk⟶pk>0. But, this is a contradiction. Hence, wi≥0for allxi∈0,1.Theorem 4.
The finite difference operatorLN in equation (13) is stable for ax+bx>0,i.e.qx>0, if wi is any mesh function, then wi≤Kmaxw0,maxxi∈0,1Lwi, for some constant K≥1.Proof.
The proof is similar to Theorem2. This proves the scheme for the case of oscillatory behavior.Definition 1.
(uniform convergence). Lety be a solution of equations (1) and (2). Consider a difference scheme for solving equations (1) and (2). If the scheme has a numerical solution yN that satisfies y−yN≤Chp, where C>0 and p>0 are independent of ε and of mesh size h, then we say the scheme uniformly convergent to y with respect to the norm ⋅ [14].Lemma 5.
The bound for derivative of the solutionyx of the problem (1)–(3) when x∈Ω1=0,1 is given by(27)ykx≤C1+ε−kexp−αxjε,0≤k≤4,j=1,2,…,N−1,ykx≤C1+ε−kexp−α1−xjε,0≤k≤4,j=1,2,…,N−1.Proof.
For the proof, refer [15]
For any mesh functionzi, define the following difference operators:(28)D+zi=zi+1−zih,D−zi=zi−zi−1h,D0zi=zi+1−zi−12h,D+D−zi=zi+1−2zi+zi−1h2.Theorem 5.
Letyxi and yi be, respectively, the analytical solution of equations (4) and (5) and numerical solutions of equation (13). Then, for sufficiently large N, the following parameter uniform error estimate holds:(29)sup0<ε≤1yxi−yi≤C∗N−4.Proof.
Let us consider the local truncation error defined as(30)LNyxi−yi=εσρd2dx2−D+D−+h212d4dx4yxi+pxiddx−D0+h26d2dx2yxi,where(31)εσρ=ε6ρp012+ρ2p20cothp0ρ2≈p0N−12cothp0N−12ε,since ρ=N−1/ε. In our assumption, ε≤h=N−1.
By consideringN is fixed and taking the limit for ε⟶0, we obtain the following:(32)limε⟶0εσρ=limε⟶0p0N−12cothp0N−12ε=CN−1.
From Taylor series expansion, the bound for the difference becomes(33)d2dx2−D+D−yxi≤CN−3d4yxidx4,ddx−D0yxi≤CN−2d3yxidx3,where(34)dkyxidxk=supxi∈x0,xNdkyxidxk,k=3,4.
Now using the bounds and the assumptionε≤N−1, equation (30) reduces to(35)LNyxi−yi=εσρd2dx2−D+D−yxi+εσρh212yi4+pxiddx−D0yxi+h26yi2≤εσρd2dx2−D+D−yxi+εσρh212yi4+pxiddx−D0yxi+h26yi2≤CN−3d4yxidx4+CN−2d3yxidx3+CN−2d2yxidx2.
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 right layer solution in Lemma5, we obtain(36)LNyxi−yi≤CN−3d4yxidx4+CN−2d3yxidx3+CN−2d2yxidx2≤CN−31+ε−4exp−α1−xjε+CN−21+ε−3exp−α1−xjε+CN−21+ε−2exp−α1−xjε≤CN−41+ε−4exp−α1−xjε,sinceε−4≥ε−3≥ε−2.Lemma 6.
For a fixed mesh and forε⟶0, it holds(37)limε⟶0max1≤j≤N−1exp−α1−xj/εεm=0,m=1,2,3,….Proof.
Refer from [16].
By using Lemma6 into equation (36), it results in(38)Lhyxi−yi≤CN−4.
Hence, by discrete maximum principle, we obtain(39)yxi−yi≤CN−4.
Thus, result of equation (39) shows equation (29). Hence proved.Remark 1.
A similar analysis for convergence may be carried out for the left layer case.


## 4. Numerical Examples and Results

To demonstrate the applicability of the method, we implemented the method on four numerical examples, two with twin boundary layers and two with oscillatory behavior. Since those examples have no exact solution, the numerical solutions are computed using the double mesh principle. The maximum absolute errors are computed using the double mesh principle given by(40)Zh=maxiyih−yih/2,i=1,2,…,N−1,where yih is the numerical solution on the mesh xi1N at the nodal point xi and xi=x0+ih,i=1,2,…,N−1 and yih/2 is the numerical solution on a mesh, obtained by bisecting the original mesh with N number of intervals [9].Example 1.
Consider the singularly perturbed delay reaction-diffusion equation with layer behavior(41)εy″x+0.25yx−δ−yx=1,under the interval and boundary conditions(42)yx=1,−δ≤x≤0,y1=0.
The maximum absolute errors are presented in Tables1 and 2 for different values of ε and δ.
The graph of the computed solution forε=0.01 and different values of δ is also given in Figure 1.Table 1
The maximum absolute errors of Example1 for different values of δ with ε=0.1.
δ↓N=100N=200N=300N=400N=500Present method0.031.1784e ‒ 097.3658e ‒ 111.4572e ‒ 114.9336e ‒ 121.9272e ‒ 120.051.1796e ‒ 097.3727e ‒ 111.4279e ‒ 114.8657e ‒ 122.3070e ‒ 120.091.1790e ‒ 097.3690e ‒ 111.4548e ‒ 114.5787e ‒ 122.5597e ‒ 12Method in [17]0.032.1999e ‒ 031.1041e ‒ 037.3705e ‒ 045.5315e ‒ 044.4269e ‒ 040.052.2012e ‒ 031.1049e ‒ 037.3749e ‒ 045.5345e ‒ 044.4293e ‒ 040.092.1999e ‒ 031.1038e ‒ 037.3676e ‒ 047.3676e ‒ 044.4247e ‒ 04Table 2
The maximum absolute errors of Example1 for different values of ε with δ=0.5ε.
εN=16N=32N=64N=128N=256Present method2−44.5775e ‒ 062.8651e ‒ 071.7913e ‒ 081.1197e ‒ 096.9995e ‒ 112−51.6246e ‒ 051.0190e ‒ 066.3830e ‒ 083.9901e ‒ 092.4940e ‒ 102−65.9281e ‒ 053.7757e ‒ 062.3632e ‒ 071.4791e ‒ 089.2453e ‒ 102−72.2949e ‒ 041.4731e ‒ 059.2549e ‒ 075.7989e ‒ 083.6250e ‒ 092−89.1215e ‒ 045.8144e ‒ 053.6824e ‒ 062.3104e ‒ 071.4448e ‒ 082−93.5308e ‒ 032.2815e ‒ 041.4669e ‒ 059.2088e ‒ 075.7724e ‒ 082−101.1709e ‒ 029.1043e ‒ 045.8034e ‒ 053.6736e ‒ 062.3055e ‒ 07Method in [17]2−41.8632e ‒ 029.6189e ‒ 034.8865e ‒ 032.4643e ‒ 031.2376e ‒ 032−52.8161e ‒ 021.4818e ‒ 027.6255e ‒ 033.8713e ‒ 031.9509e ‒ 032−63.7958e ‒ 022.0967e ‒ 021.0977e ‒ 025.6273e ‒ 032.8498e ‒ 032−75.0640e ‒ 022.8316e ‒ 021.5267e ‒ 027.9105e ‒ 034.0287e ‒ 032−86.3580e ‒ 023.7706e ‒ 022.0984e ‒ 021.1012e ‒ 025.6555e ‒ 032−98.3843e ‒ 025.0477e ‒ 022.8297e ‒ 021.5261e ‒ 027.9111e ‒ 032−109.9137e ‒ 026.3529e ‒ 023.7660e ‒ 022.0974e ‒ 021.1011e ‒ 02Figure 1
The numerical solution of Example1 with ε=0.01 and N = 100.Example 2.
Consider the singularly perturbed delay reaction-diffusion equation with layer behavior(43)εy″x−2yx−δ−yx=1under the interval and boundary conditions(44)yx=1,−δ≤x≤0,y1=0.
The maximum absolute errors are presented in Tables3 and 4 for different values of ε and δ.
The graph of the computed solution forε=0.01 and different values of δ is also given in Figure 2.Table 3
The maximum absolute errors of Example2 for different values of δ with ε=0.1.
δ↓N=100N=200N=300N=400N=500Present method0.038.2805e ‒ 095.1758e ‒ 101.0224e ‒ 103.2348e ‒ 111.3243e ‒ 110.057.7093e ‒ 094.8202e ‒ 109.5110e ‒ 113.0120e ‒ 111.2325e ‒ 110.096.0459e ‒ 093.7806e ‒ 107.4682e ‒ 112.3630e ‒ 119.6788e ‒ 12Method in [17]0.033.1674e ‒ 031.6058e ‒ 031.0754e ‒ 038.0837e ‒ 046.4760e ‒ 040.053.1437e ‒ 031.5949e ‒ 031.0685e ‒ 038.0338e ‒ 046.4367e ‒ 040.093.0784e ‒ 031.5660e ‒ 031.0502e ‒ 037.9000e ‒ 046.3310e ‒ 04Table 4
The maximum absolute errors of Example2 for different values of ε with δ=0.5ε.
εN=16N=32N=64N=128N=256Present method2−43.0267e ‒ 051.9031e ‒ 061.1950e ‒ 077.4728e ‒ 094.6725e ‒ 102−51.1987e ‒ 047.8382e ‒ 064.9134e ‒ 073.0732e ‒ 081.9223e ‒ 092−64.9863e ‒ 043.1795e ‒ 051.9986e ‒ 061.2579e ‒ 077.8649e ‒ 092−71.9386e ‒ 031.2530e ‒ 048.1293e ‒ 065.0955e ‒ 073.1919e ‒ 082−86.4424e ‒ 035.1006e ‒ 043.2516e ‒ 052.0492e ‒ 061.2889e ‒ 072−91.7543e ‒ 021.9764e ‒ 031.2772e ‒ 048.2523e ‒ 065.1725e ‒ 072−103.8002e ‒ 026.5572e ‒ 035.1492e ‒ 043.2824e ‒ 052.0727e ‒ 06Method in [17]2−42.1118e ‒ 021.1692e ‒ 026.1941e ‒ 033.1887e ‒ 031.6178e ‒ 032−52.7872e ‒ 021.6023e ‒ 028.6367e ‒ 034.4957e ‒ 032.2948e ‒ 032−63.5711e ‒ 022.1293e ‒ 021.1869e ‒ 026.2731e ‒ 033.2240e ‒ 032−74.6679e ‒ 022.8350e ‒ 021.6107e ‒ 028.6728e ‒ 034.5120e ‒ 032−85.4895e ‒ 023.6018e ‒ 022.1373e ‒ 021.1929e ‒ 026.2847e ‒ 032−95.7371e ‒ 024.7254e ‒ 022.8581e ‒ 021.6140e ‒ 028.6961e ‒ 032−105.7878e ‒ 025.5695e ‒ 023.6153e ‒ 022.1406e ‒ 021.1956e ‒ 02Figure 2
The numerical solution of Example2 with ε=0.01 and N = 100.Example 3.
Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behavior(45)εy″x+0.25yx−δ+yx=1,under the interval and boundary conditions(46)yx=1,−δ≤x≤0,y1=0.
The maximum absolute errors are presented in Table5 for different values of δ.
The graph of the computed solution forε=0.001 and δ = 0.003 is also given in Figure 3.Table 5
The maximum absolute errors of Example3 for different values of δ with ε=0.1.
δ↓N=100N=200N=300N=400N=500Present method0.034.0323e ‒ 082.5200e ‒ 094.9716e ‒ 101.5799e ‒ 106.3758e ‒ 110.053.9610e ‒ 082.4762e ‒ 094.8871e ‒ 101.5567e ‒ 107.2701e ‒ 110.093.8377e ‒ 082.3991e ‒ 094.7227e ‒ 101.5806e ‒ 106.1942e ‒ 11Method in [17]0.032.5991e ‒ 031.2872e ‒ 038.5528e ‒ 046.4039e ‒ 045.1179e ‒ 040.052.6270e ‒ 031.3013e ‒ 038.6474e ‒ 046.4750e ‒ 045.1749e ‒ 040.092.6813e ‒ 031.3289e ‒ 038.8320e ‒ 046.6139e ‒ 045.2863e ‒ 04Figure 3
The numerical solution of Example3 with ε=0.001 and N = 100.Example 4.
Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behavior(47)εy″x+yx−δ+2yx=1,under the interval and boundary conditions(48)yx=1,−δ≤x≤0,y1=0.
The maximum absolute errors are presented in Table6 for different values of δ.
The graph of the computed solution forε=0.001 and δ = 0.003 is also given in Figure 4.
The rate of convergence (ρ): in the same way, in equation (40), one can define zh/2 by replacing h by h/2 and N−1 by 2N−1, that is,(49)zh/2=maxiyih/2−yih/4,fori=1,2,…,2N−1.
The computational rate of convergenceρ is also obtained by using double mesh principle defined as follows [9]:(50)ρ=logzh−logzh/2log2.
The following tables (Tables7 and 8) show the rate of convergence ρ of the present method for different values of the mesh size h.Table 6
The maximum absolute errors of Example4 for different values of δ with ε=0.1.
δ↓N=100N=200N=300N=400N=500Present method0.031.5160e ‒ 079.4743e ‒ 091.8715e ‒ 095.9186e ‒ 102.4306e ‒ 100.051.5697e ‒ 079.8097e ‒ 091.9379e ‒ 096.1294e ‒ 102.5132e ‒ 100.091.7120e ‒ 071.0702e ‒ 082.1140e ‒ 096.6867e ‒ 102.7376e ‒ 10Method in [17]0.031.5929e ‒ 027.4850e ‒ 034.8816e ‒ 033.6202e ‒ 032.8764e ‒ 030.051.5470e ‒ 027.2782e ‒ 034.7473e ‒ 033.5209e ‒ 032.7975e ‒ 030.092.1396e ‒ 021.0097e ‒ 026.5922e ‒ 034.8916e ‒ 033.8879e ‒ 03Figure 4
The numerical solution of Example4 with ε=0.001 and N = 100.Table 7
Rate of convergenceρ for ε=0.1 and δ=0.05.
hh/2Zhh/4Zh/2ρExample11/1001/2001.1796e ‒ 091/4007.3727e ‒ 114.00001/2001/4007.3727e ‒ 111/8004.8657e ‒ 123.92151/3001/6001.4279e ‒ 111/12009.1621e ‒ 133.9621Example21/1001/2007.7093e ‒ 091/4005.1758e ‒ 103.89671/2001/4005.1758e ‒ 101/8003.0120e ‒ 114.10301/3001/6009.5110e ‒ 111/12006.0554e ‒ 123.9733Table 8
Rate of convergenceρ for ε=0.1 and δ=0.03.
hh/2Zhh/4Zh/2ρExample31/1001/2004.0323e ‒ 081/4002.5200e ‒ 094.00011/2001/4002.5200e ‒ 091/8001.5799e ‒ 103.99551/3001/6004.9716e ‒ 101/12003.0522e ‒ 114.0258Example41/1001/2001.5160e ‒ 071/4009.4743e ‒ 094.00011/2001/4009.4743e ‒ 091/8005.9186e ‒ 104.00071/3001/6001.8715e ‒ 091/12001.1710e ‒ 103.9984


### 4.1. The Effect of Delay Term on the Solution Profile

To analyze the effect of the delay term on the solution profile of the problem, the numerical solution of the problem for different values of the delay parameters has been given by the following graphs.


## 4.1. The Effect of Delay Term on the Solution Profile

To analyze the effect of the delay term on the solution profile of the problem, the numerical solution of the problem for different values of the delay parameters has been given by the following graphs.


## 5. Discussion and Conclusion

Fourth-order fitted operator numerical method for solving singularly perturbed reaction-diffusion with delay has been presented. To demonstrate the efficiency of the method, four model examples without exact solutions have been considered for different values of the perturbation parameterε and delay parameter δ, and also results are presented in the tables and figures. It is observed that, from the tables, the present method approximates the solution and the stability and convergence of the method is established well. The effect of the delay on the solution of singularly perturbed delay reaction-diffusion equation is showed by plotting graphs of four model examples. Two model examples of twin layers behavior and two model examples with oscillatory layers have been considered and solved for different values of perturbation parameter ε, delay parameter δ, and mesh size h. The numerical solutions are tabulated (Tables 1–6) in terms of maximum absolute errors and observed that the present method improves the findings in [17]. Also, it is significant that all of the maximum absolute errors decrease rapidly as N increases. The results presented in Tables 7 and 8 confirmed that computational rate of convergence as well as the theoretical estimate indicates that the method is a fourth-order convergent.Furthermore, to investigate the effect of delay on the solution of the problem, numerical solutions have been presented using graphs. Accordingly, when the order of the coefficient of the delay term is ofO1, the delay affects the boundary layer solution but maintains the layer behavior (Figure 1). When the delay parameter is of Oε, the solution maintains layer behavior although the coefficient of the delay term in the equation of O1 and the delay increases, the thickness of the left boundary layer decreases while that of the right boundary layer increases (Figure 2). For the oscillatory behavior case, one can conclude that the solution oscillates throughout the domain for different values of delay parameter δ (Figures 3 and 5). In a concise manner, the present method gives more accurate solution and is uniformly convergent for solving singularly perturbed delay reaction-diffusion equations with twin layer and oscillatory behavior. Also it can see that, as mesh size h decreases, the absolute errors also decrease from Figures 5–8.Figure 5
The point-wise absolute errors of Example2 for different values of mesh size h,ε=2−8, and δ=0.5ε.Figure 6
The point-wise absolute errors of Example1 for different values of mesh size h,ε=2−8, and δ=0.5ε.Figure 7
The point-wise absolute errors of Example3 for different values of mesh size h,ε=2−8, and δ=0.5ε.Figure 8
The point-wise absolute errors of Example4 for different values of mesh size h, ε=2−8, and δ=0.5ε.



---

*Source: 5768323-2020-06-17.xml*