Approximating non linear higher order ODEs by a three point block algorithm

Differential equations are commonly used to model various types of real life applications. The complexity of these models may often hinder the ability to acquire an analytical solution. To overcome this drawback, numerical methods were introduced to approximate the solutions. Initially when developing a numerical algorithm, researchers focused on the key aspect which is accuracy of the method. As numerical methods becomes more and more robust, accuracy alone is not sufficient hence begins the pursuit of efficiency which warrants the need for reducing computational cost. The current research proposes a numerical algorithm for solving initial value higher order ordinary differential equations (ODEs). The proposed algorithm is derived as a three point block multistep method, developed in an Adams type formulae (3PBCS) and will be used to solve various types of ODEs and systems of ODEs. Type of ODEs that are selected varies from linear to nonlinear, artificial and real life problems. Results will illustrate the accuracy and efficiency of the proposed three point block method. Order, stability and convergence of the method are also presented in the study.


Introduction
The multistep method was discovered by Bashforth and Adams [1] in their pursue to extend Euler's method. The idea which was then named the Adams-Bashforth method was formulated by obtaining the approximated solution for a point by way of solution values from multiple previous steps. Milne [2,3] then established a new form of multistep method, known as the predictor-corrector formulation. The modern multistep method was widely researched by authors such as [4][5][6][7][8][9].
Krogh [6], managed to revitalize the field of numerical method for solving ordinary differential equations (ODEs) which was almost discarded as robust, with a divided difference approach. The premise of his work is that the back values of any point of the derivative could be interpolated. In the same study, Krogh presented a comparison for two second order a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 problems similar to the problem provide in [5], where Gear had introduced a Nordsieck version of the multi step method.
The current study is part of research series motivated by [9]. Suleiman [9] developed an algorithm for solving stiff and non stiff higher order ordinary differential equations. The nonstiff portion of the algorithm was derived in a divided difference formulation whereas the stiff algorithm was red obtained using a backward differentiation formulation. Suleiman [9] introduced a divided difference formulation known as Direct Integration (DI) method for solving nonstiff higher order ODEs directly. Omar [10] then established a new variation of the DI method by implementing a block algorithm. The block DI method developed by Omar consists of two different red algorithms, namely explicit block method and implicit block method algorithms. The block method developed in Omar [10] was matched by a fully implicit block method established by Abdul Majid [11]. In [12], authors presented a predictor-corrector algorithm in backward difference form for solving higher ODEs. The backward difference formulation from [12] was then extended in [13]. Md Ijam [13] fitted the backward difference formulae with a two point block algorithm, reducing the number of function evaluations by half.
The current study is established to overcome drawbacks of the DI method as well to provide an efficient numerical multistep method. When establishing a numerical method, accuracy has constantly been the benchmark to determine its viability. But, as numerical methods become more and more robust, accuracy alone is longer sufficient hence, the need for efficiency. The efficiency of a numerical method is associated with computational cost which translates into accuracy per step. The current study proposes a three-point block multistep formulation (3PBCS) in backward difference form for solving higher order ODEs directly. The three-point formulation requires only a fraction of computational steps usually required by standard methods. Complemented with an Adams equivalent predictor-corrector algorithm, the 3PBCS method is able to reduce computational cost more significantly. The drawback of the 3PBCS method, is the need for establishing three set of coefficients for each blocks but, by establishing multiple recursive relationships between explicit and implicit coefficients and also for coefficients of different orders, the corrector algorithm can be written in terms of the predictor which eliminates redundancy in the form of unnecessary calculations. This issue becomes irrelevant if introduced with a parallel algorithm.

Establishing the three-point block method
Ordinary differential equations are often used to model the occurrence of natural phenomenon to man made mechanics. Among these applications include modeling two body motions, chemical reactions and engineering problems such as the bend of a thin clamped beam. These ODE models can be categorized into stiff and nonstiff. However, the study conducted here focuses mainly on initial value nonstiff ODEs of any order. We begin by considering the higher order ODE, with the initial conditions In a three point block formulation, derivation of all three points are necessary. But, to avoid unwanted redundancy, only derivation of the third point will be elaborated. Comprehensive derivation of the first and second points can be obtained from [28] and [13]. When formulating a predictor-corrector three-point block algorithm, obtaining the explicit and implicit integration coefficients is crucial. The proposed three-point block method follows similar formulation to the Adam-Basthforth-Moulton method.

The predictor
The third point predictor is derived by integrating (1) once with the limit of integration from t i to t i+3 which is then denoted as Then, we approximate ϕ(t, f, f 0 , . . ., f (n−1) ) using the Newton Gregory backward difference polynomial and substituting dt = hds, thus allowing Eq (2) to be rewritten as Next, the integral in (3) is substituted by which provides the following first order predictor formulation This is then followed by derivation of the second order predictor formulation. The second order predictor is established by integrating (1) twice, which results to the following Again, we substitute ϕ(t, f, f 0 , . . ., f (n−1) ) with the Newton Gregory backward difference polynomial then replacing (t − t i+3 ) by h(3 − s) where s is as previously defined, yields Then, replacing the integral gives For the n th order integration, (1) is integrated n fold, Implementing similar step as the previous order, yields The corrector Derivation of the corrector also begins with integrating (1) once, which introduces the corrector as follows Again, approximating ϕ(t, f, f 0 , . . ., f (n−1) ) using the Newton Gregory backward difference polynomial with some subtle differences, and substituting dt = hds, the corrector in the following form Let denote b c 3;1;j ¼ ðÀ 1Þ then by amending (7) we have, Because derivation of subsequent orders follow the same sequence, we skip ahead to the n th order derivation. The derivation continues with integrating (1) n number of times thus yielding after applying similar process as the preceding orders.

Integration coefficients
In a three point block predictor-corrector method, it is highly beneficial to obtain the relationship between integration coefficients. In this section, firstly the explicit and implicit integration coefficients are derived. Then a recursive relationship between explicit and implicit integration coefficients and also coefficient of different orders are established.

Explicit coefficients
To derive the explicit integration coefficient, we first denote the first order generating function, G p 3;1 ðtÞ as By replacing b p 3;1;j t j as defined in (4), G p 3;1 ðtÞ takes the following form and by solving the integral establishes the following generating function By way of (12), the generating functions can be written in terms of b p 3;1;j

Implicit coefficients
The implicit integration is derived by first defining the first order implicit generating function hence extracting the following set of coefficients The second order implicit generating function is defined by which is then translated into the following set of coefficients The n th order generating function can then be mathematically deduced as which subsequently produces the n th order implicit coefficients

Recursive relationship
Calculating both predictor and corrector can be expensive, especially when the calculation requires large number of integration. By obtaining a recursive relationship between explicit and implicit coefficients, the corrector can be written in terms of predictor which will reduce the need for extensive calculation to obtain the corrected value. We begin by establishing the recursive interrelationship between explicit and implicit coefficients of the first order. From (13) and (18), the first order explicit and implicit generating functions are given respectively as Next, consider rearranging the implicit first order generating function G c 3;1 ðtÞ which then can be denoted as and then simplified as Then substituting the set of first order explicit and implicit coefficients obtained in (11) and (17) into (21) establishes the following which can reformulated yielding the following explicit-implicit relationship.
We then continue with the second order explicit and implicit coefficients beginning with the generating functions which can be obtained respectively from (16) and (19) as Using the relationship obtained from (21), G c 3;2 ðtÞ can be written as and by way of (16), then denoted as which yields a similar relationship as (22) Finally via similar approach, the n th order generating function can be mathematically deduced as with the recursive relationship for n th order explicit-implicit integration coefficient denoted as

Order, stability and convergence
Conditions that satisfy order and zero stability of block methods have been researched by authors such as [29][30][31]. Order and stability of the proposed method presented in this study uses techniques provided in [30]. Before establishing order and stability of the method, here are preliminary definitions that are required.

Preliminaries
This section provides necessary definitions to determine issue of stability, convergence and order of the method. We begin with the definition of the general linear multistep method.

Definition 1 The general linear multistep method is denoted by
Next, lets define the linear differential operator for the general linear multistep method as Definition 2 The linear differential operator L associated with the linear multistep method is defined Let g(t) be a function that is q times differentiable. Next, expand g(t + ih) and g 0 (t + ih) about t and arrange as Definition 3 The linear multistep method and associated difference operator L are said to be of order p if,

Order of the method
In the current section, order of both predictor and corrector will be substantiated for k = 6 number of back values. For simplicity, the example selected is of f 0 . Beginning with the predictor, f i+b for b = 1, 2, 3 of the three block method can be expressed as The order of the method follows Definition 1, where (23) to (25) are defined as matrices which satisfy where the matrices obtained are as follows The coefficients matrices are then split into sets of vector columns By way of Definition 3, from the set of vector columns yield Next, the order method for the corrector begins by denoting the corrector as follows Again by Definition 3 and derived in similar fashion as the predictor, it is established that the corrector formula is of order 6 with the error constant

Zero stability
Zero stability governs certain conditions which dictates the viability of a linear multistep method. The stability of the three point block method can be established by following similar conditions as presented in [7]. To establish whether the method is zero stable, derivation begins with the predictor. By way of the standard linear test problem f 0 ¼ lf :

PLOS ONE
and applying Eqs (23) to (25) and (26) to (28), the predictor and corrector can reformulated as and 1 0 0 respectively. Next, consider the stability polynomial which is denoted by,

PLOS ONE
By way of the stability polynomial, (29) and (30) can consequently be expressed as For zero stability, we let Λ = 0 which yields the stability polynomial for both predictor and corrector as r p ðt; 0Þ ¼ r c ðt; 0Þ ¼ t 6 À t 5 with roots t 0 = . . . = t 5 = 0, t 6 = 1. Thus, Definition 4 dictates that both predictor and corrector formulae are zero stable.

Convergence of the backward difference method
Works of [32] highlights that certain conditions are required in order to validate the convergence of a backward difference formulation. Those conditions are governed by the following. (6) and (10) to be convergent are i. Methods must be consistent

ii. Methods must be zero stable
For proof of the theorem, we refer readers to works of [32]. The order of the method was established in earlier subsection as following: The predictor

PLOS ONE
The corrector Hence, by Definition 5 the predictor and corrector are consistent of orders 5 and 6 respectively. And as shown in the previous subsection, both methods are zero stable ergo satisfying the necessary conditions for convergence. Finally, we proceed with the subsequent section, the long awaited numerical results.

Results and discussion
Real life science and engineering problems are often modeled in the form of ODEs. Every so often, these problems are not able to be solved analytically, thus requires numerical approximations. These numerical approximations are often used as alternative solutions due to the absence of analytical solution. The current section provides results for selected higher order ODEs using the proposed three point block method (3PBCS) with constant step size. Problems selected varies from linear to nonlinear artificial and real life problems which are in the form of single and systems of equations. The problems selected also consist of various difficulty levels in order to provide a holistic overview of the 3PBCS method's capability. For a more just comparison, the 3PBCS method is compared against Direct Integration (DI), one point block (1PBCS) and two point block (2PBCS) methods. Results will be compared in terms of accuracy, function evaluations and efficiency. The error solution used in this study adopts techniques suggested in [12] which is given by where the n th component of the exact solution is denoted by (f i ) n and the n th component of the approximated solution for f is denoted by f(t i ) n . The error Err n is used to define three types of error test used this study, absolute error (A = 1, B = 0), relative error (A = 0, B = 1) and mixed error (A = 1, B = 1) tests. Before continuing with the numerical results, these are a few terms that will be used throughout the section.

Higher order problems
This section begins with higher order differential problems which are artificial in nature and of different orders then followed with well known differential equations. Examples 1-3 are 4 th to 5 th order problems obtained from various sources.
with the initial condtions f| 1 = 1, f 0 | 1 = −1, f@| 1 = 2, f´´´| 1 = −6f (iv) | 1 = 24 and given analytical solution f ðtÞ ¼ 1 t . (Source: [16]) Table 1 presents numerical results of the 1PBCS, 2PBCS and 3PBCS methods for Examples 1 -3. The 1PBCS approximates single steps of equidistant, the 2PBCS approximates two-points of equidistant simultaneously and the 3PBCS approximates three-points of equidistant simultaneously. Results will compare accuracy and total steps approximated by all three methods for step sizes 10 −1 , 10 −2 , 10 −3 , 10 −4 and 10 −5 . Tstep represents the total steps required and Err denotes the maximum error estimated for each method in the interval T 0 � t � T n . The 1PBCS and 2PBCS methods were selected to compare efficiency of the 3PBCS when solving higher order ODEs with other backward difference methods to provide a fair analysis. The analysis of Examples 1 to 3 begins with results provided in Table 1. As shown in Table 1, the 3PBCS method obviously requires less calculation steps compared to the latter even though some loss of accuracy is expected. For Example 2, the 3PBCS is able to match the level of accuracy of 1PBCS and 2PBCS methods but slightly under-performed in Examples 1 and 3. Whereas Table 2 provides computational time (Time) which is recorded in seconds for Examples 1-3. As exhibited in Table 2, 3PBCS requires less computational time compared to 1PBCS and 2PBCS methods for almost all step sizes, except at H = 1 × 10 −1 . This is due to the initial calculation of the three-point block integration coefficients. The small amount of Tsteps required to calculate the solutions from beginning to end outweighs the efficiency of the 3PBCS method. But as a smaller H is selected, more Tstep is required which redistributes calculation time of the initial 3PBCS integration coefficients thus leads to a lesser overall computational time. This loss is barely significant and can be overcome in a variable order stepsize algorithm. Figs 1-3 compares the efficiency of the methods. As defined in [12], the efficiency

PLOS ONE
of the methods is the accuracy per step and by way of the under most curve, graphs depicted in Figs 1-3 show the 3PBCS method to be slightly more efficient than latter methods. Example 4: A 3 rd order systems of ODE, with the initial conditions f 1 | 0 = 1, f 0  Table 3 compares results between the DI and 3PBCS methods. Opposed to the 3PBCS backward difference formulation, the DI method is derived with a divided difference formulation. The DI method was chosen to test the 3PBCS with a competing predictor-corrector multistep

PLOS ONE
method. Results presented highlights maximum error (ErrEQ1, ErrEQ2, ErrEQ3) for each of equation, f 1 , f 2 and f 3 and the overall maximum error for each step size, ErrMax. The error is obtained by comparing approximated solution with the exact solution, j ðf i Þ n À f ðt i Þ n AþBðf ðt i Þ n Þ j with conditions established in (31). The DI was selected because of its divided difference formulation. The purpose was to test the performance of the 3PBCS method against an alternative multistep method using similar parameter. In the current case, both methods were set using similar order methods (same number of back values). Results of the current example clearly illustrates the advantage of 3PBCS over the DI method.
Example 5: The problem presented is a artificial second order ODE with periodic solution. The second order ODE with the initial condtions f| 0 = −1, f 0 | 0 = 50 and given analytical solution f(t) = sin(50t) − cos (60t) provides an unique oscillatory solution. Table 4

Application problems
Test problems in the current section are practical ODEs that are found in real applications including two body motion, thin flow and electrical circuits. Lets begin with the renowned two body motion. Example 6: Newton's equation of the two body motion problem refers to the movement of two particle interacting which each other due to gravitational pull, neglecting any third body the two do not collide with (see Fig 8)  In the current research, we consider the following formulation of the two body problem f @ 1 ¼ À f 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi with the initial conditions f 1 j 0 ¼ 1 f 0 1 j 0 ¼ 0; f 2 j 0 ¼ 0; f 0 2 j 0 ¼ 1 and given analytical solution f 1 (t) = cos(t), f 2 (t) = sin(t). (Source: [36]) Data provided in Table 5 are numerical approximations for Example 6 using DI and 3PBCS method. Results exhibits that when using a larger step size, both methods are comparable but as smaller step size are used, 3PBCS begins showing an obvious advantage. Results produced in the table shows that the accuracy of 3PBCS improves at a rate of H 2 compared to e the DI method which barely equals the current step size. To further validate the capability of the 3PBCS method, Figs 9 dan 10 show graphical approximation of f 1 and f 2 by NDSOLVER and 3PBCS methods. It is clear that both methods practically overlap each other, point by point. Fig 11 is presented to illustrate the orbit of the two body motion as approximated by 3PBCS. Example 6, with the current initial conditions presents a circular orbit (orbit 1), where as, by changing the initial conditions to f 1 j 0 ¼ 0: gives Example 6 a Kepler-like solution which produces an elliptic orbit (orbit 2) and can be verified from the work of Shampine and Gordon [37].

PLOS ONE
These thin film flow problems can be found in various form like the the infamous drainage dry surface given by the function Consider a thin film pre-wetted surface with a remotely small thickness, ω > 0. This changes the function to Among the main concern is the tension surface effects when dealing with a flow of thin film of viscous fluid with free surface. As discussed in [38] f 000 ¼ f À k with initial values f| 0 = , f 0 | 0 = 1, f 0 | 0 = 1. (Source: [39]) exemplify the dynamic balance against the strength of a thin fluid layer, disregarding gravity.
Because the thin flow problem have no exact or analytical solution, the solution by 3PBCS is compared against the solution obtained by the DI. Table 6 presents the approximated solution obtained by 3PBCS method of f, f 0 and f@ for the points 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0. As shown in Table 6, both methods are comparable with similar approximation up to at least 5 decimal points. The objective is to find the capacitor charge at any time, t > 0 given the initial current and initial charge of the capacitor are both zero. By Kirchhoff law, the following second order differential with the initial conditions f| 0 = 0, f 0 | 0 = 0 and given analytical solution f ðtÞ ¼ 6e À 10t 61 ð6 sin ð50tÞ þ 5 cos ð50tÞÞ À 5 61 ð5 sin ð60tÞ þ 6 cos ð60tÞÞ. (Source: [40]) The following Table 7 contains numerical results of 3PBCS method for example. The approximated values for end points (T n ) 0.5, 1.0, 1.3 and 2.0 are provided with the corresponding exact solutions. The maximum error of all points, ErrMax is defined by T 0 � |f(t n+b ) − f(t)| � T n , b = 1, 2, 3 is also provided for a clearer view of 3PBCS abilities. As demonstrated in Table 6, 3PBCS is consistently accurate to the seventh decimal point when applying a stepsize of H = 1 × 10 −3 .
is defined as the displacement of the periodic solution and f| 0 is the amplitude of the oscillation (Source: [41]). The classical nonlinear dynamics of a self sustained oscillation is commonly modeled when ξ > 0. Parameters used for this problem are ξ = 0.025 for 0 � t � 40 with initial conditions f| 0 = 0, f 0 | 0 = 0.5.
Example 9 is originally a problem obtained from [42]. The Van Der Pol equation selected does not have any known analytical solution. Table 8

PLOS ONE
is to present a 3D illustration of f, f 0 and f@ which corresponds to similar steps points by the 3PBCS method.

Conclusion
Results provided in the current work validate that 3PBCS is a viable numerical method for solving ODEs. The three-point element of the method allows for lesser computational cost and provides better efficiency compared to its rival counterparts. By simply using a constant step size algorithm decreases computational cost and increases efficiency. Tables 1, 3 and 5 validate the convergence of the method. When a smaller H is selected, the method becomes more accurate and proven to provide a consistent set of accuracy for every point in the interval as indicated in Tables 4 and 7. For future works, the 3PBCS can be fitted with a variable order step size and parallel computing algorithm. The effectiveness of a variable order step size algorithm depends on the selection of a suitable order step size criterion. Authors a currently refining appropriate conditions to provide a better approximation for a variable order step size algorithm.