# Numerical Methods for Engineers

## Numerical Methods for Engineers

PART SEVEN

699

PT7.1 MOTIVATION

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)]:

dy

dt 5 g 2

c

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

equation,

m d 2x

dt 2 1 c

dx

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

dt (PT7.3)

which itself can be differentiated to yield

dy

dt 5

d 2x

dt2 (PT7.4)

ORDINARY DIFFERENTIAL EQUATIONS

700 ORDINARY DIFFERENTIAL EQUATIONS

Equations (PT7.3) and (PT7.4) can then be substituted into Eq. (PT7.2) to give

m dy

dt 1 cy 1 kx 5 0 (PT7.5)

or

dy

dt 5 2

cy 1 kx

m (PT7.6)

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

an(x)y (n)

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):

d 2u

dt 2 1

g

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

d 2u

dt 2 1

g

l u 5 0 (PT7.11)

We have, therefore, transformed Eq. (PT7.9) into a linear form that is easy to solve

analytically.

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

engineers.

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.

u

l

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)

dv

dt 5

F

m

q 5 2k¿ dT

dx

J 5 2D dc

dx

¢VL 5 L di

dt

F = ma

Analytical Numerical

v = (1 – e– (c/m)t) gm

c vi + 1 = vi + (g – vi) t

c

m

= g – vdv

dt

c

m

Physical law

Solution

ODE

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

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:

dy

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.

y

4

(a)

x3

dy/dx

8

(b)

x

– 8

3

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.

y

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

eigenvalues.

PT7.3 ORIENTATION

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.

CHAPTER 25

Runge-Kutta

Methods

PART 7

Ordinary

Differential

Equations

CHAPTER 26

Stiffness/

Multistep

Methods

CHAPTER 27

Boundary Value

and Eigenvalue

Problems

CHAPTER 28

Case Studies

EPILOGUE

26.2 Multistep methods

26.1 Stiffness

PT 7.2 Mathematical background

PT 7.5 Important formulas

28.4 Mechanical engineering

28.3 Electrical

engineering

28.2 Civil

engineering

28.1 Chemical

engineering

27.1 Boundary-

value problems

27.3 Software packages

27.2 Eigenvalues

PT 7.3 Orientation

PT 7.1 Motivation

25.2 Heun and midpoint methods

25.1 Euler’s method

25.3 Runge-Kutta

25.4 Systems of

ODEs

methods

PT7.3 ORIENTATION 707

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

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

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

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

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

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

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

the truncation error of the computation.

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

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

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

problem.

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

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

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

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

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

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

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

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

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

ies for solution of ODEs and eigenvalues.

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

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

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

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

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

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

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

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

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

problem solving.

708 ORDINARY DIFFERENTIAL EQUATIONS

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

25

709

C H A P T E R 25

Runge-Kutta Methods

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

dy

dx 5 f(x, y)

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

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

New value 5 old value 1 slope 3 step size

or, in mathematical terms,

yi11 5 yi 1 fh (25.1)

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

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

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

FIGURE 25.1 Graphical depiction of a one- step method.

y

x

Step size = h

Slope =

xi xi + 1

yi + 1 = yi + h

710 RUNGE-KUTTA METHODS

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

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

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

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

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

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

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

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

25.1 EULER’S METHOD

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

f 5 f(xi, yi)

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

substituted into Eq. (25.1):

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

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

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

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

EXAMPLE 25.1 Euler’s Method

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

dy

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

y

xxi + 1

error

Predicted

True

xi

h

FIGURE 25.2 Euler’s method.

25.1 EULER’S METHOD 711

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

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

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

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

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

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

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

Therefore,

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

The true solution at x 5 0.5 is

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

Thus, the error is

Et 5 true 2 approximate 5 3.21875 2 5.25 5 22.03125

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

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

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

5 5.875

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

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

TABLE 25.1 Comparison of true and approximate values of the integral of y9 5 22×3 1 12×2 2 20x 1 8.5, with the initial condition that y 5 1 at x 5 0. The approximate values were computed using Euler’s method with a step size of 0.5. The local error refers to the error incurred over a single step. It is calculated with a Taylor series expansion as in Example 25.2. The global error is the total discrepancy due to past as well as present steps.

Percent Relative Error

x ytrue yEuler Global Local

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

712 RUNGE-KUTTA METHODS

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

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

a smaller step size.

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

facilitate the error analyses that follow. Thus,

dy

dx 5 f(x)

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

x and y,

dy

dx 5 f(x, y)

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

that depend on both the independent and the dependent variables.

25.1.1 Error Analysis for Euler’s Method

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

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

to approximate values of y.

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.

y

4

0 x4

True solution

h = 0.5

20

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

y(n)i

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

or

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

expansion.

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

this purpose:

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-

order method.

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).

y

4

0 x4

True solution

h = 0.5

h = 0.25

2

(a)

0

y

– 0.5

0 x42

True

Estimated

(b)

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.

1

10

100

0.1 0.1

Step size

A b

s o

lu te

p e rc

e n

t re

la ti

v e e

rr o

r

0.01 0.0011

50

Steps

500 50005

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

‘initialize variables

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

END DO

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

DO

xend 5 x 1 xout

IF (xend . xf) THEN xend 5 xf

h 5 dx

CALL Integrator (x, y, h, xend)

m 5 m 1 1

xpm 5 x

ypm 5 y

IF (x \$ xf) EXIT

END DO

DISPLAY RESULTS

END

(b) Routine to Take One Output Step

SUB Integrator (x, y, h, xend)

DO

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

CALL Euler (x, y, h, ynew)

y 5 ynew

IF (x \$ xend) EXIT

END DO

END SUB

(c) Euler’s Method for a Single ODE

SUB Euler (x, y, h, ynew)

CALL Derivs (x, y, dydx)

ynew 5 y 1 dydx * h

x 5 x 1 h

END SUB

(d) Routine to Determine Derivative

SUB Derivs (x, y, dydx)

dydx 5 . . .

END SUB

720 RUNGE-KUTTA METHODS

EXAMPLE 25.4 Solving ODEs with the Computer

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

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

based on Newton’s second law in the form

dy

dt 5 g 2

c

m y (E25.4.1)

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

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

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

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

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

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

dy

dt 5 g 2

c

m cy 1 a a y

ymax bb d (E25.4.2)

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

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

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

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

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

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

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

solution.

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

60

y

40

20

0 t15

Nonlinear

Linear

5 100

25.2 IMPROVEMENTS OF EULER’S METHOD 721

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

solution of the linear model for comparison purposes.

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

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

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

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

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

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

aspects of the problem rather than to tedious manual computations.

25.1.3 Higher-Order Taylor Series Methods

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

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

term from Eq. (25.6) yields

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

2! h2 (25.11)

with a local truncation error of

Ea 5 f –(xi, yi)

6 h3

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

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

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

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

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

0x 1 0f(x, y)

0y dy

dx

The second derivative is

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

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

0y dy

dx

Higher-order derivatives become increasingly more complicated.

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

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

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

25.2 IMPROVEMENTS OF EULER’S METHOD

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

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

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

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

722 RUNGE-KUTTA METHODS

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

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

25.2.1 Heun’s Method

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

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

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

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

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

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

is used to extrapolate linearly to yi11:

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

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

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

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

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

y

xxi + 1xi

(a)

Slope = f (xi, yi)

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

y

xxi + 1xi

(b)

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

2

0

25.2 IMPROVEMENTS OF EULER’S METHOD 723

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

mated slope at the end of the interval:

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

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

slope for the interval:

y¿ 5 y¿i 1 y¿i11

2 5

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

2

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

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

0 i11)

2 h

which is called a corrector equation.

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

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

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

expressed concisely as

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

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

0 i11)

2 h (25.16)

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

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

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

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

f(x i , y

i ) + f (x

i+ 1 , y

i+ 1 )

y i + h

y i +

1

2

0

724 RUNGE-KUTTA METHODS

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

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

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

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

Zea Z 5 ` y j i11 2 y

j21 i11

y j i11

` 100% (25.17) where y

j21 i11 and y

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

tor, respectively.

EXAMPLE 25.5 Heun’s Method

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

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

y 5 4

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

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

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

y¿0 5 4e 0

2 0.5(2) 5 3

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

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

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

estimate of y at 1.0:

y01 5 2 1 3(1) 5 5

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

Iterations of Heun’s Method

1 15

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

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

25.2 IMPROVEMENTS OF EULER’S METHOD 725

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

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

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

the end of the interval

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

0.8(1) 2 0.5(5) 5 6.402164

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

from x 5 0 to 1

y¿ 5 3 1 6.402164

2 5 4.701082

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

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

y1 5 2 1 4.701082(1) 5 6.701082

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

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

as compared with Euler’s method.

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

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

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

2 1 5 6.275811

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

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

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

2 1 5 6.382129

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

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

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

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

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

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

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

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

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

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

is expressed concisely as

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

2 h (25.18)

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

dy

dx 5 f(x)

This equation can be solved for y by integration:

# yi11

yi

dy 5 #

xi11

xi

f(x) dx (25.19)

which yields

yi11 2 yi 5 # xi11

xi

f(x) dx (25.20)

or

yi11 5 yi 1 # xi11

xi

f(x) dx (25.21)

Now, recall from Chap. 21 that the trapezoidal rule [Eq. (21.3)] is de! ned as

# xi11

xi

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

2 (25.25)

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.

y

5

x

True solution

Euler’s method

Heun’s method

3

y

xxi + 1xi

y Slope = f (xi + 1/2, yi + 1/2)

xxi + 1/2xi

(b)

(a)

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

# b

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

# xi11

xi

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.

25.2.4 Summary

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

techniques.

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

END SUB

(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

END SUB

(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

DO

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

END DO

ynew 5 ye

x 5 x 1 h

END SUB

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)

where

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

2 (25.32)

a2q11 5 1

2 (25.33)

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

2a2 (25.35)

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)

where

k1 5 f (xi, yi) (B25.1.2)

and

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

[Eq. (25.11)]

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

(Sec. 25.1.3):

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

0x 1 0f (x, y)

0y dy

dx (B25.1.5)

Substituting Eq. (B25.1.5) into (B25.1.4) gives

yi11 5 yi 1 f (xi, yi)h 1 a 0f 0x

1 0f

0y dy

dx b h2

2! (B25.1.6)

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

Eq. (4.26)]

g(x 1 r, y 1 s) 5 g(x, y) 1 r 0g

0x 1 s 0g

0y 1p

Applying this method to expand Eq. (B25.1.3) gives

f (xi 1 p1h, yi 1 q11k1h) 5 f (xi, yi) 1 p1h 0f

0x

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

0f

0x

1 a2q11h 2 f (xi, yi)

0f

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)

(B25.1.7)

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

must hold:

a1 1 a2 5 1

a2 p1 5 1

2

a2q11 5 1

2

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)

where

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)

where

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)

where

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)

to predict

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.

y

4

0 x420

Analytical Euler Heun Midpoint Ralston

734 RUNGE-KUTTA METHODS

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

k2 5 22(0.375) 3

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

The average slope is computed by

f 5 1

3 (8.5) 1

2

3 (2.58203125) 5 4.5546875

which can be used to predict

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

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

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

25.3.2 Third-Order Runge-Kutta Methods

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

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

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

parameters. One common version that results is

yi11 5 yi 1 1

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

where

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

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

Second-Order Heun Midpoint Ralston RK

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

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

25.3 RUNGE-KUTTA METHODS 735

k2 5 f axi 1 1 2

h, yi 1 1

2 k1hb (25.39b)

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

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

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

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

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

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

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

25.3.3 Fourth-Order Runge-Kutta Methods

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

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

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

yi11 5 yi 1 1

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

where

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

k2 5 f axi 1 1 2

h, yi 1 1

2 k1hb (25.40b)

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

y

xxi+1/2

h

xi

k2

k1

k3

k3

k2

k1

k4

xi+1

736 RUNGE-KUTTA METHODS

k3 5 f axi 1 1 2

h, yi 1 1

2 k2hb (25.40c)

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

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

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

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

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

of these to arrive at the improved slope.

EXAMPLE 25.7 Classical Fourth-Order RK Method

Problem Statement.

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

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

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

(b) Similarly, integrate

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

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

Solution.

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

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

y(0.5) 5 1 1 e 1 6

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

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

order method gives an exact result.

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

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

2 0.5(2) 5 3

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

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

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

2 0.5(2.75) 5 3.510611

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

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

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

2 0.5(2.877653) 5 3.446785

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

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

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

2 0.5(3.723392) 5 4.105603

25.3 RUNGE-KUTTA METHODS 737

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

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

f 5 1

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

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

which compares favorably with the true solution of 3.751521.

25.3.4 Higher-Order Runge-Kutta Methods

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

recommended:

yi11 5 yi 1 1

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

where

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

k2 5 f axi 1 1 4

h, yi 1 1

4 k1hb (25.41b)

k3 5 f axi 1 1 4

h, yi 1 1

8 k1h 1

1

8 k2hb (25.41c)

k4 5 f axi 1 1 2

h, yi 2 1

2 k2h 1 k3hb (25.41d)

k5 5 f axi 1 3 4

h, yi 1 3

16 k1h 1

9

16 k4hb (25.41e)

k6 5 f axi 1 h, yi 2 3 7

k1h 1 2

7 k2h 1

12

7 k3h 2

12

7 k4h 1

8

7 k5hb (25.41f)

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

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

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

EXAMPLE 25.8 Comparison of Runge-Kutta Methods

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

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

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

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

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

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

738 RUNGE-KUTTA METHODS

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

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

Effort 5 nf b 2 a

h (E25.8.1)

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

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

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

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

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

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

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

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

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

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

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

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

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

100

1

10– 2

10– 4

10– 6

Euler

Heun

RK–3d

RK–4th

Butcher

Effort

P e rc

e n

t re

la ti

v e e

rr o

r

25.4 SYSTEMS OF EQUATIONS 739

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

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

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

25.3.5 Computer Algorithms for Runge-Kutta Methods

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

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

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

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

25.4 SYSTEMS OF EQUATIONS

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

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

may be represented generally as

dy1

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

dy2

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

#

#

#

dyn

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

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

value of x.

SUB RK4 (x, y, h, ynew)

CALL Derivs(x, y, k1)

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

CALL Derivs(x 1 h, ye, k4)

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

x 5 x 1 h

END SUB

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

740 RUNGE-KUTTA METHODS

25.4.1 Euler’s Method

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

system shown above. Engineering applications can involve thousands of simultaneous

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

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

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

EXAMPLE 25.9 Solving Systems of ODEs Using Euler’s Method

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

of 0.5.

dy1

dx 5 20.5y1

dy2

dx 5 4 2 0.3y2 2 0.1y1

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

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

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

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

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

x y1 y2

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

25.4.2 Runge-Kutta Methods

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

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

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

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

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

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

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

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

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

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

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

following example illustrates the approach.

25.4 SYSTEMS OF EQUATIONS 741

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

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

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

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

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

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

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

y1 1 k1,1 h

2 5 4 1 (22)

0.5

2 5 3.5

y2 1 k1, 2 h

2 5 6 1 (1.8)

0.5

2 5 6.45

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

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

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

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

y1 1 k2,1 h

2 5 4 1 (21.75)

0.5

2 5 3.5625

y2 1 k2, 2 h

2 5 6 1 (1.715)

0.5

2 5 6.42875

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

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

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

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

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

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

which can be used to compute the endpoint slopes,

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

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

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

y1(0.5) 5 4 1 1

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

y2(0.5) 5 6 1 1

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

742 RUNGE-KUTTA METHODS

Proceeding in a like manner for the remaining steps yields

x y1 y2

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

25.4.3 Computer Algorithm for Solving Systems of ODEs

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

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

1. Inputting the number of equations, n.

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

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

variables.

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

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

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

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

to the fact that

1. There are n equations.

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

EXAMPLE 25.11 Solving Systems of ODEs with the Computer

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

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

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

dy1

dx 5 y2

dy2

dx 5 216.1y1

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

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

dy3

dx 5 y4

dy4

dx 5 216.1 sin(y3)

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

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

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

25.4 SYSTEMS OF EQUATIONS 743

(a) Main or “Driver” Program

Assign values for

n 5 number of equations

yi 5 initial values of n dependent

variables

xi 5 initial value independent

variable

xf 5 final value independent variable

dx 5 calculation step size

xout 5 output interval

x 5 xi

m 5 0

xpm 5 x

DOFOR i 5 1, n

ypi,m 5 yii

yi 5 yii

END DO

DO

xend 5 x 1 xout

IF (xend . xf) THEN xend 5 xf

h 5 dx

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

m 5 m 1 1

xpm 5 x

DOFOR i 5 1, n

ypi,m 5 yi

END DO

IF (x \$ xf) EXIT

END DO

DISPLAY RESULTS

END

(b) Routine to Take One Output Step

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

DO

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

CALL RK4 (x, y, n, h)

IF (x \$ xend) EXIT

END DO

END SUB

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

SUB RK4 (x, y, n, h)

CALL Derivs (x, y, k1)

DOFOR i 5 1, n

ymi 5 yi 1 k1i * h / 2

END DO

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

DOFOR i 5 1, n

ymi 5 yi 1 k2i * h / 2

END DO

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

DOFOR i 5 1, n

yei 5 yi 1 k3i * h

END DO

CALL Derivs (x 1 h, ye, k4)

DOFOR i 5 1, n

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

yi 5 yi 1 slopei * h

END DO

x 5 x 1 h

END SUB

(d ) Routine to Determine Derivatives

SUB Derivs (x, y, dy)

dy1 5 …

dy2 5 …

END SUB

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

744 RUNGE-KUTTA METHODS

Solution.

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

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

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

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

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

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

4

2

0y

0 321

x

– 4

– 2

4

y1, y3

y2, y4

(a)

4

2

0y

0 2 31

x

– 4

– 2

4

y2 y4

y3

y1

(b)

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

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

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

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

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

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

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

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

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

I 5 # b

a f(x) dx

is equivalent to solving the differential equation

dy

dx 5 f(x)

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

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

smooth but exhibit regions of abrupt change.

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

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

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

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

1

0 1 2 3

y

x

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 ¢

15 (25.44)

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

and

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.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

125

594 k4 1

512

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

13,525

55,296 k4 1

277

14,336 k5 1

1

4 k6b h (25.46)

where

k1 5 f(xi, yi)

k2 5 f axi 1 1 5

h, yi 1 1

5 k1hb

k3 5 f axi 1 3 10

h, yi 1 3

40 k1h 1

9

40 k2hb

k4 5 f axi 1 3 5

h, yi 1 3

10 k1h 2

9

10 k2h 1

6

5 k3hb

k5 5 f axi 1 h, yi 2 11 54

k1h 1 5

2 k2h 2

70

27 k3h 1

35

27 k4hb

k6 5 f axi 1 7 8

h, yi 1 1631

55,296 k1h 1

175

512 k2h 1

575

13,824 k3h 1

44,275

110,592 k4h

1 253

4096 k5hb

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

125

594 6.832587 1

512

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

13,525

55,296 6.832587

1 227

14,336 12.09831 1

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.

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

SUBROUTINE RKkc (y,dy,x,h,yout,yerr)

PARAMETER (a250.2,a350.3,a450.6,a551.,a650.875,

b2150.2,b3153.y40.,b3259.y40.,b4150.3,b42520.9,

b4351.2,b515211.y54.,b5252.5,b535270.y27.,

b54535.y27.,b6151631.y55296.,b625175.y512.,

b635575.y13824.,b64544275.y110592.,b655253.y4096.,

c1537.y378.,c35250.y621.,c45125.y594.,

c65512.y1771.,dc15c122825.y27648.,

dc35c3218575.y48384.,dc45c4213525.y55296.,

dc552277.y14336.,dc65c620.25)

ytemp5y1b21*h*dy

CALL Derivs (x1a2*h,ytemp,k2)

ytemp5y1h*(b31*dy1b32*k2)

CALL Derivs(x1a3*h,ytemp,k3)

ytemp5y1h*(b41*dy1b42*k21b43*k3)

CALL Derivs(x1a4*h,ytemp,k4)

ytemp5y1h*(b51*dy1b52*k21b53*k31b54*k4)

CALL Derivs(x1a5*h,ytemp,k5)

ytemp5y1h*(b61*dy1b62*k21b63*k31b64*k41b65*k5)

CALL Derivs(x1a6*h,ytemp,k6)

yout5y1h*(c1*dy1c3*k31c4*k41c6*k6)

yerr5h*(dc1*dy1dc3*k31dc4*k41dc5*k51dc6*k6)

END RKkc

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

dy

dx 1 0.6y 5 10e2(x22)

2y[2(0.075)2] (E25.14.1)

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

maxstep5100

hi5.5; tiny 5 1. 3 10230

eps50.00005

print *, xi,yi

x5xi

y5yi

h5hi

istep50

DO

IF (istep . maxstep AND x # xf) EXIT

istep5istep11

CALL Derivs(x,y,dy)

yscal5ABS(y)1ABS(h*dy)1tiny

IF (x1h.xf) THEN h5xf2x

PRINT x,y

h5hnxt

END DO

END

PARAMETER (safety50.9, econ51.89e24)

h5htry

DO

CALL RKkc (y,dy,x,h,ytemp,yerr)

emax5abs(yerr/yscal/eps)

IF emax # 1 EXIT

htemp5safety*h*emax 20.25

h5max(abs(htemp),0.25*abs(h))

xnew5x1h

IF xnew5x THEN pause

END DO

IF emax . econ THEN

hnxt5safety*emax 2.2 *h

ELSE

hnxt54.*h

END IF

x5x1h

y5ytemp

FIGURE 25.22 Pseudocode for a (a) driver program and an (b) adaptive step routine to solve a single ODE.

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

1

2

0 2 4 x

(b)

0

5

10

0 2 4 x

(a)

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

PROBLEMS

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.

dy

dt 5 yt 2 2 1.1y

(a) Analytically.

(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.

dy

dt 5 (1 1 4t)1y

(a) Analytically.

(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

solve

d 2y

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:

d 2y

dx 2 1 0.6

dy

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.

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:

dy

dt 5 y sin3(t)  y(0) 5 1

25.6 Solve the following problem numerically from t 5 0 to 3:

dy

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

dy

dt 5 22y 1 5e2t

dz

dt 5 2

yz2

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

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

Table 25.2.

25.14 Develop a user-friendly computer program for the classical

fourth-order RK method. Test the program by duplicating Exam-

ple 25.7.

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

dx

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.

FIGURE P25.16

k

c

x

m

PROBLEMS 753

25.21 The logistic model is used to simulate population as in

dp

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,

dy

dt 5 2g(0)

R2

(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

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:

dy

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

equation:

d 2x

dt 2 1 (5x)

dx

dt 1 (x 1 7) sin(vt) 5 0

where

dx

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:

dy

dt 5 g 2

cd

m y2

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.

H

r

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:

d 2y

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:

dy

dt 5 y

dv

dt 5 2g 2

cd

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

2 rACd

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

755

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-

size control.

26.1 STIFFNESS

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

ODE is

dy

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),

dy

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

dt h

Substituting Eq. (26.3) gives

yi11 5 yi 2 ayih

or

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.

3

y

2

1

0 42 t0

1

0 0.020.010

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

dt h

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

dy

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.

Solution.

(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

2ti11

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.

1.5

y

1

0.5

0 0.0060.004

h = 0.0015

h = 0.0005

Exact

(a)

t0 0.002

2

y

1

0 0.40.3

Exact

h = 0.05

(b)

t0 0.20.1

Systems of ODEs can also be stiff. An example is

dy1

dt 5 25y1 1 3y2 (26.6a)

dy2

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-

tion complexity.

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

differential equations.

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

0 i11)

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 )

2 h

(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.

y

xi

(a)

xxi + 1

y

xi

(b)

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

j21 i11

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.

y

xxi+1

xi–1

xi

(a)

(b)

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

y

xxi+1xi

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

2

0

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

dy

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:

# yi11

yi

dy 5 #

xi11

xi

f(x, y) dx

The left side can be integrated and evaluated using [recall Eq. (25.21)]:

yi11 5 yi 1 # xi11

xi

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

# xi11

xi

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)

2 h

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

1

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:

# yi11

yi21

dy 5 #

xi11

xi21

f(x, y) dx

which can be integrated and rearranged to yield

yi11 5 yi21 1 # xi11

xi21

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

# xi11

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

1

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

section.

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

5

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

1

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

m i11

5 (26.26)

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,

respectively.

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

20.4583148.

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

be decreased.

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

estimate further:

ymi11d y m i11 2

ymi11 2 y 0 i11

5 (26.27)

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

4

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

modifier.

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.

Predictor:

y0i11 5 yi m 21 1 f (xi, y i

m)2h

(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

4

5 (y mi,u 2 y

0 i,u)

Corrector:

y ji11 5 y m i 1

f (xi, y m i ) 1 f (xi11, y

j21 i11)

2 h (for j 5 1 to maximum iterations m)

Error Check:

Zea Z 5 ` y j i11 2 y

j21 i11

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

m 11.)

Corrector Error Estimate:

Ec 5 2 1

5 (y mi11,u 2 y

0 i11,u)

(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.

y

x

Interpolation

(a)

y

x

(b)

y

x

(c)

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

xi2n

fn(x) dx (26.31)

y

xi + 1 xxixi – 1

(a)

xi – 2

yi + 1 = yi – 2 +

xi + 1

xi – 2 f (x, y) dx

y

xi + 1 xxixi – 1

(b)

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

xi2n11

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