Abstract
An alternative block method for solving fifth-order initial value problems (IVPs) is proposed with an adaptive strategy of implementing variable step size. The derived method is designed to compute four solutions simultaneously without reducing the problem to a system of first-order IVPs. To validate the proposed method, the consistency and zero stability are also discussed. The improved performance of the developed method is demonstrated by comparing it with the existing methods and the results showed that the 4-point block method is suitable for solving fifth-order IVPs.
1. Introduction
Many natural processes or real-world problems can be translated into the language of mathematics [1–4]. The mathematical formulation of physical phenomena in science and engineering often leads to a differential equation, which can be categorized as an ordinary differential equation (ODE) and a partial differential equation (PDE). This formulation will explain the behavior of the phenomenon in detail. The search for solutions of real-world problems requires solving ODEs and thus has been an important aspect of mathematical study. For many interesting applications, an exact solution may be unattainable, or it may not give the answer in a convenient form. The reliability of numerical approximation techniques in solving such problems has been proven by many researchers as the role of numerical methods in engineering problems solving has increased dramatically in recent years. Thus a numerical approach has been chosen as an alternative tool for approximating the solutions consistent with the advancement in technology.
Commonly, the formulation of real-world problems will take the form of a higher order differential equation associated with its initial or boundary conditions [4]. In the literature, a mathematical model in the form of a fifth-order differential equation, known as Korteweg-de Vries (KdV) equation, has been used to describe several wave phenomena depending on the values of its parameters [2, 3, 5, 6]. The KdV equation is a PDE and researchers have tackled the problem analytically and numerically. It is also noted that in certain cases by using different approaches the KdV might be transformed into a higher order ODE [7]. To date, there are a number of studies that have proposed solving fifth-order ODE directly [8, 9]. Hence, the purpose of the present paper is to solve directly the fifth-order IVPs with the implementation of a variable step size strategy. The fifth-order IVP with its initial conditions is defined as
Conventionally, (1) will be converted to a system of first-order ODEs by a simple change of variables. However, it will increase the computational cost in terms of function evaluation and thus will affect the computational time. This drawback is obviously seen when dealing with a higher order problem. Furthermore, [10] also has remarked that the block method is far more cost-effective when it is implemented in direct integration. Hence, several researchers [11–16] have shown an interest in the development of direct integration methods. A direct integration method of variable order and step size for solving systems of nonstiff higher order ODEs has been discussed in [11] whereby [12] has proposed an algorithm based on collocation of the differential system at selected grid points for direct solution of general second-order ODEs. In addition, [13] has used the Gaussian method in order to solve fourth-order differential equations directly. However, it requires a tedious computation as well, since it consists of higher order partial derivatives of Taylor series algorithm which supplies the starting values. Jator and Li [15] have proposed the linear multistep method (LMM) for solving general second-order IVPs directly. The method is self-starting, so it involves less computational time by avoiding incorporating subroutines to supply the starting values.
Thus far, a number of researchers have concerned themselves with developing a numerical method based on block features, and the characteristic feature of the block method is that in each application it generates a set of solutions concurrently [10]. Rosser [10] also has remarked that the implementation of block method in numerical computation will reduce the computational cost by reducing the number of function evaluations. Shampine and Watts [17] have constructed an -stable implicit one-step block method and Cash [18] has studied block methods based upon the Runge-Kutta method for the numerical solution of nonstiff IVPs. Furthermore [19] has used the self-starting LMM to solve second-order ODEs in a block-by-block fashion and recently [20] has constructed a predictor-corrector scheme 3-point block method with the implementation of variable step size. This research is an extension of the work in [20] in which the solution is computed at three points concurrently and it shows the satisfactory numerical results obtained when solving general higher order ODEs.
An increasing amount of literature is devoted to variable step size implementations of numerical methods [11, 21, 22]. The practicality of varying the step size for block method has been justified by [10]. This strategy is an attempt to reduce the computational cost as well as maintaining the accuracy. The Falkner method with variable step size implementation for the numerical solution of second-order IVPs has been employed in [21]. Although the implementation of the method involves varying the step size and solving directly, the computation is still tedious since the coefficients of the formulae must be calculated every time the step size is changed. On the contrary, the present work will store all the integration coefficients in the code in order to avoid the tedious calculations of the divided differences.
2. Methodology
2.1. Derivation of 4-Point Block Method
The basic approach of numerical methods for integration is performed by subdividing the interval of integration into certain subintervals. The proposed method was based on concurrent computation; hence the closed finite interval was subdivided into a series of blocks and each block contains four equal subintervals as illustrated in Figure 1.
Initially, (1) was integrated five times over the corresponding interval: , , , for first, second, third, and fourth point, respectively. The integration was started by replacing the function with the interpolating function which was generated from Lagrange polynomials. A set of points was interpolated for deriving predictor and corrector formulae, respectively. Let the Lagrange polynomial, , be written aswhereNine points were interpolated in (2) with set to be eight for deriving the corrector and thus one point less for the predictor formula. Then, the integration process was proceeded by substituting and in (2). Consistent with the number of interpolation points involved in deriving the formulae, predictor and corrector formulae were obtained in terms of variables and . The variables and refer to the distance ratio between current and previous point as a result of implementation variable step size strategy in the proposed method.
In this work, the selection of the next step size could be increased by a factor of or maintained by and for halving the current step size. This step size changing technique was limited in order to minimize the number of coefficients to be stored as well as reducing the computational storage. The increment of step size was also limited to doubling in order to control the accuracy of the computation [10]. The compact form of the 4-point block method is presented in where are the coefficients of the formulae to be calculated, is the number of points , is the number of times (1) will be integrated, and is the number of terms when the equation is integrated. The values of , and , were considered for deriving the predictor and corrector formulae, respectively. After further simplification, the associated corrector formulae of the 4-point block method when are represented below. Integrate once: Integrate twice: Integrate thrice: Integrate four times: Integrate five times:
2.2. Order and Convergence of the Method
The matrix differential equation of the derived method is given aswhere , and are the coefficients of the developed method. Consequently, the order of 4-point block method can be determined using the following formulae:As a result, a 4-point block method of order nine is developed with and the error constant obtained is given asHence, the consistency of 4-point block method is proven according to the definition in [23]. The analysis of zero stability for the developed method is tested using a similar approach as presented in [24] and the first characteristic polynomial obtained is . It is clearly seen that the roots are and . Thus from Theorem 2.1 in [23], the convergence of the proposed method is asserted.
3. Implementation
Throughout this section, the implementation of 4-point block method for solving fifth-order IVPs will be explained in detail. The code starts by finding the values of in the initial block using Euler method. However, it should be noted that Euler method will act only as a fundamental building block. Then the 4-point block method will be applied until the end of the interval. As stated earlier, the proposed method is implemented in the mode of predicting and correcting. In order to preserve the accuracy, the step will succeed if the local truncation error (LTE) is less than the specified error tolerance (TOL) such that where and are the corrector value of at the last point for each block with iterations. If (13) is satisfied, the new step size will be calculated via the step size increment formula. Otherwise, the current step size will be reduced by half. The step size increment formula is defined aswhere and denote the current and previous step size, respectively, with value of 0.5 for safety factor and is the order of corrector formulae. To show the accuracy and efficiency of the proposed method, the computational errors will be reported, equal toDifferent values of and represent the type of error test which will be considered; namely, , are for the absolute error test; , represent the mixed error test; and , correspond to the relative error test. Here we will use the mixed error test and the maximum error is calculated by
4. Numerical Results
To illustrate the technique proposed in the preceding sections, two test problems are solved and the results obtained compared with the method proposed in [8, 9, 25] and the ODE solver in MATLAB ode45. The work done by [25] involved the solving of the first-order ODEs using 4-point block method using variable step size. This means that (1) needs to be reduced into a system of first-order IVPs whereby the methods proposed by [8, 9] are for solving (1) directly with the implementation of constant step size. The notations used in Tables 1 and 2 are listed below: TOL: error tolerance limit; MTD: method used; TS: total steps taken; MAXE: maximum error of the computed solution; AVERR: average error of the computed solution; FC: total function calls; ode45: Runge-Kutta-Dormand-Prince ODE solver; 4PHODE: implementation of the 4-point block method in this research; 4P1FI: numerical results in [25]; KAYODE (a): numerical results in [8]; KAYODE (b): numerical results in [9].
Problem 1. Consider the following:
Solution is as follows: .
Source is [8].
Problem 2. Consider the following:
Solution is as follows: .
Source is [8].
5. Discussion
In this section the performances of 4PHODE, ode45, 4P1FI, KAYODE (a), and KAYODE (b) are discussed in terms of three parameters, namely, total steps taken, accuracy, and total function evaluations. It is apparent from these tables, mostly at tolerances , , and , that 4PHODE gives better accuracy compared to ode45, whereas at tolerances and , ode45 is one decimal more accurate compared to 4PHODE for both problems. As both tables show, the total steps taken by the 4PHODE reduce by nearly half to ode45 at tolerance . Although, at other tolerance, ode45 requires lesser steps to compute the solution, this result may be explained by the fact that the initial step size generated by 4PHODE is extremely small in order for the method control of the accuracy. From the data in Tables 1 and 2, it is apparent that the number of function calls is likely to be related to the number of steps taken.
In the comparison of block-by-block method, the numerical results obtained in Tables 1 and 2 demonstrate that the proposed method 4PHODE always requires fewer steps in converging to the exact solution compared to 4P1FI. This gain becomes more obvious as the tolerance decreases. However, the 4P1FI achieves better accuracy than 4PHODE. Even so, the maximum error for 4PHODE is still within the specified tolerance. Another issue that had to be addressed was the total number of function evaluations used by each method. The 4PHODE solved all the problems at much lower cost than 4P1FI and this superiority is most apparent at finer tolerance. This result is in agreement with [10] which states the drawback of using the reduction approach in the implementation of simultaneous computations.
Comparing the results with the method proposed by KAYODE (a) and KAYODE (b), for Problem 1, 4PHODE outperformed both KAYODE (a) and KAYODE (b) in terms of accuracy and total steps taken. While, for Problem 2, KAYODE (b) has lesser total steps taken, however, 4PHODE still has superiority in terms of accuracy.
6. Conclusion
The overall performance revealed that the 4-point block method is best to be implemented in a direct integration approach as it required much less storage than the reduction method while still maintaining an acceptable accuracy. Besides that, the results of this study also indicate that the developed method has better accuracy compared to the existing methods. Hence it can be said that 4PHODE is one of the alternative methods that can be used for solving fifth-order ODEs.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
Grateful acknowledgement is made to the Ministry of Higher Education, Malaysia, for financial support with Grant no. RAGS/2013/FKEKK/SG04/01/B00037.