Implicit Three-Point Block Numerical Algorithm for Solving Third Order Initial Value Problem Directly with Applications

: Recently, direct methods that involve higher derivatives to numerically approximate higher order initial value problems (IVPs) have been explored, which aim to construct numerical methods with higher order and very high precision of the solutions. This article aims to construct a fourth and ﬁfth derivative, three-point implicit block method to tackle general third-order ordinary differential equations directly. As a consequence of the increase in order acquired via the implicit block method of higher derivatives, a signiﬁcant improvement in efﬁciency has been observed. The new method is derived in a block mode to simultaneously evaluate the approximations at three points. The derivation of the new method can be easily implemented. We established the proposed method’s characteristics, including order, zero-stability, and convergence. Numerical experiments are used to conﬁrm the superiority of the method. Applications to problems in physics and engineering are given to assess the signiﬁcance of the method.


Introduction
A wide variety of real life situations are represented by mathematical models as third order ordinary differential equations (ODEs), such as chemical engineering, biology, electromagnetic waves, quantum mechanics, the motion of rocket, and thin film flow [1][2][3][4]. Nevertheless, the theoretical solutions for most of these equations are undefined; therefore, third-order ODEs have gained significant attention and the need to develop numerical methods with more accurate approximations is eminent [5][6][7][8]. In the classical way, solving higher order ODEs is done by reducing the equation into an equivalent system of first-order ODEs, but this process is too rigorous compared to the direct methods [9][10][11]. Not only that, but it is also found that the implementation process of the direct methods is simpler and more accurate than the process of reduction [11]. In order to avoid the reduction effort, many researchers have proposed different methods to solve initial value problems (IVPs) of the ODEs directly [12][13][14][15][16]. To enhance the efficacy of numerical methods, many researchers developed block methods by producing the r-point of the approximate solutions simultaneously. Kuboye and Omar have presented a direct seven-step block method for solving third-order ODEs by using a multistep collocation technique [17]. Awoyemi et al. [18] developed a direct continuous five step collocation method for solving the general third order IVPs. An implicit continuous linear multistep methods using the interpolation and collocation for solving the general third order IVPs was proposed by [6]. Normally, direct methods are constructed by interpolation and collocation strategies, which is complicated due to the process of determining the coefficients of the method. Where the points need to be collocated and interpolated after which a system of linear equations must be resolved. Moreover, scholars have looked lately at methods that involved derivatives in solving ODEs, which lead to a more accurate numerical results and to increase the order of the methods [8,[19][20][21]. In fact, the accuracy of a method increases with the increase in the order of the method. However, the idea of incorporation of higher derivatives of the solution in the process leads to higher and better accuracy and is achieved without a corresponding increase in the order of the method. Therefore, in this research, our main concern is to propose an implicit three-point block method with fourth and fifth derivatives of the solution by using a technique that can be implemented in a straightforward manner for directly solving both linear and nonlinear problems of general third-order ODEs in the form of y = f (t, y, y , y ) y(a) = y 0 , y (a) = y 0 , y (a) = y 0 , We assume that f in Equation (1) is differentiable to a desired order in region R and f (t, y, y , y ) satisfies the Lipchitz condition in its second, third and fourth terms as for all points (t, y i , y , y ), (t, y, y i , y ), and (t, y, y , y i ); i = 1, 2 in the region R. Then the IVPs in Equation (1) have a unique solution in R (see [22,23] ).
In the upcoming section we will derive the three-point implicit block method of order nine (ITPBO9) and will discuss the basic idea of how the block method works. In Section 3 the main properties of the suggested method are analyzed. The implementation of the method is presented in Section 4. Section 5 presents the discussion on the numerical experiments as well as on applications to equations of fluid flow, such as problems in thin film flow, the boundary layer equation and the nonlinear Genesio equation. Lastly, the conclusion of the research is provided in Section 6.

Methodology
The three-point block method generates three approximate values, y n+1 , y n+2 and y n+3 concurrently at t n+1 , t n+2 and t n+3 respectively, using one earlier block, where t n becomes the starting point and t n+3 is the last point in the block with step size 3h. The method is derived by applying numerical integration thrice to Equation (1) to acquire the approximate formula of y n+1 , y n+2 and y n+3 .
Integrating the first, second and third point once gives: Integrating the first, second and third point thrice gives: In order to derive the formula, we need to approximate f (t, y, y , y ) in Equation (1) using Hermite interpolation P n (t) [24]: where n is the degree of the Hermite polynomial.
For the first point, y n+1 , let s = t−t n+1 h and dt = hds be substituted into (2), (5) and (8). By evaluating the integral from −3 to −2 using MAPLE gives the following where g and q denote the fourth and fifth derivatives of the solution, respectively. Next, we introduce s = t−t n+2 h and dt = hds for the second point, y n+2 , into (3), (6) and (9) and apply a similar technique by evaluating the integral from −2 to −1 using MAPLE gives the following: Then, for the third point, we introduce s = t−t n+3 h and dt = hds into (4), (7) and (10) and apply a similar technique by evaluating the integral from −1 to 0 using MAPLE, which gives the following

Order and Error Constant
We can write the Equations (12)-(20) in a matrix difference equation to recognize the order of the proposed method as where α i ; i = 1, 2, ..., 6 defined as, The linear operator related to Equation (21) can be defined as By expanding Equation (22) in the Taylor series yields T . Therefore, we concluded that the new method has order 9. As the method's order is p ≥ 1, then, it can be said that the method is consistent (see Lambart and Fatunla [25,26]).

Zero Stability
To check the zero-stability of the implicit three-point method, we rewrite the formulas into a matrix form as below (1) and S (0) are constant coefficients and A (0) = 9 × 9 identity matrix y n−2 y n−2 y n−2 y n−1 y n−1 y n−1 y n y n y n Then, the first characteristic polynomial can be written as By solving (23), we obtain the roots R = 0, 0, 0, 0, 0, 0, 1, 1, 1. Since |R| ≤ 1, the three-point implicit block method is zero stable. Along with the consistency of the method, this property implies convergence of the new method (see Ackleh et al. [27]).

Implementation
The implementation of the three-point implicit block method on general third-order ODEs is carried out in a straightforward manner by applying the predictor-corrector schemes. The implementation begins by using Taylor's method as the predictor to compute the starting values. When the starting values are obtained, ITPBO9 will be implemented as the corrector to estimate the approximate solutions of y , y and y at t. In the iteration, we evaluate functions f , g, q at each point, which will be used to compute the approximate solutions at the next point. The procedure to solve the third-order ODEs by using the new method can be observed in the following Algorithm 1.

Algorithm 1
The procedure to solve the third-order ODEs by using the new method.  13: If t i+3 < b, then repeat step 6. Else, go to step 14. 14: Evaluate the maximum error, which is defined as MAXE = max 1≤i≤N | y(t i ) − y i | . 15: Execute the results. Complete.

Results and Discussion
In this section, some well-known single and system of linear and nonlinear IVPs as well as applications of the IVPs are presented, which aim to assess the efficiency and accuracy of the new ITPBO9 method of order nine compared to other direct block methods. C + + programming codes have been developed and applied to solve IVPs in ODE of the form (1) based on the proposed ITPBO9 method. The results are compared with other existing methods with similar characteristics and order to give an idea of how well the new method performs and to clearly display the efficiency of the ITPBO9 method. The following abbreviations will be used in the tables: ITPBO9: Implicit three-point block direct method introduced in this paper of order nine. HCD: Block hybrid collocation direct method of order six [28]. ABAM: Adams Bashforth-Adams Moulton method of order four. FSM: Five-step direct method of order nine [18]. ILMM: Implicit linear multistep direct method of order six [6]. ISHD: Three-step hybrid direct method of order nine [5]. h : Step size. NS: Number of steps. AE: Absolute error at the point considered.

MAXE:
Maximum absolute error on the grid points at the interval.
Exact solution: Problem 4. Consider the nonlinear system y 1 = 1 2 e 4t y 3 y 2 , y 1 (0) = 1, y 1 (0) = −1, y 1 (0) = 1, y 2 = 8 3 e 2t y 1 y 3 , y 2 (0) = 1, y 2 (0) = −2, y 2 (0) = 4, Exact solution:  [18] of order nine with regards to the accuracy. Table 2 shows a direct comparison between the new ITPBO9 method with HCD [28] and the well-known fourth order Adams Bashforth-Adams Moulton (ABAM) method in terms of the number of steps and accuracy at different step sizes. ITPBO9 reduces the number of steps to one third compared to the ABAM method and requires the same number of steps compared to HCD since ITPBO9 and HCD compute three points simultaneously; nevertheless, Figure 1 displayed the best performances and efficiency of the new ITPBO9 method.
In addition, we have solved the linear system in Problem 2 in order to compare the proposed ITPBO9 method with ISHD [5] of order nine and ILMM [6] of order six, which are presented in Tables 3 and 4, respectively. In Table 3, we have considered the maximum absolute errors for step size h = 1 2 j , j = 2, 3, 4, 5. While in Table 4, the numerical results have been obtained by considering the maximum absolute errors along the interval using a different number of steps where NS refers to the number of steps. It can be observed that increasing the number of steps leads to a decrease in the error. Based on the efficiency curves in Figures 1 and 2, the new ITPBO9 method is more efficient where it achieved a smaller maximum error compared to the other methods when the step size h decreases as well as when the number of steps taken is the same. Therefore, Tables 3 and 4 and Figures 1 and 2 clearly point out how ITPBO9 is superior in terms of accuracy and efficiency compared to the existing methods.
However, we have considered in Problem 3 and Problem 4 single and systems of nonlinear IVPs. Tables 5-8 show the computed solutions compared to the exact solutions and the absolute errors at different points. It is remarkable that the computed solutions at each point t agree very well with the exact solutions up to 16 digits at least.
Overall, from all the tables and figures, the convergence and high precision of the new method is clearly observed. The method is more efficient than the existing methods, in which the order is either nearly equal or identical to that of the new method at each point t and each step size h.

Application to Thin Film Flow Problem
We also consider the well-known engineering and physical problem, which is the thin film flow of a liquid on a surface. This problem was studied by several authors [5,7,[29][30][31]. They studied the fluid motion on a plane surface where the motion along the plane is in the same flow direction. The fluid dynamics problem is governed by the third-order ODEs where y(t) moves with the fluid in a coordinate frame and f (y(t)) can vary depending on the physical context, such as is a formula for a fluid draining problem on a dry surface.
is a formula for a fluid draining problem on a wet surface where ξ is the film thickness and ξ > 0. Problems regarding the flow of thin films with a free surface of viscous fluid in which surface tension effects play a role commonly leads to third order ordinary ODEs governing the shape of the free surface of the fluid, which is given by subject to y(t 0 ) = a 1 , where a 1 , a 2 and a 3 are constants. In the literature, several researchers solve the problem (27) with the initial conditions y(0) = y (0) = y (0) = 1, for the case µ = 2. A parametric representation for the exact solution of (27) was introduced by [30]. For comparison purposes, we will apply the new method to solve (27) directly. The results are in Table 9 and Figure 3 clearly shows that the method is applicable and the solution that was obtained by ITPBO9 agrees very well with the exact solutions [30].  Figure 3. Response curve concerning Equation (27) with µ = 2, h = 0.01.

Application to Boundary Layer Equation
A boundary layer in physics and fluid mechanics is the fluid layer in the immediate vicinity of a bounding surface in which the viscosity has significant effects. The equation of the nonlinear boundary layer is a nonlinear differential equation of third order defined as subject to It is a well known equation as the Blasius equation, which describes a boundary layer flow over a flat plate. This equation was already considered in [32,33]. We applied our method to solve the Blasius equation and to determine the shear stress at the plate. Figure 4 depicts that the solution to Equation (29) of the proposed method ITBPO9 in the interval t ∈ [0, 10] with h = 0.1 agree very well with approximations found by the Mathematica built-in package NDSolve.

Application to Nonlinear Genesio Equation
Consider the following nonlinear Genesio equation, which was introduced as a chaotic system by Genesio [34] y + αy + βy − f (y(t)) = 0, where f (y(t) = −γy(t) + y(t) 2 , subject to y(t 0 ) = 0.2, where α, β and γ are positive constants satisfying αβ < γ. The theoretical solution for this problem is unknown. This problem was studied by some researchers such as Bataineh et al. [35], which included the behavior of this system. Table 10 shows the computed solutions and the number of steps by the proposed ITPBO9 method, HCD [28] and the NDSolve at different b as well as different step sizes. We apply the new method to solve the nonlinear Genesio equation when α = 1.2, β = 2.92, γ = 6. It can be observed that ITPBO9 is applicable to solve Equation (29) with an advantage of fewer total steps compared to NDSolve. Figure 5 illustrates the numerical approximations for Equation (29) with h = 0.1 in t ∈ [0, 4]. It is obvious that the solutions obtained by ITPBO9 agree very well with approximations found by the Mathematica built-in package NDSolve that evaluate the efficiency of the new method.

Conclusions
In this article, we proposed a three-point implicit block method using the fourth and fifth derivatives of the solution, which aim to solve linear and nonlinear single as well as system initial value problems of the general third-order ODEs directly. The method is also applicable to solve the physical and engineering problems of the general third-order ODEs directly. The idea of incorporation of higher derivatives of the solution in the process, is that higher and better accuracy can be achieved without a corresponding increase in the order of the method. This new method is uncomplicated to implement and satisfies the property of convergence, which is indicated by the significant improvement with regards to accuracy in the numerical results. Therefore, we suggest the new ITPBO9 method as a suitable tool for solving general third-order ODEs directly with high precision and easy implementation.

Conflicts of Interest:
The authors declare no conflict of interest.