Numerical Methods for Engineers
In the rst chapter of this book, we derived the following equation based on Newton’s
second law to compute the velocity y of a falling parachutist as a function of time t
[recall Eq. (1.9)]:
dt 5 g 2
m y (PT7.1)
where g is the gravitational constant, m is the mass, and c is a drag coef cient. Such
equations, which are composed of an unknown function and its derivatives, are called
differential equations. Equation (PT7.1) is sometimes referred to as a rate equation
because it expresses the rate of change of a variable as a function of variables and pa-
rameters. Such equations play a fundamental role in engineering because many physical
phenomena are best formulated mathematically in terms of their rate of change.
In Eq. (PT7.1), the quantity being differentiated, y, is called the dependent variable.
The quantity with respect to which y is differentiated, t, is called the independent vari-
able. When the function involves one independent variable, the equation is called an
ordinary differential equation (or ODE). This is in contrast to a partial differential equa-
tion (or PDE) that involves two or more independent variables.
Differential equations are also classi ed as to their order. For example, Eq. (PT7.1)
is called a ! rst-order equation because the highest derivative is a rst derivative. A
second-order equation would include a second derivative. For example, the equation
describing the position x of a mass-spring system with damping is the second-order
m d 2x
dt 2 1 c
dt 1 kx 5 0 (PT7.2)
where c is a damping coef cient and k is a spring constant. Similarly, an nth-order equa-
tion would include an nth derivative.
Higher-order equations can be reduced to a system of rst-order equations. For Eq.
(PT7.2), this is done by de ning a new variable y, where
y 5 dx
which itself can be differentiated to yield
ORDINARY DIFFERENTIAL EQUATIONS
700 ORDINARY DIFFERENTIAL EQUATIONS
Equations (PT7.3) and (PT7.4) can then be substituted into Eq. (PT7.2) to give
dt 1 cy 1 kx 5 0 (PT7.5)
dt 5 2
cy 1 kx
Thus, Eqs. (PT7.3) and (PT7.6) are a pair of rst-order equations that are equivalent to
the original second-order equation. Because other nth-order differential equations can be
similarly reduced, this part of our book focuses on the solution of rst-order equations.
Some of the engineering applications in Chap. 28 deal with the solution of second-order
ODEs by reduction to a pair of rst-order equations.
PT7.1.1 Noncomputer Methods for Solving ODEs
Without computers, ODEs are usually solved with analytical integration techniques. For
example, Eq. (PT7.1) could be multiplied by dt and integrated to yield
y 5 # ag 2 cm yb dt (PT7.7) The right-hand side of this equation is called an inde! nite integral because the limits of
integration are unspeci ed. This is in contrast to the de nite integrals discussed previously
in Part Six [compare Eq. (PT7.7) with Eq. (PT6.6)].
An analytical solution for Eq. (PT7.7) is obtained if the inde nite integral can be
evaluated exactly in equation form. For example, recall that for the falling parachutist
problem, Eq. (PT7.7) was solved analytically by Eq. (1.10) (assuming y 5 0 at t 5 0):
y(t) 5 gm
c (1 2 e2(cym)t) (1.10)
The mechanics of deriving such analytical solutions will be discussed in Sec. PT7.2. For
the time being, the important fact is that exact solutions for many ODEs of practical
importance are not available. As is true for most situations discussed in other parts of
this book, numerical methods offer the only viable alternative for these cases. Because
these numerical methods usually require computers, engineers in the precomputer era
were somewhat limited in the scope of their investigations.
One very important method that engineers and applied mathematicians developed to
overcome this dilemma was linearization. A linear ordinary differential equation is one
that ts the general form
1 p 1 a1(x)y¿ 1 a0(x)y 5 f(x) (PT7.8)
where y(n) is the nth derivative of y with respect to x and the a’s and f ’s are speci ed
functions of x. This equation is called linear because there are no products or nonlinear
functions of the dependent variable y and its derivatives. The practical importance of
linear ODEs is that they can be solved analytically. In contrast, most nonlinear equations
PT7.1 MOTIVATION 701
cannot be solved exactly. Thus, in the precomputer era, one tactic for solving nonlinear
equations was to linearize them.
A simple example is the application of ODEs to predict the motion of a swinging
pendulum (Fig. PT7.1). In a manner similar to the derivation of the falling parachutist
problem, Newton’s second law can be used to develop the following differential equation
(see Sec. 28.4 for the complete derivation):
dt 2 1
l sin u 5 0 (PT7.9)
where u is the angle of displacement of the pendulum, g is the gravitational constant,
and l is the pendulum length. This equation is nonlinear because of the term sin u. One
way to obtain an analytical solution is to realize that for small displacements of the
pendulum from equilibrium (that is, for small values of u),
sin u > u (PT7.10)
Thus, if it is assumed that we are interested only in cases where u is small, Eq. (PT7.10)
can be substituted into Eq. (PT7.9) to give
dt 2 1
l u 5 0 (PT7.11)
We have, therefore, transformed Eq. (PT7.9) into a linear form that is easy to solve
Although linearization remains a very valuable tool for engineering problem solving,
there are cases where it cannot be invoked. For example, suppose that we were interested
in studying the behavior of the pendulum for large displacements from equilibrium. In
such instances, numerical methods offer a viable option for obtaining solutions. Today,
the widespread availability of computers places this option within reach of all practicing
PT7.1.2 ODEs and Engineering Practice
The fundamental laws of physics, mechanics, electricity, and thermodynamics are usually
based on empirical observations that explain variations in physical properties and states
of systems. Rather than describing the state of physical systems directly, the laws are
usually couched in terms of spatial and temporal changes.
Several examples are listed in Table PT7.1. These laws de ne mechanisms of change.
When combined with continuity laws for energy, mass, or momentum, differential equa-
tions result. Subsequent integration of these differential equations results in mathematical
functions that describe the spatial and temporal state of a system in terms of energy,
mass, or velocity variations.
The falling parachutist problem introduced in Chap. 1 is an example of the derivation
of an ordinary differential equation from a fundamental law. Recall that Newton’s second
law was used to develop an ODE describing the rate of change of velocity of a falling
parachutist. By integrating this relationship, we obtained an equation to predict fall veloc-
ity as a function of time (Fig. PT7.2). This equation could be utilized in a number of
different ways, including design purposes.
FIGURE PT7.1 The swinging pedulum.
702 ORDINARY DIFFERENTIAL EQUATIONS
In fact, such mathematical relationships are the basis of the solution for a great
number of engineering problems. However, as described in the previous section, many
of the differential equations of practical signi cance cannot be solved using the analyti-
cal methods of calculus. Thus, the methods discussed in the following chapters are
extremely important in all elds of engineering.
TABLE PT7.1 Examples of fundamental laws that are written in terms of the rate of change of variables (t 5 time and x 5 position).
Law Mathematical Expression Variables and Parameters
Newton’s second law Velocity (v), force (F), and of motion mass (m)
Fourier’s heat law Heat fl ux (q), thermal conductivity (k9) and temperature (T)
Fick’s law of diffusion Mass fl ux ( J), diffusion coeffi cient (D), and concentration (c)
Faraday’s law Voltage drop (DV L), inductance (L), (voltage drop across and current (i) an inductor)
q 5 2k¿ dT
J 5 2D dc
¢VL 5 L di
F = ma
v = (1 – e– (c/m)t) gm
c vi + 1 = vi + (g – vi) t
= g – vdv
FIGURE PT7.2 The sequence of events in the application of ODEs for engineering problem solving. The exam- ple shown is the velocity of a falling parachutist.
PT7.2 MATHEMATICAL BACKGROUND 703
PT7.2 MATHEMATICAL BACKGROUND
A solution of an ordinary differential equation is a speci c function of the independent
variable and parameters that satis es the original differential equation. To illustrate this
concept, let us start with a given function
y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1 (PT7.12)
which is a fourth-order polynomial (Fig. PT7.3a). Now, if we differentiate Eq. (PT7.12),
we obtain an ODE:
dx 5 22×3 1 12×2 2 20x 1 8.5 (PT7.13)
This equation also describes the behavior of the polynomial, but in a manner different
from Eq. (PT7.12). Rather than explicitly representing the values of y for each value of
x, Eq. (PT7.13) gives the rate of change of y with respect to x (that is, the slope) at every
value of x. Figure PT7.3 shows both the function and the derivative plotted versus x. Notice
FIGURE PT7.3 Plots of (a) y versus x and (b) dy/dx versus x for the function y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1.
704 ORDINARY DIFFERENTIAL EQUATIONS
how the zero values of the derivatives correspond to the point at which the original func-
tion is ! at—that is, has a zero slope. Also, the maximum absolute values of the derivatives
are at the ends of the interval where the slopes of the function are greatest.
Although, as just demonstrated, we can determine a differential equation given the
original function, the object here is to determine the original function given the differ-
ential equation. The original function then represents the solution. For the present case,
we can determine this solution analytically by integrating Eq. (PT7.13):
y 5 # (22×3 1 12×2 2 20x 1 8.5) dx Applying the integration rule (recall Table PT6.2)
#un du 5 u n11
n 1 1 1 C n ? 21
to each term of the equation gives the solution
y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 C (PT7.14)
which is identical to the original function with one notable exception. In the course of
differentiating and then integrating, we lost the constant value of 1 in the original equa-
tion and gained the value C. This C is called a constant of integration. The fact that such
an arbitrary constant appears indicates that the solution is not unique. In fact, it is but
one of an in nite number of possible functions (corresponding to an in nite number of
possible values of C) that satisfy the differential equation. For example, Fig. PT7.4 shows
six possible functions that satisfy Eq. (PT7.14).
FIGURE PT7.4 Six possible solutions for the integral of 22×3 1 12×2 2 20x 1 8.5. Each conforms to a different value of the constant of integration C.
x C = 0
C = – 1
C = – 2
C = 3
C = 2
C = 1
PT7.3 ORIENTATION 705
Therefore, to specify the solution completely, a differential equation is usually ac-
companied by auxiliary conditions. For rst-order ODEs, a type of auxiliary condition
called an initial value is required to determine the constant and obtain a unique solution.
For example, Eq. (PT7.13) could be accompanied by the initial condition that at x 5 0,
y 5 1. These values could be substituted into Eq. (PT7.14):
1 5 20.5(0)4 1 4(0)3 2 10(0)2 1 8.5(0) 1 C (PT7.15)
to determine C 5 1. Therefore, the unique solution that satis es both the differential
equation and the speci ed initial condition is obtained by substituting C 5 1 into Eq.
(PT7.14) to yield
y 5 20.5×4 1 4×3 2 10×2 1 8.5x 1 1 (PT7.16)
Thus, we have “pinned down’’ Eq. (PT7.14) by forcing it to pass through the initial
condition, and in so doing, we have developed a unique solution to the ODE and have
come full circle to the original function [Eq. (PT7.12)].
Initial conditions usually have very tangible interpretations for differential equations
derived from physical problem settings. For example, in the falling parachutist problem,
the initial condition was re! ective of the physical fact that at time zero the vertical veloc-
ity was zero. If the parachutist had already been in vertical motion at time zero, the
solution would have been modi ed to account for this initial velocity.
When dealing with an nth-order differential equation, n conditions are required to
obtain a unique solution. If all conditions are speci ed at the same value of the indepen-
dent variable (for example, at x or t 5 0), then the problem is called an initial-value
problem. This is in contrast to boundary-value problems where speci cation of conditions
occurs at different values of the independent variable. Chapters 25 and 26 will focus on
initial-value problems. Boundary-value problems are covered in Chap. 27 along with
Before proceeding to numerical methods for solving ordinary differential equations, some
orientation might be helpful. The following material is intended to provide you with an
overview of the material discussed in Part Seven. In addition, we have formulated objec-
tives to focus your studies of the subject area.
PT7.3.1 Scope and Preview
Figure PT7.5 provides an overview of Part Seven. Two broad categories of numerical
methods for initial-value problems will be discussed in this part of this book. One-step
methods, which are covered in Chap. 25, permit the calculation of yi11, given the dif-
ferential equation and yi. Multistep methods, which are covered in Chap. 26, require
additional values of y other than at i.
With all but a minor exception, the one-step methods in Chap. 25 belong to what
are called Runge-Kutta techniques. Although the chapter might have been organized
around this theoretical notion, we have opted for a more graphical, intuitive approach to
introduce the methods. Thus, we begin the chapter with Euler’s method, which has a
very straightforward graphical interpretation. Then, we use visually oriented arguments
706 ORDINARY DIFFERENTIAL EQUATIONS
FIGURE PT7.5 Schematic representation of the organization of Part Seven: Ordinary Differential Equations.
26.2 Multistep methods
PT 7.2 Mathematical background
PT 7.6 Advanced methods
PT 7.5 Important formulas
28.4 Mechanical engineering
27.3 Software packages
PT 7.4 Trade-offs
PT 7.3 Orientation
PT 7.1 Motivation
25.2 Heun and midpoint methods
25.1 Euler’s method
25.4 Systems of
25.5 Adaptive RK