Numerical Methods for Engineers
In the � rst chapter of this book, we derived the following equation based on Newton’s
second law to compute the velocity y of a falling parachutist as a function of time t
[recall Eq. (1.9)]:
dt 5 g 2
m y (PT7.1)
where g is the gravitational constant, m is the mass, and c is a drag coef� cient. Such
equations, which are composed of an unknown function and its derivatives, are called
differential equations. Equation (PT7.1) is sometimes referred to as a rate equation
because it expresses the rate of change of a variable as a function of variables and pa-
rameters. Such equations play a fundamental role in engineering because many physical
phenomena are best formulated mathematically in terms of their rate of change.
In Eq. (PT7.1), the quantity being differentiated, y, is called the dependent variable.
The quantity with respect to which y is differentiated, t, is called the independent vari-
able. When the function involves one independent variable, the equation is called an
ordinary differential equation (or ODE). This is in contrast to a partial differential equa-
tion (or PDE) that involves two or more independent variables.
Differential equations are also classi� ed as to their order. For example, Eq. (PT7.1)
is called a ! rst-order equation because the highest derivative is a � rst derivative. A
second-order equation would include a second derivative. For example, the equation
describing the position x of a mass-spring system with damping is the second-order
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
which itself can be differentiated to yield
d 2 x
ORDINARY DIFFERENTIAL EQUATIONS
700 ORDINARY DIFFERENTIAL EQUATIONS
Equations (PT7.3) and (PT7.4) can then be substituted into Eq. (PT7.2) to give
dt 1 cy 1 kx 5 0 (PT7.5)
dt 5 2
cy 1 kx
Thus, Eqs. (PT7.3) and (PT7.6) are a pair of � rst-order equations that are equivalent to
the original second-order equation. Because other nth-order differential equations can be
similarly reduced, this part of our book focuses on the solution of � rst-order equations.
Some of the engineering applications in Chap. 28 deal with the solution of second-order
ODEs by reduction to a pair of � rst-order equations.
PT7.1.1 Noncomputer Methods for Solving ODEs
Without computers, ODEs are usually solved with analytical integration techniques. For
example, Eq. (PT7.1) could be multiplied by dt and integrated to yield
y 5 # ag 2 cm yb dt (PT7.7) The right-hand side of this equation is called an inde! nite integral because the limits of
integration are unspeci� ed. This is in contrast to the de� nite integrals discussed previously
in Part Six [compare Eq. (PT7.7) with Eq. (PT6.6)].
An analytical solution for Eq. (PT7.7) is obtained if the inde� nite integral can be
evaluated exactly in equation form. For example, recall that for the falling parachutist
problem, Eq. (PT7.7) was solved analytically by Eq. (1.10) (assuming y 5 0 at t 5 0):
y(t) 5 gm
c (1 2 e2(cym)t) (1.10)
The mechanics of deriving such analytical solutions will be discussed in Sec. PT7.2. For
the time being, the important fact is that exact solutions for many ODEs of practical
importance are not available. As is true for most situations discussed in other parts of
this book, numerical methods offer the only viable alternative for these cases. Because
these numerical methods usually require computers, engineers in the precomputer era
were somewhat limited in the scope of their investigations.
One very important method that engineers and applied mathematicians developed to
overcome this dilemma was linearization. A linear ordinary differential equation is one
that � ts the general form
1 p 1 a1(x)y¿ 1 a0(x)y 5 f (x) (PT7.8)
where y (n)
is the nth derivative of y with respect to x and the a’s and f ’s are speci� ed
functions of x. This equation is called linear because there are no products or nonlinear
functions of the dependent variable y and its derivatives. The practical importance of
linear ODEs is that they can be solved analytically. In contrast, most nonlinear equations
PT7.1 MOTIVATION 701
cannot be solved exactly. Thus, in the precomputer era, one tactic for solving nonlinear
equations was to linearize them.
A simple example is the application of ODEs to predict the motion of a swinging
pendulum (Fig. PT7.1). In a manner similar to the derivation of the falling parachutist
problem, Newton’s second law can be used to develop the following differential equation
(see Sec. 28.4 for the complete derivation):
d 2 u
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
l u 5 0 (PT7.11)
We have, therefore, transformed Eq. (PT7.9) into a linear form that is easy to solve
Although linearization remains a very valuable tool for engineering problem solving,
there are cases where it cannot be invoked. For example, suppose that we were interested
in studying the behavior of the pendulum for large displacements from equilibrium. In
such instances, numerical methods offer a viable option for obtaining solutions. Today,
the widespread availability of computers places this option within reach of all practicing
PT7.1.2 ODEs and Engineering Practice
The fundamental laws of physics, mechanics, electricity, and thermodynamics are usually
based on empirical observations that explain variations in physical properties and states
of systems. Rather than describing the state of physical systems directly, the laws are
usually couched in terms of spatial and temporal changes.
Several examples are listed in Table PT7.1. These laws de� ne mechanisms of change.
When combined with continuity laws for energy, mass, or momentum, differential equa-
tions result. Subsequent integration of these differential equations results in mathematical
functions that describe the spatial and temporal state of a system in terms of energy,
mass, or velocity variations.
The falling parachutist problem introduced in Chap. 1 is an example of the derivation
of an ordinary differential equation from a fundamental law. Recall that Newton’s second
law was used to develop an ODE describing the rate of change of velocity of a falling
parachutist. By integrating this relationship, we obtained an equation to predict fall veloc-
ity as a function of time (Fig. PT7.2). This equation could be utilized in a number of
different ways, including design purposes.
FIGURE PT7.1 The swinging pedulum.
702 ORDINARY DIFFERENTIAL EQUATIONS
In fact, such mathematical relationships are the basis of the solution for a great
number of engineering problems. However, as described in the previous section, many
of the differential equations of practical signi� cance cannot be solved using the analyti-
cal methods of calculus. Thus, the methods discussed in the following chapters are
extremely important in all � elds of engineering.
TABLE PT7.1 Examples of fundamental laws that are written in terms of the rate of change of variables (t 5 time and x 5 position).
Law Mathematical Expression Variables and Parameters
Newton’s second law Velocity (v), force (F), and of motion mass (m)
Fourier’s heat law Heat fl ux (q), thermal conductivity (k9) and temperature (T)
Fick’s law of diffusion Mass fl ux (J), diffusion coeffi cient (D), and concentration (c)
Faraday’s law Voltage drop (DVL), inductance (L), (voltage drop across and current (i) an inductor)
q 5 2k¿ dT
J 5 2D dc
¢VL 5 L di
F = ma
v = (1 – e– (c/m)t) gm
c vi + 1 = vi + (g – vi ) t
= g – vdv
FIGURE PT7.2 The sequence of events in the application of ODEs for engineering problem solving. The exam- ple shown is the velocity of a falling parachutist.
PT7.2 MATHEMATICAL BACKGROUND 703
PT7.2 MATHEMATICAL BACKGROUND
A solution of an ordinary differential equation is a speci� c function of the independent
variable and parameters that satis� es the original differential equation. To illustrate this
concept, let us start with a given function
y 5 20.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:
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.
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.
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
Before proceeding to numerical methods for solving ordinary differential equations, some
orientation might be helpful. The following material is intended to provide you with an
overview of the material discussed in Part Seven. In addition, we have formulated objec-
tives to focus your studies of the subject area.
PT7.3.1 Scope and Preview
Figure PT7.5 provides an overview of Part Seven. Two broad categories of numerical
methods for initial-value problems will be discussed in this part of this book. One-step
methods, which are covered in Chap. 25, permit the calculation of yi11, given the dif-
ferential equation and yi. Multistep methods, which are covered in Chap. 26, require
additional values of y other than at i.
With all but a minor exception, the one-step methods in Chap. 25 belong to what
are called Runge-Kutta techniques. Although the chapter might have been organized
around this theoretical notion, we have opted for a more graphical, intuitive approach to
introduce the methods. Thus, we begin the chapter with Euler’s method, which has a
very straightforward graphical interpretation. Then, we use visually oriented arguments
706 ORDINARY DIFFERENTIAL EQUATIONS
FIGURE PT7.5 Schematic representation of the organization of Part Seven: Ordinary Differential Equations.
26.2 Multistep methods
PT 7.2 Mathematical background
PT 7.6 Advanced methods
PT 7.5 Important formulas
28.4 Mechanical engineering
27.3 Software packages
PT 7.4 Trade-offs
PT 7.3 Orientation
PT 7.1 Motivation
25.2 Heun and midpoint methods
25.1 Euler’s method
25.4 Systems of
25.5 Adaptive RK
PT7.3 ORIENTATION 707
to develop two improved versions of Euler’s method—the Heun and the midpoint tech-
niques. After this introduction, we formally develop the concept of Runge-Kutta (or RK)
approaches and demonstrate how the foregoing techniques are actually � rst- and second-
order RK methods. This is followed by a discussion of the higher-order RK formulations
that are frequently used for engineering problem solving. In addition, we cover the ap-
plication of one-step methods to systems of ODEs. Finally, the chapter ends with a
discussion of adaptive RK methods that automatically adjust the step size in response to
the truncation error of the computation.
Chapter 26 starts with a description of stiff ODEs. These are both individual and
systems of ODEs that have both fast and slow components to their solution. We intro-
duce the idea of an implicit solution technique as one commonly used remedy for this
Next, we discuss multistep methods. These algorithms retain information of previous
steps to more effectively capture the trajectory of the solution. They also yield the trunca-
tion error estimates that can be used to implement step-size control. In this section, we
initially take a visual, intuitive approach by using a simple method—the non-self-starting
Heun—to introduce all the essential features of the multistep approaches.
In Chap. 27 we turn to boundary-value and eigenvalue problems. For the former,
we introduce both shooting and ! nite-difference methods. For the latter, we discuss sev-
eral approaches, including the polynomial and the power methods. Finally, the chapter
concludes with a description of the application of several software packages and librar-
ies for solution of ODEs and eigenvalues.
Chapter 28 is devoted to applications from all the � elds of engineering. Finally, a
short review section is included at the end of Part Seven. This epilogue summarizes and
compares the important formulas and concepts related to ODEs. The comparison includes
a discussion of trade-offs that are relevant to their implementation in engineering prac-
tice. The epilogue also summarizes important formulas and includes references for
PT7.3.2 Goals and Objectives
Study Objectives. After completing Part Seven, you should have greatly enhanced your capability to confront and solve ordinary differential equations and eigenvalue prob-
lems. General study goals should include mastering the techniques, having the capability
to assess the reliability of the answers, and being able to choose the “best’’ method (or
methods) for any particular problem. In addition to these general objectives, the speci� c
study objectives in Table PT7.2 should be mastered.
Computer Objectives. Algorithms are provided for many of the methods in Part Seven. This information will allow you to expand your software library. For example,
you may � nd it useful from a professional viewpoint to have software that employs the
fourth-order Runge-Kutta method for more than � ve equations and to solve ODEs with
an adaptive step-size approach.
In addition, one of your most important goals should be to master several of the
general-purpose software packages that are widely available. In particular, you should
become adept at using these tools to implement numerical methods for engineering
708 ORDINARY DIFFERENTIAL EQUATIONS
TABLE PT7.2 Specifi c study objectives for Part Seven.
1. Understand the visual representations of Euler’s, Heun’s, and the midpoint methods 2. Know the relationship of Euler’s method to the Taylor series expansion and the insight it provides
regarding the error of the method 3. Understand the difference between local and global truncation errors and how they relate to the
choice of a numerical method for a particular problem 4. Know the order and the step-size dependency of the global truncation errors for all the methods
described in Part Seven; understand how these errors bear on the accuracy of the techniques 5. Understand the basis of predictor-corrector methods; in particular, realize that the effi ciency of the
corrector is highly dependent on the accuracy of the predictor 6. Know the general form of the Runge-Kutta methods; understand the derivation of the second-order
RK method and how it relates to the Taylor series expansion; realize that there are an infi nite number of possible versions for second- and higher-order RK methods
7. Know how to apply any of the RK methods to systems of equations; be able to reduce an nth-order ODE to a system of n fi rst-order ODEs
8. Recognize the type of problem context where step size adjustment is important 9. Understand how adaptive step size control is integrated into a fourth-order RK method 10. Recognize how the combination of slow and fast components makes an equation or a system of
equations stiff 11. Understand the distinction between explicit and implicit solution schemes for ODEs; in particular,
recognize how the latter (1) ameliorates the stiffness problem and (2) complicates the solution mechanics
12. Understand the difference between initial-value and boundary-value problems 13. Know the difference between multistep and one-step methods; realize that all multistep methods are
predictor-correctors but that not all predictor-correctors are multistep methods 14. Understand the connection between integration formulas and predictor-corrector methods 15. Recognize the fundamental difference between Newton-Cotes and Adams integration formulas 16. Know the rationale behind the polynomial and the power methods for determining eigenvalues; in
particular, recognize their strengths and limitations 17. Understand how Hoteller’s defl ation allows the power method to be used to compute intermediate
eigenvalues 18. Know how to use software packages and/or libraries to integrate ODEs and evaluate eigenvalues
C H A P T E R 25
This chapter is devoted to solving ordinary differential equations of the form
dx 5 f (x, y)
In Chap. 1, we used a numerical method to solve such an equation for the velocity of
the falling parachutist. Recall that the method was of the general form
New value 5 old value 1 slope 3 step size
or, in mathematical terms,
yi11 5 yi 1 fh (25.1)
According to this equation, the slope estimate of f is used to extrapolate from an old value
yi to a new value yi 1 1 over a distance h (Fig. 25.1). This formula can be applied step by
step to compute out into the future and, hence, trace out the trajectory of the solution.
FIGURE 25.1 Graphical depiction of a one- step method.
Step size = h
xi xi + 1
yi + 1 = yi + h
710 RUNGE-KUTTA METHODS
All one-step methods can be expressed in this general form, with the only difference
being the manner in which the slope is estimated. As in the falling parachutist problem,
the simplest approach is to use the differential equation to estimate the slope in the form
of the ! rst derivative at xi. In other words, the slope at the beginning of the interval is
taken as an approximation of the average slope over the whole interval. This approach,
called Euler’s method, is discussed in the ! rst part of this chapter. This is followed by
other one-step methods that employ alternative slope estimates that result in more ac-
curate predictions. All these techniques are generally called Runge-Kutta methods.
25.1 EULER’S METHOD
The ! rst derivative provides a direct estimate of the slope at xi (Fig. 25.2):
f 5 f (xi, yi)
where f(xi, yi) is the differential equation evaluated at xi and yi. This estimate can be
substituted into Eq. (25.1):
yi11 5 yi 1 f (xi, yi)h (25.2)
This formula is referred to as Euler’s (or the Euler-Cauchy or the point-slope)
method. A new value of y is predicted using the slope (equal to the ! rst derivative at the
original value of x) to extrapolate linearly over the step size h (Fig. 25.2).
EXAMPLE 25.1 Euler’s Method
Problem Statement. Use Euler’s method to numerically integrate Eq. (PT7.13):
dx 5 22x
3 1 12x
2 2 20x 1 8.5
xxi + 1
FIGURE 25.2 Euler’s method.
25.1 EULER’S METHOD 711
from x 5 0 to x 5 4 with a step size of 0.5. The initial condition at x 5 0 is y 5 1.
Recall that the exact solution is given by Eq. (PT7.16):
y 5 20.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
y(0.5) 5 1.0 1 8.5(0.5) 5 5.25
The true solution at x 5 0.5 is
y 5 20.5(0.5) 4
1 4(0.5) 3
2 10(0.5) 2
1 8.5(0.5) 1 1 5 3.21875
Thus, the error is
Et 5 true 2 approximate 5 3.21875 2 5.25 5 22.03125
or, expressed as percent relative error, et 5 263.1%. For the second step,
y(1) 5 y(0.5) 1 f (0.5, 5.25)0.5
5 5.25 1 [22(0.5) 3
1 12(0.5) 2
2 20(0.5) 1 8.5]0.5
The true solution at x 5 1.0 is 3.0, and therefore, the percent relative error is 295.8%.
The computation is repeated, and the results are compiled in Table 25.1 and Fig. 25.3.
TABLE 25.1 Comparison of true and approximate values of the integral of y9 5 22×3 1 12×2 2 20x 1 8.5, with the initial condition that y 5 1 at x 5 0. The approximate values were computed using Euler’s method with a step size of 0.5. The local error refers to the error incurred over a single step. It is calculated with a Taylor series expansion as in Example 25.2. The global error is the total discrepancy due to past as well as present steps.
Percent Relative Error
x ytrue yEuler Global Local
0.0 1.00000 1.00000 0.5 3.21875 5.25000 263.1 263.1 1.0 3.00000 5.87500 295.8 228.1 1.5 2.21875 5.12500 2131.0 21.4 2.0 2.00000 4.50000 2125.0 20.3 2.5 2.71875 4.75000 274.7 17.2 3.0 4.00000 5.87500 246.9 3.9 3.5 4.71875 7.12500 251.0 211.3 4.0 3.00000 7.00000 2133.3 253.1
712 RUNGE-KUTTA METHODS
Note that, although the computation captures the general trend of the true solution, the
error is considerable. As discussed in the next section, this error can be reduced by using
a smaller step size.
The preceding example uses a simple polynomial for the differential equation to
facilitate the error analyses that follow. Thus,
dx 5 f (x)
Obviously, a more general (and more common) case involves ODEs that depend on both
x and y,
dx 5 f (x, y)
As we progress through this part of the text, our examples will increasingly involve ODEs
that depend on both the independent and the dependent variables.
25.1.1 Error Analysis for Euler’s Method
The numerical solution of ODEs involves two types of error (recall Chaps. 3 and 4):
1. Truncation, or discretization, errors caused by the nature of the techniques employed
to approximate values of y.
FIGURE 25.3 Comparison of the true solution with a numerical solution using Euler’s method for the integral of y9 5 22×3 1 12×2 2 20x 1 8.5 from x 5 0 to x 5 4 with a step size of 0.5. The initial condition at x 5 0 is y 5 1.
h = 0.5
25.1 EULER’S METHOD 713
2. Round-off errors caused by the limited numbers of signi! cant digits that can be
retained by a computer.
The truncation errors are composed of two parts. The ! rst is a local truncation error
that results from an application of the method in question over a single step. The second
is a propagated truncation error that results from the approximations produced during
the previous steps. The sum of the two is the total, or global truncation, error.
Insight into the magnitude and properties of the truncation error can be gained by
deriving Euler’s method directly from the Taylor series expansion. To do this, realize that
the differential equation being integrated will be of the general form
y¿ 5 f (x, y) (25.3)
where y9 5 dyydx and x and y are the independent and the dependent variables, respectively. If the solution—that is, the function describing the behavior of y—has continuous deriva-
tives, it can be represented by a Taylor series expansion about a starting value (xi, yi), as in
[recall Eq. (4.7)]
yi11 5 yi 1 y¿i h 1 y–i
2 1 p 1
y (n) i
n 1 Rn (25.4)
where h 5 xi11 2 xi and Rn 5 the remainder term, de! ned as
Rn 5 y
(n 1 1)! h
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 1 p 1
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 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 …