# Engineering

PT7.3 ORIENTATION 707

to develop two improved versions of Euler’s method—the Heun and the midpoint tech-

niques. After this introduction, we formally develop the concept of Runge-Kutta (or RK)

approaches and demonstrate how the foregoing techniques are actually rst- and second-

order RK methods. This is followed by a discussion of the higher-order RK formulations

that are frequently used for engineering problem solving. In addition, we cover the ap-

plication of one-step methods to systems of ODEs. Finally, the chapter ends with a

discussion of adaptive RK methods that automatically adjust the step size in response to

the truncation error of the computation.

Chapter 26 starts with a description of stiff ODEs. These are both individual and

systems of ODEs that have both fast and slow components to their solution. We intro-

duce the idea of an implicit solution technique as one commonly used remedy for this

problem.

Next, we discuss multistep methods. These algorithms retain information of previous

steps to more effectively capture the trajectory of the solution. They also yield the trunca-

tion error estimates that can be used to implement step-size control. In this section, we

initially take a visual, intuitive approach by using a simple method—the non-self-starting

Heun—to introduce all the essential features of the multistep approaches.

In Chap. 27 we turn to boundary-value and eigenvalue problems. For the former,

we introduce both shooting and ! nite-difference methods. For the latter, we discuss sev-

eral approaches, including the polynomial and the power methods. Finally, the chapter

concludes with a description of the application of several software packages and librar-

ies for solution of ODEs and eigenvalues.

Chapter 28 is devoted to applications from all the elds of engineering. Finally, a

short review section is included at the end of Part Seven. This epilogue summarizes and

compares the important formulas and concepts related to ODEs. The comparison includes

a discussion of trade-offs that are relevant to their implementation in engineering prac-

tice. The epilogue also summarizes important formulas and includes references for

advanced topics.

PT7.3.2 Goals and Objectives

Study Objectives. After completing Part Seven, you should have greatly enhanced your capability to confront and solve ordinary differential equations and eigenvalue prob-

lems. General study goals should include mastering the techniques, having the capability

to assess the reliability of the answers, and being able to choose the “best’’ method (or

methods) for any particular problem. In addition to these general objectives, the speci c

study objectives in Table PT7.2 should be mastered.

Computer Objectives. Algorithms are provided for many of the methods in Part Seven. This information will allow you to expand your software library. For example,

you may nd it useful from a professional viewpoint to have software that employs the

fourth-order Runge-Kutta method for more than ve equations and to solve ODEs with

an adaptive step-size approach.

In addition, one of your most important goals should be to master several of the

general-purpose software packages that are widely available. In particular, you should

become adept at using these tools to implement numerical methods for engineering

problem solving.

708 ORDINARY DIFFERENTIAL EQUATIONS

TABLE PT7.2 Specifi c study objectives for Part Seven.

1. Understand the visual representations of Euler’s, Heun’s, and the midpoint methods 2. Know the relationship of Euler’s method to the Taylor series expansion and the insight it provides

regarding the error of the method 3. Understand the difference between local and global truncation errors and how they relate to the

choice of a numerical method for a particular problem 4. Know the order and the step-size dependency of the global truncation errors for all the methods

described in Part Seven; understand how these errors bear on the accuracy of the techniques 5. Understand the basis of predictor-corrector methods; in particular, realize that the effi ciency of the

corrector is highly dependent on the accuracy of the predictor 6. Know the general form of the Runge-Kutta methods; understand the derivation of the second-order

RK method and how it relates to the Taylor series expansion; realize that there are an infi nite number of possible versions for second- and higher-order RK methods

7. Know how to apply any of the RK methods to systems of equations; be able to reduce an nth-order ODE to a system of n fi rst-order ODEs

8. Recognize the type of problem context where step size adjustment is important 9. Understand how adaptive step size control is integrated into a fourth-order RK method 10. Recognize how the combination of slow and fast components makes an equation or a system of

equations stiff 11. Understand the distinction between explicit and implicit solution schemes for ODEs; in particular,

recognize how the latter (1) ameliorates the stiffness problem and (2) complicates the solution mechanics

12. Understand the difference between initial-value and boundary-value problems 13. Know the difference between multistep and one-step methods; realize that all multistep methods are

predictor-correctors but that not all predictor-correctors are multistep methods 14. Understand the connection between integration formulas and predictor-corrector methods 15. Recognize the fundamental difference between Newton-Cotes and Adams integration formulas 16. Know the rationale behind the polynomial and the power methods for determining eigenvalues; in

particular, recognize their strengths and limitations 17. Understand how Hoteller’s defl ation allows the power method to be used to compute intermediate

eigenvalues 18. Know how to use software packages and/or libraries to integrate ODEs and evaluate eigenvalues

25

709

C H A P T E R 25

Runge-Kutta Methods

This chapter is devoted to solving ordinary differential equations of the form

dy

dx 5 f(x, y)

In Chap. 1, we used a numerical method to solve such an equation for the velocity of

the falling parachutist. Recall that the method was of the general form

New value 5 old value 1 slope 3 step size

or, in mathematical terms,

yi11 5 yi 1 fh (25.1)

According to this equation, the slope estimate of f is used to extrapolate from an old value

yi to a new value yi 1 1 over a distance h (Fig. 25.1). This formula can be applied step by

step to compute out into the future and, hence, trace out the trajectory of the solution.

FIGURE 25.1 Graphical depiction of a one- step method.

y

x

Step size = h

Slope =

xi xi + 1

yi + 1 = yi + h

710 RUNGE-KUTTA METHODS

All one-step methods can be expressed in this general form, with the only difference

being the manner in which the slope is estimated. As in the falling parachutist problem,

the simplest approach is to use the differential equation to estimate the slope in the form

of the ! rst derivative at xi. In other words, the slope at the beginning of the interval is

taken as an approximation of the average slope over the whole interval. This approach,

called Euler’s method, is discussed in the ! rst part of this chapter. This is followed by

other one-step methods that employ alternative slope estimates that result in more ac-

curate predictions. All these techniques are generally called Runge-Kutta methods.

25.1 EULER’S METHOD

The ! rst derivative provides a direct estimate of the slope at xi (Fig. 25.2):

f 5 f(xi, yi)

where f(xi, yi) is the differential equation evaluated at xi and yi. This estimate can be

substituted into Eq. (25.1):

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

This formula is referred to as Euler’s (or the Euler-Cauchy or the point-slope)

method. A new value of y is predicted using the slope (equal to the ! rst derivative at the

original value of x) to extrapolate linearly over the step size h (Fig. 25.2).

EXAMPLE 25.1 Euler’s Method

Problem Statement. Use Euler’s method to numerically integrate Eq. (PT7.13):

dy

dx 5 22×3 1 12×2 2 20x 1 8.5

y

xxi + 1

error

Predicted

True

xi

h

FIGURE 25.2 Euler’s method.

25.1 EULER’S METHOD 711

from x 5 0 to x 5 4 with a step size of 0.5. The initial condition at x 5 0 is y 5 1.

Recall that the exact solution is given by Eq. (PT7.16):

y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1

Solution. Equation (25.2) can be used to implement Euler’s method:

y(0.5) 5 y(0) 1 f(0, 1)0.5

where y(0) 5 1 and the slope estimate at x 5 0 is

f(0, 1) 5 22(0)3 1 12(0)2 2 20(0) 1 8.5 5 8.5

Therefore,

y(0.5) 5 1.0 1 8.5(0.5) 5 5.25

The true solution at x 5 0.5 is

y 5 20.5(0.5)4 1 4(0.5)3 2 10(0.5)2 1 8.5(0.5) 1 1 5 3.21875

Thus, the error is

Et 5 true 2 approximate 5 3.21875 2 5.25 5 22.03125

or, expressed as percent relative error, et 5 263.1%. For the second step,

y(1) 5 y(0.5) 1 f(0.5, 5.25)0.5

5 5.25 1 [22(0.5)3 1 12(0.5)2 2 20(0.5) 1 8.5]0.5

5 5.875

The true solution at x 5 1.0 is 3.0, and therefore, the percent relative error is 295.8%.

The computation is repeated, and the results are compiled in Table 25.1 and Fig. 25.3.

TABLE 25.1 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 Euler’s method with a step size of 0.5. The local error refers to the error incurred over a single step. It is calculated with a Taylor series expansion as in Example 25.2. The global error is the total discrepancy due to past as well as present steps.

Percent Relative Error

x ytrue yEuler Global Local

0.0 1.00000 1.00000 0.5 3.21875 5.25000 263.1 263.1 1.0 3.00000 5.87500 295.8 228.1 1.5 2.21875 5.12500 2131.0 21.4 2.0 2.00000 4.50000 2125.0 20.3 2.5 2.71875 4.75000 274.7 17.2 3.0 4.00000 5.87500 246.9 3.9 3.5 4.71875 7.12500 251.0 211.3 4.0 3.00000 7.00000 2133.3 253.1

712 RUNGE-KUTTA METHODS

Note that, although the computation captures the general trend of the true solution, the

error is considerable. As discussed in the next section, this error can be reduced by using

a smaller step size.

The preceding example uses a simple polynomial for the differential equation to

facilitate the error analyses that follow. Thus,

dy

dx 5 f(x)

Obviously, a more general (and more common) case involves ODEs that depend on both

x and y,

dy

dx 5 f(x, y)

As we progress through this part of the text, our examples will increasingly involve ODEs

that depend on both the independent and the dependent variables.

25.1.1 Error Analysis for Euler’s Method

The numerical solution of ODEs involves two types of error (recall Chaps. 3 and 4):

1. Truncation, or discretization, errors caused by the nature of the techniques employed

to approximate values of y.