Abstract
The search for efficient higher order methods is a constant goal in numerical analysis. In this paper, a higher order two-step hybrid block method is presented to directly solve second-order initial value problems in ordinary differential equations. In addition to the higher order, the proposed method has been formulated in variable step-size mode to extract its best performance. Comparisons with other methods in the literature show the good accuracy it can provide. Theoretical aspects such as linear stability and convergence analysis are also discussed.
1. Introduction
In literature, the numerical solution of the general second-order initial-value problem (IVP) of the formhas been on the rise. This is because it models many physical applied problems [1]. Sometimes Equation (1) can be transformed into a system of first-order ordinary differential equations (ODEs) and thus, solving it to obtain its numerical solution. A particular drawback of the latter is the high cost of CPU time due to the greater number of function evaluations [2]. The advantage of solving Equation (1) directly without considering transforming it resides in the fact that it achieves efficiency in terms of accuracy and less CPU time [3].
Few of these numerical approaches that are applied directly to solve Equation (1) include but are not limited to: Runge–Kutta methods, differential transform methods (DTM), linear multistep methods (LLMs), to mention but a few. In recent times, linear multistep block methods, credited to Milne [4], have been widely applied to solve Equation (1) directly. These methods have great advantages since they overcome the intersections of pieces of solutions and do not require any starting values provided by other methods, that is, they are self starting. For recent methods solving Equation (1) and higher order equations using linear multistep block methods, see [1, 3, 5–10]. It is worthy to mention that most of the approaches found used a constant step-size . This approach may perform poorly, especially if there are rapid and slow changes of the solution over the interval of integration. Efficient codes for solving IVPs is meant to automatically select the suitable step-size to achieve efficiency [11].
The two approaches to achieving this include: applying a scheme such that its coefficients rely on the ratios of the step sizes, or the use of a second technique to provide estimates of the local errors. The aim is to vary the step sizes so that one retains local errors smaller than a given tolerance and concurrently solving the underlying problem as efficiently as possible.
This paper aims to achieve a higher order variable step-size block hybrid integrator for solving Equation (1) directly. Though, a higher order method does not guarantee a better efficiency in terms of smaller global errors, as can be seen in comparing a method of seventh order by Jator [9], which performs better than a method of order eight by Tsitouras [12]. However, the higher order method derived here has more advantages which include, a small number of integration subintervals, smaller number of function evaluations, less costly in terms of CPU time, and small global error in terms of accuracy and its wide robust application for solving general second order IVPs. It is worth noting that the method implemented in its block form is self-starting. This means that, it does not require any external method to obtain the starting values. The novelty of the derived method is hinged on the advantages aforementioned among others.
2. Fundamentals
Definition 1. A linear -step method for solving Equation (1) is usually written in the formwhere , .
If the method is of implicit type, otherwise it is explicit. The coefficients are usually normalized assuming that .
Remark 1. The -step LMM Equation (2) is called linear because it involves only linear combinations of the and the .
Definition 2. The k-step method has first and second characteristic polynomials of the form
Definition 3. Given a continuously differentiable function , we may associate a linear difference operator to the LMM Equation (2) given by
Expanding Equation (4) in Taylor series about the point , the following is obtainedwhere are linear combinations of the coefficients .
Definition 4. The LMM Equation (2) is said to be order if , and in whichIn this case, is known as the principal error constant.
Definition 5. The difference operator is said to be consistent of order if it satisfies Equation (6) with for every sufficiently differentiable function .
Remark 2. An LMM whose associated difference operator is consistent of order is said to be consistent.
Definition 6. It is said that a polynomial satisfies the root condition if all its roots lie within or on the boundary of a unit circle, with those on the boundary being simple. That is, all roots must satisfy , and those whose modulus is unity must be simple.
Remark 3. If is a simple root of this is equivalent to say that is a factor of with multiplicity one.
Definition 7. A LMM is called be zero-stable if its first characteristic polynomial verifies the root condition.
Definition 8. The LMM Equation (2) is said to be convergent if, for all IVPs in Equation (1) having a unique solution ,holds for all.
For details, refer to the study by Lambert [13].
3. Derivation of the Method
Consider the uniform meshwith . Assume that the solution of Equation (1) is approximated by the following polynomial , that is,whose second and third derivatives are approximated bywhere are to be determined. We then consider the points , , for approximating the solution in the two-step interval . Using Equation (9) and its first derivative at the point , and the second and third derivatives in Equation (10), respectively, applied to the points , , the following system of 16 equations is obtained ; ; ; ; with , , , and .
Using the computer algebra system Mathematica, the above system is easily solved and the values of the coefficients , are obtained. These expressions are cumbersome and are not included here. After some algebraic simplification, we obtainedwhere are coefficients that depend on . By evaluating Equation (11) at the points , , we obtain the following main formulas:
Differentiating Equation (11) we obtain the approximating polynomial as follows:Evaluating Equation (18) at the points , , the following additional formulas are derived
4. Analysis of the Method
For any numerical scheme, it is very important to know the order of accuracy, the behavior of the local error, and the stability characteristics. To this end, the analysis of the proposed method is presented in this section.
For the order and truncation error of Equation (12), we consider the corresponding differential operatorswhere is sufficiently differentiable, and the coefficients , are the corresponding constant coefficients in the formulas in Equation (12). Expanding Equation (25) in Taylor series about and after collecting all the terms in the local truncation errors take the formwhere is the principal error constant and is the order of the corresponding formula.
For instance, considering the first formula in Equation (12), the formula in Equation (26) takes the form
So that, expanding Equation (27) in Taylor series, we arrive atwhere the local truncation error is as expressed in Equation (26).
Following similar pattern for other formulas in Equation (12), the local truncation errors are given as
For the formulas in Equations (19)–(24), the local truncation errors may be obtained similarly. From the above results, the order of each of the formulas in Equations (12)–(17) is . This is the same order of the formulas in Equations (19)–(24).
4.1. Zero Stability
For a block method, zero stability is such that, the roots of the characteristic equation given by must satisfy and for those with the multiplicity does not exceed , [10].
We note that the method in Equations (12)–(17) and Equations (19)–(24) may be written in matrix form as follows:whereand , are corresponding matrices of coefficients. Zero-stability implies the stability of the system in Equation(34) as . Letting , the system in Equation (34) becomes
The roots of the characteristic equation are for and . Consequently, the proposed method is zero stable.
4.2. Convergence
The convergence of a LMM is guaranteed if it is consistent (with order ) and zero stable [13, 14]. Considering the analysis shown above, the method has order , and is zero stable. Hence the proposed method is convergent.
4.3. Linear Stability Analysis
For any numerical method, the linear stability analysis is an important aspect of the theoretical analysis of the method. The Dahlquist’s test equation given byis usually used in the linear stability analysis for numerical methods for second order differential equations. Nevertheless, since the general second order differential equations are the focus and Equation (42) does not contain the first derivative, then the following linear test equation is considered [15]
For , the solutions of Equation (43) are bounded and go to zero when .
We employ the strategy by Singh and Ramos [16] to illustrate the procedure for obtaining the absolute stability region (the region in the -complex plane where the numerical method mimics the qualitative behavior of the exact solutions). The block method generated by the formulas in Equations (11) and (18) has 12 equations in which there are six different terms of derivatives: , , , , , , , and four intermediate values , , , . All these terms are eliminated from the system of equations using the Mathematica system so that the following recurrence equation with , , is obtainedwhere , andThe characteristic Equation (44) isTo determining the region of stability, the roots of this equation must have absolute values less than unity. Solving Equation (48), the roots of the characteristic equation arewhere
The plot of the region in which , is shown in Figure 1. This shows the region of stability for the method derived, whose stability interval is (0, 8.69809). This region is in the complex -plane where the roots of the characteristic Equation (48) are bounded in modulus by unity.
4.4. Formulation in Variable Step-Size Mode and Error Estimation
To gain efficiency when using the method given by Equations (12)–(17) and (19)–(24), it is convenient to have it formulated in variable step-size mode (VHM). In this sense, a lower order formula (LOF) must be considered to have an estimation of the local error (LE) at the end point of the two-step interval . The procedure is found to be less time consuming when the LOF uses values that have been previously obtained. To this end, for the block method constituted by Equations (12)–(17) and (19)–(24), we consider the LOF given bywhere denotes another approximation of . The LOF in (51) has the local truncation error
This was used to estimate the LE at the end-point . This error value gives the basis on how the step-size for the subsequent steps are determined in the course of implementation. For a given defined tolerance , the algorithm will change the step-size according to the following strategy. The Algorithm 1 below is applied for implementation of the proposed block method in VHM:
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||
4.5. Computational Procedure
The proposed block method is implemented in VHM in such a way that, on each block interval of the form and , where is a multiple of 2 so as give an integer number of iterations to reach the end of the integration interval . The system in (12)–(17) and (19)–(24) is solved by using the Newton’s method and taking the approximations provided by the Taylor formulas as starting values. These values are given as follows:for .
To apply the proposed method to a system of second order ODEs, we have the following procedurewhere
To solve this system we use again the Newton’s method. To obtain the estimations of the third derivatives of each component at and , we use the following formula
5. Numerical Examples
In this section we show the efficiency of the derived method implemented in variable and fixed step-size modes, respectively, where applicable. As mentioned earlier, our method has a higher order compared to other methods mentioned in this paper. The comparison of numerical results shall be based on the number of integration subintervals, the number of functions evaluations, accuracy in terms of global errors and CPU times, where applicable. There is no method of its order or of a higher order that has been found in the literature for comparisons. The examples considered show the robustness of the derived method in solving general second order systems of ODEs.
The following notations are used in the course of this section:(.): predefined tolerance(.): initial step-size(.): number of steps(.)EMax: maximum absolute error along the integration interval(.)(.)Time: CPU time in seconds(.)VHM: hybrid method derived in this paper using variable step size(.)FHM: hybrid method derived in this paper using fixed step size(.)NF: number of function evaluations(.)OPTBM: method derived by Singh and Ramos [16].
Problem 1. Orbital problem considered by Singh and Ramos [16].
The exact solution of the Stiefel and Bettis orbital problem is . We thus compare the proposed method named VHM with , and the methods reported by Singh and Ramos [16].
The comparison in Table 1 was done with the proposed method in VHM. The OPTBM in [16] was formulated in VHM. SCOWE(6) by Ramos and Vigo-Aguiar [3] and I3P1B by Ismail et al. [17] are methods reported by Singh and Ramos [16]. It shows the CPU times for VHM and OPTBM, with VHM performing the best. The number of function evaluations reveals that VHM performed well, being the one with the lowest number. Overall, this shows that the proposed method is the most accurate in comparison to the methods considered.
Problem 2. Consider the nonlinear Duffing equation discussed by Jator [1].
with the solution where , , , , , , and . For the sake of comparison as reported by Jator [1], the maximum global error is given in the form , where , AE is the absolute error at the endpoint of the integration interval. Table 2 shows the maximum errors obtained from our method solved in the VHM with and and fixed step-size mode, respectively, and compared to the hybrid linear multistep method (HLMM) of order seven derived by Jator [1]. It can be seen that the number of function evaluations is lower using our method.
Problem 3. Consider the Van Der Pol oscillator discussed by Allogmany and Ismail [18].
where the parameter shows the nonlinearity and the strength of damping with the value 0.005. It is a known fact that this problem does not have an exact solution. We compare the numerical solution obtained using our method implemented in the VHM and the in-built numerical scheme called up with NDSolve in Mathematica 12.0. Table 3 shows the solution obtained over the domain [0, 10] at varying points.
Figure 2 shows that the solution obtained using the VHM agrees with that obtained with NDSolve.
Problem 4. The Kepler’s problem is considered as discussed by Jator and King [19] and given by
whose exact solution is ; . The numerical results for this problem are obtained with as shown in Table 4.
For the proposed method, the initial step taken determines the number of steps . In Table 4, used and achieved are , , , , and , . It can be seen that the number of steps achieved in the proposed method and the maximum error are better than those by Jator and King [19].
Problem 5. Consider the linear system of second order by Majid et al. and Singh and Ramos [2, 16].
The exact solution is given as and . The numerical results of the problem are obtained for with as shown in Table 5.
In Table 5, the VHM with a small number of integration subintervals produced a better result with small global error than other methods. This reveals clearly the efficiency of our method as compared to VOPTBM, VJATOR, and the MATLAB ODE solvers —ode45 and ode113.
6. Conclusion
A third derivative hybrid block method for the numerical solution of general second-order initial value problems has been derived in this paper and implemented in both variable and fixed step-size modes, respectively. First, the method is derived using the collocation approach at uniform off grid points and considering up to the third derivatives. A variable step-size formulation was obtained, seeking to improve the efficiency of the method. The examples presented displayed good performance of the proposed method which provides smaller errors and less number of steps as shown when compared with other methods in the cited literature.
Data Availability
The data used to support the findings in this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Open access funding enabled and organized by CRUE-BUCLE Gold.