Download (direct link):
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 < W It is given by
yn+1 = yn + h (k- + + 2*n3 + ), (2)
kn1 = f (tn, yn)
kn2 = f(tn + 2h> yn + 2hkmh
kn3 = f(n + 2h, yn + 1 hkn2)>
kn4 = f (tn + h> yn + 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 0 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
*Carl 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
yn+1 — yn = 6[f (n) +4 f (n +h/2) + f (n +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 = $(t) of the initial value problem
y = 1 — t + 4 y, 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.
For comparison, note that both the Runge-Kutta method with h = 0.05 and the improved Euler method with h = 0.025 require 160 evaluations of f to reach t = 2. The improved Euler method yields a result at t = 2 that is in error by 1.22%. While this error may be acceptable for some purposes, it is more than 135 times the error yielded by the Runge-Kutta method with comparable computing effort. Note also that the Runge-Kutta method with h = 0.2, or 40 evaluations of f, produces a value at t = 2
Chapter 8. Numerical Methods
with an error of 1.40%, which is only slightly greater than the error in the improved Euler method with h = 0.025, or 160 evaluations of f. Thus we see again that a more accurate algorithm is more efficient; it produces better results with similar effort, or similar results with less effort.
TABLE 8.3.1 A Comparison of Results for the Numerical Solution of the Initial Value Problem y = 1 — t + 4y, y(0) = 1
t Improved Euler Runge-Kutta Exact