# Solving Nonstiff Higher-Order Ordinary Differential Equations Using 2-Point Block Method Directly

Authors: Hazizah Mohd Ijam; Mohamed Suleiman; Ahmad Fadly Nurullah Rasedee; Norazak Senu; Ali Ahmadian; Soheil Salahshour
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

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. The method computes the approximate solutions at two points simultaneously within an equidistant block. The 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). The result indicated that, for finer step sizes, this method performs better than the other two methods, that is, 1PBD and 2PBDD.


---

## Body

## 1. Introduction

In this paper, we consider the system ofdth order ODEs of the form
(1)yi(di)=fi(x,Y~),i=1,2,…,s,
with Y~(a)=η~ in the interval a≤x≤b, where
(2)Y~(x)=(y1,…,y1(d1-1),…,ys,…,ys(ds-1)),η~=(η1,…,η1(d1-1),…,ηs,…,ηs(ds-1)).
For simplicity of discussion and without loss of generality, we consider the single equation
(3)y(d)=f(x,Y~),Y~(a)=η~,
where
(4)Y~T=(y,y′,…,y(d-1)),η~T=(η,η′,…,η(d-1)).
As shown in Figure 1, here the 2-point block method, the interval [a,b], is divided into series of blocks with each block containing two points; that is, xn-1 and xn is the first block while xn+1 and xn+2 is the second block, where solutions to (3) are to be computed.Figure 1
2-point method.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] have introduced one-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.


## 2. The Formulation of the Predict-Evaluate-Correct-Evaluate (PECE) Multistep Block Method in Its Backward Difference Form (MSBBD) for Nonstiff Higher-Order ODEs

The code developed will be using the PECE mode with constant stepsize. The predictor and corrector for first and second point will have the following form.Predictor:(5)yprn+r(d-t)=∑i=0t-1hii!yn(d-t+i)+ht∑i=0k-1γr,t,i∇ifn,
where γr,t,i is coefficient for predictor for r=1,2 and t=1,2,…,d.Corrector:(6)yn+r(d-t)=∑i=0t-1hii!yn(d-t+i)+ht∑i=0kγr,t,i*∇ifn+r,
where γr,t,i* is coefficient for corrector for r=1,2 and t=1,2,…,d.We also formulate the corrector in terms of the predictor. Both pointsyn+1 and yn+2 can be written as
(7)yn+1(d-t)=yprn+1(d-t)+h[γ1,t,k∇kfn+1],(8)yn+2(d-t)=yprn+2(d-t)+h[γ2,t,k∇kfn+2-γ2,t,k-1∇k+1fn+2].
We derived the formulation for both the predictor and corrector.


## 3. Derivation for Higher-Order Explicit Integration Coefficients

### 3.1. For the First Point

The derivation for up to third-order explicit integration coefficients for the first pointyn+1 has been given by Suleiman et al. [8].


### 3.2. For the Second Point

Integrating (3) once yields
(9)y(d-1)(xn+2)=y(d-1)(xn)+∫xnxn+2f(x,y,y′,y′′,…,y(d-1))dx.
Let Pn(x) be the interpolating polynomial which interpolates the k values (xn,fn),(xn-1,fn-1),…,(xn-k+1,fn-k+1); then
(10)Pn(x)=∑i=0k-1(-1)i(-si)∇ifn.
Approximating f in (6) with Pn(x) and letting
(11)x=xn+shors=x-xnh
gives
(12)y(d-1)(xn+2)=y(d-1)(xn)+∫02∑i=0k-1(-1)i(-si)∇ifnhds
or
(13)y(d-1)(xn+2)=y(d-1)(xn)+h∑i=0k-1γ2,1,i∇ifn,
where
(14)γ2,1,i=(-1)i∫02(-si)ds.
Define the generating function G1(t) for the coefficient γ2,1,i as follows:
(15)G1(t)=∑i=0∞γ2,1,iti.
Substituting γ2,1,i in (14) into G1(t) gives
(16)G1(t)=∑i=0∞(-t)i∫02(-si)ds,G1(t)=∫02(1-t)-sds,G1(t)=∫02e-slog⁡(1-t)ds,
which leads to
(17)G1(t)=-[(1-t)-2log⁡(1-t)-1log⁡(1-t)].
Equation (17) can be written as
(18)-(∑i=0∞γ2,1,iti)log⁡(1-t)=(2-t)[t(1-t)2]
or
(19)(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯)(t+t22+t33+⋯)=(2-t)(t+2t2+3t3+⋯).
Hence, the coefficients of γ2,1,i are given by
(20)∑i=0k(γ2,1,ik-i+1)=k+2,γ2,1,k=(k+2)-∑i=0k-1γ2,1,i(k-i+1),k=1,2,…,γ2,1.0=2.
Integrating (1) twice yields
(21)y(d-2)(xn+2)=y(d-2)(xn)+hy(d-1)(xn)+h2∑i=0k-1γ2,2,i∇ifn.
Substituting x with s gives
(22)γ2,2,i=(-1)i∫02(2-s)1!(-si)ds.
The generating function of the coefficient γ2,2,i is defined as follows:
(23)G2(t)=∑i=0∞γ2,2,iti.
Substituting (22) into G2(t) above gives
(24)G2(t)=∫02(2-s)1!e-slog⁡(1-t)ds.
Substituting G1(t) into (24) yields
(25)G2(t)=11![2log⁡(1-t)-1!G1(t)log⁡(1-t)].
Equation (25) can be written as
(26)(∑i=0∞γ2,2,iti)log⁡(1-t)=11![2-1!G1(t)]
or
(27)(γ2,2,0+γ2,2,1t+γ2,2,2t2+⋯)(t+t22+t33+⋯)=11![-2+1!(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯)].Hence the coefficients ofγ2,2,k in relation to coefficients of the previous order γ2,1,k are given by
(28)∑i=0kγ2,2,ik-i+1=γ2,1,k+1,(29)γ2,2,0=γ2,1,1,γ2,2,k=γ2,1,k+1-∑i=0k-1γ2,2,ik-i+1,k=1,2,….
By using the same process previously, we note that for integrating (d-1) times yield
(30)G(d-1)(t)=∫02(2-s)(d-2)(d-2)!e-slog⁡(1-t)ds,(31)G(d-1)(t)=1(d-2)![2(d-2)log⁡(1-t)-(d-2)!G(d-2)(t)log⁡(1-t)],
and, from (29), we get
(32)γ2,(d-1),0=γ2,(d-2),1,γ2,(d-1),k=γ2,(d-2),k+1-∑i=0k-1γ2,(d-1),ik-i+1,k=1,2,….
Integrating (d) times yield
(33)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+∫xnxn+2(xn+2-x)(d-1)(d-1)!0000000×f(x,y,y′,y′′,…,y(d-1))dx
or, in the backward difference formulation, given by
(34)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+h(d)∑i=0k-1γ2,(d),i∇ifn,
where
(35)γ2,(d),i=(-1)i∫02(2-s)(d-1)(d-1)!(-si)ds.
The generating function
(36)G(d)(t)=∑i=0∞γ2,(d),iti.
Substituting (35) into G(d)(t) above yields
(37)G(d)(t)=∫02(2-s)(d-1)(d-1)!e-slog⁡(1-t)ds.
As in (30), we now substitute G(d-1)(t) in (37) giving
(38)G(d)(t)=1(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)(t)log⁡(1-t)].
Equation (38) can be written as
(39)(∑i=0∞γ2,(d),iti)log⁡(1-t)=1(d-1)![2(d-1)-(d-1)!G(d-1)(t)]
or
(40)(γ2,(d),0+γ2,(d),1t+γ2,(d),2t2+⋯)(t+t22+t33+⋯)=1(d-1)!×[-2(d-1)+(d-1)!(γ2,(d-1),0+γ2,(d-1),1t+γ2,(d-1),2t2+⋯)].
Hence the coefficients of γ2,(d),k in relation to coefficients of the previous order γ2,(d-1),k are given by
(41)∑i=0kγ2,(d),ik-i+1=γ2,(d-1),k+1,γ2,(d),0=γ2,(d-1),1,γ2,(d),k=γ2,(d-1),k+1-∑i=0k-1γ2,(d),ik-i+1,k=1,2,….


## 3.1. For the First Point

The derivation for up to third-order explicit integration coefficients for the first pointyn+1 has been given by Suleiman et al. [8].


## 3.2. For the Second Point

Integrating (3) once yields
(9)y(d-1)(xn+2)=y(d-1)(xn)+∫xnxn+2f(x,y,y′,y′′,…,y(d-1))dx.
Let Pn(x) be the interpolating polynomial which interpolates the k values (xn,fn),(xn-1,fn-1),…,(xn-k+1,fn-k+1); then
(10)Pn(x)=∑i=0k-1(-1)i(-si)∇ifn.
Approximating f in (6) with Pn(x) and letting
(11)x=xn+shors=x-xnh
gives
(12)y(d-1)(xn+2)=y(d-1)(xn)+∫02∑i=0k-1(-1)i(-si)∇ifnhds
or
(13)y(d-1)(xn+2)=y(d-1)(xn)+h∑i=0k-1γ2,1,i∇ifn,
where
(14)γ2,1,i=(-1)i∫02(-si)ds.
Define the generating function G1(t) for the coefficient γ2,1,i as follows:
(15)G1(t)=∑i=0∞γ2,1,iti.
Substituting γ2,1,i in (14) into G1(t) gives
(16)G1(t)=∑i=0∞(-t)i∫02(-si)ds,G1(t)=∫02(1-t)-sds,G1(t)=∫02e-slog⁡(1-t)ds,
which leads to
(17)G1(t)=-[(1-t)-2log⁡(1-t)-1log⁡(1-t)].
Equation (17) can be written as
(18)-(∑i=0∞γ2,1,iti)log⁡(1-t)=(2-t)[t(1-t)2]
or
(19)(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯)(t+t22+t33+⋯)=(2-t)(t+2t2+3t3+⋯).
Hence, the coefficients of γ2,1,i are given by
(20)∑i=0k(γ2,1,ik-i+1)=k+2,γ2,1,k=(k+2)-∑i=0k-1γ2,1,i(k-i+1),k=1,2,…,γ2,1.0=2.
Integrating (1) twice yields
(21)y(d-2)(xn+2)=y(d-2)(xn)+hy(d-1)(xn)+h2∑i=0k-1γ2,2,i∇ifn.
Substituting x with s gives
(22)γ2,2,i=(-1)i∫02(2-s)1!(-si)ds.
The generating function of the coefficient γ2,2,i is defined as follows:
(23)G2(t)=∑i=0∞γ2,2,iti.
Substituting (22) into G2(t) above gives
(24)G2(t)=∫02(2-s)1!e-slog⁡(1-t)ds.
Substituting G1(t) into (24) yields
(25)G2(t)=11![2log⁡(1-t)-1!G1(t)log⁡(1-t)].
Equation (25) can be written as
(26)(∑i=0∞γ2,2,iti)log⁡(1-t)=11![2-1!G1(t)]
or
(27)(γ2,2,0+γ2,2,1t+γ2,2,2t2+⋯)(t+t22+t33+⋯)=11![-2+1!(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯)].Hence the coefficients ofγ2,2,k in relation to coefficients of the previous order γ2,1,k are given by
(28)∑i=0kγ2,2,ik-i+1=γ2,1,k+1,(29)γ2,2,0=γ2,1,1,γ2,2,k=γ2,1,k+1-∑i=0k-1γ2,2,ik-i+1,k=1,2,….
By using the same process previously, we note that for integrating (d-1) times yield
(30)G(d-1)(t)=∫02(2-s)(d-2)(d-2)!e-slog⁡(1-t)ds,(31)G(d-1)(t)=1(d-2)![2(d-2)log⁡(1-t)-(d-2)!G(d-2)(t)log⁡(1-t)],
and, from (29), we get
(32)γ2,(d-1),0=γ2,(d-2),1,γ2,(d-1),k=γ2,(d-2),k+1-∑i=0k-1γ2,(d-1),ik-i+1,k=1,2,….
Integrating (d) times yield
(33)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+∫xnxn+2(xn+2-x)(d-1)(d-1)!0000000×f(x,y,y′,y′′,…,y(d-1))dx
or, in the backward difference formulation, given by
(34)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+h(d)∑i=0k-1γ2,(d),i∇ifn,
where
(35)γ2,(d),i=(-1)i∫02(2-s)(d-1)(d-1)!(-si)ds.
The generating function
(36)G(d)(t)=∑i=0∞γ2,(d),iti.
Substituting (35) into G(d)(t) above yields
(37)G(d)(t)=∫02(2-s)(d-1)(d-1)!e-slog⁡(1-t)ds.
As in (30), we now substitute G(d-1)(t) in (37) giving
(38)G(d)(t)=1(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)(t)log⁡(1-t)].
Equation (38) can be written as
(39)(∑i=0∞γ2,(d),iti)log⁡(1-t)=1(d-1)![2(d-1)-(d-1)!G(d-1)(t)]
or
(40)(γ2,(d),0+γ2,(d),1t+γ2,(d),2t2+⋯)(t+t22+t33+⋯)=1(d-1)!×[-2(d-1)+(d-1)!(γ2,(d-1),0+γ2,(d-1),1t+γ2,(d-1),2t2+⋯)].
Hence the coefficients of γ2,(d),k in relation to coefficients of the previous order γ2,(d-1),k are given by
(41)∑i=0kγ2,(d),ik-i+1=γ2,(d-1),k+1,γ2,(d),0=γ2,(d-1),1,γ2,(d),k=γ2,(d-1),k+1-∑i=0k-1γ2,(d),ik-i+1,k=1,2,….


## 4. Derivation for Higher-Order Implicit Integration Coefficients

### 4.1. For the First Point

The derivation for up to third-order implicit integration coefficients for the first pointyn+1 has been given by Suleiman et al. [8].


### 4.2. For the Second Point

Integrating (3) once yields
(42)y(d-1)(xn+2)=y(d-1)(xn)+∫xnxn+2f(x,y,y′,y′′,…,y(d-1))dx.
Let Pn(x) be the interpolating polynomial which interpolates the k values (xn,fn),(xn-1,fn-1),…,(xn-k+1,fn-k+1); then
(43)Pn(x)=∑i=0k(-1)i(-si)∇ifn+2.
As in the previous derivation, we choose
(44)x=xn+2+shors=x-xn+2h.
Replacing x by s yields
(45)y(d-1)(xn+2)=y(d-1)(xn)+∫-20∑i=0k(-1)i(-si)∇ifn+2hds.
Simplify
(46)y(d-1)(xn+2)=y(d-1)(xn)+h∑i=0kγ2,1,i*∇ifn+2,
where
(47)γ2,1,i*=(-1)i∫-20(-si)ds.
Define the generating function G1*(t) for the coefficient γ2,1,i* as follows:
(48)G1*(t)=∑i=0∞γ2,1,i*ti.
or
(49)G1*(t)=∑i=0∞(-t)i∫-20(-si)ds,(50)G1*(t)=∫-20(1-t)-sds,(51)G1*(t)=∫-20e-slog⁡(1-t)ds,
which leads to
(52)G1*(t)=-[1log⁡(1-t)-(1-t)2log⁡(1-t)].
For the case t=2, the approximate solution of y has the form
(53)y(d-2)(xn+2)=y(d-2)(xn)+hy(d-1)(xn)+∫xnxn+2(xn+2-x)(1)1!000000×f(x,y,y′,y′′,…,y(d-1))dx.
The coefficients are given by
(54)γ2,2,i*=(-1)i∫-20(-s)1!(-si)ds,
where γ2,2,i* are the coefficients of the backward difference formulation of (54) which can be represented by
(55)y(d-2)(xn+2)=y(d-2)(xn)+hy(d-1)(xn)+h2∑i=0kγ2,2,i*∇ifn+2.
Define the generating function of the coefficient γ2,2,i* as follows:
(56)G2*(t)=∑i=0∞γ2,2,i*ti.
Substituting (54) into G2*(t) above gives
(57)G2*(t)=∫-20(-s)1!e-slog⁡(1-t)ds.
Solving (57) with the substitution of (51) produces the relationship
(58)G2*(t)=11![2(1-t)2log⁡(1-t)-1!G1*(t)log⁡(1-t)].
By using the same process previously, we note that for integrating (d-1) times yield
(59)G(d-1)*(t)=∫-20(-s)(d-2)(d-2)!e-slog⁡(1-t)ds,(60)G(d-1)*(t)=1(d-2)![2(d-2)(1-t)2log⁡(1-t)-(d-2)!G(d-2)*(t)log⁡(1-t)].
Integrating (d) times yield
(61)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+∫xnxn+2(xn+2-x)(d-1)(d-1)!000000×f(x,y,y′,y′′,…,y(d-1))dx.
The coefficients are given by
(62)γ2,(d),i*=(-1)i∫-20(-s)(d-1)(d-1)!(-si)ds,
where γ2,(d),i* are the coefficients of the backward difference formulation of (62) which can be represented by
(63)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+h(d)∑i=0kγ2,(d),i*∇ifn+2.
Define the generating function G(d)*(t) of the coefficient γ2,(d),i* as follows:
(64)G(d)*(t)=∑i=0∞γ2,(d),i*ti.
Substituting (62) into G(d)*(t) above gives
(65)G(d)*(t)=∫-20(-s)(d-1)(d-1)!e-slog⁡(1-t)ds.
Solving (65) with the substitution of (59) produces the relationship
(66)G(d)*(t)=1(d-1)![2(d-1)(1-t)2log⁡(1-t)-(d-1)!G(d-1)*(t)log⁡(1-t)].


## 4.1. For the First Point

The derivation for up to third-order implicit integration coefficients for the first pointyn+1 has been given by Suleiman et al. [8].


## 4.2. For the Second Point

Integrating (3) once yields
(42)y(d-1)(xn+2)=y(d-1)(xn)+∫xnxn+2f(x,y,y′,y′′,…,y(d-1))dx.
Let Pn(x) be the interpolating polynomial which interpolates the k values (xn,fn),(xn-1,fn-1),…,(xn-k+1,fn-k+1); then
(43)Pn(x)=∑i=0k(-1)i(-si)∇ifn+2.
As in the previous derivation, we choose
(44)x=xn+2+shors=x-xn+2h.
Replacing x by s yields
(45)y(d-1)(xn+2)=y(d-1)(xn)+∫-20∑i=0k(-1)i(-si)∇ifn+2hds.
Simplify
(46)y(d-1)(xn+2)=y(d-1)(xn)+h∑i=0kγ2,1,i*∇ifn+2,
where
(47)γ2,1,i*=(-1)i∫-20(-si)ds.
Define the generating function G1*(t) for the coefficient γ2,1,i* as follows:
(48)G1*(t)=∑i=0∞γ2,1,i*ti.
or
(49)G1*(t)=∑i=0∞(-t)i∫-20(-si)ds,(50)G1*(t)=∫-20(1-t)-sds,(51)G1*(t)=∫-20e-slog⁡(1-t)ds,
which leads to
(52)G1*(t)=-[1log⁡(1-t)-(1-t)2log⁡(1-t)].
For the case t=2, the approximate solution of y has the form
(53)y(d-2)(xn+2)=y(d-2)(xn)+hy(d-1)(xn)+∫xnxn+2(xn+2-x)(1)1!000000×f(x,y,y′,y′′,…,y(d-1))dx.
The coefficients are given by
(54)γ2,2,i*=(-1)i∫-20(-s)1!(-si)ds,
where γ2,2,i* are the coefficients of the backward difference formulation of (54) which can be represented by
(55)y(d-2)(xn+2)=y(d-2)(xn)+hy(d-1)(xn)+h2∑i=0kγ2,2,i*∇ifn+2.
Define the generating function of the coefficient γ2,2,i* as follows:
(56)G2*(t)=∑i=0∞γ2,2,i*ti.
Substituting (54) into G2*(t) above gives
(57)G2*(t)=∫-20(-s)1!e-slog⁡(1-t)ds.
Solving (57) with the substitution of (51) produces the relationship
(58)G2*(t)=11![2(1-t)2log⁡(1-t)-1!G1*(t)log⁡(1-t)].
By using the same process previously, we note that for integrating (d-1) times yield
(59)G(d-1)*(t)=∫-20(-s)(d-2)(d-2)!e-slog⁡(1-t)ds,(60)G(d-1)*(t)=1(d-2)![2(d-2)(1-t)2log⁡(1-t)-(d-2)!G(d-2)*(t)log⁡(1-t)].
Integrating (d) times yield
(61)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+∫xnxn+2(xn+2-x)(d-1)(d-1)!000000×f(x,y,y′,y′′,…,y(d-1))dx.
The coefficients are given by
(62)γ2,(d),i*=(-1)i∫-20(-s)(d-1)(d-1)!(-si)ds,
where γ2,(d),i* are the coefficients of the backward difference formulation of (62) which can be represented by
(63)y(xn+2)=y(xn)+hy′(xn)+⋯+h(d-1)(d-1)!y(d-1)(xn)+h(d)∑i=0kγ2,(d),i*∇ifn+2.
Define the generating function G(d)*(t) of the coefficient γ2,(d),i* as follows:
(64)G(d)*(t)=∑i=0∞γ2,(d),i*ti.
Substituting (62) into G(d)*(t) above gives
(65)G(d)*(t)=∫-20(-s)(d-1)(d-1)!e-slog⁡(1-t)ds.
Solving (65) with the substitution of (59) produces the relationship
(66)G(d)*(t)=1(d-1)![2(d-1)(1-t)2log⁡(1-t)-(d-1)!G(d-1)*(t)log⁡(1-t)].


## 5. The Relationship between the Explicit and Implicit Coefficients

### 5.1. For the First Point

Calculating the integration coefficients directly is time consuming when large numbers of integration are involved. An efficient technique of computing the coefficients is by formulating a recursive relationship between them. With this recursive relationship, we are able to obtain the implicit integration coefficient with minimal effort. The relationship between the explicit and implicit coefficients for the first pointyn+1 is already given by Suleiman et al. [8].


### 5.2. For the Second Point

For first-order coefficients,(67)G1*(t)=-[1log⁡(1-t)-(1-t)2log⁡(1-t)].
It can be written as
(68)G1*(t)=-(1-t)2[1(1-t)2log⁡(1-t)-1log⁡(1-t)].
By substituting
(69)G1(t)=-[1(1-t)2log⁡(1-t)-1log⁡(1-t)]
into (68), we have
(70)G1*(t)=(1-t)2G1(t),(∑i=0∞γ2,1,i*ti)=(1-t)2(∑i=0∞γ2,1,iti).Expanding the equation yields(71)(γ2,1,0*+γ2,1,1*t+γ2,1,2*t2+⋯)=1(1+2t+3t2+⋯)(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯),(γ2,1,0*+γ2,1,1*t+γ2,1,2*t2+⋯)(1+2t+3t2+⋯)=(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯).
This gives the recursive relationship
(72)∑i=0k(k-i+1)γ2,1,i*=γ2,1,k.
For second-order coefficient,
(73)G2*(t)=11![2(1-t)2log⁡(1-t)-1!G1*(t)log⁡(1-t)].
It can be written as
(74)G2*(t)=(1-t)21![2log⁡(1-t)-1!G1*(t)(1-t)2log⁡(1-t)].
Substituting (70) into the equation above gives
(75)G2*(t)=(1-t)21![2log⁡(1-t)-1!(1-t)2G1(t)(1-t)2log⁡(1-t)]
or
(76)G2*(t)=(1-t)21![2log⁡(1-t)-1!G1(t)log⁡(1-t)].
Substituting (25) into (76) gives
(77)G2*(t)=(1-t)2G2(t),(∑i=0∞γ2,2,i*ti)=(1-t)2(∑i=0∞γ2,2,iti).
Expanding the equation yields
(78)(γ2,2,0*+γ2,2,1*t+γ2,2,2*t2+⋯)=1(1+2t+3t2+⋯)(γ2,2,0+γ2,2,1t+γ2,2,2t2+⋯),(γ2,2,0*+γ2,2,1*t+γ2,2,2*t2+⋯)(1+2t+3t2+⋯)=(γ2,2,0+γ2,2,1t+γ2,2,2t2+⋯).
This gives the recursive relationship
(79)∑i=0k(k-i+1)γ2,2,i*=γ2,2,k.
By using the same process previously, we note that, for (d-1)-order coefficient, we have
(80)G(d-1)*(t)=(1-t)2G(d-1)(t),
which leads to a recursive relationship
(81)∑i=0k(k-i+1)γ2,(d-1),i*=γ2,(d-1),k.
For (d)-order coefficient, we have
(82)G(d)*(t)=1(d-1)![2(d-1)(1-t)2log⁡(1-t)-(d-1)!G(d-1)*(t)log⁡(1-t)].
It can be written as
(83)G(d)*(t)=(1-t)2(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)*(t)(1-t)2log⁡(1-t)].
Substituting (80) into (83) gives
(84)G(d)*(t)=(1-t)2(d-1)![2(d-1)log⁡(1-t)-(d-1)!(1-t)2G(d-1)(t)(1-t)2log⁡(1-t)]
or
(85)G(d)*(t)=(1-t)2(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)(t)log⁡(1-t)].
Substituting
(86)G(d)(t)=1(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)(t)log⁡(1-t)]
into (85) leads to
(87)G(d)*(t)=(1-t)2G(d)(t),(∑i=0∞γ2,(d),i*ti)=(1-t)2(∑i=0∞γ2,(d),iti).
Expanding the equation yields
(88)(γ2,(d),0*+γ2,(d),1*t+γ2,(d),2*t2+⋯)=1(1+2t+3t2+⋯)×(γ2,(d),0+γ2,(d),1t+γ2,(d),2t2+⋯),(γ2,(d),0*+γ2,(d),1*t+γ2,(d),2*t2+⋯)(1+2t+3t2+⋯)=(γ2,(d),0+γ2,(d),1t+γ2,(d),2t2+⋯),
which leads to a recursive relationship
(89)∑i=0k(k-i+1)γ2,(d),i*=γ2,(d),k.
Tables 1 and 2 are a few examples of the explicit and implicit integration coefficients.Table 1
The explicit integration coefficients fork from 0 to 6 (for yn+2).
k
0
1
2
3
4
5
6
γ
2,1
,
k
2
2
7/3
8/3
269/90
33/10
13613/3780
γ
2,2
,
k
2
4/3
4/3
62/45
43/30
94/63
1466/945
γ
2,3
,
k
4/3
2/3
3/5
26/45
359/630
179/315
16159/28350Table 2
The implicit integration coefficients fork from 0 to 6 (for yn+2).
k
0
1
2
3
4
5
6
γ
2,1
,
k
*
2
-
2
1/3
0
−1/90
−1/90
−8/945
γ
2,2
,
k
*
2
-
8/3
2/3
4/90
1/90
1/315
1/1890
γ
2,3
,
k
*
4/3
-
2
3/5
2/45
1/70
2/315
47/14175


## 5.1. For the First Point

Calculating the integration coefficients directly is time consuming when large numbers of integration are involved. An efficient technique of computing the coefficients is by formulating a recursive relationship between them. With this recursive relationship, we are able to obtain the implicit integration coefficient with minimal effort. The relationship between the explicit and implicit coefficients for the first pointyn+1 is already given by Suleiman et al. [8].


## 5.2. For the Second Point

For first-order coefficients,(67)G1*(t)=-[1log⁡(1-t)-(1-t)2log⁡(1-t)].
It can be written as
(68)G1*(t)=-(1-t)2[1(1-t)2log⁡(1-t)-1log⁡(1-t)].
By substituting
(69)G1(t)=-[1(1-t)2log⁡(1-t)-1log⁡(1-t)]
into (68), we have
(70)G1*(t)=(1-t)2G1(t),(∑i=0∞γ2,1,i*ti)=(1-t)2(∑i=0∞γ2,1,iti).Expanding the equation yields(71)(γ2,1,0*+γ2,1,1*t+γ2,1,2*t2+⋯)=1(1+2t+3t2+⋯)(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯),(γ2,1,0*+γ2,1,1*t+γ2,1,2*t2+⋯)(1+2t+3t2+⋯)=(γ2,1,0+γ2,1,1t+γ2,1,2t2+⋯).
This gives the recursive relationship
(72)∑i=0k(k-i+1)γ2,1,i*=γ2,1,k.
For second-order coefficient,
(73)G2*(t)=11![2(1-t)2log⁡(1-t)-1!G1*(t)log⁡(1-t)].
It can be written as
(74)G2*(t)=(1-t)21![2log⁡(1-t)-1!G1*(t)(1-t)2log⁡(1-t)].
Substituting (70) into the equation above gives
(75)G2*(t)=(1-t)21![2log⁡(1-t)-1!(1-t)2G1(t)(1-t)2log⁡(1-t)]
or
(76)G2*(t)=(1-t)21![2log⁡(1-t)-1!G1(t)log⁡(1-t)].
Substituting (25) into (76) gives
(77)G2*(t)=(1-t)2G2(t),(∑i=0∞γ2,2,i*ti)=(1-t)2(∑i=0∞γ2,2,iti).
Expanding the equation yields
(78)(γ2,2,0*+γ2,2,1*t+γ2,2,2*t2+⋯)=1(1+2t+3t2+⋯)(γ2,2,0+γ2,2,1t+γ2,2,2t2+⋯),(γ2,2,0*+γ2,2,1*t+γ2,2,2*t2+⋯)(1+2t+3t2+⋯)=(γ2,2,0+γ2,2,1t+γ2,2,2t2+⋯).
This gives the recursive relationship
(79)∑i=0k(k-i+1)γ2,2,i*=γ2,2,k.
By using the same process previously, we note that, for (d-1)-order coefficient, we have
(80)G(d-1)*(t)=(1-t)2G(d-1)(t),
which leads to a recursive relationship
(81)∑i=0k(k-i+1)γ2,(d-1),i*=γ2,(d-1),k.
For (d)-order coefficient, we have
(82)G(d)*(t)=1(d-1)![2(d-1)(1-t)2log⁡(1-t)-(d-1)!G(d-1)*(t)log⁡(1-t)].
It can be written as
(83)G(d)*(t)=(1-t)2(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)*(t)(1-t)2log⁡(1-t)].
Substituting (80) into (83) gives
(84)G(d)*(t)=(1-t)2(d-1)![2(d-1)log⁡(1-t)-(d-1)!(1-t)2G(d-1)(t)(1-t)2log⁡(1-t)]
or
(85)G(d)*(t)=(1-t)2(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)(t)log⁡(1-t)].
Substituting
(86)G(d)(t)=1(d-1)![2(d-1)log⁡(1-t)-(d-1)!G(d-1)(t)log⁡(1-t)]
into (85) leads to
(87)G(d)*(t)=(1-t)2G(d)(t),(∑i=0∞γ2,(d),i*ti)=(1-t)2(∑i=0∞γ2,(d),iti).
Expanding the equation yields
(88)(γ2,(d),0*+γ2,(d),1*t+γ2,(d),2*t2+⋯)=1(1+2t+3t2+⋯)×(γ2,(d),0+γ2,(d),1t+γ2,(d),2t2+⋯),(γ2,(d),0*+γ2,(d),1*t+γ2,(d),2*t2+⋯)(1+2t+3t2+⋯)=(γ2,(d),0+γ2,(d),1t+γ2,(d),2t2+⋯),
which leads to a recursive relationship
(89)∑i=0k(k-i+1)γ2,(d),i*=γ2,(d),k.
Tables 1 and 2 are a few examples of the explicit and implicit integration coefficients.Table 1
The explicit integration coefficients fork from 0 to 6 (for yn+2).
k
0
1
2
3
4
5
6
γ
2,1
,
k
2
2
7/3
8/3
269/90
33/10
13613/3780
γ
2,2
,
k
2
4/3
4/3
62/45
43/30
94/63
1466/945
γ
2,3
,
k
4/3
2/3
3/5
26/45
359/630
179/315
16159/28350Table 2
The implicit integration coefficients fork from 0 to 6 (for yn+2).
k
0
1
2
3
4
5
6
γ
2,1
,
k
*
2
-
2
1/3
0
−1/90
−1/90
−8/945
γ
2,2
,
k
*
2
-
8/3
2/3
4/90
1/90
1/315
1/1890
γ
2,3
,
k
*
4/3
-
2
3/5
2/45
1/70
2/315
47/14175


## 6. Problem Tested

The problems shown in Table3 are used to test the performance of the method.Table 3
List of test problems.
Problem
Initial value
Interval
1
y
′
=
y Exact solution: y(x)=ex Source: artificial problem
y
(
0
)
=
1
0
≤
x
≤
20
2
y
′
′
=
-
2
y
′
+
3
y Exact solution: y(x)=ex+e-3x Source: Suleiman [9]
y
(
0
)
=
2 y′(0)=-2
0
≤
x
≤
64
3
y
1
′
′
=
-
y
1
h
h
g
h
g
r
3 y2′′=-y2jhhbfr3 r=(y12+y22)2 Exact solution:  y1(x)=cos⁡x  y2(x)=sin⁡⁡x Source: Shampine and Gordon [10]
y
1
(
0
)
=
1 y1′(0)=0 y2(0)=0 y2′(0)=1
0
≤
x
≤
16
π
4
y
1
′
′
=
-
2
y
1
′
-
5
y
2
+
3 y2′=y1′+2y2 Exact solution:  y1(x)=2cos⁡x+6sin⁡⁡x-6x-2  y1(x)=2sin⁡⁡x-2cos⁡x+3 Source: Suleiman [9]
y
1
(
0
)
=
0 y1′(0)=0 y2(0)=1
0
≤
x
≤
16
π
5
y
′
′
′
=
-
1
x
y
′
′
+
1
x
2
y
′
+
1
j
g
h
f
x  Exact solution: y(x)=x28(2ln⁡(x2)-(3313)-23ln⁡(2))+(13-2621ln⁡(x2))ln⁡(2)+3326 Source: Russell and Shampine [11]
y
(
1
)
=
26
21
ln
2
⁡
(
2
)
+
99
104 y′(1)=-4021ln⁡(2)⁡-513 y′′(1)=326+47ln⁡(2)
1
≤
x
≤
50


## 7. Numerical Result

Tables4, 5, 6, 7, and 8 give the numerical results for problems given in the previous section. The results for the 2PBBD are compared with those of 2PBDD and 1PBD according to Omar [5] and Suleiman et al. [8], respectively. Also given are graphs, where Log10⁡(MAXE) is plotted against Log10⁡(H) and Log10⁡(Time). The following notations are used in the tables:H: step size,2PBBD: 2-point block backward difference method,2PBDD: 2-point block divided difference method,1PBD: 1-point backward difference method,NS: total number of steps,MAXE: maximum error,TIME: total execution times (in microsecond).Two sets of scaled graphs were plotted, namely, (i) Log10⁡(MAXE) against Log10⁡(H) and (ii) Log10⁡(MAXE) against Log10⁡(TIME). For a particular abscissa, the lowest value of the ordinate is considered to be the more efficient at the abscissa considered. Hence, for the first set of graphs, that is, Log10⁡(MAXE) against Log10⁡(H), the method 2PBBD is better when Log10⁡(H)<-5, and loses out for value of Log10⁡(H)>-5. For the second set of graphs, as the time increases, the 2PBBD is the method of choice since it is lowest for all five sets of problems (see Figures 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11). It gives us the impression of stability, where the errors grow most slowly compared with the other methods, 2PBDD and 1PBD.Table 4
Numerical result for Problem 1.
H
Method
NS
Log
10
⁡
(MAXE)
Time
10
-
2
2PBBD
1000
-
5.40549
6914
2PBDD
1000
-
5.88584
9509
1PBD
2000
-
6.77769
6912
10
-
3
2PBBD
10000
-
7.34272
49146
2PBDD
10000
-
8.87606
54025
1PBD
20000
-
9.77432
50095
10
-
4
2PBBD
100000
-
9.40044
195037
2PBDD
100000
-
10.67248
256339
1PBD
200000
-
10.48025
197800
10
-
5
2PBBD
1000000
-
9.37408
1759047
2PBDD
1000000
-
9.37799
2055758
1PBD
2000000
-
9.23059
1638426
10
-
6
2PBBD
10000000
-
8.81572
17500121
2PBDD
10000000
-
8.81572
19917218
1PBD
20000000
-
8.28209
14653209
10
-
7
2PBBD
100000000
-
7.78104
176146062
2PBDD
100000000
-
7.78104
199668362
1PBD
200000000
-
7.39425
145626635Table 5
Numerical result for Problem 2.
H
Method
NS
Log
10
⁡
MAXE
Time
10
-
2
2PBBD
3200
-
4.87654
21516
2PBDD
3200
-
4.58643
23350
1PBD
6400
-
5.83463
22890
10
-
3
2PBBD
32000
-
7.86398
86824
2PBDD
32000
-
7.56556
111374
1PBD
64000
-
8.80976
87426
10
-
4
2PBBD
320000
-
9.96338
776833
2PBDD
320000
-
9.82447
828900
1PBD
640000
-
8.99888
725855
10
-
5
2PBBD
3200000
-
8.68100
7769459
2PBDD
3200000
-
8.68099
8208180
1PBD
6400000
-
8.04024
6474481
10
-
6
2PBBD
32000000
-
7.64880
77008100
2PBDD
32000000
-
7.64879
81918617
1PBD
64000000
-
7.15759
63364050
10
-
7
2PBBD
320000000
-
6.53269
773213500
2PBDD
320000000
-
6.53269
817874121
1PBD
640000000
-
6.28394
616850084Table 6
Numerical result for Problem 3.
H
Method
NS
Log
10
⁡
MAXE
Time
10
-
2
2PBBD
2513
-
6.17611
24465
2PBDD
2513
-
5.87541
24221
1PBD
5026
-
7.07922
28522
10
-
3
2PBBD
25133
-
9.17248
152099
2PBDD
25133
-
8.87081
159389
1PBD
50265
-
10.01264
139073
10
-
4
2PBBD
251328
-
10.04346
1418931
2PBDD
251328
-
10.0434
1366837
1PBD
502655
-
9.26354
1219971
10
-
5
2PBBD
2513274
-
8.80560
14315026
2PBDD
2513274
-
8.80560
13192255
1PBD
5026548
-
8.32281
11525418
10
-
6
2PBBD
25132742
-
7.99962
141938590
2PBDD
25132742
-
7.99962
130102816
1PBD
50265482
-
7.45728
106896913
10
-
7
2PBBD
251327412
-
6.87690
1412455700
2PBDD
251327412
-
6.87690
1300828490
1PBD
502654824
-
6.44359
1063913703Table 7
Numerical result for Problem 4.
H
Method
NS
Log
10
⁡
MAXE
Time
10
-
2
2PBBD
2513
-
5.29159
23909
2PBDD
2513
-
5.69715
26684
1PBD
5026
-
6.58585
31772
10
-
3
2PBBD
25133
-
7.74946
121769
2PBDD
25133
-
8.68299
131297
1PBD
50265
-
9.58464
121536
10
-
4
2PBBD
251328
-
9.39201
1102741
2PBDD
251328
-
10.04970
1068239
1PBD
502655
-
9.26367
1046930
10
-
5
2PBBD
2513274
-
8.80463
10004344
2PBDD
2513274
-
8.80548
9407995
1PBD
5026548
-
8.32278
10087395
10
-
6
2PBBD
25132742
-
7.99900
99147571
2PBDD
25132742
-
7.99900
91521335
1PBD
50265482
-
7.45728
93207026
10
-
7
2PBBD
251327412
-
6.87692
988919951
2PBDD
251327412
-
6.87692
911824508
1PBD
502654824
-
6.44359
922126253Table 8
Numerical result for Problem 5.
H
Method
NS
Log
10
⁡
MAXE
Time
10
-
2
2PBBD
2450
-
3.65914
25345
2PBDD
2450
-
3.37200
26172
1PBD
4900
-
5.08070
27139
10
-
3
2PBBD
24500
-
6.62483
106625
2PBDD
24500
-
6.32538
116848
1PBD
49000
-
8.06583
108650
10
-
4
2PBBD
245000
-
9.64707
1021740
2PBDD
245000
-
9.30862
996446
1PBD
490000
-
10.55259
888813
10
-
5
2PBBD
2450000
-
10.06216
9797882
2PBDD
2450000
-
10.06228
8683164
1PBD
4900000
-
9.42795
8211650
10
-
6
2PBBD
24500000
-
8.99818
97582403
2PBDD
24500000
-
8.99818
85958740
1PBD
49000000
-
8.73189
80890328
10
-
7
2PBBD
245000000
-
8.10589
978631746
2PBDD
245000000
-
8.10589
844019109
1PBD
490000000
-
7.70882
790057441Figure 2
Graph ofLog10⁡(MAXE) plotted against Log10⁡(H) for Problem 1.Figure 3
Graph ofLog10⁡(MAXE) plotted against Log10⁡(TIME) for Problem 1.Figure 4
Graph ofLog10⁡(MAXE) plotted against Log10⁡(H) for Problem 2.Figure 5
Graph ofLog10⁡(MAXE) plotted against Log10⁡(TIME) for Problem 2.Figure 6
Graph ofLog10⁡(MAXE) plotted against Log10⁡(H) for Problem 3.Figure 7
Graph ofLog10⁡(MAXE) plotted against Log10⁡(TIME) for Problem 3.Figure 8
Graph ofLog10⁡(MAXE) plotted against Log10⁡(H) for Problem 4.Figure 9
Graph ofLog10⁡(MAXE) plotted against Log10⁡(TIME) for Problem 4.Figure 10
Graph ofLog10⁡(MAXE) plotted against Log10⁡(H) for Problem 5.Figure 11
Graph ofLog10⁡(MAXE) plotted against Log10⁡(TIME) for Problem 5.


## 8. Conclusion

Of the 3 methods, 2PBBD is therefore preferred as a general code and should be included as a collection of methods, as a code for parallelization purposes, as an assembly of codes to be tested and studied, and as a code for solving ODEs.



---

*Source: 867095-2014-09-17.xml*