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
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
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
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
C H A P T E R 25
This chapter is devoted to solving ordinary differential equations of the form
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.
Step size = h
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):
dx 5 22×3 1 12×2 2 20x 1 8.5
xxi + 1
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
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
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,
dx 5 f(x)
Obviously, a more general (and more common) case involves ODEs that depend on both
x and y,
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.