# Solving Nonstiff Higher Order Odes Using Variable Order Step Size Backward Difference Directly

Authors: Ahmad Fadly Nurullah Rasedee; Mohamed bin Suleiman; Zarina Bibi Ibrahim
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The current numerical techniques for solving a system of higher order ordinary differential equations (ODEs) directly calculate the integration coefficients at every step. Here, we propose a method to solve higher order ODEs directly by calculating the integration coefficients only once at the beginning of the integration and if required once more at the end. The formulae will be derived in terms of backward difference in a constant step size formulation. The method developed will be validated by solving some higher order ODEs directly using variable order step size. To simplify the evaluations of the integration coefficients, we find the relationship between various orders. The results presented confirmed our hypothesis.


---

## Body

## 1. Introduction

In this paper, we will focus only on nonstiff ODEs of the form(1)
y
(
d
)
=
f
(
x
,
Y
~
)
given
(2)
Y
~
=
(
y
,
y
′
,
y
′′
,
…
,
y
(
d
-
1
)
)
,
η
~
=
(
η
,
η
′
,
η
′′
,
…
,
η
(
d
-
1
)
)
,
where Y
~
(
a
)
=
η
~ in the interval a
≤
x
≤
b, d is the order of the ODE. We impose the condition that Y
~ is continuous. This implies that solution for Y
~ exists.Many science and engineering problems are in the form of higher order initial value problems (IVPs) ODEs. Authors such as Gear [1], Suleiman [2, 3], Hall and Suleiman [4], and Omar [5] suggested a new approach for solving higher order ODEs directly. An algorithm was designed by Suleiman [2] to solve stiff and nonstiff higher order ODEs directly without reducing the order of the problems to first order. He called it the direct integration method. The drawbacks to methods described by Omar [5], Suleiman [3], Gear [1], and Lambert [6] are the tedious calculations of the divided differences and recurrence relations in computing the integration coefficients.Here, a new and efficient algorithm for solving IVPs of higher order ODEs directly using variable order step size method in its backward difference formulation (MSBD method) is developed. In contrast to the integration coefficients used in Suleiman’s [3] direct integration (DI) method, the code MSBD calculates the integration coefficients once at the start and once when calculating the last point. The simpler nature of the backward difference compared to the divided difference gives an added advantage to the MSBD over the DI method when formulating and coding the method. This produces more elegant error formulae. The DI method in Suleiman [2] has been proven to converge. We note that the DI method is formulated based on divided difference while the proposed MSBD technique is based on backward difference. Therefore, in a similar way, we can prove the convergence of the MSBD. A brief explanation of the DI code is given below.LetP
k
,
n
(
x
) be the integrating polynomial with degree k
-
1; interpolating k points is denoted by
(3)
P
k
,
n
(
x
)
=
f
n
+
(
x
-
x
n
)
f
[
n
,
n
-
1
]
+
⋯
+
(
x
-
x
n
)
⋯
(
x
-
x
n
-
k
+
2
)
f
[
n
,
n
-
1
,
…
,
n
-
k
+
1
]
,
and obtained through divided difference.We defineg
i
,
t
,
t
>
0 to be the t-fold integral and the integration coefficients obtained are denoted as follows:
(4)
g
i
,
t
=
(
x
n
+
1
-
x
n
-
i
+
1
)
g
i
-
1
,
t
-
t
g
i
-
1
,
t
+
1
.
For the case i
=
0,
(5)
g
0
,
t
=
h
n
+
1
t
t
!
,
where h
r
=
x
r
-
x
r
-
1
,
r
=
1,2
,
3
,
….Lete be given explicitly by
(6)
e
=
f
(
x
n
+
1
,
P
_
n
+
1
)
-
P
n
+
1
(
d
)
,
P
n
+
1
(
d
-
t
)
,
t
=
1,2
,
…
,
d
,
which are the predicted derivatives.Using theg
i
,
t coefficients, the predictor and corrector are given by
(7)
P
n
+
1
(
r
)
=
∑
i
=
0
d
-
1
-
r
h
i
i
!
y
n
(
r
+
i
)
+
∑
i
=
0
k
-
1
g
i
,
d
-
r
f
[
n
,
n
-
1
,
…
,
n
-
i
]
y
n
+
1
(
d
-
t
)
=
P
n
+
1
(
d
-
t
)
+
g
k
,
t
g
k
,
0
e
,
t
=
1,2
,
…
,
d
.


## 2. Derivation ofdth Order Explicit and Implicit Backward Differences Method

In this section, the integration coefficients of the backward difference formulae are derived in order to obtain a relationship between the explicit and implicit formulae. This will make the computation easier.Integrating (1) once yields
(8)
y
(
x
n
+
1
)
=
y
(
x
n
)
+
∫
x
n
x
n
+
1
f
(
x
,
y
,
y
′
,
y
′′
)
d
x
.
Let P
n
(
x
) be the interpolating polynomial which interpolates the k values (
x
n
,
f
n
)
,
(
x
n
-
1
,
f
n
-
1
)
,
…
,
(
x
n
-
k
+
1
,
f
n
-
k
+
1
); then
(9)
P
n
(
x
)
=
∑
i
=
0
k
-
1
(
-
1
)
i
(
-
s
i
)
∇
i
f
n
.
Next, approximating f in (8) with P
n
(
x
) and letting
(10)
x
=
x
n
+
s
h
give us
(11)
y
(
x
n
+
1
)
=
y
(
x
n
)
+
∫
0
1
∑
i
=
0
k
-
1
(
-
1
)
i
(
-
s
i
)
∇
i
f
n
d
s
,
where
(12)
y
1
,
i
=
(
-
1
)
i
∫
0
1
(
-
s
i
)
d
s
.
Let the generating function G
1
(
t
) for the coefficients γ
1
,
i be defined as follows:
(13)
G
1
(
t
)
=
∑
i
=
0
∞
γ
1
,
i
t
i
.
Substituting (12) in G
1
(
t
) gives
(14)
G
1
(
t
)
=
∫
0
1
e
-
s
log
⁡
(
1
-
t
)
d
s
which leads to
(15)
G
1
(
t
)
=
-
[
(
1
-
t
)
-
1
log
⁡
(
1
-
t
)
-
1
log
⁡
(
1
-
t
)
]
.
Hence, the coefficients of the backward difference formulation, γ
1
,
k, are given by
(16)
∑
i
=
0
k
(
γ
1
,
i
k
-
i
+
1
)
=
1
.
γ
1
,
k
=
1
-
∑
i
=
0
k
-
1
(
γ
1
,
i
k
-
i
+
1
)
,
k
=
1,2
,
…
,
γ
1,0
=
1
.
To derive the dth order generating function, we integrate (1) d number of times and mathematical induction gives
(17)
G
(
d
)
(
t
)
=
1
(
d
-
1
)
!
[
1
(
d
-
1
)
log
⁡
(
1
-
t
)
-
(
d
-
1
)
!
G
(
d
-
1
)
(
t
)
log
⁡
(
1
-
t
)
]
.
This gives the generalized explicit coefficients which are denoted by
(18)
γ
(
d
)
,
0
=
γ
(
d
-
1
)
,
1
γ
(
d
)
,
k
=
γ
(
d
-
1
)
,
k
+
1
-
∑
i
=
0
k
-
1
γ
(
d
)
,
i
k
-
i
+
1
k
=
0,1
,
2
…
.
Similarly, deriving the generalized implicit formulae gives the following relationship for the generating functions:
(19)
G
(
d
)
*
(
t
)
=
1
(
d
-
1
)
!
[
(
1
-
t
)
log
⁡
(
1
-
t
)
-
(
d
-
1
)
!
G
(
d
-
1
)
*
(
t
)
log
⁡
(
1
-
t
)
]
,
where the coefficients are given by
(20)
γ
(
d
)
,
0
*
=
∑
i
=
0
1
γ
(
d
-
1
)
,
i
*
γ
(
d
)
,
k
*
=
∑
i
=
0
k
+
1
γ
(
d
-
1
)
,
i
*
-
∑
i
=
0
k
-
1
γ
(
d
)
,
i
*
I
1
,
k
+
1
-
i
,
k
=
1,2
,
3
…
,
and I
1
,
k
+
1
-
i is the Lagrange coefficient.


## 3. The Relationship between the Explicit and Implicit Integration Coefficients

Calculating the integration coefficients directly is time consuming if large numbers of integrations are involved. By obtaining a recursive relationship between the coefficients, we are able to obtain the implicit integration coefficient more efficiently. The relationship between the explicit and implicit coefficients is discussed below.For first order coefficients,(21)
G
1
*
(
t
)
=
-
[
1
log
⁡
⁡
(
1
-
t
)
-
1
-
t
log
⁡
⁡
(
1
-
t
)
]
.
It can be written as
(22)
G
1
*
(
t
)
=
-
(
1
-
t
)
[
1
(
1
-
t
)
log
⁡
⁡
(
1
-
t
)
-
1
log
⁡
⁡
(
1
-
t
)
]
.
By substituting
(23)
G
1
(
t
)
=
1
(
1
-
t
)
log
⁡
⁡
(
1
-
t
)
-
1
log
⁡
⁡
(
1
-
t
)
into (22), we have
(24)
G
1
*
(
t
)
=
(
1
-
t
)
G
1
(
t
)
(25)
(
∑
i
=
0
∞
γ
1
,
i
*
t
i
)
=
(
1
-
t
)
(
∑
i
=
0
∞
γ
1
,
i
t
i
)
.
Solving the equation above gives the recursive relationship
(26)
∑
i
=
0
k
γ
1
,
i
*
=
γ
1
,
k
.
For second order coefficient,
(27)
G
2
*
(
t
)
=
-
1
1
!
[
1
log
⁡
⁡
(
1
-
t
)
-
1
!
G
1
*
(
t
)
log
⁡
⁡
(
1
-
t
)
]
.
It can be written as
(28)
G
2
*
(
t
)
=
(
1
-
t
)
1
!
[
1
log
⁡
⁡
(
1
-
t
)
-
1
!
G
1
*
(
t
)
(
1
-
t
)
log
⁡
⁡
(
1
-
t
)
]
.
Substituting (24) into the equation above gives
(29)
G
2
*
(
t
)
=
(
1
-
t
)
1
!
[
1
log
⁡
⁡
(
1
-
t
)
-
1
!
(
1
-
t
)
G
1
(
t
)
(
1
-
t
)
log
⁡
⁡
(
1
-
t
)
]
or
(30)
G
2
*
(
t
)
=
(
1
-
t
)
1
!
[
1
log
⁡
⁡
(
1
-
t
)
-
1
!
G
1
(
t
)
log
⁡
⁡
(
1
-
t
)
]
which can be simplified as
(31)
G
2
*
(
t
)
=
(
1
-
t
)
G
2
(
t
)
,
(
∑
i
=
0
∞
γ
2
,
i
*
t
i
)
=
(
1
-
t
)
(
∑
i
=
0
∞
γ
2
,
i
t
i
)
.
Solving the equation above gives the recursive relationship
(32)
∑
i
=
0
k
γ
2
,
i
*
=
γ
2
,
k
.
For dth order coefficient, using mathematical induction, we have
(33)
G
(
d
)
*
(
t
)
=
(
1
-
t
)
G
(
d
)
(
t
)
(
∑
i
=
0
∞
γ
(
d
)
,
i
*
t
i
)
=
(
1
-
t
)
(
∑
i
=
0
∞
γ
(
d
)
,
i
t
i
)
.
In similar manner to the case of the second order coefficients, we obtain the general dth order coefficients as follows:
(34)
∑
i
=
0
k
γ
(
d
)
,
i
*
=
γ
(
d
)
,
k
.


## 4. Order and Step Size Selection

An important decision to be made is whether to accept the result of an integration step. Although the efficiency of an algorithm is affected by the strategy for selecting the order and step size, its reliability is affected by the acceptance criteria used.The literature on variable order step size codes show numerous strategies of order and step size selection. Varying the order in a multistep method is very easy. The order depends directly on the back values stored. After the completion of any step, the order can be increased by one if none of the back values used in the previous step are discarded. The order can also be decreased to any value desired by discarding the appropriate number of back values. Through experience, it is suggested that the order strategies which are not biased to selecting a lower order when implemented in an Adams code are efficient for nonstiff problems. After considering certain criteria, we adopted the order strategies similar to those used by Shampine and Gordon [7].Leth be the calculated step size and h
new the final step size. See Suleiman [3]. In practice, we multiply h by a safety factor R such that h
new
=
R
h, in order to give a more conservative estimate for h
new and hence reduce the numbers of steps rejected. For our code, we take the safety factor as 0.8. Shampine and Gordon [7] discuss a practical approach to the result of the convergence and stability when using a variable step size. Their experience suggests that their ratio of successive step sizes will need to be restricted in order to ensure stability. For this reason and reasons of economy discussed in the section of generating the algorithm, it would appear that the use of a constant step size is desirable. Opposing this view, however, is the fact that the use of the maximum possible step size minimizes the number of steps to be taken and hence the number of derivative evaluations.


## 5. Error Estimation

Adopting the same strategies from Hall and Watt [8] for higher order ODEs, the estimations for the local errors on each step of the integration are shown below.The predictor takes the form of(35)
y
pr
n
+
1
(
d
)
=
∑
i
=
0
k
-
1
γ
d
,
i
∇
i
f
n
y
pr
n
+
1
(
d
-
1
)
=
y
n
(
d
-
1
)
+
h
∑
i
=
0
k
-
1
γ
d
,
i
∇
i
f
n
⋮
y
pr
n
+
1
=
∑
i
=
0
d
-
1
h
i
i
!
y
n
i
+
h
d
∑
i
=
0
k
-
1
γ
d
,
i
∇
i
f
n
,
where the constant γ
d
,
i is independent of k. In P
k
E
C
k
+
1
E algorithm, the corrector can be written as follows:
(36)
y
n
+
1
(
d
)
=
∑
i
=
0
k
-
1
γ
d
,
i
*
∇
pr
i
f
n
+
1
y
n
+
1
(
d
-
1
)
=
y
n
(
d
-
1
)
+
h
∑
i
=
0
k
-
1
γ
d
,
i
*
∇
pr
i
f
n
+
1
⋮
y
n
+
1
=
∑
i
=
0
d
-
1
h
i
i
!
y
n
i
+
h
d
∑
i
=
0
k
-
1
γ
d
,
i
*
∇
pr
i
f
n
+
1
,
where ∇
pr
i denotes the ith backward difference using f
(
x
n
+
1
,
Y
~
n
+
1
pr
) for f
n
+
1. It can be shown that
(37)
γ
1,0
*
=
γ
1,0
=
1
,
γ
d
,
0
*
=
γ
d
,
0
,
∑
i
=
0
k
γ
d
,
i
*
=
γ
d
,
k
,
where (36), can be simplified for computation to
(38)
y
n
+
1
(
d
)
=
y
pr
n
+
1
(
d
)
+
γ
0
,
k
∇
pr
k
f
n
+
1
y
n
+
1
(
d
-
1
)
=
y
pr
n
+
1
(
d
-
1
)
+
h
γ
1
,
k
∇
pr
k
f
n
+
1
⋮
y
n
+
1
=
y
pr
n
+
1
+
h
d
γ
d
,
k
∇
pr
k
f
n
+
1
.The Milne error estimate mentions that the local truncation error (LTE) is simply the difference between the two possibilities ofk and k
-
1 giving
(39)
E
~
k
=
γ
0
,
k
*
∇
pr
k
f
n
+
1
E
~
k
(
1
)
=
h
γ
1
,
k
*
∇
pr
k
f
n
+
1
⋮
E
~
k
(
d
)
=
h
d
γ
d
,
k
*
∇
pr
k
f
n
+
1
.
The selection of an appropriate p, for E
~
k
(
d
-
p
), to control the order and step size is as mentioned in Suleiman [3].When reducing the order, the estimated error would be(40)
E
k
-
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
-
1
*
×
[
∇
pr
k
-
1
f
n
+
1
+
f
(
x
n
+
1
,
Y
~
n
+
1
pr
-
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
-
1
∇
k
-
1
f
n
)
-
f
(
x
n
+
1
,
Y
~
n
+
1
pr
)
]
.
Since a different predictor, of order one less, would have been used, from the mean value theorem, we are able to write
(41)
E
k
-
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
-
1
*
×
[
∇
p
r
k
-
1
f
n
+
1
-
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
-
1
∇
k
-
1
f
n
·
∂
f
∂
y
]
.
In practice, to avoid the extra derivative evaluation required for (40), the decision to reduce the order is based on the estimate
(42)
E
~
k
-
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
-
1
*
∇
pr
k
-
1
f
n
+
1
,
which is immediately available.To consider the increase in order, the estimated error is(43)
E
k
+
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
+
1
*
∇
pr
k
+
1
f
(
x
n
+
1
,
Y
~
n
+
1
pr
+
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
∇
k
f
n
)
,
where the predicted difference is based on a predictor with an extra term included. In estimating E
k
+
1
(
d
-
p
), we obtain
(44)
E
k
+
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
+
1
*
×
[
∇
pr
k
+
1
f
n
+
1
+
f
(
x
n
+
1
,
Y
~
n
+
1
pr
+
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
∇
k
f
n
)
-
f
(
x
n
+
1
,
Y
~
n
+
1
pr
+
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
∇
k
f
n
+
1
)
]
leading to
(45)
E
k
+
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
+
1
*
×
[
∇
pr
k
+
1
f
n
+
1
-
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
∇
k
+
1
f
n
+
1
·
∂
f
∂
y
]
,
which establishes the asymptotic validity of using
(46)
E
~
k
+
1
(
d
-
p
)
=
h
(
d
-
p
)
γ
(
d
-
p
)
,
k
+
1
*
∇
k
+
1
f
n
+
1
.


## 6. Changing the Step Size

We implement the step size changing technique as mentioned in Lambert [6]. Implementing the Adams-Bashforth methods as predictors and Adams-Moulton methods as correctors or which is commonly known as ABM methods in PECE mode in backwards difference form, we then adopt the algorithm for doubling or halving the step size from Krogh [9] which was derived for the ABM method.


## 7. Last Step for Solving Higher Order ODEs

As we mentioned previously, the method we proposed calculates the integrating coefficients only once. In the case of the last step, generally, the final step size will not always fit into our current step size strategy ofh
,
2
h, and (
1
/
2
)
h. Here, we need to calculate the coefficients once more for the last step when the step size is in the form of r
h, r
>
0. The following are the generating functions for the explicit coefficients:
(47)
G
(
d
)
(
t
)
=
1
(
d
-
1
)
!
[
(
r
)
(
d
-
1
)
log
⁡
(
1
-
t
)
-
(
d
-
1
)
!
G
(
d
-
1
)
(
t
)
log
⁡
(
1
-
t
)
]
and implicit coefficients:
(48)
G
(
d
)
*
=
1
(
d
-
1
)
!
[
(
1
-
t
)
r
(
r
)
(
d
-
1
)
log
⁡
(
1
-
t
)
-
(
d
-
1
)
!
G
(
d
-
1
)
*
(
t
)
log
⁡
(
1
-
t
)
]
.
The following yields the relationship between the explicit and implicit coefficients:
(49)
G
(
d
)
*
(
t
)
=
(
1
-
t
)
r
G
(
d
)
(
t
)
(50)
γ
(
d
)
,
k
*
=
γ
(
d
)
,
k
+
∑
i
=
0
k
-
1
(
γ
(
d
)
,
k
-
(
1
+
i
)
(
-
1
)
i
+
1
(
i
+
1
)
!
∏
j
=
0
i
(
-
j
+
r
)
)
.Tables 1 and 2 contain values for the coefficients when r
=
8
/
10.Table 1
r
=
0.8 explicit integration coefficients.
k
0
1
2
3
4
5
6
γ
1
,
k
4
5
8
25
92
375
392
1875
26234
140625
13344
78125
11736472
73828125
γ
2
,
k
8
25
32
375
112
1875
6784
140625
9712
234375
906544
24609375
12355736
369140625
γ
3
,
k
32
375
32
1875
176
15625
6176
703125
181064
24609375
112736
17578125
31750448
5537109375
γ
4
,
k
32
1875
128
46875
1216
703125
32384
24609375
399776
369140625
5155328
5537109375
113924096
138427734375
γ
5
,
k
128
46875
256
703125
5504
24609375
61696
369140625
751808
5537109375
2284928
19775390625
2313234752
22840576171875Table 2
r
=
0.8 implicit integration coefficients.
k
0
1
2
3
4
5
6
γ
1
,
k
*
4
5
-
8
25
-
28
375
24
625
-
3466
140625
-
12416
703125
-
997928
73828125
γ
2
,
k
*
8
25
-
64
375
-
64
1875
-
2336
140625
-
1448
140625
-
176944
24609375
-
1990648
369140625
γ
3
,
k
*
32
375
-
32
625
-
144
15625
-
608
140625
-
1432
546875
-
221152
123046875
-
7362256
5537109375
γ
4
,
k
*
32
1875
-
512
46875
-
256
140625
-
4096
4921875
-
36608
73828125
-
265984
791015625
-
34103936
138427734375
γ
5
,
k
*
128
46875
-
256
140625
-
1408
4921875
-
9472
73828125
-
83648
1107421875
-
7028096
138427734375
-
168729536
4568115234375


## 8. The MSBD Algorithm

The previous write-up was devoted to the calculation of the integration coefficients for the explicit formulae which can form the basis of the predicted values and also the implicit coefficients which are obtained from the explicit formulae.For the sake of clarity, we present the algorithm for the MSBD.Step
1. The integration coefficients are calculated from the algorithm in (16), (18), and (34).Step
2. Use the k back values to obtain the predictor in (30).Step
3. Calculate the corrected values in terms of the predicted values given in (35).Step
4. The errors E
~
k
-
1
(
d
-
p
)
,
E
k
(
d
-
p
), and E
~
k
+
1
(
d
-
p
) are obtained in (39), (42), and (46).Step
5. Determine whether E
k satisfies the local accuracy requirements which we take to be θ
n
+
1
pr
|
E
k
(
d
-
p
)
|
< TOL, where θ
n
+
1
pr
=
1
/
(
A
+
B
+
P
n
) with A
=
1
,
B
=
0 give the absolute error test; A
=
0
,
B
=
1 which gives a relative error test, while A
=
B
=
1 give a mixed error test.Step
6. The order is lowered by one if, for k
>
2, max
⁡
(
|
E
~
k
-
1
(
d
-
p
)
|
,
|
E
~
k
-
2
(
d
-
p
)
|
)
≤
|
E
k
(
d
-
p
)
| and k
=
2
,
|
E
~
k
-
1
(
d
-
p
)
|
≤
0.5
|
E
k
(
d
-
p
)
|. If we have E
~
k
+
1
(
d
-
p
) available, we lower the order if k
>
1 and |
E
~
k
-
1
(
d
-
p
)
|
≤
min
⁡
(
|
E
k
(
d
-
p
)
|
,
|
E
~
k
+
1
(
d
-
p
)
|
). The order is raised by one only after k
+
1 successful step at constant step size such that, for k
>
1
,
|
E
~
k
+
1
(
d
-
p
)
|
<
|
E
k
(
d
-
p
)
|
<
max
⁡
(
|
E
~
k
-
1
(
d
-
p
)
|
,
|
E
~
k
-
2
(
d
-
p
)
|
) and, for k
=
1
,
|
E
~
k
+
1
(
d
-
p
)
|
<
0.5
|
E
k
(
d
-
p
)
|. We restrict k in the range 1
≤
k
≤
12.Step
7. Step size selection is also determined based on the local accuracy requirements. The step size algorithm for doubling ∇
(
D
)
i
I
(
x
)
≔
∇
(
D
)
i
-
1
I
(
x
)
-
∇
(
D
)
i
-
1
I
(
x
-
2
h
) and halving ∇
(
H
)
i
I
(
x
)
≔
∇
(
H
)
i
-
1
I
(
x
)
-
∇
(
H
)
i
-
1
I
(
x
-
h
/
2
)
,
i
=
2,3
,
…
,
k
-
1 is derived by Krogh [9].Step
8. If h
<
|
x
end
-
x
|, where h is the current step size and x
end denotes the end of the interval, repeat Steps 2–7. If not, determine r where r
=
|
x
end
-
x
|
/
h. We then calculate the new integration coefficients as given in (47), (48), and (50). Using the new coefficients, repeat Steps 2–7 and exit the program.


## 9. Test Problems and Numerical Results

### 9.1. Numerical Results

Tables3, 4, 5, and 6 show the numerical results for Problems 1–4, when solved with direct integration, backwards difference, and also reducing to first order systems. Problems 1 and 2 are well-behaved problems whereas Problems 3 and 4 are without exact solutions. For the first two problems, we evaluate the maximum and average values of the error in the computed solution y. The definition of the error is as defined above. The following notation will indicateTIME: the execution time taken in microseconds,
FS: the number of failed steps,
STEPS: total steps,
MTD: the name of the method,
SUCSTEP: the number of successful step,
MAXE: the maximum error,
AVER: the average error,
TOL: the tolerance used,
DI: direct integration,
D1: direct integration with the reduction to first order system,
1PBDVSO: backward difference,
B1: backward difference with the reduction to first order system.The errors calculated are defined as
(51)
(
e
i
)
t
=
|
(
y
i
)
t
-
(
y
(
x
i
)
)
t
A
+
B
(
y
(
x
i
)
)
t
|
with (
y
)
t as the tth component of y
~. A
=
1
,
B
=
0 correspond to the absolute error test, A
=
1
,
B
=
1 correspond to the mixed error test, and A
=
0
,
B
=
1 correspond to the relative error test. The mixed error test is used for all test problems. The maximum and average errors are given by
(52)
MAXE
=
max
⁡
1
<
i
<
SUCSTEP
(
max
⁡
1
<
i
<
N
⁡
(
e
i
)
t
)
,
AVER
=
∑
i
=
1
SUCSTEP
∑
t
=
1
N
(
e
i
)
t
(
N
)
(
SUCSTEP
)
,
where N is the number of equations in the systems.Table 3
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem1.
TOL
MTD
STEPS
FS
MAXE
AVER
TIME (total)
10−2
D1B1DI1PBDVSO
113897671
5310
8.80309 (−2)2.07488 (−1)8.78700 (−2)1.17774 (−1)
1.19400 (−1)6.87734 (−2)2.56159 (−2)1.70700 (−2)
20781628 969 969
10−3
D1B1DI1PBDVSO
1461118579
7021
2.68984 (−2)3.05138 (−2)2.12907 (−2)1.81903 (−2)
5.25240 (−2)6.28779 (−3)9.00022 (−3)6.05764 (−3)
25602071 1135 1063
10−4
D1B1DI1PBDVSO
15117094149
2110
4.30130 (−3)5.37539 (−4)4.25614 (−3)1.97342 (−5)
5.45887 (−3)1.06165 (−4)1.53472 (−3)5.05592 (−6)
30243169 1284 1951
10−5
D1B1DI1PBDVSO
218200161160
1700
1.84329 (−5)4.40121 (−4)6.68588 (−4)1.83864 (−5)
1.33438 (−5)1.37204 (−4)2.70629 (−4)6.47215 (−6)
38603652 2094 2080
10−6
D1B1DI1PBDVSO
275121179176
3010
1.80685 (−6)1.38184 (−5)3.15166 (−4)5.88468 (−6)
1.81711 (−6)2.65555 (−6)1.29162 (−4)2.26518 (−6)
52274077 2294 2272
10−7
D1B1DI1PBDVSO
219330193197
9010
4.47927 (−6)1.06409 (−7)1.44784 (−4)2.58941 (−7)
7.17250 (−6)6.01321 (−8)5.87351 (−5)9.00261 (−8)
54136257 2502 2513
10−8
D1B1DI1PBDVSO
348287205209
2200
3.49871 (−8)4.59227 (−7)6.69118 (−5)6.96499 (−8)
2.75009 (−8)1.19307 (−7)2.76364 (−5)2.35911 (−8)
63105367 8665 2664
10−9
D1B1DI1PBDVSO
260424226225
7000
2.68246 (−8)4.48526 (−10)3.11980 (−5)1.87994 (−9)
1.79939 (−8)1.44218 (−10)1.25997 (−5)4.25589 (−10)
86338271 29442851
10−10
D1B1DI1PBDVSO
526475370248
91000
3.45300 (−9)7.88269 (−9)1.44822 (−7)4.13390 (−9)
2.61398 (−9)1.39352 (−9)5.81042 (−8)8.63144 (−10)
98668822 4639 3109Table 4
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem2.
TOL
MTD
STEPS
FS
MAXE
AVER
TIME (total)
10−2
D1B1DI1PBDVSO
13513496173
0000
4.83764 (−2)2.63945 (−2)4.99681 (−3)1.01436 (−4)
3.21705 (−2)2.14665 (−2)2.47698 (−3)8.58874 (−5)
66214358 1961 3218
10−3
D1B1DI1PBDVSO
175121157162
1000
8.81849 (−3)3.90401 (−2)6.91170 (−5)2.40484 (−3)
5.48179 (−3)1.65274 (−2)4.87625 (−5)1.71894 (−3)
82943896 2890 3005
10−4
D1B1DI1PBDVSO
233173141142
2000
2.85719 (−3)1.44401 (−3)1.79264 (−4)2.60404 (−4)
1.39202 (−3)9.91992 (−4)1.19952 (−4)1.89056 (−4)
116845520 2622 2635
10−5
D1B1DI1PBDVSO
226263238237
1200
3.05429 (−4)1.97450 (−4)5.14212 (−6)1.03862 (−5)
1.22401 (−4)5.02997 (−5)3.89969 (−6)8.18350 (−6)
110958492 4137 4268
10−6
D1B1DI1PBDVSO
337262214218
1000
1.12897 (−5)6.52399 (−5)1.18826 (−5)1.56415 (−6)
5.53620 (−6)2.34931 (−5)9.01852 (−6)1.09680 (−6)
160608429 3766 3941
10−7
D1B1DI1PBDVSO
421337368370
2000
4.00671 (−7)1.70966 (−6)1.79221 (−7)2.75078 (−7)
4.00671 (−7)1.05119 (−6)1.35202 (−7)1.91969 (−7)
1973510899 6213 6584
10−8
D1B1DI1PBDVSO
610523335336
4310
3.05239 (−7)4.90321 (−8)3.11508 (−8)3.64286 (−8)
2.22380 (−7)3.52515 (−8)2.58590 (−8)1.21667 (−8)
2906016840 5641 5996
10−9
D1B1DI1PBDVSO
534555578575
2510
4.96195 (−8)9.16995 (−8)1.20539 (−8)3.80099 (−9)
1.71145 (−8)5.35093 (−8)9.25523 (−9)2.95718 (−9)
2523417610 9347 10114
10−10
D1B1DI1PBDVSO
800671516517
2400
2.12904 (−10)8.63683 (−10)8.37294 (−10)3.30834 (−9)
1.21736 (−10)3.89440 (−10)6.51148 (−10)2.52921 (−9)
3823521617 8541 9011Table 5
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem3.
TOL
MTD
y(1)
10−2
1PBDVSODIB1D1
1.9788520539e + 0001.9163232164e + 0001.8780198255e + 000−1.#IND000000e + 000
10−4
1PBDVSODIB1D1
1.8706114418e + 0001.8694497531e + 0001.8694030590e + 000−1.#IND000000e + 000
10−6
1PBDVSODIB1D1
1.8694735377e + 0001.8694497531e + 0001.8694368915e + 000−1.#IND000000e + 000
10−8
1PBDVSODIB1D1
1.8694389431e + 0001.8694388843e + 0001.8694402726e + 000−1.#IND000000e + 000
10−10
1PBDVSODIB1D1
1.8694388662e + 0001.8694388834e + 000Fail to compute−1.#IND000000e + 000Table 6
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem4.
TOL
MTD
y(1)
10−4
1PBDVSODIB1D1
−1.5171246226e  −  009−2.5152953310e  −  0089.9810188187e  −  009Fail to compute
10−6
1PBDVSODIB1D1
1.0035458019e  −  0081.0609214996e  −  0089.9998786959e  −  009Fail to compute
10−8
1PBDVSODIB1D1
1.0004172282e  −  0081.0004031492e  −  0081.0000000079e  −  008Fail to compute
10−10
1PBDVSODIB1D1
9.9999571459e  −  0099.9999855296e  −  0099.9999999555e  −  009Fail to computeThe result of the backwards difference method (1PBDVSO) is measured against D1, B1, and DI method for solving the same problems. The variable step size order codes are implemented for solving the problems. We begin to compute the solution using a higher tolerance (TOL
*
10
-
2).Problem 1.
Consider(53)
y
1
′′
=
-
y
1
r
3
,
y
1
(
0
)
=
1
,
y
1
′
(
0
)
=
0
y
2
′′
=
-
y
2
r
3
,
y
2
(
0
)
=
0
,
y
2
′
(
0
)
=
1
r
=
(
y
1
2
+
y
2
2
)
1
/
2
0
≤
x
≤
16
π
.
Solution is as follows:
(54)
y
1
(
x
)
=
cos
⁡
⁡
x
,
y
2
(
x
)
=
sin
⁡
x
.
First order system is as follows:
(55)
y
1
′
=
y
3
,
y
2
′
=
y
4
,
y
3
′
=
-
y
1
r
3
,
y
4
′
=
-
y
2
r
3
y
1
(
0
)
=
1
,
y
2
(
0
)
=
0
,
y
3
(
0
)
=
0
,
y
4
(
0
)
=
1
.
Solution is as follows:
(56)
y
1
(
x
)
=
cos
⁡
⁡
x
,
y
2
(
x
)
=
sin
⁡
x
,
y
3
(
x
)
=
-
sin
⁡
x
,
y
4
(
x
)
=
cos
⁡
⁡
x
.
For source, see Shampine and Gordon [7].Problem 2.
Consider(57)
y
(
8
)
=
y
,
y
′
(
0
)
=
y
′′
(
0
)
=
⋯
=
y
(
7
)
(
0
)
=
1
,
0
≤
x
≤
100
.
Solution is as follows:
(58)
y
=
e
x
.
First order system is as follows:
(59)
y
1
′
=
y
2
,
y
2
′
=
y
3
,
y
3
′
=
y
4
,
y
4
′
=
y
5
,
y
5
′
=
y
6
,
y
6
′
=
y
7
,
y
7
′
=
y
8
,
y
8
′
=
y
1
y
1
(
0
)
=
y
2
(
0
)
=
⋯
=
y
8
(
0
)
=
1
.
Solution is as follows:
(60)
y
1
(
x
)
=
y
2
(
x
)
=
y
3
(
x
)
=
⋯
=
y
8
(
x
)
=
e
x
.
For source, see Suleiman [2].Problem 3.
Van der Pol’s equation is as follows:(61)
y
′′
=
5
(
1
-
y
2
)
y
′
-
y
y
(
0
)
=
2
,
y
′
(
0
)
=
0
,
0
≤
x
≤
1
.
First order system is as follows:
(62)
y
1
′
=
y
2
,
y
2
′
=
5
(
1
-
y
1
2
)
y
2
-
y
1
,
y
1
(
0
)
=
2
,
y
2
(
0
)
=
0
.For source, see Suleiman [2].Problem 4.
Control theory is as follows:(63)
y
(
i
v
)
=
(
y
2
-
sin
⁡
y
-
100
4
)
y
+
(
y
′
y
′′
1
(
y
2
+
1
)
-
4
×
100
3
)
y
′
+
(
1
-
6
×
100
2
)
y
′′
+
(
10
e
y
′′′
2
-
4
×
100
)
y
′′′
+
1
y
(
0
)
=
y
′
(
0
)
=
y
′′
(
0
)
=
y
′′′
(
0
)
=
0
,
0
≤
x
≤
1
.
First order system is as follows:
(64)
y
1
′
=
y
2
,
y
2
′
=
y
3
,
y
3
′
=
y
4
,
y
4
′
=
(
y
1
2
-
sin
⁡
y
-
100
4
)
y
1
+
(
y
2
y
3
1
(
y
1
2
+
1
)
-
4
×
100
3
)
y
2
+
(
1
-
6
×
100
2
)
y
3
+
(
10
e
y
4
2
-
4
×
100
)
y
4
+
1
y
1
(
0
)
=
y
2
(
0
)
=
y
3
(
0
)
=
y
4
(
0
)
=
0
.For source, see Enright et al. [10].


## 9.1. Numerical Results

Tables3, 4, 5, and 6 show the numerical results for Problems 1–4, when solved with direct integration, backwards difference, and also reducing to first order systems. Problems 1 and 2 are well-behaved problems whereas Problems 3 and 4 are without exact solutions. For the first two problems, we evaluate the maximum and average values of the error in the computed solution y. The definition of the error is as defined above. The following notation will indicateTIME: the execution time taken in microseconds,
FS: the number of failed steps,
STEPS: total steps,
MTD: the name of the method,
SUCSTEP: the number of successful step,
MAXE: the maximum error,
AVER: the average error,
TOL: the tolerance used,
DI: direct integration,
D1: direct integration with the reduction to first order system,
1PBDVSO: backward difference,
B1: backward difference with the reduction to first order system.The errors calculated are defined as
(51)
(
e
i
)
t
=
|
(
y
i
)
t
-
(
y
(
x
i
)
)
t
A
+
B
(
y
(
x
i
)
)
t
|
with (
y
)
t as the tth component of y
~. A
=
1
,
B
=
0 correspond to the absolute error test, A
=
1
,
B
=
1 correspond to the mixed error test, and A
=
0
,
B
=
1 correspond to the relative error test. The mixed error test is used for all test problems. The maximum and average errors are given by
(52)
MAXE
=
max
⁡
1
<
i
<
SUCSTEP
(
max
⁡
1
<
i
<
N
⁡
(
e
i
)
t
)
,
AVER
=
∑
i
=
1
SUCSTEP
∑
t
=
1
N
(
e
i
)
t
(
N
)
(
SUCSTEP
)
,
where N is the number of equations in the systems.Table 3
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem1.
TOL
MTD
STEPS
FS
MAXE
AVER
TIME (total)
10−2
D1B1DI1PBDVSO
113897671
5310
8.80309 (−2)2.07488 (−1)8.78700 (−2)1.17774 (−1)
1.19400 (−1)6.87734 (−2)2.56159 (−2)1.70700 (−2)
20781628 969 969
10−3
D1B1DI1PBDVSO
1461118579
7021
2.68984 (−2)3.05138 (−2)2.12907 (−2)1.81903 (−2)
5.25240 (−2)6.28779 (−3)9.00022 (−3)6.05764 (−3)
25602071 1135 1063
10−4
D1B1DI1PBDVSO
15117094149
2110
4.30130 (−3)5.37539 (−4)4.25614 (−3)1.97342 (−5)
5.45887 (−3)1.06165 (−4)1.53472 (−3)5.05592 (−6)
30243169 1284 1951
10−5
D1B1DI1PBDVSO
218200161160
1700
1.84329 (−5)4.40121 (−4)6.68588 (−4)1.83864 (−5)
1.33438 (−5)1.37204 (−4)2.70629 (−4)6.47215 (−6)
38603652 2094 2080
10−6
D1B1DI1PBDVSO
275121179176
3010
1.80685 (−6)1.38184 (−5)3.15166 (−4)5.88468 (−6)
1.81711 (−6)2.65555 (−6)1.29162 (−4)2.26518 (−6)
52274077 2294 2272
10−7
D1B1DI1PBDVSO
219330193197
9010
4.47927 (−6)1.06409 (−7)1.44784 (−4)2.58941 (−7)
7.17250 (−6)6.01321 (−8)5.87351 (−5)9.00261 (−8)
54136257 2502 2513
10−8
D1B1DI1PBDVSO
348287205209
2200
3.49871 (−8)4.59227 (−7)6.69118 (−5)6.96499 (−8)
2.75009 (−8)1.19307 (−7)2.76364 (−5)2.35911 (−8)
63105367 8665 2664
10−9
D1B1DI1PBDVSO
260424226225
7000
2.68246 (−8)4.48526 (−10)3.11980 (−5)1.87994 (−9)
1.79939 (−8)1.44218 (−10)1.25997 (−5)4.25589 (−10)
86338271 29442851
10−10
D1B1DI1PBDVSO
526475370248
91000
3.45300 (−9)7.88269 (−9)1.44822 (−7)4.13390 (−9)
2.61398 (−9)1.39352 (−9)5.81042 (−8)8.63144 (−10)
98668822 4639 3109Table 4
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem2.
TOL
MTD
STEPS
FS
MAXE
AVER
TIME (total)
10−2
D1B1DI1PBDVSO
13513496173
0000
4.83764 (−2)2.63945 (−2)4.99681 (−3)1.01436 (−4)
3.21705 (−2)2.14665 (−2)2.47698 (−3)8.58874 (−5)
66214358 1961 3218
10−3
D1B1DI1PBDVSO
175121157162
1000
8.81849 (−3)3.90401 (−2)6.91170 (−5)2.40484 (−3)
5.48179 (−3)1.65274 (−2)4.87625 (−5)1.71894 (−3)
82943896 2890 3005
10−4
D1B1DI1PBDVSO
233173141142
2000
2.85719 (−3)1.44401 (−3)1.79264 (−4)2.60404 (−4)
1.39202 (−3)9.91992 (−4)1.19952 (−4)1.89056 (−4)
116845520 2622 2635
10−5
D1B1DI1PBDVSO
226263238237
1200
3.05429 (−4)1.97450 (−4)5.14212 (−6)1.03862 (−5)
1.22401 (−4)5.02997 (−5)3.89969 (−6)8.18350 (−6)
110958492 4137 4268
10−6
D1B1DI1PBDVSO
337262214218
1000
1.12897 (−5)6.52399 (−5)1.18826 (−5)1.56415 (−6)
5.53620 (−6)2.34931 (−5)9.01852 (−6)1.09680 (−6)
160608429 3766 3941
10−7
D1B1DI1PBDVSO
421337368370
2000
4.00671 (−7)1.70966 (−6)1.79221 (−7)2.75078 (−7)
4.00671 (−7)1.05119 (−6)1.35202 (−7)1.91969 (−7)
1973510899 6213 6584
10−8
D1B1DI1PBDVSO
610523335336
4310
3.05239 (−7)4.90321 (−8)3.11508 (−8)3.64286 (−8)
2.22380 (−7)3.52515 (−8)2.58590 (−8)1.21667 (−8)
2906016840 5641 5996
10−9
D1B1DI1PBDVSO
534555578575
2510
4.96195 (−8)9.16995 (−8)1.20539 (−8)3.80099 (−9)
1.71145 (−8)5.35093 (−8)9.25523 (−9)2.95718 (−9)
2523417610 9347 10114
10−10
D1B1DI1PBDVSO
800671516517
2400
2.12904 (−10)8.63683 (−10)8.37294 (−10)3.30834 (−9)
1.21736 (−10)3.89440 (−10)6.51148 (−10)2.52921 (−9)
3823521617 8541 9011Table 5
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem3.
TOL
MTD
y(1)
10−2
1PBDVSODIB1D1
1.9788520539e + 0001.9163232164e + 0001.8780198255e + 000−1.#IND000000e + 000
10−4
1PBDVSODIB1D1
1.8706114418e + 0001.8694497531e + 0001.8694030590e + 000−1.#IND000000e + 000
10−6
1PBDVSODIB1D1
1.8694735377e + 0001.8694497531e + 0001.8694368915e + 000−1.#IND000000e + 000
10−8
1PBDVSODIB1D1
1.8694389431e + 0001.8694388843e + 0001.8694402726e + 000−1.#IND000000e + 000
10−10
1PBDVSODIB1D1
1.8694388662e + 0001.8694388834e + 000Fail to compute−1.#IND000000e + 000Table 6
Comparison between the 1PBDVSO, DI, B1, and D1 methods for solving Problem4.
TOL
MTD
y(1)
10−4
1PBDVSODIB1D1
−1.5171246226e  −  009−2.5152953310e  −  0089.9810188187e  −  009Fail to compute
10−6
1PBDVSODIB1D1
1.0035458019e  −  0081.0609214996e  −  0089.9998786959e  −  009Fail to compute
10−8
1PBDVSODIB1D1
1.0004172282e  −  0081.0004031492e  −  0081.0000000079e  −  008Fail to compute
10−10
1PBDVSODIB1D1
9.9999571459e  −  0099.9999855296e  −  0099.9999999555e  −  009Fail to computeThe result of the backwards difference method (1PBDVSO) is measured against D1, B1, and DI method for solving the same problems. The variable step size order codes are implemented for solving the problems. We begin to compute the solution using a higher tolerance (TOL
*
10
-
2).Problem 1.
Consider(53)
y
1
′′
=
-
y
1
r
3
,
y
1
(
0
)
=
1
,
y
1
′
(
0
)
=
0
y
2
′′
=
-
y
2
r
3
,
y
2
(
0
)
=
0
,
y
2
′
(
0
)
=
1
r
=
(
y
1
2
+
y
2
2
)
1
/
2
0
≤
x
≤
16
π
.
Solution is as follows:
(54)
y
1
(
x
)
=
cos
⁡
⁡
x
,
y
2
(
x
)
=
sin
⁡
x
.
First order system is as follows:
(55)
y
1
′
=
y
3
,
y
2
′
=
y
4
,
y
3
′
=
-
y
1
r
3
,
y
4
′
=
-
y
2
r
3
y
1
(
0
)
=
1
,
y
2
(
0
)
=
0
,
y
3
(
0
)
=
0
,
y
4
(
0
)
=
1
.
Solution is as follows:
(56)
y
1
(
x
)
=
cos
⁡
⁡
x
,
y
2
(
x
)
=
sin
⁡
x
,
y
3
(
x
)
=
-
sin
⁡
x
,
y
4
(
x
)
=
cos
⁡
⁡
x
.
For source, see Shampine and Gordon [7].Problem 2.
Consider(57)
y
(
8
)
=
y
,
y
′
(
0
)
=
y
′′
(
0
)
=
⋯
=
y
(
7
)
(
0
)
=
1
,
0
≤
x
≤
100
.
Solution is as follows:
(58)
y
=
e
x
.
First order system is as follows:
(59)
y
1
′
=
y
2
,
y
2
′
=
y
3
,
y
3
′
=
y
4
,
y
4
′
=
y
5
,
y
5
′
=
y
6
,
y
6
′
=
y
7
,
y
7
′
=
y
8
,
y
8
′
=
y
1
y
1
(
0
)
=
y
2
(
0
)
=
⋯
=
y
8
(
0
)
=
1
.
Solution is as follows:
(60)
y
1
(
x
)
=
y
2
(
x
)
=
y
3
(
x
)
=
⋯
=
y
8
(
x
)
=
e
x
.
For source, see Suleiman [2].Problem 3.
Van der Pol’s equation is as follows:(61)
y
′′
=
5
(
1
-
y
2
)
y
′
-
y
y
(
0
)
=
2
,
y
′
(
0
)
=
0
,
0
≤
x
≤
1
.
First order system is as follows:
(62)
y
1
′
=
y
2
,
y
2
′
=
5
(
1
-
y
1
2
)
y
2
-
y
1
,
y
1
(
0
)
=
2
,
y
2
(
0
)
=
0
.For source, see Suleiman [2].Problem 4.
Control theory is as follows:(63)
y
(
i
v
)
=
(
y
2
-
sin
⁡
y
-
100
4
)
y
+
(
y
′
y
′′
1
(
y
2
+
1
)
-
4
×
100
3
)
y
′
+
(
1
-
6
×
100
2
)
y
′′
+
(
10
e
y
′′′
2
-
4
×
100
)
y
′′′
+
1
y
(
0
)
=
y
′
(
0
)
=
y
′′
(
0
)
=
y
′′′
(
0
)
=
0
,
0
≤
x
≤
1
.
First order system is as follows:
(64)
y
1
′
=
y
2
,
y
2
′
=
y
3
,
y
3
′
=
y
4
,
y
4
′
=
(
y
1
2
-
sin
⁡
y
-
100
4
)
y
1
+
(
y
2
y
3
1
(
y
1
2
+
1
)
-
4
×
100
3
)
y
2
+
(
1
-
6
×
100
2
)
y
3
+
(
10
e
y
4
2
-
4
×
100
)
y
4
+
1
y
1
(
0
)
=
y
2
(
0
)
=
y
3
(
0
)
=
y
4
(
0
)
=
0
.For source, see Enright et al. [10].


## 10. Comments on Numerical Results and Conclusion

Results for the first problem in Table3 whose solution has the trigonometric form show that the 1PBDVSO is quicker in terms of execution times compared to the other three methods. This is not surprising because the divided difference has an element of division whereas backward difference is without one. And since the solution involves trigonometric functions, the divided difference can be small and can magnify round-off errors. This can clearly be seen as we decrease the tolerance. The solution to Problem 2 has the solution in the positive exponential form and, therefore, the effect in terms of accumulated errors between divided difference and backward difference is small. This is reflected in the graph of accuracy versus execution times in Figure 2, where 1PBDVSO and DI methods almost overlapped each other.In Problem3, the reduction to first order system using divided difference cannot solve the Van der Pol equation while the backward difference could. However, the DI direct method and 1PBDVSO are successful in solving this equation. This conclusion is judged by the closeness of the solution at x
=
1. This also indicates the instability of the problem when reduced to first order and solved using divided difference. For Problem 4, which has no given solution, similar behavior pattern occurs in solving the problem. This indicates the preferability of the high order problems being solved directly rather than reduction to first order. Overall, from the tables, we see that the accuracy of 1PBDVSO is better compared to the other methods. In order to see this more clearly, we present the graphs of execution times against accuracy. To determine the more accurate relationship, we take any particular values of the abscissae (the execution times) and the ordinates (accuracy) which are lowest in values representing more accurate points. For graph in Figure 1, the 1PBDVSO is the most accurate method since, for all the abscissae chosen from the graphs, they are below the values of the corresponding ordinates.Figure 1
Comparison of efficiency between the 1PBDVSO, DI, B1, and D1 methods for Problem1.Figure 2
Comparison of efficiency between the 1PBDVSO, DI, B1, and D1 methods for Problem2.The numerical results show that the 1PBDVSO is better than DI in terms of accuracy, competitive in terms of execution times, and more stable when tested with more demanding problems. We are inclined to choose 1PBDVSO because, for some problems, they are more accurate. It can be seen from the figures above, where the undermost curve is more efficient. Using this simple performance index, the point which is lowest at a particular time is more accurate. These results show that the variable order step size method gives acceptable solutions and is suitable for solving higher order ODEs directly.



---

*Source: 565137-2014-08-19.xml*