Engineering

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 ¢

15 (25.44)

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

and

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

125

594 k4 1

512

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

13,525

55,296 k4 1

277

14,336 k5 1

4 k6b h (25.46)

where

k1 5 f(xi, yi)

k2 5 f axi 1 1 5

h, yi 1 1

5 k1hb

k3 5 f axi 1 3 10

h, yi 1 3

40 k1h 1

40 k2hb

k4 5 f axi 1 3 5

h, yi 1 3

10 k1h 2

10 k2h 1

5 k3hb

k5 5 f axi 1 h, yi 2 11 54

k1h 1 5

2 k2h 2

27 k3h 1

27 k4hb

k6 5 f axi 1 7 8

h, yi 1 1631

55,296 k1h 1

175

512 k2h 1

575

13,824 k3h 1

44,275

110,592 k4h

1 253

4096 k5hb

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

125

594 6.832587 1

512

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

13,525

55,296 6.832587

1 227

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)

PARAMETER (a250.2,a350.3,a450.6,a551.,a650.875,

b2150.2,b3153.y40.,b3259.y40.,b4150.3,b42520.9,

b4351.2,b515211.y54.,b5252.5,b535270.y27.,

b54535.y27.,b6151631.y55296.,b625175.y512.,

b635575.y13824.,b64544275.y110592.,b655253.y4096.,

c1537.y378.,c35250.y621.,c45125.y594.,

c65512.y1771.,dc15c122825.y27648.,

dc35c3218575.y48384.,dc45c4213525.y55296.,

dc552277.y14336.,dc65c620.25)

ytemp5y1b21*h*dy

CALL Derivs (x1a2*h,ytemp,k2)

ytemp5y1h*(b31*dy1b32*k2)

CALL Derivs(x1a3*h,ytemp,k3)

ytemp5y1h*(b41*dy1b42*k21b43*k3)

CALL Derivs(x1a4*h,ytemp,k4)

ytemp5y1h*(b51*dy1b52*k21b53*k31b54*k4)

CALL Derivs(x1a5*h,ytemp,k5)

ytemp5y1h*(b61*dy1b62*k21b63*k31b64*k41b65*k5)

CALL Derivs(x1a6*h,ytemp,k6)

yout5y1h*(c1*dy1c3*k31c4*k41c6*k6)

yerr5h*(dc1*dy1dc3*k31dc4*k41dc5*k51dc6*k6)

END RKkc

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)

2y[2(0.075)2] (E25.14.1)

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

maxstep5100

hi5.5; tiny 5 1. 3 10230

eps50.00005

print *, xi,yi

x5xi

y5yi

h5hi

istep50

IF (istep . maxstep AND x # xf) EXIT

istep5istep11

CALL Derivs(x,y,dy)

yscal5ABS(y)1ABS(h*dy)1tiny

IF (x1h.xf) THEN h5xf2x

CALL Adapt (x,y,dy,h,yscal,eps,hnxt)

PRINT x,y

h5hnxt

END DO

END

(b) Adaptive Step Routine

SUB Adapt (x,y,dy,htry,yscal,eps,hnxt)

PARAMETER (safety50.9, econ51.89e24)

h5htry

CALL RKkc (y,dy,x,h,ytemp,yerr)

emax5abs(yerr/yscal/eps)

IF emax # 1 EXIT

htemp5safety*h*emax 20.25

h5max(abs(htemp),0.25*abs(h))

xnew5x1h

IF xnew5x THEN pause

END DO

IF emax . econ THEN

hnxt5safety*emax 2.2 *h

ELSE

hnxt54.*h

END IF

x5x1h

y5ytemp

END Adapt

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

(b)

0 2 4 x

(a)

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

PROBLEMS

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

(a) Analytically.

(b) Euler’s method with h 5 0.5 and 0.25.

(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

(a) Analytically.

(b) Euler’s method.

(d) Ralston’s method.

(e) Fourth-order RK method.

25.3 Use the (a) Euler and (b) Heun (without iteration) methods to

solve

d 2y

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:

d 2y

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

yz2

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

Table 25.2.

25.14 Develop a user-friendly computer program for the classical

fourth-order RK method. Test the program by duplicating Exam-

ple 25.7.

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.

FIGURE P25.16

PROBLEMS 753

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

equation:

d 2x

dt 2 1 (5x)

dt 1 (x 1 7) sin(vt) 5 0

where

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

m y2

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.

Engineering

Summary of Our Reviews

Why Us