Hybrid Numerical Method with Block Extension for Direct Solution of Third Order Ordinary Differential Equations

Abstract

This paper focuses on the development of a hybrid method with block extension for direct solution of initial value problems (IVPs) of general third-order ordinary differential equations. Power series was used as the basis function for the solution of the IVP. An approximate solution from the basis function was interpolated at some selected off-grid points while the third derivative of the approximate solution was collocated at all grid and off-grid points to generate a system of linear equations for the determination of the unknown parameters. The derived method was tested for consistency, zero stability, convergence and absolute stability. The method was implemented with five test problems including the Genesio equation to confirm its accuracy and usability. The rate of convergence (ROC) reveals that the method is consistent with the theoretical order of the proposed method. Comparison of the results with some existing methods shows the superiority of the accuracy of the method.

Share and Cite:

Duromola, M. and Momoh, A. (2019) Hybrid Numerical Method with Block Extension for Direct Solution of Third Order Ordinary Differential Equations. American Journal of Computational Mathematics, 9, 68-80. doi: 10.4236/ajcm.2019.92006.

1. Introduction

The focus of this article is to find an approximate solution on a given interval to third order initial value problems (IVP) of the type

y ( x ) = f ( x , y ( x ) , y ( x ) ) , y ( x ) = α a , y ( x ) = α b , y ( x ) = α c (1)

where x [ a , b ] and y ( x ) , f ( x , y ( x ) , y ( x ) , y ( x ) ) n . In recent time, direct numerical solution of (1) without reduction to equivalent first-order initial value problems (see: [1] [2] [3] [4] [5] ) has become subject of research by several authors. This method was extensively discussed in ( [6] [7] [8] [9] [10] ) to mention but a few, they developed Linear Multistep Method (LMM) which mode of implementation is Predictor-Corrector form for the solution of initial value problems of ordinary differential equations of the type (1). As reported by [10] [11] , the major drawback of this approach of implementation is that the methods are not self-starting and thus required the development of predictors which are usually of lower order, hence reducing the accuracy of the methods.

Recently, in order to remove the difficulties usually encountered by adopting this mode of solution, researchers ( [1] [11] - [22] ) have proposed direct methods other than Predictor-Corrector methods whose modes of implementation are in block-by-block manner which was first introduced by Milne [23] as a starting step for predictor-corrector. This monumental success has greatly removed the burden of developing predictors and hence resulted in methods of uniform orders that yielded more accurate results. The block-by-block technique has also made it easier to handle the general type of (1) which has been a major concern in the past years.

The quest for numerical methods with better accuracy has also led to the introduction of hybrid linear multistep methods which have recorded high success since its introduction. These successes motivated us to propose a hybrid method with block extension for the solution of (1).

In the next section, we discuss in detail the derivation of the proposed method with its implementation in block mode, followed by analysis of the proposed method to establish the numerical stability, numerical example to demonstrate the efficiency advantages of the proposed method and subsequently. Conclusion was drawn on the performance of the proposed method when applied to solve the numerical examples.

2. Mathematical Formulation

In order to obtain a numerical formula for the approximate solution of (1), the function

y ( x ) = v ( c + i ) 1 g v x v (2)

is considered as the basis where x is continuous within the interval [ a , b ] , c and i denote collocation and interpolation points respectively. Variable g v ’s are coefficients to be determined. The third derivative of (2) equated to (1) is given as

v ( c + i ) 1 v ( v 1 ) ( v 2 ) g v x v 3 = f ( x , y , y , y ) (3)

Evaluating (2) at x = x n + v , v = 2 8 , 4 8 , 6 8 , (3) at x = x n + v , v = 0 ( 2 8 ) 1 using τ = x x n + k h yield the following interpolation and collocation matrix

X A = B (4)

where

A = ( g 0 g 1 g 2 g 3 g 4 g 5 g 6 g 7 ) , B = ( y n + 2 8 y n + 4 8 y n + 6 8 f n f n + 2 8 f n + 4 8 f n + 6 8 f n + 1 ) ,

X = ( 1 1 4 1 16 1 64 1 256 1 1024 1 4096 1 16384 1 1 2 1 4 1 8 1 16 1 32 1 64 1 128 1 3 4 9 16 27 64 81 256 243 1024 729 4096 2187 16384 0 0 0 6 0 0 0 0 0 0 0 6 6 15 4 15 8 105 128 0 0 0 6 12 15 15 105 8 0 0 0 6 18 135 4 405 8 8505 128 0 0 0 6 24 60 120 210 ) .

where f n + v = f ( x n + v , y n + v , y n + v ) , y n + v y ( x n + v ) . Solving the matrix Equation (4) for coefficients g v ’s and substituting into (2) yields after some simplification the continuous method

y ¯ ( x ) = ζ 1 4 y n + 1 4 + ζ 1 2 y n + 1 2 + ζ 3 4 y n + 3 4 + h 3 ( Θ 0 ( x ) f n + Θ 1 4 ( x ) f n + 1 4 + Θ 1 2 ( x ) f n + 1 2 + Θ 3 4 ( x ) f n + 3 4 + Θ 1 ( x ) f n + 1 ) (5)

with the following coefficients:

ζ 1 4 = 8 τ 2 10 τ + 3

ζ 1 2 = 16 τ 2 + 16 τ 3

ζ 3 4 = 8 τ 2 6 τ + 1

Θ 0 = 1 322560 ( 4 τ 1 ) ( 2 τ 1 ) ( 4 τ 3 ) ( 512 τ 4 1472 τ 3 + 1360 τ 2 400 τ + 7 )

Θ 1 4 = 1 80640 ( 4 τ 3 ) ( 2 τ 1 ) ( 4 τ 1 ) ( 512 τ 4 1248 τ 3 + 688 τ 2 + 258 τ 203 )

Θ 1 2 = 1 53760 ( 4 τ 3 ) ( 2 τ 1 ) ( 4 τ 1 ) ( 512 τ 4 1024 τ 3 + 240 τ 2 + 272 τ + 147 )

Θ 3 4 = 1 80640 ( 4 τ 3 ) ( 2 τ 1 ) ( 4 τ 1 ) ( 512 τ 4 800 τ 3 + 16 τ 2 + 62 τ + 7 )

Θ 1 = 1 322560 ( 4 τ 1 ) ( 2 τ 1 ) ( 4 τ 3 ) ( 512 τ 4 576 τ 3 + 16 τ 2 + 48 τ + 7 )

Evaluating the continuous scheme (5) at τ = 0 , 1 and its first and second derivatives at τ = 0 yield two discrete, one first and second derivatives schemes. These can be represented in a block matrix finite difference form as

ϒ 0 Y m , 0 = ϒ 1 Y m , 1 + h ϒ 2 Y m , 2 + h 2 ϒ 3 Y m , 3 + h 3 H ¯ F m , 0 + h 3 H F m , 1 (6)

where T denotes the transpose,

H = ( 107 64512 103 107520 43 107520 47 645120 83 5040 1 168 13 5040 19 40320 1863 35840 243 35840 45 7168 81 71680 34 315 1 210 2 105 1 504 ) ;

H ¯ = ( 0 0 0 113 71680 0 0 0 331 40320 0 0 0 331 40320 0 0 0 31 840 ) ; ϒ 0 = ( 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ) ;

ϒ 1 = ( 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 ) ; ϒ 2 = ( 0 0 0 1 4 0 0 0 1 2 0 0 0 3 4 0 0 0 1 ) ; ϒ 3 = ( 0 0 0 1 32 0 0 0 1 8 0 0 0 9 32 0 0 0 1 2 )

Vectors Y m ,0 , Y m ,1 , Y m ,2 , Y m ,3 , F m ,0 and F m ,1 are

Y m , 0 = ( y n + 1 4 , y n + 1 2 , y n + 3 4 , y n + 1 ) T , Y m , 1 = ( y n 1 4 , y n 1 2 , y n 3 4 , y n ) T , Y m , 2 = ( y n 1 4 , y n 1 2 , y n 3 4 , y n ) T , Y m , 3 = ( y n 1 4 , y n 1 2 , y n 3 4 , y n ) T , F m , 0 = ( f n 1 4 , f n 1 2 , f n 3 4 , f n ) T , F m , 1 = ( f n + 1 4 , f n + 1 2 , f n + 3 4 , f n + 1 ) T .

3. Analysis of the Proposed Method

This section presents analysis of the proposed method (6) vis-a-vis the order, consistency, zero-stability and convergence.

The linear operator associated with the block method (6) is

L [ y ( x ) ; h ] = ϒ 0 Y m , 0 ( ϒ 1 Y m , 1 + h ϒ 2 Y m , 2 + h 2 ϒ 3 Y m , 3 + h 3 H ¯ F m , 0 + h 3 H F m , 1 ) (7)

where y ( x ) is an arbitrary function which is continuously differentiable on [ a , b ] . Following Lambert [5] and Fatunla [20] , the term in (6) can be written as a Taylor series expansion about the point x to obtain the expression,

L [ y ( x ) ; h ] = c ¯ 0 y ( x ) + c ¯ 1 h y ( x ) + c ¯ 2 h 2 y ( x ) + + c ¯ p h p y p ( x ) + , (8)

where the constant coefficients c ¯ p , p = 0 , 1 , 2 , are given as follows:

c ¯ 0 = ζ k + ζ u + ζ v + ζ w

c ¯ 1 = k ζ k + u ζ u + v ζ v + w ζ w

c ¯ 2 = 1 2 ! ( k 2 ζ k + u 2 ζ u + v 2 ζ v + w 2 ζ w )

c ¯ 3 = 1 3 ! ( k 3 ζ k + u 3 ζ u + v 3 ζ v + w 3 ζ w ) 1 0 ! ( Θ 0 + Θ k + Θ u + Θ v + Θ w )

c ¯ 4 = 1 p ! ( k 4 ζ k + u 4 ζ u + v 4 ζ v + w 4 ζ w ) 1 ( p 3 ) ! ( k Θ k + u Θ u + v Θ v + w Θ w )

c ¯ p = 1 p ! ( k p ζ k + u p ζ u + v p ζ v + w p ζ w ) 1 ( p 3 ) ! ( k p 3 Θ k + u p 3 Θ u + v p 3 Θ v + w p 3 Θ w ) , p = 4,5,

Going by Lambert [5] , the mutistep collocation method (6) has order p if

L [ y ( x ) ; h ] = 0 ( h p + 1 ) , c ¯ 0 = c ¯ 1 = = c ¯ p = 0 , c ¯ p + 3 0

Therefore c ¯ p + 3 is the error constant and c ¯ p + 3 h p + 3 y p + 3 ( x n ) is the principal local truncation error at point x n . The order of the proposed method (6) and the corresponding error constant are as reported in Table 1.

Definition 1 (consistency).

The proposed method (6) is said to be consistent if the order of method is greater than or equal to one, that is if p 1 . In addition to

1) ρ ( 1 ) = 0 and

Table 1. Order of accuracy and error constant of the proposed method.

2) ρ ( 1 ) = σ ( 1 ) where ρ ( z ) and σ ( z ) are 1st and 2nd characteristics polynomial respectively.

Definition 2 (Zero-stability).

The block method (6) is said to be zero-stable if the roots

ρ ( z ) = det [ i = 0 k ϒ ( i ) z k i ] (9)

satisfies | z i | 1 , i = 1 , , k and the roots with | z i | = 1 , the multiplicity must not exceed one. Applying (9) to the proposed method (6) yields the following

ρ ( z ) = | [ z 0 0 0 0 z 0 0 0 0 z 0 0 0 0 z ] [ 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 ] | = | [ z 0 0 1 0 z 0 1 0 0 z 1 0 0 0 z 1 ] | = z 3 ( z 1 ) (10)

This result shows that the method is zero-stable.

Definition 3 (convergence).

The necessary and sufficient condition for the proposed method (6) to be convergent are that it must be consistent and zero-stable according to Dahlquist see [5]. Hence, by definitions 4 and 5 the method is convergent.

Stability Domain of the Proposed Method

In order to study the stability domain of the proposed method (6), the test equations

y = λ y (11)

y = λ 2 y (12)

and

y = λ 3 y (13)

are applied to the block method (6) with z = λ y and η represents the roots of the first characteristic polynomial of the block method (6). This is then reformulated as a general linear method as discussed in [24]. The partition ( s 1 + s 2 ) × ( s 1 + s 2 ) matrix is expressed in the form

[ Y Y i + 1 ] = [ A V B U ] [ h 3 f ( y ) Y i 1 ] (14)

where

A = [ 0 0 0 0 0 113 71680 107 64512 103 107520 43 107520 47 645120 331 40320 83 5040 1 168 13 5040 19 40320 331 40320 1863 35840 243 35840 45 7168 81 71680 31 840 34 315 1 210 2 105 1 504 ] ,

B = [ 113 71680 107 64512 103 107520 43 107520 47 645120 31 840 34 315 1 210 2 105 1 504 ]

U = [ 0 1 0 1 ] , V = [ 0 0 0 0 0 1 1 1 1 1 ] T , f ( y ) = [ f n f n + 1 4 f n + 1 2 f n + 3 4 f n + 1 ] T ,

Y i + 1 = [ y n + 1 4 y n + 1 ] , Y i 1 = [ y n + 1 4 y n ]

By solving the stability function

p ( η , z ) = | η I ( v + z B m ( I z A m ) 1 U m ) | (15)

yields the polynomial

p ( η , z ) = η 315 [ 2835 η z 4 + 4737600 η z 3 + 50591 z 4 + 417312000 η z 2 134158720 z 3 145871596800 z 2 + 124853944320000 η 20808990720000 z 124853944320000 9 z 4 + 15040 z 3 + 1324800 z 2 + 396361728000 ] (16)

(16) and its derivatives are then plotted in the MATLAB environment given the stability region displayed in Figure 1.

Definition 4 (Lambert and Watson [25] ).

Method (6) is P-stable if its interval of periodicity is ( 0, ) . It is clearly

Figure 1. This figure depicts the region of absolute stability of the proposed method generated by plotting Equation (16).

shown in Figure 1 that the block method (6) is P-stable.

4. Numerical Example

Example I

The first numerical example to be considered is the oscillatory problem

y = 3 sin x , y ( 0 ) = 1 , y ( 0 ) = 0 , y ( 0 ) = 2 , h = 0.1

with the theoretical solution

y ( x ) = 3 cos x + x 2 2 2.

This example was solved by [11] [15] [18] [26]. The numerical solution, exact solution and absolute error generated by the proposed method when applied to example I are as presented in Table 2. The last column of the Table shows the errors generated by method in [26] when applied to example I. It is obvious from the table that the proposed method is better in term of accuracy when compared with the method in [26].

Example II

The second example considered is the special third order problem

y = e x , y ( 0 ) = 3 , y ( 0 ) = 1 , y ( 0 ) = 5 , h = 0.1

with the theoretical solution

y ( x ) = y ( x ) = 2 + 2 x 2 + e x

Source: [26]. The solution curve is shown in Figure 2.

Example III

Another example considered is a general third order problem

y + 2 y 9 y 18 y = 18 x 2 18 x + 22 , y ( 0 ) = 2 , y ( 0 ) = 8 , y ( 0 ) = 12 , h = 0.1

with the theoretical solution

Table 2. Results of Example I solved with the proposed method.

Figure 2. Solution curve obtained by our method and the exact solution of example II.

y ( x ) = 2 e 3 x + e 2 x + x 2 1.

The theoretical solution at x = 1 is y ( 1 ) 40.0357385631387227899630 . The errors were obtained at x = 1 using our method at a fixed step-size h = 0.1 ; 0.05 ; 0.025 ; 0.0125 ; 0.00625 . The numerical results are compared with those of [27]. For this example, the maximum error was compared with those reported in [18] in Table 3 for h = 0.01 and it was observed that our method perform better. The ROC, computed solutions and maximum error of the proposed method are reported in Table 4. The Table also shows the performance of our method as compared with method in [27].

Example IV

General nonlinear third order equation

y = y ( 2 x y + y ) , y ( 0 ) = 1 , y ( 0 ) = 1 / 2 , y ( 0 ) = 0.1 , h = 0.1

with the theoretical solution

y ( x ) = 1 + 1 2 log ( 2 + x 2 x )

is also considered. Source: [28]. Figure 3 is the graph of the solution of this problem.

Application to solve nonlinear Genesio equation

The chaotic Genesio equation reported in [17] given as

y = α y β y + f ( y (x))

with

f ( y ( x ) ) = γ y ( x ) + y 2 (x)

y ( 0 ) = 0.2 , y ( 0 ) = 0.3 , y ( 0 ) = 0.1 , x [ a , b ]

where α = 1.2 , β = 2.92 and γ = 6 are the positive constants that satisfied

α β < γ

Table 3. Results of Example III solved with the proposed method.

Table 4. Results of Example III solved with the proposed method.

Figure 3. Solution curve of example IV obtained by the proposed method.

for the solution to exist. The solution of this problem is presented in Figure 4(a) and Figure 4(b) in the intervals [ 0,10 ] and [ 0,100 ] respectively.

5. Conclusion

In this work, hybrid method with block extension for the direct solution of third order ordinary differential equations has been proposed. Numerical examples are considered to demonstrate the efficiency advantage of the method especially the Genesio equation which is chaotic in nature. The analysis, stability and numerical examples revealed that the proposed method is efficient for direct solution of third order ordinary differential equations.

Figure 4. Solution curve of the Genesio equation for α = 1.2 , β = 2.92 and γ = 6 for step size h = 0.1 within [0, N]: (a) N = 10 and (b) N = 100.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Henrici, P. (1962) Discrete Variable Methods in Ordinary Differential Equations. John Wiley & Sons, New York.
[2] Butcher, J.C. (1965) Modified Multistep Method for the Numerical Integration of Ordinary Differential Equations. Journal of the ACM, 12, 124-135.
https://doi.org/10.1145/321250.321261
[3] Gear, C.W. (1965) Hybrid Multistep Method for Initial Value in Ordinary Differential Equations. SIAM Journal on Numerical Analysis, 2, 69-86.
[4] Gragg, W.B. and Stetter, H.J. (1964) Generalized Multistep Predictor-Corrector Methods. Journal of the ACM, 11, 188-209. https://doi.org/10.1145/321217.321223
[5] Lambert, J.D. (1973) Computational Methods in Ordinary Differential Equations. John Willy, and Sons, New York.
[6] Awoyemi, D.O. (1996) An Efficient Two Step Numerical Integrator for General Second Order Ordinary Differential Equations. Abacus Journal of the Mathematical Association of Nigeria, 24, 31-43.
[7] Awoyemi, D.O. (1999) A Class of Continuous Methods for General Second Order Initial Value Problems in Ordinary Differential Equations. International Journal of Computer Mathematics, 72, 29-39. https://doi.org/10.1080/00207169908804832
[8] Awoyemi, D.O. and Kayode, S.J. (2005) A Maximal Order Collocation Method for Direct Solution of Initial Value Problems of General Second Order Ordinary Differential Equations. In: Proceedings of the Conference, the National Mathematical Center, Abuja.
[9] Kayode, S.J. (2005) An Improved Numerov Method for Direct Solution of General Second Order Initial Value Problems of Ordinary Differential Equations. In: Proceedings of the Seminar, the Mathematical Centre, Abuja.
[10] Kayode, S.J. (2008) An Efficient Zero Stable Numerical Method for Fourth Order Differential Equation. International Journal of Mathematics and Mathematical Sciences, 2008, Article ID: 364021. https://doi.org/10.1155/2008/364021
[11] Olabode, B.T. (2009) An Accurate Scheme by Block Method for Third Order Ordinary Differential Equations. Pacific Journal of Science and Technology, 10, 136-142.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.