Research Article
Solving Nonstiff Higher-Order Ordinary Differential Equations
Using 2-Point Block Method Directly
Hazizah Mohd Ijam,1Mohamed Suleiman,2Ahmad Fadly Nurullah Rasedee,2
Norazak Senu,1Ali Ahmadian,1,2 and Soheil Salahshour3
1Department of Mathematics, Faculty of Science, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia
2Institute for Mathematical Research, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia
3Department of Computer Engineering, Mashhad Branch, Islamic Azad University, Mashhad, Iran
Correspondence should be addressed to Ali Ahmadian; ahmadian.hosseini@gmail.com
Received 18 July 2014; Accepted 23 August 2014; Published 17 September 2014
Academic Editor: Dumitru Baleanu
Copyright © 2014 Hazizah Mohd Ijam et al. is 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.
We describe the development of a 2-point block backward difference method (2PBBD) for solving system of nonstiff higher-order
ordinary differential equations (ODEs) directly. e method computes the approximate solutions at two points simultaneously
within an equidistant block. e integration coefficients that are used in the method are obtained only once at the start of
the integration. Numerical results are presented to compare the performances of the method developed with 1-point backward
difference method (1PBD) and 2-point block divided difference method (2PBDD). e result indicated that, for finer step sizes, this
method performs better than the other two methods, that is, 1PBD and 2PBDD.
1. Introduction
In this paper, we consider the system of 𝑑th order ODEs of
the form
𝑦𝑖(𝑑𝑖)=𝑓
𝑖(𝑥,̃
𝑌), 𝑖=1,2,…,𝑠, (1)
with ̃
𝑌(𝑎)=̃
𝜂in the interval 𝑎≤𝑥≤𝑏,where
̃
𝑌(𝑥)=(𝑦1,…,𝑦(𝑑1−1)
1,…,𝑦𝑠,…,𝑦(𝑑𝑠−1)
𝑠),
̃
𝜂=(𝜂1,…,𝜂(𝑑1−1)
1,…,𝜂𝑠,…,𝜂(𝑑𝑠−1)
𝑠). (2)
For simplicity of discussion and without loss of generality, we
consider the single equation
𝑦(𝑑) =𝑓(𝑥,̃
𝑌),̃
𝑌(𝑎)=̃
𝜂, (3)
where
̃
𝑌𝑇=(𝑦,𝑦,…,𝑦(𝑑−1)),̃
𝜂𝑇=(𝜂,𝜂,…,𝜂(𝑑−1)).(4)
As shown in Figure 1, here the 2-point block method, the
interval [𝑎,𝑏], is divided into series of blocks with each block
containing two points; that is, 𝑥𝑛−1 and 𝑥𝑛is the first block
while 𝑥𝑛+1 and 𝑥𝑛+2 is the second block, where solutions to
(3)aretobecomputed.
Previous works on block method for solving (3)directly
are given by Milne [1], Rosser [2], Shampine and Watts [3],
and Chu and Hamilton [4]. According to Omar [5], both
implicit and explicit block Adams methods in their divided
difference form are developed for the solution of higher-
order ODEs. Majid [6] has derived a code based on the
variable step size and order of fully implicit block method
to solve nonstiff higher-order ODEs directly. Ibrahim [7]
has developed a new block backward differentiation formula
method of variable step size for solving first- and second-
order ODEs directly. Suleiman et al. [8]haveintroducedone–
point backward difference methods for solving higher-order
ODEs. Hence, this motivates us to extend the method to block
method in solving nonstiff higher-order ODEs.
Hindawi Publishing Corporation
Abstract and Applied Analysis
Volume 2014, Article ID 867095, 13 pages
http://dx.doi.org/10.1155/2014/867095
2Abstract and Applied Analysis
xn−2 xn−1 xnxn+2
xn+1
Figure 1: 2-point method.
2. The Formulation of the Predict-Evaluate-
Correct-Evaluate (PECE) Multistep Block
Method in Its Backward Difference Form
(MSBBD) for Nonstiff Higher-Order ODEs
e code developed will be using the PECE mode with con-
stant stepsize. e predictor and corrector for first and second
point will have the following form.
Predictor:
pr𝑦(𝑑−𝑡)
𝑛+𝑟 =𝑡−1
∑
𝑖=0 ℎ𝑖
𝑖!𝑦(𝑑−𝑡+𝑖)
𝑛+ℎ𝑡𝑘−1
∑
𝑖=0𝛾𝑟,𝑡,𝑖∇𝑖𝑓𝑛,(5)
where 𝛾𝑟,𝑡,𝑖 is coefficient for predictor for 𝑟=1,2and 𝑡=
1,2,…,𝑑.
Corrector:
𝑦(𝑑−𝑡)
𝑛+𝑟 =𝑡−1
∑
𝑖=0 ℎ𝑖
𝑖!𝑦(𝑑−𝑡+𝑖)
𝑛+ℎ𝑡𝑘
∑
𝑖=0𝛾∗
𝑟,𝑡,𝑖∇𝑖𝑓𝑛+𝑟,(6)
where 𝛾∗
𝑟,𝑡,𝑖 is coefficient for corrector for 𝑟=1,2and 𝑡=
1,2,…,𝑑.
We also formulate the corrector in terms of the predictor.
Both points 𝑦𝑛+1 and 𝑦𝑛+2 canbewrittenas
𝑦(𝑑−𝑡)
𝑛+1 =pr𝑦(𝑑−𝑡)
𝑛+1 +ℎ[𝛾1,𝑡,𝑘∇𝑘𝑓𝑛+1], (7)
𝑦(𝑑−𝑡)
𝑛+2 =pr𝑦(𝑑−𝑡)
𝑛+2 +ℎ[𝛾2,𝑡,𝑘∇𝑘𝑓𝑛+2 −𝛾2,𝑡,𝑘−1∇𝑘+1𝑓𝑛+2]. (8)
We derived the formulation for both the predictor and cor-
rector.
3. Derivation for Higher-Order Explicit
Integration Coefficients
3.1. For the First Point. e derivation for up to third-order
explicit integration coefficients for the first point 𝑦𝑛+1 has
been given by Suleiman et al. [8].
3.2. For the Second Point. Integrating (3)onceyields
𝑦(𝑑−1) (𝑥𝑛+2)=𝑦(𝑑−1) (𝑥𝑛)
+∫𝑥𝑛+2
𝑥𝑛𝑓(𝑥,𝑦,𝑦,𝑦,…,𝑦(𝑑−1))𝑑𝑥. (9)
Let 𝑃𝑛(𝑥)be the interpolating polynomial which interpolates
the 𝑘values (𝑥𝑛,𝑓𝑛),(𝑥𝑛−1,𝑓𝑛−1),…,(𝑥𝑛−𝑘+1,𝑓𝑛−𝑘+1);then
𝑃𝑛(𝑥)=𝑘−1
∑
𝑖=0(−1)𝑖(−𝑠
𝑖)∇𝑖𝑓𝑛.(10)
Approximating 𝑓in (6)with𝑃𝑛(𝑥)and letting
𝑥=𝑥𝑛+𝑠ℎ or 𝑠=𝑥−𝑥𝑛
ℎ(11)
gives
𝑦(𝑑−1) (𝑥𝑛+2)=𝑦(𝑑−1) (𝑥𝑛)+∫2
0
𝑘−1
∑
𝑖=0(−1)𝑖(−𝑠
𝑖)∇𝑖𝑓𝑛ℎ𝑑𝑠
(12)
or
𝑦(𝑑−1) (𝑥𝑛+2)=𝑦(𝑑−1) (𝑥𝑛)+ℎ𝑘−1
∑
𝑖=0𝛾2,1,𝑖∇𝑖𝑓𝑛,(13)
where
𝛾2,1,𝑖 =(−1)𝑖∫2
0(−𝑠
𝑖)𝑑𝑠. (14)
Define the generating function 𝐺1(𝑡)for the coefficient 𝛾2,1,𝑖
as follows:
𝐺1(𝑡)=∞
∑
𝑖=0𝛾2,1,𝑖𝑡𝑖.(15)
Substituting 𝛾2,1,𝑖 in (14)into𝐺1(𝑡)gives
𝐺1(𝑡)=∞
∑
𝑖=0(−𝑡)𝑖∫2
0(−𝑠
𝑖)𝑑𝑠,
𝐺1(𝑡)=∫2
0(1−𝑡)−𝑠𝑑𝑠,
𝐺1(𝑡)=∫2
0𝑒−𝑠log(1−𝑡)𝑑𝑠,
(16)
which leads to
𝐺1(𝑡)=−[(1−𝑡)−2
log (1−𝑡)−1
log (1−𝑡)]. (17)
Equation (17)canbewrittenas
−(∞
∑
𝑖=0𝛾2,1,𝑖𝑡𝑖)log (1−𝑡)=(2−𝑡)[𝑡
(1−𝑡)2](18)
or
(𝛾2,1,0 +𝛾2,1,1𝑡+𝛾2,1,2𝑡2+⋅⋅⋅)(𝑡+𝑡2
2+𝑡3
3+⋅⋅⋅)
=(2−𝑡)(𝑡+2𝑡2+3𝑡3+⋅⋅⋅). (19)
Abstract and Applied Analysis 3
Hence, the coefficients of 𝛾2,1,𝑖 are given by
𝑘
∑
𝑖=0 (𝛾2,1,𝑖
𝑘−𝑖+1)=𝑘+2,
𝛾2,1,𝑘 =(𝑘+2)−𝑘−1
∑
𝑖=0 𝛾2,1,𝑖
(𝑘−𝑖+1), 𝑘=1,2,…, 𝛾2,1.0 =2.
(20)
Integrating (1)twiceyields
𝑦(𝑑−2) (𝑥𝑛+2)=𝑦(𝑑−2) (𝑥𝑛)+ℎ𝑦(𝑑−1) (𝑥𝑛)+ℎ2𝑘−1
∑
𝑖=0𝛾2,2,𝑖∇𝑖𝑓𝑛.
(21)
Substituting 𝑥with 𝑠gives
𝛾2,2,𝑖 =(−1)𝑖∫2
0(2−𝑠)
1! (−𝑠
𝑖)𝑑𝑠. (22)
e generating function of the coefficient 𝛾2,2,𝑖 is defined as
follows:
𝐺2(𝑡)=∞
∑
𝑖=0𝛾2,2,𝑖𝑡𝑖.(23)
Substituting (22)into𝐺2(𝑡)above gives
𝐺2(𝑡)=∫2
0(2−𝑠)
1! 𝑒−𝑠log(1−𝑡)𝑑𝑠. (24)
Substituting 𝐺1(𝑡)into (24)yields
𝐺2(𝑡)=1
1![2
log (1−𝑡)−1!𝐺1(𝑡)
log (1−𝑡)]. (25)
Equation (25)canbewrittenas
(∞
∑
𝑖=0𝛾2,2,𝑖𝑡𝑖)log (1−𝑡)=1
1
or
(𝛾2,2,0 +𝛾2,2,1𝑡+𝛾2,2,2𝑡2+⋅⋅⋅)(𝑡+𝑡2
2+𝑡3
3+⋅⋅⋅)
=1
1![−2+1!(𝛾2,1,0 +𝛾2,1,1𝑡+𝛾2,1,2𝑡2+⋅⋅⋅)]. (27)
Hence the coefficients of 𝛾2,2,𝑘 in relation to coefficients of
the previous order 𝛾2,1,𝑘 are given by
𝑘
∑
𝑖=0 𝛾2,2,𝑖
𝑘−𝑖+1=𝛾2,1,𝑘+1,(28)
𝛾2,2,0 =𝛾2,1,1,
𝛾2,2,𝑘 =𝛾2,1,𝑘+1 −𝑘−1
∑
𝑖=0 𝛾2,2,𝑖
𝑘−𝑖+1, 𝑘=1,2,…. (29)
By using the same process previously, we note that for inte-
grating (𝑑−1)times yield
𝐺(𝑑−1) (𝑡)=∫2
0(2−𝑠)(𝑑−2)
(𝑑−2)!𝑒−𝑠log(1−𝑡)𝑑𝑠, (30)
𝐺(𝑑−1) (𝑡)=1
(𝑑−2)![2(𝑑−2)
log (1−𝑡)−(𝑑−2)!𝐺(𝑑−2) (𝑡)
log (1−𝑡)],
(31)
and, from (29), we get
𝛾2,(𝑑−1),0 =𝛾2,(𝑑−2),1,
𝛾2,(𝑑−1),𝑘 =𝛾2,(𝑑−2),𝑘+1 −𝑘−1
∑
𝑖=0 𝛾2,(𝑑−1),𝑖
𝑘−𝑖+1, 𝑘=1,2,…. (32)
Integrating (𝑑)times yield
𝑦(𝑥𝑛+2)=𝑦(𝑥𝑛)+ℎ𝑦(𝑥𝑛)+⋅⋅⋅+ ℎ(𝑑−1)
(𝑑−1)!𝑦(𝑑−1) (𝑥𝑛)
+∫𝑥𝑛+2
𝑥𝑛
(𝑥𝑛+2 −𝑥)(𝑑−1)
(𝑑−1)!
×𝑓(𝑥,𝑦,𝑦,𝑦,…,𝑦(𝑑−1))𝑑𝑥 (33)
or, in the backward difference formulation, given by
𝑦(𝑥𝑛+2)=𝑦(𝑥𝑛)+ℎ𝑦(𝑥𝑛)+⋅⋅⋅+ ℎ(𝑑−1)
(𝑑−1)!𝑦(𝑑−1) (𝑥𝑛)
+ℎ(𝑑)𝑘−1
∑
𝑖=0𝛾2,(𝑑),𝑖∇𝑖𝑓𝑛,
(34)
where
𝛾2,(𝑑),𝑖 =(−1)𝑖∫2
0(2−𝑠)(𝑑−1)
(𝑑−1)!(−𝑠
𝑖)𝑑𝑠. (35)
e generating function
𝐺(𝑑) (𝑡)=∞
∑
𝑖=0𝛾2,(𝑑),𝑖𝑡𝑖.(36)
Substituting (35)into𝐺(𝑑)(𝑡)above yields
𝐺(𝑑) (𝑡)=∫2
0(2−𝑠)(𝑑−1)
(𝑑−1)!𝑒−𝑠log(1−𝑡)𝑑𝑠. (37)
4Abstract and Applied Analysis
As in (30), we now substitute 𝐺(𝑑−1)(𝑡)in (37)giving
𝐺(𝑑) (𝑡)=1
(𝑑−1)![2(𝑑−1)
log (1−𝑡)−(𝑑−1)!𝐺(𝑑−1) (𝑡)
log (1−𝑡)]. (38)
Equation (38)canbewrittenas
(∞
∑
𝑖=0𝛾2,(𝑑),𝑖𝑡𝑖)log (1−𝑡)
=1
(𝑑−1)![2(𝑑−1) −(𝑑−1)!𝐺(𝑑−1) (𝑡)]
(39)
or
(𝛾2,(𝑑),0 +𝛾2,(𝑑),1𝑡+𝛾2,(𝑑),2𝑡2+⋅⋅⋅)(𝑡+𝑡2
2+𝑡3
3+⋅⋅⋅)
=1
(𝑑−1)!
×[ −2(𝑑−1)
+(𝑑−1)!(𝛾2,(𝑑−1),0 +𝛾2,(𝑑−1),1𝑡+𝛾2,(𝑑−1),2𝑡2+⋅⋅⋅)].
(40)
Hence the coefficients of 𝛾2,(𝑑),𝑘 in relation to coefficients of
the previous order 𝛾2,(𝑑−1),𝑘 are given by
𝑘
∑
𝑘−𝑖+1=𝛾2,(𝑑−1),𝑘+1,
As in the previous derivation, we choose
𝑥=𝑥𝑛+2 +𝑠ℎ or 𝑠=𝑥−𝑥𝑛+2
ℎ.(44)
Replacing 𝑥by 𝑠yields
𝑦(𝑑−1) (𝑥𝑛+2)=𝑦(𝑑−1) (𝑥𝑛)+∫0
−2
𝑘
∑
𝑖=0(−1)𝑖(−𝑠
𝑖)∇𝑖𝑓𝑛+2ℎ𝑑𝑠.
(45)
Simplify
𝑦(𝑑−1) (𝑥𝑛+2)=𝑦(𝑑−1) (𝑥𝑛)+ℎ𝑘
∑
𝑖=0𝛾∗
2,1,𝑖∇𝑖𝑓𝑛+2,(46)
where
𝛾∗
2,1,𝑖 =(−1)𝑖∫0
−2 (−𝑠
𝑖)𝑑𝑠. (47)
Define the generating function 𝐺∗
1(𝑡)for the coefficient 𝛾∗
2,1,𝑖
as follows:
𝐺∗
1(𝑡)=∞
∑
𝑖=0𝛾∗
2,1,𝑖𝑡𝑖.(48)
or
𝐺∗
1(𝑡)=∞
∑
𝑖=0(−𝑡)𝑖∫0
−2 (−𝑠
𝑖)𝑑𝑠, (49)
