3-Point Implicit Block Multistep Method for the Solution of First Order ODEs

MALAYSIANMATHEMATICAL

SCIENCESSOCIETY http://math.usm.my/bulletin

3-Point Implicit Block Multistep Method for the Solution of First Order ODEs

1S. MEHRKANOON,²Z.A. MAJID,³M. SULEIMAN 4K.I. OTHMAN AND5Z.B. IBRAHIM

1,2,3,5Department of Mathematics, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia

4Department of Mathematics, Faculty of Information Technology and Quantitative Science, Universiti Technologi MARA, 40450 Shah Alam, Selangor, Malaysia

1[email protected]

Abstract. A new 3-point three step method is developed for solving system of first order ordinary differential equations (ODEs). This method at each step approximates the solution at three points simultaneously using variable step size. The method is in a simple form as Adams Moulton method with the specific aim of gaining efficiency. The Gauss-Seidel style is used for the implementation of the proposed method in PE(CE) mode. The stability regions of the method are discussed. Numerical results show that the proposed method is more efficient than some existent block method, in terms of accuracy, total number of steps and function calls and execution times

2010 Mathematics Subject Classification: Primary: 65L05, 65L06

Keywords and phrases: 3-point block method, ordinary differential equation, numerical method.

1. Introduction

In this paper, the initial value problems (IVPs) in the following form will be considered (1.1) y^′=f(x,y) y(a) =y₀ , a≤x≤b.

Block methods are one of the numerical methods which have been introduced by several au- thors such as [2, 7, 11, 12, 13]. One of the advantages of block methods is that at each step of the application we can obtain approximate solution in more than one point. The number of points depends on the construction of the block method. Primary study on this research area may be found in [4, 14]. The authors [7] have introduced 3-point diagonally implicit block method (3P1DI) using the Jacobi iteration approach during the implementation of the 3P1DI code. Zanariah et al. in [8, 9] have used Gauss-Seidel approach for the implementa- tion of 2-point block one step and 3-point block one step methods. In this paper, the same approach as in [8, 9] will be considered for the 3-point implicit block multistep method.

The proposed block method will approximate the solutions at three points simultaneously at each step using variable step size.

Received: June 22, 2009; Revised: October 4, 2010.

The method is derived by using Lagrange interpolation polynomial and the closest point in the interval will be considered for obtaining the corrector and predictor formula. There- fore, the approximation values of y_n+1, y_n+2and y_n+3are obtained by integrating (1.1) over the interval[xn,x_n+1],[x_n+1,x_n+2]and[x_n+2,x_n+3]respectively.

2. Derivation of 3PG method

In Figure 1, the approximations of y_n+1, y_n+2and y_n+3with step size h at the points x_n+1, x_n+2and x_n+3respectively are approximated simultaneously using four back values at the points x_n, x_n−1, x_n−2and x_n−3of the previous three steps where each interval has step size rh. The set of points{xn−5, . . . ,x_n}are used for deriving the predictor formula which its order is one less than the order of corrector formula.

The interpolation points involved for obtaining the corrector formula for y_n+1, y_n+2, y_n+3,

h h h

rh rh rh

qh qh

xn xn+1 xn+2 xn+3

xn−1

xn−2

xn−3

xn−4

xn−5

Figure 1. 3-point block multistep method

are{(xn−3,f_n−3), . . . ,(xn+3,f_n+3)}. The first point, yn+1, is derived by integrating (1.1) as follow,

Z x_n+1 xn

y^′dx= Z x_n+1

xn

f(x,y)dx.

Then

(2.1) y(x_n+1) =y(xn) +

Z _xn+1 xn

f(x,y)dx.

The function f(x,y)in (2.1) is approximated by the Lagrange polynomial which interpolates the set of mentioned points. Evaluating the integral using MAPLE gives the formula for the first point in terms of r as follows:

y_n+1=y_n+h

− 66+357r+532r²

7560r³(1+r)(1+3r)(2+3r)f_n−3+ 33+238r+399r²

840r³(1+r)(1+2r)(3+2r)f_n−2

− 66+595r+1596r²

840r³(1+r)(2+r)(3+r)f_n−1+33+357r+1463r²+2835r³

7560r³ f_n

+214+1680r+4389r²+3990r³

840(1+r)(1+2r)(1+3r) f_n+1−16+147r+462r²+525r³ 840(1+r)(2+r)(2+3r) f_n+2 +18+168r+539r²+630r³

7560(1+r)(3+r)(3+2r)f_n+3

.

The approximate value for second point, y_n+2, is derived by integrating (1.1) over the interval[xn+1,x_n+2]and approximating f using the same Lagrange polynomial. Then eval- uating the integral using MAPLE the formula of the second point in terms of r is obtained

as follows:

y_n+2=y_n+1+h

354+693r+308r²

7560r³(1+r)(1+3r)(2+3r)f_n−3− 59+154r+77r²

280r³(1+r)(1+2r)(3+2r)f_n−2 + 118+385r+308r²

280r³(1+r)(2+r)(3+r)f_n−1−177+693r+847r²+315r³

7560r³ f_n

+398+1680r+2233r²+910r³

280(1+r)(1+2r)(1+3r) f_n+1+368+1281r+1386r²+455r³ 280(1+r)(2+r)(2+3r) f_n+2

−402+1512r+1771r²+630r³ 7560(1+r)(3+r)(3+2r) f_n+3

.

For deriving the approximate value for the third point, y_n+3, interval[xn+2,x_n+3]is applied and after evaluating the integral, the formula for the third point is obtained as follows:

y_n+3=y_n+2+h

− 1746+2037r+532r²

7560r³(1+r)(1+3r)(2+3r)f_n−3+ 873+1358r+399r²

840r³(1+r)(1+2r)(3+2r)f_n−2

− 1746+3395r+1596r²

840r³(1+r)(2+r)(3+r)f_n−1+873+2037r+1463r²+315r³

7560r³ f_n

−2866+6720r+4851r²+1050r³

840(1+r)(1+2r)(1+3r) f_n+1+4744+11613r+8778r²+1995r³ 840(1+r)(2+r)(2+3r) f_n+2 +19338+42168r+28259r²+5670r³

7560(1+r)(3+r)(3+2r) f_n+3

.

The predictor formula can be derived similarly, but the interpolation points involved are {(x_n−5,f_n−5), . . . ,(xn,f_n)}.

The step size strategy in the code is the same as [7], the choice for next step size will be limited to half, double or as same as the current step size. In the developed code, when the next successful step size is doubled, the ratio r is 0.5 and q can be 0.5 or 0.25. If the next successful step size remain constant, r is 1 and q can be 1, 2 or 0.5. In case of step size failure, r is 2 and q is 2. Namely taking r=1 in the above mentioned formulas, we will have the following corrector formulas for the first, second and third point respectively:

yn+1=yn+h(−191 fn−3+1608 fn−2−6771 fn−1+37504 fn+30819 fn+1−2760 fn+2+271 fn+3)

60480 ,

yn+2=yn+1+h(271 fn−3−2088 fn−2+7299 fn−1−16256 fn+46989 fn+1+25128 fn+2−863 fn+3)

60480 ,

y_n+3=y_n+2

+h(−863 fn−3+6312 f_n−2−20211 f_n−1+37504 fn−46461 f_n+1+65112 f_n+2+19087 f_n+3)

60480 .

It can be seen that the derived corrector formulas when r=1, belong to the Generalized Linear Methods (GLM), see [1, 3]. In particular, the formula for approximating the third point, y_n+3coincides with the 6-step Adams method. In our code, an estimation of local truncation error is obtained by comparing the derived corrector formula of seventh order to the the same corrector formula of sixth order at the third point.

3. Implementation of 3PG method

The code starts by finding the initial points in the starting block for the method. Initially we use the Euler method to find the initial six starting points{x_n−j}⁰_j=5using constant h. Once we find the initial points for the first starting block, then we could implement the proposed block method until the end of the interval. The method is implemented in PE(CE)^t mode where P stands for an application of a predictor, E stands for an evaluation of a function f,and C stands for an application of a corrector. During the implementation, the iteration will involve the Gauss Seidel style. The t^thiteration process for the corrector formula is as follows:

C : y^t_n+1=y_n+

∑

αn+jf_n+j E : f(xn+1,y^t_n+1)

C : y^t_n+2=y^t_n+1+

∑

βn+jf_n+j E : f(xn+2,y^t_n+2)

C : y^t_n+3=y^t_n+2+

∑

γn+jf_n+j

E : f(xn+3,y^t_n+3) t=1,2, ...

where{αn+j}³_j=−3,{βn+j}³_j=−3and{γn+j}³_j=−3coincide with the coefficients of the above mentioned corrector formula for the first, second and third point respectively. The estimated value of y^t_n+1does not use the approach of Gauss-Seidel iteration. While for obtaining the approximate value of{y^t_n+i}³_i=2at points{xn+i}³_i=2respectively the Gauss Seidel iteration is involved.

In our code, the choice for next step size is limited to half, double or the same as the current step size. If the approximated solution at step i, has desired accuracy, i.e. it is acceptable, therefore the choice for next step will be doubled or the same as current step size which may be specified by step size controller. Otherwise the step size controller will allow the step size to become half. After a successful step, the step size increment is given by:

h_new=τ×h_old× T OL

LT E_k ¹_k

,

if h_new≥2×h_old then h_new=2×h_old else h_new=h_old

whereτ=0.8 is a safety factor. The purpose of having the safety factor is to increase the probability that the next step will be accepted. The algorithm when the step failure occurs is

h_new=1 2×h_old.

In the developed code, when the next step size is doubled, the ratio r is 0.5 and q can be 0.5 or 0.25, but if the next step size remains constant, r is 1 and q can be 1 or 2 or 0.5. In case of step size failure, r is 2, and q is 2.

4. Absolute stability

The absolute stability of 3PG method using a linear first order test problem

(4.1) y^′=f =λ^y

is discussed. The stability region is plotted when the step size ratio is constant, doubled and halved for the method. The test Equation (4.1) is substituted into the corrector formula of the 3PG method. Setting the determinant of the corrector formula written in matrix form to zero will give the stability polynomial. The stability polynomials of 3PG method at r=1,0.5,2 are as follow:

For r=1,

−547411 ¯h³

6048000 +63859 ¯h²

120960 −11947 ¯h 10080 +1

t⁶+

−796417 ¯h³

756000 −24337 ¯h² 25200 −659 ¯h

315 −1

t⁵

+

7879 ¯h³

63000 +13063 ¯h² 403200 +559 ¯h

2016

t⁴+

593 ¯h³

756000+101 ¯h² 75600

t³− ¯h³t² 6048000=0.

For r=2,

−7024939 ¯h³

55566000 +23795 ¯h² 37044 −507 ¯h

392 +1

t⁶

+

−4281161 ¯h³

7408800 −706621 ¯h²

604800 −163693 ¯h 94080 −1

t⁵+

523433 ¯h³

50803200+455029 ¯h²

13547520+3133 ¯h 94080

t⁴

+

5851 ¯h³

948326400+ 5809 ¯h² 474163200

t³− ¯h³t²

14224896000=0.

For r=0.5,

−3545993 ¯h³

55566000 +79099 ¯h² 185220 −843 ¯h

784 +1

t⁶+

−9633853 ¯h³

3704400 +306023 ¯h²

264600 −46723 ¯h 11760 −1

t⁵

+

252761 ¯h³

158760 +3685 ¯h²

1323 +3011 ¯h 1470

t⁴+

15412 ¯h³

231525 +21248 ¯h² 231525

t³− 512 ¯h³t² 3472875=0.

where ¯h=hλ and the stability regions are plotted in Figure 2.

In Figure 2, the stability region is inside the boundary of the dotted points. The stability region is larger when the step size is halved(r=2)compared to the step size being doubled (r=0.5)or constant(r=1). This is expected since the region should get larger with smaller step size.

Since the method is implemented in PE(CE)^tmode, we are interested to show the stabil- ity region of the predictor and corrector formula that involved for the most cases i.e constant and doubled step size. Substituting the predictor formula into the corrector formula and then applying (4.1) the stability polynomial will be obtained. Figure 3 (a) and (b), illustrate the stability region of the method for t=1 and 2 when r=q=0.5 and r=q=1 respectively where t is the number of iterations. Through our numerical experiences, the number of it- erations is at most 2. The stability region of method is inside the boundary of the dotted

r= 2

r= 1

r= 0.5

Figure 2. Stability region of 3PG method

points. It is observed from Figure 3, that the stability region is larger as the number of iter- ations is increased. The method developed is suitable for the numerical integration of non stiff and mildly stiff differential equations.

(a) (b)

t= 1 t= 2

t= 1 t= 2

Figure 3. Stability region of 3PG method when both predictor and corrector formula are taken into account, for t=1 and 2

5. Numerical results

In order to show the efficiency of our presented method, we consider three given problems to compare our computed solutions with the solutions obtained in [7]. The following notation is used in the tables:

Problem 1. Negative exponential problem. (Mildly stiff)

y^′=−0.5y, y(0) =1, [a,b] = [0,20].

The exact solution is given by y(x) =e^−0.5x.

Problem 2.( Hairer et al. [5]) A two-body orbit problem (Mildly stiff) y^′₁=y₃, y^′₂=y₄, y^′₃=−y₁

r³, y^′₄=−y₂ r³, r=

q y²₁+y²₂

TOL Tolerance

MTD Method Employed TS Total Successful Steps FS Total Failure Steps

MAXE Absolute value of the maximum error of the computed solution FN Total Function Calls

3PG Implementation of the three point block method using half Gauss Seidel iteration

3P1DI Implementation of the three point one block diagonally implicit method using Jacobi iteration in [7]

TIME The execution time taken in microsecond

y₁(0) =1−ε, y₂(0) =0, y₃(x) =0, y₄(x) = r1+ε

1−ε^{, [a,b] = [0,20].}

The values ofεtaken in solving the problem 2, are 0 and 0.5.

Problem 3.(Johnson and Barney [6]) Nonlinear non stiff Krogh’s problem

y^′_i=−βiy_i+y²_i, y_i(0) =−1, [a,b] = [0,20], i=1,2,3,4.

Case (I):β1=β2=0.2,β3=0.3,β4=0.4.

Case (II):β1=2,β2=3,β3=4,β4=5.

The exact solution is given by y_i(x) =βi/(1+cie^βix), c_i=−(1+βi).

It should be noted that in Problem 3, we consider two cases (I) and (II). Case (I) repre- sents Problem 3 with the same initial values as in [6] and the obtained numerical results are compared with the recorded results in [7]. Case (II) points out the artificial initial values, with emphasis onβ16=β2.

The error calculated are defined as (ei)t=

(yi)t−(y(xi)t) A+B(y(xi))t

where(y)t, is the t-th component of the approximate y. A=1, B=0 corresponds to the absolute error test, A=0, B=1 corresponds to the relative error test and finally A=1, B=1 corresponds to the mixed error test (see [10]). The mixed error test is utilized for all tested problems. The maximum error is defined as follows:

MAXE= max

1≤i≤P(max

1≤t≤N(ei)t)

where N is the number of equation in the system and P is the number of points which at each step their solutions will be estimated. In the code, we iterate the corrector to converge using the convergence criteria:

y^t+1_n+3−y^t_n+3

<0.1×TOL and t is the number of iterations.

The code is written in C language and executed on DYNIX/ptx operating system. The numerical results of the above tested problems are tabulated in Table 1–3.

Tables 1–3, show that the total number of steps and function calls taken by 3PG method are less than those in 3P1DI method in all tested problems. At larger tolerances, the maximum

Table 1. Numerical results of 3PG and 3P1DI methods for solving problem 1 TOL MTD TS FS MAXE FN TIME 10⁻² 3PG 18 0 1.10×10⁻⁴ 143 246

3P1DI 29 0 1.55×10⁻⁴ 215 525 10⁻⁴ 3PG 24 0 9.57×10⁻⁷ 206 184 3P1DI 38 0 1.16×10⁻⁵ 281 236 10⁻⁶ 3PG 32 0 6.64×10⁻⁹ 293 206 3P1DI 49 0 7.72×10⁻⁸ 398 293 10⁻⁸ 3PG 50 0 7.46×10⁻¹¹ 443 274 3P1DI 70 0 2.51×10⁻⁹ 620 419 10⁻¹⁰ 3PG 80 0 5.64×10⁻¹³ 749 442 3P1DI 95 0 1.42×10⁻¹⁰ 959 590

Table 2. Numerical results of 3PG and 3P1DI method for solving problem 2 TOL ε MTD TS FS MAXE FN TIME

10⁻² 0 3PG 24 0 4.17×10⁻² 293 877

3P1DI 45 0 5.17×10⁻² 353 1331

0.5 3PG 56 3 8.61×10⁻¹ 581 1423

3P1DI - - - - -

10⁻⁴ 3PG 45 0 5.16×10⁻⁵ 467 1206

3P1DI 57 0 8.57×10⁻⁵ 626 1559 0.5 3PG 95 6 2.16×10⁻³ 1055 2398

3P1DI - - - - -

10⁻⁶ 3PG 56 0 1.08×10⁻⁵ 701 1595

3P1DI 71 0 8.43×10⁻⁶ 983 2277 0.5 3PG 187 6 3.95×10⁻⁵ 1829 4851

3P1DI - - - - -

10⁻⁸ 3PG 119 0 1.40×10⁻⁸ 1325 3117

3P1DI 136 0 3.36×10⁻⁸ 1730 3978 0.5 3PG 330 3 1.68×10⁻⁷ 3032 8143

3P1DI - - - - -

10⁻¹⁰ 3PG 274 0 2.10×10⁻¹¹ 2414 6244 3P1DI 293 0 6.11×10⁻¹⁰ 3611 8422 0.5 3PG 614 6 6.53×10⁻⁹ 5645 1267

3P1DI - - - - -

Table 3. Numerical results of 3PG and 3P1DI methods for solving problem 3 TOL Case MTD TS FS MAXE FN TIME

10⁻² (I) 3PG 19 0 2.45×10⁻⁴ 179 560

3P1DI 28 0 1.38×10⁻⁴ 227 889 (II) 3PG 70 0 1.54×10⁻⁴ 722 1402

3P1DI - - - - -

10⁻⁴ (I) 3PG 28 0 1.14×10⁻⁶ 254 582

3P1DI 42 0 1.17×10⁻⁶ 320 754 (II) 3PG 75 0 1.95×10⁻⁶ 746 1448

3P1DI - - - - -

10⁻⁶ (I) 3PG 42 0 8.31×10⁻⁹ 392 763

3P1DI 59 0 3.99×10⁻⁸ 494 1062 (II) 3PG 97 0 1.72×10⁻⁸ 974 1795

3P1DI - - - - -

10⁻⁸ (I) 3PG 69 0 7.79×10⁻¹¹ 650 1250 3P1DI 88 0 1.02×10⁻⁹ 797 1598 (II) 3PG 115 0 2.43×10⁻¹⁰ 1073 2083

3P1DI - - - - -

10⁻¹⁰ (I) 3PG 118 0 9.82×10⁻¹³ 1112 2102 3P1DI 139 0 3.32×10⁻¹¹ 1316 2535 (II) 3PG 184 0 1.94×10⁻¹² 1691 3254

3P1DI - - - - -

error of 3PG is comparable or better compared to 3P1DI. The 3PG method has better accu- racy compared to 3P1DI method specially for finer tolerance. It should be noted that, the corrector formula of the 3P1DI for the first, second and third point are of the fourth, sixth and seventh order respectively. While in the new proposed method, all the corrector formula are of the same order i.e. seventh order. It is obvious that 3PG method is faster than 3P1DI method in all tolerance when solving three given problems.

6. Conclusion

In this paper, we proposed a 3-point three-step method for solving system of initial value problems of ODEs. From the numerical results we can conclude that the method has supe- riority in terms of total number of steps, maximum error and function calls over the 3P1DI code in [7].

References

[1] L. Brugnano and D. Trigiante, Solving Differential Problems by Multistep Initial and Boundary Value Meth- ods, Stability and Control: Theory, Methods and Applications, 6, Gordon and Breach, Amsterdam, 1998.

[2] K. Burrage, Efficient block predictor-corrector methods with a small number of corrections, J. Comput. Appl.

Math. 45 (1993), no. 1-2, 139–150.

[3] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, 2003.

[4] O.S. Fatunla, Block methods for second order ODEs, Intern. J. Computer Math. 41 (1990), 55–63.

[5] E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations. I, second edition, Springer Series in Computational Mathematics, 8, Springer, Berlin, 1993.

[6] A. I. Johnson and J. R. Barney, Numerical solution of large systems of stiff ordinary differential equations in a modular simulation framework, in Numerical Methods for Differential Systems, 97–124, Academic Press, New York, 1976.

[7] Z. A. Majid and M. Suleiman, Variable steps three point diagonally implicit block method for solving ordinary differential equations, Int. J. Appl. Math. Anal. Appl. 3 (2008), no. 1, 79–90.

[8] Z.A. Majid, M. Suleiman, F. Ismail and M. Othman, 2-point implicit block one-step method half Gauss-Seidel for solving first order ordinary differential equations, Mathematika 19 (2003), no. 2, 91–100.

[9] Z. Abdul Majid, M. B. Suleiman and Z. Omar, 3-point implicit block method for solving ordinary differential equations, Bull. Malays. Math. Sci. Soc. (2) 29 (2006), no. 1, 23–31.

[10] F. Mazzia and C. Magherini, Test Set for Initial Value Problem Solvers, Department of Mathematics, Univer- sity of Bari, February 2008. Available at: http://www.dm.uniba.it/∼testset

[11] W. E. Milne, Numerical Solution of Differential Equations, Wiley, New York, 1953.

[12] J. B. Rosser, A Runge-Kutta for all seasons, SIAM Rev. 9 (1967), 417–452.

[13] L. F. Shampine and H. A. Watts, Block implicit one-step methods, Math. Comp. 23 (1969), 731–740.

[14] P. B. Worland, Parallel methods for the numerical solution of ordinary differential equations, IEEE Trans.

Computers C-25 (1976), no. 10, 1045–1048.