Higher order Runge–Kutta methods

For real-life problems first and second order methods are not applicable, be-cause within reasonable computational time we cannot achieve the numerical solution with some prescribed and required (by the application) accuracy. Our aim is to construct higher order accurate methods on the basis of formula (2.112).

We re-consider the numerical method written in the form (2.116)-(2.117).

As we have shown, this method is exact at most in second order. (The method has second order accuracy when the condition (2.115) is satisfied for its param-eters.) Therefore, to get higher order accurate numerical methods, we have to introduce new parameters and give some conditions for their choice. Based on the formulas (2.116)-(2.117), the following generalization seems to be natural:

k₁ =f(t_i, y_i),

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

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

(2.126)

and the approximation on the new mesh-point is defined as

y_i+1 =y_i+h(σ₁k₁+σ₂k₂+σ₃k₃). (2.127) The parameters of this method, according to the table (2.113), can be written as follows

Our aim is to define the parameters in (2.126) in such a way that the cor-responding numerical method is of third order accuracy. To get this condition, we have to define again the local approximation error, as it was done before.

After some long (but not difficult) calculation we obtain the following result.

Theorem 2.4.5. The numerical method (2.126) has third order accuracy if and only if the conditions

Clearly, (2.129) yields six equations (conditions) for eight unknown values.

From the possible solution we give two cases.

• The method with the parameters 0 1/3 1/3 2/3 0 2/3

1/4 0 3/4

(2.130)

is very popular in the different applications.

• The third order method 0 1/2 1/2

1 −1 2

1/6 2/3 1/6

(2.131)

is also often used. We note that this method has accuracy O(h⁵) for problems with f(t, u) = f(t) with the Simpson formula. Therefore, this method is recommended also for problems where the partial derivative

∂₂f is close to zero.

When we want a higher than third order method (p > 3), then we need some more generalization. This can be done in the following form, which is the natural generalization of the methods (2.116)-(2.117) and (2.126).

Letm≥1 be some given integer. We define the following, so called m-stage explicit Runge–Kutta method:

k₁ =f(t_i, y_i),

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

k3 =f(ti +a3h, yi+hb31k1+hb32k2), . . . .

k_m =f(t_i +a_mh, y_i+hb_m1k₁+hb_m2k₂+. . .+hbm,m−1km−1)

(2.132)

y_i+1 =y_i+h(σ₁k₁+σ₂k₂+. . .+σ_mk_m). (2.133)

By specifying the values of the parameters in this formula, we define the con-crete numerical method. As before, we can write the parameters in a table:

0 a₂ b₂₁ a3 b31 b32

. . . .

a_m b_m1 b_m2 . . . bm,m−1

σ₁ σ₂ . . . σ_m

(2.134)

To write the method in a compact form, we introduce some notations. Let the symbols σ,a ∈ R^m denote the row-vectors with coordinates σ_i and a_i, respectively, (where we always assume that a₁ = 0). Let B ∈ R^m×m denote the matrix with the elements bij, which is a strictly lower triangular matrix, i.e.,

B_ij =

b_ij, for i > j, 0, for i≤j.

Definition 2.4.6. For some explicit Runge–Kutta method the corresponding table of the parameters in the form

a^> B

σ (2.135)

is called Butcher tableau.

In the sequel, when we specify some explicit Runge–Kutta method, then we list the lower triangular part of the matrix B only.

For some fixed explicit Runge–Kutta method the order of the consistency can be defined by some simple, however usually cumbersome computation.

This can be done in the following steps.

• First on the right side of the formulas (2.132)-(2.133) we replace y_i by the values u(t_i).

• Then on the left side of this formula we replace yi+1 byu(ti+h).

• We compute the difference between the two sides, and we define its order in h.

Executing these steps for an explicit Runge–Kutta method, we get the condi-tion by which the given explicit Runge–Kutta method is of p-th order.

Remark 2.4.5. Before giving the condition of consistency for a general ex-plicit Runge–Kutta method, we note the following. In the condition of the second order (2.115) and the third order (2.129), we can observe that any i-th row-sum of the matrix B is equal to the value a_i. Because this is true for any consistent method, this implies that in the Butcher tableau only the matrix B and the vector σ can be chosen arbitrarily, and the vector a is defined from the relation a=Be, where e= [1,1, . . . ,1]∈R^m.

Let us write the condition of consistency of explicit Runge–Kutta methods in a compact way. We will use the notations

aⁿ = [aⁿ₁, aⁿ₂, . . . , aⁿ_m,]∈R^m, A=diag[a₁, a₂, . . . , a_m]∈R^m×m.

Then the condition ofp-th order consistency of an explicit Runge–Kutta method forp= 1,2,3,4 imposes the following requirements for the elements of the ma-trix B and the vector σ:¹⁵

Hence, the following statement is true.

Theorem 2.4.7. The explicit Runge–Kutta method defined by the Butcher tableau (2.135) is consistent if and only if the conditions

Be=a; σ·e^>= 1 (2.137)

are satisfied, i.e., it is true under the conditions

m

The most popular explicit Runge–Kutta method is the method given in Table 2.5, which is a fourth order method.

The algorithm of this method (i.e., how y_i+1, the approximation at the mesh-point ti+1 can be defined from the known approximation yi at the mesh-point t_i) is the following:

15These conditions can be obtained by some cumbersome computations see e.g. in [4].

0

1/2 1/2

1/2 0 1/2

1 0 0 1

1/6 1/3 1/3 1/6

Table 2.5: Butcher tableau of the fourth order explicit Runge–Kutta method.

• By the formulas

k₁ =f(t_i, y_i)

k₂ =f(t_i+ 0.5h, y_i+ 0.5hk₁) k₃ =f(t_i+ 0.5h, y_i+ 0.5hk₂) k4 =f(ti+h, yi+hk3)

(2.139)

we define the values k₁, k₂, k₃ and k₄.

• By the formula

yi+1 =yi+h

6(k1+ 2k2+ 2k3+k4) (2.140) we define the new approximation.

Remark 2.4.6. What is the connection between the number of the stages (m) and the order of consistency (p) for an explicit Runge–Kutta method? As we have seen, for the one-stage explicit Euler method the method is of first order, for the two-stage methods (Heun method, modified Euler method) the methods are of second order, for the three-stage method (2.130)-(2.131) it has third order, and for the four-stage method (2.139)-(2.140) it is consistent in fourth order. This means that for m = 1,2,3,4 the stage and the order of consistency are equal. However, for the cases m ≥5 this is not true anymore:

the order of the method is less than the number of stages, i.e., p < m. The relation between them (up to m = 10) is the following:

m 1,2,3,4 5,6,7 8,9,10

p(m) m m−1 m−2 (2.141)

We enclose an interactive animation, called onestep exp.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.1.

Alternatively, this program suggests three initial-value problems as test problems. The chosen numerical methods can be the explicit Euler method,

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

the trapezoidal method, 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 increasingnthe order of the convergence is also shown.

Some further details of explicit Runge–Kutta methods can be found under the link

http://math.fullerton.edu/mathews/n2003/RungeKuttaMod.html Here one can also see the behavior of the above investigated numerical methods as applied to the differential equation u⁰ = 1−t√³

u.

In document Numerical Methods for Ordinary Differential Equations (Pldal 45-50)