In this paper, an implicit 2-point Block Backward Differentiation formula (BBDF) method was 

considered for solving Delay Differential Equations (DDEs). The method was implemented by using a 

constant stepsize via Newton Iteration. This implicit block method was expected to produce two points 

simultaneously. The efficiency of the method was compared with the existing classical 1-point Backward 

Differentiation Formula (BDF) in terms of execution time and accurac

Heng, S. C.

Ibrahim, Zarina Bibi

Ismail, Fudziah

Suleiman, Mohamed

Universiti Putra Malaysia Institutional Repository

Pertanika J. Sci. & Technol. 21 (1): 37 - 44 (2013)SCIENCE & TECHNOLOGYJournal homepage: http://www.pertanika.upm.edu.my/ISSN: 0128-7680  © 2013 Universiti Putra Malaysia Press.1 0'( ) ( , , ( ),... ( ),     ny t f t y y t y t t tτ τ= − − ≥      (1)0( ) ( ),    y t t t tφ= ≤                                        (2)where ( )tφ is the initial function, and iτ  is the delay function with 1, 2,..., .i n=  According to Bellen and Zennaro (2003), the delays iτ  or Solving Delay Differential Equations by Using Implicit 2-Point Block Backward Differentiation FormulaHeng, S. C.*, Ibrahim, Z. B., Suleiman, M. and Ismail, F.Mathematics Department, Faculty of Science, Universiti Putra Malaysia, 43400 Serdang, Selangor, MalaysiaABSTRACTIn this paper, an implicit 2-point Block Backward Differentiation formula (BBDF) method was considered for solving Delay Differential Equations (DDEs). The method was implemented by using a constant stepsize via Newton Iteration. This implicit block method was expected to produce two points simultaneously. The efficiency of the method was compared with the existing classical 1-point Backward Differentiation Formula (BDF) in terms of execution time and accuracy.Keywords: Block Backward Differentiation formula, Delay Differential Equations, InterpolationsINTRODUCTIONIn the recent years, the Mathematics society is gradually shifting their interest into numerical treatment of DDEs due to its ability in Mathematical modelling of processes in various applications as it provides the best, if not the only realistic simulation of the observed phenomena (Ismail et al., 2002). This is because many physical systems possess the feature of having a delayed respond in input conditions so that the rate at which processes occur depends not only on the current state of the systems but also the past states (Ismail 1999).Some of the known application areas of DDEs include the fields of Engineering, Biology and Economy, such as in mixing of liquid, population growth, prey-predator population model and electrodynamics (Driver 1976). The general form of the first order DDE is as follows:Article history:Received: 18 March 2011Accepted: 13 September 2011E-mail addresses: cerylyn@yahoo.com (Heng, S. C.), zarina@math.upm.edu.my (Ibrahim, Z. B.), msuleiman@science.upm.edu.my (Suleiman, M.) fudziah@science.upm.edu.my (Ismail, F.)*Corresponding AuthorHeng, S. C., Ibrahim, Z. B., Suleiman, M. and Ismail, F.38 Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)lags as known to some, and which are always non-negative, can be divided into three different situations according to the complexity of the phenomenon. The delays may be just constant, which is referred to as the constant delay case, or the function of t, where iτ  = iτ (t) which is referred to as the variable or time dependant delay case, or even the function of t and y itself, where iτ  = ( )( ),i t y tτ , which is referred to as the state dependant delay case.Currently, most of the numerical methods for solving ordinary differential equations (ODEs) can be adapted to give corresponding techniques in solving DDEs. For example, stiff DDEs are recognized through the behaviour of ODEs. This is because of the solutions in the systems of DDEs, which may contain strongly damped components that rapidly approach equilibrium states (Roth, 1980). In general, the largest source of stiff DDEs is the stiffness in the ODEs components that are without the delay term. The stiffness of a linear system of ODEs ¢ = y Ay  is expressed as ratio ( )max | Re | min |Re |i i iil l  where il  are eigenvalues of A. For this reason, there are actually many numerical methods proposed for solving DDEs. The range of the methods comprises of one-step methods, as in Euler method and Runge-Kutta methods, as well as multistep methods and block implicit methods (Ismail, 1999). Among the proposed methods, the family of Runge-Kutta is the most famous numerical method used to solve DDEs.Research has shown various families of Runge-Kutta methods and interpolation techniques which are used to solve DDEs. Among other, Oberle and Pesch (1981) developed two numerical methods known as the Runge-Kutta-Fehlberg methods of orders 4 and 7 to solve DDEs. Hermite interpolation is used to approximate the delay term. Ismail et al. (2002) solved delay differential equations by using embedded Singly Diagonally Implicit Runge–Kutta method and the delay term was obtained by using Hermite interpolation.While most numerical methods to determine the solution for DDEs are based on the Runge-Kutta formulas, some opted for Adams and backward differentiation formula (BDF) methods to solve them. For instance, Bocharov et al. (1996) considered the application of the linear multistep methods for the numerical solution of initial value problem for stiff delay differential equations with several constant delays. As for the approximation of delayed variables, Bocharov et al. (1996) used Nordsieck’s interpolation technique.BLOCK BACKWARD DIFFERENTIATION FORMULAThe purpose of this paper is to solve stiff DDEs using 2-point block backward differentiation formula (BBDF) method derived by Ibrahim et al. (2007). The corrector formulas are formulated as follows:1 1 2 11 22 23 3+ - + +=- + - +n n n n ny y y y hf                                     (3)2 1 1 22 9 18 611 11 11 11+ - + += - + +n n n n ny y y y hf                                   (4)Solving Delay Differential Equations by Using Implicit 2-Point Block Backward Differentiation Formula39Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)The numerical results were then compared to the classical 1-point backward differentiation formula (BDF). The corrector formula (Lambert, 1993) is given as:1 2 1 12 9 18 611 11 11 11+ - - += - + +n n n n ny y y y hf                                  (5)In general, most numerical methods for solving differential equations produce only one new approximation value. However, the implicit block method can produce two new values simultaneously at each step. Thus, it is logically safe to presume that this particular method can reduce the timing of the calculation and its computational cost as well.IMPLEMENTATION OF THE METHODNewton iteration is implemented in the method to evaluate the values of 1+ny  and 2+ny  for solving equation [1]. The values for both 1+ny  and 2+ny  at ( 1)+ith iterative are given as ( 1)1++iny  and ( 1)2++iny . According to Ibrahim et al. (2007), the Newton iteration takes the form of( )( )1 1)( 1) ( )1 1 '1 1( )( )+++ ++= -ini in n inF yy yF y                                         (6)and( )( )2 2)( 1) ( )2 2 '2 2( )( )+++ ++= -ini in n inF yy yF y                                           (7)Hence,( ) ( )( )1 ( ) ( ) ( )1 1 1 1 2 1 1121 2 23+++ + + + ++æ ö¶ ÷ç ÷- - =- - + +ç ÷ç ÷ç ¶è øi i i i inn n n n nnfh y y y y hfyV           (8)( ) ( )( )1 ( ) ( ) ( )2 2 2 1 2 2 226 18 6111 11 11+++ + + + ++æ ö¶ ÷ç ÷- - = - + +ç ÷ç ÷ç ¶è øi i i i inn n n n nnfh y y y y hfyV           (9)where ¶¶fy  is a Jacobian and 1V , 2V  are the back values.Let( 1) ( 1) ( )1 1 1+ ++ + += -i i in n ne y y                                             (10) ( 1) ( 1) ( )2 2 2+ ++ + += -i i in n ne y y                                            (11)Heng, S. C., Ibrahim, Z. B., Suleiman, M. and Ismail, F.40 Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)Then, equations [8] and [9] can be rewritten in the matrix form, as follows:1( 1) ( ) ( )1 11 1 1( 1) ( ) ( )22 2 2 222 21 2 ( ) 2 013 3618 018 6 11 ( ) 111111 11ni i in n n ni i in n n nnfhy e y fhf e y fhyVV+++ + + +++ + + ++é ù¶ é ùê ú- é ùê ú- -ê ú é ù é ù é ù é ù¶ ê úê úê ú ê ú ê ú ê ú ê úê ú= + +ê úê ú ê ú ê ú ê ú ê úê ú¶ ê ú ë ûê ú ë û ë û ë û- ê ú- - ê ú ë ûê ú ê ú¶ ë ûê úë û(12)In order to solve the above matrix equation, LU decomposition, which is a matrix decomposition that rewrites the matrix in the upper triangular matrix and lower triangular matrix was used (William et al., 2007). Let112221 2318 6111 11nnnnfhyfhAy++++é ùæ ö¶ ÷çê ú÷- ç ÷ê úç ÷ç¶è øê úê úæ öê ú¶ ÷ç ÷ê ú- - ç=÷çê ú÷ç¶è øë û                                      (13)( )( )( )( )1 1 122 22 2 013618 01 1111i in ni in ny fB hy fVV+ ++ +é ù é ùê ú- - é ù é ù é ùê úê ú ê ú ê ú ê úê ú= + +ê ú ê ú ê ú ê úê úê ú ê ú ê ú ë ûë û ë û- ê úê ú ë ûê úë û                                    (14)( 1)1( 1)2inineEe++++é ùê ú= ê úë û                                                           (15)A E B× =                                                            (16)Suppose that matrix A is a product of the two matrices,L U A× =                                                             (17)where L  is the lower triangular and U  is the upper triangular. Thus,( ) .L U E B× × =                                                         (18)Let .X U E= ×                                                                 (19)Hence,.L X B× =                                                                 (20)Solving Delay Differential Equations by Using Implicit 2-Point Block Backward Differentiation Formula41Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)and then solving,U E X× =                                                            (21)to get the matrix value of E . Finally, the approximation values of 1ny +  and 2ny +  can be obtained. Furthermore, the existence of the delay term is the main difference between ODEs and DDEs. Thus, in order to get the numerical solutions for DDEs, the delay term needs to be solved in advance. For a situation where 0tα ≤ , the initial function ( )tφ  is used to calculate the delay term, ( )y α . Otherwise, the delay term is calculated by using Newton Divided Difference interpolation.Meanwhile, the formula to calculate the delay term (Richard et al., 1993) can be written as:( ) [ ] [ ]( ) ( )10 0 1 1, , ,nn k o kkP x z x z x x x x x x x -== + ¼ - ¼ -å                        (22)where [ ]0 0( )z x z x= .Equation [22] is known as the Newton’s interpolatory divided-difference formula. Normally, the delay terms are obtained by interpolation on the values of the chosen support points. In order to obtain these support points, the delay argument should be close to the centre of the points. However, if the delay argument under consideration is outside the points, it will be treated as a special case and it can only be obtained through extrapolation.PROBLEM TESTEDThe followings are some of the tested problems to get the numerical results. All the tested problems are stiff DDEs. The numerical results are obtained by using C language with constant stepsize, H.Problem 1:( ) ( ) ( )( )' 1000 (ln 1000 1 ) , 0 3,    y t y t y t t=- + - - £ £( )       , . 0ty t e t-= £Exact solution is ( ) .ty t e-=Problem 2:( ) ( ) ( ) ( )' 3 31000 997 1 1000 99           7 , 0 3,y t y t e y t e t- -=- + - + - £ £( ) ( )1 exp 3 , 0.y t t t= + - £Heng, S. C., Ibrahim, Z. B., Suleiman, M. and Ismail, F.42 Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)Exact solution is ( ) ( )1 exp 3 .y t t= + -Problem 3:( ) ( ) ( ) ( )25' 24 1          , 0 3,y t y t e y t t-=- - - £ £( ) ( )25 , 0.ty t e t-= £Exact solution is ( ) .ty t e-=NUMERICAL RESULTS AND DISCUSSIONIn this section, the numerical results showed the performance of block method in solving the first order stiff DDEs. All the results were obtained using the Microsoft Visual Studio 2010 software on a personal computer with Intel Core2 Duo CPU. The abbreviations used are stated in Table 1.TABLE 1: List of Abbreviations Notation DescriptionH Stepsize1BDF 1-point BDF method2BBDF 2-point BBDF methodTS total number of stepsMAXE Maximum errorTIME CPU time in secondsThe errors are calculated by using the following formulas:i i it exact approximateerr y y= -                                              (23)The maximum error, MAXE, is defined as:( )1MAXE max .iti TS err£ £=                                                   (24)The numerical results are given in Table 2 to Table 4. In the results, the notation 9.75159e-004 represents the value of 9.75159 X 10-4. From the numerical results, it is clear that the total number of steps in solving 2-point BBDF is half of the steps needed to solve the 1-point BDF. The execution time for the 2-point BBDF is reduced up to 77% as compared to the 1-point BDF. Meanwhile, the maximum errors for 2-point BBDF method are better than the 1-point BDF for all the tested problems.Solving Delay Differential Equations by Using Implicit 2-Point Block Backward Differentiation Formula43Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)TABLE 2: Numerical results for Problem 11H METHOD TS MAXE TIME(seconds)10-2 1BDF2BBDF3001509.75159e-0043.80236e-0040.008367480.0019086310-3 1BDF2BBDF300015003.22629e-0072.61239e-0070.011032560.0045113410-4 1BDF2BBDF30000150001.45149e-0081.33568e-0080.037421930.0307998510-5 1BDF2BBDF3000001500001.73433e-0101.57914e-0100.294731640.278207631The time of the execution of the 2-point BBDF for Problem 1 is about 77%, 59%, 17%, and 5% faster for the step-sizes of 10-2, 10-3, 10-4, and 10-5, respectively.  The smaller step-size obtained a small difference in the execution time.  This is due to the total calculation of LU decomposition that requires matrix in 2BBDF.TABLE 3: Numerical results for Problem 22H METHOD TS MAXE TIME(seconds)10-2 1BDF2BBDF3001508.73499e-0033.40594e-0030.008435240.0029511010-3 1BDF2BBDF300015002.88750e-0062.33921e-0060.011679290.0071567810-4 1BDF2BBDF30000150001.30592e-0071.20176e-0070.042872570.0356464910-5 1BDF2BBDF3000001500001.56085e-0091.42118e-0090.326525070.306108352The time of the execution of the 2 point BBDF for Problem 2 is about 65%, 38%, 16%, and 6% faster for the respective step-sizes of 10-2, 10-3, 10-4, and 10-5TABLE 4: Numerical results for Problem 33H METHOD TS MAXE TIME(seconds)10-2 1BDF2BBDF3001504.63438e-0024.41121e-0020.008573050.0029630310-3 1BDF2BBDF200015001.00359e-0039.28409e-0040.011457940.0055307110-4 1BDF2BBDF20000150001.10170e-0059.97473e-0060.036800990.0297635410-5 1BDF2BBDF2000001500001.11197e-0071.00464e-0070.268148250.266494463The time of the execution of the 2-point BBDF for Problem 3 is about 65%, 51%, 19%, and 1% faster for the stepsizes of 10-2, 10-3, 10-4, and 10-5, respectivelyHeng, S. C., Ibrahim, Z. B., Suleiman, M. and Ismail, F.44 Pertanika J. Sci. & Technol. 21 (1): 283 - 298 (2013)CONCLUSIONIn this paper, it was observed that the 2-point BBDF achieved better results in terms of the accuracy of the method as well as the execution time. Therefore, it can be concluded that the 2-point BBDF is suitable for solving stiff DDEs.ACKNOWLEDGEMENTSThis work was fully supported by the Institute of Mathematical Research, Universiti Putra Malaysia, under the Fundamental Scheme Research Grant (Project code: 01-09-09-681FR).REFERENCESBellen, A., & Zennaro, M. (2003). Numerical Methods for Delay Differential Equations. New York: Oxford University Press.Bocharov, G. A., Marchuk, G. I., & Romanyukha, A. A. (1996). Numerical solution by LMMs of Stiff Delay Differential systems modelling an Immnune Response. Numer. Math., 73, 131-148.Driver, R. D. (1976). Ordinary and Delay Differential Equations. New Year: Springer-Verlag.Ibrahim, Z. B., Othman, K. I., & Suleiman, M. (2007). Implicit r-point Block Backward Differentiation Formula for Solving First Order Stiff ODEs. Applied Mathematics and Computation, 186, 558-565.Ismail, F. (1999). Numerical Solution of Ordinary and Delay Differential Equations by Runge-Kutta Type Methods (Ph.D Thesis dissertation). Universiti Putra Malaysia, Malaysia.Ismail, F., Al-Khasawneh, R. A., Lwin, A. S., & Suleiman, M. (2002). Numerical Treatment of Delay Differential Equations by Runge-Kutta Method Using Hermite Interpolation. Matematika, 18(2), 79-90.Lambert, J. D. (1993). Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. Chichester, England: John Wiley and Sons.Oberle, H. J., & Pesch, H. J. (1981). Numerical Treatment of Delay Differential Equations by Hermite Interpolation. Numer. Math., 37, 235-255.Richard, L., Burden, J., & Douglas, F. (1993). Numerical Analysis. Boston: PWS-Kent Pub. Co.Roth, M. G. (1980). Difference Methods for Stiff Delay Differential Equations (Master Thesis dissertation). University of Illinois, Illinois.William, H. P., Saul, A. T., William, T. V., & Brian, P. F. (2007). Numerical Recipes: The Art of Scientific Computing. New York: Cambridge University Press.

Solving delay differential equations by using implicit 2-point block backward differentiation formula

http://psasir.upm.edu.my/30300/7/Solving%20Delay%20Differential%20Equations%20by%20Using%20Implicit%202-Point%20Block%20Backward%20Differentiation%20Formula.pdf

Solving delay differential equations by using implicit 2-point block backward differentiation formula

Abstract

Similar works

Full text

Available Versions

Universiti Putra Malaysia Institutional Repository