# Solving Delay Differential Equations of Small and Vanishing Lag Using Multistep Block Method

Authors: Nurul Huda Abdul Aziz; Zanariah Abdul Majid; Fudziah Ismail
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper considers the numerical solution of delay differential equations for solving the problem of small and vanishing lag using multistep block method. This problem arises when the size of a delay value is smaller than the step size,x-τ<h, and the delay time may even vanish when τ→0 in a current step. The proposed approach that is based on interpolation of Newton divided difference has been implemented by adapting this problem to the multistep block method. In order to achieve the required accuracy, this approach considered the appropriate degree of interpolation polynomial in approximating the solution of delay term. The developed code for solving small and vanishing lag is done using C program and we called it as DDEB5. The P-stability and Q-stability of this method are also studied. Numerical results are presented and compared to the existing method in order to illustrate the efficiency of the proposed method.


---

## Body

## 1. Introduction

Delay differential equation (DDE) is one of the mathematical models that commonly possess the result in differential equations with time delay. In general, the unknown function of this derivative equation not only depends on the current value but also depends on the past value which is called a delay term. This equation can be given in the form of(1)y'x=fx,yx,yx-τ,x∈a,b,yx=ϕx,x∈-τ,a,
where τ is a time delay or lag which is a positive constant, ϕ(x) is a prescribed initial function, and y(x-τ) is the solution of delay term.There are two families of difficulties that may exist in the numerical solution of DDEs: the occurrences of derivative discontinuity along the integration interval and the small and vanishing lag where the delay term is smaller than current step size and the delay time may even vanish asτ→0 in (x-τ), respectively. In this paper, we will focus on the study of dealing with the second difficulty in the adaptation of multistep block method which is specifically 2-point modified block method.In the previous work, there are several authors that have investigated the numerical solution of DDEs with the problem of vanishing delay and small delay. For instance, Enright and Hu [1] developed an approach which combines an iteration scheme and interpolation technique that is based on two time steps for Runge-Kutta methods to handle the vanishing delay. Their main idea is to use all the information from the last step including the early stages of the current step in the interpolation technique.Karoui and Vaillancourt [2] developed a SYSDEL code which is based on the numerical method of Runge-Kutta for the formula pair of order (5,6) for solving the case of vanishing lag and asymptotically vanishing lag of (1). In their study of solving vanishing lag case, they have used interpolation or extrapolation of 3-point Hermite polynomial to approximate the solution of the delay term by depending on the location of delay time fall in the history queue or in a neighborhood of vanishing lag point, respectively. Then for the asymptotically vanishing lag case, they have used the approach of the first case as starter and then the delay equation is approximated by solving it as an ODE.The latest code for the numerical scheme of small, vanishing, and asymptotically vanishing delay in DDEs has come out in Yagoub et al. [3] in the name of HBODDE. This code is based on a hybrid variable-step variable-order 3-stage Hermite-Birkhoff-Obrechkoff ODE solver that has been adapted in DDEs. The delay values are computed by Hermite interpolation and the small delay deal with extrapolation.Hayashi [4] handled small and vanishing delay by proposing three algorithms of iterative scheme which are extrapolation, special interpolant, and iteration procedure with the adaptation of continuous Runge-Kutta method, while Neves and Thompson [5] handled small and vanishing delay by restricting the step size to be smaller and using extrapolation, respectively. All of the approaches that have been proposed are due to the one-step integration method where when the delay time falls in the current step, there is no any available approximant information that can be used for computing the delay term, (x-τ). This is lead to the use of extrapolation in their way to handle the small or vanishing lag.In the past eight years, the study of DDE in the numerical solution of multistep block method has gained some attention among the researchers. These methods were initially investigated in solving ODE and have shown the advantages of less computational works and manage to obtain good accuracy. For the previous work please see [6–8]. In the study of DDE involving block method, Ishak et al. [9] has developed the two-point block for solving delay differential equations by using six points of Lagrange interpolation to evaluate the delay solution, while San et al. [10] has investigated a coupled block method that can integrate the solution within two different blocks which are namely two-point two-step block and three-point two-step block method. In both works, they only solved a typical retarded DDE problem without any difficulty that may arise in delay differential equations. These situations allow us to fill the gap taking into account the second difficulty with the adaptation of 2-point modified block method.In this paper, we have developed a new DDEB5 code of C program for handling small and vanishing lag using the current approximate information provided from the approximation ofPECEs mode to evaluate the small and vanishing lag using Newton divided difference interpolation. The detail of this numerical scheme is described in Section 3 in conjunction with the summary of DDEB5 code.


## 2. Formulation of the Method

The multistep block method based onk-step predictor-corrector pair can be defined as
(2)Predictor:∑q=0k-1αq*yn+q=h∑q=0k-1β-q*fn-q,Corrector:∑q=0kαqyn+q=h∑q=0kβ2-qfn+2-q,
where the step number from q=0 to (k-1) and k refers to the order of the predictor and corrector, respectively.There are many types of methods that can be generated from (2). Specifically, in this paper, we will only focus on the 2-point modified block method. According to Figure 1, this 2-point modified block method will approximate the solution of yn+1 and yn+2 at two points xn+1 and xn+2 simultaneously using the information in the previous block at points xn-2,xn-1, and xn. This process will continue on the next block until it reaches the end of the interval. In particular, the computed block has step size 2h, while the previous block has the step size 2rh where the use of r and q in this method is for variable step size implementation.Figure 1
2-point modified block method.Basically, this 2-point modified block method is differing from the block method that has been proposed by Abdul Majid and Suleiman [6]. In this method, the approximation of these two solutions yn+1 and yn+2 will be integrated over the interval [xn,xn+1] and [xn+1,xn+2], respectively. While in Abdul Majid and Suleiman [6], the approximation of the solution is obtained by integrating over the interval [xn,xn+1] and [xn,xn+2], respectively. Our considered approach is called Gauss-Seidel method which is an iterative procedure that is based on a modification of the Jacobi method. In Gauss-Seidel method, the first approximation is computed in the same manner as in the Jacobi method. However, in computing the second approximation, the Gauss-Seidel method assumes that the new solution values are a better approximant to the solution than the initial values. In other words, to approximate the second point yn+2 in our algorithm, the most recently calculated approximation which is yn+1 is used instead of the initial approximation of yn.The formula of predictor and corrector in 2-point modified block method can be obtained by the derivation of Lagrange interpolation polynomial. For corrector formula, the interpolation points involved foryn-2,yn-1,yn,yn+1, and yn+2 are xn-2,fn-2,xn-1,fn-1,xn,fn,xn+1,fn+1,xn+2,fn+2 and by usingMAPLE, the formula in terms of r can be written in the matrix form as the follows:(3)1001yn+1yn+2=0010yn+1yn+2+0100yn-1yn+h240(672r3+144r2+940r4+320r5)(r+1)(r+2)(2r+1)r2(-21r3-3r2-50r4-40r5)(r+1)(r+2)(2r+1)r2-(-1632r3-1300r4-320r5-624r2)(r+1)(r+2)(2r+1)r2-(-549r3-610r4-200r5-147r2)(r+1)(r+2)(2r+1)r2×fn+1fn+2+h240(-176r-28-240r2)(r+1)(r+2)(2r+1)r2(1029r3+564r2+139r+14+790r4+200r5)(r+1)(r+2)(2r+1)r2-(-92-304r-240r2)(r+1)(r+2)(2r+1)r2-(46+230r4+40r5+501r3+516r2+251r)(r+1)(r+2)(2r+1)r2×fn-1fn+h2400(37r+14+15r2)(r+1)(r+2)(2r+1)r20(46+53r+15r2)(r+1)(r+2)(2r+1)r2fn-3fn-2.Letting us considerr=1 and substituting this value in (3) will produce the following first and second point of the corrector formula:
(4)r=1,1001yn+1yn+2=0010yn+1yn+2+0100yn-1yn+h720346-19646251fn+1fn+2+h720-74456106-264fn-1fn+h7200110-19fn-3fn-2.
The coefficients of the formula (3) are then recalculated whenever the step size changes in each integration of steps. These approaches avoid the storing of the coefficients at the start of the code that may give too many subroutines in the algorithm. The same way in getting the predictor formula is employed by involving the interpolation points xn-3,fn-3,xn-2,fn-2,xn-1,fn-1,andxn,fn.Now, we consider the implementation ofPECEs mode where P denotes the application of the predictor of order (k-1), C denotes the application of the corrector of order k, E denotes the evaluation of the function f, and s denotes the number of iterations that is needed in a convergence test. The considered implementation can be described as follows:
(5)P:ypixn+1=yxn+h∑q=03β-qzxn-q,ypixn+2=ypixn+1+h∑q=03β-qzxn-q,E:zpi(xn+1)=fxn+1,ypixn+1,ypixn+1-τ,zpi(xn+2)=fxn+2,ypixn+2,ypixn+2-τ,C:[j]ycxn+1=yxn+h∑q=04β2-qzxn+2-q,[j]ycxn+2=[j]ycxn+1+h∑q=04β2-qzxn+2-q,E:[j]zcxn+1=fxn+1,[j]ycxn+1,[j]ycxn+1-τ,[j]zcxn+2=fxn+2,[j]ycxn+2,[j]ycxn+2-τ,
for i=0 and j=1,2,…,s.


## 3. Handling Small and Vanishing Lag with Newton Divided Difference Interpolation

In solving delay differential equations, one must be aware of the presence of the delay value in order to achieve the smooth solution and desired accuracy in the numerical solution. During the integration of DDEs, the delay time may even lie in the previous step, current step, or in the next step. The main difficulty in this numerical solution may arise when the delay falls in the current step where no approximation information can be used to evaluate the delay term; in particular, when one-step integration method is applied. But it might be slightly different with 2-point modified block method, where the approximate solution at the current step has been obtained from the application of predictor inPECEs mode. Therefore, in this section, we will compute the solution of small and vanishing lag using Newton divided difference interpolation by taking five points to interpolate in the code.


### 3.1. Small Lag

The small lag occurs when the delay time falls in the current step,[xn,xn+1] which is caused when the delay value is smaller than the step size, ψn+1<hn+1, where ψn+1=xn+1-τn+1 as illustrated in Figure 2.Figure 2
Small lag in 2-point modified block method.


### 3.2. Vanishing Lag

The vanishing lag may occur when the delay time vanishes asτn+1→0 at the current step ψn+1* as illustrated in Figure 3.Figure 3
Vanishing lag in 2-point modified block method.Detailed information of vanishing lag can be described by considering the case of problem 1 in Section7 as follows:
(6)y'x=1-yexp⁡1-1x,x∈0.1,10.
In this problem, the vanishing lag occurs when the lag at x=1. It can be seen that when we substitute the lag in (6), it yields
(7)y'x=1-yexp⁡1-11.y'x=1-y1,
where the delay term ψn+1*=exp⁡1-(1/1)=1=xn+1. Basically, if this happens, there is no approximate solution that can be used in the interpolation to evaluate the delay term because the delay falls in the current mesh point.In the implementation of 2-point modified block method, the location of the delay term,(x-τ), is sought first to specify whether the delay lies as the small lag or vanishing lag. This DDEB5 code will detect and handle these cases automatically and hence identifying the points that would be involved in the interpolation. In order to get an accurate approximation in delay solution, the number of interpolation points must be chosen properly as the order of interpolation is one order higher than or equal to the order of integration method. Since our method is of order five, we choose five-point number of approximate solution that has been stored closest to the delay value to do the interpolation. The approach of this case can be summarize in Algorithm 1.Algorithm 1
POINT = number of points to be computed simultaneously in a block,
EQN = number of equations in a system,
DDE(n+d,X,Y,DX,DY) = subroutine function to approximate the delay term,
FN(n+d,X,XD,YD,XDP,YDP,Y,K) =  subroutine function of function evaluation.
(1)for  d=1  to POINT do
(2)for  k=1  to EQN do
(3)P: if d=1, then [i]yp(xn+d)=y(xn)+h∑q=03β-qz(xn-q)
ifd=2, then [i]yp(xn+d)=y(xn+1)+h∑q=03β-qz(xn-q)
(4)   end for
(5)  end for
(6)for  d=1  to POINT do
(7)E: subroutine DDE(n+d,X,Y,DX,DY)
if (-τ≤ψn+d≤a)
then y(ψn+d)=ϕ(ψn+d)
else
y(ψn+d)=Newton_Divided_Difference_Interpolation
evaluate FN(n+d,X,XD,YD,XDP,YDP,Y,K)
(8)  end for
(9)for  d=1  to POINT do
(10)for  k=1  to EQN do
(11)C: if d=1, then [j]yc(xn+d)=y(xn)+h∑q=04β2-qz(xn+2-q),
ifd=2, then [j]yc(xn+d)=y(xn+1)+h∑q=04β2-qz(xn+2-q),
(12)  end for
(13) end for
(14)for  d=1  to POINT do
(15)E: subroutine DDE(n+d,X,Y,DX,DY)
if(-τ≤ψn+d≤a)
then y(ψn+d)=ϕ(ψn+d)
else
y(ψn+d)=Newton_Divided_Difference_Interpolation
evaluate FN(n+d,X,XD,YD,XDP,YDP,Y,K)
(16) end forFor the case if the small and vanishing lags occur at starting point of this 2-point modified block method, the approach of linear extrapolation will be used by determining the solution at previous point,xn-hn. Then, the solution of delay term can be obtained by extrapolating the points in the interval xn-h,xn.


## 3.1. Small Lag

The small lag occurs when the delay time falls in the current step,[xn,xn+1] which is caused when the delay value is smaller than the step size, ψn+1<hn+1, where ψn+1=xn+1-τn+1 as illustrated in Figure 2.Figure 2
Small lag in 2-point modified block method.


## 3.2. Vanishing Lag

The vanishing lag may occur when the delay time vanishes asτn+1→0 at the current step ψn+1* as illustrated in Figure 3.Figure 3
Vanishing lag in 2-point modified block method.Detailed information of vanishing lag can be described by considering the case of problem 1 in Section7 as follows:
(6)y'x=1-yexp⁡1-1x,x∈0.1,10.
In this problem, the vanishing lag occurs when the lag at x=1. It can be seen that when we substitute the lag in (6), it yields
(7)y'x=1-yexp⁡1-11.y'x=1-y1,
where the delay term ψn+1*=exp⁡1-(1/1)=1=xn+1. Basically, if this happens, there is no approximate solution that can be used in the interpolation to evaluate the delay term because the delay falls in the current mesh point.In the implementation of 2-point modified block method, the location of the delay term,(x-τ), is sought first to specify whether the delay lies as the small lag or vanishing lag. This DDEB5 code will detect and handle these cases automatically and hence identifying the points that would be involved in the interpolation. In order to get an accurate approximation in delay solution, the number of interpolation points must be chosen properly as the order of interpolation is one order higher than or equal to the order of integration method. Since our method is of order five, we choose five-point number of approximate solution that has been stored closest to the delay value to do the interpolation. The approach of this case can be summarize in Algorithm 1.Algorithm 1
POINT = number of points to be computed simultaneously in a block,
EQN = number of equations in a system,
DDE(n+d,X,Y,DX,DY) = subroutine function to approximate the delay term,
FN(n+d,X,XD,YD,XDP,YDP,Y,K) =  subroutine function of function evaluation.
(1)for  d=1  to POINT do
(2)for  k=1  to EQN do
(3)P: if d=1, then [i]yp(xn+d)=y(xn)+h∑q=03β-qz(xn-q)
ifd=2, then [i]yp(xn+d)=y(xn+1)+h∑q=03β-qz(xn-q)
(4)   end for
(5)  end for
(6)for  d=1  to POINT do
(7)E: subroutine DDE(n+d,X,Y,DX,DY)
if (-τ≤ψn+d≤a)
then y(ψn+d)=ϕ(ψn+d)
else
y(ψn+d)=Newton_Divided_Difference_Interpolation
evaluate FN(n+d,X,XD,YD,XDP,YDP,Y,K)
(8)  end for
(9)for  d=1  to POINT do
(10)for  k=1  to EQN do
(11)C: if d=1, then [j]yc(xn+d)=y(xn)+h∑q=04β2-qz(xn+2-q),
ifd=2, then [j]yc(xn+d)=y(xn+1)+h∑q=04β2-qz(xn+2-q),
(12)  end for
(13) end for
(14)for  d=1  to POINT do
(15)E: subroutine DDE(n+d,X,Y,DX,DY)
if(-τ≤ψn+d≤a)
then y(ψn+d)=ϕ(ψn+d)
else
y(ψn+d)=Newton_Divided_Difference_Interpolation
evaluate FN(n+d,X,XD,YD,XDP,YDP,Y,K)
(16) end forFor the case if the small and vanishing lags occur at starting point of this 2-point modified block method, the approach of linear extrapolation will be used by determining the solution at previous point,xn-hn. Then, the solution of delay term can be obtained by extrapolating the points in the interval xn-h,xn.


## 4. Stability of the Method

The general linear test equation for DDE is(8)y'x=λyx+μyx-τ,x≥x0,y(x)=ϕ(x),-τ≤x≤x0,
where λ and μ are complex and ϕ(x) is a continuous function.Definition 1.
For the step sizeh consider the following.(i)
Ifλ and μ are real in (8), the region RP in the (H1,H2)-plane is called the P-stability region if for any (H1,H2)∈RP the numerical solution of (8) satisfies y(xn)→0 as xn→∞. The test equation for P-stability is
(9)y'x=λyx+μyx-τ,x≥x0,yx=ϕx,-τ≤x≤x0.(ii)
Ifλ=0 and μ is complex in (8), the region RQ in the (H2)-plane is called the Q-stability region if for any (H2)∈RQ the numerical solution of (8) satisfies y(xn)→0 as xn→∞. The test equation for Q-stability is
(10)y'x=μyx-τ,x≥x0,yx=ϕx,-τ≤x≤x0.The corrector formulae of 2-point modified block method in (4) can be written as follows:
(11)A2YN+2=A1YN+1+h∑i=02BiFN+i.Applying (11) to the test equations in (9) and rearranging the equation to be equal to zero will give
(12)A2YN+2=A1YN+1+h∑i=02BiλYN+i+μYN+i-m.A2-H1B2YN+2-A1+H1B1YN+1-H1B0YN-H2∑i=02BiYN+i-m=0,where(13)A2=10-11,A1=0100,YN+2=yn+1yn+2,YN+1=yn-1yn,YN=yn-3yn-2,FN+2=fn+1fn+2,FN+1=fn-1fn,FN=fn-3fn-2,and the matrices ofB0, B1, and B2 are depending on the step size ratio r which has been formulated in corrector formula in (4).Forr=1,
(14)B2=1720346-19646251,B1=1720-74456106-264,B0=17200110-19.
Thus, the P-stability polynomial, π(ς) of 2-point modified block method, is given by
(15)πς=det⁡-H1B0tm-H2∑i=02BitiA2-H1B2t2+m-A1+H1B1t1+m-H1B0tm-H2∑i=02Biti.
When the same approaches are applied to the test equation (10), it will give the Q-stability polynomial, ψ(ς), as follows:
(16)ψ(ς)=det⁡A1t2+m-A0t1+m-H2∑i=02Biti.
By solving π(ς)=0 and ψ(ς)=0, therefore the P- and Q-stability regions are as shown in Figures 4 and 5.Figure 4
P-stability of 2-point modified block method.Figure 5
Q-stability of 2-point modified block method.


## 5. Order and Error Constant

From (2), we associate that the difference operator L is defined by
(17)Ly(x);h=∑q=0kαqyn+q-hβ2-qfn+2-q=∑q=0kαqy(x+qh)-hβ2-qy'x+2-qh,
where y(x)∈C1a,b is an arbitrary function. Expanding the functions y(x+qh) and y'(x+2-qh) as Taylor series about x gives
(18)yx+qh=yx+qhy'x+qh22!y′′x+qh33!y′′′x+⋯,y'x+2-qh=y'x+2-qhy′′x+2-qh22!y′′′x+2-qh33!y4x+⋯.
Substitute into (17)(19)Lyx;h=∑q=0kαqyx+qhy'x+qh22!y′′x+qh33!y′′′x+⋯-hβ2-q+2-qh22!y'''xy'x+2-qhy′′x+2-qh22!y′′′x+2-qh33!y4x+⋯.
Collecting terms
(20)Lyx;h=∑q=0k+h3q3αq3!-2-q2β2-q2!y'''x+⋯αqyx+hqαq-β2-qy'x+h2q2αq2!-2-qβ2-qy′′x+h3q3αq3!-2-q2β2-q2!y′′′x+⋯,
where
(21)C0=∑q=0kαq,C1=∑q=0kqαq-β2-q,C2=∑q=0kq2αq2!-2-qβ2-q,⋮Cr=∑q=0kqrαqr!-2-qr-1β2-qr-1!.Definition 2.
The multistep block method in (2) and the difference operator (17) are said to be of order P, if C0=C1=⋯=CP=0 and CP+1≠0.
By applying the formulae in (4), we obtained
(22)C0=C1=C2=C3=C4=C5=0,C6=yn+1yn+2=111440-3160≠0.
From Definition 2, we can conclude that the 2-point modified block method is of order five with the error constant
(23)C6=111440-3160T≠0.


## 6. Variable Step Size Strategy

The developed algorithm starts by finding the values of starting point atxn-2,xn-1, and xn by using the Euler method. The initial block in 2-point modified block method can be obtained by the approximation of the solutions yn+1 and yn+2 with the starting step size ratio r=1 and q=1. For obtaining the next block, the above PECEs approaches will be repeated until it reaches the end of the interval.In order to achieve the desired accuracy and the most optimal total steps in the whole interval, we use the Runge-Kutta Fehlberg variable step size strategy which has been introduced by [11]. If the integration step is successful, then the new step size in the next step should be determined by using the following formula:
(24)hnew=C×TOLEk-11/4,
where C=0.5 is a safety factor. The code then will recalculate the coefficient of the formula whenever the step size changes by using the step size ratio r and q as follows:
(25)r=holdh,q=roldholdh,
where hold and rold are the step size and the step size ratio in the previous block, respectively, and h is the current step size that is used in the computed block. In the case of failure step, the step size will be half of the hold.The convergence test is done by the iteration of corrector in the second point as given below:(26)yn+2(s+1)-yn+2(s)<0.1×TOL.
By comparing the corrector formula of the method of order k and the same corrector formulae of order k-1 at the second point xn+2, the local truncation error, Ek-1, can be estimated as
(27)Ek-1=yn+2(k)-yn+2(k-1).
Finally, the numerical results are compared with the exact solutions and the maximum error of the mixed test can be defined as follows:
(28)MAXE=max⁡1≤i≤SSTEPmax⁡1≤i≤Nyix-yxixA+Byxix,xi∈a,b,
where (yi)x is the xth component of the approximate y, y(xi) is the exact solution, N is the number of equations in the system, SSTEP is the number of successful steps, and A=1,B=1.


## 7. Numerical Results

In this section, we have tested three problems of vanishing lag delay differential equations in order to show the efficiency and reliability on the performance of 2-point modified block method where all these implementations have been tested by using the C program. The following notations are used in the notations section.Problem 1 ((vanishing lag atx=1), see [3]).
Consider(29)y'x=1-yexp⁡1-1x,x∈0.1,10,yx=ln⁡x,x∈0,0.1.
Exact solution is
(30)yx=ln⁡x.Problem 2 ((state dependent delay with vanishing lag atx=1), see [4]).
Consider(31)y1′x=y2(x),x∈[0.1,5],y2′x=-y2e1-y2xy22xe1-y2(x),x∈[0.1,5],y1x=ln⁡x,x∈0,0.1,y2x=1x,x∈0,0.1.
Exact solution:
(32)y1x=ln⁡x,x∈0.1,5,y2x=1x,x∈0.1,5.Problem 3 (delay equation with vanishing lag, see [2]).
Consider(33)y'x=x4-3x5+xln⁡⁡x-x-3+x-x-3-3yx-x-3,x∈2,30,yx=ln⁡x+x-3,x∈1.5,2.
Exact solution is
(34)yx=ln⁡x+x-3.
From the numerical results that are tabulated in Table 1 to Table 3, it can be observed that DDEB5 code gave better accuracy of maximum error with less number of function calls in the prescribed tolerances. This is also depicted clearly in the graphs shown in Figures 6, 7, and 8. For example at tolerance 10−6 in Table 3, the DDEB5 code only needs 157 total function evaluations and obtains a good accuracy of 7.17E-7, while SYSDEL code needs 483 total function calls and obtains maximum error 1.20E-4.Table 1
Numerical results for Problem1.
TOL
MTD
TS
FS
MAXE
AVERR
FNC
10−4
DDEB5
27
0
6.33
E
-
05
3.13
E
-
05
117
SYSDEL
—
—
4.05
E
-
03
—
498
10−6
DDEB5
59
0
5.61
E
-
07
2.83
E
-
07
245
SYSDEL
—
—
6.69
E
-
05
—
568
10−8
DDEB5
136
0
4.86
E
-
09
2.51
E
-
09
553
SYSDEL
—
—
1.07
E
-
07
—
631
10−10
DDEB5
329
0
4.18
E
-
11
2.27
E
-
11
1325
SYSDEL
—
—
1.83
E
-
09
—
1030
10−12
DDEB5
810
0
3.67
E
-
13
2.12
E
-
13
3249
SYSDEL
—
—
5.81
E
-
12
—
1947Figure 6
Comparison of total function calls versus maximum error for Problem1.Figure 7
Comparison of total function calls versus maximum error for Problem2.Figure 8
Comparison of total function calls versus maximum error for Problem3.At tolerance 10−6 in Table 2, it can be seen that both codes have a comparable maximum error, but DDEB5 required less total function evaluations in the integration. At tolerance 10−8 to 10−12, SYSDEL achieved better accuracy, but it needs more function evaluations at each tolerance compared to DDEB5.Table 2
Numerical results for Problem2.
TOL
MTD
TS
FS
MAXE
AVERR
FNC
10−4
DDEB5
34
0
2.04
E
-
04
9.10
E
-
05
145
SYSDEL
—
—
—
—
—
10−6
DDEB5
68
0
1.89
E
-
06
7.61
E
-
07
281
SYSDEL
—
—
1.85
E
-
06
—
553
10−8
DDEB5
160
0
2.00
E
-
08
8.33
E
-
09
649
SYSDEL
—
—
5.05
E
-
09
—
959
10−10
DDEB5
387
0
2.07
E
-
10
8.78
E
-
11
1557
SYSDEL
—
—
2.85
E
-
11
—
1946
10−12
DDEB5
957
0
2.10
E
-
12
9.03
E
-
13
3837
SYSDEL
—
—
1.11
E
-
13
—
4214Table 3
Numerical results for Problem3.
TOL
MTD
TS
FS
MAXE
AVERR
FNC
10−4
DDEB5
22
0
3.21
E
-
05
7.82
E
-
06
97
SYSDEL
—
—
—
—
—
10−6
DDEB5
37
0
7.17
E
-
07
2.59
E
-
07
157
SYSDEL
—
—
1.20
E
-
04
—
483
10−8
DDEB5
80
0
1.16
E
-
08
3.69
E
-
09
329
SYSDEL
—
—
9.39
E
-
07
—
588
10−10
DDEB5
188
0
1.23
E
-
10
3.85
E
-
11
761
SYSDEL
—
—
8.23
E
-
09
—
1099
10−12
DDEB5
459
0
1.05
E
-
12
3.04
E
-
13
1845
SYSDEL
—
—
7.96
E
-
11
—
2233Overall, it can be conclude that DDEB5 code has its own advantages in handling small and vanishing lag. This could be justified by the fact that the approach of Newton divided difference interpolation has given better accuracy in approximating delay value compared to extrapolation. The lesser number of function calls also has shown that DDEB5 code does not need more iteration in order to check the convergence of the solutions.


## 8. Conclusion

The approach of Newton divided difference interpolation in 2-point modified block method is proposed for solving the case of small and vanishing lag of delay differential equation. The numerical results proved that our developed DDEB5 code is reliable to handle the related difficulty that may exist in DDE. The DDEB5 code also has its own advantages when it can detect and automatically treat the small and vanishing lag in every step of integration without requiring any user guidance. It also has shown the capability in solving a system of state dependent vanishing lag as well.



---

*Source: 348912-2014-10-29.xml*