Mathematical Problems in Engineering

Volume 2019, Article ID 4850872, 8 pages

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

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

School of Physics & Electronics, Henan University, 475004 Kaifeng, China

Correspondence should be addressed to L. Zhang; nc.ude.uneh.piv@ielgnahz

Received 13 June 2019; Revised 12 August 2019; Accepted 29 August 2019; Published 10 October 2019

Academic Editor: Luigi Rodino

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.

#### 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.