Direct Two-Point Block One-Step Method for Solving General Second-Order Ordinary Differential Equations

Abstract

A direct two-point block one-step method for solving general second-order ordinary differential equations (ODEs) directly is presented in this paper. The one-step block method will solve the second-order ODEs without reducing to first-order equations. The direct solutions of the general second-order ODEs will be calculated at two points simultaneously using variable step size. The method is formulated using the linear multistep method, but the new method possesses the desirable feature of the one-step method. The implementation is based on the predictor and corrector formulas in the 𝑃𝐸(𝐶𝐸)^𝑚 mode. The stability and precision of this method will also be analyzed and deliberated. Numerical results are given to show the efficiency of the proposed method and will be compared with the existing method.

1. Introduction

In this paper, we are considering solving directly the general second-order initial value problems (IVPs) for systems of ODEs in the form𝑦^′′=𝑓𝑥,𝑦,𝑦^′,𝑦(𝑎)=𝑦₀,𝑦^′(𝑎)=𝑦^′₀,𝑥∈[𝑎,𝑏]. (1.1) Equation in (1.1) arises from many physical phenomena in a wide spectrum of applications especially in the science and engineering areas such as in the electric circuit, damped and undamped spring mass and some other areas of application. We will also consider solving general second order as in (1.1) using the direct block one-step method. Block methods for numerical solutions of ODEs have been proposed by several researchers such as in [1–5]. The common block methods used to solve the problems can be categorized as one-step block method and multistep block method.

One-step block method such as the implicit Runge-Kutta method is also being referred to as one previous point to obtain the solution. The multistep block method in the form of Adams type formula is presented in [5, 6]. In [7], the block backward differentiation formula (BBDF) for solving stiff ODEs has been introduced and the solutions referred to as more than one previous point. The works in [6] showed the proposed two-point four-step block method presented as in a simple form of Adams Moulton method for solving second-order ODEs directly.

The block method of Runge-Kutta type has been explored in [1], and it is suggested that a block of new approximation values is used simultaneously for solving first-order ODEs. The works in [3, 8, 9] have been considered in solving (1.1) using the block one-step method, while [3] has proposed a two-point implicit block one-step method for solving second-order ODEs directly and suggested that the method is suitable to be parallel.

In [9], Majid et al. have derived the two-point block method for solving first-order ODEs by using the closest point in the interval, that is, [𝑥_𝑛,𝑥_𝑛+1] and [𝑥_𝑛+1,𝑥_𝑛+2]. The Gauss Seidel iteration was implemented in the proposed block method. The approach in this research is to extend the idea in [9] for solving (1.1) directly without reducing system of first-order ODEs using two-point block one-step method.

2. Formulation of the Method

In order to compute the two approximation values of 𝑦_𝑛+1 and 𝑦_𝑛+2simultaneously, the interval of [𝑎,𝑏] is divided into a series of blocks with each block containing two points as shown in Figure 1.

Figure 1 

Block one-step method.

In Figure 1, we observed that 𝑘th block contains 𝑥_𝑛,𝑥_𝑛+1 and 𝑥_𝑛+2, where 𝑥_𝑛 becomes the starting point and 𝑥_𝑛+2 is the last point in the 𝑘th block with step size ℎ. The approximations values of 𝑦_𝑛+1 and 𝑦_𝑛+2 are computed simultaneously. The evaluation solution at the last point in 𝑘th block will be restored as the initial values for (𝑘+1)th block. The same procedure is used to compute the solutions for the next block until the end of the interval. The evaluation information from the previous step in a block could be used for other steps of the same block only. During the calculations of iteration, the final values of 𝑦_𝑛+2 at the point 𝑥_𝑛+2 are taken as the initial values for the next iteration.

We obtained the approximation values of 𝑦_𝑛+1 and 𝑦_𝑛+2 at the points 𝑥_𝑛+1 and 𝑥_𝑛+2 by integrating once and twice over (1.1) with respect to𝑥. We begin to evaluate the 𝑦^′_𝑛+1 and 𝑦_𝑛+1 by integrating once and twice (1.1) over the interval [𝑥_𝑛,𝑥_𝑛+1]:^𝑥𝑛+1_𝑥𝑛𝑦^′′(𝑥)𝑑𝑥=^𝑥𝑛+1_𝑥𝑛𝑓𝑥,𝑦,𝑦^′𝑑𝑥,_{𝑥𝑛+1𝑥𝑛𝑥𝑥𝑛}𝑦^′′(𝑥)𝑑𝑥𝑑𝑥=^𝑥𝑛+1_𝑥𝑛^𝑥_𝑥𝑛𝑓𝑥,𝑦,𝑦^′𝑑𝑥𝑑𝑥. (2.1) Let 𝑥_𝑛+1=𝑥_𝑛+ℎ which gives 𝑦^′𝑥_𝑛+1−𝑦^′𝑥_𝑛=^𝑥𝑛+1_𝑥𝑛𝑓𝑥,𝑦,𝑦^′𝑦𝑥𝑑𝑥,_𝑛+1−𝑦𝑥_𝑛−ℎ𝑦^′𝑥_𝑛=^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑥𝑓𝑥,𝑦,𝑦^′𝑑𝑥. (2.2) Then, 𝑓(𝑥,𝑦,𝑦^′) in (2.2) will be replaced with Lagrange interpolation polynomial that involves interpolation points in the block at (𝑥_𝑛,𝑓_𝑛),(𝑥_𝑛+1,𝑓_𝑛+1) and(𝑥_𝑛+2,𝑓_𝑛+2). Taking 𝑥=𝑥_𝑛+2+𝑠ℎ,𝑠=(𝑥−𝑥_𝑛+2)/ℎ, and 𝑑𝑥=ℎ𝑑𝑠 and replacing into (2.2). Then, the limit of integration in (2.2) will be from −2 to −1. The corrector formulae will be obtained using MATHEMATICA. The formulae of 𝑦^′_𝑛+1and 𝑦_𝑛+1 are obtained as follows:𝑦^′_𝑛+1=𝑦^′_𝑛+ℎ125𝑓_𝑛+8𝑓_𝑛+1−𝑓_𝑛+2,𝑦_𝑛+1=𝑦_𝑛+ℎ𝑦^′_𝑛+ℎ²247𝑓_𝑛+6𝑓_𝑛+1−𝑓_𝑛+2. (2.3) To approximate the value of 𝑦^′_𝑛+2 and 𝑦_𝑛+2, we take 𝑥_𝑛+2=𝑥_𝑛+2ℎ by integrating once and twice (1.1) over the interval [𝑥_𝑛+1,𝑥_𝑛+2]and apply the same process. The limit of integration will be from −1 to 0, and the following corrector formulae will be obtained:𝑦^′_𝑛+2=𝑦^′_𝑛+1+ℎ12−𝑓_𝑛+8𝑓_𝑛+1+5𝑓_𝑛+2,𝑦_𝑛+2=𝑦_𝑛+1+ℎ𝑦^′_𝑛+1+ℎ²24−𝑓_𝑛+10𝑓_𝑛+1+3𝑓_𝑛+2. (2.4) The formulae (2.3) and (2.4) may be rewritten in the form of matrix difference equation as follows:0000𝑦0−110000000−11^𝑛−1𝑦_𝑛𝑦_𝑛+1𝑦_𝑛+2𝑦=ℎ01−100100001−10010^′𝑛−1𝑦^′_𝑛𝑦^′_𝑛+1𝑦^′_𝑛+2+ℎ²0581212−1071262424−1024−1812512012−12410324𝑓24^𝑛−1𝑓_𝑛𝑓_𝑛+1𝑓_𝑛+2. (2.5) The order of this developed method is identified by referring to [10–12]. The two-point one-step block method for ODEs can be written in a matrix difference equation as follows:𝛼𝑌_𝑚=ℎ𝛽𝑌^′_𝑚+ℎ²𝛾𝐹_𝑚, (2.6) where 𝛼,𝛽, and 𝛾 are the coefficients with the m-vector 𝑌_𝑚, 𝑌^′_𝑚, and 𝐹_𝑚 be defined as𝑌_𝑚=𝑦_𝑛−1,𝑦_𝑛,𝑦_𝑛+1,𝑦_𝑛+2^𝑇,𝑌^′_𝑚=𝑦^′_𝑛−1,𝑦^′_𝑛,𝑦^′_𝑛+1,𝑦^′_𝑛+2^𝑇,𝐹_𝑚=𝑓_𝑛−1,𝑓_𝑛,𝑓_𝑛+1,𝑓_𝑛+2^𝑇. (2.7)

By applying the formulae for the constants 𝐶_𝑞, in [10] the formulae is defined as𝐶₀=^𝑘_𝑗=0𝛼_𝑗,𝐶₁=^𝑘_𝑗=0𝑗𝛼_𝑗−𝛽_𝑗,𝐶₂=^𝑘_𝑗=0𝑗²𝛼2!^𝑗−𝑗𝛽_𝑗−𝛾_𝑗,⋮𝐶_𝑞=^𝑘_𝑗=0𝑗^𝑞𝛼𝑞!^𝑗−𝑗^𝑞−1!𝛽𝑞−1^𝑗−𝑗^𝑞−2!𝛾𝑞−2^𝑗,where𝑞=3,4,5,…. (2.8) Therefore, the order and error constant of the two-point one-step block method can be obtained by using equation in (2.8).

For 𝑞=0,𝐶₀=³_𝑗=01𝑗0!⁰𝛼_𝑗=𝛼₀+𝛼₁+𝛼₂+𝛼₃=0000+000+010+0001=0000−1−1. (2.9)

For 𝑞=1,𝐶₁=³_𝑗=01𝑗1!¹𝛼_𝑗−³_𝑗=01𝑗0!⁰𝛽_𝑗=1¹⋅𝛼₁+2¹⋅𝛼₂+3¹⋅𝛼₃−𝛽₀+𝛽₁+𝛽₂+𝛽₃0000100001−0000+1100+011+000=0000=1.−1+2−1+3−1−1 (2.10) For 𝑞=2,𝐶₂=³_𝑗=01𝑗2!²𝛼_𝑗−³_𝑗=01𝑗1!¹𝛽_𝑗−³_𝑗=01𝑗0!⁰𝛾_𝑗=121²⋅𝛼₁+2²⋅𝛼₂+3²⋅𝛼₃−1¹⋅𝛽₁+2¹⋅𝛽₂+3¹⋅𝛽₃−𝛾₀+𝛾₁+𝛾₂+𝛾₃=1210000100001−11100011000−0000+5−1+4−1+9+2−1+3−171224−112−1+8246128241210+24−112−1524312=0000.24 (2.11)

For 𝑞=3,𝐶₃=³_𝑗=01𝑗3!³𝛼_𝑗−³_𝑗=01𝑗2!²𝛽_𝑗−³_𝑗=01𝑗1!¹𝛾_𝑗=0000. (2.12)

For 𝑞=4,𝐶₄=³_𝑗=01𝑗4!⁴𝛼_𝑗−³_𝑗=01𝑗3!³𝛽_𝑗−³_𝑗=01𝑗2!²𝛾_𝑗=0000. (2.13)

For 𝑞=5,𝐶₅=³_𝑗=01𝑗5!⁵𝛼_𝑗−³_𝑗=01𝑗4!⁴𝛽_𝑗−³_𝑗=01𝑗3!³𝛾_𝑗=1124−145−724≠0000360. (2.14) The method is order 𝑝 if 𝐶₀=𝐶₁=⋯=𝐶_𝑝+1=0 and 𝐶_𝑝+2≠0 is the error constant. Thus, we conclude that the method in (2.3) and (2.4) is of the order 3 and the error constant is 𝐶_𝑝+2=𝐶₅=1,124145,−724,−360^T. (2.15)

3. Implementation of the Method

The initial starting point at each block are obtained by using Euler method as predictor and the initial ℎ used is as follows:ℎ_initial=TOL2^1/(𝑝+1), (3.1) where TOL is the tolerance and 𝑝 is the order of the method. Then, the calculations are corrected using the corrector formulae in (2.3) and (2.4). For the next block, the same techniques are repeated to compute the approximation values of 𝑦_𝑛+1 and 𝑦_𝑛+2 simultaneously until the end of the interval. We use 𝑟 function evaluations per step at the corrector formulae, that is, 𝑃𝐸(𝐶𝐸)^𝑟, where P, E, and 𝐶 indicate the predictor, evaluate the function 𝑓, and the corrector, respectively.

Algorithm 3.1. Computing approximating 𝑦^′(𝑝)_𝑛+1,0,𝑦^(𝑝)_𝑛+1,0,𝑦^′(𝑝)_𝑛+2,0, and 𝑦^(𝑝)_𝑛+2,0 using the predictor formulas is as follows: for𝑖=1 to 2 do𝑃∶𝑦^′(𝑝)_𝑛+𝑖,0=𝑦₍₎^′(𝑝)_{𝑛+𝑖−1}+ℎ𝑓₍₎^(𝑝)_{𝑛+𝑖−1},𝑦^(𝑝)_𝑛+𝑖,0=𝑦₍₎^(𝑝)_{𝑛+𝑖−1}+ℎ𝑦₍₎^′(𝑝)_{𝑛+𝑖−1}+ℎ²2𝑓₍₎^(𝑝)_{𝑛+𝑖−1},𝐸∶𝑓𝑥_𝑛+𝑖,𝑦^(𝑝)_𝑛+𝑖,0,𝑦^′(𝑝)_𝑛+𝑖,0 (3.2)end forComputing approximations 𝑦^′(𝑐)_𝑛+1,0,𝑦^(𝑐)_𝑛+1,0,𝑦^′(𝑐)_𝑛+2,0, and𝑦^(𝑐)_𝑛+2,0 using corrector formulas is as follows: for 𝑟=0do𝐶∶𝑦^′(𝑐)_{𝑛+1,𝑟+1}=𝑦_𝑛^′(𝑐)+ℎ125𝑓_𝑛^(𝑐)+8𝑓^(𝑝)_𝑛+1,0−𝑓^(𝑝)_𝑛+2,0,𝑦^(𝑐)_{𝑛+1,𝑟+1}=𝑦_𝑛^(𝑐)+ℎ𝑦_𝑛^′(𝑐)+ℎ²247𝑓_𝑛^(𝑐)+6𝑓^(𝑝)_𝑛+1,0−𝑓^(𝑝)_𝑛+2,0,𝑦^′(𝑐)_{𝑛+2,𝑟+1}=𝑦^′(𝑐)_{𝑛+1,𝑟+1}+ℎ12−𝑓_𝑛^(𝑐)+8𝑓^(𝑝)_𝑛+1,0+5𝑓^(𝑝)_𝑛+2,0,𝑦^(𝑐)_{𝑛+2,𝑟+1}=𝑦^(𝑐)_{𝑛+1,𝑟+1}+ℎ𝑦^′(𝑐)_{𝑛+1,𝑟+1}+ℎ²24−𝑓_𝑛^(𝑐)+10𝑓^(𝑝)_𝑛+1,0+3𝑓^(𝑝)_𝑛+2,0, (3.3)    for𝑖=1 to 2 do𝐸∶𝑓𝑥_𝑛+𝑖,𝑦^(𝑐)_{𝑛+𝑖,𝑟+1},𝑦^′(𝑐)_{𝑛+𝑖,𝑟+1} (3.4)   end forend forfor 𝑟=1,2,…until convergence do𝐶∶𝑦^′(𝑐)_{𝑛+1,𝑟+1}=𝑦_𝑛^′(𝑐)+ℎ125𝑓_𝑛^(𝑐)+8𝑓^(𝑐)_𝑛+1,0−𝑓^(𝑐)_𝑛+2,0,𝑦^(𝑐)_{𝑛+1,𝑟+1}=𝑦_𝑛^(𝑐)+ℎ𝑦_𝑛^′(𝑐)+ℎ²247𝑓_𝑛^(𝑐)+6𝑓^(𝑐)_𝑛+1,0−𝑓^(𝑐)_𝑛+2,0,𝑦^′(𝑐)_{𝑛+2,𝑟+1}=𝑦^′(𝑐)_{𝑛+1,𝑟+1}+ℎ12−𝑓_𝑛^(𝑐)+8𝑓^(𝑐)_𝑛+1,0+5𝑓^(𝑐)_𝑛+2,0,𝑦^(𝑐)_{𝑛+2,𝑟+1}=𝑦^(𝑐)_{𝑛+1,𝑟+1}+ℎ𝑦^′(𝑐)_{𝑛+1,𝑟+1}+ℎ²24−𝑓_𝑛^(𝑐)+10𝑓^(𝑐)_𝑛+1,0+3𝑓^(𝑐)_𝑛+2,0, (3.5)for𝑖=1 to 2 do𝐸∶𝑓𝑥_𝑛+𝑖,𝑦^(𝑐)_{𝑛+𝑖,𝑟+1},𝑦^′(𝑐)_{𝑛+𝑖,𝑟+1} (3.6)end for end for

The convergence test:𝑦_{𝑛+2,𝑟+1𝑡}−𝑦_{𝑛+2,𝑟𝑡}𝐴+𝐵𝑦_{𝑛+2,𝑟+1𝑡}<0.1×TOL, (3.7) where 𝑟 is the number of iterations and (𝑦)_𝑡is the tth component of the approximate. 𝐴=1, 𝐵=0 correspond to the absolute error test. 𝐴=1, 𝐵=1 correspond to the mixed test and finally 𝐴=0, 𝐵=1 correspond to the relative error test. The mixed error test is used in all tested problems. If the convergence test is satisfied, then we control the error in the current block by performing the local truncation error as follows:LTE=𝑦^(𝑘)_𝑛+2−𝑦^(𝑘−1)_𝑛+2, (3.8) where 𝑘 is the order of corrector formula. The errors calculated in the code is defined as𝐸_𝑖𝑡=𝑦_𝑖𝑡−𝑦𝑥_𝑖𝑡𝐴+𝐵𝑦𝑥_𝑖𝑡. (3.9) For the evaluation of maximum error, it is defined as follow:MAXE=max_{1≤𝑖≤TS}max_{1≤𝑖≤𝑁}𝐸_𝑖𝑡, (3.10) where 𝑁 is the number of equation in the system and TS is the number of successful steps.

In order to make the selection of the next step size, we follow again the techniques used in [9]. If LTE ≤ TOL, then the next step size remains constant or double, otherwise the step size will be half. The step size when the integration steps are successful is given byℎ_new=𝐶×ℎ_old×TOLLTE^1/(𝑘+1),ifℎ_new≥2×ℎ_old,thenℎ_next=2×ℎ_old,elseℎ_next=ℎ_old, (3.11) where 𝐶= 0.5 is a safety factor. The algorithm when the step failure occurs isℎ_next=12×ℎ_old. (3.12)

4. Stability of the Method

In this section, we will discuss the stability of the proposed method derived in the previous section on a linear general second-order problem:𝑦^′′=𝑓=𝜃𝑦^′+𝜆𝑦. (4.1)

Firstly, the test equation in (4.1) is substituted into the predictor formula (3.2). The evaluation of 𝑓^𝑝_𝑛+1,0and 𝑓^𝑝_𝑛+2,0 will be substituted into the right-hand side of (3.3) when 𝑟=0 as shown below:𝑓^𝑝_𝑛+1,0=𝜃𝑦^′𝑝_𝑛+1,0+𝜆𝑦^𝑝_𝑛+1,0,𝑓^𝑝_𝑛+2,0=𝜃𝑦^′𝑝_𝑛+2,0+𝜆𝑦^𝑝_𝑛+2,0, (4.2) where 𝑦^′𝑝_𝑛+1,0, 𝑦^𝑝_𝑛+1,0, 𝑦^′𝑝_𝑛+2,0, and 𝑦^𝑝_𝑛+2,0in (4.2) are obtained in (3.2) after substituted the test equation (4.1). Then, the formulae are written in the matrix form and setting the determinant of the matrix to zero. Hence, the stability polynomial is obtained:𝑟=0: 𝑡²1+25ℎ𝜃+129ℎ²𝜃²4+5ℎ³𝜃³6+ℎ⁴𝜆⁴3+ℎ²𝜆8−ℎ³𝜃𝜆2+ℎ⁴𝜃²𝜆−ℎ12₅𝜃³𝜆−1237ℎ⁴𝜆²−ℎ144⁵𝜃𝜆²9−ℎ⁶𝜃²𝜆²9−ℎ⁶𝜆³+ℎ36⁷𝜃𝜆³+ℎ16⁸𝜆⁴36+𝑡³−2−13ℎ𝜃6−7ℎ²𝜃²3−5ℎ³𝜃³6−ℎ⁴𝜆⁴3−17ℎ²𝜆4−47ℎ³𝜃𝜆−127ℎ⁴𝜃²𝜆4−7ℎ⁵𝜃³𝜆−12115ℎ⁴𝜆²−7219ℎ⁵𝜃𝜆²−187ℎ⁶𝜃²𝜆²−187ℎ⁶𝜆³−3617ℎ⁷𝜃𝜆³−ℎ144⁸𝜆⁴72+𝑡⁴1+ℎ𝜃+ℎ12²𝜃²+ℎ12²𝜆8+ℎ³𝜃𝜆+ℎ12⁴𝜆²48=0. (4.3)Figure 2 showed the stability region of the direct block one-step method when 𝑟=0.

Figure 2 

Stability region for direct block one-step method when 𝑟=0.

In case 𝑟=1, the stability polynomial is obtained as follows:𝑟=1:𝑡²1+2ℎ𝜃+25ℎ²𝜃²+1219ℎ³𝜃³+ℎ12⁴𝜃⁴3+5ℎ⁵𝜃⁵+ℎ18⁶𝜃⁶9+7ℎ48³𝜃𝜆−13ℎ48⁴𝜃²1𝜆+4ℎ⁵𝜃³𝜆+17ℎ72⁶𝜃⁴1𝜆+ℎ18⁷𝜃⁵𝜆+11ℎ⁴𝜆²−28883ℎ144⁵𝜃𝜆²+101ℎ432⁶𝜃²𝜆²+13ℎ432⁷𝜃³𝜆²−17ℎ432⁸𝜃⁴𝜆²−199ℎ⁶𝜆³−172853ℎ288⁷𝜃𝜆³−7ℎ216⁸𝜃²𝜆³−5ℎ432⁹𝜃³𝜆³−29ℎ⁸𝜆⁴−32417ℎ⁹𝜃𝜆⁴+129697ℎ¹⁰𝜃²𝜆⁴−ℎ5184¹⁰𝜆⁵+1648ℎ96¹¹𝜃𝜆⁵+ℎ¹²𝜆⁶648+𝑡³−2−2ℎ𝜃−13ℎ²𝜃²6−5ℎ³𝜃³3−ℎ⁴𝜃⁴3−5ℎ⁵𝜃⁵−ℎ18⁶𝜃⁶9−4ℎ²𝜆−103ℎ24³−𝜃𝜆17ℎ⁴𝜃²𝜆6−11ℎ⁵𝜃³𝜆−1219ℎ⁶𝜃⁴𝜆−245ℎ⁷𝜃⁵𝜆−18203ℎ⁴𝜆²−144199ℎ144⁵𝜃𝜆²−365ℎ⁶𝜃²𝜆²−432361ℎ432⁷𝜃³𝜆²−127ℎ⁸𝜃⁴𝜆²−432133ℎ⁶𝜆³−864269ℎ⁷𝜃𝜆³−86423ℎ⁸𝜃²𝜆³−5473ℎ⁹𝜃³𝜆³−43255ℎ⁸𝜆⁴−1296139ℎ⁹𝜃𝜆⁴−1296289ℎ¹⁰𝜃²𝜆⁴−51847ℎ¹⁰𝜆⁵−64813ℎ¹¹𝜃𝜆⁵−ℎ1296¹²𝜆⁶1296+𝑡⁴ℎ1+²𝜃²+ℎ12³𝜃³+712ℎ48³5𝜃𝜆+ℎ48⁴𝜃²𝜆+11ℎ⁴𝜆²+1288ℎ24⁵𝜃𝜆²+ℎ⁶𝜆³192=0. (4.4) Figure 3 showed the stability region of the direct block one-step method when 𝑟=1.

Figure 3 

Stability region for direct block one-step method when 𝑟=1.

The stability region is plotted using MATHEMATICA, and the shaded region inside the boundary in Figures 2 and 3 demonstrate the stability region for the proposed method.

5. Numerical Results and Discussion

We have tested the performance of the proposed method on four problems. For the first three problems, a comparison is made between the solutions obtained in [9].

Problem 1. We have 𝑦₁^′′=−𝑦₂+sin𝜋𝑥,𝑦₁(0)=0,𝑦^′₁(0)𝑦=−1,_2′′=−𝑦₁+1−𝜋²sin𝜋𝑥,𝑦₂(0)=1,𝑦^′₂(0)=1+𝜋,[0,10]. (5.1) Solution: 𝑦₁(𝑥)=1−𝑒^𝑥,𝑦₂(𝑥)=𝑒^𝑥+sin𝜋𝑥.
First-order systems: 𝑦^′₁=𝑦₃,𝑦^′₂=𝑦₄,𝑦^′₃=−𝑦₂+sin𝜋𝑥,𝑦^′₄=−𝑦₁+1−𝜋²𝑦sin𝜋𝑥,₁(0)=0,𝑦₂(0)=1,𝑦₃(0)=−1,𝑦₄(0)=1+𝜋. (5.2) Solution:𝑦₁(𝑥)=1−𝑒^𝑥, 𝑦₂(𝑥)=𝑒^𝑥+sin𝜋𝑥, 𝑦₃(𝑥)=−𝑒^𝑥, 𝑦₄(𝑥)=𝑒^𝑥+𝜋cos𝜋𝑥.

Problem 2. We have 𝑦₁^′′=−𝑒^−𝑥𝑦₂,𝑦₁(0)=1,𝑦^′₁(0)𝑦=0,_2′′=2𝑒^𝑥𝑦^′₁,𝑦₂(0)=1,𝑦^′₂(0)=1,[0,𝜋]. (5.3) Solution: 𝑦₁(𝑥)=cos𝑥,𝑦₂(𝑥)=𝑒^𝑥cos𝑥.
First-order systems: 𝑦^′₁=𝑦₃,𝑦^′₂=𝑦₄,𝑦^′₃=−𝑒^−𝑥𝑦₂,𝑦^′₄=2𝑒^𝑥𝑦₃,𝑦₁(0)=1,𝑦₂(0)=1,𝑦₃(0)=0,𝑦₄(0)=1. (5.4) Solution: 𝑦₁(𝑥)=cos𝑥,𝑦₂(𝑥)=𝑒^𝑥cos𝑥,𝑦₃(𝑥)=−sin𝑥,𝑦₄(𝑥)=𝑒^𝑥cos𝑥−𝑒^𝑥sin𝑥.

Problem 3. We have 𝑦₁^′′=−𝑦^′₂,𝑦₁(0)=0,𝑦^′₁(0)=11−𝑒⁻¹,𝑦₂^′′=−𝑦^′₁,𝑦₂(0)=1,𝑦^′₂(0)=11−𝑒⁻¹,[0,10]. (5.5) Solution: 𝑦₁(𝑥)=(1−𝑒^−𝑥)/(1−𝑒⁻¹), 𝑦₂(𝑥)=(2−𝑒⁻¹−𝑒^−𝑥)/(1−𝑒⁻¹).
First-order systems: 𝑦^′₁=𝑦₃,𝑦^′₂=𝑦₄,𝑦^′₃=−𝑦₄,𝑦^′₄=−𝑦₃,𝑦₁(0)=0,𝑦₂(0)=1,𝑦₃(0)=11−𝑒⁻¹,𝑦₄(0)=11−𝑒⁻¹. (5.6) Solution:𝑦₁(𝑥)=(1−𝑒^−𝑥)/(1−𝑒⁻¹), 𝑦₂(𝑥)=(2−𝑒⁻¹−𝑒^−𝑥)/(1−𝑒⁻¹), 𝑦₃(𝑥)=(𝑒^−𝑥)/(1−𝑒⁻¹), 𝑦₄(𝑥)=𝑒^−𝑥/(1−𝑒⁻¹).

Problem 4 (Van Der Pol oscillator). We have 𝑦^′′−2𝜉1−𝑦²𝑦^′+𝑦=0,𝑦(0)=0,𝑦^′(0)=0.5,[0,10], (5.7) where we take 𝜉=0.025.

The codes are written in C language and executed on DYNIX/ptx operating system. The numerical results for Problems 1–3 in Tables 2–4 are solved using the proposed method and the method in [9]. The results in terms of total steps and execution times for solving Problems 1–3 are presented in histograms and graph lines in Figures 4, 5 and 6. The solutions of 2P(B) for solving Problem 4 are plotted in Figure 7 at tolerance 10⁻⁸and compared with the solutions obtained by the MATLAB built-in solver ode45.

Figure 4 

Results of total steps and time for Problem 1.

Figure 5 

Results of total steps and time for Problem 2.

Figure 6 

Results of total steps and time for Problem 3.

Figure 7 

Solution of Problem 4 when TOL = 10−8.

In Figures 4, 5 and 6, it is obvious that method 2P(B) requires less number of total steps as compared to method 2P(A) when solving the same given problems. It is also observed that the execution times of 2P(B) are faster than 2P(A) at all tested tolerances. This is expected since 2P(B) has less number of function calls; therefore, it has affected the computation time of 2P(B).

In Table 1, the RSTEP and RTIME are greater than 1.00 which shows that 2P(B) is more efficient compared to 2P(A). In fact, in some cases, the ratios are greater than 3.00, which indicates a clear advantage of method 2P(B) over 2P(A).

Table 1 

The ratios execution times and steps for solving Problems 1 to 3.

Table 2 

Comparison between 2P(A) and 2P(B) methods for solving Problem 1.

Table 3 

Comparison between 2P(A) and 2P(B) methods for solving Problem 2.

Table 4 

Comparison between 2P(A) and 2P(B) methods for solving Problem 3.

In Table 2, it can be observed that the maximum error of 2P(B) is one or two order larger than 2P(A) but still acceptable as it is within the given tolerance. This is expected since the code 2P(B) solved the given problem directly without reducing to system of first-order differential equations. We could observe that the RSTEP for Problem 1 is greater than 2 in Table 1.

In Tables 3 and 4, it is observed that the maximum error of 2P(B) is two or three order larger as compared to 2P(A) in solving Problems 2 and 3 but it is still within the given tolerances. In Figure 7, it is obvious that the 2P(B) solutions agree very well with the solutions obtained by the MATLAB built-in solver ode45 when solving Problem 4.

6. Conclusion

In this paper, we have constructed the direct two-point block one-step method which is efficient and suitable for solving general second-order ODEs directly. The block method has shown acceptable solutions and managed to solve the second-order ODE faster compared to the existing method.

Notations

TOL:	Tolerance
MTD:	Method employed
TS:	Total steps taken
FS:	Total failure steps
MAXE:	Magnitude of the maximum error of the computed solution
AVERR:	The average error
FCN:	Total function calls
TIME:	The execution time taken in microseconds
2P(A):	Implementation of the direct block method in [9] by reducing the problem to system of first-order ODEs
2P(B):	Implementation of the direct two-point block one-step method by solving the problem directly
RSTEP:	The ratio steps, TS2P(A)/TS2P(B)
RTIME:	The ratio execution times, TIME2P(A)/TIME2P(B)

Acknowledgment

The authors would like to thank the Universiti Putra Malaysia for providing financial support through Graduate Research Fellowship (GRF) during the study period.