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, , and the delay time may even vanish when 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.
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 where is a time delay or lag which is a positive constant, is a prescribed initial function, and 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 in , 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 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, . 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 of 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 on -step predictor-corrector pair can be defined as where the step number from to and 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 and at two points and simultaneously using the information in the previous block at points , and . This process will continue on the next block until it reaches the end of the interval. In particular, the computed block has step size , while the previous block has the step size where the use of and in this method is for variable step size implementation.
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 and will be integrated over the interval and , respectively. While in Abdul Majid and Suleiman [6], the approximation of the solution is obtained by integrating over the interval and , 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 in our algorithm, the most recently calculated approximation which is is used instead of the initial approximation of .
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 for , and are and by using MAPLE, the formula in terms of can be written in the matrix form as the follows:
Letting us consider and substituting this value in (3) will produce the following first and second point of the corrector formula: 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 .
Now, we consider the implementation of mode where denotes the application of the predictor of order , denotes the application of the corrector of order , denotes the evaluation of the function , and denotes the number of iterations that is needed in a convergence test. The considered implementation can be described as follows: for and .
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 in 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, which is caused when the delay value is smaller than the step size, , where as illustrated in Figure 2.
3.2. Vanishing Lag
The vanishing lag may occur when the delay time vanishes as at the current step as illustrated in Figure 3.
Detailed information of vanishing lag can be described by considering the case of problem 1 in Section 7 as follows: In this problem, the vanishing lag occurs when the lag at . It can be seen that when we substitute the lag in (6), it yields where the delay term . 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, , 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.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
For 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, . Then, the solution of delay term can be obtained by extrapolating the points in the interval .
4. Stability of the Method
The general linear test equation for DDE is where and are complex and is a continuous function.
Definition 1. For the step size consider the following.(i)If and are real in (8), the region in the -plane is called the -stability region if for any the numerical solution of (8) satisfies as . The test equation for -stability is (ii)If and is complex in (8), the region in the -plane is called the -stability region if for any the numerical solution of (8) satisfies as . The test equation for -stability is
The corrector formulae of 2-point modified block method in (4) can be written as follows:
Applying (11) to the test equations in (9) and rearranging the equation to be equal to zero will give
where
and the matrices of , , and are depending on the step size ratio which has been formulated in corrector formula in (4).
For , Thus, the -stability polynomial, of 2-point modified block method, is given by When the same approaches are applied to the test equation (10), it will give the -stability polynomial, , as follows: By solving and , therefore the - and -stability regions are as shown in Figures 4 and 5.
5. Order and Error Constant
From (2), we associate that the difference operator is defined by where is an arbitrary function. Expanding the functions and as Taylor series about gives Substitute into (17) Collecting terms where
Definition 2. The multistep block method in (2) and the difference operator (17) are said to be of order , if and .
By applying the formulae in (4), we obtained
From Definition 2, we can conclude that the 2-point modified block method is of order five with the error constant
6. Variable Step Size Strategy
The developed algorithm starts by finding the values of starting point at , and by using the Euler method. The initial block in 2-point modified block method can be obtained by the approximation of the solutions and with the starting step size ratio and . For obtaining the next block, the above 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: where 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 and as follows: where and are the step size and the step size ratio in the previous block, respectively, and 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 .
The convergence test is done by the iteration of corrector in the second point as given below: By comparing the corrector formula of the method of order and the same corrector formulae of order at the second point , the local truncation error, , can be estimated as Finally, the numerical results are compared with the exact solutions and the maximum error of the mixed test can be defined as follows: where is the th component of the approximate , is the exact solution, is the number of equations in the system, SSTEP is the number of successful steps, and .
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 at ), see [3]). Consider Exact solution is
Problem 2 ((state dependent delay with vanishing lag at ), see [4]). Consider Exact solution:
Problem 3 (delay equation with vanishing lag, see [2]). Consider Exact solution is 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 , while SYSDEL code needs 483 total function calls and obtains maximum error .
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.
Overall, 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.
Notations
| TOL: | The prescribed tolerance |
| MTD: | Method employed |
| TS: | The total number of steps |
| FS: | The total number of failure steps |
| MAXE: | Maximum value of mixed error test of the computed solution |
| FNC: | Total function calls |
| DDEB5: | A developed code in this paper, which handled the small and vanishing lag by using 5-point Newton divided difference interpolation; the integration method is 2-point modified block method of order five with the variable step size implementation |
| SYSDEL: | A code developed by Karoui and Vaillancourt [2], which handled the vanishing lag by using the extrapolation of 3-point Hermite polynomial; the integration method is Runge-Kutta formula pair of order with the variable step size implementation as described in [11] |
| —: | No data provided in the reference. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge the financial support of Putra Grant (GP-IPS/2013/9390100) from Universiti Putra Malaysia and the scholarship of Academic Training Scheme (SLAI) received from Ministry of Education Malaysia.