- Original research
- Open access
- Published:
Seventh order hybrid block method for solution of first order stiff systems of initial value problems
Journal of the Egyptian Mathematical Society volume 28, Article number: 34 (2020)
Abstract
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
where \(f : \Re \times \Re ^{2s} \longrightarrow \Re ^{s}\); y, \(y_{0}\in \Re ^{s}\), and s is the dimension of the system. f satisfies the Lipschitz condition, and the Jacobian \(\left (\frac {\partial f}{\partial y}\right)\) with the negative real parts of its eigenvalues varies slowly (see [2]). These types of equations are called stiff equations, and they arise frequently in sciences and engineering. It is well known that many problems of this type do not have analytic solutions; therefore, numerical methods remain important. Recall that Eq. (1) can be expeditiously solved using methods that are A-stable, and to obtain high accuracy, methods with higher order are preferred. But according to Dahlquist theorem (Theorem 2.2, see [3]), high order linear multistep methods (LMMs) for solving (1) cannot have order greater than 2. To overcome this barrier theorem of 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 predictor-corrector 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\in \left [t_{0}, T_{n}\right ]\), mj 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.
where n is the grid index, \(y_{n+i} = Y\left (t_{n+i}\right)\) is the numerical estimation of the analytical solution y(tn+i), and \(g_{n+i} = \frac {df}{dt}{\left (t_{n+i}, Y(t_{n+i})\right).}\)
Equations (3), (4), and (5) will provide a system of eight equations whose solutions generate the coefficients mj which are replaced in Eq. (2). After some algebraic manipulation, the continuous form of the hybrid second derivative formula is obtained as
where \(\alpha _{\frac {i}{2}(t)}, i=0,1,2,\ldots,5, \beta _{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\left (t_{n+\frac {ih}{2}}\right)\) is the numerical estimation of the analytical solution \(y\left (t_{n+\frac {ih}{2}}\right), y'_{n+k}= f\left (t_{n+k}, y_{n+k}\right)\) is an approximation to \(y'\left (t_{n+k}\right),\) and \(g_{n+k} = \frac {df}{dt}{(t, y(t))}|_{t_{n+k}}\)
The main method is generated from (6), by interpolating at \(t=\left (t_{n+3}\right)\) to obtain
While the complementary methods are obtained by differentiating (6) with respect to t to obtain
and interpolating (8) at the points \(t=\left (t_{n+\frac {i}{2}}\right), i=1,2,\ldots,5\) to yield
The hybrid methods are implemented by applying (7), (9)–(13), as a single block method to provide the approximate solution \( y_{n+\frac {1}{2}}, y_{n+1},y_{n+\frac {3}{2}}\), \(y_{n+2}, y_{n+\frac {5}{2}},y_{n+3}\), for n=0,3,…,N−3,N on the partition \(\left [t_{0},t_{3},t_{6},\cdots,t_{N-3},t_{N}\right ]\).
Analysis of HBSDBDF
This section is better discussed by first representing Eqs. (7), (9)–(13), with a matrix difference equation as
where
for σ=0,1,2,… and n=0,3,6,…,N−3.
Zero stability
The zero stability of HBSDBDF can be determined from the difference system (14), as h tends 0. Thus, as \(h\rightarrow 0\), the method (14) tends to the difference system (15)
where W1 and W0 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 |Lj|≤1, j=1,…,6, and for those roots \(\left |L_{j}\right | = 1\), the multiplicity does not exceed 1.
Local truncation error
Theorem
The HBSDBDF has a local truncation error (LTE) of \(C_{8}h^{8}y^{8}(t_{n})+\bigcirc \left (h^{9}\right)\).
Proof
Suppose \(y\left (t_{n}\right)\) is a sufficiently differentiable function, and consider the Taylor series expansions of \(y\left (t_{n}+\frac {ih}{2}\right)\), i=0,1,…,5, \(y'\left (t_{n}+3h\right)\) and \(y''\left (t_{n}+3h\right)\). We assume that yn+i=y(tn+ih), fn+3=y′(tn+3h), gn+3=y″(tn+3h) and substitute the coefficients αi/2,i=0,1,…,5, β3 and γ3 into the corresponding formula in (7); after simplifying, we obtain
Following the same process, the local truncation error of each of the additional formulas of HBSDBDF is obtained.
It is therefore established that LTE of the HBSDBDF in (14) have order p=(7,7,7,7,7,7)τ because from the local truncation error for all the methods 7, 9–13 after expansion of each method by Taylor’s series the coefficient C0 to C7=0 and C8≠0 is the error constant given by
where τ is transpose.
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.
Stability
Applying the HBSDBDF to the test problems y′=λy, y″=λ2y, λ<0 yields
where the amplification matrix Q(z) is given by \(Q(z)=\left (W_{1}-z U_{1}-z^{2} D_{1}\right)^{-1}(W_{0})\). The matrix Q(z) has eigenvalues \(\left \{\phi _{1},\phi _{2},\phi _{3},\phi _{4},\phi _{5},\phi _{6}\right \}=\{ 0, 0,0,0,0, \phi _{6}\}\). ϕ6 is the dominant eigenvalue called the stability function with real coefficient given by
The stability domain of the method is defined as \(\delta =\{z \in C: \Re (z) \leq 1\}\). To determine the stability of the HBSDBDF, we state the following definitions:
-
Dahlquist [3]: A numerical method whose stability region contains the entire left half plane is said to be A-stable.
-
Ehle [21]: A method that is A-stable and \(Lim[R(z)] \rightarrow 0\) as \(z \rightarrow -\infty \) 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\rightarrow -\infty } 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 [1]. This was due to the fact that the new method converges with correct digit of 13 from h=0.5 to h=0.1250 as shown in Table 1.
Example 2
Our second example is a nearly sinusoidal problem given in ([10], Example 5.1, p 636) in the range 0≤t≤10
with exact solution
,
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} \frac {|y_{i} - y(t_{i})|}{|(1 + y(t_{i}))|}\) and the ROC are calculated with the formula \(ROC=log_{2}\left (\frac {e^{2 h}}{e^{h}}\right)\), where eh is the greatest absolute error for h shown in Table 2. Clearly, the ROC is consistent with the order of the new block scheme.
Example 3
Consider the highly stiff system, see [19].
The eigenvalues of the Jacobian of the system are approximately λ1=−1.000000000562500×106 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 y1 in blue, y2 in brown, and y3 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 0≤t≤10
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.
Availability of data and materials
Not applicable.
References
Fatunla, S. O.: Block methods for second order IVPs. Intern. J. Comput. Math. 41, 55–63 (1991).
Jackson, L. W., Kenue, S. K.: A fourth order exponentially fitted method. SIAM J. Numer. Anal. 11, 965–978 (1974).
Dahlquist, G.: A special stability problem for linear multistep methods. BIT. 3, 27–43 (1963).
Enright, W. H.: Continuous numerical methods for ODEs with defect control. J. Comput. Appl. Math. 125(1-2), 159–170 (2000).
Onumanyi, P., Awoyemi, D. O., Jator, S. N., Sirisena, U. W.: New linear mutlistep methods with continuous coefficients for first order initial value problems. J. Nig. Math. Soc. 13, 37–51 (1994).
Akinfenwa, O. A., Jator, S. N., Yao, N. M.: On the seven step backwrard differentiation formula with continuous coefficients for stiff systems. Far East J. Math. Sci. 56(1), 23–41 (2011).
Akinfenwa, O. A., Jator, S. N., Okunuga, S. A., Sofoluwe, A. B.: A one step hybrid multistep method with two off-step points for solution of second order ordinary differential equations. Int. J. Comput. Appl. Math. 7(3), 235–247 (2012).
Akinfenwa, O. A., Jator, S. N., Yao, N. M.: A continuous hybrid method for solving parabolic PDEs. AMSE J. Adv. Modeling A Gen. Math. 48, 17–27 (2011).
Akinfenwa, O. A., Jator, S. N., Yao, N. M.: Continuous block backward differentiation formula for solving stiff ordinary differential equations. Comput. Math. Appl. 65(7), 996–1005 (2013).
Akinfenwa, O. A.: Third derivative hybrid block integrator for solution of stiff systems of initial value problems. Afrika Matematika. 28(3-4), 629–641 (2016).
Akinfenwa, O., Jator, S., Yoa, N.: Eight order backward differentiation formula with continuous coeffcients for stiff ordinary differential equations. Int. J. Math. Comput. Sci.l7(4), 172–176 (2011).
Jator, S. N., Agyingi, E.: Block hybrid k-step backward differentiation formulas for large stiff systems. Int. J. Comput. Math.Article ID 162103, 8 (2014). doi:10.1155/2014/162103.
Jator, S. N.: On the hybrid method with three off-step points for initial value problems. Int. J. Math. Educ. Sci. Technol. 41(1), 110–118 (2010).
Jator, S. N., Swindell, S., French, R.: Trigonometrically fitted block Numerov type method for y″=f(t,y,y′). Numer. Algoritm. 62(1), 13–26 (2012).
Chartier, P: L-Stable parallel one-block methods for ordinary differential equations. SIAM J. Numer. Anal. 31(2), 552–571 (1994).
Shampine, L. F., Watts, H. A.: Block implicit one-step methods. Math. Comp. 23, 731–740 (1969).
Chu, M. T., Hamilton, H.: Parallel solution of ODE’s by multi-block methods, SIAM Journal of Science Statistical Computing. SIAM J. Numer. Anal. 8, 342–353 (1987).
Suleiman, M. B., Musa, H., Ismail, F., Senu, N.: A new variable step size block backward differentiation formula for solving stiff initial value problems. Int. J. Comput. Math. 90(11), 2391–2408 (2013). https://doi.org/10.1080/00207160.3013.776677.
Yakubu, D. G., Markus, S.: Second derivative of high-order accuracy methods for the numerical integration of stiff initial value problems. Afr. Mat. 27(5-6), 963–977 (2016). https://doi.org/10.1007/s13370-015-0389-5.
Henrici, P.: Discrete Variable Methods in ODEs. John Wiley, New York (1962).
Ehle, B. L.: On padé approximations to the exponential function and A-stable methods for the solution of initial value problems, Thesis, Univ. of Waterloo, Waterloo, Ontario, Canada (1969).
Gear, C. W.: The automatic integration of stiff ordinary differential equations in information Processing. Informations processing 68 which is the proceedings of the International Federation for Information Processing IFIP. 68, 187–193 (1969).
Cash, J. R.: Second derivative extended backward differentiation formulas for numerical integration of stiff systems. SIAM J. Numer. Anal. 18(1), 21–36 (1981).
Acknowledgements
The authors would like to thank Mr Azeez Faisal Abiola of the Centre for Information Technology And Systems (C.I.T.S) University of Lagos who provides Turnitin platform for similarity checks.
Funding
There is no funding available to the authors on this work.
Author information
Authors and Affiliations
Contributions
Akinfenwa O.A proposed and derived the hybrid second derivative block method, wrote the abstract and introduction, coordinated the research work all through, and wrote the “Development of the method” and “Conclusion” sections. Abdulganiy R. I. did the analysis of the numerical method, crosschecked the computation, and wrote the “Analysis of HBSDBDF” section. Akinnukawe B.I. implemented the numerical method on the examples, verified the analysis, and wrote the “Numerical examples” and “Results and discussion” sections. Okunga S.A. read through the work for correctness and provided the leadership required to ensure the correctness of the mathematical presentation of the work. The authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors state that there are no competing interest.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Akinfenwa, O.A., Abdulganiy, R.I., Akinnukawe, B.I. et al. Seventh order hybrid block method for solution of first order stiff systems of initial value problems. J Egypt Math Soc 28, 34 (2020). https://doi.org/10.1186/s42787-020-00095-3
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s42787-020-00095-3