Research Article  Open Access
M. Y. Liu, L. Zhang, C. F. Zhang, "Study on Banded Implicit Runge–Kutta Methods for Solving Stiff Differential Equations", Mathematical Problems in Engineering, vol. 2019, Article ID 4850872, 8 pages, 2019. https://doi.org/10.1155/2019/4850872
Study on Banded Implicit Runge–Kutta Methods for Solving Stiff Differential Equations
Abstract
The implicit Runge–Kutta method with Astability is suitable for solving stiff differential equations. However, the fully implicit Runge–Kutta method is very expensive in solving large system problems. Although some implicit Runge–Kutta methods can reduce the cost of computation, their accuracy and stability are also adversely affected. Therefore, an effective banded implicit Runge–Kutta method with high accuracy and high stability is proposed, which reduces the computation cost by changing the Jacobian matrix from a full coefficient matrix to a banded matrix. Numerical solutions and results of stiff equations obtained by the methods involved are compared, and the results show that the banded implicit Runge–Kutta method is advantageous to solve large stiff problems and conducive to the development of simulation.
1. Introduction
A stiff equation is a differential equation for which certain numerical methods for solving the equation are numerically unstable, unless the step size is taken to be extremely small. Differential equations in chemical reactions, automatic control, electronic networks, and biology are often expressed as stiff equations [1]. For general numerical methods, stiff equations must take small steps to obtain satisfactory results, which increases the computational complexity and is unfavourable to simulation. In the meantime, with the increase of step size, the roundoff error becomes larger and larger, which affects the accuracy of the simulation. Therefore, for the stiff problem, the general numerical method is not competent, and it is necessary to use a method with better stability to solve stiff differential equations.
The Runge–Kutta method is a widely used numerical method because of its high accuracy and high stability [2–6]. The Runge–Kutta method is divided into two types: explicit Runge–Kutta method and implicit Runge–Kutta (IRK) method. The IRK method has wider stable region and higher accuracy compared with explicit Runge–Kutta methods [7]. Therefore, the IRK method is suitable for solving stiff differential equations, and abundant results have been mentioned in the literature [8–11].
However, unlike explicit Runge–Kutta methods, the IRK method needs to solve a system of equations in each step of calculation, which increases the cost of calculation. The aim of reducing the computational cost for IRK methods leads to the singly diagonally implicit Runge–Kutta (SDIRK) method, whose coefficient matrix is lower triangular with same diagonal elements instead of using a full coefficient matrix [12–16]. However, it also has some drawbacks that its stability and accuracy are adversely affected by the simplification of the coefficient matrix . Another classic transformed method whose Runge–Kutta matrix has just one real fold eigenvalue is called a singly implicit Runge–Kutta (SIRK) method [17, 18]. Although this method has obvious computational advantages for large systems over fully IRK method, it has a complex transformation [19]. Apart from transformation, the stability and accuracy are also reduced compared to traditional IRK methods. Therefore, the search for high accuracy and high stability while retaining computational advantages leads to the banded implicit Runge–Kutta (BIRK) method.
2. The Problem of Implementing IRK Methods
For initial value problems of an ordinary differential equation with dimension ,
The approximate solution obtained by the stage implicit Runge–Kutta method is of the formwhere and is the step size of each integration. is the node and is the weight coefficient. The coefficient matrix is in the form of and .
A complicated iteration scheme for the solutions of the internal stages is necessary for IRK methods. In order to solve (2), let
If the coefficient matrix of the IRK method is nonsingular, then (4) can be written as
To solve (5), the Newton–Raphson iteration method can be used. The system of the Newton method can be written aswhere is the Kronecker product, is the identity matrix, and is the identity matrix. Let , and is the matrix.
For the IRK method, the coefficient matrix is required for each iteration.
Let be the modified approximation after a single iteration, be the Newton update, and . Then, and are the matrix because the dimensions of and are .
We can see that the computation of the stage values is typically the most expensive component in their implementation because of the iteration. For each iteration, we need to do calculations about evaluation of , evaluation of , factorization of matrix , and back substitution to the Newton update . In addition, the computation cost grows rapidly as the stage of methods or the dimension of the system increases. Especially, for simulation of large system problems the cost is increasingly expensive. Therefore, it is needed to reduce the cost of calculation by the following method.
3. Methods
3.1. SDIRK and SIRK
In order to reduce the computation cost of the IRK method with a full coefficient matrix , the SDIRK method is introduced. SDIRK is a method whose coefficient matrix is lower triangular with same diagonal elements . For this method, will be changed from (9) to a matrix of this form:
In this case, it may only be necessary to use repeatedly the factorization of to obtain the upper triangular matrix of the factorization of the matrix . So, the cost of factorization of the matrix is decreased greatly, and the calculation cost of back substitution to the Newton update is reduced.
Although the SDIRK method has obvious computational advantages over the fully IRK method, the method has some drawbacks that its stability and accuracy are affected by the simplification of the coefficient matrix A. The maximum attainable order for stage with the Astable IRK method is , whereas the stage SDIRK method has the maximum attainable order only when [20].
Simultaneously, the SIRK method can achieve the same effect of reducing the cost of calculation as the SDIRK method through transformation. Let be a nonsingular matrix such that is the Jordan canonical form of with the same diagonal elements.where
By using some transformation matrices, we can still have the iteration matrix for each stage iteration.
The other fact is that it needs a complex transformation, and it is obvious from (12)–(15). It is required to find a transformation matrix which can make a Jordan canonical form of the coefficient matrix and then transform , , , and into , , , and , respectively. This transformation cost is relatively high, which is even poor for lowdimensional systems. In addition, not all SIRK methods are Astable, and the maximum attainable order is reduced. The SIRK method has the maximum attainable order [21]. The effective method is to reduce the computation cost on the basis of maintaining high accuracy and good stability. Therefore, our BIRK method is proposed as follows.
3.2. Construction of BIRK Methods
The construction of IRK methods relies on the simplifying assumptions
Theorem 1 (see [22]). If the coefficient , , and of a Runge–Kutta method satisfy , , and with and , then the method is of order .
According to Theorem 1, the modified coefficients , , and will affect the order of the method. Because of the modified order conditions which have to be solved, it is difficult to construct the SDIRK method or SIRK method with high stages and high order. Different from the SDIRK method and SIRK method, we do not change the coefficient matrix , node coefficient , and weight coefficient of IRK methods, which will save a lot of work. We just need to find a mature IRK method with good stability and high accuracy and then try to reduce the computation cost of the method. The method we consider is the Gauss implicit Runge–Kutta (GaussIRK) method which was introduced by Butcher [22].
Assume the numerical solution is equal to the exact solution at the point . Then,where is the local truncation error. The accuracy of a numerical solution can be interpreted in terms of the local truncation error. A numerical solution method is said to be of order as in (19), where is an integer and .
The GaussIRK method is a collocation method based on the Gaussian quadrature formulas. Since the Gaussian quadrature formulas for point have the algebraic accuracy of , the truncation error is , and it satisfies the order conditions in Theorem 1, so the order the implicit method can be obtained is . The order of numerical methods is an important indicator to measure the accuracy of the method. Generally, the higher the order, the higher the accuracy of the numerical method. The stage GaussIRK method can be obtained is , and the SDIRK and SIRK methods have the maximum attainable order . It is clear that the GaussIRK method has higher order and higher accuracy, which is the main reason why we choose this method. Even so, when this method is applied to a large system simulation, it is hesitant because it is fully implicit with a large computational cost. This is one disadvantage that contributes to the expensive computation cost by the GaussIRK method with a full coefficient matrix , and the computation cost grows increasingly expensive for high stage methods and high dimensional systems. So, it is necessary to reduce its computational cost. Since the transformation method of SIRK is relatively complex, we proposed a simple method to reduce the computation cost in this paper.
Firstly, the method introduced in this paper retains the coefficients , , and of the GaussIRK method. This avoids the same result as the SIRK method, which changes the order condition and reduces the order of the method due to the change in the method coefficients. Then, we let the two sides of equation (6) be premultiplied by in the process of the method implementation, and this operation has no effect on the equation results. Equation (6) becomeswherewhere is the adjugate matrix of the coefficient matrix , and .
This is similar in form to the transformation method of SIRK, and we can see that a simple transformation method is not as complicated as the transformation of SIRK. can be easily found by , and there is no need to transform vector . Newton iterative update also can be obtained directly from equation (8) without repeated transform. Furthermore, will be a matrix as the following form:
is
The main achievement is that the matrix is a banded diagonal matrix with bandwidth rather than a full coefficient matrix of . So we call it the banded implicit Runge–Kutta (BIRK) method.
Theorem 2 (see [23, 24]). If the order banded matrix with the upper bandwidth and the lower bandwidth has a factorization , then is the lower triangular matrix with lower bandwidth , and is the upper triangular matrix with upper bandwidth .
According to Theorem 2, the matrix is decomposed into an upper triangular matrix with bandwidth and lower triangular matrix with bandwidth by making full use of the characteristic of banded matrix factorization. Therefore, compared with the traditional GaussIRK method, the BIRK method reduces the computational complexity of the factorization and back substitution to the Newton update in iteration. The larger the system dimension , the more effective the method is. In spite of this, the stage BIRK method can obtain the maximum order because the coefficients of the GaussIRK method are retained. Therefore, the accuracy of the BIRK method are not reduced compared with the GaussIRK method. Furthermore, the BIRK method is easier to implement programmatically than SDIRK and SIRK.
3.3. Stability of BIRK Methods
When each IRK method is applied to the test equation,
A straightforward computation
Solving this and inserting equations (2) and (3) leads towhere , , and . is called the stability function of the IRK method.
The equation (26) shows that the stability function of the IRK method is related to the coefficients of the method, and the transformation method in the process of implementation has no effect on its stability function. Because the BIRK method retains the coefficients of the GaussIRK method, and the stability function of the BIRK method is the same as that of the GaussIRK method. Ehle shows that the stability function of the stage order GaussIRK method is (s, s)Padé approximation and the method is Astable [25]. Therefore, the stage BIRK method of order is Astable.
4. Experiments and Results
This section describes the implementation of some IRK methods on the MATLAB R2017b in the computer with a CPU i54590, 3.3 GHz, and RAM 6 GB. Take stiff systems (27)–(29) as three simulation problems:
Several IRK methods with stage size and Astable are used to solve the stiff systems (27)–(29). These three stiff systems have different dimensions. The test results are shown in Tables 1–3.



In the Tables 1–2, the test results of systems (27) and (28) are presented. When the step size , the maximum absolute error of the BIRK method and GaussIRK method is the same and the maximum absolute error of the SDIRK method and SIRK method is the same, but the maximum absolute error of the BIRK method and GaussIRK method is smaller than that of the SDIRK method and SIRK method. When step size , the maximum absolute error of the BIRK method and GIRK method is similar to that of the SDIRK method and SIRK method in the case of . When step size , the maximum absolute error of the BIRK method and GaussIRK method is similar to that of the SDIRK method and SIRK method in the case of . In addition, the SDIRK method has the shortest computation time, although it does not have the least number of Newton iterations among the four IRK methods. This is because the SDIRK method has the smallest computation cost in factorization of the matrix and back substitution to the Newton update. Although the SIRK method has the same computation cost in factorization of the matrix and back substitution to the Newton update, it needs a complex transformation process. It has no advantage for threedimensional systems (27) and fourdimensional systems (28), and the SIRK method has the largest number of Newton iterations, so it always needs the longest computation time. The Newton iteration times of the BIRK method and GaussIRK method differ little or sometimes the same, but they are less than the number of Newton iterations of the SDIRK method and the SIRK method. The calculation time of the BIRK method is slightly longer than that of the GaussIRK method, but shorter than that of the SIRK method. This is because the BIRK method also needs a simple transformation. For the threedimensional system (27) and the fourdimensional system (28), the reduction of computing in the factorization of the matrix and back substitution to the Newton update is not obvious, so the computation time is slightly more than that of the GaussIRK method.
Table 3 shows the test results of system (29). The same situation as Tables 1–2 is that the maximum absolute error of the BIRK method and GaussIRK method is the same and the maximum absolute error of the SDIRK method and SIRK method is the same, and the maximum absolute error of the BIRK method and GaussIRK method is smaller than that of the SDIRK method and SIRK method as the step size . When step size , the maximum absolute error of the BIRK method and GaussIRK method is similar to that of the SDIRK method and SIRK method in the case of . In Table 3, the SDIRK method still has the minimum computation time. Due to the complexity of the transformation process of the SIRK method, there is no obvious computing advantage for the 6dimensional system (29), and the number of Newton iterations is the largest, so the SIRK method also needs the most computing time compared with the other three methods. It is possible that the SIRK method needs to be presented in a system of larger dimensions to show that it has the advantage of reducing the cost of computation. However, the BIRK method has the least number of Newton iterations compared with the other three methods, and the BIRK method also takes less computing time than the GaussIRK method in the case of step size . When , only one Newton iteration is needed for each step, so the BIRK method and GaussIRK method have the same number of iterations. Although the number of Newton iterations is the same, the calculation time of the BIRK method is still shorter than that of the GaussIRK method, which indicates that the BIRK method has the advantage of reducing the cost of computation.
To compare the accuracy of the four IRK methods, we tested the maximum absolute error of the three stiff systems (27)–(29) as the step size increased from 0.001 to 0.01. The results are shown in Figures 1–3.
As can be seen from Figures 1–3, the maximum absolute error of the stiff systems (27)–(29) also gradually increases with the increase in integration step size. When the integral step size is the same, the maximum absolute error of the BIRK method is the same as the maximum absolute error of the GIRK method, and the maximum absolute error of the SDIRK method is the same as the maximum absolute error of the SIRK method, but the maximum absolute error of the BIRK method and the GIRK method is always smaller than the maximum absolute error of the SDIRK method and the SIRK method. Therefore, the BIRK method and the GaussIRK method have a higher accuracy when the step size is fixed.
We know that the order of the BIRK method and the GaussIRK method for 2 stage and Astable method are 4, but the SIRK method and the SDIRK are 3. The accuracy of the numerical method with higher order should be better for the same stage IRK methods. Since the highest order obtainable by SIRK and SDIRK methods is and the highest order obtainable by BIRK and GaussIRK is , if we add the same stages to these numerical methods, the order of SIRK and SDIRK methods will be much lower than that of BIRK and GaussIRK methods, which means that the accuracy of SIRK and SDIRK methods will be much lower than that of BIRK and GaussIRK methods. In the other words, for the same stage SDIRK method and the SIRK method, if they want to achieve the accuracy similar to the BIRK method and the GaussIRK method, they need to greatly reduce the step size, which will increase the calculation cost or may cost more computation time than the BIRK method and the GaussIRK method. For another, compared with the GaussIRK method, the computation cost in the factorization and back substitution to the Newton update in iteration of the BIRK method are less as mentioned in the previous section. Therefore, the BIRK method is more effective than the other three methods, with the advantages of high accuracy, high stability, and low computational cost.
5. Discussion and Conclusions
Although the SDIRK method and the SIRK method have lower accuracy, the cost of factorization and back substitution in the Newton iteration is greatly reduced. So, the SDIRK method and the SIRK method are suitable for the simulation of large systems with low precision requirements when the number of Newton iterations is similar. The traditional GaussIRK method is suitable for lower dimensional systems. Because the computational cost for lowdimensional systems is small, the GaussIRK method has higher accuracy than the SDIRK method and SIRK, and the transformation process of BIRK and SIRK is omitted. However, in some practical systems such as chemical reactions, electronic networks, and especially automatic control systems, systems are often large and complex with many factors. The BIRK method has higher accuracy than SDIRK and SIRK and less cost of computation than the GaussIRK method. Therefore, the BIRK method not only satisfies the requirements of high accuracy and high stability for solving the stiff differential equations in large system simulation but also has lower computation cost. It becomes more effective while the dimension of the system is increased, and to some extent, unnecessary costs are avoided. So, the BIRK method has broad application prospects.
In practical applications, we certainly hope to combine the characteristics of specific problems, make full use of the differences between different numerical methods, choose a more appropriate method, and strive to get the best possible results.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under grant no. 61773008.
References
 E. Hairer, G. Wanner, and S. P. Nørsett, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, SpringerVerlag, Berlin, China, 1996.
 V. Nirmala, V. Parimala, and P. Rajarajeswari, “Application of RungeKutta method for finding multiple numerical solutions to intuitionistic fuzzy differential equations,” Journal of Physics: Conference Series, vol. 1139, no. 1, Article ID 012012, 2018. View at: Publisher Site  Google Scholar
 B. Su, X. Tuo, and L. Xu, “Two modified symplectic partitioned RungeKutta methods for solving the elastic wave equation,” Journal of Geophysics and Engineering, vol. 14, no. 4, pp. 811–821, 2017. View at: Publisher Site  Google Scholar
 I. E. Athanasakis, E. P. Papadopoulou, and Y. G. Saridakis, “RungeKutta and Hermite Collocation for a biological invasion problem modeled by a generalized Fisher equation,” Journal of Physics: Conference Series, vol. 490, no. 1, Article ID 012133, 2014. View at: Publisher Site  Google Scholar
 J. H. Tian, J. Zou, Y. S. Wang, J. Liu, J. Yuan, and Y. Zhou, ““Simulation of bipolar charge transport with trapping and recombination in polymeric insulators using Runge–Kutta discontinuous Galerkin method,” Journal of Physics D: Applied Physics, vol. 41, no. 19, Article ID 195416, 2008. View at: Publisher Site  Google Scholar
 S. Chowdhury and P. K. G. Das, Numerical Solutions of Initial Value Problems Using Mathematica, Morgan & Claypool Publishers, SaintRaphaël, France, 2018.
 A. Samsudin, N. Rosli, and A. N. Ariffin, “Stability analysis of explicit and implicit stochastic RungeKutta methods for stochastic differential equations,” Journal of Physics: Conference Series, vol. 890, no. 1, Article ID 012084, 2017. View at: Publisher Site  Google Scholar
 V. KazemiKamyab, A. H. Van Zuijlen, and H. Bijl, “Analysis and application of high order implicit Runge–Kutta schemes to collocated finite volume discretization of the incompressible Navier–Stokes equations,” Computers & Fluids, vol. 108, no. 15, pp. 107–115, 2015. View at: Publisher Site  Google Scholar
 A. Jameson, “Evaluation of fully implicit Runge Kutta schemes for unsteady flow calculations,” Journal of Scientific Computing, vol. 73, no. 23, pp. 819–852, 2017. View at: Publisher Site  Google Scholar
 K. Liu, X. Liao, and Y. Li, “Parallel simulation of power systems transient stability based on implicit RungeKutta methods and Wtransformation,” Electric Power Components and Systems, vol. 45, no. 20, pp. 2246–2256, 2017. View at: Publisher Site  Google Scholar
 V. Antohe and I. Gladwell, “Performance of Gauss implicit RungeKutta methods on separable Hamiltonian systems,” Computers & Mathematics with Applications, vol. 45, no. 1–3, pp. 481–501, 2003. View at: Publisher Site  Google Scholar
 J. M. Franco, I. Gómez, and L. Rández, “SDIRK methods for stiff ODEs with oscillating solutions,” Journal of Computational and Applied Mathematics, vol. 81, no. 2, pp. 197–209, 1997. View at: Publisher Site  Google Scholar
 J. M. Franco and I. Gomez, “Two threeparallel and threeprocessor SDIRK methods for stiff initialvalue problems,” Journal of Computational and Applied Mathematics, vol. 87, no. 1, pp. 119–134, 1997. View at: Publisher Site  Google Scholar
 W. Liao and Y. Yan, “Singly diagonally implicit RungeKutta method for timedependent reactiondiffusion equation,” Numerical Methods for Partial Differential Equations, vol. 27, no. 6, pp. 1423–1441, 2011. View at: Publisher Site  Google Scholar
 M. K. Wing, N. Senu, and M. Suleiman, “A fivestage singly diagonally implicit RungeKuttaNyström method with reduced phaselag,” AIP Conference Proceedings, vol. 1482, no. 1, pp. 315–320, 2012. View at: Publisher Site  Google Scholar
 O. Y. Ababneh and R. Ahmad, ““Construction of thirdorder diagonal implicit Runge–Kutta methods for stiff problems,” Chinese Physics Letters, vol. 26, no. 8, Article ID 080503, 2009. View at: Publisher Site  Google Scholar
 J. C. Butcher and P. Chartier, “A generalization of singlyimplicit RungeKutta methods,” Applied Numerical Mathematics, vol. 24, no. 23, pp. 343–350, 1997. View at: Publisher Site  Google Scholar
 K. Burrage, J. C. Butcher, and F. H. Chipman, “An implementation of singlyimplicit RungeKutta methods,” BIT Numerical Mathematics, vol. 20, no. 3, pp. 326–340, 1980. View at: Publisher Site  Google Scholar
 K. Burrage, “A special family of RungeKutta methods for solving stiff differential equations,” BIT Numerical Mathematics, vol. 18, no. 1, pp. 22–41, 1978. View at: Publisher Site  Google Scholar
 R. Alexande, “Diagonally implicit RungeKutta methods for stiff O. D. E.’s,” SIAM Journal on Numerical Analysis, vol. 14, no. 6, pp. 1006–1021, 1977. View at: Publisher Site  Google Scholar
 D. J. L. Chen, “The efficiency of Singlyimplicit RungeKutta methods for stiff differential equations,” Numerical Algorithms, vol. 65, no. 3, pp. 533–554, 2014. View at: Publisher Site  Google Scholar
 J. C. Butcher, “Implicit RungeKutta processes,” Mathematics of Computation, vol. 18, no. 85, pp. 50–64, 1964. View at: Publisher Site  Google Scholar
 G. H. Golub and C. F. Van Loan, Matrix Computations, vol. 3, John Hopkins University Press, London, UK, 2012.
 J. W. Demmel, Applied numerical linear algebra, vol. 56, SIAM, Philadelphia, PA, USA, 1997.
 B. L. Ehle, “On padé approximations to the exponential function and astable methods for the numerical solution of initial value problems,” University of Waterloo, Waterloo, Canada, 1969, Doctoral dissertation. View at: Google Scholar
Copyright
Copyright © 2019 M. Y. Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.