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′, …, f⁽ⁿ⁻¹⁾) 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′, …, f⁽ⁿ⁻¹⁾) 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 where

The corrector

Derivation of the corrector also begins with integrating (1) once, which introduces the corrector as follows

Again, approximating ϕ(t, f, f′, …, f⁽ⁿ⁻¹⁾) using the Newton Gregory backward difference polynomial with some subtle differences, and substituting dt = hds, the corrector in the following form Let denote 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.

Explicit coefficients

To derive the explicit integration coefficient, we first denote the first order generating function, $$ G_{3,1}^p\left(t\right)$$ as

By replacing $$ {\beta }_{3,1,j}^p{t}^j$$ as defined in (4), $$ G_{3,1}^p\left(t\right)$$ 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 $$ {\beta }_{3,1,j}^p$$ Expanding the functions in (14) then gives the following set of first order explicit coefficients, Next, the second order generating function has similar form as (11), where and by substituting $$ {\beta }_{3,2,j}^p{t}^j$$ with (5) we can rewrite $$ G_{3,2}^p\left(t\right)$$ as Then, solving (15) yields thus establishing By mathematical induction, the n^th order explicit generating function and set of coefficients can be constructed as and

Implicit coefficients

The implicit integration is derived by first defining the first order implicit generating function

Subsequent to defining the generating function, $$ G_{3,1}^c\left(t\right)$$ is then equated as with the solution hence extracting the following set of coefficients The second order implicit generating function is defined by where $$ {\beta }_{3,2,j}^c$$ is replaced with thus re-defining the generating function as Next, $$ G_{3,2}^c\left(t\right)$$ is expressed in the form of 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

and Next, consider rearranging the implicit first order generating function $$ G_{3,1}^c\left(t\right)$$ 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 and Using the relationship obtained from (21), $$ G_{3,2}^c\left(t\right)$$ 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

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

where g(t) is an arbitrary function in C¹[a, b].

Let g(t) be a function that is q times differentiable. Next, expand g(t + ih) and g′(t + ih) about t and arrange as

where C_i, i = 0, 1.…, q, … are constants.

Definition 3 The linear multistep method and associated difference operator L are said to be of order p if, C₀ = C₁ = … = C_q = 0, C_q+1 ≠ 0 where

Definition 4 The block method is zero stable if the roots r_j, j = 1, …k of the first characteristic polynomial ρ(r) denoted by satisfies the conditions |r_j| ≤ 1 and the roots with |r_j| = 1 where the multiplicity does not exceed 2.

Definition 5 A Linear Multistep Method is said to be consistent if the LLM is of order p ≥ 1.

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′. 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 and The coefficients matrices are then split into sets of vector columns

By way of Definition 3, from the set of vector columns yield

thus, establishing the predictor is of order 5 with a corresponding error constant of 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

and applying Eqs (23) to (25) and (26) to (28), the predictor and corrector can reformulated as and respectively. Next, consider the stability polynomial which is denoted by, By way of the stability polynomial, (29) and (30) can consequently be expressed as and For zero stability, we let Λ = 0 which yields the stability polynomial for both predictor and corrector as with roots t₀ = … = t₅ = 0, t₆ = 1. Thus, Definition 4 dictates that both predictor and corrector formulae are zero stable.

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.

Example 1: A non homogeneous 4^th order ODE,

with the initial conditions f|₀ = 0, f′|₀ = 1, f″|₀ = 0, f‴|₀ = −1 and given analytical solution f(t) = sin(t). (Source: [33])

Example 2: A 4^th order non linear ODE

with the initial condtions f|₀ = 1, f′|₀ = −0.1, f″|₀ = 0.02, f‴|₀ = −0.006 and given analytical solution $$ f\left(t\right)={\displaystyle \frac{10}{10+t}}$$. (Source: [34])

Example 3: A 5^th order non linear ODE

with the initial condtions f|₁ = 1, f′|₁ = −1, f″|₁ = 2, f‴|₁ = −6f^(iv)|₁ = 24 and given analytical solution $$ f\left(t\right)={\displaystyle \frac{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⁻¹, 10⁻², 10⁻³, 10⁻⁴ and 10⁻⁵. Tstep represents the total steps required and Err denotes the maximum error estimated for each method in the interval T₀ ≤ 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⁻¹. 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 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, $$ f_1^{\prime }{|}_0=-1$$, $$ f_1^{″}{|}_0=1$$, f₂|₀ = 1, $$ f_2^{\prime }{|}_0=-2$$, $$ f_2^{″}{|}_0=4$$, f₃|₀ = 1, $$ f_3^{\prime }{|}_0=-3$$ and $$ f_3^{″}{|}_0=9$$ and given analytical solutions f₁(t) = e^−t, f₂(t) = e^−2t and f₃(t) = e^−3t. (Source: [35])

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 method. Results presented highlights maximum error (ErrEQ1, ErrEQ2, ErrEQ3) for each of equation, f₁, f₂ and f₃ and the overall maximum error for each step size, ErrMax. The error is obtained by comparing approximated solution with the exact solution, $$ \left|{\displaystyle \frac{{\left(f_i\right)}_n-f{\left(t_i\right)}_n}{A+B\left(f{\left(t_i\right)}_n\right)}}\right|$$ 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|₀ = −1, f′|₀ = 50 and given analytical solution f(t) = sin(50t) − cos(60t) provides an unique oscillatory solution.

Table 4 provides numerical approximation by the 3PBCS and the provided exact solution for Example 5. Due to the difficulty of the current example, the step size of order H = 10⁻⁵ is selected. Results show a consistent approximation error of 1×10⁻⁶ for every point. To further validate the accuracy of 3PBCS as suitable method for the approximation of a non linear ODE with periodic solution, 3PBCS is compared against NDSOLVER, a preset Euler method equipped in Mathematica 12. The plotted graph displayed in Fig 4 shows that the 3PBCS method matches the solution obtained by NDSOLVER thus, confirming the accuracy of the method. Figs 5 and 6 are 2D and 3D parametric plots of Example 5 by the 3PBCS. For a clear understanding of the difficulty level for Example 5, readers may also refer to a 3D representation of the periodic solutions as illustrated in Fig 7.

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

with the initial conditions $$ f_1{|}_0=1f_1^{\prime }{|}_0=0,f_2{{|}_0=0,f_2^{\prime }|}_0=1$$ and given analytical solution f₁(t) = cos(t), f₂(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² 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₁ and f₂ 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{|}_0=0.4f_1^{\prime }{|}_0=0,f_2{{|}_0=0,f_2^{\prime }|}_0=2$$, 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].

Example 7: Many problems in engineering and physics are based on the famous thin flow of liquid which is a third order ODEs with the ideal draining flow as shown in Fig 12.

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] with initial values f|₀ = , f′|₀ = 1, f′|₀ = 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′ 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. Fig 13 compares solutions by NDSOLVER and 3PBCS. For the current example, 3PBCS method approximates the solution using a step size of H = 1 × 10³. Even at the current stepsize, the figure shows 3PBCS ability to match step by step solution of the NDSOLVER. Whereas, Fig 14 illustrates the difference of solution between f, f′ and f″.

Example 8: Consider the second order RLC circuit in Fig 15 below, with resistance of 2 ohm, capacitor at $$ {\displaystyle \frac{1}{260}}$$ and inductor at 0.1 Henry and electrical force of 100 sin 60tV. 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 equation can be derived,

with the initial conditions f|₀ = 0, f′|₀ = 0 and given analytical solution $$ f\left(t\right)={\displaystyle \frac{6e^{-10t}}{61}}(6\text{sin}\left(50t\right)+5\text{cos}\left(50t\right))-{\displaystyle \frac{5}{61}}(5\text{sin}\left(60t\right)+6\text{cos}\left(60t\right))$$. (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₀ ≤ |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⁻³. Fig 16 is the approximate solution for f_n by NDSOLVER and 3PBCS. The graphical overlap shown in the figure exhibits the similar accuracy of both methods, whereas Fig 17 is parametric plot of Example 8 presented to express the level of difficulty of the problem.

Example 9: The Van der Pol equation is an ODE that is derived from Rayleigh differential equation. In a Van der Pol equation,

f(t) is defined as the displacement of the periodic solution and f|₀ 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, f′|₀ = 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 provides approximated values of f, f′ and f″ for 5 different end points by the DI and 3PBCS methods. Numerical approximation shows that both methods obtained similar values for f, f′ and f″ up to 1 × 10⁻⁵. Figs 18 and 19 illustrate the approximated solution, f by the 3PBCS method against NDSOLVER for the interval 0 ≤ t ≤ 40. Results in Fig 18 shows similar approximation for both 3PBCS and NDSOLVER where Fig 19 illustrates a graphical representation approximated values of t against f, f′ and f″. Fig 20 is to present a 3D illustration of f, f′ and f″ which corresponds to similar steps points by the 3PBCS method.