# Engineering

FIGURE 25.14 Comparison of the true solution with numerical solutions using three second-order RK methods and Euler’s method.

y

4

0 x420

Analytical Euler Heun Midpoint Ralston

734 RUNGE-KUTTA METHODS

For Ralston’s method, k1 for the ! rst interval also equals 8.5 and [Eq. (25.38b)]

k2 5 22(0.375) 3

1 12(0.375)2 2 20(0.375) 1 8.5 5 2.58203125

The average slope is computed by

f 5 1

3 (8.5) 1

2

3 (2.58203125) 5 4.5546875

which can be used to predict

y(0.5) 5 1 1 4.5546875(0.5) 5 3.27734375 et 5 21.82%

The computation is repeated, and the results are summarized in Fig. 25.14 and Table

25.3. Notice how all the second-order RK methods are superior to Euler’s method.

25.3.2 Third-Order Runge-Kutta Methods

For n 5 3, a derivation similar to the one for the second-order method can be performed.

The result of this derivation is six equations with eight unknowns. Therefore, values for

two of the unknowns must be speci! ed a priori in order to determine the remaining

parameters. One common version that results is

yi11 5 yi 1 1

6 (k1 1 4k2 1 k3)h (25.39)

where

k1 5 f(xi, yi) (25.39a)

TABLE 25.3 Comparison of true and approximate values of the integral of y9 5 22×3 1 12×2 2 20x 1 8.5, with the initial condition that y 5 1 at x 5 0. The approximate values were computed using three versions of second-order RK methods with a step size of 0.5.

Second-Order Heun Midpoint Ralston RK

x ytrue y |Et| (%) y |Et| (%) y |Et| (%)

0.0 1.00000 1.00000 0 1.00000 0 1.00000 0 0.5 3.21875 3.43750 6.8 3.109375 3.4 3.277344 1.8 1.0 3.00000 3.37500 12.5 2.81250 6.3 3.101563 3.4 1.5 2.21875 2.68750 21.1 1.984375 10.6 2.347656 5.8 2.0 2.00000 2.50000 25.0 1.75 12.5 2.140625 7.0 2.5 2.71875 3.18750 17.2 2.484375 8.6 2.855469 5.0 3.0 4.00000 4.37500 9.4 3.81250 4.7 4.117188 2.9 3.5 4.71875 4.93750 4.6 4.609375 2.3 4.800781 1.7 4.0 3.00000 3.00000 0 3 0 3.031250 1.0

25.3 RUNGE-KUTTA METHODS 735

k2 5 f axi 1 1 2

h, yi 1 1

2 k1hb (25.39b)

k3 5 f(xi 1 h, yi 2 k1h 1 2k2h) (25.39c)

Note that if the derivative is a function of x only, this third-order method reduces to

Simpson’s 1y3 rule. Ralston (1962) and Ralston and Rabinowitz (1978) have developed an alternative version that provides a minimum bound on the truncation error. In any

case, the third-order RK methods have local and global errors of O(h4) and O(h3), re-

spectively, and yield exact results when the solution is a cubic. When dealing with

polynomials, Eq. (25.39) will also be exact when the differential equation is cubic and

the solution is quartic. This is because Simpson’s 1y3 rule provides exact integral esti- mates for cubics (recall Box 21.3).

25.3.3 Fourth-Order Runge-Kutta Methods

The most popular RK methods are fourth order. As with the second-order approaches,

there are an in! nite number of versions. The following is the most commonly used form,

and we therefore call it the classical fourth-order RK method:

yi11 5 yi 1 1

6 (k1 1 2k2 1 2k3 1 k4)h (25.40)

where

k1 5 f(xi, yi) (25.40a)

k2 5 f axi 1 1 2

h, yi 1 1

2 k1hb (25.40b)

FIGURE 25.15 Graphical depiction of the slope estimates comprising the fourth-order RK method.

y

xxi+1/2

h

xi

k2

k1

k3

k3

k2

k1

k4

xi+1

736 RUNGE-KUTTA METHODS

k3 5 f axi 1 1 2

h, yi 1 1

2 k2hb (25.40c)

k4 5 f(xi 1 h, yi 1 k3h) (25.40d)

Notice that for ODEs that are a function of x alone, the classical fourth-order RK

method is similar to Simpson’s 1y3 rule. In addition, the fourth-order RK method is similar to the Heun approach in that multiple estimates of the slope are developed in order

to come up with an improved average slope for the interval. As depicted in Fig. 25.15,

each of the k’s represents a slope. Equation (25.40) then represents a weighted average

of these to arrive at the improved slope.

EXAMPLE 25.7 Classical Fourth-Order RK Method

Problem Statement.

(a) Use the classical fourth-order RK method [Eq. (25.40)] to integrate

f(x, y) 5 22×3 1 12×2 2 20x 1 8.5

using a step size of h 5 0.5 and an initial condition of y 5 1 at x 5 0.

(b) Similarly, integrate

f(x, y) 5 4e0.8x 2 0.5y

using h 5 0.5 with y(0) 5 2 from x 5 0 to 0.5.

Solution.

(a) Equations (25.40a) through (25.40d) are used to compute k1 5 8.5, k2 5 4.21875,

k3 5 4.21875 and k4 5 1.25, which are substituted into Eq. (25.40) to yield

y(0.5) 5 1 1 e 1 6

[8.5 1 2(4.21875) 1 2(4.21875) 1 1.25] f 0.5 5 3.21875

which is exact. Thus, because the true solution is a quartic [Eq. (PT7.16)], the fourth-

order method gives an exact result.

(b) For this case, the slope at the beginning of the interval is computed as

k1 5 f(0, 2) 5 4e 0.8(0)

2 0.5(2) 5 3

This value is used to compute a value of y and a slope at the midpoint,

y(0.25) 5 2 1 3(0.25) 5 2.75

k2 5 f(0.25, 2.75) 5 4e 0.8(0.25)

2 0.5(2.75) 5 3.510611

This slope in turn is used to compute another value of y and another slope at the midpoint,

y(0.25) 5 2 1 3.510611(0.25) 5 2.877653

k3 5 f(0.25, 2.877653) 5 4e 0.8(0.25)

2 0.5(2.877653) 5 3.446785

Next, this slope is used to compute a value of y and a slope at the end of the interval,

y(0.5) 5 2 1 3.071785(0.5) 5 3.723392

k4 5 f(0.5, 3.723392) 5 4e 0.8(0.5)

2 0.5(3.723392) 5 4.105603

25.3 RUNGE-KUTTA METHODS 737

Finally, the four slope estimates are combined to yield an average slope. This average

slope is then used to make the ! nal prediction at the end of the interval.

f 5 1

6 [3 1 2(3.510611) 1 2(3.446785) 1 4.105603] 5 3.503399

y(0.5) 5 2 1 3.503399(0.5) 5 3.751699

which compares favorably with the true solution of 3.751521.

25.3.4 Higher-Order Runge-Kutta Methods

Where more accurate results are required, Butcher’s (1964) fth-order RK method is

recommended:

yi11 5 yi 1 1

90 (7k1 1 32k3 1 12k4 1 32k5 1 7k6)h (25.41)

where

k1 5 f(xi, yi) (25.41a)

k2 5 f axi 1 1 4

h, yi 1 1

4 k1hb (25.41b)

k3 5 f axi 1 1 4

h, yi 1 1

8 k1h 1

1

8 k2hb (25.41c)

k4 5 f axi 1 1 2

h, yi 2 1

2 k2h 1 k3hb (25.41d)

k5 5 f axi 1 3 4

h, yi 1 3

16 k1h 1

9

16 k4hb (25.41e)

k6 5 f axi 1 h, yi 2 3 7

k1h 1 2

7 k2h 1

12

7 k3h 2

12

7 k4h 1

8

7 k5hb (25.41f)

Note the similarity between Butcher’s method and Boole’s rule in Table 21.2. Higher-order

RK formulas such as Butcher’s method are available, but in general, beyond fourth-order

methods the gain in accuracy is offset by the added computational effort and complexity.

EXAMPLE 25.8 Comparison of Runge-Kutta Methods

Problem Statement. Use ! rst- through ! fth-order RK methods to solve

f(x, y) 5 4e0.8x 2 0.5y

with y(0) 5 2 from x 5 0 to x 5 4 with various step sizes. Compare the accuracy of

the various methods for the result at x 5 4 based on the exact answer of y(4) 5 75.33896.

Solution. The computation is performed using Euler’s, the noniterative Heun, the third- order RK [Eq. (25.39)], the classical fourth-order RK, and Butcher’s ! fth-order RK

methods. The results are presented in Fig. 25.16, where we have plotted the absolute

738 RUNGE-KUTTA METHODS

value of the percent relative error versus the computational effort. This latter quantity is

equivalent to the number of function evaluations required to attain the result, as in

Effort 5 nf b 2 a

h (E25.8.1)

where nf 5 the number of function evaluations involved in the particular RK computa-

tion. For orders # 4, nf is equal to the order of the method. However, note that Butcher’s

! fth-order technique requires six function evaluations [Eq. (25.41a) through (25.41f)].

The quantity (b 2 a)yh is the total integration interval divided by the step size—that is, it is the number of applications of the RK technique required to obtain the result. Thus,

because the function evaluations are usually the primary time-consuming steps, Eq. (E25.8.1)

provides a rough measure of the run time required to attain the answer.

Inspection of Fig. 25.16 leads to a number of conclusions: ! rst, that the higher-order

methods attain better accuracy for the same computational effort and, second, that the

gain in accuracy for the additional effort tends to diminish after a point. (Notice that the

curves drop rapidly at ! rst and then tend to level off.)

Example 25.8 and Fig. 25.16 might lead one to conclude that higher-order RK tech-

niques are always the preferred methods. However, other factors such as programming

FIGURE 25.16 Comparison of percent relative error versus computational effort for fi rst- through fi fth-order RK methods.

100

1

10– 2

10– 4

10– 6

Euler

Heun

RK–3d

RK–4th

Butcher

Effort

P e rc

e n

t re

la ti

v e e

rr o

r

25.4 SYSTEMS OF EQUATIONS 739

costs and the accuracy requirements of the problem also must be considered when choos-

ing a solution technique. Such trade-offs will be explored in detail in the engineering

applications in Chap. 28 and in the epilogue for Part Seven.

25.3.5 Computer Algorithms for Runge-Kutta Methods

As with all the methods covered in this chapter, the RK techniques ! t nicely into the

general algorithm embodied in Fig. 25.7. Figure 25.17 presents pseudocode to determine

the slope of the classic fourth-order RK method [Eq. (25.40)]. Subroutines to compute

slopes for all the other versions can be easily programmed in a similar fashion.

25.4 SYSTEMS OF EQUATIONS

Many practical problems in engineering and science require the solution of a system of

simultaneous ordinary differential equations rather than a single equation. Such systems

may be represented generally as

dy1

dx 5 f1(x, y1, y2, p , yn)

dy2

dx 5 f2(x, y1, y2, p , yn)

#

#

#

dyn

dx 5 fn(x, y1, y2, p , yn) (25.42)

The solution of such a system requires that n initial conditions be known at the starting

value of x.

SUB RK4 (x, y, h, ynew)

CALL Derivs(x, y, k1)

ym 5 y 1 k1 ? hy2 CALL Derivs(x 1 hy2, ym, k2) ym 5 y 1 k2 ? hy2 CALL Derivs(x 1 hy2, ym, k3) ye 5 y 1 k3 ? h

CALL Derivs(x 1 h, ye, k4)

slope 5 (k1 1 2(k2 1 k3) 1 k4)y6 ynew 5 y 1 slope ? h

x 5 x 1 h

END SUB

FIGURE 25.17 Pseudocode to determine a single step of the fourth-order RK method.

740 RUNGE-KUTTA METHODS

25.4.1 Euler’s Method

All the methods discussed in this chapter for single equations can be extended to the

system shown above. Engineering applications can involve thousands of simultaneous

equations. In each case, the procedure for solving a system of equations simply involves

applying the one-step technique for every equation at each step before proceeding to the

next step. This is best illustrated by the following example for the simple Euler’s method.

EXAMPLE 25.9 Solving Systems of ODEs Using Euler’s Method

Problem Statement. Solve the following set of differential equations using Euler’s method, assuming that at x 5 0, y1 5 4, and y2 5 6. Integrate to x 5 2 with a step size

of 0.5.

dy1

dx 5 20.5y1

dy2

dx 5 4 2 0.3y2 2 0.1y1

Solution. Euler’s method is implemented for each variable as in Eq. (25.2):

y1(0.5) 5 4 1 [20.5(4)]0.5 5 3

y2(0.5) 5 6 1 [4 2 0.3(6) 2 0.1(4)]0.5 5 6.9

Note that y1(0) 5 4 is used in the second equation rather than the y1(0.5) 5 3 computed

with the ! rst equation. Proceeding in a like manner gives

x y1 y2

0 4 6 0.5 3 6.9 1.0 2.25 7.715 1.5 1.6875 8.44525 2.0 1.265625 9.094087

25.4.2 Runge-Kutta Methods

Note that any of the higher-order RK methods in this chapter can be applied to systems

of equations. However, care must be taken in determining the slopes. Figure 25.15 is help-

ful in visualizing the proper way to do this for the fourth-order method. That is, we ! rst

develop slopes for all variables at the initial value. These slopes (a set of k1’s) are then

used to make predictions of the dependent variable at the midpoint of the interval. These

midpoint values are in turn used to compute a set of slopes at the midpoint (the k2’s). These

new slopes are then taken back to the starting point to make another set of midpoint pre-

dictions that lead to new slope predictions at the midpoint (the k3’s). These are then em-

ployed to make predictions at the end of the interval that are used to develop slopes at the

end of the interval (the k4’s). Finally, the k’s are combined into a set of increment functions

[as in Eq. (25.40)] and brought back to the beginning to make the ! nal prediction. The

following example illustrates the approach.

25.4 SYSTEMS OF EQUATIONS 741

EXAMPLE 25.10 Solving Systems of ODEs Using the Fourth-Order RK Method

Problem Statement. Use the fourth-order RK method to solve the ODEs from Ex- ample 25.9.

Solution. First, we must solve for all the slopes at the beginning of the interval:

k1,1 5 f1(0, 4, 6) 5 20.5(4) 5 22

k1, 2 5 f2(0, 4, 6) 5 4 2 0.3(6) 2 0.1(4) 5 1.8

where ki, j is the ith value of k for the jth dependent variable. Next, we must calculate

the ! rst values of y1 and y2 at the midpoint:

y1 1 k1,1 h

2 5 4 1 (22)

0.5

2 5 3.5

y2 1 k1, 2 h

2 5 6 1 (1.8)

0.5

2 5 6.45

which can be used to compute the ! rst set of midpoint slopes,

k2, 1 5 f1(0.25, 3.5, 6.45) 5 21.75

k2, 2 5 f2(0.25, 3.5, 6.45) 5 1.715

These are used to determine the second set of midpoint predictions,

y1 1 k2,1 h

2 5 4 1 (21.75)

0.5

2 5 3.5625

y2 1 k2, 2 h

2 5 6 1 (1.715)

0.5

2 5 6.42875

which can be used to compute the second set of midpoint slopes,

k3, 1 5 f1(0.25, 3.5625, 6.42875) 5 21.78125

k3, 2 5 f2(0.25, 3.5625, 6.42875) 5 1.715125

These are used to determine the predictions at the end of the interval

y1 1 k3,1h 5 4 1 (21.78125)(0.5) 5 3.109375

y2 1 k3, 2h 5 6 1 (1.715125)(0.5) 5 6.857563

which can be used to compute the endpoint slopes,

k4,1 5 f1(0.5, 3.109375, 6.857563) 5 21.554688

k4, 2 5 f2(0.5, 3.109375, 6.857563) 5 1.631794

The values of k can then be used to compute [Eq. (25.40)]:

y1(0.5) 5 4 1 1

6 [22 1 2(21.75 2 1.78125) 2 1.554688]0.5 5 3.115234

y2(0.5) 5 6 1 1

6 [1.8 1 2(1.715 1 1.715125) 1 1.631794]0.5 5 6.857670

742 RUNGE-KUTTA METHODS

Proceeding in a like manner for the remaining steps yields

x y1 y2

0 4 6 0.5 3.115234 6.857670 1.0 2.426171 7.632106 1.5 1.889523 8.326886 2.0 1.471577 8.946865

25.4.3 Computer Algorithm for Solving Systems of ODEs

The computer code for solving a single ODE with Euler’s method (Fig. 25.7) can be

easily extended to systems of equations. The modi! cations include:

1. Inputting the number of equations, n.

2. Inputting the initial values for each of the n dependent variables.

3. Modifying the algorithm so that it computes slopes for each of the dependent

variables.

4. Including additional equations to compute derivative values for each of the ODEs.

5. Including loops to compute a new value for each dependent variable.

Such an algorithm is outlined in Fig. 25.18 for the fourth-order RK method. Notice

how similar it is in structure and organization to Fig. 25.7. Most of the differences relate

to the fact that

1. There are n equations.

2. The added detail of the fourth-order RK method.

EXAMPLE 25.11 Solving Systems of ODEs with the Computer

Problem Statement. A computer program to implement the fourth-order RK method for systems can be easily developed based on Fig. 25.18. Such software makes it con-

venient to compare different models of a physical system. For example, a linear model

for a swinging pendulum is given by [recall Eq. (PT7.11)]

dy1

dx 5 y2

dy2

dx 5 216.1y1

where y1 and y2 5 angular displacement and velocity. A nonlinear model of the same

system is [recall Eq. (PT7.9)]

dy3

dx 5 y4

dy4

dx 5 216.1 sin(y3)

where y3 and y4 5 angular displacement and velocity for the nonlinear case. Solve these

systems for two cases: (a) a small initial displacement (y1 5 y3 5 0.1 radians; y2 5 y4 5 0)

and (b) a large displacement (y1 5 y3 5 py4 5 0.785398 radians; y2 5 y4 5 0).

25.4 SYSTEMS OF EQUATIONS 743

(a) Main or “Driver” Program

Assign values for

n 5 number of equations

yi 5 initial values of n dependent

variables

xi 5 initial value independent

variable

xf 5 final value independent variable

dx 5 calculation step size

xout 5 output interval

x 5 xi

m 5 0

xpm 5 x

DOFOR i 5 1, n

ypi,m 5 yii

yi 5 yii

END DO

DO

xend 5 x 1 xout

IF (xend . xf) THEN xend 5 xf

h 5 dx

CALL Integrator (x, y, n, h, xend)

m 5 m 1 1

xpm 5 x

DOFOR i 5 1, n

ypi,m 5 yi

END DO

IF (x $ xf) EXIT

END DO

DISPLAY RESULTS

END

(b) Routine to Take One Output Step

SUB Integrator (x, y, n, h, xend)

DO

IF (xend 2 x , h) THEN h 5 xend 2 x

CALL RK4 (x, y, n, h)

IF (x $ xend) EXIT

END DO

END SUB

(c) Fourth-Order RK Method for a System of ODEs

SUB RK4 (x, y, n, h)

CALL Derivs (x, y, k1)

DOFOR i 5 1, n

ymi 5 yi 1 k1i * h / 2

END DO

CALL Derivs (x 1 h / 2, ym, k2)

DOFOR i 5 1, n

ymi 5 yi 1 k2i * h / 2

END DO

CALL Derivs (x 1 h / 2, ym, k3)

DOFOR i 5 1, n

yei 5 yi 1 k3i * h

END DO

CALL Derivs (x 1 h, ye, k4)

DOFOR i 5 1, n

slopei 5 (k1i 1 2*(k2i1k3i)1k4i)/6

yi 5 yi 1 slopei * h

END DO

x 5 x 1 h

END SUB

(d ) Routine to Determine Derivatives

SUB Derivs (x, y, dy)

dy1 5 …

dy2 5 …

END SUB

FIGURE 25.18 Pseudocode for the fourth-order RK method for systems.

744 RUNGE-KUTTA METHODS

Solution.

(a) The calculated results for the linear and nonlinear models are almost identical

(Fig. 25.19a). This is as expected because when the initial displacement is small,

sin (u) > u. (b) When the initial displacement is py4 5 0.785398, the solutions are much different

and the difference is magni! ed as time becomes larger and larger (Fig. 25.19b). This

is expected because the assumption that sin (u) 5 u is poor when theta is large.

25.5 ADAPTIVE RUNGE-KUTTA METHODS

To this point, we have presented methods for solving ODEs that employ a constant step

size. For a signi! cant number of problems, this can represent a serious limitation. For

example, suppose that we are integrating an ODE with a solution of the type depicted

in Fig. 25.20. For most of the range, the solution changes gradually. Such behavior sug-

gests that a fairly large step size could be employed to obtain adequate results. However,

for a localized region from x 5 1.75 to x 5 2.25, the solution undergoes an abrupt change.

The practical consequence of dealing with such functions is that a very small step size

would be required to accurately capture the impulsive behavior. If a constant step-size al-

gorithm were employed, the smaller step size required for the region of abrupt change would

have to be applied to the entire computation. As a consequence, a much smaller step size

than necessary—and, therefore, many more calculations—would be wasted on the regions

of gradual change.

4

2

0y

0 321

x

– 4

– 2

4

y1, y3

y2, y4

(a)

4

2

0y

0 2 31

x

– 4

– 2

4

y2 y4

y3

y1

(b)

FIGURE 25.19 Solutions obtained with a computer program for the fourth-order RK method. The plots represent solutions for both linear and nonlinear pendulums with (a) small and (b) large initial displacements.

25.5 ADAPTIVE RUNGE-KUTTA METHODS 745

Algorithms that automatically adjust the step size can avoid such overkill and hence

be of great advantage. Because they “adapt” to the solution’s trajectory, they are said to

have adaptive step-size control. Implementation of such approaches requires that an es-

timate of the local truncation error be obtained at each step. This error estimate can then

serve as a basis for either lengthening or decreasing the step size.

Before proceeding, we should mention that aside from solving ODEs, the methods

described in this chapter can also be used to evaluate de! nite integrals. As mentioned

previously in the introduction to Part Six, the evaluation of the integral

I 5 # b

a f(x) dx

is equivalent to solving the differential equation

dy

dx 5 f(x)

for y(b) given the initial condition y(a) 5 0. Thus, the following techniques can be em-

ployed to ef! ciently evaluate de! nite integrals involving functions that are generally

smooth but exhibit regions of abrupt change.

There are two primary approaches to incorporate adaptive step-size control into one-

step methods. In the ! rst, the error is estimated as the difference between two predictions

using the same-order RK method but with different step sizes. In the second, the local

FIGURE 25.20 An example of a solution of an ODE that exhibits an abrupt change. Automatic step-size adjustment has great advantages for such cases.

1

0 1 2 3

y

x