Implicit Runge–Kutta method

As we have seen, the explicit Runge–Kutta methods are defined by the special (strictly lower triangular) matrixB ∈R^m×m and the vectorσ ∈R^m. A further and quite natural generalization is obtained by neglecting the requirement for the matrix B. Namely, we introduce the following notion.

Definition 2.4.8. LetB∈R^m×m be an arbitrary matrix, σ ∈R^m some given vector, and a=Be. The numerical method defined by these parameters, i.e.,

the method

k₁ =f(t_i+a₁h, y_i+h[b₁₁k₁+b₁₂k₂+· · ·+b_1mk_m]), k₂ =f(t_i+a₂h, y_i+h[b₂₁k₁+b₂₂k₂+· · ·+b_2mk_m]), . . . . km =f(ti+amh, yi+h[bm1k1+bm2k2+· · ·+bmmkm]),

(2.142)

y_i+1 =y_i+h(σ₁k₁+σ₂k₂+. . .+σ_mk_m), (2.143) is called Runge–Kutta method. When B is not strictly lower triangular, i.e., there exists an elementb_ij such thati≥j andb_ij 6= 0, then the method (2.142)-(2.143) is called implicit Runge–Kutta method. When Bis not a strictly lower triangular matrix, but it is lower triangular, i.e., b_ij = 0 for any i > j, then the method (2.142)- (2.143) is called diagonally implicit Runge–Kutta method.

Remark 2.4.7. For the diagonally implicit Runge–Kutta method the algo-rithm is the following:

k₁ =f(t_i+a₁h, y_i+hb₁₁k₁),

k₂ =f(t_i+a₂h, y_i+h[b₂₁k₁+b₂₂k₂]),

. . . . k_m =f(t_i+a_mh, y_i+h[b_m1k₁+b_m2k₂+· · ·+b_mmk_m]),

(2.144)

yi+1 =yi+h(σ1k1+σ2k2+. . .+σmkm), (2.145) Remark 2.4.8. The main difference between them is as follows. For the explicit Runge–Kutta method we can compute the intermediate values k_i di-rectly, by substituting the previously defined values k₁, k₂, . . . , ki−1. For the diagonally implicit Runge–Kutta method we cannot compute k_i directly, be-cause, according to (2.144), k_i is present on the left side, too. Therefore, for its computation we have to solve an algebraic equation. Typically, for the implicit Runge–Kutta method the situation is the same, however, the computation of k_i becomes more difficult, because it is contained not only in the formula for the definition of k_i but also in all equations of (2.142). This means that for the definition of the values k_i we have to solve a system of nonlinear algebraic equations formunknown values. Hence the implicit methods (DIRK and IRK) require much more computational work. However, these methods have good stability properties, which make them very beneficial for solving some spe-cial problems (like stiff ones). The other benefit of the implicit Runge–Kutta methods is that usually they are more accurate than the explicit ones.

Remark 2.4.9. The implicit Runge–Kutta method can be written as

For the diagonally implicit Runge–Kutta method the only difference is that the formula (2.146) is changed into the formula

k_j =f(t_i+a_jh, y_i+h

j

X

s=1

b_jsk_s), j = 1,2, . . . m. (2.148) In the sequel we introduce some important implicit Runge–Kutta methods (including also the special case of diagonally implicit Runge–Kutta methods) by giving their Butcher tableau. We will see that to write these methods in the usual form of one-step methods (i.e., to give the function Φ) is more difficult than for the explicit methods.

• Let the Butcher tableau of a method be given by 1 1

1 (2.149)

This is a diagonally implicit Runge–Kutta method, which yields the fol-lowing algorithm:

k₁ =f(t_i+h, y_i+hk₁)

y_i+1 =y_i+hk₁. (2.150)

The realization of this method means the following. In the first step we solve the equation for the unknown k₁¹⁶, then we substitute this value into the second formula. Let us find the function Φ for this method. From the second formula we havehk₁ =y_i+1−y_i. Therefore, by its substitution into the first equation k₁ = f(t_i+h, y_i + (y_i+1 −y_i)) = f(t_i+h, y_i+1).

From the second equation we have y_i+1 =y_i+hf(t_i+h, y_i+1), therefore, Φ(h, t_i, y_i, y_i+1) =f(t_i+h, y_i+1). Hence, the method with (2.149) is the well-known implicit Euler method. As we already know, this is a first order method.

16For the solution usually we use some iterative method (as the Newton method).

• Let the numerical method be given by 0.5 0.5

1 (2.151)

This is also a diagonally implicit Runge–Kutta method, which results in the following computational algorithm:

k₁ =f(t_i+ 0.5h, y_i+ 0.5hk₁)

y_i+1 =y_i+hk₁. (2.152)

Let us find the function Φ for this method, too. From the second equation we have hk₁ = y_i+1−y_i. By its substitution into the first equation, we obtain k₁ =f(t_i+ 0.5h, y_i+ 0.5(y_i+1−y_i)) =f(t_i+ 0.5h,0.5(y_i+y_i+1)).

Hence Φ(h, ti, yi, yi+1) =f(ti+ 0.5h,0.5(yi+yi+1)) and the corresponding one-step method (in the usual compact form) is

yi+1 =yi+hf(ti+ 0.5h,0.5(yi+yi+1)). (2.153) The diagonally implicit Runge–Kutta method (2.153) is called implicit midpoint rule. The order of the method can be defined by analyzing the expression (local approximation error) of the form

g_i =u(t_i+1)−u(t_i)−hf(t_i + 0.5h,0.5(u(t_i) +u(t_i+1))).

Let us define the order in h of g_i. By developing into Taylor series, we get

u(t_i+1)−u(t_i) =hu⁰(t_i) + h²

2 u⁰⁰(t_i) +O(h³), and similarly

f(ti+ 0.5h,0.5(u(ti) +u(ti+1)))

=f(t_i+ 0.5h, u(t_i) + 0.5hu⁰(t_i) +O(h²))

=f(t_i, u(t_i)) + 0.5h∂₁f(t_i, u(t_i)) + 0.5hu⁰(t_i)∂₂f(t_i, u(t_i)) +O(h²)

=f(t_i, u(t_i)) + h

2[∂₁f(t_i, u(t_i)) +f(t_i, u(t_i))∂₂f(t_i, u(t_i))] +O(h²).

When we substitute these value into the formula for g_i, then taking into account the second relation in (2.4) we get gi =O(h³). This means that the implicit midpoint rule is a second order method.

• We consider the method with the Butcher tableau

0 0 0

1 0.5 0.5 0.5 0.5

(2.154)

This is a two-stage implicit Runge–Kutta method with the following this expression and the first relation into the second formula, we get k₂ =f(t_i+h, y_i+ (y_i+1−y_i)) = f(t_i+h, y_i+1). Now, putting this value of k₂, and k₁ from the first equation into the third equation, we obtain Φ(h, ti, yi, yi+1) = 0.5[f(ti, yi) +f(ti +h, yi+1)]. Hence, the numerical method with the Butcher tableau (2.154) is the well-known second order trapezoidal method.

For the following two-stage methods we give only their Butcher tableau.

• We define the method, when in the numerical integration we choose the Gaussian basic points. This implies the following method:

1

The method (2.156) is a two-stage, fourth order implicit Runge–Kutta method.

• We define another two-stage implicit Runge–Kutta method as follows

1

The method (2.157) is called two-stage Radau method. We can show that the order of accuracy of this implicit Runge–Kutta method is 3.

During the consideration of the accuracy of the above implicit Runge–Kutta methods we have seen that the accuracy of the methods (p) can be greater then the number of the stages (m). For instance, the accuracy of the one-stage trapezoidal method has second order accuracy, the two-one-stage (2.156) method with Gaussian basic points has fourth order, etc. This means that the case p > m is also possible for the implicit methods. (The reason is that

the implicit Runge–Kutta methods have more free parameters.) We can show that the upper bound for the accuracy of some m-stage implicit Runge–Kutta method is p≤2m.

The methods with p = 2m, are called implicit Runge–Kutta method with maximum accuracy. Hence, the implicit Euler method, the implicit midpoint rule, and the method (2.156) with Gaussian points are all methods of maximum accuracy.

We enclose an interactive animation, called onestep imp.exe. The program can be downloaded from

http://www.cs.elte.hu/~faragois/nummod_jegyzet_prog/.

The image of this program on the screen can be seen in Figure 2.2.

Figure 2.2: The image on the screen of the interactive program for several implicit one-step methods

Alternatively, this program suggests three initial-value problems as test problems. The chosen numerical methods can be the theta-method, where the value of theta can be chosen, the 3rd order Runge–Kutta method and the 4th order Runge–Kutta method. We can select the discretization step size by giving the parameter n, which is the number of the partitions of the interval.

Pushing ”Calculate” we get the result. The result is given graphically, and the error (in maximum norm) is also indicated. By increasing n the order of the convergence is also shown.