A New family of methods for solving delay differential equations
AA New family of methods for solving delay differential equationsYogita Mahatekar, Pallavi S. Scindia
Department of Mathematics, College of Engineering Pune, Pune - 411005, India,[email protected], [email protected]
Abstract
In the present paper, we introduce a new family of θ − methods for solving delay differential equa-tions. New methods are developed using a combination of decomposition technique viz. newiterative method proposed by Daftardar Gejji and Jafari and existing implicit numerical methods.Using Butcher tableau, we observed that new methods are non Runge-Kutta methods. Further,convergence of new methods is investigated along with its stability analysis. Applications to va-riety of problems indicates that the proposed family of methods is more efficient than existingmethods. Keywords : New iterative method (NIM), Delay differential equation.
A delay differential equation (DDE) is a differential equation in which state function is given interms of value of the function at some previous times. Introduction of delay term in modelling al-lows better representation of real life phenomenon and enriches its dynamics. Due to presence ofdelay terms in the model, Delay differential equations (DDEs) are infinite dimentional and henceare difficult to analyse. Hence now a days, solving delay differential equations is an importantarea of research. Every DDE cannot be integrated analytically and hence there is a need to be de-pendant on numerical methods to solve DDEs. To develop efficient, stable and accurate numericalalgorithms is primarily important task of research.DDEs are receiving increasing importance in many areas of science and engineering like bi-ological processes, population growth and decay models, epidemiology, physiology, neural net-works etc [20], [6], [3]. The classical numerical methods such as Eulers method, trapezoidalmethod, Runge-Kutta methods are discussed in [8]. Enright and Hu [14] have developed a toolto solve DDEs with vanishing delay with iteration and interpolation technique. Karoui and Vail-lancourt [19] presneted a SYSDEL code to solve DDEs using Runge kutta methods of desiredconvergence order. In [21], a new MATLAB program dde23 has been developed to solve widerange of DDEs with constant delays. A new Adomian decomposition method is given in [15]to solve DDEs. Further in [17], two point predictor-corrector block method for solving DDEs isdescribed. Recently New iterative method (NIM) developed by Daftardar-Gejji and Jafari is usedto develop many efficient numerical methods to solve fractional differential equations [12], partialdifferential equations [9], boundary value problems [10]. In [2], a solution of DDE is acheivedby using Aboodh transformation method. Recently in [13] , a new numerical technique for solv-ing fractional generalised pantograph-delay differential equations by using fractional order hybridBessel functions is developed. 1 a r X i v : . [ m a t h . NA ] F e b n the present article, we use the powerful technique of NIM to generate new efficient numer-ical tools to solve DDEs which are reducible to solve ODEs.The paper is organized as follows. In section 2, important preliminaries like Delay differentialequations (DDEs), New iterative method (NIM) etc. are reviewed. In next section 3, we developeda new family of numerical methods to solve DDEs and system of DDEs. In section 4, error analysisof newly proposed methods is done along with its stability analysis in section 5. In section 6, someillustrative examples are solved to check the accuracy of new methods practically. In last section7, some important observations are made on the basis of theoretical stability analysis and erroranalysis of newly proposed methods. Consider a general form of an initial value problem (IVP) representing a time delay differentialequation: y (cid:48) = f ( t, y ( t ) , y ( t − τ )); t ∈ [0 , T ] , (1) y ( t ) = φ ( t ) , − τ ≤ t ≤ , (2)where φ ( t ) : [ − τ, → R n ; n ∈ N is a real valued function which represents histroy of thesolution in the past. Existence and uniqueness conditions for the solution of DDEs are given in[23]. For existence and uniqueness of the solution of the IVP (1-2), we assume that f is non-linear,bounded and continuous real valued operator defined on [0 , T ] × R n × C ( R , R n ) to R n and fulfilsLipschitz conditions with respect to the second and third arguments: | f ( t, x , u ) − f ( t, x , u ) |≤ L | x − x | , (3) | f ( t, x, u ) − f ( t, x, u ) |≤ L | u − u | , (4)where L and L are positive constants. y ( t − τ ) The approximation to the delay term y ( t n − τ ) is denoted by ν n [12].When τ is constant ( t n − τ ) may not be a grid point t n for any n. Suppose ( m + δ ) h = τ, m ∈ N and ≤ δ < . If δ = 0 , mh = τ and y ( t n − τ ) is approximated as y ( t n − τ ) ≈ ν n = (cid:26) y n − m if n ≥ m ; φ ( t n − τ ) if n < m. Daftardar-Gejji and Jafari [11] have proposed a new iterative method (NIM) for solving linear/non-linear functional equations of the form u = f + L ( u ) + N ( u ) (5)2here f is a known function, L is a linear operator and N a non linear operator. This decom-position technique is used by many researchers for solving variety of problems such as fractionaldifferential equations [12], boundary value problems [10], system of non-linear functional equa-tions [5] etc. Recently NIM is used to develop efficient numerical algorithms to solve ordinarydifferential equations and time delay fractional differential equations too [22],[18].In this method, We assume that eq.(5) has a series solution of the form u = ∞ (cid:88) i =0 u i . Since L isa linear operator, we have L (cid:32) ∞ (cid:88) i =0 u i (cid:33) = ∞ (cid:88) i =0 L ( u i ) and non linear operator N is decomposed as N ( u ) = N ( u ) + [ N ( u + u ) − N ( u )] + [ N ( u + u + u ) − N ( u + u )] + . . . Let G = N ( u ) and G i = N (cid:32) i (cid:88) n =0 u n (cid:33) − N (cid:32) i − (cid:88) n =0 u n (cid:33) , i = 1 , , , . . . .Observe that N ( u ) = ∞ (cid:88) i =0 G i .Putting series solution u = ∞ (cid:88) i =0 u i in eq.(5), we get ∞ (cid:88) i =0 u i = f + ∞ (cid:88) i =0 L ( u i ) + ∞ (cid:88) i =0 G i . Taking u = f , and u n = L ( u n − ) + G n − , n = 1 , , , . . . Note that u = u + u + u + · · · = f + L ( u )+ N ( u )+ L ( u )+[ N ( u + u ) − N ( u )]+ · · · = f + L ( u )+ N ( u ) . Hence u satisfies the functional eq.(5). k -term NIM solution is given by u = k − (cid:88) i =0 u i . In this section, we represent a family of new numerical methods based on NIM and implicit nu-merical methods to solve delay differential equations. In [22], we developed a new algorithm tosolve DDEs and ODEs by improving existing trapezoidal rule of integration. In present paper,we generalize this task by improving general θ − methods to solve delay differential equations asfollows.For solving eq.(1-2) on [0 , T ] , consider the uniform grid t n = nh, n = − m, − m +1 , ...., − , , , ....N, where m and N are integers such that N = T /h and m = τ /h. We integrate eq.(1) from the node t n to t n +1 on both sides, which gives us y ( t n +1 ) = y ( t n ) + (cid:90) t n +1 t n f ( t, y ( t ) , y ( t − τ )) dt. (6)3ow we approximate integration on R.H.S. in above equation as; (cid:90) t n +1 t n f ( t, y ( t ) , y ( t − τ )) dt ≈ h [(1 − θ ) f ( t n , y ( t n ) , ν n ) + θf ( t n +1 , y ( t n +1 ) , ν n +1 )] , (7)where θ is a parameter which lies in the closed interval [0 , and ν n is approximation to delayterm at the node t = t n . Applying the approximation from eq.(7) in eq.(6), we get y n +1 = y n + h [(1 − θ ) f ( t n , y n , ν n ) + θf ( t n +1 , y n +1 , ν n +1 )] , n = 0 , , . . . , N − . (8)This is called as a family of θ − methods to solve DDEs.In particular, when θ = 1 , eq.(8) reduces to y n +1 = y n + hf ( t n +1 , y n +1 , ν n +1 ) , n = 0 , , . . . , N − . (9)This eq.(9) is referred as Implicit Euler’s method. Note that when θ = 0 , eq.(8) yields ExplicitEuler’s method and when θ = 12 we get implicit Trapezoidal rule of integration to solve DDEs. To generate a new family of θ − methids for solving DDEs, we note that eq.(8) is of the form u = f + N ( u ) , and hence NIM can be employed as follows:We write, u = y n +1 ,f = y n + h (1 − θ ) f ( t n , y n , ν n ) ,N ( u ) = hθf ( t n +1 , y n +1 , ν n +1 ) . For simplicity we denote, f n = f ( t n , y n , ν n ) . Three term NIM solution of eq.(8) gives us u = u + u + u = u + N ( u ) + N ( u + u ) − N ( u )= u + N ( u + u )= u + N ( u + N ( u )) . That is y n +1 = y n + h (1 − θ ) f n + N ( u + u )= y n + h (1 − θ ) f n + hθf ( t n +1 , u + u , ν n +1 ) Therefore, y n +1 = y n + h (1 − θ ) f n + hθf ( t n +1 , y n + h (1 − θ ) f n + hθf ( t n +1 , y n + h (1 − θ ) f n , ν n +1 ) , ν n +1 ) (10)Eq.(10) represents new family of θ − methods to solve DDEs which can be expressed in the fol-lowing more simple form as follows: k = f ( t n , y n , ν n ) k = f ( t n +1 , y n + h (1 − θ ) k , ν n +1 ) k = f ( t n +1 , y n + h (1 − θ ) k + hθk , ν n +1 ) where, y n +1 = y n + h (1 − θ ) k + hθk . (11)4 ase 1: When θ = 1 , k = f ( t n , y n , ν n ) k = f ( t n +1 , y n , ν n +1 ) k = f ( t n +1 , y n + hk , ν n +1 ) where, y n +1 = y n + hk . (12)Eqs.(12) represents new improved implicit Eulers method to solve DDEs. Case 2:
When θ = , k = f n = f ( t n , y n , ν n ) k = f ( t n +1 , y n + hk , ν n +1 ) k = f ( t n +1 , y n + hk + hk , ν n +1 ) where, y n +1 = y n + hk + hk ., (13)Eqs.(13) represents improved trapezoidal rule to solve DDEs which is developed in [22]. Case 3:
When θ = 0 , k = f n = f ( t n , y n , ν n ) k = f ( t n +1 , y n + hk , ν n +1 ) k = f ( t n +1 , y n + hk , ν n +1 ) where, y n +1 = y n + hk ., (14)Eqs.(14) represents a method which is usual explicit Eulers method to solve DDEs. Case 4:
When θ = , k = f n = f ( t n , y n , ν n ) k = f ( t n +1 , y n + hk , ν n +1 ) k = f ( t n +1 , y n + hk + hk , ν n +1 ) where, y n +1 = y n + hk + hk ., (15)Eqs.(15) is a new method to solve DDEs.In above family of methods to solve DDEs, when delay term y ( t n − τ ) and its approximation ν n are zero then DDE (1) reduces to ODE without delay. In accordance with this, above family ofmethods gets reduced to the family of numerical methods for solving ODEs (without delay) whichare developed in [1]. General form of Runge-Kutta method is given by y n +1 = y n + h (cid:88) i =1 b i k i and k = f n = f ( t n , y n , ν n ) k = f ( t n + c h, y n + ha k , ν n ) k = f ( t n + c h, y n + h ( a k + a k ) , ν n ) (16)5n newly proposed θ − methods represented by eqs.(11), we have b ( θ ) = 1 − θ, b ( θ ) = 0 , b ( θ ) = θ, a = 1 − θ, a = 1 − θ, a = θ, c = 1 , c = 1 . Therefore new Family of θ − methods forsolving DDEs can be stated in the form of Butcher tableu as follows.01 1- θ θ θ θ θ For a Runge Kutta method it is necessary to satisfy that i − (cid:88) j =1 a ij = c i ∀ i = 2 , cf. [4]). Fromthe above tableau, for i = 2 , (cid:88) j =1 a j = c if and only if θ = 0 . This shows that, newly proposedfamily of numerical methods for solving DDEs is different from Runge Kutta methods except for θ = 0 . θ − methods for solving a system of delay differential equations The numerical algorithm presented in above section can be generalized for solving the followingsystem of DDEs: y (cid:48) ( t ) = f ( t, y ( t ) , y ( t − τ )) ,y (cid:48) ( t ) = f ( t, y ( t ) , y ( t − τ )) , ... y (cid:48) m ( t ) = f m ( t, y ( t ) , y ( t − τ )) , with initial condition y ( t ) = ( y ( t ) , y ( t ) , . . . , y m ( t )) = ( φ ( t ) , φ ( t ) , ..., φ m ( t )); − τ ≤ t ≤ . (17)We let y n be a vector of independent variables representing values of ( y , y , . . . , y m ) at node t = t n and ν n is a vector approximation to ( y ( t − τ ) , y ( t − τ ) , . . . , y m ( t − τ )) at t = t n . Weobtain a new family of θ − methods to solve a system of DDEs as follows: k ,y i = f i ( t n , y n , ν n ) k ,y i = f i ( t n +1 , y n + h (1 − θ ) k ,y i , ν n +1 ) k ,y i = f i ( t n +1 , y n + h (1 − θ ) k ,y i + hθk ,y i , ν n +1 ) Where, y n +1 = y n + h (1 − θ ) k ,y i + hθk ,y i , i = 1 , , . . . , m. Error analysis
Theorem 4.1
The new family of θ − methods given by eqs.(11) forms a second order numericalmethod for θ = and has a first order convergence for any other value of θ ∈ [0 , . Proof: Using Taylor’s series expansion, y ( t n +1 ) = y ( t n + h ) = y n + hf n + h f t + f f y + f ν ) + O ( h ) . Now consider θ − method and Taylors expansions applied in k , k .k = f ( t n , y n , ν n ) (18) k = f ( t n +1 , y n + h (1 − θ ) k , ν n +1 )= f ( t n + h, y n + h (1 − θ ) k , ν n + h ) k = f n + ( hf t + h (1 − θ ) k f y + hf ν ) +12 (cid:0) hf tt + h (1 − θ ) k f yy + h f νν + 2 h (1 − θ ) k f ty + 2 h f tν + 2 h (1 − θ ) k f yν (cid:1) + O ( h ) . (19) k = f ( t n +1 , y n + h (1 − θ ) k + hθk , ν n +1 )= f ( t n + h, y n + h (1 − θ ) k + hθk , ν n + h ) k = f n + ( hf t + h (1 − θ ) k f y + hθk f y + hf ν )+ 12 (cid:0) h f tt + ( h (1 − θ ) k + h θ k + 2 h (1 − θ ) θk k ) f yy + h f νν + (2 h (1 − θ ) k + 2 h θk ) f yt + 2 h f tν + 2( h (1 − θ ) k + h θk ) f yν (cid:1) + O ( h ) . (20)Putting (18), (19), (20) in θ − method given by eqs.(11), error e n +1 of the method is given by (21) e n +1 = y ( t n +1 ) − y n +1 = y n + hf n + h f t + h f f y h f ν O ( h ) − (cid:0) y n + h (1 − θ ) f n + hθf n h θf t + h θ (1 − θ ) k f y + h θ k f y + h θf ν (cid:1) + O ( h ) . Clearly, for θ = 0 , , / method has a linear convergence and for θ = 1 / its order of convegenceis quadratic. 7 .1 Error analysis of new improved implicit Euler’s method for solving DDEs Refer to Eqs.(12), For θ = 1 , we get new improved implicit Eulers method to solve DDEs, whichis given by the formula: y n +1 = y n + hf [ t n +1 , y n + hf ( t n +1 , y n , ν n +1 ) , ν n +1 ] . (22)By inserting analytical solution in the above equation, we get the truncation errorv T n as y ( t n +1 ) − y ( t n ) h − f [ t n +1 , y ( t n ) + hf ( t n +1 , y ( t n ) , ν n +1 ) , ν n +1 ] = T n . (23)and by eq.(23), we get y n +1 − y n h − f [ t n +1 , y n + hf ( t n +1 , y n , ν n +1 ) , ν n +1 ] = 0 . (24)Subtracting eq. (23) and eq. (24) , we get T n = e n +1 − e n h − f [ t n +1 , y ( t n )+ hf ( t n +1 , y ( t n ) , ν n +1 ) , ν n +1 ]+ f [ t n +1 , y n + hf ( t n +1 , y n , ν n +1 ) , ν n +1 ] . (25)This implies that, hT n = e n +1 − e n − h ( f [ t n +1 , y ( t n ) + hf ( t n +1 , y ( t n ) , ν n +1 ) , ν n +1 ] − f [ t n +1 , y n + hf ( t n +1 , y n , ν n +1 ) , ν n +1 ]) . Therefore, | e n +1 |≤ | e n | + hL | ( y ( t n ) + hf ( t n +1 , y ( t n ) , ν n +1 ) − ( y n + hf ( t n +1 , y n , ν n +1 ) | + hT n ≤ | e n | + hL | e n | + h L | e n | + | T n | h ≤ (1 + hL + h L ) | e n | + | T n | h ; n = 0 , , . . . N. Let T = max ≤ n ≤ ( N − | T n | Therefore, | e n +1 |≤ (1 + hL + h L ) | e n | + T h
Now by induction, | e n |≤ (1 + hL + h L ) n | e | + (cid:18) (1 + hL + h L ) n − − hL + h L ) − (cid:19) T h ≤ (1 + hL + h L ) n | e | + (cid:18) e nhL − Lh + 1) Lh (cid:19) T h ≤ e nhL | e | + (cid:18) e nhL − L (cid:19) T ≤ e ( t n − t ) L | e | + (cid:32) e ( t n − t ) L − L (cid:33) T nh = t n − t . Noting that, Truncation error T ≤ hy (cid:48)(cid:48) ( ζ ) and if y (cid:48)(cid:48) ≤ M then T ≤ hM . Therefore, | e n |≤ e ( t n − t ) L | e | + (cid:32) ( e ( t n − t ) L − hM L (cid:33) ≤ (cid:32) ( e ( t n − t ) L − hM L (cid:33) , noting that | e | is zero.hence as h → , | e n |→ . Here L is Lipschitz constant as given in eq.(3). Hence method isconvergent. Similary, for θ = , we proved the convergence of the new improved trapezoidal rulefor solving DDEs in [22] and for other values of θ convergenrce can be proved on similar lines. Definition 5.1
The numerical method for solving IVP eq.(1-2) is said to be zero-stable if smallperturbation in the initial condition of IVP do not cause the numerical approximation to divergefrom the exact solution, provided the exact solution of the IVP is bounded.
Consider the IVP eq.(1-2) and let (cid:15)y (0) = (cid:15)y be the new initial value (perturbed initialcondition) obtained by making a small change in y (0) = y . Theorem 5.1
Let y n be the solution obtained by new improved implicit Euler’s method (Case1 in section3) for solving DDE at the node t n with the initial condition y (0) = y and let (cid:15)y n be the solution obtained by the same numerical method with perturbed initial condition (cid:15)y = y + (cid:15) ; (cid:15) > . We assume that f ( t, y ) satisfies Lipschitz condition with respect to second variableand third variable with Lipschitz constant L , L then ∃ positive constants k and (cid:15) such that | y n − (cid:15)y n |≤ k(cid:15), ∀ nh ≤ T, h ∈ (0 , (cid:15) ) whenever | (cid:15) |≤ (cid:15). Proof: We prove the result by induction. | y n − (cid:15)y n | = | y n − + hf ( t n , y n − + hf ( t n , y n − , ν n ) , ν n ) − (cid:15)y n − − hf ( t n , (cid:15)y n − + hf ( t n , (cid:15)y n − , ν n ) , ν n ) |≤ | y n − − (cid:15)y n − | + hL | y n − + hf ( t n , y n − , ν n ) − (cid:15)y n − − hf ( t n , (cid:15)y n − , ν n ) |≤ | y n − − (cid:15)y n − | + hL | y n − − (cid:15)y n − | + h L | y n − − (cid:15)y n − | = (1 + hL + h L ) | y n − − (cid:15)y n − | Therefore by induction, | y n − (cid:15)y n |≤ (1 + hL + h L ) n | y − (cid:15)y |≤ e nhL | (cid:15) | = e T M L (cid:15) = k(cid:15), where k = e T M L > . This proves that the new method for ( θ = 1 ) is stable.9 Illustrative examples
To demonstrate applicability of newly proposed methods, We present here some illustrative exam-ples which are solved using Mathematica 12.
Example 6.1
Consider the delay logistic differential equation y (cid:48) ( t ) = 0 . y ( t )(1 − y ( t − , y ( t ≤
0) = 0 . . (26) (a) Solution of (26) by new method ( θ = 1 ) (b) Error in solution of (26) when solved using newmethod ( θ = 1 ) Figure 1: Dashed graph: Solution by new method ( θ = 1 ), Dotted graph: Exact solution. It isnoted that in Fig.1 (a) exact solution and approximate solution by the new method ( θ = 1 ) overlapson each other.Step length is taken as h = 0 . . Example 6.2
Consider the differential equation without delay [7] y (cid:48) ( t ) = 2 − e − t − y, y (0) = 1 , ≤ t ≤ . (27)Exact solution of the differential equation (27) is y ( t ) = 1 + e − t − e − t . (a) Solution of (27) by new method ( θ = 1 ) (b) Error in solution of (27) when solved using θ method ( θ = 1 ) Figure 2: Dashed graph: Solution by new method, Dotted graph: Exact solution. It is noted thatin Fig.2 (a) exact solution and approximate solution by the new method overlaps on each other.Step length is taken as h = 0 . . E and error in new method ( θ = 1 ) E (obtainedby taking 3-term NIM solution) while solving (27) are compared for various values of t . It isnoteworthy that E is always smaller than E in all cases. So the new method is more accuratethan implicit backward Euler’s method. t n S S S S e e e .
01 0 . . . . . . . .
02 0 . . . . . . . .
03 0 . . . . . . . .
04 0 . . . . . . . .
05 0 . . . . . . . .
06 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 1: Ex.(8.1) S : Solution by backward Euler’s method S : Solution by new method ( θ = 1 ) obtained by taking 3-term NIM solution S : Solution by new method ( θ = 1 ) obtained by taking 4-term NIM solution S : Exact solution e : Error in solution by backward Euler’s method e : Error in solution by new method ( θ = 1 ) obtained by taking 3-term NIM solution e : Error in solution by new method ( θ = 1 ) obtained by taking 4-term NIM solution Observation:
In this example, it is observed that (4-term NIM solution) new method ( θ = 1 )method and backward Euler’s method gives same error and error in these two methods is greaterthan (3-term NIM solution) new method ( θ = 1 ). Hence, new method with three term NIMsolution gives better accuracy than implicit backward euler method and new method with 4-termNIM solution. ¨ o ssler System with delay Consider the R ¨ ossler system [16] with delay given by the following system of differential equa-tions. ˙ x = − y ( t ) − z ( t ) , ˙ y = x ( t ) + ay ( t − , ˙ z = b + z ( t )( x ( t ) − c ) . (28)11et a = b = 0 . , x (0) = y (0) = z (0) = 0 . . Note that c is a control parameter. Solvingabove system by new method for ODE we obtain the x-waveforms and x-y phase potraits whichare depicted in Figs.3(a)-(b), Figs.4(a)-(b) and in Figs.5(a)-(b) for c = 2 . , c = 2 . , c = 7 . respectively. (a) x-waveform of the system (28) (b) xy-phase diagram of the system (28) Figure 3: c=2.3 in system(28) (a) x-waveform of the system (28) (b) xy-phase diagram of the system (28)
Figure 4: c=2.9 in system(28)12 a) x-waveform of the system (28) (b) xy-phase diagram of the system (28)
Figure 5: c=7.9 in system(28)
In the present work, a new family of methods for solving delay differential equations (DDEs) hasbeen proposed which are reducible to solve ordinary differential equations too (without delay).Newly proposed methods are then compared with existing methods with respect to stability andaccuracy. New methods formed found to be stable and accurate. Further error analysis and sta-bility analysis of new methods is carried out. Numerous illustrative examples are solved usingMathematica 12 to demonstrate the efficiency of the method. It is observed that, new methodsare non-Runge Kutta methods and are more accurate than existing numerical methods for solvingDDEs.
References [1] O. Y. Ababneh. New numerical methods for solving differential equations, 2019.[2] K. Aboodh, R. Farah, I. Almardy, and A. Osman. Solving delay differential equations byaboodh transformation method.
International Journal of Applied Mathematics & StatisticalSciences , 7(2):55–64, 2018.[3] C. Baker, G. Bocharov, C. Paul, and F. Rihan. Modelling and analysis of time-lags in somebasic patterns of cell proliferation.
Journal of mathematical biology , 37(4):341–371, 1998.[4] A. Bellen and M. Zennaro.
Numerical methods for delay differential equations . OxfordUniversity Press, 2013. 135] S. Bhalekar and V. Daftardar-Gejji. Solving a system of nonlinear functional equations usingrevised new iterative method.
International Journal of Mathematical and ComputationalSciences , 6(8):968–972, 2012.[6] G. A. Bocharov and F. A. Rihan. Numerical modelling in biosciences using delay differentialequations.
Journal of Computational and Applied Mathematics , 125(1-2):183–199, 2000.[7] T. Bui. Explicit and implicit methods in solving differential equations. 2010.[8] J. C. Butcher. Numerical methods for ordinary differential equations in the 20th century.
Journal of Computational and Applied Mathematics , 125(1-2):1–29, 2000.[9] V. Daftardar-Gejji and S. Bhalekar. Solving fractional diffusion-wave equations using a newiterative method.
Fractional Calculus and Applied Analysis , 11(2):193p–202p, 2008.[10] V. Daftardar-Gejji and S. Bhalekar. Solving fractional boundary value problems with dirich-let boundary conditions using a new iterative method.
Computers & Mathematics with Ap-plications , 59(5):1801–1809, 2010.[11] V. Daftardar-Gejji and H. Jafari. An iterative method for solving nonlinear functional equa-tions.
Journal of Mathematical Analysis and Applications , 316(2):753–763, 2006.[12] V. Daftardar-Gejji, Y. Sukale, and S. Bhalekar. Solving fractional delay differential equa-tions: A new approach.
Fractional Calculus and Applied Analysis , 18(2):400–418, 2015.[13] H. Dehestani, Y. Ordokhani, and M. Razzaghi. Numerical technique for solving fractionalgeneralized pantograph-delay differential equations by using fractional-order hybrid besselfunctions.
International Journal of Applied and Computational Mathematics , 6(1):1–27,2020.[14] W. H. Enright and M. Hu. Interpolating runge-kutta methods for vanishing delay differentialequations.
Computing , 55(3):223–236, 1995.[15] D. J. Evans and K. Raslan. The adomian decomposition method for solving delay differentialequation.
International Journal of Computer Mathematics , 82(1):49–54, 2005.[16] K. Ibrahim, R. Jamal, and F. Ali. Chaotic behaviour of the rossler model and its analysis byusing bifurcations of limit cycles and chaotic attractors. In
J. Phys. Conf. Ser , volume 1003,page 012099, 2018.[17] F. Ishak, M. Suleiman, and Z. Omar. Two-point predictor-corrector block method for solvingdelay differential equations.
MATEMATIKA: Malaysian Journal of Industrial and AppliedMathematics , 24:131–140, 2008.[18] A. Jhinga and V. Daftardar-Gejji. A new numerical method for solving fractional delaydifferential equations.
Computational and Applied Mathematics , 38(4):1–18, 2019.[19] A. Karoui and R. Vaillancourt. A numerical method for vanishing-lag delay differentialequations.
Applied Numerical Mathematics , 17(4):383–395, 1995.1420] F. A. Rihan, D. Abdelrahman, F. Al-Maskari, F. Ibrahim, and M. A. Abdeen. Delay differ-ential model for tumour-immune response with chemoimmunotherapy and optimal control.
Computational and mathematical methods in medicine , 2014, 2014.[21] L. F. Shampine, S. Thompson, and J. Kierzenka. Solving delay differential equations withdde23. , 2000.[22] Y. Sukale and V. Daftardar-Gejji. New numerical methods for solving differential equations.
International Journal of Applied and Computational Mathematics , 3(3):1639–1660, 2017.[23] E. Süli and D. F. Mayers.