# Engineering

FIGURE 25.7 Pseudocode for an “improved” modular version of Euler’s method.

(a) Main or “Driver” Program

Assign values for

y 5 initial value dependent variable

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

ypm 5 y

DO

xend 5 x 1 xout

IF (xend . xf) THEN xend 5 xf

h 5 dx

CALL Integrator (x, y, h, xend)

m 5 m 1 1

xpm 5 x

ypm 5 y

IF (x $ xf) EXIT

END DO

DISPLAY RESULTS

END

(b) Routine to Take One Output Step

SUB Integrator (x, y, h, xend)

DO

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

CALL Euler (x, y, h, ynew)

y 5 ynew

IF (x $ xend) EXIT

END DO

END SUB

(c) Euler’s Method for a Single ODE

SUB Euler (x, y, h, ynew)

CALL Derivs (x, y, dydx)

ynew 5 y 1 dydx * h

x 5 x 1 h

END SUB

(d) Routine to Determine Derivative

SUB Derivs (x, y, dydx)

dydx 5 . . .

END SUB

720 RUNGE-KUTTA METHODS

EXAMPLE 25.4 Solving ODEs with the Computer

Problem Statement. A computer program can be developed from the pseudocode in Fig. 25.7. We can use this software to solve another problem associated with the falling

parachutist. You recall from Part One that our mathematical model for the velocity was

based on Newton’s second law in the form

dy

dt 5 g 2

c

m y (E25.4.1)

This differential equation was solved both analytically (Example 1.1) and numerically

using Euler’s method (Example 1.2). These solutions were for the case where g 5 9.81,

c 5 12.5, m 5 68.1, and y 5 0 at t 5 0.

The objective of the present example is to repeat these numerical computations

employing a more complicated model for the velocity based on a more complete math-

ematical description of the drag force caused by wind resistance. This model is given by

dy

dt 5 g 2

c

m cy 1 a a y

ymax bb d (E25.4.2)

where g, m, and c are the same as for Eq. (E25.4.1), and a, b, and ymax are empirical

constants, which for this case are equal to 8.3, 2.2, and 46, respectively. Note that this

model is more capable of accurately ! tting empirical measurements of drag forces versus

velocity than is the simple linear model of Example 1.1. However, this increased ” exibil-

ity is gained at the expense of evaluating three coef! cients rather than one. Furthermore,

the resulting mathematical model is more dif! cult to solve analytically. In this case,

Euler’s method provides a convenient alternative to obtain an approximate numerical

solution.

FIGURE 25.8 Graphical results for the solution of the nonlinear ODE [Eq. (E25.4.2)]. Notice that the plot also shows the solution for the linear model [Eq. (E25.4.1)] for comparative purposes.

60

y

40

20

0 t15

Nonlinear

Linear

5 100

25.2 IMPROVEMENTS OF EULER’S METHOD 721

Solution. The results for both the linear and nonlinear model are displayed in Fig. 25.8 with an integration step size of 0.1 s. The plot in Fig. 25.8 also shows an overlay of the

solution of the linear model for comparison purposes.

The results of the two simulations indicate how increasing the complexity of the for-

mulation of the drag force affects the velocity of the parachutist. In this case, the terminal

velocity is lowered because of resistance caused by the higher-order terms in Eq. (E25.4.2).

Alternative models could be tested in a similar fashion. The combination of a com-

puter-generated solution makes this an easy and ef! cient task. This convenience should

allow you to devote more of your time to considering creative alternatives and holistic

aspects of the problem rather than to tedious manual computations.

25.1.3 Higher-Order Taylor Series Methods

One way to reduce the error of Euler’s method would be to include higher-order terms

of the Taylor series expansion in the solution. For example, including the second-order

term from Eq. (25.6) yields

yi11 5 yi 1 f(xi, yi)h 1 f ¿(xi, yi)

2! h2 (25.11)

with a local truncation error of

Ea 5 f –(xi, yi)

6 h3

Although the incorporation of higher-order terms is simple enough to implement for

polynomials, their inclusion is not so trivial when the ODE is more complicated. In

particular, ODEs that are a function of both the dependent and independent variable

require chain-rule differentiation. For example, the ! rst derivative of f(x, y) is

f ¿(xi, yi) 5 0f(x, y)

0x 1 0f(x, y)

0y dy

dx

The second derivative is

f –(xi, yi) 5 0[0fy0x 1 (0fy0y)(dyydx)]

0x 1 0[0fy0x 1 (0fy0y)(dyydx)]

0y dy

dx

Higher-order derivatives become increasingly more complicated.

Consequently, as described in the following sections, alternative one-step methods

have been developed. These schemes are comparable in performance to the higher-order

Taylor-series approaches but require only the calculation of ! rst derivatives.

25.2 IMPROVEMENTS OF EULER’S METHOD

A fundamental source of error in Euler’s method is that the derivative at the beginning of

the interval is assumed to apply across the entire interval. Two simple modi! cations are

available to help circumvent this shortcoming. As will be demonstrated in Sec. 25.3, both

modi! cations actually belong to a larger class of solution techniques called Runge-Kutta

722 RUNGE-KUTTA METHODS

methods. However, because they have a very straightforward graphical interpretation, we

will present them prior to their formal derivation as Runge-Kutta methods.

25.2.1 Heun’s Method

One method to improve the estimate of the slope involves the determination of two

derivatives for the interval—one at the initial point and another at the end point. The two

derivatives are then averaged to obtain an improved estimate of the slope for the entire

interval. This approach, called Heun’s method, is depicted graphically in Fig. 25.9.

Recall that in Euler’s method, the slope at the beginning of an interval

y¿i 5 f(xi, yi) (25.12)

is used to extrapolate linearly to yi11:

y0i11 5 yi 1 f(xi, yi)h (25.13)

For the standard Euler method we would stop at this point. However, in Heun’s method

the y0i11 calculated in Eq. (25.13) is not the ! nal answer, but an intermediate prediction.

This is why we have distinguished it with a superscript 0. Equation (25.13) is called a

FIGURE 25.9 Graphical depiction of Heun’s method. (a) Predictor and (b) corrector.

y

xxi + 1xi

(a)

Slope = f (xi, yi)

Slope = f (xi + 1, y 0 i + 1)

y

xxi + 1xi

(b)

Slope = f (xi, yi) + f (xi + 1, yi + 1)

2

0

25.2 IMPROVEMENTS OF EULER’S METHOD 723

predictor equation. It provides an estimate of yi11 that allows the calculation of an esti-

mated slope at the end of the interval:

y¿i11 5 f(xi11, y 0 i11) (25.14)

Thus, the two slopes [Eqs. (25.12) and (25.14)] can be combined to obtain an average

slope for the interval:

y¿ 5 y¿i 1 y¿i11

2 5

f(xi, yi) 1 f(xi11, y 0 i11)

2

This average slope is then used to extrapolate linearly from yi to yi 1 l using Euler’s method:

yi11 5 yi 1 f(xi, yi) 1 f(xi11, y

0 i11)

2 h

which is called a corrector equation.

The Heun method is a predictor-corrector approach. All the multistep methods to

be discussed subsequently in Chap. 26 are of this type. The Heun method is the only

one-step predictor-corrector method described in this book. As derived above, it can be

expressed concisely as

Predictor (Fig. 25.9a): y0i11 5 yi 1 f(xi, yi)h (25.15)

Corrector (Fig. 25.9b): yi11 5 yi 1 f(xi, yi) 1 f(xi11, y

0 i11)

2 h (25.16)

Note that because Eq. (25.16) has yi1l on both sides of the equal sign, it can be applied

in an iterative fashion. That is, an old estimate can be used repeatedly to provide an im-

proved estimate of yi1l. The process is depicted in Fig. 25.10. It should be understood that

FIGURE 25.10 Graphical representation of iterating the corrector of Heun’s method to obtain an improved estimate.

f(x i , y

i ) + f (x

i+ 1 , y

i+ 1 )

y i + h

y i +

1

2

0

724 RUNGE-KUTTA METHODS

this iterative process does not necessarily converge on the true answer but will converge

on an estimate with a ! nite truncation error, as demonstrated in the following example.

As with similar iterative methods discussed in previous sections of the book, a ter-

mination criterion for convergence of the corrector is provided by [recall Eq. (3.5)]

Zea Z 5 ` y j i11 2 y

j21 i11

y j i11

` 100% (25.17) where y

j21 i11 and y

j i11 are the result from the prior and the present iteration of the correc-

tor, respectively.

EXAMPLE 25.5 Heun’s Method

Problem Statement. Use Heun’s method to integrate y9 5 4e0.8x 2 0.5y from x 5 0 to x 5 4 with a step size of 1. The initial condition at x 5 0 is y 5 2.

Solution. Before solving the problem numerically, we can use calculus to determine the following analytical solution:

y 5 4

1.3 (e0.8x 2 e20.5x) 1 2e20.5x (E25.5.1)

This formula can be used to generate the true solution values in Table 25.2.

First, the slope at (x0, y0) is calculated as

y¿0 5 4e 0

2 0.5(2) 5 3

This result is quite different from the actual average slope for the interval from 0 to 1.0,

which is equal to 4.1946, as calculated from the differential equation using Eq. (PT6.4).

The numerical solution is obtained by using the predictor [Eq. (25.15)] to obtain an

estimate of y at 1.0:

y01 5 2 1 3(1) 5 5

TABLE 25.2 Comparison of true and approximate values of the integral of y9 5 4e0.8x 2 0.5y, with the initial condition that y 5 2 at x 5 0. The approximate values were computed using the Heun method with a step size of 1. Two cases, corresponding to different numbers of corrector iterations, are shown, along with the absolute percent relative error.

Iterations of Heun’s Method

1 15

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

0 2.0000000 2.0000000 0.00 2.0000000 0.00 1 6.1946314 6.7010819 8.18 6.3608655 2.68 2 14.8439219 16.3197819 9.94 15.3022367 3.09 3 33.6771718 37.1992489 10.46 34.7432761 3.17 4 75.3389626 83.3377674 10.62 77.7350962 3.18

25.2 IMPROVEMENTS OF EULER’S METHOD 725

Note that this is the result that would be obtained by the standard Euler method. The true

value in Table 25.2 shows that it corresponds to a percent relative error of 19.3 percent.

Now, to improve the estimate for yi11, we use the value y 0 1 to predict the slope at

the end of the interval

y¿1 5 f(x1, y 0 1) 5 4e

0.8(1) 2 0.5(5) 5 6.402164

which can be combined with the initial slope to yield an average slope over the interval

from x 5 0 to 1

y¿ 5 3 1 6.402164

2 5 4.701082

which is closer to the true average slope of 4.1946. This result can then be substituted

into the corrector [Eq. (25.16)] to give the prediction at x 5 1

y1 5 2 1 4.701082(1) 5 6.701082

which represents a percent relative error of 28.18 percent. Thus, the Heun method with-

out iteration of the corrector reduces the absolute value of the error by a factor of 2.4

as compared with Euler’s method.

Now this estimate can be used to re! ne or correct the prediction of y1 by substitut-

ing the new result back into the right-hand side of Eq. (25.16):

y1 5 2 1 [3 1 4e0.8(1) 2 0.5(6.701082) ]

2 1 5 6.275811

which represents an absolute percent relative error of 1.31 percent. This result, in turn,

can be substituted back into Eq. (25.16) to further correct:

y1 5 2 1 [3 1 4e0.8(1) 2 0.5(6.275811) ]

2 1 5 6.382129

which represents an Zet Z of 3.03%. Notice how the errors sometimes grow as the iterations proceed. Such increases can occur, especially for large step sizes, and they prevent us

from drawing the general conclusion that an additional iteration will always improve the

result. However, for a suf! ciently small step size, the iterations should eventually con-

verge on a single value. For our case, 6.360865, which represents a relative error of 2.68

percent, is attained after 15 iterations. Table 25.2 shows results for the remainder of the

computation using the method with 1 and 15 iterations per step.

In the previous example, the derivative is a function of both the dependent variable

y and the independent variable x. For cases such as polynomials, where the ODE is solely

a function of the independent variable, the predictor step [Eq. (25.16)] is not required

and the corrector is applied only once for each iteration. For such cases, the technique

is expressed concisely as

yi11 5 yi 1 f(xi) 1 f(xi11)

2 h (25.18)