746 RUNGE-KUTTA METHODS
truncation error is estimated as the difference between two predictions using different-
order RK methods.
25.5.1 Adaptive RK or Step-Halving Method
Step halving (also called adaptive RK) involves taking each step twice, once as a full
step and independently as two half steps. The difference in the two results represents an
estimate of the local truncation error. If y1 designates the single-step prediction and y2
designates the prediction using the two half steps, the error D can be represented as
¢ 5 y2 2 y1 (25.43)
In addition to providing a criterion for step-size control, Eq. (25.43) can also be used to
correct the y2 prediction. For the fourth-order RK version, the correction is
y2 d y2 1 ¢
This estimate is ! fth-order accurate.
EXAMPLE 25.12 Adaptive Fourth-Order RK Method
Problem Statement. Use the adaptive fourth-order RK method to integrate y9 5 4e0.8x 2 0.5y from x 5 0 to 2 using h 5 2 and an initial condition of y(0) 5 2. This is the same
differential equation that was solved previously in Example 25.5. Recall that the true
solutions is y(2) 5 14.84392.
Solution. The single prediction with a step of h is computed as
y(2) 5 2 1 1
6 [3 1 2(6.40216 1 4.70108) 1 14.11105]2 5 15.10584
The two half-step predictions are
y(1) 5 2 1 1
6 [3 1 2(4.21730 1 3.91297) 1 5.945681]1 5 6.20104
y(2) 5 6.20104 1 1
6 [5.80164 1 2(8.72954 1 7.99756) 1 12.71283]1 5 14.86249
Therefore, the approximate error is
Ea 5 14.86249 2 15.10584
15 5 20.01622
which compares favorably with the true error of
Et 5 14.84392 2 14.86249 5 20.01857
The error estimate can also be used to correct the prediction
y(2) 5 14.86249 2 0.01622 5 14.84627
which has an Et 5 20.00235.
25.5 ADAPTIVE RUNGE-KUTTA METHODS 747
25.5.2 Runge-Kutta Fehlberg
Aside from step halving as a strategy to adjust step size, an alternative approach for
obtaining an error estimate involves computing two RK predictions of different order.
The results can then be subtracted to obtain an estimate of the local truncation error. One
shortcoming of this approach is that it greatly increases the computational overhead. For
example, a fourth- and ! fth-order prediction amount to a total of 10 function evaluations
per step. The Runge-Kutta Fehlberg or embedded RK method cleverly circumvents this
problem by using a ! fth-order RK method that employs the function evaluations from
the accompanying fourth-order RK method. Thus, the approach yields the error estimate
on the basis of only six function evaluations!
For the present case, we use the following fourth-order estimate
yi11 5 yi 1 a 37 378
k1 1 250
621 k3 1
594 k4 1
1771 k6b h (25.45)
along with the ! fth-order formula:
yi11 5 yi 1 a 2825 27,648
k1 1 18,575
48,384 k3 1
55,296 k4 1
14,336 k5 1
4 k6b h (25.46)
k1 5 f(xi, yi)
k2 5 f axi 1 1 5
h, yi 1 1
k3 5 f axi 1 3 10
h, yi 1 3
40 k1h 1
k4 5 f axi 1 3 5
h, yi 1 3
10 k1h 2
10 k2h 1
k5 5 f axi 1 h, yi 2 11 54
k1h 1 5
2 k2h 2
27 k3h 1
k6 5 f axi 1 7 8
h, yi 1 1631
55,296 k1h 1
512 k2h 1
13,824 k3h 1
Thus, the ODE can be solved with Eq. (25.46) and the error estimated as the difference
of the ! fth- and fourth-order estimates. It should be noted that the particular coef! cients
used above were developed by Cash and Karp (1990). Therefore, it is sometimes called
the Cash-Karp RK method.
EXAMPLE 25.13 Runge-Kutta Fehlberg Method
Problem Statement. Use the Cash-Karp version of the Runge-Kutta Fehlberg approach to perform the same calculation as in Example 25.12 from x 5 0 to 2 using h 5 2.
748 RUNGE-KUTTA METHODS
Solution. The calculation of the k’s can be summarized in the following table:
x y f(x, y)
k1 0 2 3 k2 0.4 3.2 3.908511 k3 0.6 4.20883 4.359883 k4 1.2 7.228398 6.832587 k5 2 15.42765 12.09831 k6 1.75 12.17686 10.13237
These can then be used to compute the fourth-order prediction
y1 5 2 1 a 37 378
3 1 250
621 4.359883 1
594 6.832587 1
1771 10.13237b 2 5 14.83192
along with a ! fth-order formula:
y1 5 2 1 a 2825 27,648
3 1 18,575
48,384 4.359883 1
14,336 12.09831 1
4 10.13237b 2 5 14.83677
The error estimate is obtained by subtracting these two equations to give
Ea 5 14.83677 2 14.83192 5 0.004842
25.5.3 Step-Size Control
Now that we have developed ways to estimate the local truncation error, it can be used
to adjust the step size. In general, the strategy is to increase the step size if the error is
too small and decrease it if the error is too large. Press et al. (2007) have suggested the
following criterion to accomplish this:
hnew 5 hpresent ` ¢new ¢present
` a (25.47) where hpresent and hnew 5 the present and the new step sizes, respectively, Dpresent 5 the
computed present accuracy, Dnew 5 the desired accuracy, and a 5 a constant power that
is equal to 0.2 when the step size is increased (that is, when Dpresent # Dnew) and 0.25
when the step size is decreased (Dpresent . Dnew).
The key parameter in Eq. (25.47) is obviously Dnew because it is your vehicle for
specifying the desired accuracy. One way to do this would be to relate Dnew to a rela-
tive error level. Although this works well when only positive values occur, it can cause
problems for solutions that pass through zero. For example, you might be simulating
an oscillating function that repeatedly passes through zero but is bounded by maximum
absolute values. For such a case, you might want these maximum values to ! gure in
the desired accuracy.
25.5 ADAPTIVE RUNGE-KUTTA METHODS 749
A more general way to handle such cases is to determine Dnew as
¢new 5 eyscale
where e 5 an overall tolerance level. Your choice of yscale will then determine how the error
is scaled. For example, if yscale 5 y, the accuracy will be couched in terms of fractional
relative errors. If you are dealing with a case where you desire constant errors relative to
a prescribed maximum bound, set yscale equal to that bound. A trick suggested by Press
et al. (2007) to obtain the constant relative errors except very near zero crossings is
yscale 5 Zy Z 1 ` h dy dx `
This is the version we will use in our algorithm.
25.5.4 Computer Algorithm
Figures 25.21 and 25.22 outline pseudocode to implement the Cash-Karp version of the
Runge-Kutta Fehlberg algorithm. This algorithm is patterned after a more detailed imple-
mentation by Press et al. (2007) for systems of ODEs.
Figure 25.21 implements a single step of the Cash-Karp routine (that is Eqs. 25.45
and 25.46). Figure 25.22 outlines a general driver program along with a subroutine that
actually adapts the step size.
SUBROUTINE RKkc (y,dy,x,h,yout,yerr)
CALL Derivs (x1a2*h,ytemp,k2)
FIGURE 25.21 Pseudocode for a single step of the Cash-Karp RK method.
750 RUNGE-KUTTA METHODS
EXAMPLE 25.14 Computer Application of an Adaptive Fourth-Order RK Scheme
Problem Statement. The adaptive RK method is well-suited for the following ordinary differential equation
dx 1 0.6y 5 10e2(x22)
Notice for the initial condition, y(0) 5 0.5, the general solution is
y 5 0.5e20.6x (E25.14.2)
which is a smooth curve that gradually approaches zero as x increases. In contrast, the
particular solution undergoes an abrupt transition in the vicinity of x 5 2 due to the nature
of the forcing function (Fig. 25.23a). Use a standard fourth-order RK scheme to solve
Eq. (E25.14.1) from x 5 0 to 4. Then employ the adaptive scheme described in this sec-
tion to perform the same computation.
Solution. First, the classical fourth-order scheme is used to compute the solid curve in Fig. 25.23b. For this computation, a step size of 0.1 is used so that 4y(0.1) 5 40 applica- tions of the technique are made. Then, the calculation is repeated with a step size of 0.05
for a total of 80 applications. The major discrepancy between the two results occurs in the
region from 1.8 to 2.0. The magnitude of the discrepancy is about 0.1 to 0.2 percent.
(a) Driver Program
INPUT xi, xf, yi
hi5.5; tiny 5 1. 3 10230
print *, xi,yi
IF (istep . maxstep AND x # xf) EXIT
IF (x1h.xf) THEN h5xf2x
CALL Adapt (x,y,dy,h,yscal,eps,hnxt)
(b) Adaptive Step Routine
SUB Adapt (x,y,dy,htry,yscal,eps,hnxt)
PARAMETER (safety50.9, econ51.89e24)
CALL RKkc (y,dy,x,h,ytemp,yerr)
IF emax # 1 EXIT
IF xnew5x THEN pause
IF emax . econ THEN
hnxt5safety*emax 2.2 *h
FIGURE 25.22 Pseudocode for a (a) driver program and an (b) adaptive step routine to solve a single ODE.
25.5 ADAPTIVE RUNGE-KUTTA METHODS 751
Next, the algorithm in Figs. 25.21 and 25.22 is developed into a computer program
and used to solve the same problem. An initial step size of 0.5 and an e 5 0.00005 were
chosen. The results were superimposed on Fig. 25.23b. Notice how large steps are taken
in the regions of gradual change. Then, in the vicinity of x 5 2, the steps are decreased
to accommodate the abrupt nature of the forcing function.
FIGURE 25.23 (a) A bell-shaped forcing function that induces an abrupt change in the solution of an ODE [Eq. (E25.14.1)]. (b) The solution. The points indicate the predictions of an adaptive step-size routine.
0 2 4 x
0 2 4 x
The utility of an adaptive integration scheme obviously depends on the nature of the
functions being modeled. It is particularly advantageous for those solutions with long
smooth stretches and short regions of abrupt change. In addition, it has utility in those
situations where the correct step size is not known a priori. For these cases, an adaptive
routine will “feel” its way through the solution while keeping the results within the
desired tolerance. Thus, it will tiptoe through the regions of abrupt change and step out
briskly when the variations become more gradual.
752 RUNGE-KUTTA METHODS
25.1 Solve the following initial value problem over the interval from
t 5 0 to 2 where y(0) 5 1. Display all your results on the same graph.
dt 5 yt 2 2 1.1y
(b) Euler’s method with h 5 0.5 and 0.25.
(c) Midpoint method with h 5 0.5.
(d) Fourth-order RK method with h 5 0.5.
25.2 Solve the following problem over the interval from x 5 0 to 1
using a step size of 0.25 where y(0) 5 1. Display all your results on
the same graph.
dt 5 (1 1 4t)1y
(b) Euler’s method.
(c) Heun’s method without iteration.
(d) Ralston’s method.
(e) Fourth-order RK method.
25.3 Use the (a) Euler and (b) Heun (without iteration) methods to
dt 2 2 0.5t 1 y 5 0
where y(0) 5 2 and y9(0) 5 0. Solve from x 5 0 to 4 using h 5 0.1.
Compare the methods by plotting the solutions.
25.4 Solve the following problem with the fourth-order RK method:
dx 2 1 0.6
dx 1 8y 5 0
where y(0) 5 4 and y9(0) 5 0. Solve from x 5 0 to 5 with h 5 0.5.
Plot your results.
25.5 Solve from t 5 0 to 3 with h 5 0.1 using (a) Heun (without
corrector) and (b) Ralston’s second-order RK method:
dt 5 y sin3(t) y(0) 5 1
25.6 Solve the following problem numerically from t 5 0 to 3:
dt 5 22y 1 t2 y(0) 5 1
Use the third-order RK method with a step size of 0.5.
25.7 Use (a) Euler’s and (b) the fourth-order RK method to solve
dt 5 22y 1 5e2t
dt 5 2
over the range t 5 0 to 0.4 using a step size of 0.1 with y(0) 5 2 and
z(0) 5 4.
25.8 Compute the ” rst step of Example 25.14 using the adaptive
fourth-order RK method with h 5 0.5. Verify whether step-size
adjustment is in order.
25.9 If e 5 0.001, determine whether step size adjustment is re-
quired for Example 25.12.
25.10 Use the RK-Fehlberg approach to perform the same calcula-
tion as in Example 25.12 from x 5 0 to 1 with h 5 1.
25.11 Write a computer program based on Fig. 25.7. Among other
things, place documentation statements throughout the program to
identify what each section is intended to accomplish.
25.12 Test the program you developed in Prob. 25.11 by duplicat-
ing the computations from Examples 25.1 and 25.4.
25.13 Develop a user-friendly program for the Heun method with
an iterative corrector. Test the program by duplicating the results in
25.14 Develop a user-friendly computer program for the classical
fourth-order RK method. Test the program by duplicating Exam-
25.15 Develop a user-friendly computer program for systems of
equations using the fourth-order RK method. Use this program to
duplicate the computation in Example 25.10.
25.16 The motion of a damped spring-mass system (Fig. P25.16)
is described by the following ordinary differential equation:
m d 2x
dt2 1 c
dt 1 kx 5 0
where x 5 displacement from equilibrium position (m), t 5 time
(s), m 5 20-kg mass, and c 5 the damping coef” cient (N ? s/m).
The damping coef” cient c takes on three values of 5 (under-
damped), 40 (critically damped), and 200 (overdamped). The
spring constant k 5 20 N/m. The initial velocity is zero, and the
initial displacement x 5 1 m. Solve this equation using a numerical
method over the time period 0 # t # 15 s. Plot the displacement
versus time for each of the three values of the damping coef” cient
on the same curve.
25.21 The logistic model is used to simulate population as in
dt 5 kgm(1 2 pypmax)p
where p 5 population, kgm 5 the maximum growth rate under un-
limited conditions, and pmax 5 the carrying capacity. Simulate the
world’s population from 1950 to 2000 using one of the numerical
methods described in this chapter. Employ the following initial
conditions and parameter values for your simulation: p0 (in 1950) 5
2555 million people, kgm 5 0.026/yr, and pmax 5 12,000 million
people. Have the function generate output corresponding to the
dates for the following measured population data. Develop a plot of
your simulation along with these data.
t 1950 1960 1970 1980 1990 2000
p 2555 3040 3708 4454 5276 6079
25.22 Suppose that a projectile is launched upward from the
earth’s surface. Assume that the only force acting on the object is
the downward force of gravity. Under these conditions, a force
balance can be used to derive,
dt 5 2g(0)
(R 1 x)2
where y 5 upward velocity (m/s), t 5 time (s), x 5 altitude (m)
measured upwards from the earth’s surface, g(0) 5 the gravita-
tional acceleration at the earth’s surface (> 9.81 m/s2), and R 5 the earth’s radius (> 6.37 3 106 m). Recognizing that dx/dt 5 y, use Euler’s method to determine the maximum height that would be
obtained if y(t 5 0) 5 1500 m/s.
25.23 The following function exhibits both # at and steep regions
over a relatively short x region:
f (x) 5 1
(x 2 0.3)2 1 0.01 1
(x 2 0.9)2 1 0.04 2 6
Determine the value of the de” nite integral of this function between
x 5 0 and 1 using an adaptive RK method.
25.17 If water is drained from a vertical cylindrical tank by open-
ing a valve at the base, the water will # ow fast when the tank is full
and slow down as it continues to drain. As it turns out, the rate at
which the water level drops is:
dt 5 2k1y
where k is a constant depending on the shape of the hole and the
cross-sectional area of the tank and drain hole. The depth of the
water y is measured in meters and the time t in minutes. If k 5 0.06,
determine how long it takes the tank to drain if the # uid level is
initially 3 m. Solve by applying Euler’s equation and writing a
computer program or using Excel. Use a step of 0.5 minutes.
25.18 The following is an initial value, second-order differential
dt 2 1 (5x)
dt 1 (x 1 7) sin(vt) 5 0
dt (0) 5 1.5 and x(0) 5 6
Note that v 5 1. Decompose the equation into two ” rst-order dif-
ferential equations. After the decomposition, solve the system from
t 5 0 to 15 and plot the results.
25.19 Assuming that drag is proportional to the square of velocity,
we can model the velocity of a falling object like a parachutist with
the following differential equation:
dt 5 g 2
where y is velocity (m/s), t 5 time (s), g is the acceleration due to
gravity (9.81 m/s2), cd 5 a second-order drag coef” cient (kg/m),
and m 5 mass (kg). Solve for the velocity and distance fallen by a
90-kg object with a drag coef” cient of 0.225 kg/m. If the initial
height is 1 km, determine when it hits the ground. Obtain your solu-
tion with (a) Euler’s method and (b) the fourth-order RK method.
25.20 A spherical tank has a circular ori” ce in its bottom through
which the liquid # ows out (Fig. P25.20). The # ow rate through the
hole can be estimated as
Qout 5 CA12gH where Qout 5 out# ow (m
3/s), C 5 an empirically-derived coef” –
cient, A 5 the area of the ori” ce (m2), g 5 the gravitational con-
stant (5 9.81 m/s2), and H 5 the depth of liquid in the tank. Use
one of the numerical methods described in this chapter to determine
how long it will take for the water to # ow out of a 3-m-diameter
tank with an initial height of 2.75 m. Note that the ori” ce has a di-
ameter of 3 cm and C 5 0.55.