Numerical Methods for Engineers
In the rst chapter of this book, we derived the following equation based on Newton’s
second law to compute the velocity y of a falling parachutist as a function of time t
[recall Eq. (1.9)]:
dt 5 g 2
m y (PT7.1)
where g is the gravitational constant, m is the mass, and c is a drag coef cient. Such
equations, which are composed of an unknown function and its derivatives, are called
differential equations. Equation (PT7.1) is sometimes referred to as a rate equation
because it expresses the rate of change of a variable as a function of variables and pa-
rameters. Such equations play a fundamental role in engineering because many physical
phenomena are best formulated mathematically in terms of their rate of change.
In Eq. (PT7.1), the quantity being differentiated, y, is called the dependent variable.
The quantity with respect to which y is differentiated, t, is called the independent vari-
able. When the function involves one independent variable, the equation is called an
ordinary differential equation (or ODE). This is in contrast to a partial differential equa-
tion (or PDE) that involves two or more independent variables.
Differential equations are also classi ed as to their order. For example, Eq. (PT7.1)
is called a ! rst-order equation because the highest derivative is a rst derivative. A
second-order equation would include a second derivative. For example, the equation
describing the position x of a mass-spring system with damping is the second-order
m d 2x
dt 2 1 c
dt 1 kx 5 0 (PT7.2)
where c is a damping coef cient and k is a spring constant. Similarly, an nth-order equa-
tion would include an nth derivative.
Higher-order equations can be reduced to a system of rst-order equations. For Eq.
(PT7.2), this is done by de ning a new variable y, where
y 5 dx
which itself can be differentiated to yield
ORDINARY DIFFERENTIAL EQUATIONS
700 ORDINARY DIFFERENTIAL EQUATIONS
Equations (PT7.3) and (PT7.4) can then be substituted into Eq. (PT7.2) to give
dt 1 cy 1 kx 5 0 (PT7.5)
dt 5 2
cy 1 kx
Thus, Eqs. (PT7.3) and (PT7.6) are a pair of rst-order equations that are equivalent to
the original second-order equation. Because other nth-order differential equations can be
similarly reduced, this part of our book focuses on the solution of rst-order equations.
Some of the engineering applications in Chap. 28 deal with the solution of second-order
ODEs by reduction to a pair of rst-order equations.
PT7.1.1 Noncomputer Methods for Solving ODEs
Without computers, ODEs are usually solved with analytical integration techniques. For
example, Eq. (PT7.1) could be multiplied by dt and integrated to yield
y 5 # ag 2 cm yb dt (PT7.7) The right-hand side of this equation is called an inde! nite integral because the limits of
integration are unspeci ed. This is in contrast to the de nite integrals discussed previously
in Part Six [compare Eq. (PT7.7) with Eq. (PT6.6)].
An analytical solution for Eq. (PT7.7) is obtained if the inde nite integral can be
evaluated exactly in equation form. For example, recall that for the falling parachutist
problem, Eq. (PT7.7) was solved analytically by Eq. (1.10) (assuming y 5 0 at t 5 0):
y(t) 5 gm
c (1 2 e2(cym)t) (1.10)
The mechanics of deriving such analytical solutions will be discussed in Sec. PT7.2. For
the time being, the important fact is that exact solutions for many ODEs of practical
importance are not available. As is true for most situations discussed in other parts of
this book, numerical methods offer the only viable alternative for these cases. Because
these numerical methods usually require computers, engineers in the precomputer era
were somewhat limited in the scope of their investigations.
One very important method that engineers and applied mathematicians developed to
overcome this dilemma was linearization. A linear ordinary differential equation is one
that ts the general form
1 p 1 a1(x)y¿ 1 a0(x)y 5 f(x) (PT7.8)
where y(n) is the nth derivative of y with respect to x and the a’s and f ’s are speci ed
functions of x. This equation is called linear because there are no products or nonlinear
functions of the dependent variable y and its derivatives. The practical importance of
linear ODEs is that they can be solved analytically. In contrast, most nonlinear equations
PT7.1 MOTIVATION 701
cannot be solved exactly. Thus, in the precomputer era, one tactic for solving nonlinear
equations was to linearize them.
A simple example is the application of ODEs to predict the motion of a swinging
pendulum (Fig. PT7.1). In a manner similar to the derivation of the falling parachutist
problem, Newton’s second law can be used to develop the following differential equation
(see Sec. 28.4 for the complete derivation):
dt 2 1
l sin u 5 0 (PT7.9)
where u is the angle of displacement of the pendulum, g is the gravitational constant,
and l is the pendulum length. This equation is nonlinear because of the term sin u. One
way to obtain an analytical solution is to realize that for small displacements of the
pendulum from equilibrium (that is, for small values of u),
sin u > u (PT7.10)
Thus, if it is assumed that we are interested only in cases where u is small, Eq. (PT7.10)
can be substituted into Eq. (PT7.9) to give
dt 2 1
l u 5 0 (PT7.11)
We have, therefore, transformed Eq. (PT7.9) into a linear form that is easy to solve
Although linearization remains a very valuable tool for engineering problem solving,
there are cases where it cannot be invoked. For example, suppose that we were interested
in studying the behavior of the pendulum for large displacements from equilibrium. In
such instances, numerical methods offer a viable option for obtaining solutions. Today,
the widespread availability of computers places this option within reach of all practicing
PT7.1.2 ODEs and Engineering Practice
The fundamental laws of physics, mechanics, electricity, and thermodynamics are usually
based on empirical observations that explain variations in physical properties and states
of systems. Rather than describing the state of physical systems directly, the laws are
usually couched in terms of spatial and temporal changes.
Several examples are listed in Table PT7.1. These laws de ne mechanisms of change.
When combined with continuity laws for energy, mass, or momentum, differential equa-
tions result. Subsequent integration of these differential equations results in mathematical
functions that describe the spatial and temporal state of a system in terms of energy,
mass, or velocity variations.
The falling parachutist problem introduced in Chap. 1 is an example of the derivation
of an ordinary differential equation from a fundamental law. Recall that Newton’s second
law was used to develop an ODE describing the rate of change of velocity of a falling
parachutist. By integrating this relationship, we obtained an equation to predict fall veloc-
ity as a function of time (Fig. PT7.2). This equation could be utilized in a number of
different ways, including design purposes.
FIGURE PT7.1 The swinging pedulum.
702 ORDINARY DIFFERENTIAL EQUATIONS
In fact, such mathematical relationships are the basis of the solution for a great
number of engineering problems. However, as described in the previous section, many
of the differential equations of practical signi cance cannot be solved using the analyti-
cal methods of calculus. Thus, the methods discussed in the following chapters are
extremely important in all elds of engineering.
TABLE PT7.1 Examples of fundamental laws that are written in terms of the rate of change of variables (t 5 time and x 5 position).
Law Mathematical Expression Variables and Parameters
Newton’s second law Velocity (v), force (F), and of motion mass (m)
Fourier’s heat law Heat fl ux (q), thermal conductivity (k9) and temperature (T)
Fick’s law of diffusion Mass fl ux ( J), diffusion coeffi cient (D), and concentration (c)
Faraday’s law Voltage drop (DV L), inductance (L), (voltage drop across and current (i) an inductor)
q 5 2k¿ dT
J 5 2D dc
¢VL 5 L di
F = ma
v = (1 – e– (c/m)t) gm
c vi + 1 = vi + (g – vi) t
= g – vdv
FIGURE PT7.2 The sequence of events in the application of ODEs for engineering problem solving. The exam- ple shown is the velocity of a falling parachutist.
PT7.2 MATHEMATICAL BACKGROUND 703
PT7.2 MATHEMATICAL BACKGROUND
A solution of an ordinary differential equation is a speci c function of the independent
variable and parameters that satis es the original differential equation. To illustrate this
concept, let us start with a given function
y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1 (PT7.12)
which is a fourth-order polynomial (Fig. PT7.3a). Now, if we differentiate Eq. (PT7.12),
we obtain an ODE:
dx 5 22×3 1 12×2 2 20x 1 8.5 (PT7.13)
This equation also describes the behavior of the polynomial, but in a manner different
from Eq. (PT7.12). Rather than explicitly representing the values of y for each value of
x, Eq. (PT7.13) gives the rate of change of y with respect to x (that is, the slope) at every
value of x. Figure PT7.3 shows both the function and the derivative plotted versus x. Notice
FIGURE PT7.3 Plots of (a) y versus x and (b) dy/dx versus x for the function y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1.
704 ORDINARY DIFFERENTIAL EQUATIONS
how the zero values of the derivatives correspond to the point at which the original func-
tion is ! at—that is, has a zero slope. Also, the maximum absolute values of the derivatives
are at the ends of the interval where the slopes of the function are greatest.
Although, as just demonstrated, we can determine a differential equation given the
original function, the object here is to determine the original function given the differ-
ential equation. The original function then represents the solution. For the present case,
we can determine this solution analytically by integrating Eq. (PT7.13):
y 5 # (22×3 1 12×2 2 20x 1 8.5) dx Applying the integration rule (recall Table PT6.2)
#un du 5 u n11
n 1 1 1 C n ? 21
to each term of the equation gives the solution
y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 C (PT7.14)
which is identical to the original function with one notable exception. In the course of
differentiating and then integrating, we lost the constant value of 1 in the original equa-
tion and gained the value C. This C is called a constant of integration. The fact that such
an arbitrary constant appears indicates that the solution is not unique. In fact, it is but
one of an in nite number of possible functions (corresponding to an in nite number of
possible values of C) that satisfy the differential equation. For example, Fig. PT7.4 shows
six possible functions that satisfy Eq. (PT7.14).
FIGURE PT7.4 Six possible solutions for the integral of 22×3 1 12×2 2 20x 1 8.5. Each conforms to a different value of the constant of integration C.
x C = 0
C = – 1
C = – 2
C = 3
C = 2
C = 1
PT7.3 ORIENTATION 705
Therefore, to specify the solution completely, a differential equation is usually ac-
companied by auxiliary conditions. For rst-order ODEs, a type of auxiliary condition
called an initial value is required to determine the constant and obtain a unique solution.
For example, Eq. (PT7.13) could be accompanied by the initial condition that at x 5 0,
y 5 1. These values could be substituted into Eq. (PT7.14):
1 5 20.5(0)4 1 4(0)3 2 10(0)2 1 8.5(0) 1 C (PT7.15)
to determine C 5 1. Therefore, the unique solution that satis es both the differential
equation and the speci ed initial condition is obtained by substituting C 5 1 into Eq.
(PT7.14) to yield
y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1 (PT7.16)
Thus, we have “pinned down’’ Eq. (PT7.14) by forcing it to pass through the initial
condition, and in so doing, we have developed a unique solution to the ODE and have
come full circle to the original function [Eq. (PT7.12)].
Initial conditions usually have very tangible interpretations for differential equations
derived from physical problem settings. For example, in the falling parachutist problem,
the initial condition was re! ective of the physical fact that at time zero the vertical veloc-
ity was zero. If the parachutist had already been in vertical motion at time zero, the
solution would have been modi ed to account for this initial velocity.
When dealing with an nth-order differential equation, n conditions are required to
obtain a unique solution. If all conditions are speci ed at the same value of the indepen-
dent variable (for example, at x or t 5 0), then the problem is called an initial-value
problem. This is in contrast to boundary-value problems where speci cation of conditions
occurs at different values of the independent variable. Chapters 25 and 26 will focus on
initial-value problems. Boundary-value problems are covered in Chap. 27 along with
Before proceeding to numerical methods for solving ordinary differential equations, some
orientation might be helpful. The following material is intended to provide you with an
overview of the material discussed in Part Seven. In addition, we have formulated objec-
tives to focus your studies of the subject area.
PT7.3.1 Scope and Preview
Figure PT7.5 provides an overview of Part Seven. Two broad categories of numerical
methods for initial-value problems will be discussed in this part of this book. One-step
methods, which are covered in Chap. 25, permit the calculation of yi11, given the dif-
ferential equation and yi. Multistep methods, which are covered in Chap. 26, require
additional values of y other than at i.
With all but a minor exception, the one-step methods in Chap. 25 belong to what
are called Runge-Kutta techniques. Although the chapter might have been organized
around this theoretical notion, we have opted for a more graphical, intuitive approach to
introduce the methods. Thus, we begin the chapter with Euler’s method, which has a
very straightforward graphical interpretation. Then, we use visually oriented arguments
706 ORDINARY DIFFERENTIAL EQUATIONS
FIGURE PT7.5 Schematic representation of the organization of Part Seven: Ordinary Differential Equations.
26.2 Multistep methods
PT 7.2 Mathematical background
PT 7.6 Advanced methods
PT 7.5 Important formulas
28.4 Mechanical engineering
27.3 Software packages
PT 7.4 Trade-offs
PT 7.3 Orientation
PT 7.1 Motivation
25.2 Heun and midpoint methods
25.1 Euler’s method
25.4 Systems of
25.5 Adaptive RK
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.
FIGURE 25.3 Comparison of the true solution with a numerical solution using Euler’s method for the integral of y9 5 22×3 1 12×2 2 20x 1 8.5 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.
h = 0.5
25.1 EULER’S METHOD 713
2. Round-off errors caused by the limited numbers of signi! cant digits that can be
retained by a computer.
The truncation errors are composed of two parts. The ! rst is a local truncation error
that results from an application of the method in question over a single step. The second
is a propagated truncation error that results from the approximations produced during
the previous steps. The sum of the two is the total, or global truncation, error.
Insight into the magnitude and properties of the truncation error can be gained by
deriving Euler’s method directly from the Taylor series expansion. To do this, realize that
the differential equation being integrated will be of the general form
y¿ 5 f(x, y) (25.3)
where y9 5 dyydx and x and y are the independent and the dependent variables, respectively. If the solution—that is, the function describing the behavior of y—has continuous deriva-
tives, it can be represented by a Taylor series expansion about a starting value (xi, yi), as in
[recall Eq. (4.7)]
yi11 5 yi 1 y¿i h 1 y–i
2! h2 1 p 1
n! hn 1 Rn (25.4)
where h 5 xi11 2 xi and Rn 5 the remainder term, de! ned as
Rn 5 y(n11)(j)
(n 1 1)! hn11 (25.5)
where j lies somewhere in the interval from xi to xi11. An alternative form can be de-
veloped by substituting Eq. (25.3) into Eqs. (25.4) and (25.5) to yield
yi11 5 yi 1 f(xi, yi)h 1 f ¿(xi, yi)
2! h2 1 p 1
f (n21)(xi, yi)
n! hn 1 O(hn11) (25.6)
where O(hn11) speci! es that the local truncation error is proportional to the step size
raised to the (n 1 1)th power.
By comparing Eqs. (25.2) and (25.6), it can be seen that Euler’s method corresponds
to the Taylor series up to and including the term f(xi, yi)h. Additionally, the comparison
indicates that a truncation error occurs because we approximate the true solution using
a ! nite number of terms from the Taylor series. We thus truncate, or leave out, a part of
the true solution. For example, the truncation error in Euler’s method is attributable to
the remaining terms in the Taylor series expansion that were not included in Eq. (25.2).
Subtracting Eq. (25.2) from Eq. (25.6) yields
Et 5 f ¿(xi, yi)
2! h2 1 p 1 O(hn11) (25.7)
where Et 5 the true local truncation error. For suf! ciently small h, the errors in the terms
in Eq. (25.7) usually decrease as the order increases (recall Example 4.2 and the ac-
companying discussion), and the result is often represented as
Ea 5 f ¿(xi, yi)
2! h2 (25.8)
714 RUNGE-KUTTA METHODS
Ea 5 O(h 2) (25.9)
where Ea 5 the approximate local truncation error.
EXAMPLE 25.2 Taylor Series Estimate for the Error of Euler’s Method
Problem Statement. Use Eq. (25.7) to estimate the error of the ! rst step of Example 25.1. Also use it to determine the error due to each higher-order term of the Taylor series
Solution. Because we are dealing with a polynomial, we can use the Taylor series to obtain exact estimates of the errors in Euler’s method. Equation (25.7) can be written as
Et 5 f ¿(xi, yi)
2! h2 1
f –(xi, yi)
3! h3 1
f (3)(xi, yi)
4! h4 (E25.2.1)
where f 9(xi, yi) 5 the ! rst derivative of the differential equation (that is, the second de-
rivative of the solution). For the present case, this is
f ¿(xi, yi) 5 26x 2
1 24x 2 20 (E25.2.2)
and f 0(xi, yi) is the second derivative of the ODE
f –(xi, yi) 5 212x 1 24 (E25.2.3)
and f (3)(xi, yi) is the third derivative of the ODE
f (3)(xi, yi) 5 212 (E25.2.4)
We can omit additional terms (that is, fourth derivatives and higher) from Eq. (E25.2.1)
because for this particular case they equal zero. It should be noted that for other func-
tions (for example, transcendental functions such as sinusoids or exponentials) this would
not necessarily be true, and higher-order terms would have nonzero values. However, for
the present case, Eqs. (E25.2.1) through (E25.2.4) completely de! ne the truncation error
for a single application of Euler’s method.
For example, the error due to truncation of the second-order term can be calculated as
Et, 2 5 26(0.0)2 1 24(0.0) 2 20
2 (0.5)2 5 22.5 (E25.2.5)
For the third-order term:
Et, 3 5 212(0.0) 1 24
6 (0.5)3 5 0.5
and the fourth-order term:
Et, 4 5 212
24 (0.5)4 5 20.03125
These three results can be added to yield the total truncation error:
Et 5 Et, 2 1 Et, 3 1 Et, 4 5 22.5 1 0.5 2 0.03125 5 22.03125
25.1 EULER’S METHOD 715
which is exactly the error that was incurred in the initial step of Example 25.1. Note
how Et, 2 . Et, 3 . Et, 4, which supports the approximation represented by Eq. (25.8).
As illustrated in Example 25.2, the Taylor series provides a means of quantifying
the error in Euler’s method. However, there are limitations associated with its use for
1. The Taylor series provides only an estimate of the local truncation error—that is, the
error created during a single step of the method. It does not provide a measure of the
propagated and, hence, the global truncation error. In Table 25.1, we have included
the local and global truncation errors for Example 25.1. The local error was computed
for each time step with Eq. (25.2) but using the true value of yi (the second column
of the table) to compute each yi1l rather than the approximate value (the third column),
as is done in the Euler method. As expected, the average absolute local truncation
error (25 percent) is less than the average global error (90 percent). The only reason
that we can make these exact error calculations is that we know the true value a
priori. Such would not be the case in an actual problem. Consequently, as discussed
below, you must usually apply techniques such as Euler’s method using a number of
different step sizes to obtain an indirect estimate of the errors involved.
2. As mentioned above, in actual problems we usually deal with functions that are more
complicated than simple polynomials. Consequently, the derivatives that are needed
to evaluate the Taylor series expansion would not always be easy to obtain.
Although these limitations preclude exact error analysis for most practical problems,
the Taylor series still provides valuable insight into the behavior of Euler’s method. Ac-
cording to Eq. (25.9), we see that the local error is proportional to the square of the step
size and the ! rst derivative of the differential equation. It can also be demonstrated that
the global truncation error is O(h), that is, it is proportional to the step size (Carnahan
et al., 1969). These observations lead to some useful conclusions:
1. The error can be reduced by decreasing the step size.
2. The method will provide error-free predictions if the underlying function (that is, the
solution of the differential equation) is linear, because for a straight line the second
derivative would be zero.
This latter conclusion makes intuitive sense because Euler’s method uses straight-line
segments to approximate the solution. Hence, Euler’s method is referred to as a rst-
It should also be noted that this general pattern holds for the higher-order one-step
methods described in the following pages. That is, an nth-order method will yield perfect
results if the underlying solution is an nth-order polynomial. Further, the local truncation
error will be O(hn11) and the global error O(hn).
EXAMPLE 25.3 Effect of Reduced Step Size on Euler’s Method
Problem Statement. Repeat the computation of Example 25.1 but use a step size of 0.25.
716 RUNGE-KUTTA METHODS
Solution. The computation is repeated, and the results are compiled in Fig. 25.4a. Halving the step size reduces the absolute value of the average global error to 40 percent
and the absolute value of the local error to 6.4 percent. This is compared to global and
local errors for Example 25.1 of 90 percent and 24.8 percent, respectively. Thus, as
expected, the local error is quartered and the global error is halved.
Also, notice how the local error changes sign for intermediate values along the range.
This is due primarily to the fact that the ! rst derivative of the differential equation is a
parabola that changes sign [recall Eq. (E25.2.2) and see Fig. 25.4b]. Because the local
error is proportional to this function, the net effect of the oscillation in sign is to keep the
global error from continuously growing as the calculation proceeds. Thus, from x 5 0 to
x 5 1.25, the local errors are all negative, and consequently, the global error increases
FIGURE 25.4 (a) Comparison of two numerical solutions with Euler’s method using step sizes of 0.5 and 0.25. (b) Comparison of true and estimated local truncation error for the case where the step size is 0.5. Note that the “estimated” error is based on Eq. (E25.2.5).
h = 0.5
h = 0.25
25.1 EULER’S METHOD 717
over this interval. In the intermediate section of the range, positive local errors begin to
reduce the global error. Near the end, the process is reversed and the global error again
in” ates. If the local error continuously changes sign over the computation interval, the net
effect is usually to reduce the global error. However, where the local errors are of the
same sign, the numerical solution may diverge farther and farther from the true solution
as the computation proceeds. Such results are said to be unstable.
The effect of further step-size reductions on the global truncation error of Euler’s
method is illustrated in Fig. 25.5. This plot shows the absolute percent relative error at
x 5 5 as a function of step size for the problem we have been examining in Examples
25.1 through 25.3. Notice that even when h is reduced to 0.001, the error still exceeds
0.1 percent. Because this step size corresponds to 5000 steps to proceed from x 5 0 to
x 5 5, the plot suggests that a ! rst-order technique such as Euler’s method demands
great computational effort to obtain acceptable error levels. Later in this chapter, we
present higher-order techniques that attain much better accuracy for the same computa-
tional effort. However, it should be noted that, despite its inef! ciency, the simplicity of
Euler’s method makes it an extremely attractive option for many engineering problems.
Because it is very easy to program, the technique is particularly useful for quick analy-
ses. In the next section, a computer algorithm for Euler’s method is developed.
FIGURE 25.5 Effect of step size on the global truncation error of Euler’s method for the integral of y9 5 22×3 1 12×2 2 20x 1 8.5. The plot shows the absolute percent relative global error at x 5 5 as a function of step size.
p e rc
v e e
718 RUNGE-KUTTA METHODS
25.1.2 Algorithm for Euler’s Method
Algorithms for one-step techniques such as Euler’s method are extremely simple to
program. As speci! ed previously at the beginning of this chapter, all one-step methods
have the general form
New value 5 old value 1 slope 3 step size (25.10)
The only way in which the methods differ is in the calculation of the slope.
Suppose that you want to perform the simple calculation outlined in Table 25.1. That
is, you would like to use Euler’s method to integrate y9 5 22×3 1 12×2 2 20x 1 8.5,
with the initial condition that y 5 1 at x 5 0. You would like to integrate out to x 5 4
using a step size of 0.5, and display all the results. A simple pseudocode to accomplish
this task could be written as in Fig. 25.6.
Although this program will “do the job” of duplicating the results of Table 25.1, it
is not very well designed. First, and foremost, it is not very modular. Although this is
not very important for such a small program, it would be critical if we desired to mod-
ify and improve the algorithm.
Further, there are a number of issues related to the way we have set up the iterations.
For example, suppose that the step size were to be made very small to obtain better ac-
curacy. In such cases, because every computed value is displayed, the number of output
values might be very large. Further, the algorithm is predicated on the assumption that
the calculation interval is evenly divisible by the step size. Finally, the accumulation of
x in the line x 5 x 1 dx can be subject to quantizing errors of the sort previously dis-
FIGURE 25.6 Pseudocode for a “dumb” version of Euler’s method.
‘set integration range
xi 5 0
xf 5 4
x 5 xi
y 5 1
‘set step size and determine
‘number of calculation steps
dx 5 0.5
nc 5 (xf 2 xi)/dx
‘output initial condition
PRINT x, y
‘loop to implement Euler’s method
‘and display results
DOFOR i 5 1, nc
dydx 5 22×3 1 12×2 2 20x 1 8.5
y 5 y 1 dydx ? dx
x 5 x 1 dx
PRINT x, y
25.1 EULER’S METHOD 719
cussed in Sec. 3.4.1. For example, if dx were changed to 0.01 and standard IEEE ” oat-
ing point representation were used (about seven signi! cant digits), the result at the end
of the calculation would be 3.999997 rather than 4. For dx 5 0.001, it would be 3.999892!
A much more modular algorithm that avoids these dif! culties is displayed in Fig. 25.7.
The algorithm does not output all calculated values. Rather, the user speci! es an output
interval, xout, that dictates the interval at which calculated results are stored in arrays, xpm
and ypm. These values are stored in arrays so that they can be output in a variety of ways
after the computation is completed (for example, printed, graphed, or written to a ! le).
The Driver Program takes big output steps and calls an Integrator routine that takes
! ner calculation steps. Note that the loops controlling both large and small steps exit on
logical conditions. Thus, the intervals do not have to be evenly divisible by the step sizes.
The Integrator routine then calls an Euler routine that takes a single step with Euler’s
method. The Euler routine calls a Derivative routine that calculates the derivative value.
Whereas such modularization might seem like overkill for the present case, it will
greatly facilitate modifying the program in later sections. For example, although the
program in Fig. 25.7 is speci! cally designed to implement Euler’s method, the Euler
module is the only part that is method-speci! c. Thus, all that is required to apply this
algorithm to the other one-step methods is to modify this routine.
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
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
(b) Routine to Take One Output Step
SUB Integrator (x, y, h, xend)
IF (xend 2 x , h) THEN h 5 xend 2 x
CALL Euler (x, y, h, ynew)
y 5 ynew
IF (x $ xend) EXIT
(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
(d) Routine to Determine Derivative
SUB Derivs (x, y, dydx)
dydx 5 . . .
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
dt 5 g 2
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
dt 5 g 2
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
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.
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)
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)
The second derivative is
f –(xi, yi) 5 0[0fy0x 1 (0fy0y)(dyydx)]
0x 1 0[0fy0x 1 (0fy0y)(dyydx)]
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.
xxi + 1xi
Slope = f (xi, yi)
Slope = f (xi + 1, y 0 i + 1)
xxi + 1xi
Slope = f (xi, yi) + f (xi + 1, yi + 1)
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
f(xi, yi) 1 f(xi11, y 0 i11)
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
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
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 +
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
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-
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
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)
726 RUNGE-KUTTA METHODS
Notice the similarity between the right-hand side of Eq. (25.18) and the trapezoidal rule
[Eq. (21.3)]. The connection between the two methods can be formally demonstrated by
starting with the ordinary differential equation
dx 5 f(x)
This equation can be solved for y by integration:
dy 5 #
f(x) dx (25.19)
yi11 2 yi 5 # xi11
f(x) dx (25.20)
yi11 5 yi 1 # xi11
f(x) dx (25.21)
Now, recall from Chap. 21 that the trapezoidal rule [Eq. (21.3)] is de! ned as
f(x) dx >
f(xi) 1 f(xi11)
2 h (25.22)
where h 5 xi11 2 xi. Substituting Eq. (25.22) into Eq. (25.21) yields
yi11 5 yi 1 f(xi) 1 f(xi11)
2 h (25.23)
which is equivalent to Eq. (25.18).
Because Eq. (25.23) is a direct expression of the trapezoidal rule, the local truncation
error is given by [recall Eq. (21.6)]
Et 5 2 f –(j)
12 h3 (25.24)
where j is between xi and xi1l. Thus, the method is second order because the second de-
rivative of the ODE is zero when the true solution is a quadratic. In addition, the local and
global errors are O(h3) and O(h2), respectively. Therefore, decreasing the step size decreases
the error at a faster rate than for Euler’s method. Figure 25.11, which shows the result of
using Heun’s method to solve the polynomial from Example 25.1 demonstrates this behavior.
25.2.2 The Midpoint (or Improved Polygon) Method
Figure 25.12 illustrates another simple modi! cation of Euler’s method. Called the mid-
point method (or the improved polygon or the modi ed Euler), this technique uses Euler’s
method to predict a value of y at the midpoint of the interval (Fig. 25.12a):
yi11y2 5 yi 1 f(xi, yi) h
25.2 IMPROVEMENTS OF EULER’S METHOD 727
FIGURE 25.11 Comparison of the true solution with a numerical solution using Euler’s and Heun’s methods for the integral of y9 5 22×3 1 12×2 2 20x 1 8.5.
xxi + 1xi
y Slope = f (xi + 1/2, yi + 1/2)
xxi + 1/2xi
Slope = f (xi + 1/2, yi + 1/2)
FIGURE 25.12 Graphical depiction of the midpoint method. (a) Eq. (25.25) and (b) Eq. (25.27).
728 RUNGE-KUTTA METHODS
Then, this predicted value is used to calculate a slope at the midpoint:
y¿i11y2 5 f(xi11y2, yi11y2) (25.26)
which is assumed to represent a valid approximation of the average slope for the entire
interval. This slope is then used to extrapolate linearly from xi to xi1l (Fig. 25.12b):
yi11 5 yi 1 f(xi11y2, yi11y2)h (25.27)
Observe that because yi1l is not on both sides, the corrector [Eq. (25.27)] cannot be ap-
plied iteratively to improve the solution.
As in the previous section, this approach can also be linked to Newton-Cotes inte-
gration formulas. Recall from Table 21.4, that the simplest Newton-Cotes open integra-
tion formula, which is called the midpoint method, can be represented as
a f(x) dx > (b 2 a) f(x1)
where x1 is the midpoint of the interval (a, b). Using the nomenclature for the present
case, it can be expressed as
f(x) dx > h f(xi11y2)
Substitution of this formula into Eq. (25.21) yields Eq. (25.27). Thus, just as the Heun
method can be called the trapezoidal rule, the midpoint method gets its name from the
underlying integration formula upon which it is based.
The midpoint method is superior to Euler’s method because it utilizes a slope estimate
at the midpoint of the prediction interval. Recall from our discussion of numerical differ-
entiation in Sec. 4.1.3 that centered ! nite divided differences are better approximations of
derivatives than either forward or backward versions. In the same sense, a centered ap-
proximation such as Eq. (25.26) has a local truncation error of O(h2) in comparison with
the forward approximation of Euler’s method, which has an error of O(h). Consequently,
the local and global errors of the midpoint method are O(h3) and O(h2), respectively.
25.2.3 Computer Algorithms for Heun and Midpoint Methods
Both the Heun method with a single corrector and the midpoint method can be easily
programmed using the general structure depicted in Fig. 25.7. As in Fig. 25.13a and b,
simple routines can be written to take the place of the Euler routine in Fig. 25.7.
However, when the iterative version of the Heun method is to be implemented, the
modi! cations are a bit more involved (although they are still localized within a single
module). We have developed pseudocode for this purpose in Fig. 25.13c. This algorithm
can be combined with Fig. 25.7 to develop software for the iterative Heun method.
By tinkering with Euler’s method, we have derived two new second-order techniques. Even
though these versions require more computational effort to determine the slope, the accom-
panying reduction in error will allow us to conclude in a subsequent section (Sec. 25.3.4)
25.3 RUNGE-KUTTA METHODS 729
that the improved accuracy is usually worth the effort. Although there are certain cases where
easily programmable techniques such as Euler’s method can be applied to advantage, the
Heun and midpoint methods are generally superior and should be implemented if they are
consistent with the problem objectives.
As noted at the beginning of this section, the Heun (without iterations), the midpoint
method, and in fact, the Euler technique itself are versions of a broader class of one-step
approaches called Runge-Kutta methods. We now turn to a formal derivation of these
25.3 RUNGE-KUTTA METHODS
Runge-Kutta (RK) methods achieve the accuracy of a Taylor series approach without
requiring the calculation of higher derivatives. Many variations exist but all can be cast
in the generalized form of Eq. (25.1):
yi11 5 yi 1 f(xi, yi, h)h (25.28)
where f(xi, yi, h) is called an increment function, which can be interpreted as a represen-
tative slope over the interval. The increment function can be written in general form as
f 5 a1k1 1 a2k2 1 p 1 ankn (25.29)
(a) Simple Heun without Iteration
SUB Heun (x, y, h, ynew)
CALL Derivs (x, y, dy1dx)
ye 5 y 1 dy1dx ? h
CALL Derivs(x 1 h, ye, dy2dx)
Slope 5 (dy1dx 1 dy2dx)y2 ynew 5 y 1 Slope ? h
x 5 x 1 h
(b) Midpoint Method
SUB Midpoint (x, y, h, ynew)
CALL Derivs(x, y, dydx)
ym 5 y 1 dydx ? hy2 CALL Derivs (x 1 hy2, ym, dymdx) ynew 5 y 1 dymdx ? h
x 5 x 1 h
(c) Heun with Iteration
SUB HeunIter (x, y, h, ynew)
es 5 0.01
maxit 5 20
CALL Derivs(x, y, dy1dx)
ye 5 y 1 dy1dx ? h
iter 5 0
yeold 5 ye
CALL Derivs(x 1 h, ye, dy2dx)
slope 5 (dy1dx 1 dy2dx)y2 ye 5 y 1 slope ? h
iter 5 iter 1 1
ea 5 ` ye 2 yeold ye
` 100% IF (ea # es OR iter . maxit) EXIT
ynew 5 ye
x 5 x 1 h
FIGURE 25.13 Pseudocode to implement the (a) simple Heun, (b) midpoint, and (c) iterative Heun methods.
730 RUNGE-KUTTA METHODS
where the a’s are constants and the k’s are
k1 5 f(xi, yi) (25.29a)
k2 5 f(xi 1 p1h, yi 1 q11k1h) (25.29b)
k3 5 f(xi 1 p2h, yi 1 q21k1h 1 q22k2h) (25.29c)
kn 5 f(xi 1 pn21h, yi 1 qn21, 1k1h 1 qn21, 2k2h 1 p 1 qn21, n21kn21h) (25.29d)
where the p’s and q’s are constants. Notice that the k’s are recurrence relationships. That
is, k1 appears in the equation for k2, which appears in the equation for k3, and so forth.
Because each k is a functional evaluation, this recurrence makes RK methods ef! cient
for computer calculations.
Various types of Runge-Kutta methods can be devised by employing different num-
bers of terms in the increment function as speci! ed by n. Note that the ! rst-order RK
method with n 5 1 is, in fact, Euler’s method. Once n is chosen, values for the a’s, p’s,
and q’s are evaluated by setting Eq. (25.28) equal to terms in a Taylor series expansion
(Box 25.1). Thus, at least for the lower-order versions, the number of terms, n, usually
represents the order of the approach. For example, in the next section, second-order RK
methods use an increment function with two terms (n 5 2). These second-order methods
will be exact if the solution to the differential equation is quadratic. In addition, because
terms with h3 and higher are dropped during the derivation, the local truncation error is
O(h3) and the global error is O(h2). In subsequent sections, the third- and fourth-order
RK methods (n 5 3 and 4, respectively) are developed. For these cases, the global trun-
cation errors are O(h3) and O(h4), respectively.
25.3.1 Second-Order Runge-Kutta Methods
The second-order version of Eq. (25.28) is
yi11 5 yi 1 (a1k1 1 a2k2)h (25.30)
k1 5 f(xi, yi) (25.30a)
k2 5 f(xi 1 p1h, yi 1 q11k1h) (25.30b)
As described in Box 25.1, values for al, a2, p1, and q11 are evaluated by setting Eq. (25.30)
equal to a Taylor series expansion to the second-order term. By doing this, we derive three
equations to evaluate the four unknown constants. The three equations are
a1 1 a2 5 1 (25.31)
a2 p1 5 1
a2q11 5 1
25.3 RUNGE-KUTTA METHODS 731
Because we have three equations with four unknowns, we must assume a value of one
of the unknowns to determine the other three. Suppose that we specify a value for a2. Then
Eqs. (25.31) through (25.33) can be solved simultaneously for
a1 5 1 2 a2 (25.34)
p1 5 q11 5 1
Box 25.1 Derivation of the Second-Order Runge-Kutta Methods
The second-order version of Eq. (25.28) is
yi11 5 yi 1 (a1k1 1 a2k2)h (B25.1.1)
k1 5 f (xi, yi) (B25.1.2)
k2 5 f (xi 1 p1h, yi 1 q11k1h) (B25.1.3)
To use Eq. (B25.1.1) we have to determine values for the
constants a1, a2, p1, and q11. To do this, we recall that the second-
order Taylor series for yi11 in terms of yi and f(xi, yi) is written as
yi11 5 yi 1 f (xi, yi)h 1 f ¿(xi, yi)
2! h2 (B25.1.4)
where f 9(xi, yi) must be determined by chain-rule differentiation
f ¿(xi, yi) 5 0f (x, y)
0x 1 0f (x, y)
Substituting Eq. (B25.1.5) into (B25.1.4) gives
yi11 5 yi 1 f (xi, yi)h 1 a 0f 0x
dx b h2
The basic strategy underlying Runge-Kutta methods is to use alge-
braic manipulations to solve for values of a1, a2, p1, and q11 that
make Eqs. (B25.1.1) and (B25.1.6) equivalent.
To do this, we ! rst use a Taylor series to expand Eq. (25.1.3).
The Taylor series for a two-variable function is de! ned as [recall
g(x 1 r, y 1 s) 5 g(x, y) 1 r 0g
0x 1 s 0g
Applying this method to expand Eq. (B25.1.3) gives
f (xi 1 p1h, yi 1 q11k1h) 5 f (xi, yi) 1 p1h 0f
1 q11k1h 0f
0y 1 O(h2)
This result can be substituted along with Eq. (B25.1.2) into Eq.
(B25.1.1) to yield
yi11 5 yi 1 a1h f (xi, yi) 1 a2h f(xi, yi) 1 a2 p1h 2
1 a2q11h 2 f (xi, yi)
0y 1 O(h3)
or, by collecting terms,
yi11 5 yi 1 [a1 f (xi, yi) 1 a2 f (xi, yi) ]h
1 c a2 p1 0f 0x
1 a2q11 f (xi, yi) 0f
0y d h2 1 O(h3)
Now, comparing like terms in Eqs. (B25.1.6) and (B25.1.7), we
determine that for the two equations to be equivalent, the following
a1 1 a2 5 1
a2 p1 5 1
a2q11 5 1
These three simultaneous equations contain the four unknown con-
stants. Because there is one more unknown than the number of
equations, there is no unique set of constants that satisfy the equa-
tions. However, by assuming a value for one of the constants, we
can determine the other three. Consequently, there is a family of
second-order methods rather than a single version.
732 RUNGE-KUTTA METHODS
Because we can choose an in! nite number of values for a2, there are an in! nite
number of second-order RK methods. Every version would yield exactly the same results
if the solution to the ODE were quadratic, linear, or a constant. However, they yield
different results when (as is typically the case) the solution is more complicated. We
present three of the most commonly used and preferred versions:
Heun Method with a Single Corrector (a2 5 1y2). If a2 is assumed to be 1y2, Eqs. (25.34) and (25.35) can be solved for a1 5 1y2 and pl 5 q11 5 1. These parameters, when substituted into Eq. (25.30), yield
yi11 5 yi 1 a1 2
k1 1 1
2 k2b h (25.36)
k1 5 f(xi, yi) (25.36a)
k2 5 f(xi 1 h, yi 1 k1h) (25.36b)
Note that k1 is the slope at the beginning of the interval and k2 is the slope at the end
of the interval. Consequently, this second-order Runge-Kutta method is actually Heun’s
technique without iteration.
The Midpoint Method (a2 5 1). If a2 is assumed to be 1, then a1 5 0, p1 5 q11 5 1y2, and Eq. (25.30) becomes
yi11 5 yi 1 k 2h (25.37)
k1 5 f(xi, yi) (25.37a)
k2 5 f axi 1 1 2
h, yi 1 1
2 k1hb (25.37b)
This is the midpoint method.
Ralston’s Method (a2 5 2y3). Ralston (1962) and Ralston and Rabinowitz (1978) determined that choosing a2 5 2y3 provides a minimum bound on the truncation error for the second-order RK algorithms. For this version, a1 5 1y3 and p1 5 q11 5 3y4 and yields
yi11 5 yi 1 a1 3
k1 1 2
3 k2b h (25.38)
k1 5 f(xi, yi) (25.38a)
k2 5 f axi 1 3 4
h, yi 1 3
4 k1hb (25.38b)
25.3 RUNGE-KUTTA METHODS 733
EXAMPLE 25.6 Comparison of Various Second-Order RK Schemes
Problem Statement. Use the midpoint method [Eq. (25.37)] and Ralston’s method [Eq. (25.38)] to numerically integrate Eq. (PT7.13)
f(x, y) 5 22×3 1 12×2 2 20x 1 8.5
from x 5 0 to x 5 4 using a step size of 0.5. The initial condition at x 5 0 is y 5 1.
Compare the results with the values obtained using another second-order RK algorithm,
that is, the Heun method without corrector iteration (Table 25.3).
Solution. The ! rst step in the midpoint method is to use Eq. (25.37a) to compute
k1 5 22(0) 3
1 12(0)2 2 20(0) 1 8.5 5 8.5
However, because the ODE is a function of x only, this result has no bearing on the
second step—the use of Eq. (25.37b) to compute
k2 5 22(0.25) 3
1 12(0.25)2 2 20(0.25) 1 8.5 5 4.21875
Notice that this estimate of the slope is much closer to the average value for the interval
(4.4375) than the slope at the beginning of the interval (8.5) that would have been used
for Euler’s approach. The slope at the midpoint can then be substituted into Eq. (25.37)
y(0.5) 5 1 1 4.21875(0.5) 5 3.109375 et 5 3.4%
The computation is repeated, and the results are summarized in Fig. 25.14 and Table 25.3.
FIGURE 25.14 Comparison of the true solution with numerical solutions using three second-order RK methods and Euler’s method.
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
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)
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)
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.
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
(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.
(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
yi11 5 yi 1 1
90 (7k1 1 32k3 1 12k4 1 32k5 1 7k6)h (25.41)
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
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
16 k4hb (25.41e)
k6 5 f axi 1 h, yi 2 3 7
k1h 1 2
7 k2h 1
7 k3h 2
7 k4h 1
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
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.
P e rc
v e e
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
dx 5 f1(x, y1, y2, p , yn)
dx 5 f2(x, y1, y2, p , yn)
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
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
dx 5 20.5y1
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)
2 5 3.5
y2 1 k1, 2 h
2 5 6 1 (1.8)
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)
2 5 3.5625
y2 1 k2, 2 h
2 5 6 1 (1.715)
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
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)]
dx 5 y2
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)]
dx 5 y4
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
xi 5 initial value independent
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
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
IF (x $ xf) EXIT
(b) Routine to Take One Output Step
SUB Integrator (x, y, n, h, xend)
IF (xend 2 x , h) THEN h 5 xend 2 x
CALL RK4 (x, y, n, h)
IF (x $ xend) EXIT
(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
CALL Derivs (x 1 h / 2, ym, k2)
DOFOR i 5 1, n
ymi 5 yi 1 k2i * h / 2
CALL Derivs (x 1 h / 2, ym, k3)
DOFOR i 5 1, n
yei 5 yi 1 k3i * h
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
x 5 x 1 h
(d ) Routine to Determine Derivatives
SUB Derivs (x, y, dy)
dy1 5 …
dy2 5 …
FIGURE 25.18 Pseudocode for the fourth-order RK method for systems.
744 RUNGE-KUTTA METHODS
(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.
0 2 31
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
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.
0 1 2 3
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 ¢
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
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
594 k4 1
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
55,296 k4 1
14,336 k5 1
4 k6b h (25.46)
k1 5 f(xi, yi)
k2 5 f axi 1 1 5
h, yi 1 1
k3 5 f axi 1 3 10
h, yi 1 3
40 k1h 1
k4 5 f axi 1 3 5
h, yi 1 3
10 k1h 2
10 k2h 1
k5 5 f axi 1 h, yi 2 11 54
k1h 1 5
2 k2h 2
27 k3h 1
k6 5 f axi 1 7 8
h, yi 1 1631
55,296 k1h 1
512 k2h 1
13,824 k3h 1
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
594 6.832587 1
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
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)
CALL Derivs (x1a2*h,ytemp,k2)
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)
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
hi5.5; tiny 5 1. 3 10230
print *, xi,yi
IF (istep . maxstep AND x # xf) EXIT
IF (x1h.xf) THEN h5xf2x
CALL Adapt (x,y,dy,h,yscal,eps,hnxt)
(b) Adaptive Step Routine
SUB Adapt (x,y,dy,htry,yscal,eps,hnxt)
PARAMETER (safety50.9, econ51.89e24)
CALL RKkc (y,dy,x,h,ytemp,yerr)
IF emax # 1 EXIT
IF xnew5x THEN pause
IF emax . econ THEN
hnxt5safety*emax 2.2 *h
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
0 2 4 x
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
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
(b) Euler’s method with h 5 0.5 and 0.25.
(c) Midpoint method with h 5 0.5.
(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
(b) Euler’s method.
(c) Heun’s method without iteration.
(d) Ralston’s method.
(e) Fourth-order RK method.
25.3 Use the (a) Euler and (b) Heun (without iteration) methods to
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:
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
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
25.14 Develop a user-friendly computer program for the classical
fourth-order RK method. Test the program by duplicating Exam-
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.
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
dt 2 1 (5x)
dt 1 (x 1 7) sin(vt) 5 0
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
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.
FIGURE P25.20 A spherical tank.
754 RUNGE-KUTTA METHODS
from its equilibrium position (m), and g 5 gravitational acceleration
(9.81 m/s2). Solve these equations for the positions and velocities of
the three jumpers given the initial conditions that all positions and
velocities are zero at t 5 0. Use the following parameters for your
calculations: m1 5 60 kg, m2 5 70 kg, m3 5 80 kg, k1 5 k3 5 50,
and k2 5 100 (N/m).
25.24 Given the initial conditions, y(0) 5 1 and y9(0) 5 0, solve
the following initial-value problem from t 5 0 to 4:
dt 2 1 4y 5 0
Obtain your solutions with (a) Euler’s method and (b) the fourth-
order RK method. In both cases, use a step size of 0.125. Plot both
solutions on the same graph along with the exact solution y 5 cos 2t.
25.25 Use the following differential equations to compute the
velocity and position of a soccer ball that is kicked straight up in the
air with an initial velocity of 40 m/s:
dt 5 y
dt 5 2g 2
m y Zy Z
where y 5 upward distance (m), t 5 time (s), y 5 upward velocity
(m/s), g 5 gravitational constant (5 9.81 m/s2), cd 5 drag coef” –
cient (kg/m), and m 5 mass (kg). Note that the drag coef” cient is
related to more fundamental parameters by
cd 5 1
where r 5 air density (kg/m3), A 5 area (m2), and Cd 5 the di-
mensionless drag coef” cient. Use the following parameter values
for your calculation: d 5 22 cm, m 5 0.4 kg, r 5 1.3 kg/m3, and
Cd 5 0.52.
25.26 Three linked bungee jumpers are depicted in Fig. P25.26. If
the bungee cords are idealized as linear springs (i.e., governed by
Hooke’s law), the following differential equations based on force
balances can be developed
m1 d 2×1
dt 2 5 m1g 1 k2(x2 2 x1) 2 k1x1
m2 d 2×2
dt 2 5 m2g 1 k3(x3 2 x2) 1 k2(x1 2 x2)
m3 d 2×3
dt 2 5 m3g 1 k3(x2 2 x3)
where mi 5 the mass of jumper i (kg), kj 5 the spring constant for
cord j (N/m), xi 5 the displacement of jumper i measured downward
FIGURE P25.26 Three individuals connected by bungee cords.
x1 = 0
(a) Unstretched (b) Stretched
x2 = 0
x3 = 0
26 C H A P T E R 26
Stiffness and Multistep Methods
This chapter covers two areas. First, we describe stiff ODEs. These are both indi-
vidual and systems of ODEs that have both fast and slow components to their solution.
We introduce the idea of an implicit solution technique as one commonly used remedy
for this problem. Then we discuss multistep methods. These algorithms retain informa-
tion of previous steps to more effectively capture the trajectory of the solution. They
also yield the truncation error estimates that can be used to implement adaptive step-
Stiffness is a special problem that can arise in the solution of ordinary differential equa-
tions. A stiff system is one involving rapidly changing components together with slowly
changing ones. In many cases, the rapidly varying components are ephemeral transients
that die away quickly, after which the solution becomes dominated by the slowly varying
components. Although the transient phenomena exist for only a short part of the integra-
tion interval, they can dictate the time step for the entire solution.
Both individual and systems of ODEs can be stiff. An example of a single stiff
dt 5 21000y 1 3000 2 2000e2t (26.1)
If y(0) 5 0, the analytical solution can be developed as
y 5 3 2 0.998e21000t 2 2.002e2t (26.2)
As in Fig. 26.1, the solution is initially dominated by the fast exponential term
(e21000t). After a short period (t , 0.005), this transient dies out and the solution becomes
dictated by the slow exponential (e2t).
Insight into the step size required for stability of such a solution can be gained by
examining the homogeneous part of Eq. (26.1),
dt 5 2ay (26.3)
756 STIFFNESS AND MULTISTEP METHODS
If y(0) 5 y0, calculus can be used to determine the solution as
y 5 y0e 2at
Thus, the solution starts at y0 and asymptotically approaches zero.
Euler’s method can be used to solve the same problem numerically:
yi11 5 yi 1 dyi
Substituting Eq. (26.3) gives
yi11 5 yi 2 ayih
yi11 5 yi(1 2 ah) (26.4)
The stability of this formula clearly depends on the step size h. That is, 01 2 ah 0 must be less than 1. Thus, if h . 2ya, 0 yi 0 n q as i n q. For the fast transient part of Eq. (26.2), this criterion can be used to show that the step
size to maintain stability must be , 2y1000 5 0.002. In addition, it should be noted that, whereas this criterion maintains stability (that is, a bounded solution), an even smaller step
size would be required to obtain an accurate solution. Thus, although the transient occurs for
only a small fraction of the integration interval, it controls the maximum allowable step size.
Super! cially, you might suppose that the adaptive step-size routines described at the
end of the last chapter might offer a solution for this dilemma. You might think that they
would use small steps during the rapid transients and large steps otherwise. However,
this is not the case, because the stability requirement will still necessitate using very
small steps throughout the entire solution.
FIGURE 26.1 Plot of a stiff solution of a single ODE. Although the solution appears to start at 1, there is actually a fast transient from y 5 0 to 1 that occurs in less than 0.005 time unit. This transient is perceptible only when the response is viewed on the fi ner timescale in the inset.
0 42 t0
26.1 STIFFNESS 757
Rather than using explicit approaches, implicit methods offer an alternative remedy.
Such representations are called implicit because the unknown appears on both sides of
the equation. An implicit form of Euler’s method can be developed by evaluating the
derivative at the future time,
yi11 5 yi 1 dyi11
This is called the backward, or implicit, Euler’s method. Substituting Eq. (26.3) yields
yi11 5 yi 2 ayi11 h
which can be solved for
yi11 5 yi
1 1 ah (26.5)
For this case, regardless of the size of the step, 0 yi 0 n 0 as i n q. Hence, the approach is called unconditionally stable.
EXAMPLE 26.1 Explicit and Implicit Euler
Problem Statement. Use both the explicit and implicit Euler methods to solve
dt 5 21000y 1 3000 2 2000e2t
where y(0) 5 0. (a) Use the explicit Euler with step sizes of 0.0005 and 0.0015 to solve
for y between t 5 0 and 0.006. (b) Use the implicit Euler with a step size of 0.05 to
solve for y between 0 and 0.4.
(a) For this problem, the explicit Euler’s method is
yi11 5 yi 1 (21000yi 1 3000 2 2000e 2ti)h
The result for h 5 0.0005 is displayed in Fig. 26.2a along with the analytical solu-
tion. Although it exhibits some truncation error, the result captures the general shape
of the analytical solution. In contrast, when the step size is increased to a value just
below the stability limit (h 5 0.0015), the solution manifests oscillations. Using
h . 0.002 would result in a totally unstable solution, that is, it would go in! nite
as the solution progressed.
(b) The implicit Euler’s method is
yi11 5 yi 1 (21000yi11 1 3000 2 2000e 2ti11)h
Now because the ODE is linear, we can rearrange this equation so that yi11 is isolated
on the left-hand side,
yi11 5 yi 1 3000h 2 2000he
1 1 1000h
The result for h 5 0.05 is displayed in Fig. 26.2b along with the analytical solution.
Notice that even though we have used a much bigger step size than the one that
758 STIFFNESS AND MULTISTEP METHODS
induced instability for the explicit Euler, the numerical solution tracks nicely on
the analytical result.
FIGURE 26.2 Solution of a “stiff” ODE with (a) the explicit and (b) implicit Euler methods.
h = 0.0015
h = 0.0005
h = 0.05
Systems of ODEs can also be stiff. An example is
dt 5 25y1 1 3y2 (26.6a)
dt 5 100y1 2 301y2 (26.6b)
For the initial conditions y1(0) 5 52.29 and y2(0) 5 83.82, the exact solution is
y1 5 52.96e 23.9899t
2 0.67e2302.0101t (26.7a)
y2 5 17.83e 23.9899t
1 65.99e2302.0101t (26.7b)
Note that the exponents are negative and differ by about 2 orders of magnitude. As with
the single equation, it is the large exponents that respond rapidly and are at the heart of
the system’s stiffness.
26.2 MULTISTEP METHODS 759
An implicit Euler’s method for systems can be formulated for the present example as
y1, i11 5 y1, i 1 (25y1, i11 1 3y2, i11)h (26.8a)
y2, i11 5 y2, i 1 (100y1, i11 2 301y2, i11)h (26.8b)
Collecting terms gives
(1 1 5h)y1, i11 2 3hy2, i11 5 y1, i (26.9a)
2100hy1, i11 1 (1 1 301h)y2, i11 5 y2, i (26.9b)
Thus, we can see that the problem consists of solving a set of simultaneous equations
for each time step.
For nonlinear ODEs, the solution becomes even more dif! cult since it involves
solving a system of nonlinear simultaneous equations (recall Sec. 6.6). Thus, although
stability is gained through implicit approaches, a price is paid in the form of added solu-
The implicit Euler method is unconditionally stable and only ! rst-order accurate. It
is also possible to develop in a similar manner a second-order accurate implicit trapezoi-
dal rule integration scheme for stiff systems. It is usually desirable to have higher-order
methods. The Adams-Moulton formulas described later in this chapter can also be used
to devise higher-order implicit methods. However, the stability limits of such approaches
are very stringent when applied to stiff systems. Gear (1971) developed a special series
of implicit schemes that have much larger stability limits based on backward difference
formulas. Extensive efforts have been made to develop software to ef! ciently implement
Gear’s methods. As a result, this is probably the most widely used method to solve stiff
systems. In addition, Rosenbrock and others (see Press et al., 2007) have proposed
implicit Runge-Kutta algorithms where the k terms appear implicitly. These methods have
good stability characteristics and are quite suitable for solving systems of stiff ordinary
26.2 MULTISTEP METHODS
The one-step methods described in the previous sections utilize information at a single
point xi to predict a value of the dependent variable yi11 at a future point xi11 (Fig. 26.3a).
Alternative approaches, called multistep methods (Fig. 26.3b), are based on the insight
that, once the computation has begun, valuable information from previous points is at
our command. The curvature of the lines connecting these previous values provides
information regarding the trajectory of the solution. The multistep methods explored in
this chapter exploit this information to solve ODEs. Before describing the higher-order
versions, we will present a simple second-order method that serves to demonstrate the
general characteristics of multistep approaches.
26.2.1 The Non-Self-Starting Heun Method
Recall that the Heun approach uses Euler’s method as a predictor [Eq. (25.15)]:
y0i11 5 yi 1 f(xi, yi)h (26.10)
760 STIFFNESS AND MULTISTEP METHODS
and the trapezoidal rule as a corrector [Eq. (25.16)]:
yi11 5 yi 1 f(xi, yi) 1 f(xi11, y
2 h (26.11)
Thus, the predictor and the corrector have local truncation errors of O(h2) and O(h3),
respectively. This suggests that the predictor is the weak link in the method because it
has the greatest error. This weakness is signi! cant because the ef! ciency of the iterative
corrector step depends on the accuracy of the initial prediction. Consequently, one way
to improve Heun’s method is to develop a predictor that has a local error of O(h3). This
can be accomplished by using Euler’s method and the slope at yi, and extra information
from a previous point yi21, as in
y0i11 5 yi21 1 f(xi, yi)2h (26.12)
Notice that Eq. (26.12) attains O(h3) at the expense of employing a larger step size, 2h. In
addition, note that Eq. (26.12) is not self-starting because it involves a previous value of the
dependent variable yi 2 1. Such a value would not be available in a typical initial-value problem.
Because of this fact, Eqs. (26.11) and (26.12) are called the non-self-starting Heun method.
As depicted in Fig. 26.4, the derivative estimate in Eq. (26.12) is now located at the
midpoint rather than at the beginning of the interval over which the prediction is made.
As demonstrated subsequently, this centering improves the error of the predictor to O(h3).
However, before proceeding to a formal derivation of the non-self-starting Heun, we will
summarize the method and express it using a slightly modi! ed nomenclature:
Predictor: y0i11 5 y m i21 1 f(xi, y
m i )2h (26.13)
Corrector: y j i11 5 y
m i 1
f(xi, y m i ) 1 f(xi11, y
j21 i11 )
(for j 5 1, 2, p , m) (26.14)
FIGURE 26.3 Graphical depiction of the fundamental difference between (a) one-step and (b) multistep methods for solving ODEs.
xxi + 1
xxi + 1xi – 1xi – 2
26.2 MULTISTEP METHODS 761
where the superscripts have been added to denote that the corrector is applied iteratively
from j 5 1 to m to obtain re! ned solutions. Note that ymi and y m i21 are the ! nal results
of the corrector iterations at the previous time steps. The iterations are terminated at any
time step on the basis of the stopping criterion
Zea Z 5 ` y j i11 2 y
y j i11
` 100% (26.15) When ea is less than a prespeci! ed error tolerance es, the iterations are terminated. At this
point, j 5 m. The use of Eqs. (26.13) through (26.15) to solve an ODE is demonstrated
in the following example.
EXAMPLE 26.2 Non-Self-Starting Heun Method
Problem Statement. Use the non-self-starting Heun method to perform the same com- putations as were performed previously in Example 25.5 using Heun’s method. That is,
FIGURE 26.4 A graphical depiction of the non-self-starting Heun method. (a) The midpoint method that is used as a predictor. (b) The trapezoidal rule that is employed as a corrector.
Slope = f (xi+1, yi+1) 0
Slope = f (xi, yi) + f (xi+1, yi+1)
762 STIFFNESS AND MULTISTEP METHODS
integrate y9 5 4e0.8x 2 0.5y from x 5 0 to x 5 4 using a step size of 1.0. As with Example
25.5, the initial condition at x 5 0 is y 5 2. However, because we are now dealing with a
multistep method, we require the additional information that y is equal to 20.3929953 at
x 5 21.
Solution. The predictor [Eq. (26.13)] is used to extrapolate linearly from x 5 21 to x 5 1.
y01 5 20.3929953 1 [4e 0.8(0)
2 0.5(2)] 2 5 5.607005
The corrector [Eq. (26.14)] is then used to compute the value:
y11 5 2 1 4e0.8(0) 2 0.5(2) 1 4e0.8(1) 2 0.5(5.607005)
2 1 5 6.549331
which represents a percent relative error of 25.73 percent (true value 5 6.194631). This
error is somewhat smaller than the value of 28.18 percent incurred in the self-starting Heun.
Now, Eq. (26.14) can be applied iteratively to improve the solution:
y21 5 2 1 3 1 4e0.8(1) 2 0.5(6.549331)
2 1 5 6.313749
which represents an et of 21.92%. An approximate estimate of the error can also be
determined using Eq. (26.15):
0 ea 0 5 ` 6.313749 2 6.549331 6.313749
` 100% 5 3.7% Equation (26.14) can be applied iteratively until ea falls below a prespeci! ed value of
es. As was the case with the Heun method (recall Example 25.5), the iterations converge
on a value of 6.360865 (et 5 22.68%). However, because the initial predictor value is
more accurate, the multistep method converges at a somewhat faster rate.
For the second step, the predictor is
y02 5 2 1 [4e 0.8(1)
2 0.5(6.360865) ] 2 5 13.44346 et 5 9.43%
which is superior to the prediction of 12.08260 (et 5 18%) that was computed with the
original Heun method. The ! rst corrector yields 15.76693 (et 5 6.8%), and subsequent
iterations converge on the same result as was obtained with the self-starting Heun method:
15.30224 (et 5 23.1%). As with the previous step, the rate of convergence of the corrector
is somewhat improved because of the better initial prediction.
Derivation and Error Analysis of Predictor-Corrector Formulas. We have just em- ployed graphical concepts to derive the non-self-starting Heun. We will now show how
the same equations can be derived mathematically. This derivation is particularly interest-
ing because it ties together ideas from curve ! tting, numerical integration, and ODEs.
The exercise is also useful because it provides a simple procedure for developing higher-
order multistep methods and estimating their errors.
The derivation is based on solving the general ODE
dx 5 f(x, y)
26.2 MULTISTEP METHODS 763
This equation can be solved by multiplying both sides by dx and integrating between
limits at i and i 1 1:
dy 5 #
f(x, y) dx
The left side can be integrated and evaluated using [recall Eq. (25.21)]:
yi11 5 yi 1 # xi11
f(x, y) dx (26.16)
Equation (26.16) represents a solution to the ODE if the integral can be evaluated.
That is, it provides a means to compute a new value of the dependent variable yi11 on
the basis of a prior value yi and the differential equation.
Numerical integration formulas such as those developed in Chap. 21 provide one
way to make this evaluation. For example, the trapezoidal rule [Eq. (21.3)] can be used
to evaluate the integral, as in
f(x, y) dx 5
f(xi, yi) 1 f(xi11, yi11)
2 h (26.17)
where h 5 xi1 1 2 xi is the step size. Substituting Eq. (26.17) into Eq. (26.16) yields
yi11 5 yi 1 f(xi, yi) 1 f(xi11, yi11)
which is the corrector equation for the Heun method. Because this equation is based on
the trapezoidal rule, the truncation error can be taken directly from Table 21.2,
Ec 5 2 1
12 h3y(3)(jc) 5 2
12 h3f –(jc) (26.18)
where the subscript c designates that this is the error of the corrector.
A similar approach can be used to derive the predictor. For this case, the integration
limits are from i 2 1 to i 1 1:
dy 5 #
f(x, y) dx
which can be integrated and rearranged to yield
yi11 5 yi21 1 # xi11
f(x, y) dx (26.19)
Now, rather than using a closed formula from Table 21.2, the ! rst Newton-Cotes open
integration formula (see Table 21.4) can be used to evaluate the integral, as in
xi 2 1
f(x, y) dx 5 2h f(xi, yi) (26.20)
which is called the midpoint method. Substituting Eq. (26.20) into Eq. (26.19) yields
yi11 5 yi21 1 2h f(xi, yi)
764 STIFFNESS AND MULTISTEP METHODS
which is the predictor for the non-self-starting Heun. As with the corrector, the local
truncation error can be taken directly from Table 21.4:
Ep 5 1
3 h3y(3)(jp) 5
3 h3f –(jp) (26.21)
where the subscript p designates that this is the error of the predictor.
Thus, the predictor and the corrector for the non-self-starting Heun method have
truncation errors of the same order. Aside from upgrading the accuracy of the predic-
tor, this fact has additional bene! ts related to error analysis, as elaborated in the next
Error Estimates. If the predictor and the corrector of a multistep method are of the same order, the local truncation error may be estimated during the course of a computa-
tion. This is a tremendous advantage because it establishes a criterion for adjustment of
the step size.
The local truncation error for the predictor is estimated by Eq. (26.21). This error
estimate can be combined with the estimate of yi1l from the predictor step to yield [recall
our basic de! nition of Eq. (3.1)]
True value 5 y0i11 1 1
3 h3y(3)(jp) (26.22)
Using a similar approach, the error estimate for the corrector [Eq. (26.18)] can be com-
bined with the corrector result yi1l to give
True value 5 y mi11 2 1
12 h3y(3)(jc) (26.23)
Equation (26.22) can be subtracted from Eq. (26.23) to yield
0 5 ymi11 2 y 0 i11 2
12 h3y(3)(j) (26.24)
where j is now between xi2l and xi1l. Now, dividing Eq. (26.24) by 5 and rearranging
the result gives
y0i11 2 y m i11
5 5 2
12 h3y(3)(j) (26.25)
Notice that the right-hand sides of Eqs. (26.18) and (26.25) are identical, with the excep-
tion of the argument of the third derivative. If the third derivative does not vary appre-
ciably over the interval in question, we can assume that the right-hand sides are equal,
and therefore, the left-hand sides should also be equivalent, as in
Ec 5 2 y0i11 2 y
Thus, we have arrived at a relationship that can be used to estimate the per-step truncation
error on the basis of two quantities—the predictor (y0i11) and the corrector (y m i11)—that
are routine by-products of the computation.
26.2 MULTISTEP METHODS 765
EXAMPLE 26.3 Estimate of Per-Step Truncation Error
Problem Statement. Use Eq. (26.26) to estimate the per-step truncation error of Example 26.2. Note that the true values at x 5 1 and 2 are 6.194631 and 14.84392,
Solution. At xi1l 5 1, the predictor gives 5.607005 and the corrector yields 6.360865. These values can be substituted into Eq. (26.26) to give
Ec 5 2 6.360865 2 5.607005
5 5 20.1507722
which compares well with the exact error,
Et 5 6.194631 2 6.360865 5 20.1662341
At xi1l 5 2, the predictor gives 13.44346 and the corrector yields 15.30224, which
can be used to compute
Ec 5 2 15.30224 2 13.44346
5 5 20.3717550
which also compares favorably with the exact error, Et 5 14.84392 2 15.30224 5
The ease with which the error can be estimated using Eq. (26.26) provides a ratio-
nal basis for step-size adjustment during the course of a computation. For example, if
Eq. (26.26) indicates that the error is greater than an acceptable level, the step size could
Modifi ers. Before discussing computer algorithms, we must note two other ways in which the non-self-starting Heun method can be made more accurate and ef! cient. First,
you should realize that besides providing a criterion for step-size adjustment, Eq. (26.26)
represents a numerical estimate of the discrepancy between the ! nal corrected value at
each step yi11 and the true value. Thus, it can be added directly to yi11 to re! ne the
ymi11d y m i11 2
ymi11 2 y 0 i11
Equation (26.27) is called a corrector modi! er. (The symbol m is read “is replaced by.”)
The left-hand side is the modi! ed value of ymi11.
A second improvement, one that relates more to program ef! ciency, is a predictor
modi! er, which is designed to adjust the predictor result so that it is closer to the ! nal
convergent value of the corrector. This is advantageous because, as noted previously at
the beginning of this section, the number of iterations of the corrector is highly dependent
on the accuracy of the initial prediction. Consequently, if the prediction is modi! ed
properly, we might reduce the number of iterations required to converge on the ultimate
value of the corrector.
766 STIFFNESS AND MULTISTEP METHODS
Such a modi! er can be derived simply by assuming that the third derivative is
relatively constant from step to step. Therefore, using the result of the previous step at
i, Eq. (26.25) can be solved for
h3y(3)(j) 5 2 12
5 (y0i 2 y
m i ) (26.28)
which, assuming that y(3) (j) > y(3) (jp), can be substituted into Eq. (26.21) to give
Ep 5 4
5 (ymi 2 y
0 i ) (26.29)
which can then be used to modify the predictor result:
y0i11d y 0 i11 1
5 (ymi 2 y
0 i ) (26.30)
EXAMPLE 26.4 Effect of Modifi ers on Predictor-Corrector Results
Problem Statement. Recompute Example 26.3 using both modi! ers.
Solution. As in Example 26.3, the initial predictor result is 5.607005. Because the predictor modi! er [Eq. (26.30)] requires values from a previous iteration, it cannot be
employed to improve this initial result. However, Eq. (26.27) can be used to modify the
corrected value of 6.360865 (et 5 22.684%), as in
ym1 5 6.360865 2 6.360865 2 5.607005
5 5 6.210093
which represents an et 5 20.25%. Thus, the error is reduced over an order of magnitude.
For the next iteration, the predictor [Eq. (26.13)] is used to compute
y02 5 2 1 [4e 0.8(0)
2 0.5(6.210093) ] 2 5 13.59423 et 5 8.42%
which is about half the error of the predictor for the second iteration of Example 26.3,
which was et 5 18.6%. This improvement occurs because we are using a superior
estimate of y (6.210093 as opposed to 6.360865) in the predictor. In other words,
the propagated and global errors are reduced by the inclusion of the corrector
Now because we have information from the prior iteration, Eq. (26.30) can be em-
ployed to modify the predictor, as in
y02 5 13.59423 1 4
5 (6.360865 2 5.607005) 5 14.19732 et 5 24.36%
which, again, halves the error.
This modi! cation has no effect on the ! nal outcome of the subsequent corrector
step. Regardless of whether the unmodi! ed or modi! ed predictors are used, the correc-
tor will ultimately converge on the same answer. However, because the rate or ef! ciency
of convergence depends on the accuracy of the initial prediction, the modi! cation can
reduce the number of iterations required for convergence.
26.2 MULTISTEP METHODS 767
Implementing the corrector yields a result of 15.21178 (et 5 22.48%), which rep-
resents an improvement over Example 26.3 because of the reduction of global error.
Finally, this result can be modi! ed using Eq. (26.27):
ym2 5 15.21178 2 15.21178 2 13.59423
5 5 14.88827 et 5 20.30%
Again, the error has been reduced an order of magnitude.
As in the previous example, the addition of the modi! ers increases both the ef! –
ciency and accuracy of multistep methods. In particular, the corrector modi! er effectively
increases the order of the technique. Thus, the non-self-starting Heun with modi! ers is
third order rather than second order as is the case for the unmodi! ed version. However,
it should be noted that there are situations where the corrector modi! er will affect the
stability of the corrector iteration process. As a consequence, the modi! er is not included
in the algorithm for the non-self-starting Heun delineated in Fig. 26.5. Nevertheless, the
corrector modi! er can still have utility for step-size control, as discussed next.
FIGURE 26.5 The sequence of formulas used to implement the non-self-starting Heun method. Note that the corrector error estimates can be used to modify the corrector. However, because this can affect the corrector’s stability, the modifi er is not included in this algorithm. The corrector error estimate is included because of its utility for step-size adjustment.
y0i11 5 yi m 21 1 f (xi, y i
(Save result as y0i11,u 5 y 0 i11 where the subscript u designates that the variable is unmodifi ed.)
Predictor Modifi er:
y0i11d y 0 i11,u 1
5 (y mi,u 2 y
y ji11 5 y m i 1
f (xi, y m i ) 1 f (xi11, y
2 h (for j 5 1 to maximum iterations m)
Zea Z 5 ` y j i11 2 y
y ji11 ` 100%
(If |ea| . error criterion, set j 5 j 11 and repeat corrector; if ea # error criterion, save result as y i
m 11,u 5 y i
Corrector Error Estimate:
Ec 5 2 1
5 (y mi11,u 2 y
(If computation is to continue, set i 5 i 1 1 and return to predictor.)
768 STIFFNESS AND MULTISTEP METHODS
26.2.2 Step-Size Control and Computer Programs
Constant Step Size. It is relatively simple to develop a constant step-size version of the non-self-starting Heun method. About the only complication is that a one-step method
is required to generate the extra point to start the computation.
Additionally, because a constant step size is employed, a value for h must be chosen
prior to the computation. In general, experience indicates that an optimal step size should
be small enough to ensure convergence within two iterations of the corrector (Hull and
Creemer, 1963). In addition, it must be small enough to yield a suf! ciently small trunca-
tion error. At the same time, the step size should be as large as possible to minimize
run-time cost and round-off error. As with other methods for ODEs, the only practical
way to assess the magnitude of the global error is to compare the results for the same
problem but with a halved step size.
Variable Step Size. Two criteria are typically used to decide whether a change in step size is warranted. First, if Eq. (26.26) is greater than some prespeci! ed error criterion,
the step size is decreased. Second, the step size is chosen so that the convergence criterion
of the corrector is satis! ed in two iterations. This criterion is intended to account for the
trade-off between the rate of convergence and the total number of steps in the calculation.
For smaller values of h, convergence will be more rapid but more steps are required. For
larger h, convergence is slower but fewer steps result. Experience (Hull and Creemer,
1963) suggests that the total steps will be minimized if h is chosen so that the corrector
converges within two iterations. Therefore, if over two iterations are required, the step
size is decreased, and if less than two iterations are required, the step size is increased.
Although the above strategy speci! es when step size modi! cations are in order, it
does not indicate how they should be changed. This is a critical question because mul-
tistep methods by de! nition require several points to compute a new point. Once the step
size is changed, a new set of points must be determined. One approach is to restart the
computation and use the one-step method to generate a new set of starting points.
A more ef! cient strategy that makes use of existing information is to increase and
decrease by doubling and halving the step size. As depicted in Fig. 26.6b, if a suf! cient
number of previous values have been generated, increasing the step size by doubling is
a relatively straightforward task (Fig. 26.6c). All that is necessary is to keep track of
subscripts so that old values of x and y become the appropriate new values. Halving the
step size is somewhat more complicated because some of the new values will be unavail-
able (Fig. 26.6a). However, interpolating polynomials of the type developed in Chap. 18
can be used to determine these intermediate values.
In any event, the decision to incorporate step-size control represents a trade-off
between initial investment in program complexity versus the long-term return because
of increased ef! ciency. Obviously, the magnitude and importance of the problem itself
will have a strong bearing on this trade-off. Fortunately, several software packages and
libraries have multistep routines that you can use to obtain solutions without having to
program them from scratch. We will mention some of these when we review packages
and libraries at the end of Chap. 27.
26.2.3 Integration Formulas
The non-self-starting Heun method is characteristic of most multistep methods. It em-
ploys an open integration formula (the midpoint method) to make an initial estimate.
26.2 MULTISTEP METHODS 769
This predictor step requires a previous data point. Then, a closed integration formula (the
trapezoidal rule) is applied iteratively to improve the solution.
It should be obvious that a strategy for improving multistep methods would be to use
higher-order integration formulas as predictors and correctors. For example, the higher-
order Newton-Cotes formulas developed in Chap. 21 could be used for this purpose.
Before describing these higher-order methods, we will review the most common inte-
gration formulas upon which they are based. As mentioned above, the ! rst of these are the
Newton-Cotes formulas. However, there is a second class called the Adams formulas that
we will also review and that are often preferred. As depicted in Fig. 26.7, the fundamental
difference between the Newton-Cotes and Adams formulas relates to the manner in which
the integral is applied to obtain the solution. As depicted in Fig. 26.7a, the Newton-Cotes
formulas estimate the integral over an interval spanning several points. This integral is then
used to project from the beginning of the interval to the end. In contrast, the Adams for-
mulas (Fig. 26.7b) use a set of points from an interval to estimate the integral solely for
the last segment in the interval. This integral is then used to project across this last segment.
FIGURE 26.6 A plot indicating how a halving-doubling strategy allows the use of (b) previously calculated val- ues for a third-order multistep method. (a) Halving; (c) doubling.
770 STIFFNESS AND MULTISTEP METHODS
Newton-Cotes Formulas. Some of the most common formulas for solving ordinary differential equations are based on ! tting an nth-degree interpolating polynomial to n 1 1
known values of y and then using this equation to compute the integral. As discussed
previously in Chap. 21, the Newton-Cotes integration formulas are based on such an
approach. These formulas are of two types: open and closed forms.
Open Formulas. For n equally spaced data points, the open formulas can be expressed
in the form of a solution of an ODE, as was done previously for Eq. (26.19). The general
equation for this purpose is
yi11 5 yi2n 1 # xi11
fn(x) dx (26.31)
xi + 1 xxixi – 1
xi – 2
yi + 1 = yi – 2 +
xi + 1
xi – 2 f (x, y) dx
xi + 1 xxixi – 1
xi – 2
yi + 1 = yi +
xi + 1
xi f (x, y) dx
FIGURE 26.7 Illustration of the fundamental difference between the Newton-Cotes and Adams integration for- mulas. (a) The Newton-Cotes formulas use a series of points to obtain an integral estimate over a number of segments. The estimate is then used to project across the entire range. (b) The Adams formulas use a series of points to obtain an integral estimate for a single segment. The estimate is then used to project across the segment.
26.2 MULTISTEP METHODS 771
where fn(x) is an nth-order interpolating polynomial. The evaluation of the integral em-
ploys the nth-order Newton-Cotes open integration formula (Table 21.4). For example,
if n 5 1,
yi11 5 yi21 1 2h fi (26.32)
where fi is an abbreviation for f(xi, yi)—that is, the differential equation evaluated at xi
and yi. Equation (26.32) is referred to as the midpoint method and was used previously
as the predictor in the non-self-starting Heun method. For n 5 2,
yi11 5 yi22 1 3h
2 ( fi 1 fi21)
and for n 5 3,
yi11 5 yi23 1 4h
3 (2 fi 2 fi21 1 2 fi22) (26.33)
Equation (26.33) is depicted graphically in Fig. 26.8a.
Closed Formulas. The closed form can be expressed generally as
yi11 5 yi2n11 1 # xi11
fn(x) dx (26.34)
where the integral is approximated by an nth-order Newton-Cotes closed integration
formula (Table 21.2). For example, for n 5 1,
yi11 5 yi 1 h
2 ( fi 1 fi11)
which is equivalent to the trapezoidal rule. For n 5 2,
yi11 5 yi21 1 h
3 ( fi21 1 4fi 1 fi11) (26.35)
which is equivalent to Simpson’s 1y3 rule. Equation (26.35) is depicted in Fig. 26.8b.
Adams Formulas. The other types of integration formulas that can be used to solve ODEs are the Adams formulas. Many popular computer algorithms for the multistep solution
of ODEs are based on these methods.
Open Formulas (Adams-Bashforth). The Adams formulas can be derived in a variety
of ways. One technique is to write a forward Taylor series expansion around xi:
yi11 5 yi 1 fi h 1 f ¿i