Derivation and Application of a Linear Multistep Numerical Scheme
Abdulrahman NDANUSA
Mathematics/Computer Science Department, Federal University of Technology, Minna, Nigeria
Abstract
A great many physical occurrences give rise to problems that often result in differential equations. When we solve a differential equation, we are in effect solving the physical problem it represents. Traditionally, solutions to differential equations were derived using analytical or exact methods. These solutions are often useful as they provide excellent insight into the behavior of some systems. However, certain differential equations are very difficult to solve by any means other than an approximate solution by the application of numerical methods. These methods can be classified into two thus: Onestep and multistep methods. However, in this work, our focus is on a class of multistep methods known as the Linear Multistep Methods. We thus use Taylor series expansion to derive a linear multistep method of order eight. The method is tested for consistency and zerostability in order to establish its convergence. Also, provided are some examples of problems solved using the new scheme, and the results compared with exact solution.
Introduction
We consider the initial value problem (IVP):
y^{'} = f(t,y), y(t_{0}) = y_{0}, t > t_{0} 
(1) 
Let y_{n} be an approximation to the theoretical solution at t_{n}, that is, to y(t_{n}), and let f_{n}=f(t_{n},y_{n}). Then, we say a linear multistep method (LMM) of step number k, or a linear k step method is a computational method for determining, the sequence {y_{n}} that takes the form of a linear relationship between y_{n+j}, f_{n+j}, j = 0,1, …, k.
Thus the general lmm may be written:
_{} 
(2) 
where α_{j} and β_{j} are constants; we assume α_{j} = 1 and that not both α_{0} and β_{0} are zero. We say that the method is explicit if β_{k} = 0, and implicit if β_{k} ≠ 0 [1].
Convergence is a minimal property which any acceptable linear multistep method must possess. Thus, the necessary and sufficient conditions for a linear multistep method to be convergent are that it be consistent and zerostable [1].
Consistency demands that,
ρ(1) = 0 and ρ^{'}(1) = σ(1) 
(3) 
where,
_{} 
(4) 
_{} 
(5) 
are called the first and second characteristic polynomials of (2) respectively.
The lmm (2) is said to be zerostable if no root of (4) has modulus greater than one, and if every root with modulus one is simple [1].
Derivation of the Lmm of Order 8
We use the method of Taylor series expansion described by [1] as follows:
Let _{} be the linear difference operator defined by:
_{} 
(6) 
where y(t) is an arbitrary function, continuously differentiable on [a , b].
Suppose we choose to expand y(t+jh) and y^{’}(t+jh) about t+rh; where r need not necessarily be an integer. We obtain:
_{} 
(7) 
The formulae for the constants D_{q} expressed in terms of α_{j}, β_{j} are:
_{} 
(8) 
q = 2, 3, …
In this research work, we wish to derive an optimal 6step method. Therefore, all the roots of the first characteristic polynomial _{} must be on the unit circle. We know that _{} is a polynomial of degree 6. Hence, by consistency, it has one real root at +1 and another real root at 1. The four remaining roots must be complex.
Hence we have:
ξ_{1} = +1 , ξ_{2} = 1, ξ_{3} = e^{iθ1}, ξ_{3} = e^{iθ1}, ξ_{3} = e^{iθ2}, ξ_{3} = e^{iθ2} 

Hence:
ρ(ξ) = (ξ  1)(ξ + 1)(ξ  e^{iθ1})(ξ  e^{iθ1})(ξ  e^{iθ2}) (ξe^{iθ2}) 

α_{6} = +1, α_{5} = 2(a+b), α_{4} = (4ab + 1), α_{3} = 0, α_{2} = (4ab+1), α_{1} = 2(a+b), α_{0} = 1 

We require the method to have order 8. We now state the order requirement in terms of the coefficients D_{q}_{ .}
From (8) we have the following:
D_{0} = α_{0} + α_{1} + α_{2} + α_{3} + α_{4} + α_{5} + α_{6} D_{1} = rα_{0} + (1  r)α_{1} + (2  r)α_{2} + (3  r)α_{3} + (4  r)α_{4} + (5  r)α_{5} + (6  r)α_{6}   (β_{0} + β_{1} + β_{2} + β_{3} + β_{4} + β_{5} + β_{6}) D_{2} = 1/2![( r)^{2}α_{0} + (1  r)^{2}α_{1} + (2  r)^{2}α_{2} + (3  r)^{2}α_{3} + (4  r)^{2}α_{4} + (5  r)^{2}α_{5} + + (6  r)^{2}α_{6}]  [rβ_{0} + (1  r)β_{1} + (2 r)β_{2} + (3  r)β_{3} + (4  r)β_{4} + (5  r)β_{5} + + (6  r)β_{6}] . . . D_{8} = 1/8![( r)^{8}α_{0} + (1  r)^{8}α_{1} + (2  r)^{8}α_{2} + (3  r)^{8}α_{3} + (4  r)^{8}α_{4} + (5  r)^{8}α_{5} + + (6  r)^{8}α_{6}]  1/7![(r)^{ 7}β_{0} + (1  r)^{7}β_{1} + (2 r)^{ 7}β_{2} + (3  r)^{ 7}β_{3} + (4r)^{ 7}β_{4} + + (5  r)^{ 7}β_{5} + (6  r)^{ 7}β_{6}] D_{9} = 1/9![( r)^{9}α_{0} + (1  r)^{9}α_{1} + (2  r)^{9}α_{2} + (3  r)^{9}α_{3} + (4  r)^{9}α_{4} + (5  r)^{9}α_{5} + + (6  r)^{9}α_{6}]  1/8![(r)^{8}β_{0} + (1  r)^{8}β_{1} + (2 r)^{8}β_{2} + (3  r)^{8}β_{3} + (4  r)^{8}β_{4} + + (5  r)^{8}β_{5} + (6  r)^{8}β_{6}] 

Setting r = 3 and D_{q} = 0, q = 2, 3, 4, 5, 6, 7, 8 we have:
D_{2} = 1/2![3^{2}α_{0} +2^{2}α_{1} + α_{2} + α_{4} + 2^{2}α_{5} + 3^{2}α_{6}]   [3β_{0}  2β_{1}  β_{2} + β_{4} + 2β_{5} + 3β_{6}] = 0 D_{3} = 1/3![3^{3}α_{0}  2^{3}α_{1}  α_{2} + α_{4} + 2^{3}α_{5} + 3^{3}α_{6}]   1/2![3^{2}β_{0} + 2^{2}β_{1} + β_{2} + β_{4} + 2^{2}β_{5} + 3^{2}β_{6}] = 0 D_{4} = 1/4![3^{4}α_{0} + 2^{4}α_{1} + α_{2} + α_{4} + 2^{4}α_{5} + 3^{4}α_{6}]   1/3![3^{3}β_{0}  2^{3}β_{1}  β_{2} + β_{4} + 2^{3}β_{5} + 3^{3}β_{6}] = 0 D_{5} = 1/5![3^{5}α_{0}  2^{5}α_{1}  α_{2} + α_{4} + 2^{5}α_{5} + 3^{5}α_{6}]   1/4![3^{4}β_{0} + 2^{4}β_{1} + β_{2} + β_{4} + 2^{4}β_{5} + 3^{4}β_{6}] = 0 D_{6} = 1/6![3^{6}α_{0} + 2^{6}α_{1} + α_{2} + α_{4} + 2^{6}α_{5} + 3^{6}α_{6}]   1/5![3^{5}β_{0}  2^{5}β_{1}  β_{2} + β_{4} + 2^{5}β_{5} + 3^{5}β_{6}] = 0 D_{7} = 1/7![3^{7}α_{0}  2^{7}α_{1}  α_{2} + α_{4} + 2^{7}α_{5} + 3^{7}α_{6}]   1/6![3^{6}β_{0} + 2^{6}β_{1} + β_{2} + β_{4} + 2^{6}β_{5} + 3^{6}β_{6}] = 0 D_{8} = 1/8![3^{8}α_{0} + 2^{8}α_{1} + α_{2} + α_{4} + 2^{8}α_{5} + 3^{8}α_{6}]   1/7![3^{7}β_{0}  2^{7}β_{1}  β_{2} + β_{4} + 2^{7}β_{5} + 3^{7}β_{6}] = 0 

However, on inserting the values we have obtained for the α_{j} into these equations we have:
3β_{0}  2β_{1}  β_{2} + β_{4} + 2β_{5} + 3β_{6}] = 0 3^{2}β_{0} + 2^{2}β_{1} + β_{2} + β_{4} + 2^{2}β_{5} + 3^{2}β_{6} = 2/3[28 + 4ab  16(a + b)] 3^{3}β_{0}  2^{3}β_{1}  β_{2} + β_{4} + 2^{3}β_{5} + 3^{3}β_{6} = 0 3^{4}β_{0} + 2^{4}β_{1} + β_{2} + β_{4} + 2^{4}β_{5} + 3^{4}β_{6} = 2/5[244 + 4ab  64(a + b)] 3^{5}β_{0}  2^{5}β_{1}  β_{2} + β_{4} + 2^{5}β_{5} + 3^{5}β_{6} = 0 3^{6}β_{0} + 2^{6}β_{1} + β_{2} + β_{4} + 2^{6}β_{5} + 3^{6}β_{6} = 2/7[2188 + 4ab  256 (a + b)] 3^{7}β_{0}  2^{7}β_{1}  β_{2} + β_{4} + 2^{7}β_{5} + 3^{7}β_{6} = 0 
(9) 
We can satisfy the first, third, fifth and seventh of these equations if we choose: β_{2} = β_{4}, β_{1} = β_{5}, β_{0} = β_{6}. The remaining three equations give:
3^{2} β_{0} + 2^{2} β_{1} + β_{2} = 1/3[28 + 4ab  16(a+b)] 
(10) 
3^{4} β_{0} + 2^{4} β_{1} + β_{2} = 1/5[244 + 4ab  64(a+b)] 
(11) 
3^{6} β_{0} + 2^{6} β_{1} + β_{2} = 1/7[2188 + 4ab  256(a+b)] 
(12) 
From (10) we have:
β_{2} = 1/3[28 + 4ab  16(a+b)]  9β_{0}  4 β_{1} 
(13) 
Substituting (13) into (11):
81β_{0} + 16β_{1} + 1/3[28 + 4ab  16(a+b)]  9β_{0}  4 β_{1} = 1/5[244 + 4ab  64(a+b)] 72β_{0} + 12β_{1} = 1/5[244 + 4ab  64(a+b)]  1/3[28 + 4ab  16(a+b)] 72β_{0} + 12β_{1} = 1/15[592  8ab  112(a+b)] 
(14) 
Substituting (13) into (12):
729β_{0} + 64β_{1} + 1/3[28 + 4ab  16(a+b)]  9β_{0}  4 β_{1} = 1/7[2188 + 4ab  256(a+b)] 720β_{0} + 60β_{1} = 1/7[2188 + 4ab  256(a+b)]  1/3[28 + 4ab  16(a+b)] 720β_{0} + 60β_{1} = 1/21[6368 + 4ab  656(a+b)] 
(15) 
From (14) and (15), solving simultaneously to have:
β_{0} = 1/945[278 + 5ab + 16(a+b)] = β_{6} 
(16) 
Substituting (16) into (14), we have:
β_{1} = 1/150[160  8ab  76(a+b)] = β_{5} 
(17) 
Substituting (16) and (17) into (13) gives:
β_{2} = 1/105[62 + 167ab  272(a+b)] = β_{4} 
(18) 
Finally, solving D_{1 }= 0 gives:
β_{ 3} = 1/945[3008 + 5688ab  1328(a+b)] 
(19) 
We solve for the error constant, D_{9}:
D_{9} = 1/9![3^{9}α_{0}  2^{9}α_{1}  α_{2} + α_{4} + 2^{9}α_{5} + 3^{9}α_{6}]   1/8![3^{8}β_{0} + 2^{8}β_{1} + β_{2} + β_{4} + 2^{8}β_{5} + 3^{8}β_{6}] D_{9} = 1/907200[6016 + 736ab  8576(a+b)] 
(20) 
Since a = cosθ_{1}, b = cosθ_{2}, 0 < θ_{1} < π, 0 < θ_{2} < π, a and b are restricted to the range 1 < a < 1 and 1 < b < 1. Our choice of values for a and b is guided by the fact that we like to minimize the error constant as well as the need to develop a method that makes computation easier by reducing the number of operations involved.
The following values are therefore, assigned to the variables: a = 3/4, b = 1/3.
Hence the following values are obtained for the coefficients α_{i}, β_{i}:
α_{6} = + 1, α_{5} =  5/6, α_{4} = 0, α_{3} = 0, α_{2} = 0, α_{1} = 5/6, α_{0} =  1 β_{0} = 3401/11340 = β_{6}, β_{1} = 391/315 = β_{5}, β_{2} = 1117/1260 = β_{4}, β_{1} = 3848/2835 

Substituting the values of a and b into (20), the Error constant is: 0.002489711924.
And finally, we have the scheme:
y_{n+6} 5/6 y_{n+5} +5/6 y_{n+1}  y_{n} = = h[(3401/11340)f_{n+6} + (391/315)f_{n+5}  (1117/1260)f_{n+4} + (3848/2835)f_{n+3}   (1117/1260)f_{n+2} + (391/315)f_{n+1} + (3401/11340)f_{n}] 
(21) 
Test for Convergence
To prove that the scheme converges, it is sufficient for us to show that it is consistent as well as zerostable.
To prove consistency, we have from (3), (4) and (5):
ρ(1) = 1  5/6 + 5/6  1 = 0 
(22) 
ρ^{'}(1) = 6(1)  5(5/6) + 1(5/6)  0(0) = 2.666666667 
(23) 
σ(1) = 2.666666667 
(24) 
Equations (22)  (24) establish consistency.
For zerostability, we consider the following equation that represents the first characteristic polynomial of the scheme:
ρ(ξ) = ξ^{6}  5/6 ξ^{5} 5 5/6 ξ  1 
(25) 
The following are the roots of (25):
ξ_{1} = + 1, ξ_{2} =  1, ξ_{3} = [(3 + Ö7)/4]i, ξ_{4} = [(3  Ö7)/4]i, ξ_{5} = [(1 +2Ö2)/3]i, ξ_{6} = [(1  2Ö2)/3]i 

It is obvious that  ξ_{1} ≤ 1, I = 1, 2, 3, 4, 5, 6. Thus ζ_{i} = 1, 2, …, 6 satisfy the zero stability condition. Hence, we conclude that scheme 1 is convergent.
Numerical Examples
We use the scheme to solve some differential equations as follows. However, as with all kstep methods (k < 1) we need to use another method to calculate additional starting values. In this work, we decide to use the Fourth Order RungeKutta method to evaluate the starting values y_{n}, n = 0, 1, …, 5, since the RungeKutta methods constitute the most efficient method for generating starting values for linear multistep methods.
_{ The Fourth order RungeKutta method is given below:}
_{yn+1 = yn + h/6(k1 + 2k2 + 2k3 + k4)} _{k1 = f(xn, yn), k2 = f(xn +1/2h, yn + 1/2hk1), k3 = f(xn +1/2h, yn + 1/2hk2),} k_{4} = f(x_{n} +h, y_{n} + hk_{1}) 

We choose as our predictor, the fourth order AdamsBashforth method:
y_{n + 4} = y_{n + 3} +h/24(55f_{n+3}  59f_{n+2} + 37f_{n+1}  9f_{n}) 

Discussion of Results
In Table 1 the scheme produced errors as shown, within acceptable limits. The scheme is very accurate as exhibited in Table 2; there is no error at all. This is understandable because the solution of the differential equation is a polynomial of degree three, i.e., not greater than six, the step number.
Conclusion
Since we have used the numerical method to solve various differential equations we conclude that the method is accurate as it produced results which are comparable with those produced by other similar methods.
Table 1. Problem: F=X+Y; Y(0)=1; h=0.1; exact solution: Y(X)=2*EXP(X)X1
_{} 
Exact 
_{} 
Error 
0.0 
1.0000000000 
1.0000000000 
0.0000000000E+00 
0.1 
1.1103418362 
1.1103416667 
1.6948462878E07 
0.2 
1.2428055163 
1.2428051417 
3.7461895075E07 
0.3 
1.3997176152 
1.3997169941 
6.2102693121E07 
0.4 
1.5836493953 
1.5836484802 
9.1512116951E07 
0.5 
1.7974425414 
1.7974412772 
1.2642065803E06 
0.6 
2.0442376008 
2.0442361876 
1.4132161703E06 
0.7 
2.3275054149 
2.3275025204 
2.8945628872E06 
0.8 
2.6510818570 
2.6510783589 
3.4980860688E06 
0.9 
3.0192062223 
3.0192004614 
5.7608878645E06 
1.0 
3.4365636569 
3.4365566462 
7.0107547572E06 
Table 2. Problem: F=3X^{2}6X+5; Y(0)=1; h=0.1; Exact Solution: Y(X)=X^{3}3X^{2}+5X+1
x 
Exact 
Y(X) 
Error 
0.0 
1.0000000000 
1.0000000000 
0.0000000000E+00 
0.1 
1.4710000000 
1.4710000000 
0.0000000000E+00 
0.2 
1.8880000000 
1.8880000000 
0.0000000000E+00 
0.3 
2.2570000000 
2.2570000000 
0.0000000000E+00 
0.4 
2.5840000000 
2.5840000000 
0.0000000000E+00 
0.5 
2.8750000000 
2.8750000000 
0.0000000000E+00 
0.6 
3.1360000000 
3.1360000000 
0.0000000000E+00 
0.7 
3.3730000000 
3.3730000000 
0.0000000000E+00 
0.8 
3.5920000000 
3.5920000000 
0.0000000000E+00 
0.9 
3.7990000000 
3.7990000000 
0.0000000000E+00 
1.0 
4.0000000000 
4.0000000000 
0.0000000000E+00 
Reference
1. Lambert j. d., Computational Methods in Ordinary Differential Equations, John Wiley, 1973, p. 12 & 33& 38.