A finite-difference method for a singularly perturbed delay integro-differential equation

Accepted Manuscript A finite-difference method for a singularly perturbed delay integro-differential equation Mustafa Kudu, Ilhame Amirali, Gabil M. Amiraliyev PII: DOI: Reference:

S0377-0427(16)30290-4 http://dx.doi.org/10.1016/j.cam.2016.06.018 CAM 10685

To appear in:

Journal of Computational and Applied Mathematics

Received date: 2 September 2015 Revised date: 14 April 2016 Please cite this article as: M. Kudu, I. Amirali, G.M. Amiraliyev, A finite-difference method for a singularly perturbed delay integro-differential equation, Journal of Computational and Applied Mathematics (2016), http://dx.doi.org/10.1016/j.cam.2016.06.018 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Revised Manuscript Click here to view linked References

A FINITE-DIFFERENCE METHOD FOR A SINGULARLY PERTURBED DELAY INTEGRO-DIFFERENTIAL EQUATION MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA Abstract. We consider the singularly perturbed initial value problem for a linear first order Volterra integro-differential equation with delay. Our purpose is to construct and analyse a numerical method with uniform convergence in the perturbation parameter. The numerical solution of this problem is discretised using an implicit difference rules for differential part and the composite numerical quadrature rules for integral part. On a layer-adapted mesh error estimations for the approximate solution are established. Numerical examples supporting the theory are presented.

1. Introduction In this paper, we consider the following singularly perturbed delay Volterra integro-differential equation in the interval I = [0, T ] :

Lu := εu (t) + a(t)u(t) + b(t)u(t − r) + ′

= f (t), t ∈ I,

Zt

0

{K(t, s)u(s) + L(t, s)u(s − r)}ds (1.1)

u(t) = ψ(t), t ∈ I0 ,

(1.2)

m

where I = (0, T ] = ∪ Ip , Ip = {t : rp−1 < t ≤ rp }, 1 ≤ p ≤ m and rs = sr, for p=1

0 ≤ s ≤ m and I0 = [−r, 0]. ε ∈ (0, 1] is the perturbation parameter and r is a constant delay, which is independent of ε. We assume that a(t) ≥ α > 0, b(t), f (t) (t ∈ I), ψ(t) (t ∈ I0 ), K(t, s) and L(t, s) ((t, s) ∈ I × I ) are given sufficiently smooth functions satisfying certain regularity conditions to be specified. The initial value problem (1.1) has in general boundary layers on the right side of each point t = rs (0 ≤ s ≤ m − 1) for small values of ε (see Section 2). Volterra delay-integro-differential equations (VDIDE’s) arise widely in scientific fields such as biology, ecology, medicine and physics[7 − 9, 15, 21]. This class of equations plays an important role in modelling diverse problems of engineering and natural science, and hence has led researchers to develop a theory and numerical analysis for VDIDE’s. 2000 Mathematics Subject Classification. 65L11, 65L12, 65L20, 65R05. Key words and phrases. Delay Difference Scheme; Layer Adapted Mesh; Uniform Convergence; Singular Perturbation. 1

2

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

Differential equations with a small parameter ε multiplying the highest order derivative terms are said to be singularly perturbed and normally boundary layers occur in their solutions. These equations play an important role in today’s advanced scientific computations. Many mathematical models starting from fluid dynamics to the problems in mathematical biology are modelled by singularly perturbed problems. Typical examples include high Reynold’s number flow in the fluid dynamics, heat transport problem etc. For more details on singular perturbation, one can refer the books [10, 11, 22, 24, 25] and the references therein. It is well known that, for small values of ε, standard numerical methods for solving such problems are unstable and do not give accurate results. Therefore, it is important to develop suitable numerical methods for solving these problems, whose accuracy does not depend on the parameter value ε, i.e., methods that are convergent εuniformly. These include fitted finite difference methods, finite element methods using special elements such as exponential elements, and methods which use a priori refined or special non-uniform grids which condense in the boundary layers in a special manner. The various approaches to the design and analysis of appropriate numerical methods for singularly perturbed differential equations can be found in [1 − 3, 11, 16, 22, 23, 27] (see also references cited in them). One of the simplest ways to derive parameter uniform methods consists of using a class of special piecewise uniform meshes (a Shishkin mesh, see e.g.[11, 27]), which are constructed a priori and depend on the perturbation parameter, the problem data, and the number of corresponding mesh points. For a survey of early results in the theoretical analysis of singularly perturbed Volterra integro-differential equations (VIDE’s) and in the numerical analysis and implementation of various techniques for these problems we refer to the book [17]. Various approximating aspects for singularly perturbed VIDE’s have also been investigated in [2, 5, 6, 18, 20, 26, 28, 31]. Recently, there has been a growing interest in the numerical solution of VDIDE’s. For example, Koto [19] studied stability of Runge-Kutta method for VDIDE’s with a constant delay. The qualititative behaviour of numerical approximations to a nonlinear VDIDE with unbounded delay investigated in [30]. Zang and Vanderwalle [33] gave a numerical method based on the combination of the general linear methods with compound quadrature rules for VDIDE’s. Gan[12] studied the analytic and numerical dissipativity of θ-methods for nonlinear VDIDE’s. The adaptation of linear multistep methods for VDIDE’s has been discussed in [14]. Shakourifar and Enright [29] considered standard software based on the collocation method for solving VDIDE’s. The numerical stability of linear VDIDE’s with real coefficients has been discussed by Zhao et al. [35]. Bellour and Bousselsal [4] used the Taylor polynomial method for approximating VDIDE’s. Zhao et al.[34] constructed a methodology based on the sinc collocation technique to approximate pantograph VDIDE’s. The above mentioned papers, related with VDIDE’s were only concerned with the regular cases. Also singularly perturbed Volterra delay-integro-differential equations (SPVDIDE’s) frequently arise in many scientific applications also .Wu and Gan [32] investigated error behaviour of linear multistep methods applied to SPVDIDE’s and derived global error estimates A(α)− stable linear multistep methods with convergent quadrature rule. He and Xu [13] discussed the exponential stability of impulsive SPVDIDE’s. Amiraliyev and Yilmaz [3] gave an exponentially

3

fitted difference method on a uniform mesh for (1.1)-(1.2) except for a delay term in differential part and shown that the method is first-order convergent uniformly in ǫ. The aim of this paper is discretizing (1.1)-(1.2) using a numerical method, which is composed of an implicit finite difference scheme on piecewise uniform Shishkinmeshes on each time-subinterval. The scheme is constructed by the method based on using appropriate quadrature rules with the weight and remainder terms in integral form. This method of approximation has the advantage that the schemes can be effectively applied also in the case when the original problem has a solution with certain singularities (presence of boundary layer, nonsmooth solutions, etc.). But we need to have a prior information about the location and width of the layers to generate the mesh. It is proved under some additional conditions that our scheme is stable uniformly with respect to ε and that its solution approximates the solution of (1)-(2) with the error O(N −1 ln N ),where N is the mesh parameter. The approach to construct discrete problem and error analysis for approximate solution is similar to those one’s from [1 − 3]. The structure of this paper is as follows. In Section 2, we state some important properties of the exact solution. In Section 3, we describe the finite difference discretization and introduce the piecewise uniform grid. In Section 4, we present the error analysis for the approximate solution. Uniform convergence is proved in the discrete maximum norm. Numerical results are given in Section 5 to support the predicted theory. The paper ends with a summary of the main conclusions. Notation. Throughout the paper, C denotes a generic positive constant independent of ε and the mesh parameter. Some specific, fixed constants of this kind are indicated by subscripting C. For any continuous function g(t), kgk∞ denotes a continuous maximum norm on the corresponding closed interval, in particular we shall use kgk∞,I = kgk∞,p = max |g(t)| , 0 ≤ p ≤ m. We also use p

Ip∗

= {t : 0 < t ≤ rp }, 1 ≤ p ≤ m.

t∈Ip

2. Asymptotic Behavior of the Exact Solution In this section, we give a priori estimates for the solution of the problem (1.1) -(1.2), which indicate the asymptotic behavior of the solution and its first derivative in respect to perturbation parameter. These estimates are unimprovable in terms of the view of behavior in ε and will be used in order to analyse of the numerical solution. Lemma 2.1. ([3]) Consider the following initial value problem εv ′ (t) + a(t)v(t) = F (t), 0 < t < T,

(2.1)

v(0) = A.

(2.2)

Let a(t) ≥ α > 0, F (t) ∈ C(I), |F (t)| ≤ F (t) and the continuous function F (t) is nondecreasing. Then the solution of the problem (2.1) -(2.2) satisfies the following estimate

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

4

|v(t)| ≤ |A| + α−1 F (t), 0 ≤ t ≤ T.

(2.3)

Lemma 2.2. Let a, b, f ∈ C(I), ψ ∈ C(I0 ) and K, L ∈ C(I × I).Then for the solution of the problem (1.1)-(1.2) the following estimate holds kuk∞,I ∗ ≤ Cp , 1 ≤ p ≤ m,

(2.4)

p

where p  −1 Cp = kψk∞,0 α−1 kbk∞ eα (K+L)T   Z0   −1 + |ψ(0)| + α−1 kf k∞ + α−1 L |ψ(t)| dt eα (K+L)T   −r

p  p−s X −1 , α−1 kbk∞ eα (K+L)T . s=1

K = max |K(t, s)| , L = max |L(t, s)| ; p = 1, 2, ..., m. I×I

I×I

In addition if a, b, f ∈ C 1 (I), ψ ∈ C 1 (I0 ) and K, L ∈ C11 (I × I) then   (t − rp−1 )p−1 − α(t−rp−1 ) ε , (1 ≤ p ≤ m). e |u (t)| ≤ C 1 + εp ′

(2.5)

Proof. First we show the validity of (2.4). The equation(1.1) can be rewritten in the form εu′ (t) + a(t)u(t) = F (t), t ∈ Ip∗ ,

(2.6)

where

F (t) = f (t) − b(t)u(t − r) −

Zt

0

{K(t, s)u(s) + L(t, s)u(s − r)}ds.

Then

|F (t)| ≤ kf k

∗ ∞,I p

+ kbk

∗ ∞,I p

kuk

∗ ∞,I p−1

+

Zt

0

{K |u(s)| + L |u(s − r)|}ds

and based on Lemma 2.1 from (2.6) it follows that

5

o n |u(t)| ≤ |ψ(0)| + α−1 kf k∞,I ∗ + kbk∞,I ∗ kuk∞,I ∗ p

+ α−1 K

Zt

0

p

p−1

Zt |u(s)| ds + α−1 L |u(s − r)| ds. 0

o n , after replacing Denoting now δp = |ψ(0)| + α−1 kf k∞,I ∗ + kbk∞,I ∗ kuk∞,I ∗ p p p−1 s − r = η we arrive at |u(t)| ≤ δp + α

−1

Zt

|u(s)| ds + α

L

Z0

|ψ(s)| ds + α

  Zt |u(s)| ds. K +L

K

0

≤ δp + α

−1

L

−1

−r

−1

−r

Z0

|ψ(η)| dη + α

−1

L

t−r Z 0

0

Applying now the Gronwall’s inequality leads to 

|u(t)| ≤ δp + α−1 L

Z0

−r



|ψ(t)| dt ea

−1

(K+L)t

,

which implies the inequality for the solution u(t)   n kuk∞,I ∗ ≤ |ψ(0)| + α−1 kf k∞,I ∗ + kbk∞,I ∗ kuk∞,I ∗ p p p p−1  Z0  −1 +α−1 L |ψ(t)| dt ea (K+L)t .  −r

From here, the following first order difference inequality follows

with

wp ≤ qwp−1 + ρ

wp = kuk∞,I ∗ , p

−1

q = α kbk∞ ea (K+L)T ,   0 Z  −1  ρ = |ψ(0)| + α−1 kf k∞ + α−1 L |ψ(t)| dt ea (K+L)T ,   −1

−r

which yields the estimate wp ≤ w0 q + p

p X s=1

q

p−s

!

ρ,

|u(η)| dη

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

6

so we arrive at (2.4). Now we prove (2.5) by induction. For p = 1 it is known that [see,e.g.,10,11]   1 αt |u′ (t)| ≤ C 1 + e− ε . ε

Now, let the inequality (2.5) holds true for p = k |u (t)| ≤ C ′

(

(t − rk−1 )k−1 − α(t−rk−1 ) ε e 1+ εk

)

, t ∈ Ik (1 ≤ k ≤ m).

(2.7)

Differentiating the (1.1) we have the relation for p = k + 1 εu′′ (t) + a(t)u′ (t) = F (t), t ∈ Ik+1 ,

(2.8)

where

F (t) = f ′ (t) − a′ (t)u(t) − b(t)u′ (t − r) − b′ (t)u(t − r) − K(t, t)u(t) −

tZ

∂ K(t, s)u(s)ds − L(t, t)u(t − r) − ∂t

0

Zt

0

∂ L(t, s)u(s − r)ds. ∂t

Then from (2.8) we have the following relation for u′ (t)

u (t) = u (rk )e ′

′

− 1ε

Rt

rk

a(η)dη

1 + ε

Zt

F (ξ)e

Rt − 1ε a(η)dη ξ

(2.9)

dξ.

rk

Using the estimate (2.7) for t = rk , we have u (rk ) ≤ C ′

(

k−1

(rk − rk−1 ) 1+ εk

e

−

α(rk −rk−1 ) ε

)

  rk−1 αr ≤ C 1 + k e− ε . ε

Applying the inequality xk e−x ≤ Ce−γx , 0 < γ < 1, x ∈ [0, ∞) to last relation we deduce |u′ (rk )| ≤ C, k ≥ 1.

(2.10)

Furthermore, using again (2.7) we get the following estimate for F (t) |F (t)| ≤ |f ′ (t)| + |a′ (t)| |u(t)| + |b(t)| |u′ (t − r)| + |b′ (t)| |u(t − r)| + |K(t, t)| |u(t)| Zt Zt ∂ ∂ + K(t, s) |u(s)| ds + |L(t, t)| |u(t − r)| + L(t, s) |u(s − r)| ds ∂t ∂t 0

≤ C (1 + |u′ (t − r)|)   k−1 ) (t − r − rk−1 )k−1 − α(t−r−r ε e ≤C 1+ εk

0

7

  (t − rk )k−1 − α(t−rk ) ε e ≤C 1+ εk

(2.11)

Taking into account (2.10) and (2.11) in (2.9) we have

|u (t)| ≤ Ce ′

−

α(t−rk ) ε

C + ε

Zt 

1+

rk

(ξ − rk )k−1 − α(ξ−rk ) ε e εk



e−

α(t−ξ) ε

dξ

Zt

Zt C C (ξ − rk )k−1 − α(t−rk ) − α(t−ξ) ε ε ≤C+ dξ + dξ e e ε ε εk rk rk   (t − rk )k − α(t−rk ) ε , t ∈ Ik+1 , e ≤C 1+ εk+1 

which proves (2.5).

3. Discrete Problem Let ω N0 be any nonuniform mesh on I : ω N0 = {0 = t0 < t1 < ... < tN0 −1 < tN0 = T } with step-sizes τi = ti − ti−1 , which contains by N mesh points at each subinterval Ip (1 ≤ p ≤ m), ωN,p = {ti : (p − 1)N + 1 ≤ i ≤ pN }, ωN,p = {ti : (p − 1)N + 1 ≤ i ≤ pN }, 1 ≤ p ≤ m, and consequently

ωN0 =

m [

ωN,p .

p=1

Also we denote ∗ ωN,p = {ti : 1 ≤ i ≤ pN }, (1 ≤ p ≤ m).

To simplify the notation, we set gi = g(ti ) for any function g(t) and moreover yi denotes an approximation of u(t) at ti . For any mesh function gi defined on ωN0 we use the backward difference gt.i := (gi − gi−1 )/τi and norms kgk∞,N,p = max |gi | , kgk∞,ω∗ := max |gi | , 1 ≤ p ≤ m. kgk∞,ωN,p := (p−1)N ≤i≤pN

N,p

0≤i≤pN

For the difference approximation to (1.1), we integrate (1.1) over (ti−1 , ti ):

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

8

Zti

ετi−1

u (t)dt + ′

+ τi−1

ti−1

= τi−1

Zti

Zti

a(t)u(t)dt +

τi−1

Zti

ti−1

ti−1

ti−1

Zti

τi−1

b(t)u(t − r)dt

t  Z  {K(t, s)u(s) + L(t, s)u(s − r)}ds dt 0

f (t)dt.

(3.1)

ti−1

Applying right side rectangle rule yields the relation: εut,i + ai ui + bi ui−N  t Z Zti  {K(t, s)u(s) + L(t, s)u(s − r)}ds dt + Ri(1) + τi−1 ti−1

0

= fi , 1 ≤ i ≤ N0

(3.2)

with the truncation error

(1) Ri

=

−τi−1

Zti

ti−1

(t − ti−1 )

d [a(t)u(t) + b(t)u(t − r) − f (t)] dt. dt

(3.3)

Next, using right side rectangle rule in (3.2), we have t  Zti Z  {K(t, s)u(s) + L(t, s)u(s − r)}ds dt τi−1 ti−1

0

Zti (2) = τi−1 (ti − ti−1 ) {K(ti , s)u(s) + L(ti , s)u(s − r)}ds + Ri ,

(3.4)

0

where

(2) Ri

=

−τi−1

Zti

ti−1

d [(t − ti−1 ) dt

Zt

0

{K(t, s)u(s) + L(t, s)u(s − r)}ds]dt.

(3.5)

Furthermore applying the composite left side rectangle rule to the right-hand side integral term in (3.4), we obtain Zti 0

where

{K(ti , s)u(s)+L(ti , s)u(s−r)}ds =

i−1 X j=0

(3)

τj+1 [K(ti , tj )uj + L(ti , tj )uj−N ]+Ri ,

9

tj

(3) Ri

=

i Z X j=1 t

(tj − s)

d [K(ti , s)u(s) + L(ti , s)u(s − r)]ds. ds

(3.6)

j−1

Consequently we have the exact relation for u(ti ) εut,i + ai ui + bi ui−N +

i−1 X

τj+1 [K(ti , tj )uj + L(ti , tj )uj−N ] + Ri

j=0

= fi , i = 1, 2, ..., N0 with remainder term (1)

Ri = Ri

(2)

+ Ri

(3)

(3.7)

+ Ri .

As a consequence of this relation, we propose the following difference scheme for approximating (1.1) and (1.2)

LN yi := εyt,i + ai yi + bi yi−N +

i−1 X

τj+1 [K(ti , tj )yj + L(ti , tj )yj−N ]

(3.8)

j=0

= fi , i = 1, 2, ..., N0 , yi = ψi , − N ≤ i ≤ 0.

(3.9)

The difference scheme (3.8) and (3.9), in order to be ε−uniform convergent, we will use the special nonuniform mesh (Shishkin type piecewise uniform mesh). Let the discretisation parameter N be an even positive integer and let σp denote a mesh transition parameter defined by σp = rp−1 + min{r/2, α−1 θε ln N }

where θ > 1 is a constant. We consider the piecewise uniform mesh ωN,p which divides each of the interval [rp−1 , σp ] (layer region) and [σp , rp ] (outside layer) into N/2 equidistant subintervals and so, if we denote the stepsizes in these subintervals (1) (2) with τp and τp respectively, we have:

so

τp(1) = 2(σp − rp−1 )N −1 , τp(2) = 2(rp − σp )N −1 , 1 ≤ p ≤ m

ω N,p =

(

(1)

ti = rp−1 + (i − (p − 1)N )τp , i = (p − 1)N, ..., (p − 1/2)N (2) ti = σp + (i − (p − 1/2)N )τp , i = (p − 1/2)N + 1, ..., pN

where (1 ≤ p ≤ m). In the rest of the paper we only consider this type mesh.

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

10

4. Stability Bound and Convergence To investigate the convergence of this method, note that the error function zi = yi − ui , i = 0, 1, 2, ..., N0 is the solution of the discrete problem LN zi = Ri , i = 0, 1, 2, ..., N0 ,

(4.1)

zi = 0, −N ≤ i ≤ 0,

(4.2)

where the truncation error Ri is given by (3.7). Lemma 4.1. Consider the following difference problem ℓN vi := εvt,i + ai vi = Fi , i = 0, 1, 2, ..., N0,

(4.3)

v0 = A.

(4.4)

Let |Fi | ≤ Fi and Fi be nondecreasing function. Then the solution of (4.3)-(4.4) satisfies |vi | ≤ |A| + α−1 Fi , 0, 1, 2, ..., N0.

(4.5)

Proof. Consider the barrier function = ±vi + |A| + α Fi and take into consideration that Ft,i = (Fi − Fi−1 )/τ ≥ 0 because of Fi being nondecreasing function, thereby it can be seen clearly that φ± i

and

−1

−1 ℓN φ± Fi ≥ Fi + Fi ≥ 0 i ≥ ±Fi + ai α

−1 φ± Fi ≥ 0. 0 = ±A + A + α

So that accordingly to the maximum principle φ± 0 ≥ 0,



which proves (4.5).

Lemma 4.2. Let yi be the solution of (3.8) and (3.9). Then the following estimate holds kyk∞,ω∗

−1

N,p

≤ kψk∞,ωN,0 Qp + eα −1

+α−1 eα

(K+L)T

p P

s=1

(K+L)T

(|ψ(0)| + L kψk1 )

p P

Qp−s

s=1

kf k∞,ω∗ Qp−s , 1 ≤ p ≤ m N,s

(4.6)

11

where

Qp−k =

  1,  e

p=k

p Q

α−1 (K+L)T

−1

α

s=k+1

0≤k ≤p−1

kbk∞,ωN,s ,

Proof. (3.8) can be rewritten in the form εyt,i + ai yi = Fi , where i−1 X

Fi = fi − bi yi−N −

j=0

τj+1 {K(ti , tj )yj + L(ti , tj )yj−N }.

From here the following estimate holds for Fi |Fi | ≤ |fi | + |bi | |yi−N | + ≤ kf k∞,ω∗

N,p

i−1 X j=0

τj+1 {|K(ti , tj )| |yj | + |L(ti , tj )| |yj−N |}

+ kbk∞,ω∗

N,p

max |yi−N | +

0≤i≤pN

Further applying Lemma 4.1, we have |yi | ≤ |ψ(0)| + α

−1

+

 kf k∞,ω∗

N,p

i−1 X j=0

τj+1 {K |yj | + L |yj−N |}.

+ kbk∞,ω∗ max |yi−N | N,p 0≤i≤pN 

i−1 X τj+1 {K |yj | + L |yj−N |} j=0

≤ δ p + α−1 where

i X j=1

τj K |yj−1 | + α−1

 δ p = kψk∞,ωN,0 + α−1 kf k∞,ω∗

N,p

i X j=1

(4.7)

τj L |yj−N −1 | ,

+ kbk∞,ω∗

N,p

kyk∞,ω∗

N,p−1

After replacing j − N = q in (4.7) we obtain |yi | ≤ δ p + α−1 Denote that kψk1 = a) for i ≤ N

i X j=1

τj K |yj−1 | + α−1

0 P

q=1−N

i−N X

q=1−N



.

τq+N L |yq−1 | .

τq+N |ψq−1 | and from here we arrive at

12

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

|yi | ≤ δ p + α−1 L kψk1 + α−1 b) for i > N |yi | ≤ δ p + α−1 L kψk1 + α−1 which implies |yi | ≤ δ p + α−1 L kψk1 + α−1

i X

τj K |yj−1 | ,

i X

τj K |yj−1 | + α−1

i X

(K + L)τj |yj−1 | .

j=1

j=1

j=1

(4.8)

i X q=1

τq+N L |yq−1 | ,

(4.9)

From the inequalities (4.8) and(4.9) by virtue of the difference analogue of Gronwall ‘s inequality with delay, it follows that   −1  −1  |yi | ≤ δ p + α−1 L kψk1 eα (K+L)tj ≤ δ p + α−1 L kψk1 eα (K+L)T ,

thereby

kyk∞,ω∗ N,p  n ≤ |ψ(0)| + α−1 kf k∞,ω∗

N,p

(4.10) + kbk∞,ω∗

N,p

kyk∞,ω∗

N,p−1

+ L kψk1

Solving this first order difference inequality we arrive at (4.6).

o

−1

eα

(K+L)T



Lemma 4.3. Suppose that zi is the solution of the problem (4.3)-(4.4). Then the following estimate holds

kzk∞,ω∗

N,p

≤C

p X

k=1

kRk∞,ωN,k .

Proof. It evidently follows from (4.6) by taking ψ ≡ 0 and f ≡ R.

(4.11) 

Lemma 4.4. Under the above assumptions of Lemma 2.1, for the truncation error R, the following estimate holds kRk∞,N,p ≤ CN −1 ln N, 1 ≤ p ≤ m. Proof. From (3.3), taking also into consideration (2.4), on an arbitrary mesh we obtain   Zti   (1) (|u′ (t)| + |u′ (t − r)|) dt , i = 1, 2, ..., N0 , Ri ≤ C τi +   ti−1

which, in view of (2.5), yields

.

13

  Zti   αt 1 (1) e− ε dt , i = 1, 2, ..., N, Ri ≤ C τi +   ε

(4.12)

ti−1

and

 Zti  (t − rp−1 )p−1 − α(t−rp−1 ) (1) ε dt e Ri ≤ C τi +  εp

(4.13)

ti−1

+

Zti

ti−1

 (t − rp−1 )p−2 − α(t−rp−1 )  ε e dt  εp−1

where ti ∈ ωN,p (p > 1). Applying the inequality xk e−x ≤ Ce−γx , 0 < γ < 1, x ∈ [0, ∞) to (4.13) we deduce   Zti   αt 1 (1) e− θε dt , ti ∈ ωN,p , θ > 1 , i = 1, 2, ..., N. Ri ≤ C τi +   ε

(4.14)

ti−1

Combining (4.13) and (4.14) we can write

  Zti   α(t−r ) p−1 1 (1) dt , ti ∈ ωN,p , (p = 1, 2, ..., m) , θ > 1 . e− θε Ri ≤ C τi +   ε ti−1

(4.15)

For (3.5) we write

Zti Zt d (2) Ri ≤ dt {K(t, s)u(s) + L(t, s)u(s − r)}ds dt ti−1

0

and after applying Leibnitz rule

Zti (2) {|K(t, t)| |u(t)| + |L(t, t)| |u(t − r)|} dt Ri ≤ ti−1

Zti Zt ∂ dt + {K(t, s)u(s) + L(t, s)u(s − r)} ds ∂t ti−1

0

and taking into account (2.4) and (2.5), we have (2) Ri ≤ Cτi , i = 1, 2, ..., N0 .

Last for (3.6) we have

(4.16)

14

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

t   i Zj X ∂ ∂ (3) K(ti , s)u(s) + L(ti , s)u(s − r) ds (tj − s) Ri ≤ ∂s ∂s j=1 tj−1 i Ztj X (tj − s) {K(ti , s)u′ (s) + L(ti , s)u′ (s − r)} ds + j=1 tj−1 t  t i i Z Z  ≤ Cτi (|u(s)| + |u(s − r)|) ds + (|u′ (s)| + |u′ (s − r)|) ds   0

0

≤ Cτi , i = 1, 2, ..., N0 .

(4.17)

At each submesh ωN,p we estimate the truncation error Ri on Shishkin mesh as follows. We consider first the case σp = rp−1 + r/2, and so r/2 ≤ α−1 θε ln N and (1) (2) τp = τp = τp = r/N. Hereby, since 1 ε

Zti

ti−1

e−

α(t−rp−1 ) θε

dt ≤ =

ti − ti−1 τp = ε ε r ≤ 2α−1 θN −1 ln N ≡ O(N −1 ln N ) Nε

and τi = O(N −1 ) it follows from (4.15) that (1) Ri ≤ CN −1 ln N.

(4.18)

For (4.16) and (4.17) we obtain the following estimates (2) Ri ≤ Cτp ≤ CN −1

(3) Ri ≤ Cτp ≤ CN −1 .

(4.19)

(4.20)

Thus from (4.18), (4.19) and (4.20) we can establish the estimate |Ri | ≤ CN −1 ln N , (p − 1)N ≤ i ≤ pN , 1 ≤ p ≤ m.

(4.21)

We now consider the case σp = rp−1 + α−1 θε ln N and estimate Ri on [rp−1 , σp ] and [σp , rp ] separately. In the layer region [rp−1 , σp ], the inequalities (4.15), (4.16), (4.17) reduce to, respectively

15

(1) Ri ≤ C =



(

(1)

τp σp − rp−1 + 2 N ε

)

2α−1 θε ln N 2α−1 θε ln N + N N

≤ CN −1 ln N,



≤ CN −1 2α−1 θε ln N (4.22)

2 (σp − rp−1 ) (2) = CN −1 2α−1 θε ln N ≤ CN −1 ln N , Ri ≤ Cτp(1) = C N (4.23) (3) Ri ≤ CN −1 ln N

(4.24)

which obviously leads to

|Ri | ≤ CN −1 ln N , (p − 1)N ≤ i ≤ (p − 1/2)N, 1 ≤ p ≤ m.

(4.25)

Next we estimate Ri for (p − 1/2)N + 1, ..., pN . In this case, recalling that ti = σp + (i − (p − 1/2)N )τp(2) = rp−1 + α−1 θε ln N + (i − (p − 1/2)N )τp(2)

we obtain from (4.15)

   α(ti−1 −rp−1 ) α(ti −rp−1 ) (1) (2) −1 − − θε θε −e Ri ≤ C τp + α θ e   α(ti−1 −rp−1 ) −1 −1 − θε ≤ C 2(rp − σp )N + α θe

and similarly

≤ CN −1 ,

(2) Ri ≤ Cτp(2) = 2C(rp − σp )N −1 ≤ CN −1 ,

(3) Ri ≤ CN −1 .

(4.26)

(4.27) (4.28)

Hence from (4.26), (4.27) and (4.28) also on [σp , rp ] the following estimate holds |Ri | ≤ CN −1 , (p − 1/2)N ≤ i ≤ pN , 1 ≤ p ≤ m.

(4.29)

So, in view of (4.21) and (4.29) we can write on the interval Ip (3) Ri ≤ CN −1 ln N.

Thus, the proof is completed.



16

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

Theorem 4.5. Suppose that the conditions of Lemma 2.2 are satisfied. Then the solution of the problems (3.8)-(3.9) converges to the solution of (1.1)-(1.2) and satisfies the following ε−uniform error estimate ky − uk∞,ω∗N,p ≤ CN −1 ln N, 1 ≤ p ≤ m.

(4.30)

Proof. This follows immediately by combining the previous lemmas.



5. Numerical Results We now look at computational results as a particular problem

εu (t) + u(t) + u(t − 1) + ′

Zt

0

u(s − 1)ds = 1, t ∈ (0, 2]

u(t) = e−t , −1 ≤ t ≤ 0. The exact solution for 0 ≤ t ≤ 2 is given by  −t   e ε, t−1 t 1 1 u(t) = e−1 − εe− ε + εe− ε + (e−1 + εe− ε )e− ε  2(t−1) 1 1  +(e− ε + εe− ε )e− ε ,

t ∈ [0, 1] t ∈ (1, 2].

We define the exact error eN,p and the computed ε-uniform maximum pointwise ε error eN,p as follows eN,p = ky − uk∞,ωN,p , p = 1, 2, ε

eN,p = maxeN,p ε , p = 1, 2, ε

where y is the numerical approximation to u for various values of N and ε. We also define the computed parameter-uniform rate of convergence to be rN,p = ln(eN,p /e2N,p)/ ln 2, p = 1, 2. The values of ε and N for which we solve the test problem are ε = 2−i , i = 2, 4, ..., 16; N = 32, 64, 128, 256, 512, 1024.We employed the numerical method of steps. That is, on submesh ωN,2 , we used the previously computed values of y(t) on ωN,1 . The results of numerical experiment are displayed in Tables 1 and 2. From these tables we see that the rates of convergence rN,p are monotonically increasing towards one and numerical results are in good agreement with our theoretical findings.

17 N,1 Table 1. Exact errors eN,1 and ε , computed ε-uniform errors e N,1 convergence rates r on ωN,1

ε 2−2 2−4 2−6 2−8 2−10 2−12 2−14 2−16 eN pN

N = 32

N = 64

N = 128

N = 256

N = 512

N = 1024

0.099760 0.98 0.126678 0.96 0.134667 0.96 0.139987 0.95 0.147376 0.95 0.147378 0.95 0.147378 0.95 0.147378 0.95 0.147378 0.95

0.050448 0.99 0.0650013 0.97 0.068946 0.97 0.072398 0.96 0.076173 0.96 0.076162 0.96 0.076162 0.96 0.076162 0.96 0.076162 0.96

0.025399 0.99 0.033135 0.98 0.0351430 0.97 0.0371888 0.97 0.0391234 0.96 0.0391255 0.96 0.0391255 0.96 0.0391255 0.96 0.0391255 0.96

0.012787 1.00 0.016771 0.99 0.017726 0.98 0.018968 0.98 0.020094 0.97 0.0200123 0.97 0.0200123 0.97 0.0200123 0.97 0.0200123 0.97

6.350016e-003 1.00 8.430045e-003 1.00 9.073430e-003 0.99 9.608980e-003 0.98 0.0102491 0.97 0.0102491 0.97 0.0102491 0.97 0.0102491 0.97 0.0102491 0.97

3.175007e-003 4.186296e-003 4.559858e-003 4.871485e-003 5.223687e-003 5.173687e-003 5.173687e-003 5.173687e-003 5.173687e-003

N,2 Table 2. Exact errors eN,2 and ε , computed ε-uniform errors e N,2 convergence rates r on ω

ε

N = 32

N = 64

N = 128

N = 256 N = 512

N = 1024

2−2

0.124589 0.98 0.129886 0.98 0.131633 0.97 0.132799 0.97 0.134885 0.96 0.134891 0.96 0.134891 0.96 0.134891 0.96 0.134891 0.96

0.0063109 0.98 0.065611 0.98 0.067141 0.97 0.067735 0.97 0.069226 0.96 0.069237 0.96 0.069237 0.96 0.069237 0.96 0.069237 0.96

0.0031943 0.98 0.033214 0.98 0.0334246 0.98 0.033457 0.97 0.0355608 0.96 0.0355746 0.96 0.0355746 0.96 0.0355746 0.96 0.0355746 0.96

0.016046 0.99 0.016824 0.99 0.017331 0.98 0.017610 0.98 0.018277 0.97 0.018392 0.97 0.018392 0.97 0.018392 0.97 0.018392 0.97

4.060531e-003

2−4 2−6 2−8 2−10 2−12 2−14 2−16

eN pN

8.071963e-003 0.99 8.463149e-003 0.99 8.785166e-003 0.98 8.894040e-003 0.98 9.329126e-003 0.97 9.343173e-003 0.97 9.343173e-003 0.97 9.343173e-003 0.97 9.343173e-003 0.97

4.257298e-003 4.445991e-003 4.505179e-003 4.7584653e-003 4.7645349e-003 4.7645349e-003 4.7645349e-003 4.7645349e-003

18

MUSTAFA KUDUA,∗ , ILHAME AMIRALIB AND GABIL M. AMIRALIYEVA

6. Conclusion A singularly perturbed delay integro-differential equation for a first-order initial value problem is considered. The difference scheme is constructed by the method of integral identities with the use of interpolating quadrature rules with the weight and remainder terms in integral form. The numerical method presented here comprises a standard backward difference operators on a non-uniform mesh which consists of the special piecewise uniform meshes on each subinterval. It is shown that the method displays uniform convergence with respect to the perturbation parameter. Numerical results confirm our theoretical analysis. The main lines for the analysis of the uniform convergence carried out here can be used for the study of more complicated linear as well as nonlinear singularly perturbed Volterra delay-integrodifferential equations. References [1] Amiraliyev, G.M., Erdo˘ gan, F.: Uniform Numerical Method for Singularly Perturbed Delay Differential Equations, Comput. Math. Appl. 53,1251-1259(2007). [2] Amiraliyev, G.M., S ¸ evgin,S.: Uniform Difference Method for Singularly Perturbed Volterra Integro- Differential Equations, Appl. Math. Comput.179, 731-741(2006). [3] Amiraliyev, G.M. Yilmaz, B.: Fitted Difference Method for a Singularly Perturbed Initial Value Problem, J. Comput. Appl. Math. 22, 1-10(2014). [4] Belloura, A. and Bousselsal, M.: Numerical solution of delay integro-differential equations by using Taylor collocation method, Math. Meth. Appl. Sci. 37, 1491–1506(2014). [5] Bijura, A.M.: Singularly Perturbed Volterra Integro-differential Equations, Quaestiones Mathematicae, 25:2, 229-248(2002), DOI: 10.2989/16073600209486011 [6] Bobodzhanov, A. A. and Safonov, V. F.: Singularly Perturbed Integro-Di erential Equations with Diagonal Degeneration of the Kernel in Reverse Time, Differential Equations, 40(1),120– 1272004. [7] Bocharov, G.A., Rihan, F.A., Numerical modelling in biosciences with delay differential equations, J. Comput. Appl. Math. 125, 183–199(2000) [8] Brunner, H., van der Houwen, P.J. The Numerical Solution of Volterra Equations, in: CWI Monographs, Vol. 3, North-Holland, Amsterdam, 1986. [9] De Gaetano, A., Arino, O.: Mathematical modelling of the intravenous glucose tolerance test. J.Math.Biol. 40, 136–168 (2000) [10] Doolan, E.R., Miller, J.J.H., Schilders, W.H.A.: Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [11] Farrel, P.A., Hegarty, A.F., Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Robust Computational Techniques for Boundary Layers, Chapman Hall/CRC, New York, 2000. [12] Gan, S.: Dissipativity of θ-Methods for Nonlinear Volterra Delay Integro-Differential Equations, Comput. Math. Appl. 206, 898-907(2007). [13] He, D. and Xu,L.:Integrodifferential Inequality for Stability of Singularly Perturbed Impulsive Delay Integrodifferential Equations, Journal of Inequalities and Applications, 2009, ID 369185, 11, doi:10.1155/2009/369185. [14] Huang, C.: Stability of linear multistep methods for delay integro-differential equations, Computers and Mathematics with Applications 55, 2830–2838(2008). [15] Jerri, A. Introduction to Integral Equations with Applications, Wiley, New York, 1999. [16] Kadalbajoo, M.K., Gupta, V.: A brief survey on numerical methods for solving singularly perturbed problems, Appl. Math. Comput. 217, 3641-3716(2010). [17] Kauthen, J.P.: A survey on singularly perturbed Volterra equations, Appl. Numer. Math. 24, 95-114(1997). [18] Khater, A. H., Shamardan A. B., Callebaut D. K., Sakran, M. R. A.: Numerical solutions of integral and integro-differential equations using Legendre polynomials, Numer Algor 46,195– 218(2007), DOI 10.1007/s11075-007-9130-2 [19] Koto, T.: Stability of Runge–Kutta methods for delay integro-differential equations, J. Comput. Appl. Math.145 ,483–492(2002).

19

[20] Lodge, A.S., McLeod, J.B.: Nohel, J.A.: A nonlinear singularly perturbed Volterra integrodifferential equation occurring in polymer rheology, Proc. Roy. Soc. Edinburgh Sect. A 80,99-137(1978). [21] Marino, S., Beretta, E., Kirschner, D.E.: The role of delays in innate and adaptive immunity to intracellular bacterial infection. Math. Biosci. Eng. 4, 261–288(2007). [22] Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [23] Mishra, H. K., Saini S.: Various Numerical Methods for Singularly Perturbed Boundary Value problems, American J. of Appl. Math. and Statistics, 2, 129-142(2014). [24] Nayfeh, A.H.: Introduction to Perturbation Techniques, Wiley, New York, 1993. [25] O’Malley, R.E.: Singular Perturbation Methods for Ordinary Differential Equations, Springer Verlag, New York, 1991. [26] Ramos, J.I.:Exponential techniques and implicit Runge Kutta method for singularly perturbed Volterra integro differential equations, Neural Parallel Sci. Comput.16, 387-404(2008). [27] Roos, H.G., Stynes, M., Tobiska, L.: Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [28] Salama, A.A., Bakr, S.A.: Difference schemes of exponential type for singularly perturbed Volterra integro-differential problems, Appl. Math. Model. 31, 866-879(2007). [29] Shakourifar, M., Enright, W.: Superconvergent interpolants for collocation methods applied to Volterra integro-differential equations with delay, BIT Numer Math 52, 725–740(2012),DOI 10.1007/s10543-012-0373-5. [30] Song, Y., Baker, C. T.H.: Qualitative behaviour of numerical approximations to Volterra integro-differential equations, J. Comput. Appl. Math. 172 , 101–115 (2004). [31] S ¸ evgin,S.:Numerical solution of a singularly perturbed Volterra integro-differential equation, Adv. Difference Equ. 2014,171(2014). [32] Wu, S., Gan, S.: Errors of linear multistep methods for singularly perturbed Volterra delayintegro-differential equations, Math. Comput. Simulat. 79, 3148-3159(2009). [33] Zhang, C., Vandewalle, S.: Stability Analysis of Volterra Delay- Integro- Differential Equations and their Backward Differentiation Time Discretization, J. Comput. Appl. Math. 164165, 797-814 (2004). [34] Zhao, J., Cao, Y. and Xu, Y.: Sinc numerical solution for pantograph Volterra delayintegro-differential equation, International Journal of Computer Mathematics, (2016), DOI: 10.1080/00207160.2016.1149577. [35] Zhao, J., Fan, Y., Xu, Y.: Delay-dependent stability analysis of symmetric boundary value methods for linear delay integro-differential equations, Numer Algor 65, 125–151(2014), DOI 10.1007/s11075-013-9698-7 a Department of Mathematics, Faculty of Arts and Sciences, Erzincan University, b Department of Mathematics, Faculty of Arts and Sci24100, Erzincan, Turkey, ¨ ¨ zce, Turkey ences, Duzce University, 816200 Du ∗ [email protected], [email protected], [email protected] E-mail address: