Seventh order hybrid block method for solution of first order stiff systems of initial value problems

A hybrid second derivative three-step method of order 7 is proposed for solving first order stiff differential equations. The complementary and main methods are generated from a single continuous scheme through interpolation and collocation procedures. The continuous scheme makes it easy to interpolate at off-grid and grid points. The consistency, stability, and convergence properties of the block formula are presented. The hybrid second derivative block backward differentiation formula is concurrently applied to the first order stiff systems to generate the numerical solution that do not coincide in time over a given interval. The numerical results show that the new method compares favorably with some known methods in the literature.


Introduction
The numerical solutions of stiff systems have been one of the major worries for numerical analysts. A numerical method that is potentially good for solving systems of stiff ODEs must have some reliablity in terms of its region of absolute stability and good accuracy [1]. Consider the system of initial value problems (IVPs) given by y = f (t, y), y(t 0 ) = y 0 (1) where f : × 2s −→ s ; y, y 0 ∈ s , and s is the dimension of the system. f satisfies the Lipschitz condition, and the Jacobian (2020) 28:34 Page 2 of 11 Dahlquist, several linear multistep methods have been developed including continuous ones. Continuous multistep methods have been the subject of growing interest due to the fact that continuous methods enjoy certain advantages, such as they have the potential to provide defect control (see Enright [4]) as well as they are able to generate complementary methods, which are applied together as a single block scheme (see Onumanyi et al. [5], Akinfenwa et al. [6]- [11], and Jator [12]- [14]). Most of the block methods proposed in the literature such as Chartier [15], Shampine and Watts [16], Chu and Hamilton [17], and Suleiman et al. [18] to mention but a few are usually implemented in the predictorcorrector mode. In this paper, we propose an Hybrid Block Second Derivative Backward Differentiation Formula (HBSDBDF) which simultaneously provides the solution of (1) in each block without the use of predictors (see Akinfenwa et al. [6]- [11], Jator et al. [12], [13], and Yakubu and Markus [19]).

Development of the method
We seek the numerical estimation of the analytic solution y(t) by assuming an approximate solution Y (t) of the form where t ∈ [t 0 , T n ], m j are undetermined coefficients that must be obtained and ϕ j (t) are the basis polynomial functions of degree 7. It is required that the eight equations below must be satisfied. 7 j=0 m j t j = y n+i , i = 0, where n is the grid index, y n+i = Y (t n+i ) is the numerical estimation of the analytical solution y(t n+i ) , and g n+i = df dt (t n+i , Y (t n+i )) . Equations (3), (4), and (5) will provide a system of eight equations whose solutions generate the coefficients m j which are replaced in Eq. (2). After some algebraic manipulation, the continuous form of the hybrid second derivative formula is obtained as where α i 2 (t) , i = 0, 1, 2, . . . , 5, β k (t), and γ k (t) are continuous coefficients; k = 3 is the step number; and h is the step length. We assume that y n+i/2 = Y t n+ ih 2 is the numerical estimation of the analytical solution y t n+ ih 2 , y n+k = f t n+k , y n+k is an approximation to y t n+k , and g n+k = df dt (t, y(t))| t n+k The main method is generated from (6), by interpolating at t = (t n+3 ) to obtain While the complementary methods are obtained by differentiating (6) with respect to t to obtain and interpolating (8) The hybrid methods are implemented by applying (7), (9)-(13), as a single block method to provide the approximate solution y n+ 1 2 , y n+1 , y n+ 3 2 , y n+2 , y n+ 5 2 , y n+3 , for n = 0, 3, . . . , N − 3, N on the partition [t 0 , t 3 , t 6 , · · · , t N−3 , t N ].

Zero stability
The zero stability of HBSDBDF can be determined from the difference system (14), as h tends 0. Thus, as h → 0, the method (14) tends to the difference system (15) where W 1 and W 0 are 6 by 6 constant matrices. Hence, from (15), we obtain the first characteristic polynomial π(L) given by The HBSDBDF is zero-stable since π(L) = 0 satisfies that the roots |L j | ≤ 1 , j = 1, . . . , 6, and for those roots L j = 1, the multiplicity does not exceed 1.

Local truncation error Theorem
The HBSDBDF has a local truncation error

Consistency and convergence
The HBSDBDF (14) is consistent as each of the schemes has order p > 1. Following Henrici [20], since the HBSDBDF is zero-stable and consistent, then the HBSDBDF is convergent.
The stability domain of the method is defined as δ = {z ∈ C : (z) ≤ 1}. To determine the stability of the HBSDBDF, we state the following definitions: (2020) 28:34 Page 6 of 11 (i) Dahlquist [3]: A numerical method whose stability region contains the entire left half plane is said to be A-stable. (ii) Ehle [21]: A method that is A-stable and Lim[ R(z)] → 0 as z → −∞ is said to be L-stable.
The region of absolute stability (RAS) of the method is plotted using the boundary locus technique. Figure 1 depicts the stability region for the HBSDBDF of the dominant eigenvalue R(z). It can be seen that the HBSDBDF is A-stable because the stability region contains the whole left half complex plane whose interval [ −3.4, 0]. Also the HBSDBDF is A-stable and as in Ehle [21] the requirement that lim z→−∞ R(z) = 0 is satisfied. Thus, it is L-stable.

Numerical examples
Example 1 This example is a linear system on the range 0 ≤ t ≤ 10 (see [6]).
For the problem, the maximum absolute error was computed on the given interval. It was found that the HBSDBDF of order 7 is more accurate when compared to the BBDF in [6] of the same order. Although when h = 0.1250, the max error is greater than the error in Our aim here is to demonstrate the accuracy, rate of convergence (ROC), and good stability properties of the HBSDBDF . For different step sizes h, the relative error max i  Table 2. Clearly, the ROC is consistent with the order of the new block scheme. Example 3 Consider the highly stiff system, see [19]. y 1 = −10 7 y 1 + 0.075y 2 , y 1 (0) = 1 The eigenvalues of the Jacobian of the system are approximately λ 1 = −1.000000000562500 × 10 6 and λ 2 = −0.0743749995813. The result of the HBSDBDF is compared with that of Yakubu and Markus [19] using second derivative method and shown as displayed in Table 3.   Example 4 Next consider the chemistry problem in Gear [22] , Cash [23], and Yakubu [19], We solve this problem in the interval 0 ≤ t ≤ 50 using the HBSDBDF. The result is y 1 in blue, y 2 in brown, and y 3 in red as shown in Fig. 2 with the numerical values for h = 0.001 at the point T = 10, 20, 30, 40, and 50. See Table 4.

Example 5 Consider the well-known nonlinear problem (Kaps problem) in the range
This problem is solved using method in [10] and the new HBSDBDF so as to show the advantage of the new method over that in [10]. The absolute error at the end point t = 10, NFE the number of function evaluations, h the step size, and N the number of computation steps are displayed in Table 5.

Results and discussion
The stability of the newly derived method was obtained by using the boundary locus approach. The technique involves finding the roots of the stability function which is a  rational complex function and by plotting the imaginary root against the real root. The interval of stability read from the plot of the region of absolute stability gives [−3.4,0]. The result obtained in Example 1 showed the accuracy and stability of the method. However, when h = 0.1250, the max error is greater then the error in [1]. This was due the fact that the new method converges with correct digit of 13 from h = 0.5 to h = 0.1250. The second example was solved to show that the rate of convergence of the method is in agreement with the order of the method. The third example is a highly stiff system, and it is solved to show the effectiveness of the method. Despite the fact that the method is of order 7, it was compared with those [19] of orders 8 and 11. Table 3 shows the superiority of the new method over those in [19]. The fourth problem solved is also a standard chemistry problem and the result plotted in Fig. 2 and the numerical solution shown in Table 4 is in agreement with those in the literature. The advantage of the HBSDBDF can be seen in Table 5 where the new method converges even for very large step size. The low number of function evaluation shows that the new method can save computer memory with reduced computation time.

Conclusion
In this article, a new hybrid second derivative block backward differentiation formula for solving stiff systems of first order initial value problems is reported. The stability analysis has been conducted based on the boundary locus technique to obtain the region of absolute stability. The HBSDBDF is implemented without the need for predictors or starting values, and therefore, subroutines that are sometimes complicated are not necessary. Five standard numerical examples, both linear and non-linear stiff systems, have been solved to show the accuracy and efficiency of the methods. From the results obtained, the rate of convergence confirmed the order of the method. Detailed results are displayed in Tables  1, 2, 3, 4 and 5. The results have shown that HBSDBDF is suitable for solution of stiff problems and converges accurately even for large step size h. The advantages of the new method are that it is more accurate than those in [10] and [19] in the manuscript. It has less number of function evaluation when compared with [10] and [19], hence reduced the time of computation and reduced use in computer memory. Another advantage is that it converges accurately with large step size h while those in [10] and [19] are less accurate as evident in Tables 3 and 5. This is the goal of numerical analysts. The disadvantage however is that the new method will converge with very fewer digit/s of accuracy when compared with the method in [10] for problems using very small h.