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

Authors: M. Y. Liu; L. Zhang; C. F. Zhang
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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 matrixA 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 A. Another classic transformed method whose Runge–Kutta matrix has just one real s 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 dimensionN,(1)y′=fyx,yx0=y0.The approximate solution obtained by thes stage implicit Runge–Kutta method is of the form(2)Yi=yn+h∑j=1saijfxn+cjh,Yj,i=1,…,s,(3)yn+1=yn+h∑j=1sbjfxn+cjh,Yj,where xn=x0+nhn=0,1,…,h>0 and h is the step size of each integration. ci is the node and bi is the weight coefficient. The coefficient matrix A=aij i,j=1,2,…,s is in the form of s×s and ci=∑jsaij.A complicated iteration scheme for the solutions of the internal stages is necessary for IRK methods. In order to solve (2), let(4)zi=Yi−yn.If the coefficient matrix of the IRK method is nonsingular, then (4) can be written as(5)z1⋮zs=Ahftn+c1h,yn+z1⋮hftn+csh,yn+zs.To solve (5), the Newton–Raphson iteration method can be used. The system of the Newton method can be written as(6)MΔZk=−Zk+hA⊗INFZk,(7)M=IS⊗IN−hA⊗J,(8)Zk+1=Zk+ΔZk,where ⊗ is the Kronecker product, Is is the s×s identity matrix, and IN is the N×N identity matrix. Let J=∂f/∂y, and M is the sN×sN matrix.(9)M=IN−ha11J−ha12J⋯−ha1sJ−ha21JIN−ha22J⋯−ha2sJ⋮⋮⋱⋮−has1J−has2J⋯IN−ha22J.For the IRK method, the coefficient matrixM is required for each iteration.LetZk+1=z1k+1,…,zsk+1T be the modified approximation after a single iteration, ΔZk=Δz1k,…,ΔzskT be the Newton update, and FZk=fxn+c1h,yn+z1k,…,fxn+csh,yn+zskT. Then, F and Z are the s×N matrix because the dimensions of Yi and fxn+cjh,Yj are N.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 ofF, evaluation of −Zk+hA⊗INFZk, LU factorization of matrix M, and back substitution to the Newton update ΔZk. In addition, the computation cost grows rapidly as the stage of methods s or the dimension of the system N 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 matrixA, the SDIRK method is introduced. SDIRK is a method whose coefficient matrix A is lower triangular with same diagonal elements λ. For this method, M will be changed from (9) to a matrix of this form:(10)MD=IN−hλJ0⋯0−ha12ININ−hλJ⋯0⋮⋮⋱⋮−has1INas2IN⋯IN−hλJ.In this case, it may only be necessary to use repeatedly theLU factorization of IN−hλJ to obtain the upper triangular matrix of the LU factorization of the matrix MD. So, the cost of LU factorization of the matrix MD 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 matrixA. The maximum attainable order for s stage with the A-stable IRK method is 2s, whereas the s stage SDIRK method has the maximum attainable order s+1 only when s=1,2,3,5 [20].Simultaneously, the SIRK method can achieve the same effect of reducing the cost of calculation as the SDIRK method through transformation. LetT be a nonsingular matrix such that A¯=T−1AT is the Jordan canonical form of A with the same diagonal elements.(11)M¯ΔZ¯k=−Z¯k+hA¯⊗INF¯Zk,where(12)M¯=IS⊗IN−hA¯⊗J,(13)Z¯k+1=T−1⊗INZk,(14)F¯Zk=T−1⊗INFZk,(15)ΔZk=T⊗INΔZ¯k.By using some transformation matrices, we can still have the iteration matrixIN−hλJ 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 T which can make a Jordan canonical T−1AT form of the coefficient matrix A and then transform M, Zk, ΔZk, and FZk into M¯, Z¯k+1, ΔZ¯k, and F¯Zk, 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 s+1 [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(16)Bp:∑i=1sbiciq−1=1q,q=1,…,p;(17)Cη:∑j=1saijcjq−1=ciqq,i=1,…,s,q=1,…,η;(18)Dζ:∑i=1sbiciq−1aij=bjq1−cjq,j=1,…,s,q=1,…,ζ.Theorem 1 (see [22]).
If the coefficientbi, ci, and aij of a Runge–Kutta method satisfy Bp, Cη, and Dζ with p≤η+ζ+1 and p≤2η+2, then the method is of order p.According to Theorem 1, the modified coefficientsbi, ci, and aij 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 A, node coefficient ci, and weight coefficient bi 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 solutionyn is equal to the exact solution yxn at the point xn. Then,(19)ln+1=yxn+1−yn+1=οhp+1,where ln+1 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 p as in (19), where p is an integer and p>0.The Gauss-IRK method is a collocation method based on the Gaussian quadrature formulas. Since the Gaussian quadrature formulas for points have the algebraic accuracy of 2s−1, the truncation error is ln+1=yxn+1−yn+1=οh2s+1, and it satisfies the order conditions in Theorem 1, so the order the implicit method can be obtained is 2s. 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 s stage Gauss-IRK method can be obtained is 2s, and the SDIRK and SIRK methods have the maximum attainable order s+1. 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 A, 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 coefficientsbi, ci, and aij 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 A∗⊗IN in the process of the method implementation, and this operation has no effect on the equation results. Equation (6) becomes(20)MBΔZk=−A∗⊗INZk+hA∗A⊗INFZk,where(21)MB=A∗⊗IN−hA∗A⊗J,where A∗=Aij i,j=1,2,…,s is the adjugate matrix of the coefficient matrix A, and A∗A=AIS.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.A∗ can be easily found by A, and there is no need to transform vector FZk. Newton iterative update ΔZk also can be obtained directly from equation (8) without repeated transform. Furthermore, MB will be a matrix as the following form:(22)A11IN−hAJA21IN⋯AS1INA12INA22IN−hAJ⋯AS2IN⋮⋮⋱⋮A1SINA2SIN⋯ASSIN−hAJ.A
i
j
I
N is(23)Aij0⋯00Aij⋯0⋮⋮⋱⋮00⋯Aij.The main achievement is that the matrixMB is a banded diagonal matrix with bandwidth Ns−1 rather than a full coefficient matrix of sN×sN. So we call it the banded implicit Runge–Kutta (BIRK) method.Theorem 2 (see [23, 24]).
If then order banded matrix M with the upper bandwidth q and the lower bandwidth p has a LU factorization M=LU, then L is the lower triangular matrix with lower bandwidth p, and U is the upper triangular matrix with upper bandwidth q.According to Theorem 2, the matrixMB is decomposed into an upper triangular matrix with bandwidth Ns−1 and lower triangular matrix with bandwidth Ns−1 by making full use of the characteristic of banded matrix LU factorization. Therefore, compared with the traditional Gauss-IRK method, the BIRK method reduces the computational complexity of the LU factorization and back substitution to the Newton update in iteration. The larger the system dimension N, the more effective the method is. In spite of this, the s stage BIRK method can obtain the maximum order 2s 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,(24)y′=λy.A straightforward computation(25)y1=Rhλy0.Solving this and inserting equations (2) and (3) leads to(26)Rz=1+zbTI−zA−1e,where z=hλ, bT=b1,…,bs, and e=1,…,1T. Rz 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 s stage 2s order Gauss-IRK method is (s, s)-Padé approximation and the method is A-stable [25]. Therefore, the s stage BIRK method of order 2s is A-stable.


## 3.1. SDIRK and SIRK

In order to reduce the computation cost of the IRK method with a full coefficient matrixA, the SDIRK method is introduced. SDIRK is a method whose coefficient matrix A is lower triangular with same diagonal elements λ. For this method, M will be changed from (9) to a matrix of this form:(10)MD=IN−hλJ0⋯0−ha12ININ−hλJ⋯0⋮⋮⋱⋮−has1INas2IN⋯IN−hλJ.In this case, it may only be necessary to use repeatedly theLU factorization of IN−hλJ to obtain the upper triangular matrix of the LU factorization of the matrix MD. So, the cost of LU factorization of the matrix MD 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 matrixA. The maximum attainable order for s stage with the A-stable IRK method is 2s, whereas the s stage SDIRK method has the maximum attainable order s+1 only when s=1,2,3,5 [20].Simultaneously, the SIRK method can achieve the same effect of reducing the cost of calculation as the SDIRK method through transformation. LetT be a nonsingular matrix such that A¯=T−1AT is the Jordan canonical form of A with the same diagonal elements.(11)M¯ΔZ¯k=−Z¯k+hA¯⊗INF¯Zk,where(12)M¯=IS⊗IN−hA¯⊗J,(13)Z¯k+1=T−1⊗INZk,(14)F¯Zk=T−1⊗INFZk,(15)ΔZk=T⊗INΔZ¯k.By using some transformation matrices, we can still have the iteration matrixIN−hλJ 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 T which can make a Jordan canonical T−1AT form of the coefficient matrix A and then transform M, Zk, ΔZk, and FZk into M¯, Z¯k+1, ΔZ¯k, and F¯Zk, 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 s+1 [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(16)Bp:∑i=1sbiciq−1=1q,q=1,…,p;(17)Cη:∑j=1saijcjq−1=ciqq,i=1,…,s,q=1,…,η;(18)Dζ:∑i=1sbiciq−1aij=bjq1−cjq,j=1,…,s,q=1,…,ζ.Theorem 1 (see [22]).
If the coefficientbi, ci, and aij of a Runge–Kutta method satisfy Bp, Cη, and Dζ with p≤η+ζ+1 and p≤2η+2, then the method is of order p.According to Theorem 1, the modified coefficientsbi, ci, and aij 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 A, node coefficient ci, and weight coefficient bi 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 solutionyn is equal to the exact solution yxn at the point xn. Then,(19)ln+1=yxn+1−yn+1=οhp+1,where ln+1 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 p as in (19), where p is an integer and p>0.The Gauss-IRK method is a collocation method based on the Gaussian quadrature formulas. Since the Gaussian quadrature formulas for points have the algebraic accuracy of 2s−1, the truncation error is ln+1=yxn+1−yn+1=οh2s+1, and it satisfies the order conditions in Theorem 1, so the order the implicit method can be obtained is 2s. 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 s stage Gauss-IRK method can be obtained is 2s, and the SDIRK and SIRK methods have the maximum attainable order s+1. 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 A, 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 coefficientsbi, ci, and aij 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 A∗⊗IN in the process of the method implementation, and this operation has no effect on the equation results. Equation (6) becomes(20)MBΔZk=−A∗⊗INZk+hA∗A⊗INFZk,where(21)MB=A∗⊗IN−hA∗A⊗J,where A∗=Aij i,j=1,2,…,s is the adjugate matrix of the coefficient matrix A, and A∗A=AIS.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.A∗ can be easily found by A, and there is no need to transform vector FZk. Newton iterative update ΔZk also can be obtained directly from equation (8) without repeated transform. Furthermore, MB will be a matrix as the following form:(22)A11IN−hAJA21IN⋯AS1INA12INA22IN−hAJ⋯AS2IN⋮⋮⋱⋮A1SINA2SIN⋯ASSIN−hAJ.A
i
j
I
N is(23)Aij0⋯00Aij⋯0⋮⋮⋱⋮00⋯Aij.The main achievement is that the matrixMB is a banded diagonal matrix with bandwidth Ns−1 rather than a full coefficient matrix of sN×sN. So we call it the banded implicit Runge–Kutta (BIRK) method.Theorem 2 (see [23, 24]).
If then order banded matrix M with the upper bandwidth q and the lower bandwidth p has a LU factorization M=LU, then L is the lower triangular matrix with lower bandwidth p, and U is the upper triangular matrix with upper bandwidth q.According to Theorem 2, the matrixMB is decomposed into an upper triangular matrix with bandwidth Ns−1 and lower triangular matrix with bandwidth Ns−1 by making full use of the characteristic of banded matrix LU factorization. Therefore, compared with the traditional Gauss-IRK method, the BIRK method reduces the computational complexity of the LU factorization and back substitution to the Newton update in iteration. The larger the system dimension N, the more effective the method is. In spite of this, the s stage BIRK method can obtain the maximum order 2s 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,(24)y′=λy.A straightforward computation(25)y1=Rhλy0.Solving this and inserting equations (2) and (3) leads to(26)Rz=1+zbTI−zA−1e,where z=hλ, bT=b1,…,bs, and e=1,…,1T. Rz 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 s stage 2s order Gauss-IRK method is (s, s)-Padé approximation and the method is A-stable [25]. Therefore, the s stage BIRK method of order 2s 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:(27)y1′=−0.1y1+49.9y2,y10=2,y2′=−40y2,y20=1,y3′=70y2−300y3,y30=2,t∈0,5,(28)y1′=−y1+y2,y10=1,y2′=−100y1−y2,y20=0,y3′=−100y3+y4,y30=1,y4′=−10000y3−100y4,y40=0,t∈0,5,(29)y1′=−10y1+100y2,y10=1,y2′=−100y1−10y2,y20=1,y3′=−4y3,y30=1,y4′=−y4,y40=1,y5′=−0.5y5,y50=1,y6′=−0.1y6,y60=1,t∈0,5.Several IRK methods with stage sizes=2 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.Table 1
The test results for system (27).
Maximum absolute error
Number of Newton iterations
Time
System (24)
Step sizeh=0.01
BIRK
4.696875980934exp-02
999
0.0818
SDIRK
2.951335592737exp-01
1002
0.0771
SIRK
2.951335592737exp-01
1098
0.1158
Gauss-IRK
4.696875980934exp-02
999
0.0801
Step sizeh=0.005
BIRK
4.632141120340exp-03
1366
0.1015
SDIRK
8.059127347433exp-02
1996
0.0941
SIRK
8.059127347433exp-02
2100
0.1867
Gauss-IRK
4.632141120340exp-03
1364
0.0949
Step sizeh=0.001
BIRK
7.163047590741exp-06
5120
0.2795
SDIRK
1.163100888497exp-03
5139
0.2143
SIRK
1.163100888497exp-03
5276
0.4993
Gauss-IRK
7.163047590741exp-06
5120
0.2769Table 2
The test results for system (28).
Maximum absolute error
Number of Newton iterations
Time
System (25)
Step sizeh=0.01
BIRK
1.899433882193exp-01
2499
0.1194
SDIRK
5.402682507296
2933
0.0929
SIRK
5.402682507296
4076
0.2534
Gauss-IRK
1.899433882193exp-01
2493
0.1068
Step sizeh=0.005
BIRK
1.231643373468exp-02
4254
0.1899
SDIRK
9.714581002096exp-01
4666
0.1325
SIRK
9.714581002096exp-01
5188
0.3267
Gauss-IRK
1.231643373468exp-02
4174
0.1725
Step sizeh=0.001
BIRK
2.003231104197exp-05
15604
0.3766
SDIRK
1.124295167250exp-02
16546
0.3154
SIRK
1.124295167250exp-02
18473
1.0362
Gauss-IRK
2.003231104197exp-05
15526
0.3651Table 3
The test results for system (29).
Maximum absolute error
Number of Newton iterations
Time
System (26)
Step sizeh=0.01
BIRK
6.745545867356exp-03
509
0.0700
SDIRK
2.431772336420exp-01
704
0.0692
SIRK
2.431772336420exp-01
837
0.1206
Gauss-IRK
6.745545867356exp-03
702
0.0829
Step sizeh=0.005
BIRK
4.451608892082exp-04
1003
0.0873
SDIRK
4.724422976707exp-02
1004
0.0834
SIRK
4.724422976707exp-02
1430
0.1654
Gauss-IRK
4.451608892082exp-04
1005
0.0990
Step sizeh=0.001
BIRK
7.399547229947exp-07
5000
0.2770
SDIRK
4.679447423339exp-04
5000
0.2294
SIRK
4.679447423341exp-04
5009
0.5155
Gauss-IRK
7.399547229947exp-07
5000
0.3121In the Tables1–2, the test results of systems (27) and (28) are presented. When the step size h=0.01,0.005,0.001, 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 h=0.01, 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 h=0.005. When step size h=0.005, 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 h=0.001. 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 LU factorization of the matrix and back substitution to the Newton update. Although the SIRK method has the same computation cost in LU 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 LU 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.Table3 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 h=0.01,0.005,0.001. When step size h=0.005, 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 h=0.001. 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 h=0.01,0.005. When h=0.001, 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.Figure 1
Relation curve of maximum absolute error of stiff system (27) with step size. Step size is from 0.001 to 0.1.Figure 2
Relation curve of maximum absolute error of stiff system (28) with step size. Step size is from 0.001 to 0.1.Figure 3
Relation curve of maximum absolute error of stiff system (29) with step size. Step size is from 0.001 to 0.1.As can be seen from Figures1–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 isp=s+1 and the highest order obtainable by BIRK and Gauss-IRK is p=2s, 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 LU 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 ofLU 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.



---

*Source: 4850872-2019-10-10.xml*