6.2 Runge–Kutta Methods

INTRODUCTION

Probably one of the more popular, as well as most accurate, numerical procedures used in obtaining approximate solutions to a first-order initial-value problem y′ = f(x, y), y(x0) = y0 is the fourth-order Runge–Kutta method devised by the German applied mathematicians Carl David Runge (1856–1927) and Martin Wilhelm Kutta (1867–1944) in the late 1890s. As the words fourth-order suggest, there are Runge–Kutta methods of different orders.

Runge–Kutta Methods

Fundamentally, all Runge–Kutta methods are generalizations of the basic Euler formula (1) of Section 6.1 in that the slope function f is replaced by a weighted average of slopes over the interval defined by xn ≤ x ≤ xn + 1. That is,

(1)

Here the weights wi, i = 1, 2, . . . , m are constants that generally satisfy w1 + w2 + ⋯ + wm = 1 and each ki, i = 1, 2, . . . , m is the function f evaluated at a selected point (x, y) for which xn ≤ x ≤ xn + 1. We shall see that the ki are defined recursively. The number m is called the order of the method. Observe that by taking m = 1, w1 = 1, and k1 = f(xn, yn) we get the familiar Euler formula yn + 1 = yn + hf(xn, yn). Hence Euler’s method is said to be a first-order Runge–Kutta method.

The average in (1) is not formed willy-nilly, but parameters are chosen so that (1) agrees with a Taylor polynomial of degree m. As we have seen in the last section, if a function y(x) possesses k + 1 derivatives that are continuous on an open interval containing a and x, then we can write

where c is some number between a and x. If we replace a by xn and x by xn + 1 = xn + h, then the foregoing formula becomes

where c is now some number between xn and xn + 1. When y(x) is a solution of y′ = f(x, y), in the case k = 1 and the remainder h2y″(c) is small, we see that a Taylor polynomial y(xn + 1) = y(xn) + hy′(xn) of degree 1 agrees with the approximation formula of Euler’s method

A Second-Order Runge–Kutta Method

To further illustrate (1), we consider now a second-order Runge–Kutta method. This consists of finding constants, or parameters, w1, w2, α, and β so that the formula

where (2)

agrees with a Taylor polynomial of degree 2. For our purposes it suffices to say that this can be done whenever the constants satisfy

(3)

This is an algebraic system of three equations in four unknowns and has infinitely many solutions:

(4)

where w2 ≠ 0. For example, the choice w2 = yields w1 = , α = 1, β = 1 and so (2) becomes

where

Since xn + h = xn+1 and yn + hk1 = yn + hf(xn, yn) the foregoing result is recognized to be the improved Euler’s method that is summarized in (3) and (4) of Section 6.1.

In view of the fact that w2 ≠ 0 can be chosen arbitrarily in (4), there are many possible second-order Runge–Kutta methods. See Problem 2 in Exercises 6.2.

We shall skip any discussion of third-order methods in order to come to the principal point of discussion in this section.

A Fourth-Order Runge–Kutta Method

A fourth-order Runge–Kutta procedure consists of finding parameters so that the formula

yn+1 = yn + h(w1k1 + w2k2 + w3k3 + w4k4),(5)

where

agrees with a Taylor polynomial of degree 4. This results in a system of 11 equations in 13 unknowns. The most commonly used set of values for the parameters yields the following result:

(6)

While other fourth-order formulas are easily derived, the algorithm summarized in (6) is so widely used and recognized as a valuable computational tool it is often referred to as the fourth-order Runge–Kutta method or the classical Runge–Kutta method. It is (6) that we have in mind, hereafter, when we use the abbreviation “the RK4 method.”

You are advised to look carefully at the formulas in (6); note that k2 depends on k1, k3 depends on k2, and k4 depends on k3. Also, k2 and k3 involve approximations to the slope at the midpoint xn + h of the interval [xn , xn+ 1].

EXAMPLE 1 RK4 Method

Use the RK4 method with h = 0.1 to obtain an approximation to y(1.5) for the solution of y′ = 2xy, y(1) = 1.

SOLUTION

For the sake of illustration, let us compute the case when n = 0. From (6) we find

and therefore

The remaining calculations are summarized in Table 6.2.1, whose entries are rounded to four decimal places. ≡

TABLE 6.2.1 RK4 Method with h = 0.1
The table captioned R K 4 Method with h equals 0.1 has five columns and six rows. The column headings are as follows: Column 1, x subscript n. Column 2, y subscript n. Column 3, Actual Value. Column 4, Absolute Error. Column 5, Percentage Relative Error. The row entries are as follows: Row 1. x subscript n, 1.00. y subscript n, 1.0000. Actual Value, 1.0000. Absolute Error, 0.0000. Percentage Relative Error, 0.00. Row 2. x subscript n, 1.10. y subscript n, 1.2337. Actual Value, 1.2337. Absolute Error, 0.0000. Percentage Relative Error, 0.00. Row 3. x subscript n, 1.20. y subscript n, 1.5527. Actual Value, 1.5527. Absolute Error, 0.0000. Percentage Relative Error, 0.00. Row 4. x subscript n, 1.30. y subscript n, 1.9937. Actual Value, 1.9937. Absolute Error, 0.0000. Percentage Relative Error, 0.00. Row 5. x subscript n, 1.40. y subscript n, 2.6116. Actual Value, 2.6117. Absolute Error, 0.0001. Percentage Relative Error, 0.00. Row 6. x subscript n, 1.50. y subscript n, 3.4902. Actual Value, 3.4904. Absolute Error, 0.0001. Percentage Relative Error, 0.00.

Inspection of Table 6.2.1 shows why the fourth-order Runge–Kutta method is so popular. If four-decimal-place accuracy is all that we desire, there is no need to use a smaller step size. Table 6.2.2 on page 314 compares the results of applying Euler’s, the improved Euler’s, and the fourth-order Runge–Kutta methods to the initial-value problem y′ = 2xy, y(1) = 1. See Tables 6.1.1–6.1.4.

TABLE 6.2.2 y′ = 2xy, y(1) = 1
The table is captioned y prime equals 2xy, y(1) equals 1. It compares the methods with h equals 0.1 and the methods with h equals 0.05. The first comparison is titled: Comparison of Numerical Methods with h equals 0.1. It has five columns and six rows. The column headings are as follows: Column 1, x subscript n. Column 2, Euler. Column 3, Improved Euler. Column 4, R K 4. Column 5, Actual Value. The row entries are as follows: Row 1. x subscript n, 1.00. Euler, 1.0000. Improved Euler, 1.0000. R K 4, 1.0000. Actual Value, 1.0000. Row 2. x subscript n, 1.10. Euler, 1.2000. Improved Euler, 1.2320. R K 4, 1.2337. Actual Value, 1.2337. Row 3. x subscript n, 1.20. Euler, 1.4640. Improved Euler, 1.5479. R K 4, 1.5527. Actual Value, 1.5527. Row 4. x subscript n, 1.30. Euler, 1.8154. Improved Euler, 1.9832. R K 4, 1.9937. Actual Value, 1.9937. Row 5. x subscript n, 1.40. Euler, 2.2874. Improved Euler, 2.5908. R K 4, 2.6116. Actual Value, 2.6117. Row 6. x subscript n, 1.50. Euler, 2.9278. Improved Euler, 3.4509. R K 4, 3.4902. Actual Value, 3.4904. The second comparison is titled: Comparison of Numerical Methods with h equals 0.05. It has five columns and 11 rows. The column headings are as follows: Column 1, x subscript n. Column 2, Euler. Column 3, Improved Euler. Column 4, R K 4. Column 5, Actual Value. The row entries are as follows: Row 1. x subscript n, 1.00. Euler, 1.0000. Improved Euler, 1.0000. R K 4, 1.0000. Actual Value, 1.0000. Row 2. x subscript n, 1.05. Euler, 1.1000. Improved Euler, 1.1077. R K 4, 1.1079. Actual Value, 1.1079. Row 3. x subscript n, 1.10. Euler, 1.2155. Improved Euler, 1.2332. R K 4, 1.2337. Actual Value, 1.2337. Row 4. x subscript n, 1.15. Euler, 1.3492. Improved Euler, 1.3798. R K 4, 1.3806. Actual Value, 1.3806. Row 5. x subscript n, 1.20. Euler, 1.5044. Improved Euler, 1.5514. R K 4, 1.5527. Actual Value, 1.5527. Row 6. x subscript n, 1.25. Euler, 1.6849. Improved Euler, 1.7531. R K 4, 1.7551. Actual Value, 1.7551. Row 7. x subscript n, 1.30. Euler, 1.8955. Improved Euler, 1.9909. R K 4, 1.9937. Actual Value, 1.9937. Row 8. x subscript n, 1.35. Euler, 2.1419. Improved Euler, 2.2721. R K 4, 2.2762. Actual Value, 2.2762. Row 9. x subscript n, 1.40. Euler, 2.4311. Improved Euler, 2.6060. R K 4, 2.6117. Actual Value, 2.6117. Row 10. x subscript n, 1.45. Euler, 2.7714. Improved Euler, 3.0038. R K 4, 3.0117. Actual Value, 3.0117. Row 11. x subscript n, 1.50. Euler, 3.1733. Improved Euler, 3.4795. R K 4, 3.4903. Actual Value, 3.4904.

Truncation Errors for the RK4 Method

In Section 6.1 we saw that global truncation errors for Euler’s method and for the improved Euler’s method are, respectively, O(h) and O(h2). Because the first equation in (6) agrees with a Taylor polynomial of degree 4, the local truncation error for this method is y(5)(c)h5/5! or O(h5), and the global truncation error is thus O(h4). It is now obvious why Euler’s method, the improved Euler’s method, and (6) are first-, second-, and fourth-order Runge–Kutta methods, respectively.

EXAMPLE 2 Bound for Local Truncation Errors

Find a bound for the local truncation errors for the RK4 method applied to y′ = 2xy, y(1) = 1.

SOLUTION

By computing the fifth derivative of the known solution y(x) = we get

y(5)(c) = (120c + 160c3 + 32c5) .(7)

Thus with c = 1.5, (7) yields a bound of 0.00028 on the local truncation error for each of the five steps when h = 0.1. Note that in Table 6.2.1 the error in y1 is much less than this bound.

Table 6.2.3 gives the approximations to the solution of the initial-value problem at x = 1.5 that are obtained from the RK4 method. By computing the value of the analytic solution at x = 1.5 we can find the error in these approximations. Because the method is so accurate, many decimal places must be used in the numerical solution to see the effect of halving the step size. Note that when h is halved, from h = 0.1 to h = 0.05, the error is divided by a factor of about 24 = 16, as expected. ≡

TABLE 6.2.3 RK4 Method
h Approx. Error
0.1 3.49021064 1.32321089 × 10−4
0.05 3.49033382 9.13776090 × 10−6

≡

Adaptive Methods

We have seen that the accuracy of a numerical method for approximating solutions of differential equations can be improved by decreasing the step size h. Of course, this enhanced accuracy is usually obtained at a cost—namely, increased computation time and greater possibility of round-off error. In general, over the interval of approximation there may be subintervals where a relatively large step size suffices and other subintervals where a smaller step size is necessary in order to keep the truncation error within a desired limit. Numerical methods that use a variable step size are called adaptive methods. One of the more popular of the adaptive routines is the Runge–Kutta–Fehlberg method. Because the German mathematician Erwin Fehlberg (1911–1990) employed two Runge–Kutta methods of differing orders, a fourth- and a fifth-order method, this algorithm is frequently denoted as the RKF45 method.*

6.2 Exercises Answers to selected odd-numbered problems begin on page ANS-14.

  1. Use the RK4 method with h = 0.1 to approximate y(0.5), where y(x) is the solution of the initial-value problem y′ = (x + y − 1)2, y(0) = 2. Compare this approximate value with the actual value obtained in Problem 11 in Exercises 6.1.
  2. Assume that w2 = in (4). Use the resulting second-order Runge–Kutta method to approximate y(0.5), where y(x) is the solution of the initial-value problem in Problem 1. Compare this approximate value with the approximate value obtained in Problem 11 in Exercises 6.1.

In Problems 3–12, use the RK4 method with h = 0.1 to obtain a four-decimal approximation to the indicated value.

  1. y′ = 2x − 3y + 1,     y(1) = 5; y(1.5)
  2. y′ = 4x − 2y,     y(0) = 2; y(0.5)
  3. y′ = 1 + y2,     y(0) = 0; y(0.5)
  4. y′ = x2 + y2,     y(0) = 1; y(0.5)
  5. y′ = e−y,     y(0) = 0; y(0.5)
  6. y′ = x + y2,     y(0) = 0; y(0.5)
  7. y′ = (x − y)2,     y(0) = 0.5; y(0.5)
  8. y′ = xy + ,     y(0) = 1; y(0.5)
  9. y′ = xy2 − ,     y(1) = 1; y(1.5)
  10. y′ = y − y2,     y(0) = 0.5; y(0.5)
  11. If air resistance is proportional to the square of the instantaneous velocity, then the velocity ν of a mass m dropped from a given height h is determined from

    m = mg − kv2, k > 0.

    Let ν(0) = 0, k = 0.125, m = 5 slugs, and g = 32 ft/s2.

    1. Use the RK4 method with h = 1 to approximate the velocity ν(5).
    2. Use a numerical solver to graph the solution of the IVP on the interval [0, 6].
    3. Use separation of variables to solve the IVP and then find the actual value ν(5).
  12. A mathematical model for the area a(in cm2) that a colony of bacteria (B. dendroides) occupies is given by

    = A(2.128 − 0.0432a).

    Suppose that the initial area is 0.24 cm2.

    1. Use the RK4 method with h = 0.5 to complete the following table.
    2. Use a numerical solver to graph the solution of the initial-value problem. Estimate the values A(1), A(2), A(3), A(4), and A(5) from the graph.
    3. Use separation of variables to solve the initial-value problem and compute the values A(1), A(2), A(3), A(4), and a(5).
  13. Consider the initial-value problem y′ = x2 + y3, y(1) = 1. See Problem 12 in Exercises 6.1.
    1. Compare the results obtained from using the RK4 method over the interval [1, 1.4] with step sizes h = 0.1 and h = 0.05.
    2. Use a numerical solver to graph the solution of the initial-value problem on the interval [1, 1.4].
  14. Consider the initial-value problem y′ = 2y, y(0) = 1. The analytic solution is y(x) = e2x.
    1. Approximate y(0.1) using one step and the fourth-order RK4 method.
    2. Find a bound for the local truncation error in y1.
    3. Compare the actual error in y1 with your error bound.
    4. Approximate y(0.1) using two steps and the RK4 method.
    5. Verify that the global truncation error for the RK4 method is O(h4) by comparing the errors in parts (a) and (d).
  15. Repeat Problem 16 using the initial-value problem y′ = −2y + x, y(0) = 1. The analytic solution is

  16. Consider the initial-value problem y′ = 2x − 3y + 1, y(1) = 5. The analytic solution is

    1. Find a formula involving c and h for the local truncation error in the nth step if the RK4 method is used.
    2. Find a bound for the local truncation error in each step if h = 0.1 is used to approximate y(1.5).
    3. Approximate y(1.5) using the RK4 method with h = 0.1 and h = 0.05. See Problem 3. You will need to carry more than six decimal places to see the effect of reducing the step size.
  17. Repeat Problem 18 for the initial-value problem y′ = e−y, y(0) = 0. The analytic solution is y(x) = ln(x + 1). Approximate y(0.5). See Problem 7.

Discussion Problem

  1. A count of the number of evaluations of the function f used in solving the initial-value problem y′ = f(x, y), y(x0) = y0 is used as a measure of the computational complexity of a numerical method. Determine the number of evaluations of f required for each step of Euler’s, the improved Euler’s, and the RK4 methods. By considering some specific examples, compare the accuracy of these methods when used with comparable computational complexities.

Computer Lab Assignment

  1. The RK4 method for solving an initial-value problem over an interval [a, b] results in a finite set of points that are supposed to approximate points on the graph of the exact solution. In order to expand this set of discrete points to an approximate solution defined at all points on the interval [a, b], we can use an interpolating function. This is a function, supported by most computer algebra systems, that agrees with the given data exactly and assumes a smooth transition between data points. These interpolating functions may be polynomials or sets of polynomials joined together smoothly. In Mathematica the command y = Interpolation[data] can be used to obtain an interpolating function through the points data = { {x0, y0}, {x1, y1}, . . . , {xn, yn} }. The interpolating function y[x] can now be treated like any other function built into the computer algebra system.
    1. Find the analytic solution of the initial-value problem y′ = −y + 10 sin 3x; y(0) = 0 on the interval [0, 2]. Graph this solution and find its positive roots.
    2. Use the RK4 method with h = 0.1 to approximate a solution of the initial-value problem in part (a). Obtain an interpolating function and graph it. Find the positive roots of the interpolating function on the interval [0, 2].

 

*The Runge–Kutta method of order four used in RKF45 is not the same as that given in (6).