International Journal of Mathematics and Mathematical Sciences
Volume 2016 (2016), Article ID 9823147, 8 pages
http://dx.doi.org/10.1155/2016/9823147

Research Article

A 4-Point Block Method for Solving Higher Order Ordinary Differential Equations Directly

Nazreen Waeleh¹ and Zanariah Abdul Majid2,3

¹Faculty of Electronic & Computer Engineering, Universiti Teknikal Malaysia Melaka (UTeM), 76100 Melaka, Malaysia
²Department of Mathematics, Faculty of Science, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia
³Institute for Mathematical Research, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia

Received 18 April 2016; Revised 23 June 2016; Accepted 30 June 2016

Academic Editor: Harvinder S. Sidhu

Copyright © 2016 Nazreen Waeleh and Zanariah Abdul Majid. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

Figure 1: 4-point block method.

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 as

where

Nine 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 as

where , 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 as

Hence, 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 as

where 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 to

Different 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].

Table 1: Numerical results for solving Problem 1.

TOL
MTD
TS
MAXE
AVERR
FC

10⁻²
ode45
11
1.743 (−6)
2.244 (−6)
67

4P1FI
23
1.640 (−4)
5.346 (−7)
210

4PHODE
20
9.031 (−9)
1.611 (−10)
186

10⁻⁴
ode45
11
1.641 (−6)
2.186 (−6)
67

4P1FI
37
9.727 (−7)
6.988 (−9)
339

4PHODE
27
7.059 (−10)
1.092 (−11)
243

10⁻⁶
ode45
18
5.234 (−7)
8.596 (−7)
109

4P1FI
50
5.582 (−9)
5.727 (−11)
447

4PHODE
35
1.009 (−9)
3.401 (−11)
307

10⁻⁸
ode45
44
5.485 (−9)
9.472 (−9)
265

4P1FI
96
7.597 (−11)
2.881 (−13)
799

4PHODE
42
1.838 (−10)
2.619 (−12)
363

10⁻¹⁰
ode45
108
5.552 (−11)
9.551 (−11)
649

4P1FI
788
1.834 (−14)
2.554 (−15)
6331

4PHODE
50
1.023 (−10)
3.443 (−12)
439

KAYODE (a)
Not stated
1.638 (−6)
Not stated
Not stated

KAYODE (b)
100
5.082 (−7)
Not stated
Not stated

Table 2: Numerical results for solving Problem 2.

TOL
MTD
TS
MAXE
AVERR
FC

10⁻²
ode45
11
5.027 (−3)
3.806 (−3)
67

4P1FI
58
2.528 (−7)
3.002 (−8)
490

4PHODE
28
5.537 (−6)
2.245 (−7)
250

10⁻⁴
ode45
13
1.967 (−4)
1.253 (−4)
79

4P1FI
68
9.819 (−9)
1.117 (−9)
571

4PHODE
36
4.193 (−7)
1.432 (−8)
315

10⁻⁶
ode45
25
4.592 (−7)
1.870 (−7)
151

4P1FI
86
1.219 (−10)
9.992 (−12)
715

4PHODE
46
1.003 (−7)
4.290 (−9)
395

10⁻⁸
ode45
58
1.861 (−9)
1.006 (−9)
349

4P1FI
135
1.363 (−12)
1.597 (−13)
1107

4PHODE
56
2.413 (−8)
1.012 (−9)
475

10⁻¹⁰
ode45
142
4.555 (−11)
2.026 (−11)
853

4P1FI
807
5.630 (−14)
1.250 (−14)
6483

4PHODE
71
1.949 (−10)
7.564 (−12)
595

KAYODE (a)
Not stated
2.225 (−7)
Not stated
Not stated

KAYODE (b)
10
3.156 (−4)
Not stated
Not stated

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.

References

M. T. Darvishi, S. Kheybari, and F. Khani, “A numerical solution of the Korteweg-de Vries equation by pseudospectral method using Darvishi's preconditionings,” Applied Mathematics and Computation, vol. 182, no. 1, pp. 98–105, 2006.

L. Jin, “Application of variational iteration method to the fifth-order KdV equation,” International Journal of Contemporary Mathematical Sciences, vol. 3, no. 5, pp. 213–221, 2008.

L. Kaur, “Generalized (G′/G)—expansion method for generalized fifth order KdV equation with time-dependent coefficients,” Mathematical Sciences Letters, vol. 3, no. 3, pp. 255–261, 2014.

M. Suleiman, Z. B. Ibrahim, and A. F. N. Bin Rasedee, “Solution of higher-order ODEs using backward difference method,” Mathematical Problems in Engineering, vol. 2011, Article ID 810324, 18 pages, 2011.

U. Goktas and W. Hereman, “Symbolic computation of conserved densities for systems of nonlinear evolution equations,” Journal of Symbolic Computation, vol. 24, no. 5, pp. 591–621, 1997.

N. Khanal, R. Sharma, J. Wu, and J.-M. Yuan, “A dual-Petrov-Galerkin method for extended fifth-order Korteweg-de Vries type equations,” Discrete and Continuous Dynamical Systems, pp. 442–450, 2009.

J. Li and Z. Qiao, “Explicit soliton solutions of the Kaup-Kupershmidt equation through the dynamical system approach,” Journal of Applied Analysis and Computation, vol. 1, no. 2, pp. 243–250, 2011.

S. J. Kayode and D. O. Awoyemi, “A multiderivative Collocation method for 5th order ordinary differential equations,” Journal of Mathematics and Statistics, vol. 6, no. 1, pp. 60–63, 2010.

S. J. Kayode, “An order seven continuous explicit method for direct solution of general fifth order ordinary differential equations,” International Journal of Differential Equations and Applications, vol. 13, no. 2, pp. 71–80, 2014.

J. B. Rosser, “A Runge-Kutta for all seasons,” SIAM Review, vol. 9, no. 3, pp. 417–452, 1967.

M. B. Suleiman, “Solving nonstiff higher order ODEs directly by the direct integration method,” Applied Mathematics and Computation, vol. 33, no. 3, pp. 197–219, 1989.

D. O. Awoyemi, “A new sixth-order algorithm for general second order ordinary differential equations,” International Journal of Computer Mathematics, vol. 77, no. 1, pp. 117–124, 2001.

S. J. Kayode, “An efficient zero-stable numerical method for fourth-order differential equations,” International Journal of Mathematics and Mathematical Sciences, vol. 2008, Article ID 364021, 10 pages, 2008.

B. T. Olabode and Y. Yusuph, “A new block method for special third order ordinary differential equations,” Journal of Mathematics and Statistics, vol. 5, no. 3, pp. 167–170, 2009.

S. N. Jator and J. Li, “A self-starting linear multistep method for a direct solution of the general second-order initial value problem,” International Journal of Computer Mathematics, vol. 86, no. 5, pp. 827–836, 2009.

N. Waeleh, Z. A. Majid, F. Ismail, and M. Suleiman, “Numerical solution of higher order ordinary differential equations by direct block code,” Journal of Mathematics and Statistics, vol. 8, no. 1, pp. 77–81, 2011.

L. F. Shampine and H. A. Watts, “Block implicit one-step methods,” Mathematics of Computation, vol. 23, pp. 731–740, 1969.

J. R. Cash, “Block Runge-Kutta methods for the numerical integration of initial value problems in ordinary differential equations. Part I. The nonstiff case,” Mathematics of Computation, vol. 40, no. 161, pp. 175–191, 1983.

S. N. Jator, “Solving second order initial value problems by a hybrid multistep method without predictors,” Applied Mathematics and Computation, vol. 217, no. 8, pp. 4036–4046, 2010.

N. Waeleh, Z. A. Majid, and F. Ismail, “A new algorithm for solving higher order IVPs of ODEs,” Applied Mathematical Sciences, vol. 5, no. 53–56, pp. 2795–2805, 2011.

J. Vigo-Aguiar and H. Ramos, “Variable stepsize implementation of multistep methods for y′′=f(x, y, y′),” Journal of Computational and Applied Mathematics, vol. 192, no. 1, pp. 114–131, 2006.

Z. A. Majid, N. A. Azmi, M. Suleiman, and Z. B. Ibrahaim, “Solving directly general third order ordinary differential equations using two-point four step block method,” Sains Malaysiana, vol. 41, no. 5, pp. 623–632, 2012.

J. D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & Sons, London, UK, 1973.

S. Ola Fatunla, “Block methods for second order ODEs,” International Journal of Computer Mathematics, vol. 41, no. 1-2, pp. 55–63, 1991.

Z. A. Majid and M. B. Suleiman, “Implementation of four-point fully implicit block method for solving ordinary differential equations,” Applied Mathematics and Computation, vol. 184, no. 2, pp. 514–522, 2007.