# An Efficient Zero-Stable Numerical Method for

Authors: S. J. Kayode
Journal: International Journal of Mathematics and Mathematical Sciences
 (2008)

---

## Abstract

The purpose of this paper is to produce an efficient zero-stable numerical
method with the same order of accuracy as that of the main starting values (predictors) for direct solution of fourth-order differential equations without reducing it to a system of first-order equations. The method of collocation of the differential system arising from the approximate solution to the problem is adopted using the power series as a basis function. The method is consistent, symmetric, and of optimal orderp
=
6.
The main predictor for the method is also consistent, symmetric, zero-stable, and of optimal order p
=
6.


---

## Body

## 1. Introduction

Many problems leading to fourth-order differential equations
of the form(1.1)
y
i
v
=
f
(
t
,
y
,
y
1
,
y
″
,
y
′
′
′
)
,


y
(
m
−
1
)
(
t
0
)
=
μ
m
−
1
,

m
=
1
,
2
,
3
,
4, are of great interest to scientists
and engineers  formulating mathematical
models for elastic bodies. Investigators, both theoretical and numerical
analysts, have concerned themselves with the
study and solutions of such equations. However, only a limited number of
analytical methods are available for solving () directly without reduction to
a first-order system of initial value problems (ivps) of type (1.1)
y
′
=
f
(
t
,
y
)
,

y
(
t
0
)
=
μ
,

f
∈
C
1
[
a
,
b
]
,
t
∈
R
,
y
∈
R
m
. are of great interest to scientists
and engineers formulating mathematical
models for elastic bodies. Investigators, both theoretical and numerical
analysts, have concerned themselves with the
study and solutions of such equations. However, only a limited number of
analytical methods are available for solving (1.1) directly without reduction to
a first-order system of initial value problems (ivps) of type (1.2)
y
(
x
)
=
∑
j
=
0
k
+
2
a
j
x
j
. The many numerical methods adopted for such higher-order
differential equations (1.1) are only for handling first-order equations (1.2),
[1–4].The reduction of such
problems of type (1.1) to
systems of first-order equations, to our knowledge,
leads to serious computational burden as well as wastage in computer time. These
setbacks had been addressed by some authors [5–7]. In recent times,
attempts have been made to formulate numerical algorithms capable of solving
problem (1.1) directly using differing methods of derivation with varying degree
of accuracy [8–12]. These attempts are not without their associated limitations, some of which
include low stepnumber k, too many function evaluations, and the combination of lower-order
predictors with the correctors in the predictor-corrector methods [13, 14]. All
these have effects on the accuracy of the methods.The proposed zero-stable order six methods
are designed to have higher step number and reduced function evaluations,
in order to address these observed
problems. The derived method is capable of handling problems wheref is either linear or nonlinear. The numerical results
of the method are compared with an existing method [3].


## 2. Description of the Method

Let the
approximate solution to problem (1.1) be a partial sum of a
power series of the form: (2.1)
∑
j
=
4
k
+
2
j
(
j
−
1
)
(
j
−
2
)
(
j
−
3
)
a
j
x
j
−
4
=
f
(
x
,
y
,
y
′
,
y
′
′
,
y
′
′
′
)
.   Taking the fourth derivative of (2.1) and
using this in (1.1) yields(2.2)
A
X
=
B
, To be able to obtain an open method
with three function evaluations, collocation points are taken in (2.2) at all odd
grid points x
=
x
n
+
i
,
i
=
1
(
2
)
k.
Equation (1.2) is interpolated at all grid points x
=
x
n
+
i
,
i
=
0
(
1
)
k
−
1
,
x
∈
R,
where R is the set of real numbers. The choice of number of interpolation
points was guided by the intended order of the method as well as the number of
collocation points.The system of equations obtained from
the collocation and interpolation above is represented by the matrix equation:(2.3)
A
m
×
m
=
[
0
0
0
0
24
120
x
n
+
i
⋯
ξ
x
n
+
1
k
−
2
0
0
0
0
24
120
x
n
+
3
⋯
ξ
x
n
+
3
k
−
2
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
0
0
0
0
24
120
x
n
+
k
⋯
ξ
x
n
+
k
k
−
2
1
x
n
x
n
2
x
n
3
x
n
4
x
n
5
⋯
x
n
k
+
2
1
x
n
+
1
x
n
+
1
2
x
n
+
1
3
x
n
+
1
4
x
n
+
1
5
⋯
x
n
+
1
k
+
2
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
1
x
n
+
k
−
1
x
n
+
k
−
1
2
x
n
+
k
−
1
3
x
n
+
k
−
1
4
⋯
⋯
x
n
+
k
−
1
k
+
2
]
,where A,
an m by m matrix, X and B are, respectively, given by(2.4)
y
(
t
)
=
∑
j
=
0
k
−
1
α
j
(
t
)
y
n
+
j
+
h
4
∑
j
=
0
k
β
j
(
t
)
f
n
+
j
,

β
k
⁢
≠
0.
X
=
[
a
0
a
1
a
2
…
a
k
+
1
a
k
+
2
]
T,
B
=
[
f
n
+
1
f
n
+
3
…
f
n
+
k
y
n
y
n
+
1
…
y
n
+
k
−
1
]
T,
ξ
=
(
k
−
1
)
k
(
k
+
1
)
(
k
+
2
)
,f
n
+
i
=
f
(
x
n
+
i
,
y
n
+
i
,
y
n
+
i
′
,
y
n
+
i
′
′
,
y
n
+
i
′
′
′
)
,
i
=
0
,
1
,
2
,
…
,
y
n
+
i
=
y
(
x
n
+
i
),  and T is the usual matrix transpose.Using the Gaussian method to solve
the matrix (2.3), the values of a
j
′
s
,
j
=
0
,
1
,
…
,
k
+
2,
for k
=
5, are obtained as shown in
Appendix A. These values of a
j
′
s are substituted into (2.1). By using the transformation x
=
t
h
+
x
n
+
4
,
t
∈
(
0
,
1
]   with h being steplength, in the resulting equation after the substitution, an
open method with variable coefficients is obtained as(2.5)
y
n
+
5
=
4
y
n
+
4
−
6
y
n
+
3
+
4
y
n
+
2
−
y
n
+
1
+
h
4
24
(
f
n
+
5
+
22
f
n
+
3
+
f
n
+
1
)
, Taking k
=
5, the coefficients α
j
(
x
) and β
j
(
x
) as well as their first, second, and third
derivatives in (2.5) are computed and, respectively, obtained as shown in
Appendix
B.The values oft could be taken in the interval I
=
(
0
,
1
] in
order to obtain a particular discrete scheme. In this work, we take t
=
1 in (2.5) to have a zero-stable method and its
derivatives as(2.6)
y
n
+
5
′
=
1
210
h
{
853
y
n
+
4
−
1767
y
n
+
3
+
1128
y
n
+
2
−
157
y
n
+
1
−
57
y
n





+
h
4
48
(
1319
f
n
+
5
+
20738
f
n
+
3
+
1679
f
n
+
1
)
}
,
(2.7)
y
n
+
5
′
′
=
1
60
h
2
{
119
y
n
+
4
−
236
y
n
+
3
+
54
y
n
+
2
+
124
y
n
+
1
−
60
y
n





+
h
4
8
(
159
f
n
+
5
+
1526
f
n
+
3
+
203
f
n
+
1
)
}
,
(2.8)
y
n
+
5
′
′
′
=
1
40
h
3
{
−
23
y
n
+
4
+
132
y
n
+
3
−
258
y
n
+
2
+
212
y
n
+
1
−
63
y
n





+
h
4
12
(
317
f
n
+
5
+
1364
f
n
+
3
+
275
f
n
+
1
)
}
,(2.9)
y
(
t
)
=
∑
j
=
0
k
−
1
α
j
(
t
)
y
n
+
j
+
∑
j
=
0
k
−
1
β
j
(
t
)
f
n
+
j
,

k
=
5
,
4. respectively.The order of the open method (2.6) and
its derivatives (2.7), (2.8), (2.9)
is the same and is p
=
6. The error constant of (2.6) is found to be C
p
+
2
≈
−
4.31
×
10
−
2.
The method is consistent and zero stable, satisfying the necessary and
sufficient conditions for convergence of linear multistep methods [3, 4].


## 3. Derivation of Starting Values for the Method

Method (2.6) above and its derivatives
are open; β
k
≠
0.
This implies that the sample discrete method (2.6) and its derivatives (2.7), (2.8),
(2.9) require additional closed starting values, β
k
=
0,
for the evaluation of f
n
+
i
,
i
=
1
,
2
,
…
,
k.
In this work, efforts are made at obtaining the main closed predictors of the
same order of accuracy as starting values for f
n
+
i
,
i
=
4
,
…
,
k.
These predictors are obtained by using the same approach described in Section 2.
The following consistent and zero-stable closed methods are obtained, for k
=
5
,
4,(3.1)
y
n
+
5
=
4
y
n
+
4
−
6
y
n
+
3
+
4
y
n
+
2
−
y
n
+
1
+
h
4
6
(
f
n
+
4
+
4
f
n
+
3
+
f
n
+
2
)
,   Taking t
=
1 and for k
=
5, we have
(3.2)
y
N
+
5
′
=
1
210
h
{
898
y
n
+
4
−
1947
y
n
+
3
+
1398
y
n
+
2
−
337
y
n
+
1
−
12
y
n





+
h
4
12
(
159
f
n
+
4
+
132
f
n
+
3
+
115
f
n
+
2
)
}
,
(3.3)
y
n
+
5
′
′
=
1
15
h
2
{
38
y
n
+
4
−
92
y
n
+
3
+
63
y
n
+
2
−
2
y
n
+
1
−
7
y
n





+
h
4
8
(
159
f
n
+
4
+
132
f
n
+
3
+
115
f
n
+
2
)
}
,
(3.4)
y
n
+
5
′
′
′
=
−
1
10
h
3
{
11
y
n
+
4
−
54
y
n
+
3
+
96
y
n
+
2
−
74
y
n
+
1
+
21
y
n





+
h
4
12
(
317
f
n
+
4
−
124
f
n
+
3
+
359
f
n
+
2
)
}
,(3.5)
y
n
+
4
=
4
y
n
+
3
−
6
y
n
+
2
+
4
y
n
+
1
−
y
n
+
h
4
6
(
f
n
+
3
+
4
f
n
+
2
+
f
n
+
1
)
, The main predictors (3.2)
and its associated derivatives (3.3), (3.4), and (3.5) are each of order 6 and
their error constants C
p
+
2 are, respectively, 1.389
×
10
−
3,  6.948
×
10
−
3, 7.0367
×
10
−
2, 3.269
×
10
−
1.Fork
=
4, we have the following:
(3.6)
y
n
+
4
(
1
)
=
1
6
h
{
26
y
n
+
3
−
57
y
n
+
2
+
42
y
n
+
1
−
11
y
n
+
h
4
60
(
185
f
n
+
3
+
452
f
n
+
2
+
113
f
n
+
1
)
}
, order p
=
6, C
p
+
2
≈
1.389
×
10
−
3;(3.7)
y
n
+
4
(
2
)
=
1
h
2
{
3
y
n
+
3
−
8
y
n
+
2
+
7
y
n
+
1
−
2
y
n
+
h
4
60
(
449
f
n
+
3
+
452
f
n
+
2
+
149
f
n
+
1
)
}
, order p
=
5, C
p
+
2
≈
9.524
×
10
−
3;(3.8)
y
n
+
4
(
3
)
=
1
h
3
{
y
n
+
3
−
3
y
n
+
2
+
3
y
n
+
1
−
y
n
+
h
4
24
(
55
f
n
+
4
8
f
n
+
2
+
13
f
n
+
1
)
}
, order p
=
5, C
p
+
2
≈
7.778
×
10
−
2;(3.9)
y
(
4
)
+
y
(
2
)
=
0
,

0
≤
x
≤
π
2
,

y
(
0
)
=
0
,
y
(
1
)
(
0
)
=
−
1.1
72
−
50
Π
,

y
(
2
)
(
0
)
=
1
144
−
100
Π
,
y
(
3
)
(
0
)
=
1.2
144
−
100
Π
.
Theoretical solution:
y
(
x
)
=
(
1
−
x
−
cos
⁢
⁡
x
−
1.2
sin
⁢
⁡
x
)
(
144
−
100
π
)
. order p
=
5, C
p
+
2
≈
0.35.Fork
<
4.To be able to have sufficient collocation and
interpolation points at the grid points that will not produce system of
underdetermined equations, the minimum stepnumberk required for the derivation
of any method for fourth-order problem is 4. Therefore, Taylor
series expansion was adopted to evaluate y
n
+
i
,
y
n
+
i
(
1
)
,
y
n
+
i
(
2
)
,
y
n
+
i
(
3
)
,
i
=
1
,
2
,
3;
see Appendix C.The initial values for the problems
are used to evaluatey
0
,
y
0
(
1
)
,
y
0
(
2
)
,
y
0
(
3
),
where y
0
(
r
)
=
y
(
r
)
(
x
0
)
,
r
=
0
,
1
,
2
,
3.


## 4. Numerical Experience

The accuracy of the proposed method (2.6)
has been tested, and the results of the following two initial value problems
are shown in the tables in Appendix C below.Test Problem 1
(4.1)
y
(
4
)
=
(
y
(
1
)
)
2
−
y
y
(
#
)
−
4
x
2
+
e
x
(
1
−
4
x
+
x
2
)
,

0
≤
x
≤
1
,
y
(
0
)
=
1
,

y
(
1
)
(
0
)
=
1
,

y
(
2
)
(
0
)
=
3
,

y
(
3
)
(
0
)
=
1.
Theoretical  solution:
y
(
x
)
=
x
2
+
e
x
.Test Problem 2
(4.2)
y
(
4
)
=
(
y
(
1
)
)
2
−
y
y
(
#
)
−
4
x
2
+
e
x
(
1
−
4
x
+
x
2
)
,
0
≤
x
≤
1
,
y
(
0
)
=
1
,
y
(
1
)
(
0
)
=
1
,
y
(
2
)
(
0
)
=
3
,
y
(
3
)
(
0
)
=
1.
Theoretical  solution:
y
(
x
)
=
x
2
+
e
x
.The test problems 1 and 2 were solved
with method (2.6). The steplength h
=
1
/
32 is used for the purpose of comparison
with the existing method Kayode [3]. The maximum errors obtained are compared with
the method in Kayode [3] as shown
in Tables 1
and 2
.Table 1
Comparison of errors in Kayode [3] and new method (2.6) for
Problem 1.
X
Y-exact
Y-computed with new method (2.6)
Errors in Kayode [3]
Errors in new method (2.6)
0.103125
0.13007993D-02
0.13007993D-02
0.49873300D-15
0.48355417D-16
0.206250
0.25317732D-02
0.25317732D-02
0.67654216D-15
0.13933299D-15
0.306250
0.36524788D-02
0.36524783D-02
0.31350790D-14
0.66893539D-15
0.406250
0.46959523D-02
0.46959523D-02
0.94360283D-14
0.20129384D-14
0.506250
0.56576413D-02
0.56576413D-02
0.22116857D-13
0.46736053D-14
0.603125
0.65077534D-02
0.65077534D-02
0.43379363D-13
0.91874598D-14
0.703125
0.72983134D-02
0.72983134D-02
0.77870869D-13
0.16069038D-13
0.803125
0.79985187D-02
0.79985187D-02
0.12863495D-12
0.25407974D-13
0.903125
0.86072451D-02
0.86072451D-02
0.19927116D-12
0.38108926D-13
1.003125
0.91242822D-02
0.91242822D-02
0.29323245D-12
0.54051538D-13Table 2
Comparison of errors in Kayode [3] and new method (2.6) for
Problem 2.
X
Y-exact
Y-computed with new method (2.6)
Errors in Kayode [3]
Errors in new method (2.6)
0.103125
0.11192647D+01
0.11192647D+01
4.6842974D-11
8.1714635D-12
0.206250
0.12715995D+01
0.12715995D+01
2.0687163D-10
1.4680923D-10
0.306250
0.1452111D+01
0.14521109D+01
9.0421981D-10
7.8813045D-10
0.406250
0.16662167D+01
0.16662169D+01
2.9137912D-09
2.6795135D-09
0.506250
0.19153471D+01
0.19153471D+01
7.5114031D-09
7.0615016D-09
0.603125
0.21915816D+01
0.21915816D+01
1.6231339D-08
1.5438411D-08
0.703125
0.25144403D+01
0.25144403D+01
3.2237145D-08
3.0900764D-08
0.803125
0.28775164D+01
0.28775164D+01
5.8891835D-08
5.6759556D-08
0.903125
0.32829362D+01
0.32829362D+01
1.0079927D-07
9.7549608D-08
1.003125
0.37330495D+01
0.37330495D+01
1.6373663D-07
1.5897083D-07


## 5. Conclusion

A predictor-corrector method whose main
predictors (3.2) and (3.6) have the same order of accuracy with the corrector (2.6)
was formulated. The method was tested to be consistent and zero stable with order
six. It is remarkable to note that the main predictors are of the same order 6 as
that of the method with comparable error constants. This, to a large extent,
reduces the effects that global error could have on the accuracy of the method.The method was used to solve both
linear and nonlinear problems of fourth-order differential equations without
reduction to system of first order equations.The methods (3.2) and (3.6) developed
as main predictors are of the same order of accuracy with the corrector (2.6). The
method in Kayode [3] has predictors of lower order of accuracy, this accounts
for an improvement of method (2.6) over
that of [3] as compared in the maximum errors shown in
Tables 1
and 2
.



---

*Source: 364021-2008-07-07.xml*