Derivation and Application of a Linear Multistep Numerical Scheme

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: One-step and multistep methods. However, in this work, our focus is on a class of multistep methods known as the Linear Multi-step Methods. We thus use Taylor series expansion to derive a linear multistep method of order eight. The method is tested for consistency and zero-stability 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₀) = y₀, t > t₀

(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 α₀ and β₀ 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 zero-stable [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 zero-stable 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 6-step 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, ξ₃ = e^iθ1, ξ₃ = e^-iθ1, ξ₃ = e^iθ2, ξ₃ = e^-iθ2

Hence:

ρ(ξ) = (ξ - 1)(ξ + 1)(ξ - e^iθ1)(ξ - e^-iθ1)(ξ - e^iθ2) (ξ-e^-iθ2)
α₆ = +1, α₅ = -2(a+b), α₄ = (4ab + 1), α₃ = 0, α₂ = -(4ab+1), α₁ = 2(a+b), α₀ = -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₀ = α₀ + α₁ + α₂ + α₃ + α₄ + α₅ + α₆

D₁ = -rα₀ + (1 - r)α₁ + (2 - r)α₂ + (3 - r)α₃ + (4 - r)α₄ + (5 - r)α₅ + (6 - r)α₆ -

- (β₀ + β₁ + β₂ + β₃ + β₄ + β₅ + β₆)

D₂ = 1/2![(- r)²α₀ + (1 - r)²α₁ + (2 - r)²α₂ + (3 - r)²α₃ + (4 - r)²α₄ + (5 - r)²α₅ +

+ (6 - r)²α₆] - [-rβ₀ + (1 - r)β₁ + (2- r)β₂ + (3 - r)β₃ + (4 - r)β₄ + (5 - r)β₅ +

+ (6 - r)β₆]

D₈ = 1/8![(- r)⁸α₀ + (1 - r)⁸α₁ + (2 - r)⁸α₂ + (3 - r)⁸α₃ + (4 - r)⁸α₄ + (5 - r)⁸α₅ +

+ (6 - r)⁸α₆] - 1/7![(-r)⁷β₀ + (1 - r)⁷β₁ + (2- r)⁷β₂ + (3 - r)⁷β₃ + (4-r)⁷β₄ +

+ (5 - r)⁷β₅ + (6 - r)⁷β₆]

D₉ = 1/9![(- r)⁹α₀ + (1 - r)⁹α₁ + (2 - r)⁹α₂ + (3 - r)⁹α₃ + (4 - r)⁹α₄ + (5 - r)⁹α₅ +

+ (6 - r)⁹α₆] - 1/8![(-r)⁸β₀ + (1 - r)⁸β₁ + (2- r)⁸β₂ + (3 - r)⁸β₃ + (4 - r)⁸β₄ +

+ (5 - r)⁸β₅ + (6 - r)⁸β₆]

Setting r = 3 and D_q = 0, q = 2, 3, 4, 5, 6, 7, 8 we have:

D₂ = 1/2![3²α₀ +2²α₁ + α₂ + α₄ + 2²α₅ + 3²α₆] -

- [-3β₀ - 2β₁ - β₂ + β₄ + 2β₅ + 3β₆] = 0

D₃ = 1/3![-3³α₀ - 2³α₁ - α₂ + α₄ + 2³α₅ + 3³α₆] -

- 1/2![3²β₀ + 2²β₁ + β₂ + β₄ + 2²β₅ + 3²β₆] = 0

D₄ = 1/4![3⁴α₀ + 2⁴α₁ + α₂ + α₄ + 2⁴α₅ + 3⁴α₆] -

- 1/3![-3³β₀ - 2³β₁ - β₂ + β₄ + 2³β₅ + 3³β₆] = 0

D₅ = 1/5![-3⁵α₀ - 2⁵α₁ - α₂ + α₄ + 2⁵α₅ + 3⁵α₆] -

- 1/4![3⁴β₀ + 2⁴β₁ + β₂ + β₄ + 2⁴β₅ + 3⁴β₆] = 0

D₆ = 1/6![3⁶α₀ + 2⁶α₁ + α₂ + α₄ + 2⁶α₅ + 3⁶α₆] -

- 1/5![-3⁵β₀ - 2⁵β₁ - β₂ + β₄ + 2⁵β₅ + 3⁵β₆] = 0

D₇ = 1/7![-3⁷α₀ - 2⁷α₁ - α₂ + α₄ + 2⁷α₅ + 3⁷α₆] -

- 1/6![3⁶β₀ + 2⁶β₁ + β₂ + β₄ + 2⁶β₅ + 3⁶β₆] = 0

D₈ = 1/8![3⁸α₀ + 2⁸α₁ + α₂ + α₄ + 2⁸α₅ + 3⁸α₆] -

- 1/7![-3⁷β₀ - 2⁷β₁ - β₂ + β₄ + 2⁷β₅ + 3⁷β₆] = 0

However, on inserting the values we have obtained for the α_j into these equations we have:

-3β₀ - 2β₁ - β₂ + β₄ + 2β₅ + 3β₆] = 0

3²β₀ + 2²β₁ + β₂ + β₄ + 2²β₅ + 3²β₆ = 2/3[28 + 4ab - 16(a + b)]

-3³β₀ - 2³β₁ - β₂ + β₄ + 2³β₅ + 3³β₆ = 0

3⁴β₀ + 2⁴β₁ + β₂ + β₄ + 2⁴β₅ + 3⁴β₆ = 2/5[244 + 4ab - 64(a + b)]

-3⁵β₀ - 2⁵β₁ - β₂ + β₄ + 2⁵β₅ + 3⁵β₆ = 0

3⁶β₀ + 2⁶β₁ + β₂ + β₄ + 2⁶β₅ + 3⁶β₆ = 2/7[2188 + 4ab - 256 (a + b)]

-3⁷β₀ - 2⁷β₁ - β₂ + β₄ + 2⁷β₅ + 3⁷β₆ = 0

(9)

We can satisfy the first, third, fifth and seventh of these equations if we choose: β₂ = β₄, β₁ = β₅, β₀ = β₆. The remaining three equations give:

3² β₀ + 2² β₁ + β₂ = 1/3[28 + 4ab - 16(a+b)]	(10)
3⁴ β₀ + 2⁴ β₁ + β₂ = 1/5[244 + 4ab - 64(a+b)]	(11)
3⁶ β₀ + 2⁶ β₁ + β₂ = 1/7[2188 + 4ab - 256(a+b)]	(12)

From (10) we have:

β₂ = 1/3[28 + 4ab - 16(a+b)] - 9β₀ - 4 β₁

(13)

Substituting (13) into (11):

81β₀ + 16β₁ + 1/3[28 + 4ab - 16(a+b)] - 9β₀ - 4 β₁ = 1/5[244 + 4ab - 64(a+b)]

72β₀ + 12β₁ = 1/5[244 + 4ab - 64(a+b)] - 1/3[28 + 4ab - 16(a+b)]

72β₀ + 12β₁ = 1/15[592 - 8ab - 112(a+b)]

(14)

Substituting (13) into (12):

729β₀ + 64β₁ + 1/3[28 + 4ab - 16(a+b)] - 9β₀ - 4 β₁ = 1/7[2188 + 4ab - 256(a+b)]

720β₀ + 60β₁ = 1/7[2188 + 4ab - 256(a+b)] - 1/3[28 + 4ab - 16(a+b)]

720β₀ + 60β₁ = 1/21[6368 + 4ab - 656(a+b)]

(15)

From (14) and (15), solving simultaneously to have:

β₀ = 1/945[278 + 5ab + 16(a+b)] = β₆

(16)

Substituting (16) into (14), we have:

β₁ = 1/150[160 - 8ab - 76(a+b)] = β₅

(17)

Substituting (16) and (17) into (13) gives:

β₂ = 1/105[62 + 167ab - 272(a+b)] = β₄

(18)

Finally, solving D₁= 0 gives:

β₃ = 1/945[3008 + 5688ab - 1328(a+b)]

(19)

We solve for the error constant, D₉:

D₉ = 1/9![-3⁹α₀ - 2⁹α₁ - α₂ + α₄ + 2⁹α₅ + 3⁹α₆] -

- 1/8![3⁸β₀ + 2⁸β₁ + β₂ + β₄ + 2⁸β₅ + 3⁸β₆]

D₉ = -1/907200[6016 + 736ab - 8576(a+b)]

(20)

Since a = cosθ₁, b = cosθ₂, 0 < θ₁ < π, 0 < θ₂ < π, 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:

α₆ = + 1, α₅ = - 5/6, α₄ = 0, α₃ = 0, α₂ = 0, α₁ = 5/6, α₀ = - 1

β₀ = 3401/11340 = β₆, β₁ = 391/315 = β₅, β₂ = -1117/1260 = β₄, β₁ = 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 zero-stable.

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 zero-stability, we consider the following equation that represents the first characteristic polynomial of the scheme:

ρ(ξ) = ξ⁶ - 5/6 ξ⁵ 5 5/6 ξ - 1

(25)

The following are the roots of (25):

ξ₁ = + 1, ξ₂ = - 1, ξ₃ = [(3 + Ö7)/4]i, ξ₄ = [(3 - Ö7)/4]i, ξ₅ = [(-1 +2Ö2)/3]i,

ξ₆ = [(-1 - 2Ö2)/3]i

It is obvious that | ξ₁| ≤ 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 k-step methods (k < 1) we need to use another method to calculate additional starting values. In this work, we decide to use the Fourth Order Runge-Kutta method to evaluate the starting values y_n, n = 0, 1, …, 5, since the Runge-Kutta methods constitute the most efficient method for generating starting values for linear multistep methods.

_{The Fourth order Runge-Kutta
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₄ = f(x_n +h, y_n + hk₁)

We choose as our predictor, the fourth order Adams-Bashforth 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)-X-1

	Exact		Error
0.0	1.0000000000	1.0000000000	0.0000000000E+00
0.1	1.1103418362	1.1103416667	1.6948462878E-07
0.2	1.2428055163	1.2428051417	3.7461895075E-07
0.3	1.3997176152	1.3997169941	6.2102693121E-07
0.4	1.5836493953	1.5836484802	9.1512116951E-07
0.5	1.7974425414	1.7974412772	1.2642065803E-06
0.6	2.0442376008	2.0442361876	1.4132161703E-06
0.7	2.3275054149	2.3275025204	2.8945628872E-06
0.8	2.6510818570	2.6510783589	3.4980860688E-06
0.9	3.0192062223	3.0192004614	5.7608878645E-06
1.0	3.4365636569	3.4365566462	7.0107547572E-06

Table 2. Problem: F=3X²-6X+5; Y(0)=1; h=0.1; Exact Solution: Y(X)=X³-3X²+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.