Solution of Higher-Order ODEs Using Backward Difference Method

Abstract

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)^𝑖¹⁰−𝑖𝑠𝑑𝑠.(3.6) The generating function 𝐺₁(𝑡) for the coefficients 𝛾_1,𝑖 is defined as follows:𝐺₁(𝑡)=^∞_𝑖=0𝛾_1,𝑖𝑡^𝑖.(3.7) Substituting 𝛾_1,𝑖 in (3.6) in 𝐺₁(𝑡) gives𝐺₁(𝑡)=^∞_𝑖=0(−𝑡)^𝑖𝐺¹⁰−𝑖𝑠𝑑𝑠,₁𝐺(𝑡)=¹⁰(1−𝑡)^(−𝑠)𝑑𝑠,₁(𝑡)=¹⁰𝑒^{−𝑠log(1−𝑡)}𝑑𝑠(3.8) which leads to𝐺₁(𝑡)=−(1−𝑡)⁻¹log(1−𝑡)−log(11−𝑡).(3.9) Equation (3.9) can be written as−^∞_𝑖=0𝛾_1,𝑖𝑡^𝑖log(1−𝑡)=(1−𝑡𝑡)(3.10)or𝛾_1,0+𝛾_1,1𝑡+𝛾_1,2𝑡²+𝛾_1,3𝑡³+⋯𝑡+𝑡2²+𝑡3³+𝑡4⁴+⋯=𝑡1+𝑡+𝑡²+𝑡³+⋯.(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)^𝑖¹⁰(1−𝑠)1!−𝑖𝑠𝑑𝑠.(4.2) The generating function 𝐺₂(𝑡) of the coefficients 𝛾_2,𝑖 is defined as follows𝐺₂(𝑡)=^∞_𝑖=0𝛾_2,𝑖𝑡^𝑖.(4.3) Substituting (4.2) into 𝐺₂(𝑡) above gives𝐺₂(𝑡)=¹⁰(1−𝑠)1!𝑒^{−𝑠log(1−𝑡)}𝑑𝑠.(4.4) Substituting 𝐺₁(𝑡) into (4.4) yields𝐺₂(𝑡)=11!log(11−𝑡)−1!𝐺¹(𝑡)log(1−𝑡).(4.5) Equation (4.5) can be written as^∞_𝑖=0𝛾_2,𝑖𝑡^𝑖log(1−𝑡)=11!1−1!𝐺¹(𝑡)(4.6)or𝛾_2,0+𝛾_2,1𝑡+𝛾_2,2𝑡²+⋯𝑡+𝑡2²+𝑡3³+⋯=11!−1+1!𝛾^1,0+𝛾_1,1𝑡+𝛾_1,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!𝑦^𝑥_𝑛+^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑥⁽²⁾(2)!𝑓𝑥,𝑦,𝑦^,𝑦^𝑑𝑥(5.1) or in the backward difference formulation, given by𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^𝑥_𝑛+ℎ²2!𝑦^𝑥_𝑛+ℎ^3𝑘_𝑖=0𝛾_3,𝑖∇^𝑖𝑓_𝑛,(5.2) where𝛾_3,𝑖=(−1)^𝑖¹⁰(1−𝑠)2!²−𝑖𝑠𝑑𝑠.(5.3) The generating function 𝐺₃(𝑡) of the coefficients 𝛾_3,𝑖 is defined as follows:𝐺₃(𝑡)=^∞_𝑖=0𝛾_3,𝑖𝑡^𝑖.(5.4)Substituting (5.3) into 𝐺₃(𝑡) above yields𝐺₃(𝑡)=¹⁰(1−𝑠)2!²𝑒^{−𝑠log(1−𝑡)}𝑑𝑠.(5.5) As in (4.4), we now substitute 𝐺₂(𝑡) in (5.5) which gives𝐺₃(𝑡)=12!log(11−𝑡)−2!𝐺²(𝑡)log(1−𝑡).(5.6) Equation (5.6) can be written as^∞_𝑖=0𝛾_3,𝑖𝑡^𝑖log(1−𝑡)=21!1−2!𝐺²(𝑡)(5.7)

or𝛾_3,0+𝛾_3,1𝑡+𝛾_3,2𝑡²+⋯𝑡+𝑡2²+𝑡3³+⋯=21!−1+2!𝛾^2,0+𝛾_2,1𝑡+𝛾_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)^𝑖⁰⁻¹−𝑖𝑠𝑑𝑠.(6.6) The generating function 𝐺^∗1(𝑡) of the coefficients𝛾^∗1,𝑖 is defined as follows:𝐺^∗1(𝑡)=^∞_𝑖=0𝛾^∗1,𝑖𝑡^𝑖(6.7)

or𝐺^∗1(𝑡)=^∞_𝑖=0(−𝑡)^𝑖𝐺⁰⁻¹−𝑖𝑠𝑑𝑠,_∗1𝐺(𝑡)=⁰⁻¹(1−𝑡)^(−𝑠)𝑑𝑠,_∗1(𝑡)=⁰⁻¹𝑒^{−𝑠log(1−𝑡)}𝑑𝑠(6.8) which leads to𝐺^∗1(𝑡)=−1log(1−𝑡)−log(1(−1−𝑡)𝑡).(6.9) For the case 𝑑=2, the approximate solution of 𝑦 has the form𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^𝑥_𝑛+^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑥(1)!⁽¹⁾𝑓𝑥,𝑦,𝑦^,𝑦^𝑑𝑥.(6.10) The coefficients are given by𝛾^∗2,𝑖=(−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(𝑡)=⁰⁻¹(−𝑠)1!𝑒^{−𝑠log(1−𝑡)}𝑑𝑠.(6.14) Solving (6.14) with the substitution of (3.8) produces the relationship𝐺^∗2(𝑡)=11!log(1(−1−𝑡)𝑡)−1!𝐺^∗1(𝑡)log(1−𝑡).(6.15) Integrating (1.1) thrice yields𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^𝑥_𝑛+ℎ²2!𝑦^𝑥_𝑛+^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑥⁽²⁾(2)!𝑓𝑥,𝑦,𝑦^,𝑦^𝑑𝑥.(6.16) The coefficients are given by𝛾^∗3,𝑖=(−1)^𝑖⁰⁻¹(−𝑠)2!²−𝑖𝑠𝑑𝑠,(6.17) where 𝛾^∗3,𝑖 are the coefficients of the backward difference formulation of (6.17) which can be represented by𝑦𝑥_𝑛+1=𝑦𝑥_𝑛+ℎ𝑦^𝑥_𝑛+ℎ²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(𝑡)=⁰⁻¹(−𝑠)2!²𝑒^{−𝑠log(1−𝑡)}𝑑𝑠.(6.20) Solving (6.11) with the substitution of (6.20) produces the relationship𝐺^∗3(𝑡)=12!log(1(−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(𝑡)=−1log(1−𝑡)−log1(−1−𝑡𝑡).(7.1)

It can be written as𝐺^∗1(𝑡)=−(1−𝑡)1(1−𝑡)log(1−𝑡)−log(11−𝑡).(7.2) By substituting 𝐺₁(𝑡)=1(1−𝑡)log(1−𝑡)−log(11−𝑡)(7.3) into (7.2), we have 𝐺^∗1(𝑡)=(1−𝑡)𝐺₁(𝑡),_∞_𝑖=0𝛾^∗1,𝑖𝑡^𝑖=(1−𝑡)^∞_𝑖=0𝛾_1,𝑖𝑡^𝑖.(7.4) Expanding the equation yields𝛾^∗1,0+𝛾^∗1,1𝑡+𝛾^∗1,2𝑡²+⋯=11+𝑡+𝑡²+⋯𝛾^1,0𝛾+𝛾_1,1𝑡+𝛾_1,2𝑡²+⋯_∗1,0,+𝛾^∗1,1𝑡+𝛾^∗1,2𝑡²+⋯1+𝑡+𝑡²+⋯=𝛾_1,0+𝛾_1,1𝑡+𝛾_1,2𝑡²+⋯.(7.5) This gives the recursive relationship^𝑘_𝑖=0𝛾^∗1,𝑖=𝛾_1,𝑘.(7.6) For second-order coefficient, 𝐺^∗2(𝑡)=−11!log(11−𝑡)−1!𝐺^∗1(𝑡)log(1−𝑡).(7.7) It can be written as𝐺^∗2(𝑡)=(1−𝑡)1!log(11−𝑡)−1!𝐺^∗1(𝑡)(1−𝑡)log(1−𝑡).(7.8) Substituting (7.4) into the equation above gives𝐺^∗2(𝑡)=(1−𝑡)1!log(11−𝑡)−1!(1−𝑡)𝐺¹(𝑡)(1−𝑡)log(1−𝑡)(7.9)or𝐺^∗2(𝑡)=(1−𝑡)1!log(11−𝑡)−1!𝐺¹(𝑡)log(1−𝑡).(7.10) Substituting (4.5) into (7.10) gives𝐺^∗2(𝑡)=(1−𝑡)𝐺₂(𝑡),_∞_𝑖=0𝛾^∗2,𝑖𝑡^𝑖=(1−𝑡)^∞_𝑖=0𝛾_2,𝑖𝑡^𝑖.(7.11) Expanding the equation, we have𝛾^∗2,0+𝛾^∗2,1𝑡+𝛾^∗2,2𝑡²+⋯=11+𝑡+𝑡²+⋯𝛾^2,0𝛾+𝛾_2,1𝑡+𝛾_2,2𝑡²+⋯_∗2,0,+𝛾^∗2,1𝑡+𝛾^∗2,2𝑡²+⋯1+𝑡+𝑡²+⋯=𝛾_2,0+𝛾_2,1𝑡+𝛾_2,2𝑡²+⋯.(7.12) The above gives the relationship^𝑘_𝑖=0𝛾^∗2,𝑖=𝛾_2,𝑘.(7.13) For third-order coefficient, we have𝐺^∗3(𝑡)=−12!log(11−𝑡)−2!𝐺^∗2(𝑡)log(1−𝑡).(7.14) It can be written as𝐺^∗3(𝑡)=(1−𝑡)2!log(11−𝑡)−2!𝐺^∗2(𝑡)(1−𝑡)log(1−𝑡).(7.15) Substituting (7.10) into (7.15) gives𝐺^∗3(𝑡)=(1−𝑡)2!log(11−𝑡)−2!(1−𝑡)𝐺²(𝑡)(1−𝑡)log(1−𝑡)(7.16)

or𝐺^∗3(𝑡)=(1−𝑡)2!log(11−𝑡)−2!𝐺²(𝑡)log(1−𝑡).(7.17) Substituting 𝐺₃(𝑡)=12!log(11−𝑡)−2!𝐺²(𝑡)log(1−𝑡)(7.18)into (6.15) leads to𝐺^∗3(𝑡)=(1−𝑡)𝐺₃(𝑡),_∞𝑘_𝑖=0𝛾^∗3,𝑖𝑡^𝑖=(1−𝑡)^∞_𝑖=0𝛾_3,𝑖𝑡^𝑖.(7.19) Expanding the equation into,𝛾^∗3,0+𝛾^∗3,1𝑡+𝛾^∗3,2𝑡²+⋯=11+𝑡+𝑡²+…𝛾^3,0𝛾+𝛾_3,1𝑡+𝛾_3,2𝑡²+⋯_∗3,0,+𝛾^∗3,1𝑡+𝛾^∗3,2𝑡²+⋯1+𝑡+𝑡²+…=𝛾_3,0+𝛾_3,1𝑡+𝛾_3,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.

Table 1 

The explicit integration coefficients for 𝑘 from 0 to 6.

Table 2 

The implicit integration coefficients for 𝑘 from 0 to 6.

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.

Table 3 

Table 4 

Table 5 

Table 6 

Table 7 

Table 8 

Table 9 

Problem 1. 𝑦^=2𝑦^−𝑦,0≤𝑦𝑥(0≤)=640.,𝑦^(0)=1,(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𝑒⁶⁴ which is considered large and therefore not suitable for absolute error test and hence the large values of the error.

Problem 2. 𝑦(0𝑦)=^=1,2𝑦𝑦^_−4,0≤𝑥≤30,(0)=2,𝑦^(0)=6.(8.4)

Solution 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−05≤𝑦₂𝑥+≤3,16𝜋𝑦,^2₁=𝑦^1+2𝑦₂,(0)=0,𝑦^1(0)=0,𝑦₂(0)=1.(8.6)

Solution 3. 𝑦₁𝑦(𝑥)=2cos𝑥+6sin𝑥−2−6𝑥,₂(𝑥)=−2cos𝑥+2sin𝑥+3.(8.7)

Source: Bronson [9].

For this problem, all error tests worked well.

Problem 4. 𝑦𝑦^1=1−𝑦₂−𝑦₃,0𝑦≤^2𝑥=≤𝑦₃4−𝜋,₁𝑦₁,𝑦^3=𝑦^1+𝑦^2,(0)=0,𝑦^1(0)=1,𝑦₂(0)=0,𝑦^2(0)=0,𝑦₃(0)=−1,𝑦^3(0)=1.(8.8)

Solution 4. 𝑦𝑦₁(𝑥)=sin𝑥,₂𝑦(𝑥)=cos𝑥−1,₂(𝑥)=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𝑥²𝑦^+1𝑦(1)=26𝑥,1≤𝑥≤50,21ln²(2)+99104,𝑦^(1)=−4021ln(2)−153,𝑦^(1)=326+47ln(2).(8.10)

Solution 5. 𝑦(𝑥)=𝑥²82ln𝑥2−3133−23ln(2)+31−2261ln𝑥2ln(2)+3236.(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=−𝑦¹𝑟³𝑦,𝑦₁(0)=1,𝑦^1(0)=0,_2=−𝑦²𝑟³𝑟,=𝑦₂𝑦(₂₁0)=0,𝑦^2(0)=1,0≤𝑥≤+𝑦1²²6𝜋.^1/2,(8.12)

Solution 6. 𝑦₁(𝑥)=cos𝑥,𝑦₂(𝑥)=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. 𝑦𝑦(^0)==2𝑦1³,,𝑦0_≤𝑥≤5,(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.

Acknowledgments

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