Study on Banded Implicit Runge–Kutta Methods for Solving Stiff Differential Equations

Liu, M. Y.; Zhang, L.; Zhang, C. F.

doi:https://doi.org/10.1155/2019/4850872

Mathematical Problems in Engineering

On this page

Abstract Introduction Methods Results Discussion and Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 4850872 | https://doi.org/10.1155/2019/4850872

Study on Banded Implicit Runge–Kutta Methods for Solving Stiff Differential Equations

M. Y. Liu,¹L. Zhang,¹and C. F. Zhang¹

Academic Editor: Luigi Rodino

Received13 Jun 2019

Revised12 Aug 2019

Accepted29 Aug 2019

Published10 Oct 2019

Abstract

The implicit Runge–Kutta method with A-stability 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 round-off 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 A-stable 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 low-dimensional systems. In addition, not all SIRK methods are A-stable, 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 (Gauss-IRK) 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 Gauss-IRK 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 Gauss-IRK method can be obtained is , and the SDIRK and SIRK methods have the maximum attainable order . It is clear that the Gauss-IRK 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 Gauss-IRK 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 Gauss-IRK 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 pre-multiplied 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 Gauss-IRK 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 Gauss-IRK method are retained. Therefore, the accuracy of the BIRK method are not reduced compared with the Gauss-IRK 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 Gauss-IRK method, and the stability function of the BIRK method is the same as that of the Gauss-IRK method. Ehle shows that the stability function of the stage order Gauss-IRK method is (s, s)-Padé approximation and the method is A-stable [25]. Therefore, the stage BIRK method of order is A-stable.

4. Experiments and Results

This section describes the implementation of some IRK methods on the MATLAB R2017b in the computer with a CPU i5-4590, 3.3 GHz, and RAM 6 GB. Take stiff systems (27)–(29) as three simulation problems:

Several IRK methods with stage size and A-stable 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 Gauss-IRK 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 Gauss-IRK 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 Gauss-IRK 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 three-dimensional systems (27) and four-dimensional 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 Gauss-IRK 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 Gauss-IRK method, but shorter than that of the SIRK method. This is because the BIRK method also needs a simple transformation. For the three-dimensional system (27) and the four-dimensional 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 Gauss-IRK 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 Gauss-IRK 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 Gauss-IRK 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 Gauss-IRK 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 6-dimensional 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 Gauss-IRK method in the case of step size . When , only one Newton iteration is needed for each step, so the BIRK method and Gauss-IRK 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 Gauss-IRK 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 Gauss-IRK method have a higher accuracy when the step size is fixed.

We know that the order of the BIRK method and the Gauss-IRK method for 2 stage and A-stable 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 Gauss-IRK 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 Gauss-IRK methods, which means that the accuracy of SIRK and SDIRK methods will be much lower than that of BIRK and Gauss-IRK 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 Gauss-IRK 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 Gauss-IRK method. For another, compared with the Gauss-IRK 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 Gauss-IRK method is suitable for lower dimensional systems. Because the computational cost for low-dimensional systems is small, the Gauss-IRK 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 Gauss-IRK 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 Differential-Algebraic Problems, Springer-Verlag, Berlin, China, 1996.
V. Nirmala, V. Parimala, and P. Rajarajeswari, “Application of Runge-Kutta 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 Runge-Kutta 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, “Runge-Kutta 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, Saint-Raphaël, France, 2018.
A. Samsudin, N. Rosli, and A. N. Ariffin, “Stability analysis of explicit and implicit stochastic Runge-Kutta 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. Kazemi-Kamyab, 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. 2-3, 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 Runge-Kutta methods and W-transformation,” 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 Runge-Kutta 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 three-parallel and three-processor SDIRK methods for stiff initial-value 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 Runge-Kutta method for time-dependent reaction-diffusion 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 five-stage singly diagonally implicit Runge-Kutta-Nyström method with reduced phase-lag,” 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 third-order 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 singly-implicit Runge-Kutta methods,” Applied Numerical Mathematics, vol. 24, no. 2-3, pp. 343–350, 1997.
View at: Publisher Site | Google Scholar
K. Burrage, J. C. Butcher, and F. H. Chipman, “An implementation of singly-implicit Runge-Kutta methods,” BIT Numerical Mathematics, vol. 20, no. 3, pp. 326–340, 1980.
View at: Publisher Site | Google Scholar
K. Burrage, “A special family of Runge-Kutta 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 Runge-Kutta 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 Singly-implicit Runge-Kutta 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 Runge-Kutta 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 a-stable 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.

PDF Download Citation

Download other formats

Order printed copies

Views

2423

Downloads

974

Citations