Iran J Sci Technol Trans Sci DOI 10.1007/s40995-017-0268-z

RESEARCH PAPER

The Numerical Solution of Volterra Integro-Differential Equations with State-Dependent Delay M. Zarebnia1 • L. Shiri1

Received: 21 April 2015 / Accepted: 4 September 2015 Ó Shiraz University 2017

Abstract This paper presents the piecewise collocation method for a class of functional integro-differential equations with state-dependent delays. A numerical analysis is applied to solve the delay integro-differential equation and detect the points that have discontinuity in the derivatives of the solution of the problem and insert them into the mesh to guarantee the required accuracy. The convergence analysis of the method is investigated and a numerical experiment is presented for clarifying the robustness of the method. Keywords Delay integro-differential equation  State-dependent delay  Collocation method  Convergence analysis

1 Introduction

(2011). Investigations of equations with state-dependent delay essentially differ from those of equations with constant or time-dependent delay. For these reasons the theory of equations with state-dependent delay has drawn the attention of researchers in recent years, see for instance Hartung (2001) and Hartung et al. (2000). The main objective of the current study is to implement the collocation method for delay Volterra integro-differential equation (DVIDE) of the form Z t 0 y ðtÞ ¼ f ðtÞ þ k1 ðt; sÞyðsÞds þ

Z

0 hðt;yðtÞÞ

k2 ðt; sÞyðsÞds; t 2 ½a; T;

ð1Þ

0

yðtÞ ¼ /ðtÞ; t 2 ½h ; a; where h ¼ min hðt; yðtÞÞ:

Delay Volterra integro-differential equations (DVIDEs) arise widely in scientific fields such as biology (Marino et al. 2007), medicine (De Gaetano and Arino 2000), control theory (Hale and Verduyn Lunel 1993; Kuang 1993; and their references therein), etc. Due to the practical application of these equations, they must be solved successfully with efficient numerical approaches. In recent years, the numerical solutions of integral equations with constant or time-dependent delay have been investigated by many authors, for example, Bellour and Bousselsal (2014), Vanani et al. (2012) and Shakourifar and Enright

t2½a;T

Let hðt; yðtÞÞ  t be defined on ½a; T. In particular, our attention is directed towards state-dependent problems where the retarding function h is allowed to vary not only as a function of t but also of yðtÞ. Our motivation comes from the fact that Eq. (1) may be viewed as a generalization of the equation encountered in the theory of a circulating fuel nuclear reactor (Hale and Verduyn Lunel 1993) Z t y0 ðtÞ ¼ kðt; sÞyðsÞds: ts

& M. Zarebnia [email protected] 1

Department of Mathematics, University of Mohaghegh Ardabili, Ardabil, Iran

A serious difficulty is the possible loss of regularity of the solution even in the presence of smooth functions f ðtÞ, k1 ðt; sÞ, k2 ðt; sÞ and hðt; yÞ in the problem (1). If y0 ðaÞ 6¼ /0 ðaÞ, the derivative solution is obviously

123

Iran J Sci Technol Trans Sci

discontinuous at a. However, even when y0 ðaÞ ¼ /0 ðaÞ there is no reason why the (right-hand) derivative y00 ðaÞ which is given by the DVIDE should be equal to the (lefthand) derivative /00 ðaÞ. This irregularity at a may propagate along the integration interval by means of the deviating argument hðt; yðtÞÞ (Feldstein and Neves 1984). Evidently, as soon as hðf; yðfÞÞ ¼ a for some f [ a, due to the fact that some derivatives of yðtÞ have a jump discontinuity at a, the function is not smooth at f. In general, this creates a further jump discontinuity in some derivatives of the solution yðtÞ at f, and so on. Such discontinuity points are referred to in the literature as breaking points. In the case h depends on t (but not on the solution yðtÞ), the breaking points can be computed in advance by solving first the equation hðf1 Þ ¼ a for f1 , then for every solution fk the equation hðfkþ1 Þ ¼ fk , etc. (Brunner 2004). This may be useful, if these equations do not have too many solutions. For an efficient computation, computed breaking points can be inserted in advance into the mesh of collocation. But in the general (so-called state-dependent) case such an a priori computation is not possible. The main purpose of this paper is to use an algorithm that automatically computes the breaking points and includes them in the collocation points. In this way the accuracy of the approximation will be improved. The discussion in this paper is organized as follows. In the next section we will review some preliminary concepts and results, which we will use throughout the paper. We explain how collocation method can be extended to problems of the form (1), and we briefly discuss the appearance of breaking points. An algorithm for solving DVIDE of the form (1) and detecting the breaking points is presented in Sect. 3. A few theoretical results (Sect. 4) and the numerical experiment (Sect. 5) conclude the discussion.

2 Supporting Results In this section, we give some basic important results, which we will use throughout the paper. At first, we consider the following definition which is followed from Neves and Feldstein (1976). Definition 2.1 The DVIDE (1) has continuity class p  1, if and only if the following hold over the appropriate domains: (i) (ii)

All of the mixed partial derivatives of k1i;j and k2i;j are continuous, for all i þ j  p: All of the mixed partial derivatives of hi;j are continuous, for all i þ j  p:

123

(iii)

f 2 C p ½a; T and / 2 C p ½h ; a, /ðaÞ ¼ yðaÞ, however /ðjÞ ðaÞ need not necessarily equal yðjÞ ðaÞ for any j  1.

In what follows, we give an explanation of several principles concerning breaking point propagation without a lengthy and detailed hypothesis. Definition 2.2 Suppose that fk be the zero of the equation hðt; yðtÞÞ ¼ fk1 , then fk1 is known as an ancestor of fk and f0 ¼ a is called a 0-level breaking point. Consequently, fk is called a k-level breaking point, and fk1 be a ðk  1Þ-level breaking point. Definition 2.3 (Bellen and Zennaro 2003) A breaking point fk is said to be of order zk if yðsÞ ðfk Þ exists for s ¼ 0; . . .; zk : Theorem 2.1 (Hauber 1997) If the problem (1) has continuity class p  1, then the following statements hold: (i)

(ii)

(iii) (iv)

(v)

Let S0 :¼ fag and for 1  k  p define Sk :¼ ff^ 2 ½a; T : f^ is a zero with odd (left handed) integer multiplicity m  p of hðt; yðtÞÞ  f, where f\f^ and S f 2 Sk1 g: Then pk¼0 Sk contains all possible breaking points in the first p þ 1 derivatives of y in ½a; T: If f 2 Sk , then y 2 C k at f for 1  k  p. If y does not belong in Cz at f 2 ½a; T for some S z  p þ 1 then f 2 z1 k¼0 Sk : The breaking points of level z  p þ 1 are not fixed points of h. Let hypothesis (H) hold for breaking points f, f^ and let m be the (left-handed) integer multiplicity of the zero f^ of hðt; yðtÞÞ  f, then m is odd and z^ mz þ 1: y 2 C pþ1 between the breaking points of level z  p þ 1.

Here we have used the hypothesis (H) of Neves and Feldstein, (1976) which guarantees that the breaking points are isolated: Hypothesis (H) Let f^ 2 ða; TÞ be a breaking point with ^ yðfÞÞ ^ 2 ½a; TÞ: Let there exist e , e^ with ancestor f ¼ hðf; f þ e  f^  e^ and z , z^ 2 f1; 2; . . .; p þ 1g such that   ^ e^Þ;  ^ ^ y 2 Dpþ1 ðf; yðfe;fþeÞ 2 Dpþ1 z ðf; eÞ and y ðf^ z^ e;fþ^ eÞ where  z1 Dpþ1 nCz ðf  e; f þ eÞ : yðfe;f z ðf; eÞ :¼ fy 2C 2 C p ; y½f;fþeÞ 2 C p g: If fk is a zero of the equation whose multiplicity is even, then fk cannot be a breaking point, that is, y 2 C pþ1 near fk . Noting that such a fk has no computational role (For

Iran J Sci Technol Trans Sci

more details see Bellen and Zennaro 2003; Neves and Feldstein 1976). So from now on, our mean with regard to the breaking points is the odd multiplicity zeros, with order of continuity less than or equal to p þ 1 and in the case z0  p þ 1, it is not necessary to compute them.

Li ðvÞ ¼

p a v  ck ; v 2 ½0; 1; c  ck k¼1;k6¼i i

and set u0 ðtn þ vhn Þ ¼

p X

Li ðvÞYn;i ; v 2 ð0; 1;

i¼1

3 Description of the Method

where

We approximate Eq. (1) with the piecewise collocation method in Sect. 3.1. By use of an appropriate strategy for detecting the breaking points which has been taken from Guglielmi and Hairer (2008), we compute and insert them into the set of mesh points in Sect. 3.2.

Let Ih :¼ ftn : a ¼ t0 \t1 \    \tN ¼ Tg be a non uniform mesh on the given interval I :¼ ½a; T and hn :¼ tnþ1  tn ; 0  n  N  1; where the diameter of the mesh is h ¼ maxn hn . For a given integer p  1, the linear space of piecewise polynomials with respect to the mesh Ih is defined by   Sð0Þ p ðIh Þ :¼ fu 2 CðIÞ : u ½tn ;tnþ1  2 pp ð0  n  N  1Þg; where pp denotes the space of all polynomials of degree not

þ

ð0Þ Sp ðIh Þ

0 hðt;yðtÞÞ

Z

So, the restriction of the collocation solution to the interval ½tn ; tnþ1 , i.e. un ðtÞ, can be expressed as un ðtÞ ¼ uðtn þ vhn Þ ¼ yn þ hn

p X

bj ðvÞYn;j ; v 2 ½0; 1;

ð3Þ

i¼1

3.1 Collocation Method

exceeding p. The collocation solution u 2 defined by the collocation equation Z t 0 k1 ðt; sÞuðsÞds u ðtÞ ¼ f ðtÞ þ

Yn;i :¼ u0 ðtn þ ci hn Þ:

k2 ðt; sÞuðsÞds; t 2 Xh ;

Rv where yn ¼ uðtn Þ and bj ðvÞ ¼ 0 Lj ðsÞds: Since hðtn;i ; uðtn;i ÞÞ  tn;i , there exists an index j  n , tj  hðtn;i ; uðtn;i ÞÞ  tjþ1 ;

ð4Þ

or hðtn;i ; uðtn;i ÞÞ  a, in this case we have Z hðt;uðtÞÞ Z hðt;uðtÞÞ k2 ðt; sÞuðsÞds ¼ k2 ðt; sÞ/ðsÞds: 0

0

We obtain the approximate equations, Yn;i ¼ f ðtn;i Þ þ hn

Z

(

ci

k1 ðtn;i ; tn þ shn Þ yn þ hn

0

can be

p X

) bj ðsÞYn;j ds

j¼1

þ Fn ðtn;i Þ þ ðVh uÞðtn;i Þði ¼ 1; . . .; pÞ;

ð5Þ where ð2Þ

Fn ðtÞ ¼

Z

tn

k1 ðt; sÞuðsÞds ¼

0

0

uðtÞ ¼ /ðtÞ; t 2 ½h ; aÞ;

n1 X

¼

where hðt; uðtÞÞ  t, for all t 2 Xh and Xh :¼ ftn;i ¼ tn þ ci hn : 0  c1 \    \cp  1; ð0  n  N  1Þg; is the set of collocation points including the collocation parameters fci g: The collocation (2) has the breaking points similar to the original equation, namely a ¼ fh;0 \fh;1 \    \fh;k \    \fh;g  T whose levels are less than p þ 1 and can be obtained from the following equation: hðt; uðtÞÞ ¼ fh;k ; 0  k  g;

n1 X

Z

l¼0

Z

(

1

1

k1 ðt; tl þ shl Þuðtl þ shl Þds

hl 0

k1 ðt; tl þ shl Þ yl þ hl

hl 0

l¼0

p X

) bj ðsÞYl;j ds;

j¼1

and ðVh uÞðtÞ ¼ Wðj1Þ ðtÞ þ n

Z

hðt;uðtÞÞ

k2 ðt; sÞuðsÞds; tk

to approximate Wðj1Þ ðtÞ we proceeded as follows: n wðj1Þ ðtÞ ¼ n

Z

tj

k2 ðt; sÞuðsÞds ¼

0

Z 1 j1 X hl k2 ðt; tl þ shl Þuðtl þ shl Þds l¼0

0

Z 1 p j1 X X hl k2 ðt; tl þ shl Þfyl þ hl bj ðsÞYl;j gds: ¼ l¼0

0

j¼1

fh;0 ¼ a: A convenient computational form of the collocation equation is provided when we employ the local Lagrange basis functions

Recall that tj  hðtn;i ; uðtn;i ÞÞ  tjþ1 , so we can write hðtn;i ; uðtn;i ÞÞ ¼ tj þ c~i hj , then

123

Iran J Sci Technol Trans Sci

hðtn;i ; uðtn;i ÞÞ  tj hj P hðtn;i ; yn þ hn pj¼1 bj ðci ÞYn;j Þ  tj ¼ ; hj

c~i ðYn Þ ¼

Also, Z

hðtn;i ;uðtn;i ÞÞ

k2 ðtn;i ; sÞuðsÞds ¼

tj

¼ hj

Z

tj þhj c~i ðYn Þ

k2 ðtn;i ; sÞuðsÞds

tj

Z

c~i ðYn Þ

k2 ðtn;i ; tj þ shj Þuðtj þ shj Þds ( ) Z c~i ðYn Þ p X ¼ hj k2 ðtn;i ; tj þ shj Þ yj þ hj bj ðsÞYj;j ds: 0

0

j¼1

The solution Yn 2 Rp is given by the following nonlinear system Yn ¼ fn þ hn yn Kn þ h2n Cn Yn þ Fn þ Qðj1Þ þ hj yj K~n ðYn Þ n 2~ þ hj Cn ðYn ÞYj ; ð6Þ and the matrices are Z ci ðCn Þi;j :¼ k1 ðtn;i ; tn þ shn Þbj ðsÞds 0

ðC~n ðYn ÞÞi;j :¼

Z

c~i ðYn Þ

k2 ðtn;i ; tn þ shn Þbj ðsÞds; 0

ðj1Þ the vectors Kn , K~n ðYn Þ, fn , F and Qn are defined by Z ci k1 ðtn;i ; tn þ shn Þds; ðKn Þi :¼ n 0

ðK~n ðYn ÞÞi :¼

Z

c~i ðYn Þ

k2 ðtn;i ; tj þ shj Þds;

computed in advance and inserted into the collocation points. In the DVIDEs under consideration, the breaking points cannot be computed in advance, therefore we should obtain the collocation solutions and the breaking points simultaneously. In our numerical algorithm, all breaking points with level less than or equal to p þ 1 are computed to get a method with order of convergence p. To compute Sk , a set of numerical approximation to Sk in Theorem 2.1, we start with S0 ¼ fag, and t1 :¼ a þ h, where h is an initial step size. Suppose that the equation is solved successfully until tn . The problem is to find the zeros of the function dk1 ðtÞ ¼ hðt; un1 ðtÞÞ  fh;k1 ; where fh;k1 2 Sk1 is a previous breaking point and un1 ðtÞ is the collocation solution of the preceding step. An approach would be to check in every interval whether at least one of the functions dk1 ðtÞ changes sign. The breaking point can then be localized and added to the set Sk : Once dk1 ðtÞ changes sign on ½tn ; tn þ h, it is expected that there exists a breaking point in this interval. We compute this breaking point and call it fh;k . After that, we set tnþ1 :¼ fh;k and hn :¼ fh;k  tn ; and then we consider the system (6) for the unknowns Yn;1 ; . . .; Yn;p : For hn the system (6) is usually solved by simplified Newton iterations, while it is necessary to choose suitable initial values for unknowns with desired accuracy (we take ½0

Yn;j :¼ Yn1;j ). Noting that in the initial point a , we use u00 ðaþ Þ as an approximation of Y0;i such that Z hða;/ðaÞÞ u00 ðaþ Þ ¼ f ðaÞ þ k2 ðt; sÞ/ðsÞds: 0

0

fn :¼ ðf ðtn;1 Þ; . . .; f ðtn;p ÞÞT ; Fn :¼ ðFn ðtn;1 Þ; . . .; Fn ðtn;p ÞÞT ; :¼ ðwðj1Þ ðtn;1 Þ; . . .; wðj1Þ ðtn;p ÞÞT : Qðj1Þ n n n This nonlinear system can be solved numerically to determine Yn : To approximate the j in (4), we propose a method based on the computation of un1 ðtn1;p Þ. Assume that un1 ðtn1;p Þ is computed in the previous interval ½tn1 ; tn . Using un1 ðtn1;p Þ, we approximate j tj  hðtn;i ; un1 ðtn1;p ÞÞ  tjþ1 : 3.2 Numerical Approximation of Breaking Points and Computation of the Solution

We summarize the above discussions in the following algorithm. Algorithm 3.1 Suppose that the problem is computed successfully until n: If dk1 ðtn Þdk1 ðtn þ hÞ\0look for zeros of the function dk1 ðtÞ for t 2 ½tn ; tn þ h: Compute the root fh;k ; then replace h with hn :¼ fh;k  tn , and consider tnþ1 :¼ fh;k (the next grid point). Then compute (6) with hn : Otherwise compute Yn;1 ; . . .; Yn;p in relation (6) by applying a simplified Newton iteration with respect to the variables Yn;1 ; . . .; Yn;p and with given h; this yields an approximation un ðtn þ vhn Þ for v 2 ½0; 1.

4 Convergence The efficiency of a numerical method for DVIDEs can be improved significantly if breaking points are computed carefully. For constant delays and for delays hðtÞ not depending on the solution, all breaking points can be

123

The convergence of the collocation method which is introduced in the previous sections is discussed in the present section. The following Lemma from Feldstein and

Iran J Sci Technol Trans Sci

Neves (1984) is useful to obtain the collocation error at the near breaking points. Lemma 4.1 If (1) has continuity class p  1 , then for all sufficiently small h [ 0, ky^  yk½f1i ;fi þh ¼ Oðhzi Þ;  ½fi h;fi  ^  ¼ Oðhzi Þ; y^  y

Now, let fh;i denote an approximate value of fi , which is generated by the introduced approach and ri is the rate of   error, i.e. fh;i  fi  ¼ Oðhri Þ. We can now prove the following theorem.

j¼1

Z

v

0

ð1Þ

Rpþ1;l ðsÞds:

Recalling the local representation (3) of the collocation solution u on ½tl ; tlþ1 , and setting el;j ¼ Zl:j  Yl;j , the collocation error e :¼ y  u may be written as eðtl þ vhl Þ ¼ eðtl Þ þ hl

bj ðvÞel;j þ hpþ1 Rpþ1;l ðvÞ; v l

j¼1

2 ½0; 1;

e0 ðtl þ vhl Þ ¼

jhðt; yÞ  hðt; zÞj  Lh jy  zj: The DVIDE (1) has continuity class p  1 .

(iii)

u 2 Sp ðIh Þ is the collocation solution to (1) defined by (3) and (6). The DVIDE (1) has a finite number of breaking points in the first p þ 1 derivatives of yðtÞ, each of which satisfies hypothesis (H).

ð0Þ

Then the estimates  ðvÞ  y  uðvÞ  ¼ Oðhminfrz;pg Þðv ¼ 0; 1Þ 1 hold for any set Xh 0  c1 \    \cp  1.

of collocation points with

Proof Let us assume that there is not any breaking point on ½tl ; tlþ1  , we conclude y 2 Cpþ1 ½tl ; tlþ1 , and hence y0 2 C pþ1 ½tl ; tlþ1 . Thus we have, using Peano’s Theorem for y0 on ½tl ; tlþ1 , p X

Lj ðvÞZl;j þ

ð1Þ hpl Rpþ1;l ðvÞ; v

2 ð0; 1;

p X

ð1Þ

Lj ðvÞel;j þ hpþ1 Rpþ1;l ðvÞ; v 2 ½0; 1: l

ð8Þ

j¼1

(ii)

ð7Þ

Now we suppose that there is a breaking point in ½tl ; tlþ1 ; since the set of collocation points includes all of the breaking points, we can write tl ¼ fh;k : Let us assume that fh;k ¼ fk and y^ is a C pþ1 extension of y^ from ½fk ; tlþ1  back to ½fh;k ; tlþ1 ; consequently using Peano’s Theorem, we have jeðtl þ vhl Þj  jyðtl þ vhl Þ  y^ðtl þ vhl Þj þ jy^ðtl þ vhl Þ  uðtl þ vhl Þj p  X   ½fh;k ;fk  bj ðvÞel;j   ky  y^k1 þjeðtl Þj þ hl þ

with Zl;j ¼ y0 ðtl;j Þ. The Peano reminder term and Peano kernel are given by Z 1 ð1Þ Kp ðv; zÞyðpþ1Þ ðtl þ zhl Þdz; Rpþ1;l ðvÞ :¼ 0





j¼1

Rpþ1;l ðvÞ: hpþ1 l

ð9Þ

Since eh is continuous in ½tl ; tlþ1 ; we also have the relation with bj ¼ bj ð1Þ; eðtl Þ ¼ eðtl1 þ hl1 Þ ¼ eðtl1 Þ þ hl1

p X

bj el1;j þ hpþ1 l1 Rpþ1;l1 ð1Þ;

j¼1

j¼1

for l ¼ 1 we have eðt1 Þ ¼ eðh0 Þ ¼ h0

p X

bj e0;j þ hpþ1 0 Rpþ1;0 ð1Þ:

j¼1

According to (9), we obtain

and (

Kp ðv; zÞ :¼

p X

While

h satisfies the Lipschitz condition

y ðtl þ vhl Þ ¼

bj ðvÞZl;j þ hpþ1 Rpþ1;l ðvÞ; v l

Assume

Theorem 4.1

0

p X

2 ½0; 1;

Rpþ1;l ðvÞ :¼

^ are where zi is the order of continuity of fi and y^ and y^ extensions of y .

(iv)

yðtl þ vhl Þ ¼ yðtl Þ þ hl

where

1

(i)

Integration of (7) leads to

1 ðv  zÞp1 þ  ðp  1Þ!

p X k¼1

) Lk ðvÞðck  zÞp1 : þ

½fh;k ;fk  þjeðtl1 Þj jeðtl Þj ¼ jeðtl1 þ hl1 Þj  ky  y^k1 p   X    bj el1;j  þ hpþ1 Rpþ1;l1 ð1Þ: þ hl1 l1 j¼1

123

Iran J Sci Technol Trans Sci

We suppose that there are k breaking points until tl and r :¼ maxi ri and z :¼ maxi zi , so we find that jeðtl Þj 

k X

ky  y^k½f1i;k ;fi  þ

i¼0

l1 X

(

p    X   bj ei;j  þ hpþ1 Rpþ1;i ð1Þ hi

i¼0

 Oðhrz Þ þ b

l1 X

)

j¼1

hi kei k1 þ Oðhp Þ;

i¼0

½fh;k ;fk  þjeðtl Þj jeðtl þ vhl Þj  ky  y^k1 p X     bj ðvÞel;j  þ hpþ1 Rpþ1;l ðvÞ þ hl l j¼1

 Oðhrz Þ þ b

l1 X

Z

1

jeðtk þ shk Þjds þ hl K1

Z

i¼0

l X

tl;i

0

k¼0

 jeðtl þ shl Þjds 0   þ K2 Lh k yk1 eðtl;i Þ ( Z 1 l1 k X X    ðK1 þ K2 Þ hk ½Oðhrz Þ þ Oðhp Þþb hj ej 1 ds þ hl

Z

ci

Z

ci

0

k¼0

"

j¼0

Oðhrz Þ þ Oðhp Þ þ b

j¼0

"

þ K2 Lh k yk1

hi ke i k1 ;

ð10Þ  R1  and kp :¼ maxv2½0;1 0 Kp ðv; zÞ

where Mpþ1 1 dz: In the rest, we find an upper bound for kei k1 , 0

Oðh Þ þ Oðh Þ þ b rz

p

0

 ¼

Z 0

þ

0

Z

tl;i

k1 ðtl;i ; sÞuðsÞds þ

0 tl;i

k1 ðtl;i ; sÞeðsÞds þ Z

Z

k2 ðtl;i ; sÞyðsÞds Z

l1 X

# k X     hj e j 1 ;

ðt; sÞ 2 ½a; T  ½0; T g

0 hðtl;i ;uðtl;i ÞÞ

k2 ðtl;i ; sÞuðsÞds

for

b

l1 X

j¼0

)

hi kei k1 ds

þ Oðhrz Þ þ Oðhp Þ

j¼0

( ) l1 l1 X X   bðK1 þ K2 Þ T hk ke k k1 þ hl hi ke i k1 j¼0

k¼0

k2 ðtl;i ; sÞuðsÞds

þ Oðhrz Þ þ Oðhp Þ  1 þ K2 Þ ¼ bðK

l1 X

ðT þ hl Þhk kek k1 þ Oðhrz Þ þ Oðhp Þ:

k¼0

þ hl

Z

k1 ðtl;i ; tk þ shk Þeðtk þ shk Þds

kel k1  C1

l1 X

ðT þ hl Þhk kek k1 þ Oðhrz Þ þ Oðhp Þ;

k¼0 ci

k1 ðtl;i ; tl 0 hðtl;i ;uðtl;i ÞÞ

þ shl Þeðtl þ shl Þds

k2 ðtl;i ; sÞeðsÞds Z

Therefore,

1

0

þ

ci

hðtl;i ;uðtl;i ÞÞ

k2 ðtl;i ; sÞyðsÞds Z

0

ds

hðtl;i ;yðtl;i ÞÞ

hk

Z

0

k¼0

Z

0

k¼0

þ

  hj ej 1

#)

 1 þ K2 Lh kyk ÞÞkel k ð1  hbðK 1 1 ( Z 1 X l1 k X    ðK1 þ K2 Þ hk hj ej 1 ds b

0

hðtl;i ;yðtl;i ÞÞ

hðtl;i ;uðtl;i ÞÞ

¼

where Ki ¼ supfjki ðt; sÞj; i ¼ 1; 2. Then, we get

þ hl

0

el;i ¼ y ðtl;i Þ  u ðtl;i Þ Z Z tl;i k1 ðtl;i ; sÞyðsÞds þ ¼

k X

j¼0

i¼0

  :¼ yðpþ1Þ 

jeðtl þ shl Þjds

  jeðsÞjds þ K2 Lh k yk1 eðtl;i Þ 0 ( Z 1 l1 X hk ¼ ðK1 þ K2 Þ jeðtk þ shk Þjds þ K2

þ hl bkel k1 þhpl Mpþ1 kp T  Oðhrz Þ þ Oðhp Þ þ b

ci

0

0

hi kei k1 þ Oðhp Þ

Z

0

k¼0

þhl

Pp   j¼1 bj . Thus

where b ¼

l1 X   el;i   K1 hk

hðtl;i ;yðtl;i ÞÞ

hðtl;i ;uðtl;i ÞÞ

Therefore,

123

k2 ðtl;i ; sÞyðsÞds:

this inequality is a generalised discrete Gronwall inequality; its solution is bounded by kel k1  ðOðhrz Þ p

þ Oðh ÞÞ exp C1

l1 X k¼0

þ Oðhp ÞÞ expð2C1 TÞ;

! ðT þ hl Þhk

 ðOðhrz Þ

Iran J Sci Technol Trans Sci

in other words, rz

p

kel k1  Oðh Þ þ Oðh Þ:

ð11Þ

jdðtk Þ  dh ðtk Þj  jdðtk Þj:

Recall now (10), it yields the estimate jeðtl þ vhl Þj  Oðhrz Þ þ Oðhp Þ þ b rz  þ Oðhp Þ þ bTðOðh Þ þ Oðhp ÞÞ

l X

hi kei k1  Oðhrz Þ

i¼0

so kek1  Oðhminfrz;pg Þ:

  For e0 , let Kp :¼ maxj Lj 1 , from (8) and (11) we can write je0 ðtl þ vhl Þj  Kp kel k1 þhp Mpþ1 Kp  Kp ðOðhrz Þ þ Oðhp ÞÞ þ hp Mpþ1 Kp :

Therefore ke0 k1  Oðhminfrz;pg Þ: h The problem of tracking the breaking points as zeros of functions of the hðt; yðtÞÞ  f form can be reduced to the (easier and more practicable) observation of sign changes of the switching functions hðt; yðtÞÞ  fh . The next lemma shows that checking the sign of hðt; yðtÞÞ  fh while using Algorithm 3.1 also detects the breaking points. If a sign change indicates a breaking point, then it can be more closely approximated by locating the zero f^h by any iterative root-finding routine. Lemma 4.2 Assume the hypothesis (H) holds for some breaking point f^ of order z^ p þ 1 and let m be the (lefthanded) integer multiplicity of the zero f^ of hðt; yðtÞÞ  f . Furthermore, let jfh  fj ¼ Oðhr Þ for some r  1 with rz  p. Then for sufficiently small h   ^  fh  f^ ¼ Oðhminfp;rg=m Þ:

ðmÞ

ðaðtÞÞ ; aðtÞ 2 ðtk ; tk þ hÞ; m!

ð12Þ

hence dðtÞ ¼ Oðhm Þ. Due to the Lipschitz continuity of h with respect to its second variable, Theorem 4.1 and rz  p, we have jdðtÞ  dh ðtÞj  Lh jyðtÞ  uðtÞj þ jf  fh j ¼ Oðhminfrz;pg Þ þ Oðhr Þ ¼ Oðhminfrz;pg Þ:

This implies that d and dh have the same sign at tk and similarly at tk þ h. Due to changing sign of d on ½tk ; tk þ h, dh also changes sign on ½tk ; tk þ h . Since dh ðf^h Þ ¼ 0 and dh is continuous on ½tk ; tk þ h, then f^h 2 ½tk ; tk þ hÞ. Con^ 6¼ 0, for sufficiently small h, sidering the fact that dðmÞ ðfÞ there exists an open interval V I such that f^h ; f^ 2 V and  ðmÞ  d ðtÞ [ 0 for all t 2 V,  such that  ðmÞ  d ðtÞ  M [ 0; 8t 2 V:  ð14Þ From (13) we get    ^  dðfh Þ ¼ Oðhminfp;rg Þ:

ð15Þ

It follows from (12), (15) and inequality (14) that    m f ðf^h Þ ^  ¼ Oðhminfp;rg Þ: fh  f^ ¼ ðmÞ jd ðaÞj h To complete this section, we give the following theorem demonstrating the high order of convergence, when the breaking points are included in the mesh points. Theorem 4.2 Assume that (1) has continuity class p  1, and there exist a finite number of breaking points with order of continuity less than p þ 1 in ½a; T. Furthermore, hypothesis (H) holds and the breaking points are approximated by the method described in Sect. 3, then ky  uk1 ¼ Oðhp Þ: Proof The proof is given in Theorem 5.4, Feldstein and Neves (1984). h

Proof We assume f^ 2 ½tk ; tk þ h. Let us set dðtÞ ¼ hðt; yðtÞÞ  f and dh ðtÞ ¼ hðt; yðtÞÞ  fh . Hence, Taylor expansion implies ^md dðtÞ ¼ ðt  fÞ

Since h is sufficiently small, using dðtk Þ ¼ Oðhm Þ and (13) we have

5 Numerical Results In this section, we illustrate the theoretical results of the previous sections by the following example. The numerical experiment is implemented in Matlab. Example Consider the delay integro-differential equation Z yðtÞ y0 ðtÞ ¼ yðsÞds; t  2; 0

ð13Þ

/ðtÞ ¼1; t  2:

123

Iran J Sci Technol Trans Sci Table 1 Values of ky  uk1

p¼2

p¼4

h0

2.7757E-2

6.9974E-5

h0 2 h0 4 h0 8 h0 16

1.9329E-3

2.7757E-2

1.6593E-3

2.9826E-7

5.0162E-4

2.6932E-8

1.2509E-4

1.6452E-9

behavior changes at the breaking point f1 ¼ 2 þ ln 2 ¼ 2:6931, as we see a jump for p ¼ 4.

6 Conclusion We propose the piecewise collocation method in order to solve a class of state-dependent delay integro-differential equations. The convergence properties of the method is analyzed and the propagated discontinuities are included in the set of the mesh points. The numeric experiment confirms the theoretical findings of the paper and the robustness of the method.

References

Fig. 1 Results for p ¼ 2

Fig. 2 Results for p ¼ 4

If we assume hðt; yðtÞÞ :¼ yðtÞ\t, we can easily see that the exact solution is as follows: 8 et2 ; <  ½2; 2 þ ln 2  8 yðtÞ ¼ 2 þ ln 2; 2 þ ln ; :  lnð4et  e2 Þ; 3 The results are shown in Table 1 for different values of p and h with h0 ¼ 14 ln 83 . Figures 1 and 2 show error in the collocation points tn , for different values of h with p ¼ 2; 4. They indicate that error decays as h decreases. We observe that the error

123

Bellen A, Zennaro M (2003) Numerical methods for delay differential equations. Oxford University Press, Oxford Bellour A, Bousselsal M (2014) A Taylor collocation method for solving delay integral equations. Numer Algorithms 65(4):843–857 Brunner H (2004) Collocation methods for Volterra integral and related functional equations. Cambridge University Press, Cambridge De Gaetano A, Arino O (2000) Mathematical modelling of the intravenous glucose tolerance test. J Math Biol 40:136–168 Feldstein A, Neves KW (1984) High order methods for statedependent delay differential equations with nonsmooth solutions. SIAM J Numer Anal 21(5):844–863 Guglielmi N, Hairer E (2008) Computing breaking points in implicit delay differential equations. Adv Comput Math. doi:10.1007/ s10444-007-9044-5 Hale JK, Verduyn Lunel SM (1993) Introduction to functional differential equations. Springer, Berlin Hartung F (2001) Parameter estimation by quasi-linearization in functional differential equations with state-dependent delays: a numerical study. Nonlinear Anal 47(7):4557–4566 Hartung F, Herdman TL, Turi J (2000) Parameter identification in classes of neutral differential equations with state-dependent delays. Nonlinear Anal 39(3):305–325 Hauber R (1997) Numerical treatment of retarded differentialalgebraic equations by collocation methods. Adv Comput Math 7:573–592 Kuang Y (1993) Delay differential equations with applications in population dynamics. Academic Press, San Diego Marino S, Beretta E, Kirschner DE (2007) The role of delays in innate and adaptive immunity to intracellular bacterial infection. Math Biosci Eng 4(2):261–286 Neves KW, Feldstein A (1976) Characterization of jump discontinuities for state-dependent delay integral equations. J Math Anal Appl 56:689–707 Shakourifar M, Enright WH (2011) Reliable approximate solution of systems of Volterra integro-differential equations with timedependent delays. SIAM J Sci Comput. doi:10.1137/100793098 Vanani SK, Yldrm A, Soleymani F, Khan M (2012) A Numerical algorithm for solving nonlinear delay Volterra integral equations by means of homotopy perturbation method. Int J Nonlinear Sci Numer Simul. doi:10.1515/IJNSNS.2011.061