Solving Second-Order Delay Differential Equations by Direct Adams-Moulton Method

Abstract

This paper will consider the implementation of fifth-order direct method in the form of Adams-Moulton method for solving directly second-order delay differential equations (DDEs). The proposed direct method approximates the solutions using constant step size. The delay differential equations will be treated in their original forms without being reduced to systems of first-order ordinary differential equations (ODEs). Numerical results are presented to show that the proposed direct method is suitable for solving second-order delay differential equations.

1. Introduction

In the recent years, there are rigorous and numerous researches undertaken in the areas of science and engineering that are skewed towards the developments of the mathematical models involving the delay differential equations. In mathematics, the DDEs are differing from ODEs in which the derivative of the unknown function at a certain time is given in terms of the values of the function at previous times. Several numerical methods have been proposed to solve first-order DDEs such as in [1–4]. However, less attention was made to solve second-order DDEs. Spline collocation methods [5] and Adomian decomposition method [6] have been proposed to solve second-order DDEs directly. There were few numerical methods that have been proposed for solving ODEs and those methods have been extended to solve DDEs with some modifications in the algorithm. These are some works carried out for solving first-order DDEs using the extended version such as in San et al. [1], Ismail et al. [2], Radzi et al. [3], and Ishak et al. [4].

In this paper, we are concerned with solving second-order delay differential equations (DDEs) as follows: 𝑦^′′(𝑡)(),𝑦=𝑓𝑡,𝑦,𝑦𝑡−𝜏𝑎≤𝑡≤𝑏,𝜏>0,_′(𝑎)𝑦(𝑡)(𝑡)||,=Ω,=𝜙,𝛼≤𝑡≤𝑎,0≤𝜏≤𝑎−𝛼 (1)(2)(3) where 𝜙(𝑡) is the initial function and 𝜏 is the delay term.

The direct Adams-Moulton methods were studied by several researchers and the methods have shown their ability to solve first-, second- and higher-order ODEs [7, 8] and boundary value problems [9] effectively and accurately. Hence, in this paper, we aim to propose the direct method of order five in the form of Adams-Moulton method for solving (1) using constant step size. The proposed direct method has the advantage to solve the second-order delay differential equations directly without reducing to system of first order. Therefore, the second-order DDEs problems will be handled in their original forms. This idea will reduce the computation cost at each step for the proposed direct method.

2. The Direct Method

2.1. Formulation of the Method

Most numerical methods for ODEs can be used to solve DDEs. In this paper, we adopted direct method proposed by Majid and Suleiman [8] to solve DDEs. The following is the derivation of one point direct Adams-Moulton method.

The point 𝑦_𝑛+1 at 𝑥_𝑛+1 can be obtained by integrating 𝑦^′′=𝑓𝑡,𝑦(𝑡),𝑦^′ (4) once and twice as follows.

Integrating once ^𝑥𝑛+1_𝑥𝑛𝑦^′′(𝑥)𝑑𝑥=^𝑥𝑛+1_𝑥𝑛𝑓𝑥,𝑦,𝑦^′𝑦𝑑𝑥,_′𝑥_𝑛+1=𝑦^′𝑥_𝑛+^𝑥𝑛+1_𝑥𝑛𝑓𝑥,𝑦,𝑦^′𝑑𝑥. (5)(6)

Integrating twice ^𝑥𝑛+1_𝑥𝑛^𝑥_𝑥𝑛𝑦^′′(𝑥)𝑑𝑥𝑑𝑥=^𝑥𝑛+1_𝑥𝑛^𝑥_𝑥𝑛𝑓𝑥,𝑦,𝑦^′𝑦𝑥𝑑𝑥𝑑𝑥,_𝑛+1−𝑦𝑥_𝑛−ℎ𝑦^′𝑥_𝑛=^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑥𝑓𝑥,𝑦,𝑦^′𝑑𝑥. (7)(8)

The function 𝑓(𝑥,𝑦,𝑦^′) in (6) and (8) will be approximated using Lagrange interpolation polynomial and the interpolation points involved are five points, that is, {𝑥_𝑛−3,𝑥_𝑛−2,𝑥_𝑛−1,𝑥_𝑛,𝑥_𝑛+1}. Taking 𝑠=(𝑥−𝑥_𝑛+1)/ℎ and replacing 𝑑𝑥=ℎ𝑑𝑠, the value of 𝑦_𝑛+1 can be obtained by integrating (6) and (8) using MAPLE. The direct method is the combination of predictor of one order less than the corrector. The following is the predictor and corrector of the direct method.

Predictor 𝑦^′𝑥_𝑛+1=𝑦^′𝑥_𝑛+ℎ2455𝑓_𝑛−59𝑓_𝑛−1+37𝑓_𝑛−2−9𝑓_𝑛−3,𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^′𝑥_𝑛+ℎ²360323𝑓_𝑛−264𝑓_𝑛−1+159𝑓_𝑛−2−38𝑓_𝑛−3. (9)

Corrector 𝑦^′𝑥_𝑛+1=𝑦^′𝑥_𝑛−ℎ×720−251𝑓_𝑛+1−646𝑓_𝑛+264𝑓_𝑛−1−106𝑓_𝑛−2+19𝑓_𝑛−3,𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^′𝑥_𝑛−ℎ²×1440−135𝑓_𝑛+1−752𝑓_𝑛+246𝑓_𝑛−1−96𝑓_𝑛−2+17𝑓_𝑛−3. (10)(11)

2.2. Order and Error Constant of the Method

The order of this developed method is calculated by referring to [10–14]. Linear 𝑘-step method can be written in the form of 𝛼𝑌_𝑚=ℎ𝛽𝑌^′_𝑚+ℎ²𝛾𝐹_𝑚, (12) where 𝛼, 𝛽, and 𝛾 are the coefficients with the 𝑚-vector. 𝑌_𝑚, 𝑌^′_𝑚 and 𝐹_𝑚 are 𝑌_𝑚=𝑦_𝑛−1,𝑦_𝑛^𝑇,𝑌_𝑚=𝑦^′_𝑛−1,𝑦^′_𝑛^𝑇,𝐹_𝑚=𝑓_𝑛−1,𝑓_𝑛^𝑇. (13)

The formulae is defined as 𝐶₀=𝛼₀+𝛼₁+⋯+𝛼_𝑘𝐶₁=𝛼₁+2𝛼₂+⋯+𝑘𝛼_𝑘−𝛽₀+𝛽₁+⋯+𝛽_𝑘𝐶₂=1𝛼2!¹+2²𝛼₂+⋯+𝑘²𝛼_𝑘−𝛽₁+2𝛽₂+⋯+𝑘𝛽_𝑘−𝛾₀+𝛾₁+⋯+𝛾_𝑘⋮𝐶_𝑞=^𝑘_𝑗=0𝑗^𝑞𝛼𝑞!^𝑗−𝑗^𝑞−1!𝛽𝑞−1^𝑗−𝑗^𝑞−2!𝛾𝑞−2^𝑗,where𝑞=3,4,5,…. (14)

The method has order 𝑝 if  𝐶₀=𝐶₁=⋯=𝐶_𝑝=𝐶_𝑝+1=0 and 𝐶_𝑝+2≠0 is the error constant.

By referring to the corrector formulae in (11) and the general multistep method in (12), we will obtain 𝛼₀=0,𝛼₁𝛼=0,₂=−1,𝛼₃=0,𝛼₄𝛽=1;₀=0,𝛽₁𝛽=0,₂=0,𝛽₃=1,𝛽₄𝛾=0;₀=−171440,𝛾₁=96,𝛾1440₂=−2461440,𝛾₃=7521440,𝛾₃=1351440 (15) and substitute those values into (14): 𝐶₀=𝛼₀+𝛼₁+𝛼₂+𝛼₃+𝛼₄𝐶=0,₁=𝛼₁+2𝛼₂+3𝛼₃+4𝛼₄−𝛽₀+𝛽₁+𝛽₂+𝛽₃+𝛽₄𝐶=0,₂=1𝛼2!¹+2²𝛼₂+3²𝛼₃+4²𝛼₄−𝛽₁+2𝛽₂+3𝛽₃+4𝛽₄−𝛾₀+𝛾₁+𝛾₂+𝛾₃+𝛾₄𝐶=0,₃=1𝛼3!¹+2²𝛼₂+3²𝛼₃+4²𝛼₄−1𝛽2!¹+2²𝛽₂+3²𝛽₃+4²𝛽₄−𝛾₁+2𝛾₂+3𝛾₃+4𝛾₄𝐶=0,₄𝐶=0,₅𝐶=0,₆𝐶=0,₇=−415040≠0. (16)

Hence, 𝐶₀=𝐶₁=⋯=𝐶₆=0 and 𝐶₇=−(41/5040)≠0 is the error constant.

The corrector of the direct method is of order five and the error constant is −(41/5040). The method is said to be consistent if it has at least one order. Since the proposed method is order five, hence the method is said to be consistent.

2.3. Stability Analysis

The method is zero stable provided the roots 𝑅_𝑗 of the first characteristic polynomial 𝜌(𝑅) specified as 𝜌(𝑅)=det[∑^𝑘_𝑖=0𝐴^(𝑖)𝑅^𝑘−𝑖]=0 and satisfy |𝑅_𝑗|≤1.

We rewrite (10) and (11) in the matrix form: 10𝑦01𝑦^′𝑛+1_𝑛+1=10𝑦01𝑦^′𝑛_𝑛𝑦+ℎ0010𝑦^′𝑛+1_𝑛+10+ℎ251𝑓72000𝑓^′𝑛+1_𝑛+10+ℎ646𝑓72000𝑓^′𝑛_𝑛+ℎ0−264𝑓72000𝑓^′𝑛−1_𝑛−10+ℎ106𝑓72000𝑓^′𝑛−2_𝑛−2+ℎ0−19𝑓72000𝑓^′𝑛−3_𝑛−3+ℎ²000135𝑓1440^′𝑛+1𝑓_𝑛+1+ℎ²000752𝑓1440^′𝑛𝑓_𝑛+ℎ²000−246𝑓1440^′𝑛−1𝑓_𝑛−1+ℎ²00096𝑓1440^′𝑛−2𝑓_𝑛−2+ℎ²000−17𝑓1440^′𝑛−3𝑓_𝑛−3. (17)

The first characteristic polynomial of the method is given as

𝜌(𝑅)=det[𝑅𝐴⁰−𝐴¹]=0, where 𝐴⁰=¹⁰₀₁ and 𝐴¹=¹⁰₀₁. Consider that 𝜌(𝑅)()=det𝑅−100𝑅−1=0,𝑅−1²=0,𝑅=1,1. (18)

Since |𝑅_𝑗|≤1, the method is said to be zero stable.

There are many concepts of stability for numerical methods when applied to DDEs, depending on the test equation as well as the delay term involved. We would like to study the stability of the method by substituting the following test equation: 𝑓=𝜆𝑦(𝑥)+𝜇𝑦(𝑥−𝜏) (19) into the proposed method (10) and (11). The method can be described in the following matrix form: 𝐴₀𝑌_𝑁+1=𝐴₁𝑌_𝑁+ℎ⁴_𝑖=0𝐵_𝑖+1𝐹_𝑁−𝑖+ℎ²⁴_𝑖=0𝐶_𝑖+1𝐹_𝑁−𝑖, (20) where 𝑌_𝑁+1=𝑦^′_𝑛+1𝑦_𝑛+1,𝑌_𝑁=𝑦^′_𝑛𝑦_𝑛,𝐹_𝑁+𝑖=𝑦^′_𝑛+𝑖𝑦_𝑛+𝑖,𝐴₀=1−2513720ℎ𝜆01−ℎ32²𝜆,𝐴₁=,𝐵1001₁=0323360ℎ𝜆+251,𝐵720ℎ𝜇10₂=0−1130ℎ𝜆+323,𝐵360ℎ𝜇00₃=053360ℎ𝜆−11,𝐵30ℎ𝜇00₄=0−19720ℎ𝜆+53,𝐵360ℎ𝜇00₅=0−19720ℎ𝜇00,𝐶₁=00047390ℎ𝜆+,𝐶32ℎ𝜇₂=000−1160ℎ𝜆+4790ℎ𝜇,𝐶₃=010015ℎ𝜆−11,𝐶60ℎ𝜇₄=000−1711440ℎ𝜆+15ℎ𝜇,𝐶₅=000−17.1440ℎ𝜇 (21)

We solve the determination of 𝑑=0, where 𝑑=𝑡⁵𝐴₀−𝑡⁴𝐴₁+ℎ𝐵₁+ℎ²𝐶₁−𝑡³ℎ𝐵₂+ℎ²𝐶₂−𝑡²ℎ𝐵₃+ℎ²𝐶₃−𝑡ℎ𝐵₄+ℎ²𝐶₁−ℎ𝐵₅+ℎ²𝐶₅. (22)

The following stability polynomial is obtained by letting 𝐻1=ℎ²𝜇 and 𝐻2=ℎ²𝜆: 31−𝑡32𝐻2¹⁰+3−2+32𝐻1−373𝑡480𝐻2⁹+1−373480𝐻1−23𝑡120𝐻2⁸+−237120𝐻1+𝑡60𝐻2⁷+760𝐻1−11𝑡160𝐻2⁶+−117160𝐻1+𝑡480𝐻2⁵+7480𝐻1=0. (23)

The boundary of the stability region in 𝐻1-𝐻2 plane is determined by substituting the values of 𝑡=1,−1, and 𝑒^𝑖𝜃 into the stability polynomial, where 0≤𝜃≤2𝜋. Figure 1 shows the stability region of the direct method and the stable region is the bounded shaded region.

Figure 1 

Stability region of the direct method.

2.4. Convergence Analysis

A basic property for an effective numerical method is that the method needs to converge. A linear multistep method is convergent if and only if it is stable and consistent [15].

By definition, a linear multistep method of the form ^𝑘_𝑗=0𝛼_𝑖𝑦_𝑛+𝑗=ℎ^2𝑘_𝑗=0𝛽_𝑖𝑓_𝑛+𝑗 (24) is said to be consistent if the LMM is of order 𝑝≥1.

The proposed direct method is of order five, where 𝑝=5 and has error constant −(3/160)−(41/5040).

Since the method is order five which is ≥1, therefore, the method is consistent according to the definition. In the previous section, it has been shown that the method is zero stable. Therefore, we can conclude that this method is convergent.

3. Implementation

In the code of PECE scheme 𝑃 stands for an application of a predictor, 𝐸 stands for an evaluation of a function 𝑓, and 𝐶 stands for an application of a corrector. The developed code starts by using Adams-Bashforth method once at the beginning to calculate the three starting initial points. Once the points are calculated, then the proposed method can be applied until the end of the interval. The values of the delay term will be stored for future use. For 𝑥_𝑛−𝜏≤0, the delay term is calculated using the initial function given, 𝜙(𝑥). Otherwise the delay term depends on the location of (𝑥−𝜏). From the location we are able to recall the value which we had stored earlier. In this paper, no interpolation is required due to the implementation of constant step sizes. The algorithms of the proposed method were developed in 𝐶 language.

In the code, the selection of step sizes is predetermined.

3.1. Algorithm of Direct Method

Step 1. Set starting value 𝑎, ending value 𝑏, and step size ℎ, given initial value and given initial function 𝑔(𝑥).

Step 2. For 𝑛=1,2,3, set 𝑥_𝑛+1=𝑎+𝑛ℎ, compute function 𝑓_𝑛, and delay term 𝑑_𝑛.
Evaluate the approximate value 𝑦_𝑛+1 with direct Adams-Bashforth method.

Step 3. While 𝑥_𝑛<𝑏, do Step 4.

Step 4. Set 𝑥_𝑛+1=𝑥_𝑛+ℎ, compute function 𝑓_𝑛, and delay term 𝑑_𝑛.
Evaluate the approximate value 𝑦_𝑛+1 with proposed method.
Computing 𝑦^′_𝑛+1 and 𝑦_𝑛+1, we use the predictor formulas as follows: 𝑦^′𝑥_𝑛+1=𝑦^′𝑥_𝑛+ℎ2455𝑓_𝑛−59𝑓_𝑛−1+37𝑓_𝑛−2−9𝑓_𝑛−3,𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^′𝑥_𝑛+ℎ²360323𝑓_𝑛−264𝑓_𝑛−1+159𝑓_𝑛−2−38𝑓_𝑛−3. (25)

Computing 𝑦^′_𝑛+1 and 𝑦_𝑛+1, we use the corrector formulas as follows: 𝑦^′𝑥_𝑛+1=𝑦^′𝑥_𝑛−ℎ×720−251𝑓_𝑛+1−646𝑓_𝑛+264𝑓_𝑛−1−106𝑓_𝑛−2+19𝑓_𝑛−3,𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^′𝑥_𝑛−ℎ²×1440−135𝑓_𝑛+1−752𝑓_𝑛+246𝑓_𝑛−1−96𝑓_𝑛−2+17𝑓_𝑛−3. (26)

Step 5. Complete.

4. Numerical Result and Discussion

In order to study the efficiency of the proposed direct method, we presented three second-order DDE problems with constant delay in the following test problems. The numerical results of the direct method when solving Problems 1, 2, and 3 will be compared with cubic spline in [5], variable multistep method (VMM) in [16], and dde23 in MATLAB solver, respectively.

Problem 1. Consider 𝑦^′′(𝑡)1=−2𝑦(𝑡)+12𝑦()[],𝑦(𝑡)(𝑡)𝑡−𝜋,𝑡∈0,𝜋=1−sin,−𝜋≤𝑡≤0. (27)

Exact solution 𝑦(𝑡)=1−sin(𝑡). (28)

Problem 2. Consider 𝑦^′′(𝑡)(𝑡)()[],𝑦(𝑡)+𝑦=𝑦𝑡−1,𝑡∈0,1=𝑡²𝑦+3𝑡+2,−1≤𝑡≤0,_′(0)=0. (29)

Exact solution 𝑦(𝑡)=4cos(𝑡)−sin(𝑡)+𝑡²+𝑡−2. (30)

Problem 3. Consider 𝑦^′′(𝑡)()[],𝑦(𝑡)(𝑡)=𝑦𝑡−𝜋,𝑡∈0,𝜋=sin,−𝜋≤𝑡≤0. (31)

Exact solution 𝑦(𝑡)=sin(𝑡). (32)

The algorithm of the C language was executed on the Microsoft Visual C++ environment. The notations at the end of the paper are used in Tables 1–3.

Table 1 

Comparison of the numerical results for solving Problem 1.

Table 2 

Comparison of the numerical results for solving Problem 2.

Table 3 

Comparison of the numerical results for solving Problem 3.

The numerical results for solving Problems 1–3 are displayed in Tables 1–3.

In Problem 1, we solved the DDEs by the direct method and compare our results with the cubic spline method in [1]. Table 1 shows that the direct method managed to obtain highly accurate results compared to the cubic spline method at the same values of ℎ. At larger step size, that is, ℎ=(𝜋/10), we observed that the maximum error obtained by the direct method and the cubic spline is 6.12𝐸−06and 1.84𝐸−02, respectively. We could also observe that the maximum error for the direct method and the cubic spline was comparable, that is, 6.12𝐸−06 and 4.52𝐸−06 at ℎ=(𝜋/10) and 𝜋/(10∗2⁶), respectively. Hence, these results show that the direct method is able to obtain comparable maximum error compared to cubic spline at larger step size and therefore less expensive. The total function calls for the direct method during the computation are also shown in Table 1. Figure 2 display the comparison of the maximum error at different values of step sizes for solving Problem 1.

Figure 2 

Comparison of maximum error between cubic spline and direct method for solving Problem 1.

Table 2 displays the results for the direct method with the range of step sizes from 0.0128 to 0.1 compared to the variable multistep size method in [10] with varies step sizes (range from 0.0016 to 0.002) at different values of 𝑡 when solving Problem 2. The direct method is clearly superior compared to variable multistep method since it is able to obtain comparable and better results at larger step sizes. Hence, the direct method has less computational cost. The total function calls are also shown in Table 2. Figure 3 shows the comparison of maximum error at the given value of 𝑡.

Figure 3 

Comparison of maximum error between variable multistep method and direct method for solving Problem 2.

Table 3 displays the numerical results for solving Problem 3 for the MATLAB solver dde23 and the present direct method. Since dde23 is a variable step size method, therefore, the user may control the efficiency and accuracy of the solutions by changing the real tolerance and the absolute tolerance. The table will only show the compatible maximum absolute error with the least function calls. We could observe that the direct method managed to obtain comparable maximum errors with lesser total function calls and total step and to have no failure step compared to dde23. Hence, the direct method has less computational cost compared to dde23. Figure 4 shows the comparison results for maximum error and total function calls.

Figure 4 

Comparison of maximum error and function calls between dde23 and direct method for solving Problem 3.

5. Conclusion

In this study, we have shown that the proposed direct Adams-Moulton method using constant step size is suitable for solving second-order DDEs directly. The proposed direct method has solved the second-order DDEs in their original forms without being reduced to first-order ODEs. This approach has given advantage in terms of computational cost to the direct method. The method has shown superiority in terms of accuracy and it has less computational cost.

Notations

ℎ:	Step size
MAXE:	Magnitude of the maximum absolute error
FCN:	Total function calls
TS:	Total steps
TFS:	Total failure steps
RT:	Real tolerance
AT:	Absolute tolerance
dde23:	MATLAB solver dde23 based on the explicit Runge-Kutta (2,3) pair [17].

Acknowledgment

The authors gratefully acknowledged the financial support of University Putra Malaysia Grant.