Download (direct link):
8.3 The Runge-Kutta Method
In preceding sections we have introduced the Euler formula, the backward Euler formula, and the improved Euler formula as ways to solve the initial value problem
y = f (t, ÓÕ y(t0) = Ó0
Chapter 8. Numerical Methods
numerically. The local truncation errors for these methods are proportional to h2, h2, and h3, respectively. The Euler and improved Euler methods belong to what is now called the Runge-Kutta1 class of methods.
In this section we discuss the method originally developed by Runge and Kutta. This method is now called the classic fourth order four-stage Runge-Kutta method, but it is often referred to simply as the Runge-Kutta method, and we will follow this practice
for brevity. This method has a local truncation error that is proportional to h5. Thus it
is two orders of magnitude more accurate than the improved Euler method and three orders of magnitude better than the Euler method. It is relatively simple to use and is sufficiently accurate to handle many problems efficiently. This is especially true of adaptive Runge-Kutta methods in which a provision is made to vary the step size as needed.
The Runge-Kutta formula involves a weighted average of values of f (t, y) at different points in the interval tn < t < Ï+1. It is given by
Óï+1 = ÓÏ + h ( kn1 + ^ + 2kn3 + kn4) , (2)
kn1 = f (tn, Óï)
kn2 = f(tn + 2h> Óï + 2hkrnh kn3 = f(tn + -2h> Óï + 1 hkn2)’ kn4 = f (tn + h> Óï + hkn3).
The sum (kn1 + 2kn2 + 2kn3 + kn4)/6 can be interpreted as an average slope. Note that kn1 is the slope at the left end of the interval, kn2 is the slope at the midpoint using the Euler formula to go from tn to tn + h/2, kn3 is a second approximation to the slope at the midpoint, and finally kn4 is the slope at tn + h using the Euler formula and the slope kn3 to go from tn to tn + h.
While in principle it is not difficult to show that Eq. (2) differs from the Taylor expansion of the solution ô by terms that are proportional to h5, the algebra is rather lengthy.2 Thus we will accept the fact that the local truncation error in using Eq. (2) is proportional to h5 and that for a finite interval the global truncation error is at most a constant times h4.
Clearly the Runge-Kutta formula, Eqs. (2) and (3), is more complicated than any of the formulas discussed previously. This is of relatively little significance, however, since it is not hard to write a computer program to implement this method. Such a program has the same structure as the algorithm for the Euler method outlined in
1Carl David Runge (1856-1927), German mathematician and physicist, worked for many years in spectroscopy. The analysis of data led him to consider problems in numerical computation, and the Runge-Kutta method originated in his paper on the numerical solution of differential equations in 1895. The method was extended to systems of equations in 1901 by M. Wilhelm Kutta (1867-1944). Kutta was a German mathematician and aerodynamicist who is also well known for his important contributions to classical airfoil theory.
2See, for example, Chapter 3 of the book by Henrici listed in the references.
8.3 The Runge-Kutta Method
Section 8.1. To be specific, the lines in Step 6 in the Euler algorithm must be replaced by the following:
Step 6. k 1 = f (t, y)
k2 = f (t + 0.5 * h, y + 0.5 * h * k 1)
k3 = f (t + 0.5 * h, y + 0.5 * h * k2)
k4 = f (t + h, y + h * k3)
y = y + (h/6) * (k 1 + 2 * k2 + 2 * k3 + k4)
t = t + h
Note that if f does not depend on y, then
kn1 = f(tn ) kn2 = kn3 = f(tn + h/2), kn4 = f(tn + h)> (4)
and Eq. (2) reduces to
Óï+1 - Óï = 6[ f (tn) + 4 f (tn + h/2) + f (tn + h)]. (5)
Equation (5) can be identified as Simpson’s (1710-1761) rule for the approximate evaluation of the integral of Ó = f(t). The fact that Simpson’s rule has an error proportional to h5 is consistent with the local truncation error in the Runge-Kutta formula.
Use the Runge-Kutta method to calculate approximate values of the solution y = ô() of the initial value problem
/ = 1 - t + 4y, y(0) = 1. (6)
Taking h = 0.2 we have
k01 = f(0, 1) = 5; hk01 = 1.0, k02 = f (0 + 0.1, 1 + 0.5) = 6.9; hk02 = 1.38, k03 = f (0 + 0.1, 1 + 0.69) = 7.66; hk03 = 1.532,
k04 = f (0 + 0.2, 1 + 1.532) = 10.928.
y1 = 1 + —[5 + 2(6.9) + 2(7.66) + 10.928]
= 1 + 1.5016 = 2.5016.
Further results using the Runge-Kutta method with h = 0.2, h = 0.1, and h = 0.05 are given in Table 8.3.1. Note that the Runge-Kutta method yields a value at t = 2 that differs from the exact solution by only 0.122% if the step size is h = 0.1, and by only 0.00903% if h = 0.05. In the latter case the error is less than one part in ten thousand, and the calculated value at t = 2 is correct to four digits.