Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2011 / Article

Research Article | Open Access

Volume 2011 |Article ID 810324 | 18 pages |

Solution of Higher-Order ODEs Using Backward Difference Method

Academic Editor: Francesco Pellicano
Received09 Nov 2010
Revised25 Mar 2011
Accepted13 May 2011
Published12 Jul 2011


The current numerical technique for solving a system of higher-order ordinary differential equations (ODEs) is to reduce it to a system of first-order equations then solving it using first-order ODE methods. Here, we propose a method to solve higher-order ODEs directly. The formulae will be derived in terms of backward difference in a constant stepsize formulation. The method developed will be validated by solving some higher-order ODEs directly with constant stepsize. To simplify the evaluations of the integration coefficients, we find the relationship between various orders. The result presented confirmed our hypothesis.

1. Introduction

Differential equations constantly arise in various branches of science and engineering. Many of these problems are in the form of higher-order ordinary differential equations (ODEs). A few examples where these problems can be found are, in the motion of projectiles, the bending of a thin clamped beam and population growth.

The popular practice for solving a system of higher-order ODEs is by reducing it to a system of first-order equations and then solving with first-order methods. These methods worked, so that methods for solving higher-order ODEs have been disregarded as robust codes. However, the work by Krogh [1], Suleiman [2], Majid and Suleiman [3], and Omar and Suleiman [4] has revived the interest in solving higher-order ODEs directly and the theoretical development of the methods.

Related works for solving higher-order ODEs can be found in Collatz [5], Gear [6], Krogh [1, 7], and Suleiman [2]. Krogh [7] proposed the direct integration (DI) method for nonstiff problems using modified divided difference while Suleiman [2] proposed the DI method using the standard divided difference. In this paper, we will derive the constant stepsize backward difference formulae of solving higher-order ODEs up to third order. The main reason for developing the constant stepsize formulae is that, in developing the theory on convergence and stability, the approach is through constant stepsize formulation. Another reason is that it is possible to use this formula in conjunction with other similar formulae as in Majid and Suleiman [3] to develop a code for variable stepsize and order.

The advantage of such a code is that the integration or differentiation constants are calculated only once at the start of the first step of integration, whereas other formulations calculate the constants at every step.

In this paper, we will focus only on nonstiff ODEs of the form š‘¦(š‘‘)ī‚€ī‚š‘Œī‚ī‚ī€·=š‘“š‘„,,(1.1)š‘Œ(š‘„)=š‘“š‘„,š‘¦,š‘¦ī…ž,š‘¦ī…žī…ž,ā€¦,š‘¦(š‘‘āˆ’1)ī€ø,ī€·Ģƒšœ‚=šœ‚,šœ‚ī…ž,šœ‚ī…žī…ž,ā€¦,šœ‚(š‘‘āˆ’1)ī€ø,(1.2) where ī‚š‘Œ(š‘Ž)=Ģƒšœ‚ in the interval š‘Žā‰¤š‘„ā‰¤š‘ and š‘‘ is the order of the ODE.

Without loss of generality, we will be considering the scalar equation in (1.1).

This paper will be organized as follows. First, the integration coefficients of the explicit constant stepsize backward difference formulation of the DI method will be derived. Then, the coefficients of the implicit method are formulated and their relationship with the explicit coefficients is shown. We start the derivation with the coefficients of the first-order system, which is given in Henrici [8]. Next, the second-order coefficients are derived and their relationship with the corresponding first-order coefficients is given, likewise the relationship of the coefficients for the second- and third-order systems. Finally, the method developed using backward difference will be validated numerically.

2. The Formulation of the Predict-Evaluate-Correct-Evaluate (PECE) Multistep Method in Its Backward Difference Form (MSBD) for Nonstiff Higher-Order ODEs

The code developed will be using the PECE mode. The predictor and corrector will have the following form:

predictor:š‘š‘Ÿš‘¦(š‘‘āˆ’š‘”)š‘›+1=š‘”āˆ’1ī“š‘–=0ā„Žš‘–š‘¦š‘–!š‘›(š‘‘āˆ’š‘”+1)+ā„Žš‘‘š‘˜āˆ’1ī“š‘–=0š›¾(š‘‘āˆ’š‘”),š‘–āˆ‡š‘–š‘“š‘›,š‘”=1,2,ā€¦,š‘‘,(2.1) corrector:š‘¦(š‘‘āˆ’š‘”)š‘›+1=š‘”āˆ’1ī“š‘–=0ā„Žš‘–š‘¦š‘–!š‘›(š‘‘āˆ’š‘”+1)+ā„Žš‘‘š‘˜āˆ’1ī“š‘–=0š›¾āˆ—(š‘‘āˆ’š‘”),š‘–āˆ‡š‘–š‘“š‘›+1,š‘”=1,2,ā€¦,š‘‘.(2.2) The corrector will be reformulated, so that it will be in terms of the predictor. The reformulated corrector can be written asš‘¦(š‘‘āˆ’š‘”)š‘›+1=š‘š‘Ÿš‘¦(š‘‘āˆ’š‘”)š‘›+1+š›¾(š‘‘āˆ’š‘”),š‘˜āˆ‡š‘˜š‘š‘Ÿī€·š‘“š‘›+1ī€ø,š‘”=1,2,ā€¦,š‘‘,(2.3) where š‘š‘Ÿš‘“š‘›+1 indicates š‘“š‘›+1 evaluated using predicted values. The integration coefficient š›¾(š‘‘āˆ’š‘”),š‘– and š›¾āˆ—(š‘‘āˆ’š‘”),š‘– will be derived using the method of generating function. Finally, the constant stepsize algorithm will be constructed and validated with some test problems and examples from physical situations.

3. Derivation up to Third-Order Explicit Integration Coefficients

Integrating (1.1) once yields š‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ī€œš‘„š‘›+1š‘„š‘›š‘“ī€·š‘„,š‘¦,š‘¦ī…ž,š‘¦ī…žī…žī€øš‘‘š‘„.(3.1) Let š‘ƒš‘›(š‘„) be the interpolating polynomial which interpolates the š‘˜ values (š‘„š‘›,š‘“š‘›),(š‘„š‘›āˆ’1,š‘“š‘›āˆ’1),ā€¦,(š‘„š‘›āˆ’š‘˜+1,š‘“š‘›āˆ’š‘˜+1), thenš‘ƒš‘›(š‘„)=š‘˜āˆ’1ī“š‘–=0(āˆ’1)š‘–īƒ©š‘–īƒŖāˆ‡āˆ’š‘ š‘–š‘“š‘›.(3.2) Next, approximating š‘“ in (3.1) with š‘ƒš‘›(š‘„) and lettingš‘„=š‘„š‘›+š‘ ā„Žorš‘ =š‘„āˆ’š‘„š‘›ā„Ž(3.3) gives usš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ī€œ10š‘˜āˆ’1ī“š‘–=0(āˆ’1)š‘–īƒ©š‘–īƒŖāˆ‡āˆ’š‘ š‘–š‘“š‘›š‘‘š‘ ,(3.4) orš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘˜āˆ’1ī“š‘–=0š›¾1,š‘–āˆ‡š‘–š‘“š‘›,(3.5) whereš›¾1,š‘–=(āˆ’1)š‘–ī€œ10īƒ©š‘–īƒŖāˆ’š‘ š‘‘š‘ .(3.6) The generating function šŗ1(š‘”) for the coefficients š›¾1,š‘– is defined as follows:šŗ1(š‘”)=āˆžī“š‘–=0š›¾1,š‘–š‘”š‘–.(3.7) Substituting š›¾1,š‘– in (3.6) in šŗ1(š‘”) givesšŗ1(š‘”)=āˆžī“š‘–=0(āˆ’š‘”)š‘–ī€œ10īƒ©š‘–īƒŖšŗāˆ’š‘ š‘‘š‘ ,1ī€œ(š‘”)=10(1āˆ’š‘”)(āˆ’š‘ )šŗš‘‘š‘ ,1ī€œ(š‘”)=10š‘’āˆ’š‘ log(1āˆ’š‘”)š‘‘š‘ (3.8) which leads tošŗ1ī‚ø(š‘”)=āˆ’(1āˆ’š‘”)āˆ’1āˆ’1log(1āˆ’š‘”)ī‚¹log(1āˆ’š‘”).(3.9) Equation (3.9) can be written asāˆ’īƒ©āˆžī“š‘–=0š›¾1,š‘–š‘”š‘–īƒŖī‚øš‘”log(1āˆ’š‘”)=ī‚¹(1āˆ’š‘”)(3.10)orī€·š›¾1,0+š›¾1,1š‘”+š›¾1,2š‘”2+š›¾1,3š‘”3ī€øī‚µš‘”+ā‹Æš‘”+22+š‘”33+š‘”44ī‚¶ī€·+ā‹Æ=š‘”1+š‘”+š‘”2+š‘”3ī€ø+ā‹Æ.(3.11) Hence, the coefficients of š›¾1,š‘˜ are given byš‘˜ī“š‘–=0ī‚€š›¾1,š‘–ī‚š›¾š‘˜āˆ’š‘–+1=1,1,š‘˜=1āˆ’š‘˜āˆ’1ī“š‘–=0ī‚€š›¾1,š‘–ī‚š‘˜āˆ’š‘–+1š‘˜=1,2,ā€¦,š›¾1,0=1.(3.12)

4. Second-Order ODE Formulae

Integrate (1.1) twice for second-order ODEs where š‘‘=2. Integrating once leads to the same coefficients as given in (3.6). Integrating twice yieldsš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ā„Ž2š‘˜ī“š‘–=0š›¾2,š‘–āˆ‡š‘–š‘“š‘›.(4.1) Substituting š‘„ with š‘  givesš›¾2,š‘–=(āˆ’1)š‘–ī€œ10(1āˆ’š‘ )īƒ©š‘–īƒŖ1!āˆ’š‘ š‘‘š‘ .(4.2) The generating function šŗ2(š‘”) of the coefficients š›¾2,š‘– is defined as followsšŗ2(š‘”)=āˆžī“š‘–=0š›¾2,š‘–š‘”š‘–.(4.3) Substituting (4.2) into šŗ2(š‘”) above givesšŗ2ī€œ(š‘”)=10(1āˆ’š‘ )š‘’1!āˆ’š‘ log(1āˆ’š‘”)š‘‘š‘ .(4.4) Substituting šŗ1(š‘”) into (4.4) yieldsšŗ21(š‘”)=ī‚ø11!āˆ’log(1āˆ’š‘”)1!šŗ1(š‘”)ī‚¹log(1āˆ’š‘”).(4.5) Equation (4.5) can be written asīƒ©āˆžī“š‘–=0š›¾2,š‘–š‘”š‘–īƒŖ1log(1āˆ’š‘”)=ī€ŗ1!1āˆ’1!šŗ1ī€»(š‘”)(4.6)orī€·š›¾2,0+š›¾2,1š‘”+š›¾2,2š‘”2ī€øī‚µš‘”+ā‹Æš‘”+22+š‘”33ī‚¶=1+ā‹Æī€ŗī€·š›¾1!āˆ’1+1!1,0+š›¾1,1š‘”+š›¾1,2š‘”2+ā‹Æī€øī€».(4.7) Hence the coefficients of š›¾2,š‘˜ in relation to coefficients of the previous order š›¾1,š‘˜ are given byš‘˜ī“š‘–=0š›¾2,š‘–š‘˜āˆ’š‘–+1=š›¾1,š‘˜+1,š›¾2,0=š›¾1,1,š›¾2,š‘˜=š›¾1,š‘˜+1āˆ’š‘˜ī“š‘–=0š›¾2,š‘–š‘˜āˆ’š‘–+1,š‘˜=1,2,ā€¦.(4.8)

5. Third-Order Formulae

Next, the case of the third-order ODE where š‘‘=3 will be considered. In the case of š‘¦ī…žī…ž,š‘¦ī…ž, the corresponding coefficients are š›¾1,š‘–,š›¾2,š‘– as in (3.6) and (4.2). For the solution š‘¦(š‘„), integrating three times yieldsš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ā„Ž2š‘¦2!ī…žī…žī€·š‘„š‘›ī€ø+ī€œš‘„š‘›+1š‘„š‘›ī€·š‘„š‘›+1ī€øāˆ’š‘„(2)š‘“ī€·(2)!š‘„,š‘¦,š‘¦ī…ž,š‘¦ī…žī…žī€øš‘‘š‘„(5.1) or in the backward difference formulation, given byš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ā„Ž2š‘¦2!ī…žī…žī€·š‘„š‘›ī€ø+ā„Ž3š‘˜ī“š‘–=0š›¾3,š‘–āˆ‡š‘–š‘“š‘›,(5.2) whereš›¾3,š‘–=(āˆ’1)š‘–ī€œ10(1āˆ’š‘ )2!2īƒ©š‘–īƒŖāˆ’š‘ š‘‘š‘ .(5.3) The generating function šŗ3(š‘”) of the coefficients š›¾3,š‘– is defined as follows:šŗ3(š‘”)=āˆžī“š‘–=0š›¾3,š‘–š‘”š‘–.(5.4)Substituting (5.3) into šŗ3(š‘”) above yieldsšŗ3ī€œ(š‘”)=10(1āˆ’š‘ )2!2š‘’āˆ’š‘ log(1āˆ’š‘”)š‘‘š‘ .(5.5) As in (4.4), we now substitute šŗ2(š‘”) in (5.5) which givesšŗ31(š‘”)=ī‚ø12!āˆ’log(1āˆ’š‘”)2!šŗ2(š‘”)ī‚¹log(1āˆ’š‘”).(5.6) Equation (5.6) can be written asīƒ©āˆžī“š‘–=0š›¾3,š‘–š‘”š‘–īƒŖ1log(1āˆ’š‘”)=ī€ŗ2!1āˆ’2!šŗ2ī€»(š‘”)(5.7)

orī€·š›¾3,0+š›¾3,1š‘”+š›¾3,2š‘”2ī€øī‚µš‘”+ā‹Æš‘”+22+š‘”33ī‚¶=1+ā‹Æī€ŗī€·š›¾2!āˆ’1+2!2,0+š›¾2,1š‘”+š›¾2,2š‘”2+ā‹Æī€øī€».(5.8) Hence, the coefficients of š›¾3,š‘˜ in relation to coefficients of the previous order š›¾2,š‘˜ are given byš‘˜ī“š‘–=0š›¾3,š‘–š‘˜āˆ’š‘–+1=š›¾2,š‘˜+1,š›¾3,0=š›¾2,1,š›¾3,š‘˜=š›¾2,š‘˜+1āˆ’š‘˜ī“š‘–=0š›¾3,š‘–š‘˜āˆ’š‘–+1,š‘˜=1,2,ā€¦.(5.9)

6. Derivation up to the Third-Order Implicit Integration Coefficients

Integrating (1.1) once yields š‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ī€œš‘„š‘›+1š‘„š‘›š‘“ī€·š‘„,š‘¦,š‘¦ī…ž,š‘¦ī…žī…žī€øš‘‘š‘„.(6.1) Let š‘ƒš‘›(š‘„) be the interpolating polynomial which interpolates the š‘˜ values (š‘„š‘›,š‘“š‘›),(š‘„š‘›āˆ’1,š‘“š‘›āˆ’1),ā€¦,(š‘„š‘›āˆ’š‘˜+1,š‘“š‘›āˆ’š‘˜+1):š‘ƒš‘›(š‘„)=š‘˜ī“š‘–=0(āˆ’1)š‘–īƒ©š‘–īƒŖāˆ‡āˆ’š‘ š‘–š‘“š‘›+1.(6.2) As in the previous derivation, we chooseš‘„=š‘„š‘›+1+š‘ ā„Ž,orš‘ =š‘„āˆ’š‘„š‘›+1ā„Ž.(6.3) Replacing š‘„by š‘  yields š‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ī€œ0š‘˜āˆ’1ī“š‘–=0(āˆ’1)š‘–īƒ©š‘–īƒŖāˆ‡āˆ’š‘ š‘–š‘“š‘›+1š‘‘š‘ .(6.4) Simplify š‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘˜ī“š‘–=0š›¾āˆ—1,š‘–āˆ‡š‘–š‘“š‘›,(6.5) whereš›¾āˆ—1,š‘–=(āˆ’1)š‘–ī€œ0āˆ’1īƒ©š‘–īƒŖāˆ’š‘ š‘‘š‘ .(6.6) The generating function šŗāˆ—1(š‘”) of the coefficientsš›¾āˆ—1,š‘– is defined as follows:šŗāˆ—1(š‘”)=āˆžī“š‘–=0š›¾āˆ—1,š‘–š‘”š‘–(6.7)

oršŗāˆ—1(š‘”)=āˆžī“š‘–=0(āˆ’š‘”)š‘–ī€œ0āˆ’1īƒ©š‘–īƒŖšŗāˆ’š‘ š‘‘š‘ ,āˆ—1ī€œ(š‘”)=0āˆ’1(1āˆ’š‘”)(āˆ’š‘ )šŗš‘‘š‘ ,āˆ—1ī€œ(š‘”)=0āˆ’1š‘’āˆ’š‘ log(1āˆ’š‘”)š‘‘š‘ (6.8) which leads tošŗāˆ—1ī‚ø1(š‘”)=āˆ’āˆ’log(1āˆ’š‘”)(1āˆ’š‘”)ī‚¹log(1āˆ’š‘”).(6.9) For the case š‘‘=2, the approximate solution of š‘¦ has the formš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ī€œš‘„š‘›+1š‘„š‘›ī€·š‘„š‘›+1ī€øāˆ’š‘„(1)!(1)š‘“ī€·š‘„,š‘¦,š‘¦ī…ž,š‘¦ī…žī…žī€øš‘‘š‘„.(6.10) The coefficients are given byš›¾āˆ—2,š‘–=(āˆ’1)š‘–ī€œ0āˆ’1(āˆ’š‘ )īƒ©š‘–īƒŖ1!āˆ’š‘ š‘‘š‘ ,(6.11) where š›¾āˆ—2,š‘– are the coefficients of the backward difference formulation of (6.11) which can be represented byš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ā„Ž2š‘˜ī“š‘–=0š›¾āˆ—2,š‘–āˆ‡š‘–š‘“š‘›.(6.12) The generating function šŗāˆ—2(š‘”) for the coefficients š›¾āˆ—2,š‘– is defined as follows:šŗāˆ—2(š‘”)=āˆžī“š‘–=0š›¾āˆ—2,š‘–š‘”š‘–.(6.13) Substituting (6.11) into šŗāˆ—2(š‘”) above,šŗāˆ—2ī€œ(š‘”)=0āˆ’1(āˆ’š‘ )š‘’1!āˆ’š‘ log(1āˆ’š‘”)š‘‘š‘ .(6.14) Solving (6.14) with the substitution of (3.8) produces the relationshipšŗāˆ—21(š‘”)=ī‚ø1!(1āˆ’š‘”)āˆ’log(1āˆ’š‘”)1!šŗāˆ—1(š‘”)ī‚¹log(1āˆ’š‘”).(6.15) Integrating (1.1) thrice yieldsš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ā„Ž2š‘¦2!ī…žī…žī€·š‘„š‘›ī€ø+ī€œš‘„š‘›+1š‘„š‘›ī€·š‘„š‘›+1ī€øāˆ’š‘„(2)š‘“ī€·(2)!š‘„,š‘¦,š‘¦ī…ž,š‘¦ī…žī…žī€øš‘‘š‘„.(6.16) The coefficients are given byš›¾āˆ—3,š‘–=(āˆ’1)š‘–ī€œ0āˆ’1(āˆ’š‘ )2!2īƒ©š‘–īƒŖāˆ’š‘ š‘‘š‘ ,(6.17) where š›¾āˆ—3,š‘– are the coefficients of the backward difference formulation of (6.17) which can be represented byš‘¦ī€·š‘„š‘›+1ī€øī€·š‘„=š‘¦š‘›ī€ø+ā„Žš‘¦ī…žī€·š‘„š‘›ī€ø+ā„Ž2š‘¦2!ī…žī…žī€·š‘„š‘›ī€ø+ā„Ž3š‘˜ī“š‘–=0š›¾āˆ—3,š‘–āˆ‡š‘–š‘“š‘›.(6.18)

The generating function šŗāˆ—3(š‘”) of the coefficients š›¾āˆ—3,š‘– is defined as follows:šŗāˆ—3(š‘”)=āˆžī“š‘–=0š›¾āˆ—3,š‘–š‘”š‘–.(6.19)

Substituting (6.17) into šŗāˆ—3(š‘”) above yieldsšŗāˆ—3ī€œ(š‘”)=0āˆ’1(āˆ’š‘ )2!2š‘’āˆ’š‘ log(1āˆ’š‘”)š‘‘š‘ .(6.20) Solving (6.11) with the substitution of (6.20) produces the relationshipšŗāˆ—31(š‘”)=ī‚ø2!(1āˆ’š‘”)āˆ’log(1āˆ’š‘”)2!šŗāˆ—2(š‘”)ī‚¹log(1āˆ’š‘”).(6.21)

7. The Relationship between the Explicit and Implicit Integration Coefficients

Calculating the integration coefficients directly is time consuming when large numbers of integration are involved. A more efficient way of computing the coefficients is by obtaining a recursive relationship between the coefficients. With this recursive relationship, we are able to obtain the implicit integration coefficient with minimal time consumption. The relationship between the explicit and implicit coefficients is expressed below.

For first-order coefficients,šŗāˆ—1ī‚ø1(š‘”)=āˆ’āˆ’log(1āˆ’š‘”)1āˆ’š‘”ī‚¹log(1āˆ’š‘”).(7.1)

It can be written asšŗāˆ—1ī‚ø1(š‘”)=āˆ’(1āˆ’š‘”)āˆ’1(1āˆ’š‘”)log(1āˆ’š‘”)ī‚¹log(1āˆ’š‘”).(7.2) By substituting šŗ11(š‘”)=āˆ’1(1āˆ’š‘”)log(1āˆ’š‘”)log(1āˆ’š‘”)(7.3) into (7.2), we have šŗāˆ—1(š‘”)=(1āˆ’š‘”)šŗ1īƒ©(š‘”),āˆžī“š‘–=0š›¾āˆ—1,š‘–š‘”š‘–īƒŖ=īƒ©(1āˆ’š‘”)āˆžī“š‘–=0š›¾1,š‘–š‘”š‘–īƒŖ.(7.4) Expanding the equation yieldsī‚€š›¾āˆ—1,0+š›¾āˆ—1,1š‘”+š›¾āˆ—1,2š‘”2ī‚=1+ā‹Æī€·1+š‘”+š‘”2ī€øī€·š›¾+ā‹Æ1,0+š›¾1,1š‘”+š›¾1,2š‘”2ī€ø,ī‚€š›¾+ā‹Æāˆ—1,0+š›¾āˆ—1,1š‘”+š›¾āˆ—1,2š‘”2ī‚ī€·+ā‹Æ1+š‘”+š‘”2ī€ø=ī€·š›¾+ā‹Æ1,0+š›¾1,1š‘”+š›¾1,2š‘”2ī€ø.+ā‹Æ(7.5) This gives the recursive relationshipš‘˜ī“š‘–=0š›¾āˆ—1,š‘–=š›¾1,š‘˜.(7.6) For second-order coefficient, šŗāˆ—21(š‘”)=āˆ’ī‚ø11!āˆ’log(1āˆ’š‘”)1!šŗāˆ—1(š‘”)ī‚¹log(1āˆ’š‘”).(7.7) It can be written asšŗāˆ—2(š‘”)=(1āˆ’š‘”)ī‚ø11!āˆ’log(1āˆ’š‘”)1!šŗāˆ—1(š‘”)ī‚¹(1āˆ’š‘”)log(1āˆ’š‘”).(7.8) Substituting (7.4) into the equation above givesšŗāˆ—2(š‘”)=(1āˆ’š‘”)ī‚ø11!āˆ’log(1āˆ’š‘”)1!(1āˆ’š‘”)šŗ1(š‘”)ī‚¹(1āˆ’š‘”)log(1āˆ’š‘”)(7.9)oršŗāˆ—2(š‘”)=(1āˆ’š‘”)ī‚ø11!āˆ’log(1āˆ’š‘”)1!šŗ1(š‘”)ī‚¹log(1āˆ’š‘”).(7.10) Substituting (4.5) into (7.10) givesšŗāˆ—2(š‘”)=(1āˆ’š‘”)šŗ2īƒ©(š‘”),āˆžī“š‘–=0š›¾āˆ—2,š‘–š‘”š‘–īƒŖ=īƒ©(1āˆ’š‘”)āˆžī“š‘–=0š›¾2,š‘–š‘”š‘–īƒŖ.(7.11) Expanding the equation, we haveī‚€š›¾āˆ—2,0+š›¾āˆ—2,1š‘”+š›¾āˆ—2,2š‘”2ī‚=1+ā‹Æī€·1+š‘”+š‘”2ī€øī€·š›¾+ā‹Æ2,0+š›¾2,1š‘”+š›¾2,2š‘”2ī€ø,ī‚€š›¾+ā‹Æāˆ—2,0+š›¾āˆ—2,1š‘”+š›¾āˆ—2,2š‘”2ī‚ī€·+ā‹Æ1+š‘”+š‘”2ī€ø=ī€·š›¾+ā‹Æ2,0+š›¾2,1š‘”+š›¾2,2š‘”2ī€ø.+ā‹Æ(7.12) The above gives the relationshipš‘˜ī“š‘–=0š›¾āˆ—2,š‘–=š›¾2,š‘˜.(7.13) For third-order coefficient, we havešŗāˆ—31(š‘”)=āˆ’ī‚ø12!āˆ’log(1āˆ’š‘”)2!šŗāˆ—2(š‘”)ī‚¹log(1āˆ’š‘”).(7.14) It can be written asšŗāˆ—3(š‘”)=(1āˆ’š‘”)ī‚ø12!āˆ’log(1āˆ’š‘”)2!šŗāˆ—2(š‘”)ī‚¹(1āˆ’š‘”)log(1āˆ’š‘”).(7.15) Substituting (7.10) into (7.15) givesšŗāˆ—3(š‘”)=(1āˆ’š‘”)ī‚ø12!āˆ’log(1āˆ’š‘”)2!(1āˆ’š‘”)šŗ2(š‘”)ī‚¹(1āˆ’š‘”)log(1āˆ’š‘”)(7.16)

oršŗāˆ—3(š‘”)=(1āˆ’š‘”)ī‚ø12!āˆ’log(1āˆ’š‘”)2!šŗ2(š‘”)ī‚¹log(1āˆ’š‘”).(7.17) Substituting šŗ31(š‘”)=ī‚ø12!āˆ’log(1āˆ’š‘”)2!šŗ2(š‘”)ī‚¹log(1āˆ’š‘”)(7.18)into (6.15) leads tošŗāˆ—3(š‘”)=(1āˆ’š‘”)šŗ3īƒ©(š‘”),āˆžš‘˜ī“š‘–=0š›¾āˆ—3,š‘–š‘”š‘–īƒŖīƒ©=(1āˆ’š‘”)āˆžī“š‘–=0š›¾3,š‘–š‘”š‘–īƒŖ.(7.19) Expanding the equation into,ī‚€š›¾āˆ—3,0+š›¾āˆ—3,1š‘”+š›¾āˆ—3,2š‘”2ī‚=1+ā‹Æī€·1+š‘”+š‘”2ī€øī€·š›¾+ā€¦3,0+š›¾3,1š‘”+š›¾3,2š‘”2ī€ø,ī‚€š›¾+ā‹Æāˆ—3,0+š›¾āˆ—3,1š‘”+š›¾āˆ—3,2š‘”2ī‚ī€·+ā‹Æ1+š‘”+š‘”2ī€ø=ī€·š›¾+ā€¦3,0+š›¾3,1š‘”+š›¾3,2š‘”2ī€ø+ā‹Æ(7.20)

which leads to a recursive relationship š‘˜ī“š‘–=0š›¾āˆ—3,š‘–=š›¾3,š‘˜.(7.21)

Tables 1 and 2 are a few examples of the explicit and implicit integration coefficients.

š‘˜ 0123456

š›¾ 1 , š‘˜ 11/25/123/8251/72095/28819087/60480
š›¾ 2 , š‘˜ 1/21/61/819/1803/32863/10080275/3456
š›¾ 3 , š‘˜ 1/61/247/24017/72041/2016731/403208563/518400

š‘˜ 0123456

š›¾ āˆ— 1 , š‘˜ 1āˆ’1/2āˆ’1/12āˆ’1/24āˆ’19/720āˆ’3/160āˆ’813/60480
š›¾ āˆ— 2 , š‘˜ 1/2āˆ’1/3āˆ’1/24āˆ’7/360āˆ’17/1440āˆ’41/5040āˆ’731/120960
š›¾ āˆ— 3 , š‘˜ 1/6āˆ’1/8āˆ’1/80āˆ’1/180āˆ’11/3360āˆ’89/40320āˆ’5849/3628800

8. Numerical Results

For error calculations, we will be using the three error tests, namely, absolute error, relative error, and mixed error tests. The error formula is given by,Error=Maxš‘„š‘›ā€–ā€–ā€–ā€–š‘¦ī€·š‘„š‘›ī€øāˆ’š‘¦š‘›ī€·š‘„š“+šµš‘¦š‘›ī€øā€–ā€–ā€–ā€–,(8.1) where š“=1, šµ=0 gives the absolute error test, š“=0, šµ=1 gives the relative error test, and š“=1, šµ=1 gives the mixed error test.

In (8.1), š‘¦(š‘„š‘›) is the exact solution for the problem considered and š‘¦š‘› the computed solution. In a general code when the exact solution is not available for the relative error, š‘¦(š‘„š‘›) is replaced by š‘¦š‘› the computed value.

When ā€–š‘¦(š‘„š‘›)āˆ’š‘¦š‘›ā€– is small, the error in (8.1) will approximate the absolute error. However, when it is large, the mixed error test will approximate the relative error. The numerical results give the three errors.

The following notations hold MAX ABS: maximum error when using absolute error test, MAX MIX: maximum error when using mixed error test, MAX REL: maximum error when using relative error test, ā„Ž: step size selected.

For the choice of problems to be tested, we choose four linear problems consisting of a second- and a third-order problem. The third problem is a mix system of second- and first-order equations and the fourth problem is a system of three second-order equations. Our reason for choosing the linear problems is that if the formulae are correct, then they should solve linear problems. The choice of system of equations is to raise the degree of difficulty of solving the problems. The rest of the problems are nonlinear, which occur in physical situations. The choices of the physical problems are those with exact solutions known. We give our comments on the numerical results right after the numerical Tables 3, 4, 5, 6, 7, 8, and 9.


1 0 āˆ’ 1 4.73734 + 0261.18857 āˆ’ 0031.18857 ā€“ 003
1 0 āˆ’ 2 4.70433 + 0231.16697 āˆ’ 0061.16697 ā€“ 006
1 0 āˆ’ 3 4.42360 + 0201.18335 āˆ’ 0091.18335 āˆ’ 009
1 0 āˆ’ 4 4.05706 + 0201.01668 āˆ’ 0091.01668 āˆ’ 009
1 0 āˆ’ 5 3.69529 + 0219.26023 āˆ’ 0099.26024 āˆ’ 009


1 0 āˆ’ 1 1.60072 + 0231.40364 āˆ’ 0031.40364 āˆ’ 003
1 0 āˆ’ 2 1.52595 + 0201.33620 āˆ’ 0061.33620 āˆ’ 006
1 0 āˆ’ 3 1.55425 + 0171.36098 āˆ’ 0091.36098 āˆ’ 009
1 0 āˆ’ 4 1.27710 + 0161.11807 āˆ’ 0101.11807 āˆ’ 010
1 0 āˆ’ 5 2.20767 + 0171.93311 āˆ’ 0091.93311 āˆ’ 009


1 0 āˆ’ 1 1.74071 āˆ’ 0021.68285 āˆ’ 0036.12179 āˆ’ 003
1 0 āˆ’ 2 2.32581 āˆ’ 0054.04838 āˆ’ 0065.34483 āˆ’ 005
1 0 āˆ’ 3 2.38274 āˆ’ 0084.39132 āˆ’ 0095.74568 āˆ’ 007
1 0 āˆ’ 4 5.33390 āˆ’ 0095.47762 āˆ’ 0109.42528 āˆ’ 009
1 0 āˆ’ 5 4.54132 āˆ’ 0084.54132 āˆ’ 0085.25218 āˆ’ 007


1 0 āˆ’ 1 1.40347 + 0011.55464 + 0006.27279 + 002
1 0 āˆ’ 2 1.78122 āˆ’ 0021.75004 āˆ’ 0021.60254 + 001
1 0 āˆ’ 3 2.32649 āˆ’ 0052.32643 āˆ’ 0059.97054 āˆ’ 001
1 0 āˆ’ 4 2.44607 āˆ’ 0082.44600 āˆ’ 0089.82644 āˆ’ 001
1 0 āˆ’ 5 9.90539 āˆ’ 0109.90531 āˆ’ 0109.99724 āˆ’ 001


1 0 āˆ’ 1 2.34537 āˆ’ 0014.96059 āˆ’ 0041.00000 + 000
1 0 āˆ’ 2 2.85286 āˆ’ 0042.86491 āˆ’ 0071.00000 + 000
1 0 āˆ’ 3 2.91084 āˆ’ 0076.19617 āˆ’ 0101.00083 + 000
1 0 āˆ’ 4 3.00856 āˆ’ 0087.74497 āˆ’ 0131.80152 + 000
1 0 āˆ’ 5 2.43305 āˆ’ 0071.09024 āˆ’ 0117.08729 + 000


1 0 āˆ’ 1 1.01336 āˆ’ 0048.29238 āˆ’ 0055.15346 āˆ’ 002
1 0 āˆ’ 2 8.33352 āˆ’ 0088.33254 āˆ’ 0089.08294 āˆ’ 004
1 0 āˆ’ 3 9.71863 āˆ’ 0119.71317 āˆ’ 0114.58972 āˆ’ 006
1 0 āˆ’ 4 5.45085 āˆ’ 0105.45075 āˆ’ 0103.45420 āˆ’ 004
1 0 āˆ’ 5 4.75544 āˆ’ 0094.75543 āˆ’ 0091.91598 āˆ’ 002


1 0 āˆ’ 1 6.95439 āˆ’ 0025.63803 āˆ’ 0022.97860 āˆ’ 001
1 0 āˆ’ 2 8.31669 āˆ’ 0057.12977 āˆ’ 0054.99583 āˆ’ 004
1 0 āˆ’ 3 8.60027 āˆ’ 0087.37166 āˆ’ 0085.16016 āˆ’ 007
1 0 āˆ’ 4 8.66052 āˆ’ 0117.42330 āˆ’ 0115.19631 āˆ’ 010
1 0 āˆ’ 5 1.86645 āˆ’ 0121.55538 āˆ’ 0129.33222 āˆ’ 012

Problem 1. š‘¦ī…žī…ž=2š‘¦ī…žāˆ’š‘¦,š‘¦(0)=0,š‘¦ī…ž(0)=1,0ā‰¤š‘„ā‰¤64.(8.2)

Solution 1. š‘¦=š‘„š‘’š‘„.(8.3)

Source: Krogh [7].

This is a linear equation used by Krogh [7] to test his code. The solution increases exponentially to a maximum value of š‘¦(64)=64š‘’64 which is considered large and therefore not suitable for absolute error test and hence the large values of the error.

Problem 2. š‘¦ī…žī…žī…ž=2š‘¦ī…žī…žāˆ’4,0ā‰¤š‘„ā‰¤30,š‘¦(0)=1,š‘¦ī…ž(0)=2,š‘¦ī…žī…ž(0)=6.(8.4)

Solution 2. š‘¦(š‘„)=š‘„2+š‘’2š‘„.(8.5)

Source: Omar and Suleiman [4].

This is a third-order problem with an exponential solution. The difference between Problems 1 and 2 is that one is third order and the other is second order. Again, absolute error test does not work for the same reason given above.

Problem 3. š‘¦1ī…žī…ž=āˆ’2š‘¦ī…ž1āˆ’5š‘¦2+3,š‘¦ī…ž2=š‘¦ī…ž1+2š‘¦2,š‘¦0ā‰¤š‘„ā‰¤16šœ‹,1(0)=0,š‘¦ī…ž1(0)=0,š‘¦2(0)=1.(8.6)

Solution 3. š‘¦1š‘¦(š‘„)=2cosš‘„+6sinš‘„āˆ’2āˆ’6š‘„,2(š‘„)=āˆ’2cosš‘„+2sinš‘„+3.(8.7)

Source: Bronson [9].

For this problem, all error tests worked well.

Problem 4. š‘¦1ī…žī…ž=1āˆ’š‘¦2āˆ’š‘¦3,š‘¦2ī…žī…ž=š‘¦3āˆ’š‘¦1,š‘¦3ī…žī…ž=š‘¦ī…ž1+š‘¦ī…ž2,š‘¦0ā‰¤š‘„ā‰¤4šœ‹,1(0)=0,š‘¦ī…ž1(0)=1,š‘¦2(0)=0,š‘¦ī…ž2(0)=0,š‘¦3(0)=āˆ’1,š‘¦ī…ž3(0)=1.(8.8)

Solution 4. š‘¦1š‘¦(š‘„)=sinš‘„,2š‘¦(š‘„)=cosš‘„āˆ’1,2(š‘„)=sinš‘„āˆ’cosš‘„.(8.9)

Source: Bronson [9].

This problem does not work for relative error test because of the small value of the solution for certain values of š‘„.

Problem 5. š‘¦ī…žī…žī…ž1=āˆ’š‘„š‘¦ī…žī…ž+1š‘„2š‘¦ī…ž+1š‘„,1ā‰¤š‘„ā‰¤50,š‘¦(1)=2621ln2(2)+99104,š‘¦ī…ž(1)=āˆ’40521ln(2)āˆ’13,š‘¦ī…žī…ž3(1)=+4267ln(2).(8.10)

Solution 5. š‘„š‘¦(š‘„)=28ī‚€ī‚€š‘„2ln2ī‚āˆ’33āˆ’2133ī‚+ī‚€1ln(2)3āˆ’26ī‚€š‘„21ln2ī‚ī‚ln(2)+3326.(8.11)

Source: Russel and Shampine [10].

This problem is the symmetrical bending of a laterally loaded circular plate.

The numerical results of this problem show the failure to control the error using relative error test. This is because the solution is zero when š‘„=2.

Problem 6. š‘¦1ī…žī…žš‘¦=āˆ’1š‘Ÿ3,š‘¦1(0)=1,š‘¦ī…ž1(š‘¦0)=0,2ī…žī…žš‘¦=āˆ’2š‘Ÿ3,š‘¦2(0)=0,š‘¦ī…ž2ī€·š‘¦(0)=1,š‘Ÿ=21+š‘¦22ī€ø1/2,0ā‰¤š‘„ā‰¤16šœ‹.(8.12)

Solution 6. š‘¦1(š‘„)=cosš‘„,š‘¦2(š‘„)=sinš‘„.(8.13)

Source: Shampine and Gordon [11].

This problem is Newtonā€™s equations of motion for the two-body problem.

Again, relative error test does not work too well for this problem because š‘¦š‘› is very small at certain points š‘„š‘›.

Problem 7. š‘¦ī…žī…ž=2š‘¦3,0ā‰¤š‘„ā‰¤5,š‘¦(0)=1,š‘¦ī…ž(0)=āˆ’1.(8.14)

Solution 7. 1š‘¦(š‘„)=š‘„+1.(8.15)

Source: Robert Jr. [12].

For this problem, all error tests worked well.

All the numerical results show that the errors in the mixed error mode give a reliable error estimate for all the problems given. The absolute error mode failed to give meaningful error results for Problems 1 and 2. This is because the value of ā€–š‘¦š‘›ā€– increases as š‘„ increases and this becomes large. Similarly, for Problems 4, 5, and 6, the relative error failed to give an acceptable result because ā€–š‘¦š‘›ā€– is small.

The research work done shows that the method developed for solving higher-order ODEs directly using the backward difference is successful. We recommend that, for multistep method, the error control procedure should use the mixed error test. This research suggests the potential of this work developing a robust code for solving higher-order ODEs directly.


This paper has been supported by the UPM Graduate Research Fellowship (GRF) and Science Fund.


  1. F. T. Krogh, ā€œAlgorithms for changing the step size,ā€ SIAM Journal on Numerical Analysis, vol. 10, pp. 949ā€“965, 1973. View at: Publisher Site | Google Scholar
  2. M. B. Suleiman, ā€œSolving nonstiff higher order ODEs directly by the direct integration method,ā€ Applied Mathematics and Computation, vol. 33, no. 3, pp. 197ā€“219, 1989. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  3. Z. A. Majid and M. B. Suleiman, ā€œDirect integration implicit variable steps method for solving higher order systems of ordinary differential equations directly,ā€ Sains Malaysiana, vol. 35, no. 2, pp. 63ā€“68, 2006. View at: Google Scholar | Zentralblatt MATH
  4. Z. B. Omar and M. B. Suleiman, ā€œParallel 2-point explicit block method for solving higher order ordinary differential equations directly,ā€ International Journal of Simulation and Process Modelling, vol. 2, no. 3, pp. 227ā€“231, 2006. View at: Publisher Site | Google Scholar
  5. L. Collatz, The Numerical Treatment of Differential Equations, Springler, New York, NY, USA, 1966. View at: Zentralblatt MATH
  6. C. W. Gear, ā€œThe numerical integration of ordinary differential equations,ā€ Mathematics of Computation, vol. 21, pp. 146ā€“156, 1967. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  7. F. T. Krogh, ā€œA variable step, variable order multistep method for the numerical solution of ordinary differential equations,ā€ in Proceedings of the IFIP Congress in Information Processing, vol. 68, pp. 194ā€“199, 1968. View at: Google Scholar
  8. P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, New York, NY, USA, 1962.
  9. R. Bronson, Modern Introductory Differential Equation: Schaumā€™s Outline Series, McGraw-Hill, USA, 1973. View at: Zentralblatt MATH
  10. R. D. Russell and L. F. Shampine, ā€œA collocation method for boundary value problems,ā€ Numerische Mathematik, vol. 19, pp. 1ā€“28, 1972. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  11. L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary Differential Equations, W. H. Freeman and Co., San Francisco, Calif, USA, 1975.
  12. C. E. Roberts Jr., Ordinary Differential Equations. A Computational Approach, Prentice-Hall, Englewood Cliffs, NJ, USA, 1979.

Copyright Ā© 2011 Mohamed Bin Suleiman et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

1835Ā Views | 711Ā Downloads | 6Ā Citations
 PDF  Download Citation  Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly and safely as possible. Any author submitting a COVID-19 paper should notify us at to ensure their research is fast-tracked and made available on a preprint server as soon as possible. We will be providing unlimited waivers of publication charges for accepted articles related to COVID-19. Sign up here as a reviewer to help fast-track new submissions.