# ENGINEERING

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

2 x

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 2 x

dt 2

(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 2 u

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 2 u

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

concept, let us start with a given function

y 5 20.5x 4

1 4x 3

2 10x 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 22x

3 1 12x

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 n 1 1

n 1 1 1 C n ? 21

to each term of the equation gives the solution

y 5 20.5x 4

1 4x 3

2 10x 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.5x 4

1 4x 3

2 10x 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.6 Advanced methods

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.4 Trade-offs

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

25.5 Adaptive RK

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

advanced topics.

PT7.3.2 Goals and Objectives

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

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

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

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

study objectives in Table PT7.2 should be mastered.

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

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

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

an adaptive step-size approach.

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

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

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

problem solving.

708 ORDINARY DIFFERENTIAL EQUATIONS

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

25

709

C H A P T E R 25

Runge-Kutta Methods

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

dy

dx 5 f (x, y)

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

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

New value 5 old value 1 slope 3 step size

or, in mathematical terms,

yi11 5 yi 1 fh (25.1)

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

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

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

FIGURE 25.1 Graphical depiction of a one- step method.

y

x

Step size = h

Slope =

xi xi + 1

yi + 1 = yi + h

710 RUNGE-KUTTA METHODS

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

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

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

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

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

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

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

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

25.1 EULER’S METHOD

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

f 5 f (xi, yi)

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

substituted into Eq. (25.1):

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

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

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

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

EXAMPLE 25.1 Euler’s Method

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

dy

dx 5 22x

3 1 12x

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.5x 4

1 4x 3

2 10x 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! h

2 1 p 1

y (n) i

n! h

n 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)! h

n11 (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! h

2 1 p 1

f (n21)

(xi, yi)

n! h

n 1 O(hn11) (25.6)

where O(h n11

) 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! h

2 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 …