6_nome_e_numero_Chapra_Canale_1998_Numerical_Differentiation_and_Integration

Prévia do material em texto

PART S 
NUMER 
Calculus is the mathematics of change. Because engineers must continuously deal with sys- 
tems and processes that change, calculus is an essential tool of our profession. Standing at the 
heart of calculus are the related mathematical concepts of differentiation and integration. 
According to the dictionary definition, to dlfferentiate means "to mark off by differ- 
ences; distinguish; . . . to perceive the difference in or between." Mathematically, the deriva- 
tive, which serves as the fundamental vehicle for differentiation, represents the rate of change 
of a dependent variable with respect to an independent variable. As depicted in Fig. PT6.1, 
the mathematical definition of the derivative begins with a difference approximation: 
where y and j ( x ) are alternative representatives for the dependent variable and x is the 
independent variable. If Ax is allowed to approach zero, as occurs in moving from 
Fig. PT6.la to c, the difference becomes a derivative 
d y - = lim f(xr + Ax) - f(xi> 
dx AX-o AX 
FlGURE PT6.1 
The graphical definition of a derivative: as Ax approaches zero in going from (a] to (c], the 
difference approximation becomes a derivative. 
568 NUMERICAL DIFFERENTIATION AND INTEGRATION 
FlGURE PT6.2 
Graphical representation of the integral of f(x\ between the limits x = a to b. The integral is 
equivalent to the area under the curve. 
where dyldx [which can also be designated as y' or f '(xi)] is the first derivative of y with re- 
spect to x evaluated at xi. As seen in the visual depiction of Fig. PT6. lc, the derivative is 
the slope of the tangent to the curve at xi. 
The inverse process to differentiation in calculus is integration. According to the dic- 
tionary definition, to integrate means "to bring together, as parts, into a whole; to unite; to 
indicate the total amount . . . ." Mathematically, integration is represented by 
which stands for the integral of the function f (x) with respect to the independent variable x, 
evaluated between the limits x = a to x = b. The function f (x ) in Eq. (PT6.2) is referred to 
as the integrand. 
As suggested by the dictionary definition, the "meaning" of Eq. (PT6.2) is the total value, 
or summation, off (x) dx over the rangex = a to b. In fact, the symbol is actually a stylized 
capital S that is intended to signify the close connection between integration and summation. 
Figure PT6.2 represents a graphical manifestation of the concept. For functions lying 
above the x axis, the integral expressed by Eq. (PT6.2) corresponds to the area under the 
curve of f (x ) between x = a and b.' 
As outlined above, the "marking off or discrimination" of differentiation and the 
"bringing together" of integration are closely linked processes that are, in fact, inversely 
'It should be noted that the process represented by Eq. (PT6.2) and Fig. PT6.2 is called definite integration. There 
is another type called indefinite integration in which the limits a and b are unspecified. As will be discussed in 
Part Seven, indefinite integration deals with determining a function whose derivative is given. 
PT6.1 MOTIVATION 569 
FIGURE PT6.3 
The contrast between (a) differ- 
entiation and (b) integration. 
related (Fig. PT6.3). For example, if we are given a function y( t ) that specifies an object's 
position as a function of time, differentiation provides a means to determine its velocity, as 
in (Fig. PT6.3a). 
Conversely, if we are provided with velocity as a function of time, integration can be used 
to determine its position (Fig. PT6.3b), 
Thus, we can make the general claim that the evaluation of the integral 
is equivalent to solving the differential equation 
for y(b) given the initial condition y(a) = 0. 
Because of this close relationship, we have opted to devote this part of the book to both 
processes. Among other things, this will provide the opportunity to highlight their similar- 
ities and differences from a numerical perspective. In addition, our discussion will have 
relevance to the next parts of the book, where we will cover differential equations. 
NUMERICAL DIFFERENTIATION AND INTEGRATION 
PT6.1.1 Noncornputer Methods For Differentiation and Infegratisn 
The function to be differentiated or integrated will typically be in one of the following three 
forms: 
1. A simple continuous function such as a polynomial, an exponential, or a trigonometric 
function. 
2. A complicated continuous function that is difficult or impossible to differentiate or in- 
tegrate directly. 
3. A tabulated function where values of x and f(x) are given at a number of discrete points, 
as is often the case with experimental or field data. 
In the first case, the derivative or integral of a simple function may be evaluated ana- 
lytically using calculus. For the second case, analytical solutions are often impractical, and 
sometimes impossible, to obtain. In these instances, as well as in the third case of discrete 
data, approximate methods must be employed. 
A noncomputer method for determining derivatives from data is called equal-area 
graphical dzfferentiation. In this method, the (x, y) data are tabulated and, for each interval, 
a simple divided difference AylAx is employed to estimate the slope. Then these values 
are plotted as a stepped curve versus x (Fig. PT6.4). Next, a smooth curve is drawn that at- 
tempts to approximate the area under the stepped curve. That is, it is drawn so that visually, 
the positive and negative areas are balanced. The rates at given values of x can then be read 
from the curve. 
In the same spirit, visually oriented approaches were employed to integrate tabulated 
data and complicated functions in the precomputer era. A simple intuitive approach is to 
plot the function on a grid (Fig. PT6.5) and count the number of boxes that approximate the 
area. This number multiplied by the area of each box provides a rough estimate of the total 
FIGURE PT6.4 
Equal-area differentiation. 
(a) Centered finite divided 
differences are used to estimate 
the derivative for each interval 
between the data points. 
(b) The derivative estimates are 
plotted as a bar graph. A 
smooth curve is superimposed 
on this plot to a the 
area under the E ar graph. This 
is accomplished by drawing the 
curve so that equal positive and 
negative areas are balanced. 
(c) Values of dyldx can then be 
read off the smooth curve. 
PT6.1 MOTIVATION 571 
FIGURE PT6.5 
The use of a grid to approxi- 
mate an integral. 
FlGURE PT6.6 
The use of rectangles, or strips, 
to approximate the integral. 
area under the curve. This estimate can be refined, at the expense of additional effort, by 
using a finer grid. 
Another commonsense approach is to divide the area into vertical segments, or strips, 
with a height equal to the function value at the midpoint of each strip (Fig. PT6.6). The area 
of the rectangles can then be calculated and summed to estimate the total area. In this 
572 NUMERICAL DIFFERENTIATION AND INTEGRATION 
approach, it is assumed that the value at the midpoint provides a valid approximation of the 
average height of the function for each strip. As with the grid method, refined estimates are 
possible by using more (and thinner) strips to approximate the integral. 
Although such simple approaches have utility for quick estimates, alternative numeri- 
cal techniques are available for the same purpose. Not surprisingly, the simplest of these 
methods is similar in spirit to the noncomputer techniques. 
For differentiation, the most fundamental numerical techniques use finite divided dif- 
ferences to estimate derivatives. For data with error, an alternative approach is to fit a 
smooth curve to the data with a technique such as least-squares regression and then differ- 
entiate this curve to obtain derivative estimates. 
In a similar spirit, numerical integration or quadrature methods are available to obtain 
integrals. These methods, which are actually easier to implementthan the grid approach, 
are similar in spirit to the strip method. That is, function heights are multiplied by strip 
widths and summed to estimate the integral. However, through clever choices of weighting 
factors, the resulting estimate can be made more accurate than that from the simple "strip 
method." 
As in the simple strip method, numerical integration and differentiation techniques uti- 
lize data at discrete points. Because tabulated information is already in such a form, it is 
naturally compatible with many of the numerical approaches. Although continuous func- 
tions are not originally in discrete form, it is usually a simple proposition to use the given 
equation to generate a table of values. As depicted in Fig. PT6.7, this table can then be eval- 
uated with a numerical method. 
PT6.1.2 Numerical Differentiation and Integration in Engineering 
The differentiation and integration of a function has so many engineering applications that 
you were required to take differential and integral calculus in your first year at college. 
Many specific examples of such applications could be given in all fields of engineering. 
Differentiation is commonplace in engineering because so much of our work involves 
characterizing the changes of variables in both time and space. In fact, many of the laws 
and other generalizations that figure so prominently in our work are based on the pre- 
dictable ways in which change manifests itself in the physical world. A prime example is 
Newton's second law, which is not couched in terms of the position of an object but rather 
in its change of position with respect to time. 
Aside from such temporal examples, numerous laws governing the spatial behavior of 
variables are expressed in terms of derivatives. Among the most common of these are those 
laws involving potentials or gradients. For example, Fourier's law of heat conduction 
quantifies the observation that heat flows from regions of high to low temperature. For the 
one-dimensional case, this can be expressed mathematically as 
dT 
Heat flux = -k' - 
dx 
Thus, the derivative provides a measure of the intensity of the temperature change, or gra- 
dient, that drives the transfer of heat. Similar laws provide workable models in many other 
areas of engineering, including the modeling of fluid dynamics, mass transfer, chemical re- 
action kinetics, and electromagnetic flux. The ability to accurately estimate derivatives is 
an important facet of our capability to work effectively in these areas. 
PT6.1 MOTIVATION 573 
FlGURE BT6.7 
Application of a nurner~cal inte- 
gration method: (a) A compli- 
cated, cont~nuous function. 
(b) Table of discrete values of 
f ix) generated from the function. 
[c) Use of a numerical method 
(the strip method here) to esti- 
mate the integral on the basis of 
the discrete points. For a tabu- 
lated function, the data is al- 
ready in tabular form (b); there- 
fore, step (a) 1s unnecessary. 
Just as accurate estimates of derivatives are important in engineering, the calculation 
of integrals is equally valuable. A number of examples relate directly to the idea of the in- 
tegral as the area under a curve. Figure PT6.8 depicts a few cases where integration is used 
for this purpose. 
Other common applications relate to the analogy between integration and summation. 
For example, a common application is to determine the mean of continuous functions. In 
Part Five, you were introduced to the concept of the mean of n discrete data points [recall 
Eq. (PT6.1)]: 
2 Yi 
/ = I Mean = - 
n 
where y; are individual measurements. The determination of the mean of discrete points is 
depicted in Fig. PT6.9a. 
574 NUMERICAL DIFFERENTIATION AND INTEGRATION 
FIGURE PT6.8 
Examples of how integration is used to evaluate areas in engineering applications. (a] A sur- 
veyor might need to know the area of a field bounded by a meandering stream and two roads 
[b] A water-resource engineer might need to know the cross-sectional area of a river. (c) A struc- 
tural engineer might need to determine the net force due to a nonuniform wind blowing against 
the side of a skyscraper. 
FlGURE BT6.9 
An illustration of the mean for 
(a) discrete and (b) continuous 
data 
PT6.1 MOTIVATION 575 
In contrast, suppose that y is a continuous function of an independent variable x, as de- 
picted in Fig. PT6.9b. For this case, there are an infinite number of values between a and 
b. Just as Eq. (PT6.3) can be applied to determine the mean of the discrete readings, you 
might also be interested in computing the mean or average of the continuous function 
y = f (x ) for the interval from a to 6. Integration is used for this purpose, as specified by the 
formula 
l b f ( x ) dx 
Mean = (PT6.4) 
b - a 
This formula has hundreds of engineering applications. For example, it is used to calculate 
the center of gravity of irregular objects in mechanical and civil engineering and to deter- 
mine the root-mean-square current in electrical engineering. 
Integrals are also employed by engineers to evaluate the total amount or quantity of a 
given physical variable. The integral may be evaluated over a line, an area, or a volume. 
For example, the total mass of chemical contained in a reactor is given as the product of the 
concentration of chemical and the reactor volume, or 
Mass = concentration x volume 
where concentration has units of mass per volume. However, suppose that concentration 
varies from location to location within the reactor. In this case, it is necessary to sum the 
products of local concentrations ci and corresponding elemental volumes AVi: 
n 
Mass = ci AVi 
1=1 
where n is the number of discrete volumes. For the continuous case, where c(x, y, z) is a 
known function and x, y, and z are independent variables designating position in Cartesian 
coordinates, integration can be used for the same purpose: 
Mass = c ( x , y , z ) d x d y dz 
which is referred to as a volume integral. Notice the strong analogy between summation 
and integration. 
Similar examples could be given in other fields of engineering. For example, the total 
rate of energy transfer across a plane where the flux (in calories per square centimeter per 
second) is a function of position is given by 
which is referred to as an areal integral, where A = area. 
576 NUMERICAL DIFFERENTIATION AND INTEGRATION 
Similarly, for the one-dimensional case, the total mass of a variable-density rod with 
constant cross-sectional area is given by 
where m = total weight (kg), L = length of the rod (m), p(x) = known density (kg/m3) as 
a function of length x (m), and A = cross-sectional area of the rod (m2). 
Finally, integrals are used to evaluate differential or rate equations. Suppose the ve- 
locity of a particle is a known continuous function of time v(t), 
The total distance y traveled by this particle over a time t is given by (Fig. PT6.3b) 
(PT6.5) 
These are just a few of the applications of differentiation and integration that you might 
face regularly in the pursuit of your profession. When the functions to be analyzed are sim- 
ple, you will normally choose to evaluate them analytically. For example, in the falling 
parachutist problem, we determined the solution for velocity as a function of time 
[Eq. (1.10)]. This relationship could be substituted into Eq. (PT6.5), which could then be 
integrated easily to determine how far the parachutist fell over a time period t. For this case, 
the integral is simple to evaluate. However, it is difficult or impossible when the function 
is complicated, as is typically the case in more realistic examples. In addition, the underly- 
ing function is often unknown and defined only by measurement at discrete points. For 
both these cases, you must have the ability to obtain approximate values for derivatives and 
integrals using numerical techniques. Several such techniques will be discussed in this part 
of the book. 
BT6.2 MATHEMATICAL BACKGROUNDIn high school or during your first year of college, you were introduced to differential and 
integral calculus. There you learned techniques to obtain analytical or exact derivatives 
and integrals. 
When we differentiate a function analytically, we generate a second function that can 
be used to compute the derivative for different values of the independent variable. General 
rules are available for this purpose. For example, in the case of the monomial 
the following simple rule applies (n # 0): 
which is the expression of the more general rule for 
PT6.2 MATHEMATICAL BACKGROUND 597 
where u = a function of x. For this equation, the derivative is computed via 
Two other useful formulas apply to the products and quotients of functions. For example, 
if the product of two functions of x(u and v) is represented as y = u v, then the derivative 
can be computed as 
For the division, y = ulv, the derivative can be computed as 
Other useful formulas are summarized in Table PT6.1. 
Similar formulas are available for definite integration, which deals with determining 
an integral between specified limits, as in 
According to the fundamental theorem of integral calculus, Eq. (PT6.6) is evaluated as 
where F(x) = the integral of f(x)-that is, any function such that F'(x) = f(x). The nomen- 
clature on the right-hand side stands for 
TABLE F"T6.1 Some commonly used derivatives. 
d 
- sin x = cos x 
d 
- cot x = -csc2 x 
dx dx 
d 
- cos x = -sin x 
d 
- sec x = sec x tan x 
dx dx 
d 
- tan x = sec2 x 
dx 
d 
- csc x = -csc x cot x 
dx 
d 1 
- logo x = -- 
dx x l n a 
578 NUMERICAL DIFFERENTIATION AND INTEGRATION 
An example of a definite integral is 
~ 0 . 8 
For this case, the function is a simple polynomial that can be integrated analytically by 
evaluating each term according to the rule 
b 
xn dx = - (PT6.9) 
where n cannot equal -1. Applying this rule to each term in Eq. (PT6.8) yields 
which can be evaluated according to Eq. (PT6.7) as I = 1.6405333. This value is equal to 
the area under the original polynomial [Eq. (PT6.8)] between x = 0 and 0.8. 
The foregoing integration depends on knowledge of the rule expressed by Eq. (PT6.9). 
Other functions follow different rules. These "rules" are all merely instances of antidiffer- 
entiation, that is, finding F(x) so that Ff(x) = f (x ) . Consequently, analytical integration de- 
pends on prior knowledge of the answer. Such knowledge is acquired by training and 
TABLE PT6.2 Some simple integrals that are used in Part Six. The a and b in this table 
are constants and should not be confused with the limits of integration 
discussed in the text. 
J u d v = u v - J v d u 
1 I cos [ax + b) dx= - sin (OX + b) + C 
dx 1 6 6 I,, = a ton-' y x + C 
PT6.3 ORIENTATION 579 
experience. Many of the rules are summarized in handbooks and in tables of integrals. We 
list some commonly encountered integrals in Table PT6.2. However, many functions of 
practical importance are too complicated to be contained in such tables. One reason why 
the techniques in the present part of the book are so valuable is that they provide a means 
to evaluate relationships such as Eq. (PT6.8) without knowledge of the rules. 
PT6.3 ORIENTATION 
Before proceeding to the numerical methods for integration, some further orientation might 
be helpful. The following is intended as an overview of the material discussed in Part Six. 
In addition, we have formulated some objectives to help focus your efforts when studying 
the material. 
PT6.3.1 Scope and Preview 
Figure PT6.10 provides an overview of Part Six. Chapter 21 is devoted to the most com- 
mon approaches for numerical integration-the Newton-Cotes ,formulas. These relation- 
ships are based on replacing a complicated function or tabulated data with a simple poly- 
nomial that is easy to integrate. Three of the most widely used Newton-Cotes formulas are 
discussed in detail: the trapezoidal rule, Simpson's 1 / 3 rule, and Simpson's 318 rule. All 
these formulas are designed for cases where the data to be integrated is evenly spaced. In 
addition, we also include a discussion of numerical integration of unequally spaced data. 
This is a very important topic because many real-world applications deal with data that is 
in this form. 
All the above material relates to closed integration, where the function values at the 
ends of the limits of integration are known. At the end of Chap. 21, we present open inte- 
gration formulas, where the integration limits extend beyond the range of the known data. 
Although they are not commonly used for definite integration, open integration formulas 
are presented here because they are utilized extensively in the solution of ordinary differ- 
ential equations in Part Seven. 
The formulations covered in Chap. 21 can be employed to analyze both tabulated data 
and equations. Chapter 22 deals with two techniques that are expressly designed to inte- 
grate equations and functions: Romherg integration and Gauss quad]-ature. Computer al- 
gorithms are provided for both of these methods. In addition, methods for evaluating im- 
proper integrals are discussed. 
In Chap. 23, we present additional information on numerical di'erentiation to supple- 
ment the introductory material from Chap. 4. Topics include high-accuracy finite-difference 
formulas, Richardson's extrapolation, and the differentiation of unequally spaced data. The 
effect of errors on both numerical differentiation and integration is discussed. Finally, the 
chapter is concluded with a description of the application of several software packages and 
libraries for integration and differentiation. 
Chapter 24 demonstrates how the methods can be applied for problem solving. As 
with other parts of the book, applications are drawn from all fields of engineering. 
A review section, or epilogue, is included at the end of Part Six. This review includes 
a discussion of trade-offs that are relevant to implementation in engineering practice. In ad- 
dition, important formulas are summarized. Finally, we present a short review of advanced 
580 NUMERICAL DIFFERENTIATION AND INTEGRATION 
FIGURE PT6.10 
Schematic of the organization of material in Part Six: Numerical Integration and Differentiation. 
methods and alternative references that will facilitate your further studies of numerical 
differentiation and integration. 
PT6.3.2 Goals and Objectives 
Study Objectives. After completing Part Six you should be able to solve many numeri- 
cal integration and differentiation problems and appreciate their application for engineer- 
PT6.3 ORIENTATION 58 l 
TABLE PT6.3 Specific study objectives for Part Six. 
Understand the derivation of the Newton-Cotes formulas; know how to derve the trapezoidal rule and 
how to set up the derivation of both of Simpson's rules; recognize that the trapezoidal and Simpson's 
1 /3 and 3/8 rules represent the areas under flrst-, second-, and third-order polynom~als, respectively 
Know the formulas and error equations for la) the trapezoidal rule, Ib) the multiple-application trape- 
zoidal rule, (c) Slmpson's 1 /3 rule, (dl Simpson's 3/8 rule, and (e) the multiple-application Simpson's 
rule Be able to choose the "best" among these formulas for any particular problem context. 
Recognize that Slrnpson's 1 /3 rule is fourth-order accurate even though lt IS based on only three points; 
realize that all the even-segment-odd-point Newton-Cotes formulas have sim~lar enhanced accuracy 
Know how to evaluate the ~ntegral and derivative of unequally spaced data. 
Recognize the difference between open and closed integration formulas. 
Understand the theoret~cal basis of Richardson extrapolation and how ~t IS applied in the Romberg ~nte- 
gration algorithm and for numerical differentiation 
Understand the fundamental difference between Newton-Cotes and Gauss quadrature formulas. 
Recognize why both Romberg integration and Gauss quadrature have utility when ~ntegrating equa- 
tions(as opposed to tabular or discrete data). 
Know how open integration formulas are employed to evaluate improper integrals, 
Understand the application of high-accuracy numerical-differentiat~on formulas 
Know how to differentiate unequally spaced data. 
Recognize the differing effects of data error on the processes of numerical integration and differentiation 
ing problem solving. You should strive to master several techniques and assess their relia- 
bility. You should understand the trade-offs involved in selecting the "best" method (or 
methods) for any particular problem. In addition to these general objectives, the specific 
concepts listed in Table PT6.3 should be assimilated and mastered. 
Computer Objectives. You have been provided with software and simple computer al- 
gorithms to implement the techniques discussed in Part Six. All have utility as learning 
tools. 
The Numerical Methods TOOLKIT personal computer software is user-friendly. It 
employs the trapezoidal rule to evaluate the integral of either continuous functions or tab- 
ular data. The graphics associated with this software will enable you to easily visualize 
your problem and the associated mathematical operations as the area between the curve 
and the x axis. The software is very easy to apply to solve many practical problems and can 
be used to check the results of any computer programs you may develop yourself. 
In addition, algorithms are provided for most of the other methods in Part Five. This 
information will allow you to expand your software library to include techniques beyond 
the trapezoidal rule. For example, you may find it useful from a professional viewpoint to 
have software to implement numerical integration and differentiation of unequally spaced 
data. You may also want to develop your own software for Simpson's rules, Romberg inte- 
gration, and Gauss quadrature, which are usually more efficient and accurate than the 
trapezoidal rule. 
Finally, 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. 
CHAPTER 2 
The Newton-Cotes formulas are the most common numerical integration schemes. They 
are based on the strategy of replacing a complicated function or tabulated data with an 
approximating function that is easy to integrate: 
where fn(x) = a polynomial of the form 
where n is the order of the polynomial. For example, in Fig. 21. la , a first-order polynomial 
(a straight line) is used as an approximation. In Fig. 21. lh, a parabola is employed for the 
same purpose. 
The integral can also be approximated using a series of polynomials applied piecewise 
to the function or data over segments of constant length. For example, in Fig. 21.2, three 
FlGURE 2 1.1 
The approximat~on of an inte- 
gral by the area under (a) a sin- 
gle straight line and (b) a single 
parabola 
NEWTON-COTES INTEGRATION FORMULAS 583 
straight-line segments are used to approximate the integral. Higher-order polynomials can 
be utilized for the same purpose. With this background, we now recognize that the "strip 
method" in Fig. PT6.6 employed a series of zero-order polynomials (that is, constants) to 
approximate the integral. 
Closed and open forms of the Newton-Cotes formulas are available. The closed 
forms are those where the data points at the beginning and end of the limits of integration 
are known (Fig. 21 .3~) . The open forms have integration limits that extend beyond the 
range of the data (Fig. 21.3b). In this sense, they are akin to extrapolation as discussed in 
Sec. 18.5. Open Newton-Cotes formulas are not generally used for definite integration. 
FIGURE 2 1.2 
The approximation of an inte- 
gral by the area under three 
straight-line segments. 
FIGURE 2 1.3 
The difference between 
(a] closed and (b) open integra- 
tion formulas. 
584 NEWTON-COTES INTEGRATION FORMULAS 
However, they are utilized for evaluating improper integrals and for the solution of ordi- 
nary differential equations. This chapter emphasizes the closed forms. However, material 
on open Newton-Cotes formulas is briefly introduced at the end of this chapter. 
2 1.1 THE TRAPEZOIDAL RULE 
The trapezoidal rule is the first of the Newton-Cotes closed integration formulas. It corre- 
sponds to the case where the polynomial in Eq. (21.1) is first-order: 
b b 
f (x) dx 2 f 1 ( x ) dx 
Recall from Chap. 18 that a straight line can be represented as [Eq. (18.2)] 
f l(x) = f(a) + f(b) - f(a) (x 
b - a 
- a) 
The area under this straight line is an estimate of the integral off (x ) between the limits a 
and b: 
f (b) - f (a) (x - ) 
b - a 
The result of the integration (see Box 21.1 for details) is 
which is called the trapezoidal rule. 
Before integration, Eq. (21.2) can be expressed as This result can be evaluated to give 
f (b ) - f W X + f ( a ) - 
f l ( x ) = 
af (b ) - af (a ) 
b - a b - a 
I = f (b) - f (a ) (b2 - a2) + bf (a ) - af (b ) (b - a ) 
b - a 2 b - a 
Grouping the last two terms gives Now, since b2 - a2 = (b - a)(b + a), 
f1(x) = f (b ) - f ( a ) x + bf (a ) - a f ( a ) - af (b ) + a f ( a ) b + a 
b - a b - a I = [ f ( b ) - f(a)l + bf(a> - a f ( b ) 
or Multiplying and collecting terms yields 
f1(x) = 
f (b ) - f ( a ) x + bf (a ) - o f @ ) 
b - a b - a 
which can be integrated between x = a and x = b to yield which is the formula for the trapezoidal rule. 
I = f (b) - f (a ) * + 
b - a 2 
21 .1 THE TRAPEZOIDAL RULE 585 
Geometrically, the trapezoidal rule is equivalent to approximating the area of the 
trapezoid under the straight line connecting f(a) and f(b) in Fig. 21.4. Recall from geome- 
try that the formula for computing the area of a trapezoid is the height times the average of 
the bases (Fig. 21.5~). In our case, the concept is the same but the trapezoid is on its side 
(Fig. 21.5b). Therefore, the integral estimate can be represented as 
I E width x average height (21.4) 
FIGURE 2 1.4 
Graphical depiction of the trapezoidal rule. 
FIGURE 2 1.5 
[a) The formula for computin the area of a trapezoid: height times the average of the bases 
[b) For the trapezoidal rule, R t e concept is the same but the trapezoid is on its side. 
586 NEWTON-COTES INTEGRATION FORMULAS 
I E (h - a ) x average height (21.5) 
where, for the trapezoidal rule, the average height is the average of the function values at 
the end points, or [ f (a) + f (b)]/2. 
All the Newton-Cotes closed formulas can be expressed in the general format of 
Eq. (21.5). In fact, they differ only with respect to the formulation of the average height. 
2 1 .I. 1 Error of the Trapezoidal Rule 
When we employ the integral under a straight-line segment to approximate the integral 
under a curve, we obviously can incur an error that may be substantial (Fig. 21.6). An 
estimate for the local truncation error of a single application of the trapezoidal rule is 
(Box. 21.2) 
where < lies somewhere in the interval from a to b. Equation (21.6) indicates that if the 
function being integrated is linear, the trapezoidal rule will be exact. Otherwise, for func- 
tions with second- and higher-order derivatives (that is, with curvature), some error can 
occur. 
FIGURE 2 1.6 
Graphical depiction of the use of a single opplication of the trapezoidal rule to approximate the 
integral of fix) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5 from x = 0 to 0.8. 
21.1 THE TRAPEZOIDAL RULE 587 
An alternative derivation of the trapezoidal rule is possible by inte- constant, this equation can be integrated: 
grating the forward Newton-Gregory interpolating polynomial. Re- 
call that for the first-order version with error term, the integral 1 = h 
would be (Box 18.2) 
a I and evaluated as 1 = Sb [ f ( a ) + Af(o)a + m a ( a - l )h2 dx (821.2.1) 2 
To simplify the analysis, realizethat because a = (x - a)/h , 
Because Af(a) = f(b) - f(a), the result can be written as dx = h da 
Inasmuch as h = b - a (for the one-segment trapezoidal rule), the 
limits of integration a and b correspond to 0 and 1, respectively. 
Therefore, Eq. (B21.2.1) can be expressed as 
I = h ' [ f ( a ) + Af(a)a + m a ( , - l )h2 d o 
2 1 le and the second is an 
approximation for the error. 
If it is assumed that, for small h, the term f"(6) is approximately 
EXAMPLE 2 1 .1 Single Application of the Trapezoidal Rule 
I Problem Statement. Use Eq. (21.3) to numerically integrate 
i 
i f ( x ) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5 
from a = 0 to b = 0.8. Recall from PT6.2 that the exact value of the integral can be deter- 
mined analytically to be 1.640533. 
Solution. The function values 
can be substituted into Eq. (21.3) to yield 
which represents an error of 
which corresponds to a percent relative error of E, = 89.5%. The reason for this large error 
is evident from the graphical depiction in Fig. 21.6. Notice that the area under the straight 
line neglects a significant portion of the integral lying above the line. 
In actual situations, we would have no foreknowledge of the true value. Therefore, 
an approximate error estimate is required. To obtain this estimate, the function's second 
588 NEWTON-COTES INTEGRATION FORMULAS 
1 derivative over the interval can be computed by differentiating the original function twice 
I to give 
i 
fl'(x) = -400 + 4 0 5 0 ~ - 10,800x2 + 8000x3 
I 
! 
, The average value of the second derivative can be computed using Eq. (PT6.4): 
0 8 
I ( -400 + 4050x - 10,800x2 + 8000x3) dx 
i f "(x) = = -60 
0.8 - 0 
i 
a which can be substituted into Eq. (21.6) to yield 
I 
1 which is of the same order of magnitude and sign as the true error. A discrepancy does exist, ' however, because of the fact that for an interval of this size, the average second derivative I . 
is not necessarily an accurate approximation off'((). Thus, we denote that the error is ap- 
proximate by using the notation E,, rather than exact by using E,. 
I 
2 1,1.2 The Multiple-Application Trapezoidal Rule 
One way to improve the accuracy of the trapezoidal rule is to divide the integration interval 
from a to h into a number of segments and apply the method to each segment (Fig. 21.7). 
The areas of individual segments can then be added to yield the integral for the entire 
interval. The resulting equations are called multiple-application, or composite, integration 
formulas. 
Figure 21.8 shows the general format and nomenclature we will use to characterize 
multiple-application integrals. There are n + 1 equally spaced base points (xo, XI , x;?, . . . , 
x,). Consequently, there are n segments of equal width: 
If a and h are designated as xo and x,, respectively, the total integral can be represen- 
ted as 
Substituting the trapezoidal rule for each integral yields 
or, grouping terms, 
n - l 
2 1 .1 THE TRAPEZOIDAL RULE 589 
FlGURE 2 1 -7 
Illustration of the multiple-appl~cation trapezoidal rule. (a) Two segments, (b) three segments, 
(c) four segments, and (d ) five segments. 
590 NEWTON-COTES INTEGRATION FORMULAS 
FIGURE 2 1.8 
The general format and nornen- 
clature for multiple-application 
integrals. 
or, using Eq. (21.7) to express Eq. (21.9) in the general form of Eq. (21.5), 
Because the summation of the coefficients of f (x ) in the numerator divided by 2n is equal 
to 1, the average height represents a weighted average of the function values. According to 
Eq. (21. l o ) , the interior points are given twice the weight of the two end points f(x0) and 
f (-4. 
An error for the multiple-application trapezoidal rule can be obtained by summing the 
individual errors for each segment to give 
2 1 .1 THE TRAPEZOIDAL RULE 5 9 1 
where f1'(Ci) is the second derivative at a point ci located in segment i . This result can be 
simplified by estimating the mean or average value of the second derivative for the entire 
interval as [Eq. (PT6.3)] 
C f " ( t i ) 
Therefore C f "(C;) ;.) n f and Eq. ( 2 1.11) can be rewritten as 
Thus, if the number of segments is doubled, the truncation error will be quartered. Note 
that Eq. (2 1 . 1 3 ) is an approximate error because of the approximate nature of Eq. (2 1.12). 
EXAMPLE 2 1.2 Multiple-Application Trapezoidal Rule 
i Problem Statement. Use the two-segment trapezoidal rule to estimate the integral of 
1 f ( x ) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5 
I 
from a = 0 to h = 0.8. Employ Eq. (21.13) to estimate the error. Recall that the correct 
a value for the integral is 1.640533. 
i Solution. n = 2(h = 0.4): 
: where -60 is the average second derivative determined previously in Example 21.1 
The results of the previous example, along with three- through ten-segment applica- 
tions of the trapezoidal rule, are summarized in Table 21.1. Notice how the error decreases 
as the number of segments increases. However, also notice that the rate of decrease is grad- 
ual. This is because the error is inversely related to the square of n [Eq. (21.13)]. Therefore, 
doubling the number of segments quarters the error. In subsequent sections we develop 
higher-order formulas that are more accurate and that converge more quickly on the true in- 
tegral as the segments are increased. However, before investigating these formulas, we will 
first discuss computer algorithms to implement the trapezoidal rule. 
2 1.1.3 Computer Algorithms For the Trapezoidal Rule 
Two simple algorithms for the trapezoidal rule are listed in Fig. 21.9. The first (Fig. 2 1 . 9 ~ ) 
is for the single-segment version. The second (Fig. 21.9h) is for the multiple-segment 
version with a constant segment width. Note that both are designed for data that is in tabu- 
592 NEWTON-COTES INTEGRATION FORMULAS 
TABLE 21 .I Results for multiple-application trapezoidal 
rule to estimate the integral of f(x) = 0.2 
+ 2 5 x - 200x2 + 675x3 - 900x4 + 400x5 
from x = 0 to 0.8. The exact value is 
1.640533. 
(a) Single-segment (b) Multiple-segment 
FliNCTiON Trap (h, &O, f l ) FUNCTION Trapm (h, n, f ) 
Trap = h ++ j f O + f ? ) / 2 sum = fo 
EWD Trap D O i = I , n - I 
sum = sum + 2 * F, 
END GO 
sum = sum + f , 
Trapm = h * s u m / 2 
END f rapm 
FIGURE 2 1.9 
Algor~thns for the ( a ] sngle-segment and (b] multiple-segment trapezo~dal rule 
lated form. A general program should have the capability to evaluate known functions or 
equations as well. We will illustrate how functions are handled in the next chapter. 
An example of a user-friendly program for the multiple-segment trapezoidal rule is in- 
cluded in the supplementary Numerical Methods TOOLKIT software associated with this 
text. This software can evaluate the integrals of either tabulated data or user-defined func- 
tions. The following example demonstrates its utility for evaluating integrals. It also pro- 
vides a good reference for assessing and testing your own software. 
EXAMPLE 21.3 Evaluating Integrals with the Computer 
Problem Statement. Use the Numerical Methods TOOLKIT software to solve a prob- 
lem related to our friend the falling parachutist. As you recall from Example 1.1, the 
velocity of the parachutist is given as the following function of time: 
21 .1 THE TRAPEZOIDAL RULE 593 
2 where v = velocity (m/s), g = the gravitational constant of 9.8 m/s , m = mass of the 
parachutist equal to 68.1 kg, and c = the drag coefficient of 12.5 kg/s. The model predicts 
the velocity of the parachutist as a function of time as described in Example 1.1. 
Suppose we would like to know how far the parachutist has fallen after a certain time t. 
This distance is given by [Eq. (PT6.5)I 
i 
I where d is the distance in meters. Substituting Eq. (E21.3. I), 
Use the Numerical Methods TOOLKIT software and your own software to determine this I integral with the multiple-segment trapezoidal rule using different numbers of segments. 
1 Notethat performing the integration analytically and substituting known parameter values 
[ results in an exact value of d = 289.43515 m. 
Solution. Click the Integrate Function button on the TOOLKIT main menu to obtain a 
blank screen similar to Fig. 21.10. This screen contains the input and output information 
needed to integrate an analytical function or tabular data. 
First, we can click the Input Function table and enter the function, 
Next click the Input Parameters block and enter values for the lower and upper limits of in- 
tegration of 0 and 10. Next, enter the value 10 for the number of segments along with plot 
dimensions, as in Fig. 21 .lo. 
1 FIGURE21.10 
1 Computer screen showing the integral of a function using the multiple-segment trapezoidal rule 
1 program from the Numerical Methods TOOLKIT. 
j 
594 NEWTON-COTES INTEGRATION FORMULAS 
After clicking the Calc and Plot buttons, a calculated integral of 288.7491 is displayed. 
Thus, we have attained the integral to three significant digits of accuracy. The integral is 
equivalent to the area under the v(t) versus t curve, as shown in Fig. 21.10. Visual inspec- 
tion confirms that the integral is the width of the interval (10 s) times the average height 
(about 29 m/s). Aplot of the function showing the segments is shown in Fig. 21.10. We can 
conveniently try other numbers of segments. By n = 500, the result is accurate to six sig- 
nificant digits. 
At this point, it is important to recognize that the Numerical Methods TOOLKIT uses 
double precision to obtain its integral estimate. We therefore repeated this problem using a 
program based on the pseudocode from Fig. 21.9b and employed single-precision numbers 
(i.e., about seven significant decimal digits of precision). The results are 
Segments Segment size Estimated d, m 
Thus, up to about 500 segments, the multiple-application trapezoidal rule attains 
excellent accuracy. However, notice how the error changes sign and begins to increase in 
absolute value beyond the 500-segment case. The 10,000-segment case actually seems to 
be diverging from the true value. This is due to the intrusion of round-off error because of 
the great number of computations for this many segments. Thus, the level of precision is 
limited, and we would never reach the exact result of 289.4351 obtained analytically. This 
limitation will be discussed in further detail in Chap. 22. 
Three major conclusions can be drawn from the Example 21.3: 
For individual applications with nicely behaved functions, the multiple-segment trape- 
zoidal rule is just fine for attaining the type of accuracy required in many engineering 
applications. 
If high accuracy is required, the multiple-segment trapezoidal rule demands a great deal 
of computational effort. Although this effort may be negligible for a single application, 
it could be very important when (a ) numerous integrals are being evaluated or (b) where 
the function itself is time consuming to evaluate. For such cases, more efficient ap- 
proaches (like those in the remainder of this chapter and the next) may be necessary. 
Finally, round-off error can represent a limitation on our ability to determine integrals. 
This is due both to the machine precision as well as to the numerous computations 
involved in simple techniques like the multiple-segment trapezoidal rule. 
21.2 SIMPSON'S RULES 595 
We now turn to one way in which efficiency is improved. That is, by using higher- 
order polynomials to approximate the integral. 
SIMPSON'S RULES 
Aside from applying the trapezoidal rule with finer segmentation, another way to obtain a 
more accurate estimate of an integral is to use higher-order polynomials to connect the 
points. For example, if there is an extra point midway between f(a) and f(b), the three 
points can be connected with a parabola (Fig. 21.11~). If there are two points equally 
spaced between f(a) and f(b), the four points can be connected with a third-order polyno- 
mial (Fig. 21.1 lb). The formulas that result from taking the integrals under these poly- 
nomials are called Simpson's rules. 
2 1.2.1 Simpson's 1 / 3 Rule 
Simpson's 1/3 rule results when a second-order interpolating polynomial is substituted 
into Eq. (21.1): 
If a and b are designated as xo and x2 and f2(x) is represented by a second-order Lagrange 
polynomial [Eq. (18.23)], the integral becomes 
FIGURE 2 8.1 1 
[a) Graphical depiction of Simp- 
son's 1 /3 rule: it consists of 
takin the area under a 
para E ola connecting three 
points. (b] Graphical depiction 
of Simpson's 318 rule: it con- 
sists of taking the area under a 
cubic equation connecting four 
points. 
596 NEWTON-COTES INTEGRATION FORMULAS 
After integration and algebraic manipulation, the following formula results: 
4(,4 i where, for thls case, h = (b - a ) / 2 . This equation is known as Slmpson's 113 rule. It is the A *.J b 
second Newton-Cotes closed integration formula. The label "1/3" stems from the fact that 
h is divided by 3 in Eq. (21.14). An alternative denvation is shown in Box 21.3 where the 
Newton-Gregory polynom~al is integrated to obtain the same formula. 
Simpson's 113 rule can also be expressed using the format of Eq. (21.5): 
3 Derivation and Error 
As was done in Box 21.2 for the trapezoidal rule, Simpson's 113 
rule can be derived by integrating the forward Newton-Gregory in- 
terpolating polynomial (Box 18.2): 
Notice that we have written the polynomial up to the fourth-order 
term rather than the third-order term as would be expected. The rea- 
son for this will be apparent shortly. Also notice that the limits of 
integration are from xo to x2. Therefore, when the simplifying sub- 
stitutions are made (recall Box 21.2), the integral is from u = 
0 to 2: 
+ - --+ - - - a 2 ) ,f *)(5)h412 ( 2 8 o 
and evaluated for the limits to give 
A2f(x0) 
2f(xo) + 2Af(xo) + --- 
3 
(B21.3.1) 
Notice the significant result that the coefficient of the third divided 
difference is zero. Because Af(xo) = f ( x l ) - f(xo) and A2f(xo) = 
f(x2) - 2f(xl) + f(xo), Eq. (B21.3.1) can be rewritten as 
/ I I 
I = - If(.r") + 4 f ( x t ) + f(sz)j - - f'"'(4)h5 
3 90 
' - 
Strnpson's 1'3 rule Truncation error 
4 u ( u - 1 - 2 - 3 ) da + -- I Thus, the first tern is Simpson's 113 rule and the second is the trun- 24 cation error. Because the third divided difference dropped out, we 
which can be integrated to yield obtain the significant result that the formula is third-order accurate. 
2 1.2 SIMPSON'S RULES 597 
where a = xo, h = x2, and xl = the point midway between a and 6, which is given by 
(b + a) /2 . Notice that, according to Eq. (21.15), the middle point is weighted by two- 
thirds and the two end points by one-sixth. 
It can be shown that a single-segment application of Simpson's 1/3 rule has a trunca- 
tion error of (Box 21.3) 
or, because h = (h - a) /2 , 
where < lies somewhere in the interval from a to b. Thus, Simpson's 113 rule is more ac- 
curate than the trapezoidal rule. However, comparison with Eq. (21.6) indicates that it is 
more accurate than expected. Rather than being proportional to the third derivative, the 
error is proportional to the fourth derivative. This is because, as shown in Box 21.3, the co- 
efficient of the third-order term goes to zero during the integration of the interpolating 
polynomial. Consequently, Simpson's 1 / 3 rule is third-order accurate even though it is 
based on only three points. In other words, it yields exact results for cubic polynomials 
even though it is derived from a parabola! 
EXAMPLE 2 1.4 Single Application of Simpson's 1 /3 Rule 
: Problem Statement. Use ~ q . (21.15) to integrate 
I 
; f ( x ) = 0.2 + 25n - 200x2 + 675x3 - 900x4 + 400x5 
I 
1 from a = 0 to b = 0.8. Recall that the exact integral is 1.640533. 
I Solution. 
i 
I 
Therefore, Eq. ( 2 1.15) can be used to compute 
j which represents an exact error of 
3 
Et = 1.640533 - 1.367467 = 0.2730667 = 16.6% 
1 
i which is approximately5 times more accurate than for a single application of the trape- 
i zoidal rule (Example 21.1). 
I The estimated error is [Eq . (21.16)] 
; where -2400 is the average fourth derivative for the interval as obtained using Eq. (PT6.4). ' As was the case in Example 21.1, the error is approximate (E,) because the average fourth 
598 NEWTON-COTES INTEGRATION FORMULAS 
r derivative is not an exact estimate off (4)(<). However, because this case deals with a fifth- 
: order polynomial, the result matches. 
L--- -- - ----- ------- --- - ---- --------------- - 
2 1.2.2 The Multiple-Application Simpson's 1 /3 Rule 
Just as with the trapezoidal rule, Simpson's rule can be improved by dividing the integra- 
tion interval into a number of segments of equal width (Fig. 21.12): 
h - a 
h = - 
n 
The total integral can be represented as 
Substituting Simpson's 113 rule for the individual integral yields 
+ . . .+ 2h f(xn-2) + 4f(xn-1> + f ( ~ n > 
6 
or, combining terms and using Eq. (21.17), 
1=1,3,5 j=2.4.h 
1 2 (b - u ) (21 18, 
3n -- / 
Width Average height 
FlGURE 21.12 
Graphical representation of the 
multiple application of Simp- 
son's 1 /3 rule Note that the 
method can be employed only 
if the number of segments is 
even. 
21.2 SIMPSON'S RULES 599 
Notice that, as illustrated in Fig. 21.12, an even number of segments must be utilized to 
implement the method. In addition, the coefficients "4" and "2" in Eq. (21.18) might seem 
peculiar at first glance. However, they follow naturally from Simpson's 1/3 rule. The odd 
points represent the middle term for each application and hence carry the weight of 4 from 
Eq. (21.15). The even points are common to adjacent applications and hence are counted 
twice. 
An error estimate for the multiple-application Simpson's rule is obtained in the same 
fashion as for the trapezoidal rule by summing the individual errors for the segments and 
averaging the derivative to yield 
where f (4) is the average fourth derivative for the interval. 
EXAMPLE 2 1.5 Multiple-Application Version of Simpson's 1 / 3 Rule 
I Problem Statement. Use Eq. (21.18) with n = 4 to estimate the integral of i 
1 f(x) = 0.2 + 25.1 - 200x2 + 675x3 - 900x4 + 400x5 
from a = 0 to h = 0.8. Recall that the exact integral is 1.640533. j 
! Solution. n = 4 (h = 0.2): 
I The estimated error [Eq. (21.19)] is 1 
The previous example illustrates that the multiple-application version of Simpson's 
113 rule yields very accurate results. For this reason, it is considered superior to the trape- 
zoidal rule for most applications. However, as mentioned previously, it is limited to cases 
where the values are equispaced. Further, it is limited to situations where there are an even 
number of segments and an odd number of points. Consequently, as discussed in the next 
section, an odd-segment-even-point formula known as Simpson's 3/8 rule is used in con- 
junction with the 1 /3 rule to permit evaluation of both even and odd numbers of segments. 
600 NEWTON-COTES INTEGRATION FORMULAS 
2 1.2.3 Simpson's 318 Rule 
In a similar manner to the derivation of the trapezoidal and Simpson's 113 rule, a third- 
order Lagrange polynomial can be fit to four points and integrated: 
to yield 
where h = (h - a ) / 3 . This equation is called Simpson's 318 rule because h is multiplied by 
318. It is the third Newton-Cotes closed integration formula. The 318 rule can also be ex- 
pressed in the form of Eq. (21.5): 
FIGURE 21.13 
lliustration of how Simpson's 
1 /3 and 3/8 rules can be 
applied in tandem to handle 
rnult~ple applications with odd 
numbers of intervals. 
21.2 SIMPSON'S RULES 60 1 
Thus, the two interior points are given weights of three-eighths, whereas the end points are 
weighted with one-eighth. Simpson's 3/8 rule has an error of 
or, because h = (b - a)/3, 
Because the denominator of Eq. (21.21) is larger than for Eq. (21.16), the 3/8 rule is some- 
what more accurate than the 1 /3 rule. 
Simpson's 113 rule is usually the method of preference because it attains third-order 
accuracy with three points rather than the four points required for the 3/8 version. How- 
ever, the 3/8 rule has utility when the number of segments is odd. For instance, in Ex- 
ample 21.5 we used Simpson's rule to integrate the function for four segments. Suppose 
that you desired an estimate for five segments. One option would be to use a multiple- 
application version of the trapezoidal rule as was done in Examples 21.2 and 21.3. This 
may not be advisable, however, because of the large truncation error associated with this 
method. An alternative would be to apply Simpson's 1/3 rule to the first two segments 
and Simpson's 3/8 rule to the last three (Fig. 21.13). In this way, we could obtain an es- 
timate with third-order accuracy across the entire interval. 
EXAMPLE 2 1.6 Simpson's 3 / 8 Rule 
Problem Statement. 
(a) Use Simpson's 3/8 rule to integrate 
from a = 0 to b = 0.8. 
(b) Use it in conjunction with Simpson's 113 rule to integrate the same function for five 
segments. 
I Solution. 
(a) A single application of Simpson's 3/8 rule requires four equally spaced points: 
,f(O) = 0.2 f(0.2667) = 1.432724 
f(0.5333) = 3.487177 f(0.8) = 0.232 
Using Eq. (21.20), 
0.2 + 3(1.432724 + 3.487177) + 0.232 = 
I 2 0.8 
8 
602 NEWTON-COTES INTEGRATION FORMULAS 
(b) The data needed for a five-segment application (h = 0.16) is 
f(0) = 0.2 f(0.16) = 1.296919 
f(0.32) = 1.743393 f(0.48) = 3.186015 
f(0.64) = 3.181929 f(O.80) = 0.232 
The integral for the first two segments is obtained using Simpson's 1 /3 rule: 
For the last three segments, the 318 rule can be used to obtain 
The total integral is computed by summing the two results: 
I = 0.3803237 + 1.264753 = 1.645077 
E, = 1.640533 - 1.645077 = -0.00454383 E, = -0.28% 
2 1.2.4 Computer Algorithms for Simpson" Rules 
Pseudocode for a number of forms of Simpson's rule are outlined in Fig. 21.14. Note that 
all are designed for data that is in tabulated form. A general program should have the capa- 
bility to evaluate known functions or equations as well. We will illustrate how functions are 
handled in the next chapter. 
Notice that the program in Fig. 21.14d is set up so that either an even or odd number 
of segments may be used. For the even case, Simpson's 113 rule is applied to each pair of 
segments, and the results are summed to compute the total integral. For the odd case, Simp- 
son's 318 rule is applied to the last three segments, and the 113 rule is applied to all the pre- 
vious segments. 
2 1.2.5 Higher-Order Nedon-Cotes Closed Formulas 
As noted previously, the trapezoidal rule and both of Simpson's rules are members of a 
family of integrating equations known as the Newton-Cotes closed integration formulas. 
Some of the formulas are summarized in Table 21.2, along with their truncation-error esti- 
mates. 
Notice that, as was the case with Simpson's 113 and 318 rules, the five- and six-point 
formulas have the same order error. This general characteristic holds for the higher-point 
formulas and leads to the result that the even-segment-odd-point formulas (for example, 
113 rule and Boole's rule) are usually the methods of preference. 
However, it must also be stressed that, in engineering practice, the higher-order (that 
is, greater than four-point) formulas are rarely used. Simpson's rules are sufficient for most 
applications. Accuracy can be improved by using the multiple-application version. Fur- 
thermore, when the function is known and high accuracy is required, methods such as 
21.2 SIMPSON'S RULES 603 
(a) 
FUNCTION Simp13 (h, fO, f?, f 2 ) 
Simp13 = 2*h* (f0+4*Fl+f2) / 6 
END Simp13 
(b) 
FUNCTION Simp38 (h, fO, f?, f2, f 3 ) 
S imp38 = 3*h* ( f0+3* ( f l t f 2 )+ f3 ) / 8 
END 5 imp38 
( ~ 9 
FUNCTION 5imp13m [h, n, f ) 
sum = f ( 0 ) 
D O ; = 1,n - 2 , 2 
s u m = s u m + 4 * f , + 2 * F , + 7 
END DO 
sum = sum + 4 * f n - , + f, 
Sirnp13m = h * sum / 3 
END SimpI3m(4 
FUNCTION Simplnt(a, b, n, f ) 
h = ( b - a ) / n 
IF n = 1 THEN 
sum = Trap(h, f,-l, f,) 
ELSE 
m = n 
odd = n / 2 - INT(n / 2 ) 
I F o d d > O A N D n > I T H E N 
sum = sum+5imp38(h,fn-sfn-2,fn-1,fn) 
m = n - 3 
END IF 
IF m > 1 THEN 
sum = sum + 5imp13m(h, m, f ) 
END IF 
END IF 
Simplnt = sum 
END Simplnt 
FBOURE 21 . I 4 
Pseudocode for Simpson's rules. (a ) Single-application Simpson's 1 / 3 rule, [b) s~ngle-application 
Simpson's 318 rule, [c) multiple-application Simpson's 1 / 3 rule, and [d) multiple-application 
Simpson's rule for both odd and even number of segments. Note that for all cases n must be > I 
TABLE 21 Newton-Cotes closed integration formulas. The formulas are presented in the 
format of Eq. (21.5) so that the weighting of the data points to estimate the 
average height is apparent. The step size i s given by h = (b - o)/n. 
Segmenb 
(4 Points Name Formula Truncation Error 
1 2 f(xoi + {(xi i Trapezoidal rule (b - a) 7 
2 3 Simpson's 1 / 3 rule (b - a] fix01 + 4fix1) + fix21 
6 
3 4 Simpson's 318 rule (b - a) fixo) + 3fIx11 + 3 f ( ~ ~ ] + Fix3] 
8 - i3/80)h5f141([) 
4 5 Boole's rule Ib - a1 7fIxol + 32fix1) + 12fIxz) + 32fix31 + 7f(x4) 
90 - (8/945)h7f161[c] 
604 NEWTON-COTES INTEGRATION FORMULAS 
Romberg integration or Gauss quadrature, described in Chap. 22, offer viable and attrac- 
tive alternatives. 
2 1.3 IlUTEGlZATlON WITH UNEQUAL SEGMENTS 
To this point, all formulas for numerical integration have been based on equally spaced 
data points. In practice, there are many situations where this assumption does not hold and 
we must deal with unequal-sized segments. For example, experimentally derived data is 
often of this type. For these cases, one method is to apply the trapezoidal rule to each seg- 
ment and sum the results: 
where hi = the width of segment i. Note that this was the same approach used for the 
multiple-application trapezoidal rule. The only difference between Eqs. (21.8) and (21.22) 
is that the h's in the former are constant. Consequently, Eq. (21.8) could be simplified 
by grouping terms to yield Eq. (21.9). Although this simplification cannot be applied to 
Eq. (21.22), a computer program can be easily developed to accommodate unequal-sized 
segments. Before describing such an algorithm, we will illustrate in the following exam- 
ple how Eq. (21.22) is applied to evaluate an integral. 
EXAMPLE 2 1.7 Trapezoidal Rule with Unequal Segments 
Problem Statement. The information in Table 21.3 was generated using the same 
polynomial employed in Example 21.1. Use Eq. (21,22) to determine the integral for this 
data. Recall that the correct answer is 1.640533. 
Solution. Applying Eq. (21.22) to the data in Table 21.3 yields 
which represents an absolute percent relative error of E, = 2.8% 
TABLE 21.3 Data for f(x) = 0.2 + 25x - 200x2 + 675x3 - 
900x4 + 400x5, with unequally spaced values 
X f(X) X f (XI 
0 0 0 200000 0.44 2.842985 
0.12 1.309729 0.54 3.507297 
0.22 1.305241 0.64 3.181929 
0.32 1 743393 0.70 2 363000 
0.36 2.074903 0.80 0.232000 
0 40 2.456000 
21.3 INTEGRATION WITH UNEQUAL SEGMENTS 605 
FIGURE 2 1.1 5 
Use of the trapezoidal rule to determine the integral of unevenly spaced data. Notice how the 
shaded segments could be evaluated with Sirnpson's rule to attain higher accuracy. 
The data from Example 21.7 is depicted in Fig. 21.15. Notice that some adjacent seg- 
ments are of equal width and, consequently, could have been evaluated using Simpson's 
rules. This usually leads to more accurate results, as illustrated by the following example. 
EXAMPLE 2 1.8 Inclusion of Simpson's Rules in the Evaluation of Uneven Data 
1 Problem Statement. Recompute the integral for the data in Table 21.3, but use 
/ Simpson's rules for those segments where they are appropriate. 
Solution. The first segment is evaluated with the trapezoidal rule: 1 
1 Because the next two segments from x = 0.12 to 0.32 are of equal length, their integral can 
be computed with Simpson's 113 rule: 
8 
1 The next three segments are also equal and, as such, may be evaluated with the 3/8 rule to 
! give I = 0.2726863. Similarly, the 1/3 rule can be applied to the two segments from 
x = 0.44 to 0.64 to yield I = 0.6684701. Finally, the last two segments, which are of un- 
I equal length, can be evaluated with the trapezoidal rule to give values of 0.1663479 and 
1 0.1297500, respectively. The area of these individual segments can be summed to yield a 
606 NEWTON-COTES INTEGRATION FORMULAS 
1 total integral of 1.603641. This represents an error of el = 2.2%, which is superior to the 
! result using the trapezoidal rule in Example 21.7. 
L 
Computer Program for Unequally Spaced Data. It is a fairly simple proposition to 
program Eq. (21.22). Such an algorithm is listed in Fig. 21.16~. 
However, as demonstrated in Example 21.8, the approach is enhanced if it implements 
Simpson's rules wherever possible. For this reason, we have developed a second algorithm 
that incorporates this capability. As depicted in Fig 21.16b, the algorithm checks the length 
of adjacent segments. If two consecutive segments are of equal length, Simpson's 113 rule 
is applied. If three are equal, the 318 rule is used. When adjacent segments are of unequal 
length, the trapezoidal rule is implemented. 
FUGURE 2 1.16 
Pseudocode for integratin unequally spaced data. (a) Trapezoidal rule and (b) combination 
Simpson's and trapezoida 9 rules. 
(4 
FUNCTION Trapun (x, y, n) 
LOCAL i, sum 
5um = 0 
D O i = l , n 
sum = sum + (xi - xi-,)~(yi-, + yi)/2 
END DO 
Trapun = sum 
END Trapun 
(6) 
FUNCTION Uneven (n, x, f ) 
h = XI - xo 
k = 1 
sum = 0. 
D O j = 1, n 
h f = xj+, - X, J 
IF A 0 5 ( h - h f ) < .000001 THEN 
IF k = 3 THEN 
sum = sum + Simp13 (h, 5-3,6-8G-I) 
k = k - 1 
ELSE 
k = k + 1 
END IF 
ELSE 
IF k = 1 THEN 
sum = sum + Trap (h, j-,,$) 
ELSE 
IF k = 2 THEN 
sum = sum + Simp13 (h, cT2, 5-,, 5) 
ELSE 
sum = sum + S i m p 3 8 (h, 5-3, 5-2, 5-,, 5) 
END IF 
k = l 
END IF 
END IF 
h = h f 
END DO 
Uneven = sum 
END Uneven 
PROBLEMS 607 
TABLE 21.4 Newton-Cotes open integration formulas. The formulas are presented in the for- 
mat of Eq. (2 1.5) so that the weighting of the data points to estimate the aver- 
age height i s apparent. The step size i s given by h = (b - a) /n . 
Segments 
(4 Points Name Formula Truncation Error 
2 1 Midpoint method lb - oi fix11 [ 1 /31h3plt) 
.S 4 
1 1 flx,! + flx2i + fjxgi + 1 1 fjx',) 
I b - ai 24 195,') 44)h5f[4lje) 
Thus, not only does it allow evaluation of unequal segment data, but if equally 
spaced information is used, it reduces to Simpson's rules. As such, it represents a basic, 
all-purpose algorithm for the determination of the integral of tabulated data. 
21.4 OPEN INTEGMTION FORMULAS 
Recall from Fig 21.3b that open integration formulas have limits that extend beyond the 
range of the data. Table 21.4 summarizes the Newton-Cotes open integration formulas. 
The formulas are expressed in the form of Eq. (21.5) so that the weighting factors are 
evident. As with the closed versions, successive pairs of the formulas have the same order 
error. The even-segment-odd-point formulas are usually the methods of preference be- 
cause they require fewer points to attain the same accuracy as the odd-segment-even-point 
formulas. 
The open formulas are not often used for definite integration. However, as discussed 
in Chap. 22, they have utility for analyzing improper integrals. In addition, they will have 
relevance to our discussion of multistep methods for solving ordinary differential equa- 
tions in Chap. 26. 
PROBLEMS 
21.1 Use analytical means to evaluate 
(a) (I - .-"I d* 
0 
21.2 Use a single application of the trapezoidal rule to evaluate the 
integrals from Prob. 21.1. 
21.3 Evaluate the integrals from Prob. 21.1 with a multiple- 
application trapezoidal rule, with n = 2,4, and 6. 
21.4 Evaluate the integrals from Prob. 21.1 with a singleapplica- 
tion of Simpson's 113 rule. 
21.5 Evaluate the integrals from Prob. 21.1 with a multiple- 
application Simpson's 113 rule, with n = 4 and 6. 
608 NEWTON-COTES INTEGRATION FORMULAS 
21.6 Evaluate the integrals from Prob. 21.1 with a single applica- 
tion of Simpson's 318 rule. 
21.7 Evaluate the integrals from Prob. 21.1, but use a multiple- 
application Simpson's rule, with n = 5. 
21.8 Integrate the following function using the trapezoidal rule, 
with n = 1, 2, 3, and 4: 
Compute percent relative errors with respect to the true value of 
4.8333 to evaluate the accuracy of the trapezoidal approximations. 
21.9 Integrate the following function both analytically and using 
Simpson's rules, with n = 4 and 5: 
Discuss the results 
21.10 Integrate the following function both analytically and nu- 
merically. Use both the trapezoidal and Simpson's 113 rules to nu- 
merically integrate the function. For both cases, use the multiple- 
application version, with n = 4. 
Compute percent relative errors for the numerical results. 
21.11 Integrate the following function both analytically and nu- 
merically. For the numerical evaluations use (a) a single applica- 
tion of the trapezoidal rule, (b) Simpson's 113 rule, (c) Simpson's 
318 rule, (d) Boole's rule, (e) the midpoint method, ( f ) the 3 
segmentl2-point open integration formula, and (g) the 4- 
segmentl3-point open integration formula. 
21.13 Integrate the following function, 
Note that the true value is 0.602297. Evaluate this integral with the 
multiple-segment trapezoidal rule. Use a sufficiently high n that 
you get 4 significant digits of accuracy. Discuss your results. 
21.14 Evaluate the integral of the following tabular data with the 
trapezoidal rule: 
21.15 Perform the same evaluation as in Prob. 21.14, but use 
Simpson's rules. 
21.16 Evaluate the integral of the following tabular data using the 
trapezoidal rule: 
x 1 -3 - 1 1 3 5 7 9 1 1 
f ix) 1 1 -4 -9 2 4 2 6 -3 
21.17 Perform the same evaluation as in Prob. 21.16, but use 
Simpson's rules. 
21.18 Determine the mean value of the function 
between x = 2 and 10 by (a) graphing the function and visually 
estimating the mean value, (b) using Eq. (PT6.4) and the analyti- 
cal evaluation of the integral, and (c) using Eq. (PT6.4) and a five- 
segment version of Simpson's rule to estimate the integral. Calcu- 
late the relative percent error. 
21.19 The function 
f(x) = e-' 
can be used to generate the following table of unequally spaced data: 
Compute percent relative errors for the numerical results. x I 0 0.1 0.3 0.5 0.7 0.95 1.2 
21.12 Integrate the following function both analytically and fix) I 1 0.9048 07408 06065 0 4966 0.3867 0.301 2 
numerically. For the numerical evaluations use (a) single applica- 
Evaluate the integral from a = 0 to b = 1.2 using (a) analytical 
tion of the trapezoidal rule: (b) Simpson's 113 rule, (c) Simpson's 
means, (b) the trapezoidal rule, and (c) a of the trape- 
318 rule, (d) multiple application of Simpson's rules (n = 5 ) , 
zoidal and Simpson's rules; employ Simpson's rules wherever pos- 
(e) Boole's rule, ( f ) the midpoint method, (g) the 3-segment/2- 
sible to obtain the highest accuracy. For (b) and (c), compute the 
point open integration formula, and (h) the 4-segmentl3-point 
percent relative error ( E ~ ) . 
open integration formula. 
21.20 Evaluate the following double integral: 
P R 
J (5 + 3 sin x) dx 
Compute percent relative errors for the numerical results. 
PROBLEMS 609 
(a) analytically. (b) using a multiple-application trapezoidal rule data. Test your program by duplicating the computation from Ex- 
(n = 2), (c) using single applications of Simpson's 113 rule. For ample 21.2. 
(b) and (c), compute the percent relative error ( E ~ ) . 21.23 Develop a user-friendly computer program for the multiple- 
21.21 Evaluate the triple integral application version of Simpson's rule based on Fig. 21.14~. Test it 
by duplicating the computations from Example 21.5. l6 /: (x3 - 2y2) d . ~ d y dr 21.24 Develop a user-friendly computer program for integrating 
unequally spaced data based on Fig. 21.16b. Test it by duplicating 
the computation from Example 2 1.8. 
(a) analytically and (b) using single applications of Simpson's 113 21-25 use the Integrate Function program on the Numerical 
rule. For (b) compute the percent relative error (s,). Methods TOOLKIT disk (or your own program from Prob. 2 1.22) to 
Develop a user-friendly computer program for the repeat (a) Prob. 21.2, (b) Prob. 21.3, (c) Prob. 21.8, (d) Prob. 21.10, 
application trapezoidal rule based on Fig. 21.9. Among (e) Prob. 21.12. Use the graphical option to help you visualize 
things, (a) add documentation statements to the code, (b) make the the concept that I = lab f ( x ) d x is the area between the f(.x) curve 
and Output user-Oriented, and (c) the program so that and the axis, Try several different step sizes for each problem, 
it is capable of evaluating given functions in addition to tabular 
CHAPTER 22 
ntegration o Equations 
In the introduction to Part Six, we noted that functions to be integrated numerically will 
typically be of two forms: a table of values or a function. The form of the data has an im- 
portant influence on the approaches that can be used to evaluate the integral. For tabulated 
information, you are limited by the number of points that are given. In contrast, if the func- 
tion is available, you can generate as many values of f ( x ) as are required to attain accept- 
able accuracy (recall Fig. PT6.7). 
This chapter is devoted to two techniques that are expressly designed to analyze cases 
where the function is given. Both capitalize on the ability to generate function values to de- 
velop efficient schemes for numerical integration. The first is based on Richardson's ex- 
trapolation, which is a method for combining two numerical integral estimates to obtain a 
third, more accurate value. The computational algorithm for implementing Richardson's 
extrapolation in a highly efficient manner is called Romherg integration. This technique 
is recursive and can be used to generate an integral estimate within a prespecified error 
tolerance. 
The second method is called Gauss quadrature. Recall that, in the last chapter, values 
of f (x ) for the Newton-Cotes formulas were determined at specified values of x. For exam- 
ple, if we used the trapezoidal rule to determine an integral, we were constrained to take the 
weighted average of f ( x ) at the ends of the interval. Gauss-quadrature formulas employ x 
values that are positioned between a and h in such a manner that a much more accurate in- 
tegral estimate results. 
In addition to these two standard techniques, we devote a final section to the evalua- 
tion of improper integrals. In this discussion, we focus on integrals with infinite limits and 
show how a change of variable and open integration formulas prove useful for such cases. 
22.1 NEWTON-COTES ALGORITHMS FOR EQUATlONS 
In Chap. 21, we presented algorithms for multiple-application versions of the trapezoidal 
rule and Simpson's rules. Although these pseudocodes can certainly be used to analyze 
equations, in our effort to make them compatible with either data or functions, they could 
not exploit the convenience of the latter. 
Figure 22.1 shows pseudocodes that are expressly designed for cases where the func- 
tion is analytical. In particular, notice that neither the independent nor the dependent 
FlGURE 22.9 
Algorithms for multiple appli- 
cations of the (a) trapezoidal 
and (b) Simpson's 113 rules, 
where the function is available. 
FIGURE 22.2 
Absolute value of the true per- 
cent relative error versus number 
of segments for the determina- 
tion of the integral of f ix) = 
0.2 + 25x - 200x2 + 
675x3 - 900x4 + 400x5, 
evaluated from a = O to 
b = 0.8 using the multiple- 
application trapezoidal rule and 
the multiple-application Sirnp- 
son's 1 /3 rule. Note that bothresults indicate that for a large 
number of segments, round-off 
errors limit precision. 
22.1 NEWTON-COTES ALGORITHMS FOR EQUATIONS 61 1 
(a) 
FUNCTION TrapEq (n, a, b) 
h = ( b - a ) / n 
x = a 
sum = f(x) 
D O i = 1 , n - 1 
x = x + h 
sum = sum + 2 a ( x ) 
END DO 
sum = sum + f (b ) 
TrapEq = ( b - a) a sum / (2 * n) 
END TrapEq 
(61 
FUNCnON 5impEq (n, a, b) 
h = ( b - a ) / n 
x = a 
sum = (x) 
D O i = ? , n - 2,2 
x = x + h 
sum = sum + 4 a(x) 
x = x + h 
sum = sum + 2 x f ( x ) 
END DO 
x = x + h 
sum = sum + 4 ( x ) 
sum = sum + f(b) 
5impEq = (b - a ) * sum / (3 n) 
END 5impEq 
61 2 INTEGRATION OF EQUATIONS 
variable values are passed into the function via its argument as was the case for the codes 
in Chap. 21. For the independent variable x, the integration interval (a , b) and the number 
of segments are passed. This information is then employed to generate equispaced values 
of x within the function. For the dependent variable, the function values in Fig. 22.1 are 
computed using calls to the function being analyzed, f(x). 
We developed single-precision programs based on these pseudocodes to analyze the 
effort involved and the errors incurred as we progressively used more segments to esti- 
mate the integral of a simple function. For an analytical function, the error equations 
[Eqs. (21.13) and (21.19)] indicate that increasing the number of segments n will result in 
more accurate integral estimates. This observation is borne out by Fig. 22.2, which is a plot 
of true error versus n for the integral of f(x) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 
400x5. Notice how the error drops as n increases. However, also notice that at large values 
of n, the error starts to increase as round-off errors begin to dominate. Also observe that a 
very large number of function evaluations (and, hence, computational effort) is required to 
attain high levels of accuracy. As a consequence of these shortcomings, the multiple- 
application trapezoidal rule and Simpson's rules are sometimes inadequate for problem 
contexts where high efficiency and low errors are needed. 
Romberg integration is one technique that is designed to attain efficient numerical integrals 
of functions. It is quite similar to the techniques discussed in Chap. 21, in the sense that it 
is based on successive application of the trapezoidal rule. However, through mathematical 
manipulations, superior results are attained for less effort. 
22.2.1 Richardson's Extrapolation 
Recall that, in Sec. 10.3.3, we used iterative refinement to improve the solution of a set of 
simultaneous linear equations. Error-correction techniques are also available to improve 
the results of numerical integration on the basis of the integral estimates themselves. Gen- 
erally called Richardson's extrapolation, these methods use two estimates of an integral to 
compute a third. more accurate approximation. 
The estimate and error associated with a multiple-application trapezoidal rule can be 
represented generally as 
where I = the exact value of the integral, I(h) = the approximation from an n-segment ap- 
plication of the trapezoidal rule with step size h = (b - a)/n, and E(h) = the truncation 
error. If we make two separate estimates using step sizes of hl and h2 and have exact val- 
ues for the error, 
Now recall that the error of the multiple-application trapezoidal rule can be represented ap- 
proximately by Eq. (21.13) [with rz = (b - a)/h]: 
h - a 
? h 2 7 (22.2) 
22.2 ROMBERG INTEGRATION 61 3 
If it is assumed that j" is constant regardless of step size, Eq. (22.2) can be used to deter- 
mine that the ratio of the two errors will be 
This calculation has the important effect of removing the term f from the computation. In 
so doing, we have made it possible to utilize the information embodied by Eq. (22.2) with- 
out prior knowledge of the function's second derivative. To do this, we rearrange Eq. (22.3) 
to give 
which can be substituted into Eq. (22.1): 
which can be solved for 
Thus, we have developed an estimate of the truncation error in terms of the integral esti- 
mates and their step sizes. This estimate can then be substituted into 
I = l ( h z ) + E(h2) 
to yield an improved estimate of the integral: 
It can be shown (Ralston and Rabinowitz, 1978) that the error of this estimate is 0(h4) . 
Thus, we have combined two trapezoidal rule estimates of 0 (h2 ) to yield a new estimate 
of O(h4). For the special case where the interval is halved (h2 = h 1 / 2 ) , this equation 
becomes 
or, collecting terms, 
EXAMPLE 22.1 Error Corrections of the Trapezoidal Rule 
Problem Statement. In the previous chapter (Example 21.1 and Table 21.1), we used a 
variety of numerical integration methods to evaluate the integral of f(x-) = 0.2 + 251- - 
200x2 + 675x3 - 900x4 + 400x5 from a = 0 to h = 0.8. For example, single and 
614 INTEGRATION OF EQUATIONS 
/ multiple applications of the trapezoidal rule yielded the following results: 
Segments h Integral E,, % 
Use this information along with Eq. (22.5) to compute improved estimates of the integral. 
Solution. The estimates for one and two segments can be combined to yield 
The error of the improved integral is E, = 1.640533 - 1.367467 = 0.273067 ( 6 , = 
16.6%), which is superior to the estimates upon which it was based. 
In the same manner, the estimates for two and four segments can be combined to give 
Equation (22.4) provides a way to combine two applications of the trapezoidal rule 
with error 0(h2) to compute a third estimate with error 0(h4). This approach is a subset of 
a more general method for combining integrals to obtain improved estimates. For in- 
stance, in Example 22.1, we computed two improved integrals of 0(h4) on the basis of 
three trapezoidal rule estimates. These two improved estimates can, in turn, be combined 
to yield an even better value with 0(h6). For the special case where the original trape- 
zoidal estimates are based on successive halving of the step size, the equation used for 
0(h6) accuracy is 
where I , and I I are the more and less accurate estimates, respectively. Similarly, two 0(h6) 
results can be combined to compute an integral that is O(hs) using 
EXAMPLE 22.2 Higher-Order Error Correction of Integral Estimates 
Problem Statement. In Example 22.1, we used Richardson's extrapolation to compute 
two integral estimates of 0(h4) . Utilize Eq. (22.6) to combine these estimates to compute 
an integral with 0(h6) . 
22.2 ROMBERG INTEGRATION 61 5 
Solution. The two integral estimates of O(h4) obtained in Example 22.1 were 1.367467 
and 1.623467. These values can be substituted into Eq. (22.6) to yield 
I which is the correct answer to the seven significant figures that are carried in this example. 
22.2.2 The Romberg integration Algorithm 
Notice that the coefficients in each of the extrapolation equations [Eqs. (22.5), (22.6), and 
(22.7)] add up to 1. Thus, they represent weighting factors that, as accuracy increases, 
place relatively greater weight on the superior integral estimate. These formulations can be 
expressed in a general form that is well-suited for computer implementation: 
where I,+,, k-1 and I,, k-1 = the more and less accurate integrals, respectively, and I,,J = the 
improved integral. The index k signifies the level of the integration, where k = I corre- 
sponds to the original trapezoidal rule estimates, k = 2 corresponds to 0(h4)), k = 3 to 
O(hb), and so forth. The index j is used to distinguish between the more ( B + 1) and the less 
( j ) accurate estimates. For example, for k = 2 and j = 1, Eq. (22.8) becomes 
which is equivalent to Eq. (22.5). 
The general form represented by Eq. (22.8) is attributed to Romberg, and its system- 
atic application to evaluate integrals is known as Romberg inregration. Figure 22.3 is a 
graphical depiction of the sequence of integral estimates generated using this approach. 
Each matrix corresponds to a singleiteration. The first column contains the trapezoidal rule 
FIGURE 22.3 
Graphical depiction of the 
O(h2) O(h4) 
sequence of integral estimates (4 0.172800 1 367467 
generated using Rornberg 1 068800 
~ntegration. 
616 INTEGRATION OF EQUATIONS 
evaluations that are designated I j ,~ , where j = 1 is for a single-segment application (step 
size is b - a), j = 2 is for a two-segment application [step size is (b - a)/2], j = 3 is for 
a four-segment application [step size is (b - a)/4], and so forth. The other columns of the 
matrix are generated by systematically applying Eq. (22.8) to obtain successively better es- 
timates of the integral. 
For example, the first iteration (Fig. 2 2 . 3 ~ ) involves computing the one- and two- 
segment trapezoidal rule estimates and 12,1). Equation (22.8) is then used to compute 
the element = 1.367467, which has an error of 0(h4). 
Now, we must check to determine whether this result is adequate for our needs. As in 
other approximate methods in this book, a termination, or stopping, criterion is required to 
assess the accuracy of the results. One method that can be employed for the present pur- 
poses is [Eq. (3.5)] 
where E, = an estimate of the percent relative error. Thus, as was done previously in other 
iterative processes, we compare the new estimate with a previous value. When the change 
between the old and new values as represented by E, is below a prespecified error criterion 
E,, the computation is terminated. For Fig. 22.3a, this evaluation indicates an 87.4 percent 
change over the course of the first iteration. 
The object of the second iteration (Fig. 22.3b) is to obtain the 0(h6) estimate-Il,s. To 
do this, an additional trapezoidal rule estimate, 13,~ = 1.4848, is determined. Then it is com- 
bined with 12,1 using Eq. (22.8) to generate 12,~ = 1.623467. The result is, in turn, combined 
FIGURE 22.4 
Pseudocode for Romberg integration that uses the equal-size-segment version of the trapezoidal 
rule from Fig. 22. 1 . 
FUNCnON Rhomberg (a, b, maxit, es) 
LOCAL 1(10,10) 
n = 1 
/?,I = TrapEq(n, a, b) 
i te r = 0 
DO 
i ter = i te r + 1 , = 2 iter 
/iter+l,l = TrapEq(n, a, b) 
DO k = 2, i ter + 1 
j = 2 + i te r - k 
1. ~ , k - - (4k-I * $+,,k-l - $,k-l) / (4k-I - 1) 
END DO 
ea = ABS((ll,itcr+l - 11,iter) / Il,iter+l) 100 
i f ( i t e r 2 maxit OR ea S es) EXIT 
END DO 
Rhomberg = /l,iter+l 
END Rhomberg 
22.3 GAUSS QUADRATURE 617 
with 11,2 to yield I l s = 1.640533. Equation (22.9) can be applied to determine that this 
result represents a change of 22.6 percent when compared with the previous result 1 1 ~ 2 . 
The third iteration (Fig. 22 .3~ ) continues the process in the same fashion. In this case, 
a trapezoidal estimate is added to the first column, and then Eq. (22.8) is applied to com- 
pute successively more accurate integrals along the lower diagonal. After only three itera- 
tions, because we are evaluating a fifth-order polynomial, the result (11,4 = 1.640533) is 
exact. 
Romberg integration is more efficient than the trapezoidal rule and Simpson's rules 
discussed in Chap. 21. For example, for determination of the integral as shown in Fig. 22.1, 
Simpson's 1/3 rule would require a 256-segment application to yield an estimate of 
1.640533. Finer approximations would not be possible because of round-off error. In con- 
trast, Romberg integration yields an exact result (to seven significant figures) based on 
combining one-, two-, four-, and eight-segment trapezoidal rules; that is, with only 15 
function evaluations! 
Figure 22.4 presents pseudocode for Romberg integration. By using loops, this algo- 
rithm impleme~lts the method in an efficient manner. Romberg integration is designed for 
cases where the function to be integrated is known. This is because knowledge of the func- 
tion permits the evaluations required for the initial implementations of the trapezoidal rule. 
Tabulated data is rarely in the form needed to make the necessary successive halvings. 
GAUSS QUADMTURE 
In Chap. 21, we studied the group of numerical integration or quadrature formulas known 
as the Newton-Cotes equations. A characteristic of these formulas (with the exception of 
the special case of Sec. 21.3) was that the integral estimate was based on evenly spaced 
function values. Consequently, the location of the base points used in these equations was 
predetermined or fixed. 
For example, as depicted in Fig. 22.5a, the trapezoidal rule is based on taking the area 
under the straight line connecting the function values at the ends of the integration interval. 
The formula that is used to compute this area is 
where a and b = the limits of integration and b - a = the width of the integration inter- 
val. Because the trapezoidal rule must pass through the end points, there are cases such as 
Fig. 22.5a where the formula results in a large error. 
Now, suppose that the constraint of fixed base points was removed and we were free 
to evaluate the area under a straight line joining any two points on the curve. By position- 
ing these points wisely, we could define a straight line that would balance the positive and 
negative errors. Hence, as in Fig. 22.56, we would arrive at an improved estimate of the 
integral. 
Gauss quadrature is the name for one class of techniques to implement such a strat- 
egy. The particular Gauss quadrature formulas described in this section are called Gauss- 
Legendl-e formulas. Before describing the approach, we will show how numerical in- 
tegration formulas such as the trapezoidal rule can be derived using the method of 
INTEGRATION OF EQUATIONS 
FlGURE 22.5 
(a) Graphical depiction of the 
trapezoidal rule as the area 
under the straight line ioining 
fixed end points. (b) An 
improved integral estimate ob- 
tained by taking the area under 
the straight line passing through 
two intermediate points. By po- 
sitioning these points wisely, the 
ositive and negative errors are 
Lalance, and an improved in- 
tegral estimate results. 
undetermined coefficients. This method will then be employed to develop the Gauss- 
Legendre formulas. 
22.3.1 Method of Undetermined Coefficients 
In Chap. 21, we derived the trapezoidal rule by integrating a linear interpolating polynomial 
and by geometrical reasoning. The method of undetermined coefficients offers a third ap- 
proach that also has utility in deriving other integration techniques such as Gauss quadrature. 
To illustrate the approach, Eq. (22.10) is expressed as 
where the c's = constants. Now realize that the trapezoidal rule should yield exact results 
when the function being integrated is a constant or a straight line. Two simple equations 
that represent these cases are y = 1 and y = x. Both are illustrated in Fig. 22.6. Thus, the 
following equalities should hold: 
C O + C ~ = 1 dx 
and 
22.3 GAUSS QUADRATURE 61 9 
FlGURE 22.6 
Two integrals that should be evaluated exactly by the trapezoidal rule: (a) a constant and (b) a 
straight line. 
or, evaluating the integrals, 
and 
b - a b - a 
-Co - +cl2=O 
2 
These are two equations with two unknowns that can be solved for 
b - a 
Co = C1 = -- 
2 
which, when substituted back into Eq. (22.1 I ) , gives 
b - a b - a 
I = - f (a) + T f ( b ) 
2 
which is equivalent to the trapezoidal rule. 
620 INTEGRATION OF EQUATIONS 
22,3.2 Derivation of the Tvvo-Point Gauss-Legendre Formula 
Just as was the case for the above derivation of the trapezoidal rule, the object of Gauss 
quadrature is to determine the coefficients of an equation of the fonn 
where the c's = the unknown coefficients. However, in contrast to the trapezoidal rule that 
used fixed end points a and b, the function arguments xo and xl are not fixed at the end 
points, but are unknowns (Fig. 22.7). Thus, we now have a total of four unknowns that 
must be evaluated, and consequently, we require four conditions to determine them 
exactly. 
Just as for the trapezoidal rule, we can obtaintwo of these conditions by assuming that 
Eq. (22.12) fits the integral of a constant and a linear function exactly. Then, to arrive at the 
other two conditions, we merely extend this reasoning by assuming that it also fits the in- 
tegral of a parabolic (y = 2) and a cubic (y = x3) function. By doing this, we determine 
all four unknowns and in the bargain derive a linear two-point integration formula that is 
exact for cubics. The four equations to be solved are 
co f(xo) + cl f(x1) = (22.13) 
co f(xo) f cl f(x1) = (22.14) 
FIGURE 22.7 
Graphical depiction of the unknown variables xo and X I for integration by Gauss quadrature 
22.3 GAUSS QUADRATURE 62 1 
Equations (22.13) through (22.16) can be solved simultaneously for 
which can be substituted into Eq. (22.12) to yield the two-point Gauss-Legendre formula 
Thus, we arrive at the interesting result that the simple addition of the function values at 
x = 1/& and -I/& yields an integral estimate that is third-order accurate. 
Notice that the integration limits in Eqs. (22.13) through (22.16) are from - 1 to 1. This 
was done to simplify the mathematics and to make the formulation as general as possible. 
A simple change of variable can be used to translate other limits of integration into this 
form. This is accomplished by assuming that a new variable xd is related to the original 
variable x in a linear fashion, as in 
X = a0 $. alxd (22.18) 
If the lower limit, x = a, corresponds to xd = -1, these values can be substituted into 
Eq. (22.18) to yield 
a = ao +al(-1) (22.19) 
Similarly, the upper limit, x = b, corresponds to xd = 1, to give 
b = ag + a l ( l ) (22.20) 
Equations (22.19) and (22.20) can be solved simultaneously for 
and 
b - a 
a1 = - 
2 
which can be substituted into Eq. (22.18) to yield 
This equation can be differentiated to give 
h - a 
dx = - 
2 dxd 
Equations (22.23) and (22.24) can be substituted for x and dx, respectively, in the equation 
to be integrated. These substitutions effectively transform the integration interval without 
22.3 GAUSS QUADRATURE 62 1 
Equations (22.13) through (22.16) can be solved simultaneously for 
which can be substituted into Eq. (22.12) to yield the two-point Gauss-Legendre formula 
Thus, we arrive at the interesting result that the simple addition of the function values at 
x = I /& and - 1 /& yields an integral estimate that is third-order accurate. 
Notice that the integration limits in Eqs. (22.13) through (22.16) are from - 1 to 1. This 
was done to simplify the mathematics and to make the formulation as general as possible. 
A simple change of variable can be used to translate other limits of integration into this 
form. This is accomplished by assuming that a new variable xd is related to the original 
variable x in a linear fashion, as in 
x = a0 + alxd (22.18) 
If the lower limit, x = a, corresponds to xd = -1, these values can be substituted into 
Eq. (22.18) to yield 
a = a o +a l ( -1 ) (22.19) 
Similarly, the upper limit. x = b, corresponds to xd = 1, to give 
b = a o f a ~ ( l ) (22.20) 
Equations (22.19) and (22.20) can be solved simultaneously for 
and 
b - a 
a ] = - 
2 
which can be substituted into Eq. (22.18) to yield 
This equation can be differentiated to give 
Equations (22.23) and (22.24) can be substituted for s and ds, respectively, in the equation 
to be integrated. These substitutions effectively transform the integration interval without 
622 INTEGRATION OF EQUATIONS 
changing the value of the integral. The following example illustrates how this is done in 
practice. 
EXAMPLE 22.3 Two-Point Gauss-Legendre Formula 
I Problem Statement. Use Eq. (22.17) to evaluate the integral of 
I 
between the limits x = 0 to 0.8. Recall that this was the same problem that we solved in 
1 Chap. 21 using a variety of Newton-Cotes formulations. The exact value of the integral is 
1 1.640533. 
Solution. Before integrating the function, we must perform a change of variable so 
; that the limits are from -1 to +1. To do this, we substitute a = 0 and b = 0.8 into 
i Eq. (22.23) to yield 
The derivative of this relationship is [Eq. (22.24)] 
Both of these can be substituted into the original equation to yield 
Therefore, the right-hand side is in the form that is suitable for evaluation using Gauss 
quadrature. The transformed function can be evaluated at - l / f i to be equal to 0.516741 
and at 1/& to be equal to 1.305837. Therefore, the integral according to Eq. (22.17) is 
which represents a percent relative error of - 11.1 percent. This result is comparable in 
magnitude to a four-segment application of the trapezoidal rule (Table 21.1) or a single ap- 
plication of Simpson's 113 and 318 rules (Examples 21.4 and 21.6). This latter result is to 
be expected because Simpson's rules are also third-order accurate. However, because of 
the clever choice of base points, Gauss quadrature attains this accuracy on the basis of only 
two function evaluations. 
22.3.3 Higher-Point Formulas 
Beyond the two-point formula described in the previous section, higher-point versions can 
be developed in the general form 
22.3 GAUSS QUADRATURE 623 
TABLE 22.1 Weighting factors c and function arguments x used in Gauss-Legendre 
formulas. 
u- 
Weighting Funcgon Truncatisn 
Points Factors Arguments Error 
2 co = 1 .OOOOOOO xo = -0.577350269 =/141(5) 
ci = 1 0000000 X I = 0.577350269 
3 co = 0.5555556 xo = -0.774596669 =fi6ll() 
cl = 0.8888889 XI = 0.0 
q = 0.5555556 x2 = 0.774596669 
4 co = 0 3478548 xo = -0.86 1 1363 12 ~f 181[() 
cj = 0.6521 452 XI = -0 33998 1044 
c:! = 0.652 1452 X>= 0339981044 
cg = 0.3478548 X, = 0.861 1363 12 
co = 0.2369269 xo = -0.906 1 79846 =fl'Ol[{) 
ci = 0.4786287 xi = -0.5384693 10 
c:, = 0.5688889 x2 = 0.0 
cg = 0.4786287 x3 = 0.538469310 
~4 = 0.2369269 x4= 0.906179846 
co = 0.171 3245 xo = -0.9324695 14 = f 1 1 2)(~) 
cl = 0.360761 6 X I = -0.661 209386 
c2 = 0.4679 1 39 x, = -0.23861 91 86 
cg = 0.4679 1 39 x3 = 0.238619186 
cq = 0.360761 6 x4 = 0.661 209386 
cg = O 1713245 x5 = 0.9324695 14 
where n = the number of points. Values for c's and x's for up to and including the six-point 
formula are summarized in Table 22.1. 
EXAMPLE 22.4 Three-Point Gauss-Legendre Formula 
i Problem Statement. Use the three-point formula from Table 22.1 to estimate the inte- i 
I gral for the same function as in Example 22.3. 
! Solution. According to Table 22.1, the three-point formula is 
1 I = 0,5555556 f(-0.7745967) + 0.8888889 f(0) + 0.5555556 f(0.7745967) 
I 
which is equal to 
I = 0.2813013 + 0.8732444 + 0.4859876 = 1.640533 
"hich is exact. i 
Because Gauss quadrature requires function evaluations at nonuniformly spaced 
points within the integration interval, it is not appropriate for cases where the function is 
unknown. Thus, it is not suited for engineering problems that deal with tabulated data. 
However, where the function is known, its efficiency can be a decided advantage. This is 
particularly true when numerous integral evaluations must be performed. 
624 INTEGRATION OF EQUATIONS 
EXAMPLE 22.5 Applying Gauss Quadrature to the Falling Parachutist Problem 
I 
Problem Statement. In Example 21.3, we used the multiple-application trapezoidal 
i where g = 9.8, c = 12.5, and m = 68.1. The exact value of the integral was determined by calculus to be 289.435 1. Recall that the best estimate obtained using a 500-segment trape- 
I zoidal rule was 289.4348 with an 2 1.15 x lo-' percent. Repeat this computation 1 using Gauss quadrature. 
/ Solution. After modifying the function, the following results are obtained: 
Two-point estimate = 290.0145 
Three-point estimate = 289.4393 
Four-point estimate = 289.4352 
Five-point estimate = 289.4351 
Six-point estimate = 289.4351 
1 Thus, the five- and six-point estimates yield results that are exact to seven significant 
t 
' figures. i 
22.3.4 Error Analysis For Gauss Quadrature 
The error for the Gauss-Legendre formulas is specifiedgenerally by (Carnahan et al., 1969) 
where n = the number of points minus one and f (2n+2)(t) = the (2n + 2)th derivative of 
the function after the change of variable with 6 located somewhere on the interval from - 1 
to 1. Comparison of Eq. (22.26) with Table 21.2 indicates the superiority of Gauss quadra- 
ture to Newton-Cotes formulas, provided the higher-order derivatives do not increase sub- 
stantially with increasing n. Problem 22.8 at the end of this chapter illustrates a case where 
the Gauss-Legendre formulas perform poorly. In these situations, the multiple-application 
Simpson's rule or Romberg integration would be preferable. However, for many functions 
confronted in engineering practice, Gauss quadrature provides an efficient means for eval- 
uating integrals. 
22.4 IMPROPER INTEGRALS 
- - - 
To this point, we have dealt exclusively with integrals having finite limits and bounded in- 
tegrands. Although these types are commonplace in engineering, there will be times when 
improper integrals must be evaluated. In this section we will focus on one type of improper 
integral-that is, one with a lower limit of -m or an upper limit of +m. 
Such integrals usually can be evaluated by making a change of variable that transforms 
the infinite range to one that is finite. The following identity serves this purpose and works 
22.4 IMPROPER INTEGRALS 625 
for any function that decreases toward zero at least as fast as 1/.x2 as x approaches infinity: 
for ab > 0. Therefore, it can be used only when a is positive and b is a, or when a is -00 
and b is negative. For cases where the limits are from -cm to a positive value or from a 
negative value to cm, the integral can be implemented in two steps. For example, 
(22.28) 
where -A is chosen as a sufficiently large negative value so that the function has begun to 
approach zero asymptotically at least as fast as 1/x2. After the integral has been divided 
into two parts, the first can be evaluated with Eq. (22.27) and the second with a Newton- 
Cotes closed formula such as Simpson's 1/3 rule. 
One problem with using Eq. (22.27) to evaluate an integral is that the transformed 
function will be singular at one of the limits. The open integration formulas can be used to 
circumvent this dilemma as they allow evaluation of the integral without employing data at 
the end points of the integration interval. To allow the maximum flexibility, a multiple- 
application version of one of the open formulas from Table 21.4 is required. 
Multiple-application versions of the open formulas can be concocted by using closed 
formulas for the interior segments and open formrllas for the ends. For example, the multiple- 
segment trapezoidal rule and the midpoint rule can be combined to give 
In addition, semiopen formulas can be developed for cases where one or the other end 
of the interval is closed. For example, a formula that is open at the lower limit and closed 
at the upper limit is given as 
Although these relationships can be used, a preferred formula is (Press et at., 1986) 
which is called the extended midpoint rule. Notice that this formula is based on limits of in- 
tegration that are h/2 after and before the first and last data points (Fig. 22.8). 
FIGURE 22.8 
Placement of data points relative to integration limits for the extended midpoint rule. 
INTEGRATION OF EQUATIONS 
EXAMPLE 22.6 Evaluation of an Improper Integral 
: Problem Statement. The cumulative normal distribution is an important formula in 
statistics (see Fig. 22.9): 
FIGURE 22.9 
(a) The normal d~strlbutlon, (b) the transformed abscissa ~n terms of the standardized normal devi- 
ate, and (c) the cumulative normal distr~bution The shaded area in (a) and the po~nt ~n (c) repre 
sent the probab~llty that a ra~dom event will be less than the mean plus one standard dev~atron 
22.4 IMPROPER INTEGRALS 621 
where x = (y - E;)/sy is called the normalized standard deviate. It represents a change of 
variable to scale the normal distribution so that it is centered on zero and the distance along 
the abscissa is measured in multiples of the standard deviation (Fig. 22.9b). 
Equation (E22.6.1) represents the probability that an event will be less than x. For ex- 
ample, if x = 1, Eq. (E22.6.1) can be used to determine that the probability that an event 
will occur that is less than one standard deviation above the mean is N(l) = 0.8413. In 
other words, if 100 events occur, approximately 84 will be less than the mean plus one stan- 
dard deviation. Because Eq. (E22.6.1) cannot be evaluated in a simple functional form, it 
is solved numerically and listed in statistical tables. Use Eq. (22.28) in conjunction with 
Simpson's 1/3 rule and the extended midpoint rule to determine N(l) numerically. 
Solution. Equation (E22.6.1) can be re-expressed in terms of Eq. (22.28) as 
The first integral can be evaluated by applying Eq. (22.27) to give 
Then the extended midpoint rule with h = 1/8 can be employed to estimate 
Simpson's 1/3 rule with h = 0.5 can be used to estimate the second integral as 
Therefore, the final result can be computed as 
which represents an error of E* = 0.046 percent. 
The foregoing computation can be improved in a number of ways. First, higher-order 
formulas could be used. For example, a Romberg integration could be employed. Second, 
more points could be used. Press et al. (1986) explore both options in depth. 
Aside from infinite limits, there are other ways in which an integral can be improper. 
Common examples include cases where the integral is singular at either the limits or at a 
point within the integral. Press et al. (1986) provide a nice discussion of ways to handle 
these situations. 
628 INTEGRATION OF EQUATIONS 
PROBLEMS 
22.1 Use Romberg integration to evaluate 
to an accuracy of E, = 0.5%. Your results should be presented in 
the form of Fig. 22.3. Use the true value of 4.8333 to determine the 
true error E, of the result obtained with Romberg integration. Check 
that E, is less than the stopping criterion E,. 
22.2 Use order of hs Romberg integration to evaluate 
Compare E, and E,. 
22.3 Use Romberg integration to evaluate 
to an accuracy of the order of h8. Your results should be presented 
in the form of Fig. 22.3. 
22.4 Obtain an estimate of the integral from Prob. 22.1, but using 
two-, three-, and four-point Gauss-Legendre formulas. Compute E, 
for each case on the basis of the analytical solution. 
22.5 Obtain an estimate of the integral from Prob. 22.2, but using 
two-, three-, and four-point Gauss-Legendre formulas. Compute E, 
for each case on the basis of the analytical solution. 
22.6 Obtain an estimate of the integral from Prob. 22.3 using six- 
point Gauss-Legendre formulas. 
22.7 Perform the computation in Examples 21.3 and 22.5 for the 
falling parachutist, but use Romberg integration (E, = 0.01%). 
22.8 Employ two- through six-point Gauss-Legendre formulas to 
solve 
Interpret your results in light of Eq. (22.26). 
22.9 Use numerical integration to evaluate the following: 
Note that (d) is the normal distribution (recall Fig. 22.9). 
22.10 Develop a user-friendly computer program for the multiple- 
segment trapezoidal and Simpson's 1/3 rule based on Fig. 22.1. 
Test it by integrating 
Use the true value of 0.602297 to compute E, for n = 4. 
22.11 Develop a user-friendly computer program for Romberg in- 
tegration based on Fig. 22.4. Test it by duplicating the results of 
Examples 22.3 and 22.4 and the function in Prob. 22.10. 
22.12 Develop a user-friendly computer program for Gauss quad- 
rature. Test it by duplicating the results of Examples 22.3 and 22.4 
and the function in Prob. 22.10. 
22.13 Use the program developed in Prob. 22.11 to solve Probs. 
(a) 22.1, (b) 22.2, and (c) 22.3. 
22.14 Use the program developed in Prob. 22.12 to solve Probs. 
(a) 22.4, (b) 22.5, and (c) 22.6. 
22.15 Develop a program to implement the extended midpoint rule 
iteratively. Start the iterations with an initial estimatebased on a 
single point and the midpoint rule from Table 21.4. Then succes- 
sively apply Eq. (22.29) with the interval divided by 3 at each 
stage; that is, h = (b - a)/3 , (b - a)/9 , (b - a)/27, etc. Note that 
this means that one-third of the function estimates will have been 
determined in the previous iteration. Develop your algorithm so 
that it capitalizes on this property. Perform iterations until an ap- 
proximate error estimate EQ falls below a prespecified stopping cri- 
teria E ~ . Test the program by evaluating Example 22.6. 
03 
(b) a p Y sin2 dr 
CHAPTER 23 
We have already introduced the notion of numerical differentiation in Chap. 4. Recall that 
we employed Taylor series expansions to derive finite-divided-difference approximations 
of derivatives. In Chap. 4, we developed forward, backward, and centered difference ap- 
proximations of first and higher derivatives. Recall that, at best, these estimates had errors 
that were 0(h2)-that is. their errors were proportional to the square of the step size. This 
level of accuracy is due to the number of terms of the Taylor series that were retained dur- 
ing the derivation of these formulas. We will now illustrate how to develop more accurate 
formulas by retaining more terms. 
HIGH-ACCURACY DIFFERENTIATION FORMULAS 
As noted above, high-accuracy divided-difference formulas can be generated by including 
additional terms from the Taylor series expansion. For example, the forward Taylor series 
expansion can be written as [Eq. (4.21)] 
which can be solved for 
In Chap. 4, we truncated this result by excluding the second- and higher-derivative 
terms and were thus left with a final result of 
In contrast to this approach, we now retain the second-derivative term by substituting 
the following approximation of the second derivative [recall Eq. (4.24)] 
6301 NUMERICAL DIFFERENTIATION 
First Derivative 
Pix,] = fix,+11 - fixd 
h 
Error 
O!hl 
Second Derivat~ve 
Third Derivative 
FIGURE 23.1 
= -3fix1+41 + 14f1xt+3) - 24f(x;+2) + 1 Bf!x,+ij - 5fixJ 
Forward finite-divided-difference 
2h3 
formulas: two versions are pre- Fourth Derivative 
sented for each derivative. The 
latter version incorporates more mxd = f(x.+4l - 4fix.+31 + 6f(x,+21 - 4 f ( ~ ~ + l / + fbtl 
terms of the Taylor series expan- h4 
sion and is, consequently, more pix,l = -2f(x,+5l + 1 1 flx,+4l - 24fix,+31 + 26f[x,+2) - 14f(x,+i] + 3f(x,l 
accurate. h4 
FIGURE 23.2 
Backward finite-divided- 
difference formulas: two ver- 
sions are resented for each de- 
rivative. T g e latter version 
incorporates more terms of the 
Taylor series expansion and is, 
consequently, more accurate. 
First Derivative 
pix) - fix,) - fix,- 1 i 
I - 
h 
Second Derivative 
= 2fIx,l - 5fix,-1l + 4fIx,-zl - fix,-3) 
h2 
Third Derivative 
Fourth Derivative 
Error 
Oihl 
Oih21 
23.1 HIGH-ACCURACY DIFFERENTIATION FORMULAS 63 1 
First Derivative 
Second Derivative 
f , , = f + - 2flxJ + fix,-li 
h2 
Pix,) = -f(x,+21 + 16f(x,+1) - 30f(x,1 + 16fix,-11 - fix,-21 
1 2h2 
Third Derivative 
- f(x,+21 - 2f(x,+,I + 2f(x,-11 - fk-21 , - 
2 h3 
Fourth Derivative 
Error 
Olh21 
a h d l 
FIGURE 23.3 
Centered finite-divided-difference foimulas: two versions are presented for each derivative. The 
latter version incorporates more terms of the Taylor series expansion and is, consequently, more 
accurate. 
into Eq. (23.2) to yield 
f(xi+l) - f(xi) f(xi+i) - 2f(*i+i) + f(xi) + O ( h 2 ) 
f '(xi) = - h 2h2 
or, by collecting terms, 
Notice that inclusion of the second-derivative term has improved the accuracy to 
0(h2). Similar improved versions can be developed for the backward and centered forrnu- 
las as well as for the approximations of the higher derivatives. The formulas are summa- 
rized in Figs. 23.1 through 23.3 along with all the results from Chap. 4. The following ex- 
ample illustrates the utility of these formulas for estimating derivatives. 
EXAMPLE 23.1 High-Accuracy Differentiation Formulas 
Problem Statement. Recall that in Example 4.4 we estimated the derivative of 
f ( x ) = -0.1x4 - 0.15x3 - 0.5x2 - 0 . 2 5 ~ + 1.2 
632 NUMERICAL DIFFERENTIATION 
' at x = 0.5 using finite divided differences and a step size of h = 0.25, i 
Forward 
O(h) 
Backward Centered 
Q(h2) 
Estimate 
E i (%I 
where the errors were computed on the basis of the true value of -0.9125. Repeat this com- 
putation, but employ the high-accuracy formulas from Figs. 23.1 through 23.3. 
Solution. The data needed for this example are f 
I 
! The forward difference of accuracy 0(h2) is computed as (Fig. 23.1) 
i 
! The backward difference of accuracy O(h2) is computed as (Fig. 23.2) 
1 
I 
I 
The centered difference of accuracy 0(h4) is computed as (Fig. 23.3) 
1 As expected, the errors for the forward and backward differences are consrderably 
j more accurate than the results from Example 3.1 3. However. surprisingly, the centered dif- 
f ference yields a perfect result. This is because the formulas based on the Taylor series are 
i 
i equivalent to passing polynomials through the data points. 
23.2 RlCNARDSON EXTRAPOLATION 
To this point, we have seen that there are two way$ to improve derivative estimates when 
employing finite divided differences: (1) decrease the step size or (2) use a higher-order 
formula that employs more points. A third approach, based on Richardson extrapolation, 
uses two derivative estimates to compute a third, more accurate approximation. 
23.2 RICHARDSON EXTRAPOLATION 633 
Recall from Sec. 22.1.1 that Richardson extrapolation provided a means to obtain an 
improved integral estimate I by the formula [Eq. (22.4)l 
1 
I Z I(h2) + [I(hz) - I(hl)I (23.6) 
( h ~ l h z ) ~ - 1 
where I(hl) and I(h2) are integral estimates using two step sizes hl and h2. Because of its 
convenience when expressed as a computer algorithm, this formula is usually written for 
the case where h2 = h1/2, as in 
In a similar fashion, Eq. (23.7) can be written for derivatives as 
For centered difference approximations with 0(h2) , the application of this formula will 
yield a new derivative estimate of 0(h4). 
EXAMPLE 23.2 Richardson Extrapolation 
Problem Statement. Using the same function as in Example 23.1, estimate the first 
derivative at x = 0.5 employing step sizes of hl = 0.5 and h2 = 0.25. Then use Eq. (23.8) 
to compute an improved estimate with Richardson extrapolation. Recall that the true 
value is -0.9125. 
Solution. The first-derivative estimates can be computed with centered differences as 
i 
The improved estimate can be determined by applying Eq. (23.8) to give 
which for the present case is a perfect result. 
The previous example yielded a perfect result because the function being analyzed 
was a fourth-order polynomial. The perfect outcome was due to the fact that Richardson 
extrapolation is actually equivalent to fitting a higher-order polynomial through the data 
and then evaluating the derivatives by centered divided differences. Thus, the present case 
matched the derivative of the fourth-order polynomial precisely. For most other functions, 
of course, this would not occur and our derivative estimate would be improved but not 
perfect. Consequently, as was the case for the application of Richardson extrapolation, the 
634 NUMERICAL DIFFERENTIATION 
approach can be applied iteratively using a Romberg algorithm until the result falls below 
an acceptable error criterion. 
DERiVATIVES OF UNEQUALLY SPACED DATA 
The approaches discussed to this point are primarily designed to determine the derivative 
of a given function. For the finite-divided-difference approximations of Sec. 23.1, the data 
had to be evenly spaced. For the Richardson extrapolation technique of Sec. 23.2, the data 
had to be evenly spaced and generated for successively halved intervals. Such control of 
data spacing is usually available only in cases where we can use a function to generate a 
table of values. 
In contrast, empiricallyderived information-that is, data from experiments or field 
studies-is often collected at unequal intervals. Such information cannot be analyzed with 
the techniques discussed to this point. 
One way to handle nonequispaced data is to fit a second-order Lagrange interpolating 
polynomial [recall Eq. (18.23)] to each set of three adjacent points. Remember that this 
polynomial does not require that the points be equispaced. The second-order polynomial 
can be differentiated analytically to give 
where x is the value at which you want to estimate the derivative. Although this equation is 
certainly more complicated than the first-derivative approximations from Figs. 23.1 
through 23.3, it has some important advantages. First, it can be used to estimate the deriv- 
ative anywhere within the range prescribed by the three points. Second, the points them- 
selves do not have to be equally spaced. Third, the derivative estimate is of the same accu- 
racy as the centered difference [Eq. (4.22)]. In fact, for equispaced points, Eq. (23.9) 
evaluated at x = xi reduces to Eq. (4.22). 
EXAMPLE 23.3 Differentiating Unequally Spaced Data 
Problem Statement. As in Fig. 23.4, a temperature gradient can be measured down into 
the soil. The heat flux at the soil-air interface can be computed with Fourier's law, 
where q = heat flux (W/m2), k = coefficient of thermal diffusivity in soil (= 3.5 x 
lo-' m2/s), p = soil density (z 1800 kg/m3), and C = soil specific heat (= 840 J/(kg. "C)). 
Note that a positive value for flux means that heat is transferred from the air to the soil. Use 
numerical differentiation to evaluate the gradient at the soil-air interface and employ this 
estimate to determine the heat flux into the ground. 
23.4 DERIVATIVES AND INTEGRALS FOR DATA WITH ERRORS 635 
i 
FIGURE 23.4 
Temperature versus depth ~nto the soil 
i 
i 
8 
Solution. Equation (23.9) can be used to calculate the derivative as 
' which can be used to compute (note that 1 W = 1 J/s), 
I 
23.4 DERIVATIVES AND INTEGRALS FOR DATA WITH ERRORS 
Aside from unequal spacing, another problem related to differentiating empirical data is 
that it usually i~lcludes measurement error. A shortcoming of numerical differentiation is 
that it tends to amplify errors in the data. Figure 23.5a shows smooth, error-free data that 
when numerically differentiated yield a smooth result (Fig. 23.53). In contrast, Fig. 2 3 . 5 ~ 
uses the same data, but with some points raised and some lowered slightly. This minor 
modification is barely apparent from Fig. 23 .5~. However, the resulting effect in Fig. 23.5d 
is significant because the process of differentiation amplifies errors. 
As might be expected, the primary approach for determining derivatives for impre- 
cise data is to use least-squares regression to fit a smooth, differentiable function to the 
636 NUMERICAL DIFFERENTIATION 
FIGURE 23.5 
Illustration of how small data 
errors are amplified by nurneri- 
cal differentiation: [a) data with 
no error, (b ] the resulting numeri- 
cal differentlation, (c] data rnod- 
ified slightly, (dl the resulting 
differentiation manifesting in- 
creased variabil i~. In contrast, 
the reverse operation of integra- 
tion [moving from (dl to [c) by 
taking the area under (dl] tends 
to attenuate or smooth data 
errors 
data. In the absence of any other information, a lower-order polynomial regression might 
be a good first choice. Obviously, if the true functional relationship between the depen- 
dent and independent variable is known, this relationship should form the basis for the 
least-squares fit. 
23.4.1 DiFFerentiation versus Integration of Uncertain Data 
Just as curve-fitting techniques like regression can be used to differentiate uncertain data, 
a similar process can be employed for integration. However, because of the difference in 
stability between differentiation and integration, this is rarely done. 
As depicted in Fig. 23.5, differentiation tends to be unstable-that is, it amplifies er- 
rors. In contrast, the fact that integration is a summing process tends to make it very for- 
giving with regard to uncertain data. In essence, as points are summed to form an integral, 
random positive and negative errors tend to cancel out. In contrast, because differentiation 
is subtractive, random positive and negative errors tend to add. 
23.5 NUMERICAL INTEOWTION/DIFFERENTIATIOIV 
WlTH LlBRARlES AND PACKAGES 
Software libraries and packages have great capabilities for numerical integration and dif- 
ferentiation. In this section, we will give you a taste of some of the more useful ones. 
23.5 NUMERICAL INTEGRATION/DIFFERENTlATlON WITH LIBRARIES, PACKAGES 637 
Mathcad has operators that perform numerical integration and differentiation. These oper- 
ators employ and look like the same traditional mathematical symbols you have used since 
high school or your first semester of college. 
The integration operator uses a sequence of trapezoidal rule evaluations of the integral 
and the Romberg algorithm. Iterations are performed until successive results vary by less 
than TOL. The derivative operator uses a similar method to compute derivatives between 
order 0 and 5. This operator creates a table of approximations based on divided-difference 
calculations of the derivative using various orders and step sizes. Extrapolation techniques 
are used to estimate values for zero step size in a manner resembling Richardson's method. 
Figure 23.6 shows a Mathcad example where f(x) is created using the definition sym- 
bol (:=), and then the integral is calculated over a range from x = 0 to x = 0.8. In this case, 
Mathcad screen to determine 
the integral of a polynomial 
with Romberg integration. NUMERICALLY CALCULATE INTEGRALS 
Enter a function: 
f(x) := 0.2+25~x-200~~+675~-900~~~+400~~ 
Enter integration interval: 
a:=O 
f (x) dx = 1.64053333 
FIGURE 23.7 
Mathcad screen to implement 
numerical differentiation. 
x ) : = 2 . ~ + 3 + c o s ( x ) ~ 
Compute third derivative 
638 NUMERICAL DIFFERENTIATION 
we used the simple polynomial we used throughout Chap. 21. Note that the range as de- 
fined by the variables a and b is input using the definition symbol. 
Figure 23.7 shows a Mathcad example where f(x) is createdusing the definition symbol 
(:=) and then first and third derivatives are calculated at a point where x = -6. Note that the 
location of the point and the order of the derivative are input using the definition symbol. 
MATLAB has a variety of built-in functions that allow functions and data to be integrated 
and differentiated. The following example illustrates how a few of them can be used. 
EXAMPLE 23.4 Using MATLAB for Integration and Differentiation 
Problem Statement. Explore how MATLAB can be employed to integrate and differ- 
' entiate the function 
i 
a 
f(x) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5 
I 
from a = 0 to h = 0.8. Recall from Chaps. 21 and 22 that the true value of the integral can 
be determined analytically to be 1.640533. 
I 
4 
1 Solution. First, we will use the MATLAB function quad to integrate the function. To 
I use quad, we must first develop an M-file to hold the function. Using a text editor, we can 
j create the following file: 
i 
1 f u n c t i o n y = f x ( x ) 
! y=0.2+25*x-200*x.A2+675*x.A3-900*x.A4+400*x.A5; 
' It can then be stored on the MATLAB directory as fx.m. 
After entering MATLAB, we can invoke quad by typing 
i 
1 where the second two entries are the integration limits. The result is 
I 
1 Thus, MATLAB provides an accurate integral estimate. 
1 Next, let's investigate how MATLAB handles integrals of tabular data. To do this, 
we will repeat Example 21.7, where we sampled the function at unequal intervals (recall 
' Table 21.3). We can generate the same information on MATLAB by first defining the val- I ues of the independent variable, 
i Then, we can generate a vector y containing the corresponding dependent variable values 
by invoking fx, 
i 
I Y= 
I 
Columns 1 t h r o u g h 7 
0.2000 1.3097 1.3052 1.7434 2.0749 2.4560 2.8430 
23.5 NUMERICAL INTEGRATIQN/DIFFERENTlATlQN WITH LIBRARIES, PACKAGES 639 
Co lumns 8 t h r o u g h 11 
3 . 5 0 7 3 3 . 1 8 1 9 2 . 3 6 3 0 0 . 2 3 2 0 
We can integrate these values by invoking the trapz function, 
i n t e g r a l = 
1 . 5 9 4 8 
Thus, as the name implies, trapz, applies the trapezoidal rule to each interval and sums the 
results to obtain the total integral. 
Finally, we can differentiate the unequally spaced data in x and y. To do this we use the 
diff function, which merely determines the differences between adjacent elements of a vec- 
tor, for example, 
a n s = 
C o l u m n s 1 t h r o u g h 7 
0 . 1 2 0 0 0 . 1 0 0 0 0 . 1 0 0 0 0 . 0 4 0 0 0 . 0 4 0 0 0 . 0 4 0 0 0 . 1 0 0 0 
C o l u m n s 8 t h r o u g h 1 0 
0 . 1 0 0 0 0 . 0 6 0 0 0 . 1 0 0 0 
The result represents the differences between each pair of elements of x. To compute 
divided-difference approximations of the derivative, we merely perform a vector division 
of the y differences by the x differences by entering 
which yields 
d = 
C o l u m n s 1 t h r o u g h 7 
9 . 2 4 7 7 - 0 . 0 4 4 9 4 . 3 8 1 5 8 . 2 8 7 7 9 . 5 2 7 4 9 . 6 7 4 6 6 . 6 4 3 1 
C o l u m n s 8 t h r o u g h 1 0 
- 3 . 2 5 3 7 - 1 3 . 6 4 8 8 - 2 1 . 3 1 0 0 
These represent crude estimates of the derivatives for each interval. This approach could be 
refined by using finer spacing. 
IMSL has several routines for integration and differentiation (Table 23.1). In the present 
discussion, we will focus on the QDAG routine. This routine integrates a function using a 
globally adaptive scheme based on Gauss-Kronrod rules. 
QDAG is implemented by the following CALL statement: 
C A L L Q D A G (F, A, B, E R R A B S , E R R R E L , I R U L E , R E S U L T , E R R E S T ) 
where F = User-supplied FUNCTION to be integrated. The form is F(X), where X is the 
independent variable. Note that F must be declared EXTERNAL in the calling 
program. 
W o NUMERICAL DIFFERENTIATION 
TABLE 23.1 IMSL routines to integrate and differentiate. 
Category Routines Capabiliw 
Unvarate 
quodrature 
QDAGS 
QDAG 
QDAGP 
QDAGI 
QDAWO 
QDAWF 
QDAWS 
QDAWC 
Q D N G 
Adaptive general-purpose end-point singularities 
Adaptive general purpose 
Adaptive general-purpose ponts of sngularlty 
Adapt~ve general-purpose infinlte interval 
Adaptive weighted oscillatory (trigonometr~c] 
Adaptive weighted Fourier (trigonometric] 
Adaptive weighted algebraic end-point singularities 
Adaptive weighted Couchy princ~pal value 
Nonadaptive general purpose 
Multidimensional 
quadrature 
TWODQ Two-dimensional quadrature [iterated integral) 
QAND Adaptive N-dimensional quadrature over a hyper-rectangle 
Gauss rules and 
three-term 
recurrences 
GQRUL Gauss quadrature rule for classical we~ghts 
GQRCF Gauss quadrature rule from recurrence coeffcients 
RECCF Recurrence coefficients for classlcol weights 
RECQR Recurrence coefficients from quadrature rule 
FQRUL Fejer quadrature rule 
Differentlation 
OtKiV Approximation to first, second, or third derrvat~ve 
A = Lower limit of integration. (Input) 
B = Upper limit of integration. (Input) 
ERRABS = Absolute accuracy desired. (Input) 
ERRREL = Relative accuracy desired. (Input) 
IRULE = Choice of quadrature rule. (Input). IRULE = 2 is recommended for 
most functions. If the function has a peak singularity, use IRULE = 1. 
If the function is oscillatory, use IRULE = 6. 
RESULT = Estimate of the integral from A to B of F. (Output) 
ERREST = Estimate of the absolute value of the error. (Output) 
EXAMPLE 23.5 Using IMSL to Integrate a Function 
1 Problem Statement. Use QDAG to determine the integral of 
f(x) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5 
from a = 0 to b = 0.8. Recall from Chaps. 21 and 22 that the exact value of the integral I can be determined analytically to be 1.640533. 
PROBLEMS 64 1 
Solution. An example of a main Fortran 90 program and function using QDAG to solve 
this problem can be written as 
PROGRAM I n t e g r a t e 
USE m i m s l 
I M P L I C I T NONE 
I N T E G E R : :i r u L e = l 
R E A L : : a = O . , b = 0 . 8 , e r r a b s = 0 . 0 . e r r r r e l = O O O 0 1 
R E A L : : e r r e s t , r e s , f 
EXTERNAL f 
CALL QDAG ( f , a , e r r a b s , e r r r e l , i r u L e , r e s , e r r e s t ) 
P R I N T ' ( I ' C o m p u t e d = " , F 8 . 4 ) ' , r e s 
P R I N T ' ( " E r r o r e s t i m a t e = " , 1 P E 1 0 . 3 ) ' , e r r e s t 
END PROGRAM 
FUNCTION f ( x ) 
I M P L I C I T NONE 
REAL: :x,f 
f = 0 . 2 + 2 5 . * X - 2 0 0 . * X * * 2 + 6 7 5 . * X * * 3 . - 9 0 0 . * X * * 4 + 4 ~ O ~ * X * * 5 
END FUNCTION 
Output: 
C o m p u t e d = 1 . 6 4 0 5 
E r r o r e s t i m a t e = 5 . 0 0 0 E - 0 5 
PROBLEMS 
23.1 Compute forward and backward difference approximations 
of O(h) and 0(h2), and central difference approximations of 0(h2) 
and O(h4) for the first derivative of y = sin x at x = n/4 using a 
value of h = n/12. Estimate the true percent relative error e, for 
each approximation. 
23.2 Repeat Prob. 23.1, but for y = log x evaluated at x = 20 with 
h = 2. 
23.3 Use centered difference approximations to estimate the first 
and second derivatives of y = eX at x = I for h = 0.1. Employ both 
0(h2) and 0(h4) formulas for your estimates. 
23.4 Use Richardson extrapolation to estimate the first derivative 
of y = sin x at x = n/4 using step sizes of hl = n/3 and h2 = n/6. 
Employ centered differences of O(h2) for the initial estimates. 
23.5 Repeat Prob. 23.4, but for the first derivative of In x at x = 4 
using hl = 2 and hZ = 1. 
23.6 Employ Eq. (23.9) to determine the first derivative of 
y = 32r4 - 7x3 - 10x - 8 at x = 0 based on values at xo = -0.5, 
xl = 1, and x2 = 2. Compare this result with the true value and 
with an estimate obtained using a centered difference approxima- 
tion based on h = 1. 
23.7 Prove that for equispaced data points, Eq. (23.9) reduces to 
Eq. (4.221 at x = x,. 
23.8 Compute the first order central difference approximations of 
O(h4) for each of the following functions at the specified location 
and for the specified step size: 
(a) y = x 3 + 4 x - 15 a t x = O , h = 0 . 5 
(b) y = x2 cos x at x = 0.5, h = 0.2 
(e) y = tan (x/3) a t x = 3 , k=O. I 
(d)y=sin(0 .5&)/x a t x = l , h = 0 . 1 
(e) y = er + x a t x = 1, h=0 .25 
W 2 NUMERICAL DIFFERENTIATION 
23.9 The following data were collected for the distance traveled and the distance traveled can be obtained by 
versus time for a rocket: 
(a) Use Mathcad to integrate Eq. (P23.12~) from t = 0 to 10. 
Use numerical differentiation to estimate the rocket's velocity and (b) Analytically integrate Eq. (P23.12b) with the initial condition 
acceleration at each time. that d = 0 at t = 0. Evaluate the result at t = 10 to confirm (a). 
23.10 Develop a user-friendly subprogram to apply a Romberg al- (c) differentiate Eq' (P23'12b) at = lo' 
gorithm to estimate the derivative of a given function. (d) Evaluate Eq. (P23.12~) at t = 10 to confirm (c). 
23.11 Develop a user-friendly subprogram to obtain first-deriva- 23.13 The normal distribution is defined as 
tive estimates for unequally spaced data. Test it with the following 1 
f ( x ) = - e - x 2 / 2 data: 6 
x 1 1.2 3 3 2 5 7 (a) Use Mathcad to integrate this function from x = - 1 to 1 and 
ijxj 1 1.807 0.7468 0.6522 0 1684 0.03192 from -2 to 2. 
(b) Use Mathcad to determine the inflection point of this function. 
where f ( x ) = Se-'x. Compare your results with the true deem- Because the function is symmetric, limit your analysis to posi- 
tives. tive x. 
23.12 Recall that for the falling parachutist problem, the velocity 23.14 Use MATLAB'S quad function to integrate Eq. (P23.12~) 
is given by from t = 0 to 10. 
23.15 The following data were generated from the normal distribution: 
(a) Use MATLAB to integrate this data from x = - 1 to 1 and -2 to2. 
(b) Use MATLAB to estimate the inflection points of this data. 
23.16 Use IMSL to integrate the normal distribution (see 
Prob. 23.13) from x = -1 to 1, from -2 to 2, and from -3 to 3. 
CHAPTER 24 
The purpose of this chapter is to apply the methods of numerical integration and differen- 
tiation discussed in Part Six to practical engineering problems. Two situations are most fre- 
quently encountered. In the first case, the function under study can be expressed in analytic 
form but is too complicated to be readily evaluated using the methods of calculus. Numer- 
ical methods are applied to situations of this type by using the analytic expression to gen- 
erate a table of arguments and function values. In the second case, the function to be eval- 
uated is inherently tabular in nature. This type of function usually represents a series of 
measurements, observations, or some other empirical information. Data for either case is 
directly compatible with several schemes discussed in this part of the book. 
Section 24.1, which deals with heat calculations from chemical engineering, involves 
equations. In this application, an analytic function is integrated numerically to determine 
the heat required to raise the temperature of a material. 
Sections 24.2 and 24.3 also involve functions that are available in equation form. Sec- 
tion 24.2, which is taken from civil engineering, uses numerical integration to determine 
the total wind force acting on the mast of a racing sailboat. Section 24.3 determines the 
root-mean-square current for an electric circuit. This example is used to demonstrate the 
utility of Romberg integration and Gauss quadrature. 
Section 24.4 focuses on the analysis of tabular information to determine the work re- 
quired to move a block. Although this application has a direct connection with mechanical 
engineering, it is germane to all other areas of engineering. Among other things, we use this 
example to illustrate the integration of unequally spaced data. 
24.1 BNTEGlDATlOlU TO DETERMINE THE TOTAL QlBARlTlTV 
OF HEAT (CHEMICAL/PETRObEUM ENGINEERING) 
Background. Heat calculations are employed routinely in chemical and petroleum engi- 
neering as well as in many other fields of engineering. This application provides a simple 
but useful example of such computations. 
One problem that is often encountered is the determination of the quantity of heat 
required to raise the temperature of a material. The characteristic that is needed to carry out 
this computation is the heat capacity r. This parameter represents the quantity of heat 
644 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
required to raise a unit mass by a unit temperature. If c is constant over the range of tem- 
peratures being examined, the required heat AH (in calories) can be calculated by 
where c has units of cal/(g."C), m = mass (g), and A T = change in temperature ("C). For 
example, the amount of heat required to raise 20 g of water from 5 to 10°C is equal to 
A H = 20(1)(10 - 5) = 100 cal 
where the heat capacity of water is approximately 1 cal/(g."C). Such a computation is ad- 
equate when A T is small. However, for large ranges of temperature, the heat capacity is not 
constant and, in fact, varies as a function of temperature. For example, the heat capacity of 
a material could increase with temperature according to a relationship such as 
In this instance you are asked to compute the heat required to raise 1000 g of this material 
from - 100 to 200°C. 
Solution. Equation (PT6.4) provides a way to calculate the average value of c ( T ) : 
which can be substituted into Eq. (24.1) to yield 
where AT = T2 - T I . NOW because, for the present case, c ( T ) is a simple quadratic, AH 
can be determined analytically. Equation (24.2) is substituted into Eq. (24.4) and the result 
integrated to yield an exact value of A H = 42,732 cal. It is useful and instructive to com- 
pare this result with the numerical methods developed in Chap. 21. To accomplish this, it 
is necessary to generate a table of values of c for various values of T: 
These points can be used in conjunction with a six-segment Simpson's 1 /3 rule to compute 
an integral estimate of 42,732. This result can be substituted into Eq. (24.4) to yield a value 
of AH = 42,732 cal, which agrees exactly with the analytical solution. This exact agree- 
ment would occur no matter how many segments were used. This is to be expected because 
c is a quadratic function and Simpson's rule is exact for polynomials of the third order or 
less (see Sec. 21.2.1). 
24.2 EFFECTIVE FORCE ON THE MAST OF A RACING SAILBOAT 645 
TABLE 24.1 Results using the trapezoidal rule with varl- 
ous step sizes. 
v= - - - *w - " IBL - -_gU_- - 
Step Size, "C AH PI (%) 
303 96,048 125 
150 43 029 C 7 
1 CO 4'2,864 0 3 
50 42,755 0 07 
2 5 42,740 0018 
10 42 733 3 <0 01 
5 42 732 3 <O 01 
1 42 732 01 <0 01 
0 05 42,732 00001 <0 01 
The results using the trapezoidal rule are listed in Table 24.1. It is seen that the trape- 
zoidal rule is also capable of estimating the total heat very accurately. However, a small 
step (< 10'C) is required for five-place accuracy. This example is a good illustration of 
why Simpson's rule is very popular. It is casy to perform with either a hand calculator or. 
better get. a personal computer. In addition, it is usually sufficiently accurate for relatively 
large step sizes and is exact for polynomials of third order or less. 
24.2 EFFECTIVE FORCE ON THE MAST OF A MCING SAILBOAT 
(CIVILiENVIRONMENTAL ENGINEERING) 
Background. A cross section of a racing sailboat is shown in Fig. 24 .1~ . Wind forces ( f ) 
exerted per foot of mast from the sails vary as a function of distance above the deck of the 
boat (z), aa in Fig. 24.lh. Calculate the tensile force T in the left mast support cable. 
FIGURE 24.1 
(a) Cross sectior, 01 a racing 
sailboaL. (D) W~nd fol-ces f ex- 
erted per foot of mast as a 
functior of dijta-rce z above ttle 
deck of the boat. 
646 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
assuming that the right support cable is completely slack and the mast joins the deck in a 
manner that transmits horizontal or vertical forces but no moments. Assume that the mast 
remains vertical. 
Solution. To proceed with this problem, it is required that the distributed force f be con- 
verted to an equivalent total force F and that its effective location above the deck d be cal- 
culated (Fig. 24.2). This computation is complicated by the fact that the force exerted per 
foot of mast varies with the distance above the deck. The total force exerted on the mast can 
be expressed as the integral of a continuous function: 
This nonlinear integral is difficult to evaluate analytically. Therefore, it is convenient to 
d = 13.05 ft employ numerical approaches such as Simpson's rule and the trapezoidal rule for this prob- i T ! lem. This is accomplished by calculating f(z) for various values of z and then using 
Eq. (21.10) or (21.18). For example, Table 24.2 has values of f ( z ) for a step size of 3 ft that 
I provide data for Simpson's 113 rule or the trapezoidal rule. Results for several step sizes 
I 
0 r - - - H are given in Table 24.3. It is observed that both methods give a value of F = 1480.6 lb as 
I 3fi the step size becomes small. In this case, step sizes of 0.05 ft for the trapezoidal rule and 
I 
I/ 0.5 ft for Simpson's rule provide good results. 
FIGURE 24.2 
Free-body diagram of the forces 
exerted on )he of a sail. TABLE 24.2 Values of f(z) For a step size of 3 ft that provide data for the trapezoidal rule 
boat . and Simpson's 1 /3 rule. 
TABLE 24.3 Values of F computed on the basis of various versions of the trapezoidal rule 
and Simpson's 1/3 rule. 
Technique Skp size, f-t Segments 
Tra~ezoidal rule 
Simpson's 1 / 3 rule 
24.3 ROOT-MEAN-SQUARE CURRENT BY NUMERICAL INTEGRATION 647 
The effective line of action of F (Fig. 24.2)can be calculated by evaluation of the 
integral 
d = (24.6) 
30 
2002 [ z / (5 + z ) ] e-2Z/30 d z 
d = 
1480.6 
This integral can be evaluated using methods similar to the above. For example, Simpson's 
1/3 rule with a step size of 0.5 gives d = 19,326.911480.6 = 13.05 ft. 
With F and d known from numerical methods, a free-body diagram is used to develop 
force and moment balance equations. This free-body diagram is shown in Fig. 24.2. Sum- 
ming forces in the horizontal and vertical direction and taking moments about point 0 gives 
C F H = O = F - T s i n 8 - H (24.8) 
where T = the tension in the cable and H and V = the unknown reactions on the mast 
transmitted by the deck. The direction, as well as the magnitude, of H and V is unknown. 
Equation (24.10) can be solved directly for V because F and dare known. 
Therefore, from Eq. (24.9), 
V 6440.6 
T = ---- - - - = 6473 Ib 
cos 8 0.995 
and from Eq. (24.8), 
H = F - T sin 8 = 1480.6 - (6473)(0.0995) = 836.54 lb 
These forces now enable you to proceed with other aspects of the structural design of the 
boat such as the cables and the deck support system for the mast. This problem illustrates 
nicely two uses of numerical integration that may be encountered during the engineering 
design of structures. It is seen that both the trapezoidal rule and Simpson's 1/3 rule are easy 
to apply and are practical problem-solving tools. Simpson's 1/3 rule is more accurate than 
the trapezoidal rule for the same step size and thus may often be preferred. 
24.3 ROOT-MEAN-SQUARE CURRENT BY NUMERICAL 
IIVTEGMTION (ELECTRICAL ENGINEERING) 
Background. The average value of an oscillating electric current over one period may be 
zero. For example, suppose that the current is described by a simple sinusoid: i ( t ) = 
sin (2n /T) , where T is the period. The average value of this function can be determined by 
648 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
FIGURE 24.3 
A varying 
current 
electric 
the following equation: 
T 
sin (F) dr 
I = - - - cos (2n) i- cos 0 
T = 0 T - 0 
Despite the fact that the net result is zero, such current is capable of performing work 
and generating heat. Therefore, electrical engineers often characterize such current by 
where i(t) = the instantaneous current. Calculate the RMS or root-mean-square current of 
the waveform shown in Fig. 24.3 using the trapezoidal rule, Simpson's 1/3 rule, Romberg 
integration, and Gauss quadrature for T = 1 s. 
Solution. Integral estimates for various applications of the trapezoidal rule and Simp- 
son's 1/3 rule are listed in Table 24.4. Notice that Simpson's rule is more accurate than the 
trapezoidal rule. 
The exact value for the integral is 15.41261. This result is obtained using a 128- 
segment trapezoidal rule or a 32-segment Simpson's rule. The same estimate is also deter- 
mined using Romberg integration (Fig. 24.4). 
In addition, Gauss quadrature can be used to make the same estimate. The determina- 
tion of the root-mean-square current involves the evaluation of the integral (T = I) 
1 /2 
I = (10e-' sin 2nt)' d t 
24.3 ROOT-MEAN-SQUARE CURRENT BY NUMERICAL INTEGRATION 649 
TABLE 24.4 Values for the integral calculated using various numerical schemes. The 
percent relative error F , is based on a true value of 15.41 261 . 
Techniaue Segments Integral 
Trapezoidal rule 
Simpson's 113 rule 
FlGLdRE 24.4 
Result of using Romberg Q(h2) 0(h4) O(hb) O(h8) 0(hl0) 0(h12) 
integration to estimate the 0 20.2 1 769 15.16503 15.41502 
RMS current. 
1541261 15.41261 
15.16327 15.48082 15.41 1 1 1 15.41262 15.41261 
15.40143 15.41547 15.41 225 15 41 261 
15.41 196 15.4 1277 15.41261 
15.4 1257 15.41 262 
15.41261 
First, a change in variable is performed by applying Eqs. (22.23) and (22.24) to yield 
These relationships can be substituted into Eq. (24.12) to yield 
For the two-point Gauss-Legendre formula, this function is evaluated at t d = -l/& 
and I/&, with the results being 7.684096 and 4.313728, respectively. These values can 
be substituted into Eq. (22.17) to yield an integral estimate of 11.99782, which represents 
an error of E~ = 22.1 %. 
The three-point formula is (Table 22.1) 
650 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
TABLE 24.5 Results of using various-point Gauss 
quadrature formulas to approximate 
the integral. 
Points Estimate 8, (%) 
The results of using the higher-point formulas are summarized in Table 24.5. 
The integral estimate of 15.41261 can be substituted into Eq. (24.12) to compute an 
IRMs of 3.925890 A. This result could then be employed to guide other aspects of the de- 
sign and operation of the circuit. 
24.4 NUMERICAL IIUTEGMTION TO COMPUTE WORK 
(MECHANICAL/AERBSPAGE ENGINEERING) 
Background. Many engineering problems involve the calculation of work. The general 
formula is 
Work = force x distance 
When you were introduced to this concept in high school physics, simple applications were 
presented using forces that remained constant throughout the displacement. For example, 
if a force of 10 lb was used to pull a block a distance of 15 ft, the work would be calculated 
as 150 ft . lb. 
Although such a simple computation is useful for introducing the concept, realistic 
problem settings are usually more complex. For example, suppose that the force varies dur- 
ing the course of the calculation. In such cases, the work equation is re-expressed as 
(24.14) 
where W = work (ft . lb), xo and x,, = the initial and final positions, respectively, and F(x) 
a force that varies as a function of position. If F(x) is easy to integrate, Eq. (24.14) can be 
evaluated analytically. However, in a realistic problem setting, the force might not be ex- 
pressed in such a manner. In fact, when analyzing measured data, the force might be avail- 
able only in tabular form. For such cases, numerical integration is the only viable option for 
the evaluation. 
Further complexity is introduced if the angle between the force and the direction of 
movement also varies as a function of position (Fig. 24.5). The work equation can be mod- 
ified further to account for this effect, as in 
(24.15) 
F(x) cos [ ~ ( x ) ] dx 
Again, if F(x) and O(x) are simple functions, Eq. (24.15) might be solved analytically. How- 
ever, as in Fig. 24.5, it is more likely that the functional relationship is complicated. For this 
situation, numerical methods provide the only alternative for determining the integral. 
a 24.4 NUMERICAL INTEGRATION TO COMPUTE WORK 65 1 
FIGURE 24.5 
The case of a variable force 
acting on a block. For this case, 
the angle, as well as the rnagni- 
tude, of the force varies. 
TABLE 24.6 Data for force F(x) and angle O(x) as a function 
of position x. 
- ~ ~ > - z ~ * * ~ * v e - ~ * * * - ~ - w ~ - ~ ~ ~ - ~ ~ ~ 
2% ff Fb), ib 8, rad F(x) eos 8 
Suppose that you have to perform the computation for the situation depicted in 
Fig. 24.5. Although the figure shows the continuous values for F(x) and Q(x), assume that, 
because of experimental constraints, you are provided with only discrete measurements at 
x = 5-ft intervals (Table 24.6). Use single- and multiple-application versions of the trape- 
zoidal rule and Simpson's 113 and 318 rules to compute work for this data. 
652 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
Solution. The results of the analysis are summarized in Table 24.7. A percent relative 
error st was computed in reference to a true value of the integral of 129.52 that was esti- 
mated on the basis of values taken from Fig. 24.5 at 1-ft intervals. 
The results are interesting because the most accurate outcome occurs for the simple 
two-segment trapezoidal rule. More refined estimates using more segments, as well as 
Simpson's rules, yield less accurate results. 
The reason for this apparently counterintuitive result is that the coarse spacing of the 
points is not adequate to capture the variations ofthe forces and angles. This is particularly 
evident in Fig. 24.6, where we have plotted the continuous curve for the product of F(x) 
and cos [O(x)]. Notice how the use of seven points to characterize the continuously varying 
function misses the two peaks at x = 2.5 and 12.5 ft. The omission of these two points 
TABLE 24.7 Estimates of work calculated using the trapezoidal rule and Simpson's rules. 
The percent relative error E, as computed in reference to a true value of the 
integral (1 29.52 ft .Ib) that was estimated on the basis of values at 1 -ft 
intervals. 
Technique Segments W r k E,, % 
Trapezoidal 1 5.31 95.9 
2 133.19 2.84 
3 124 98 3.5 1 
6 1 19.09 8.05 
Simpson's 1 /3 rule 2 175.82 -35.75 
6 117 13 9.57 
Simpson's 3/8 rule 3 139 93 -8 04 
FlGURE 24.6 
A continuous of 
Fix) cos [O(x)] versus position 
with the seven discrete points 
used to develop the numerical 
integration estimates in Table 
24.7. Notice how the use of 
seven points to characterize thi 
continuously varying function 
misses two peaks at x = 2.5 
and 12.5 ft. 
24.4 NUMERICAL INTEGRATION TO COMPUTE WORK 653 
effectively limits the accuracy of the numerical integration estimates in Table 24.7. The 
fact that the two-segment trapezoidal rule yields the most accurate result is due to the 
chance positioning of the points for this particular problem (Fig. 24.7). 
The conclusion to be drawn from Fig. 24.6 is that an adequate number of measure- 
ments must be made to accurately compute integrals. For the present case, if data were 
FIGURE 24.7 
Graphical depiction of why the 
two-segment trapezoidal rule 
yields a good estimate of the in- 
tegral for this particular case. 
By chance, the use of two 
trapezoids happens to lead to 
an even balance between posi- 
tive and negative errors. 
FIGURE 24.8 
The unequal segmentation 
scheme that results from the in- 
clusion of two additional points 
at x = 2.5 and 12.5 in the 
data in Table 24.6. The numeri- 
cal integration formulas applied 
to each set of segments are 
shown. 
654 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
available at F(2.5) cos [0(2.5)] = 4.3500 and F(12.5) cos [Q(12.5)] = 11.3600, we could 
determine an integral estimate using the algorithm for unequally spaced data described pre- 
viously in Sec. 21.3. Figure 24.8 illustrates the unequal segmentation for this case. Includ- 
ing the two additional points yields an improved integral estimate of 126.9 ( E , = 2.02%). 
Thus, the inclusion of the additional data would incorporate the peaks that were missed pre- 
viously and, as a consequence, lead to better results. 
PROBLEMS 
Chemical/Petroleum Engineering c = concentration, and x = distance (cm). An environmental engi- 
24.1 Perform the same computation as Sec. 24.1, but compute the neer measures the following concentration of a pollutant in the sed- 
amount of heat required to raise the temperature of 1500 g of the iments underlying a lake (x = 0 at the sediment-water interface and 
material from -150 to 50°C. Use Simpson's rule for your computa- increases downward): 
tion, with values of T at 50°C increments. 
24.2 Repeat Prob. 24.1, but use Romberg integration to E, = '' Crn 0 1 3 
0.01%. c, 1 o - ~ c ~ / c r n ~ 1 0 1 0 4 0 9 
24.3 Repeat Prob. 24.1, but use a two- and a three-point Gauss- 
Legendre formula. Interpret your results. 
24.4 Integration provides a means to compute how much mass en- 
ters or leaves a reactor over a specified time period, as in 
where tl and t2 = the initial and final times. This formula makes in- 
tuitive sense if you recall the analogy between integration and sum- 
mation. Thus, the integral represents the summation of the product 
of flow times concentration to give the total mass entering or leav- 
ing from tl to t2. If the flow rate is constant, Q can be moved out- 
side the integral: 
Use numerical integration to evaluate this equation for the data in 
Table P24.4. Note that Q = 4 m3/min. 
24.5 The outflow chemical concentration from a completely 
mixed reactor is measured as 
t, rnin 1 0 2 4 6 8 121620 
c, rng/rn3 I 10 20 30 40 60 72 70 50 
For an outflow of Q = 12 m3/min, estimate the mass of chemical 
that exits the reactor from t = 0 to 20 min. 
24.6 Fick'sfirst diffusion law states that 
dc 
Mass flux = -D - 
dx 
Use the best numerical differentiation technique available to esti- 
mate the derivative at x = 0. Employ this estimate in conjunction 
with Eq. (P24.6) to compute the mass flux of pollutant out of the 
sediments and into the overlying waters (D = 2 x cm2/s). For 
a lake with 3 x lo6 m2 of sediments, how much pollutant would be 
transported into the lake over a year's time? 
24.7 The following data were collected when a large oil tanker was 
loading: 
t , min I 0 15 30 45 60 90 120 
V, lo6 barrels 1 0.5 0 65 0.73 0.88 1.03 1 14 1.30 
Calculate the flow rate Q (that is, dV/dt) for each time to the order 
of h2. 
TABLE P24.4 Values of concentration measured in the exit 
pipe of a reactor. 
43 J L 
where mass flux = the quantity of mass that passes across a unit 50 34 
area per unit time (g/cm2/s), D = a diffusion coefficient (cm2/s), 
PROBLEMS 655 
FlGldRE P24.12 
A stream cross section 
24.8 You are interested in measuring the fluid velocity in a narrow 
rectangular open channel carrying petroleum waste between loca- 
tions in an oil refinery. You know that, because of bottom friction, the 
velocity varies with depth in the channel. If your technician has time 
to perform only two velocity measurements, at what depths would 
you take them to obtain the best estimate of the average velocity? 
State your recommendation in terms of the percent of total depth d 
measured from the fluid surface. For example, measuring at the top, 
would be O%d, whereas at the very bottom would be 10O%d. 
CivilIEnvironmental Engineering 
24.9 Perforrn the same computation in Sec. 24.2, but use order of 
h8 Romberg integration to evaluate the integral. 
24.10 Perform the same computation as in Sec. 24.2, but use Gauss 
quadrature to evaluate the integral. 
24.11 Compute F using the Trapezoidal rule, Simpson's 113 and 
Simpson's 318 rule. Divide the mast into five-foot intervals. 
24.12 Stream cross-sectional areas (A) are required for a number 
of tasks in water resources engineering, including flood forecasting 
and reservoir design. Unless electronic sounding devices are avail- 
able to obtain continuous profiles of the channel bottom, the engi- 
neer must rely on discrete depth measurements to compute A. An 
example of a typical stream cross section is shown in Fig. P24.12. 
The data points represent locations where a boat was anchored and 
depth readings taken. Use two trapezoidal rule applications ( h = 4 
and 2 m) and Simpson's 1/3 rule to estimate the cross-sectional 
area from this data. 
24.13 During a survey, you are required to compute the area of the 
field shown in Fig. P24.13. Use Simpson's rules to determine the 
area. 
24.14 A transportation engineering study requires the calculation 
of the total number of cars that pass through an intersection over a 
24-h period. An individual visits the intersection at various times 
during the course of a day and counts the number of cars that pass 
through the intersection in a minute. Utilize this data, summarized 
in Table P24.14, to estimate the total number of cars that pass 
through the intersection per day. (Be careful of units.) 
24.15 A wind force distributed against the side of a skyscraper is 
measured as 
Height I, m 10 30 60 90 120 150 180 210 240 
Force ,F ( / ] ,N /m 10 350 1000 1500 2600 3000 3300 3500 3600 
Compute the net force and the line of action due to this distributed 
wind. 
656 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
FlGlDRE P24.13 
A field bounded by two roads and a creek 
TABLE P24.14 Traffic flow rate (cars/min) for an intersection measured at various times 
over a 24-h period. 
WIN-*-------- 
6:00 A.M. 
PROBLEMS 657 
FIGURE P24.16Water exerting pressure on the upstream face of a dam: (a) side view showing force increasing 
linearly with depth; (b] front view showing width of dam in meters. 
24.16 Water exerts pressure on the upstream face of a dam as shown 
in Fig. P24.16. The pressure can be characterized by 
where p(z) = pressure in pascals (or N/m2) exerted at an elevation 
z meters above the reservoir bottom; p = density of water, which 
for this problem is assumed to be a constant lo3 kg/m3; g = accel- 
eration due to gravity (9.8 m/s2); and D = elevation (in m) of the 
water surface above the reservoir bottom. According to Eq. 
(P24.16). pressure increases linearly with depth, as depicted in 
Fig. P24.16~. Omitting atmospheric pressure (because it works 
against both sides of the dam face and essentially cancels out), the 
total force f can be determined by multiplying pressure times the 
area of the dam face (as shown in Fig. P24.16b). Because both 
pressure and area vary with elevation, the total force is obtained by 
evaluating 
where w(z) = width of the dam face (m) at elevation z (Fig. 
P24.16b). The line of action can also be obtained by evaluating 
Use Simpson's rule to computef, and d. Check the results with your 
computer program for the trapezoidal rule. 
24.17 To estimate the size of a new dam, you have to determine the 
total volume of water (m3) that flows down a river in a year's time. 
You have available the following long-term average data for the 
river: 
Mid- Mid Mid- Mid- Mid- Mid- Mid- Mid- Mid- 
Determine the volume. Be careful of units, and take care to make a 
proper estimate of flow at the end points. 
24.18 The data listed in the following table give hourly measure- 
ments of heat flux q (cal/cm2/h) at the surface of a solar collector. 
As an architectural engineer, you must estimate the total heat 
absorbed by a 150,000-cm2 collector panel during a 14-h period. 
The panel has an absorption efficiency e,b of 45%. The total heat 
absorbed is given by 
where A = area and q = heat flux. 
658 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
24.19 The heat flux q is the quantity of heat flowing through a unit 
area of a material per unit time. It can be computed with Fourier's 
law, 
where J has units of J/m2/s or W/m2 and k is a coefficient of ther- 
mal conductivity that parameterizes the heat-conducting properties 
of the material and has units of W/("C .m). T = temperature ("C); 
and x = distance (m) along the path of heat flow. Fourier's law is 
used routinely by architectural engineers to determine heat flow 
through walls. The following temperatures are measured into a 
stone wall: 
If the flux at x = 0 is 60 w/m2, compute k. 
Electrical Engineering 
24.20 Perform the same computation in Sec. 24.3, but for the cur- 
rent as specified by 
i(t) = 4e-'." sin 2nt for 0 5 t 5 T/2 
i(t) = 0 for T/2 < t 5 T 
where T = 1 s. Use five point Gauss Quadrature to estimate the 
integral. 
24.21 Repeat Prob. 24.20, but use a five segment Simpson's 113 
Rule. 
24.22 Repeat Prob. 24.20, but use Romberg integration to 
6, = 1%. 
24.23 Faraday's law characterizes the voltage drop across an in- 
ductor as 
where VL = voltage drop (V), L = inductance (in henrys; 1 H = 
1 V.s/A), i = current (A), and t = time (s). Determine the voltage 
drop as a function of time from the following data for an inductance 
of 4 H. 
24.24 Suppose that the current through a resistor is described by 
the function 
Compute the average voltage over t = 0 to 60 using the multiple- 
segment Simpson's 113 rule. 
Mechanical/Aerospace Engineering 
24.25 Perform the same computation as in Sec. 24.4. but use the 
following equation to compute: 
Employ the values of 8 from Table 24.6. 
24.26 Perform the same computation as in Sec. 24.4, but use the 
following equation to compute: 
Employ the equation from Prob. 24.25 for F(x). Use 4-, 8-, and 
16-segment trapezoidal rules to compute the integral. 
24.27 Repeat Prob. 24.26. but use Simpson's 113 rule. 
24.28 Repeat Prob. 24.26, but use Romberg integration to 
E, = 0.5%. 
24.29 Repeat Prob. 24.26, but use Gauss quadrature. 
24.30 The work done on an object is equal to the force times the 
distance moved in the direction of the force. The velocity of an ob- 
ject in the direction of a force is given by 
where v = m/s. Employ the multiple-application trapezoidal rule to 
determine the work if a constant force of 200 N is applied for all t. 
24.31 The rate of cooling of a body (Fig. P24.31) can be express- 
ed as 
where T = temperature of the body ("C), T, = temperature of the 
surrounding medium ("C), and k = a proportionality constant (per 
minute). Thus, this equation (called Newton's law of cooling) spec- 
ifies that the rate of cooling is proportional to the difference in the 
temperatures of the body and of the surrounding medium. If a metal 
and the resistance is a function of the current, 
R = 10i + 2i2/3 
PROBLEMS 659 
ball heated to 90°C is dropped into water that is held constant at 
T, = 25"C, the temperature of the ball changes, as in 
Time, mln 1 0 5 10 15 20 25 
T, "C / 90 49.9 33.8 28.4 26.2 25.4 
Utilize numerical differentiation to determine dT/dt at each value 
of time. Plot dT/dt versus T - T, and employ linear regression to 
evaluate k. 
24.32 A rod subject to an axial load (Fig. P24.32~) will be de- 
formed, as shown in the stress-strain curve in Fig. P24.32b. The 
area under the curve from zero stress out to the point of rupture is 
called the modulus of toughness of the material. It provides a 
measure of the energy per unit volume required to cause the mater- 
ial to rupture. As such, it is representative of the material's ability 
24.33 If the velocity distribution of a fluid flowing through a pipe 
is known (Fig. P24.33), the flow rate Q (that is, the volume of water 
passing through the pipe per unit time) can be computed by Q = 
j'u dA, where v is the velocity and A is the pipe's cross-sectional 
area. (To grasp the meaning of this relationship physically, recall 
the close connection between summation and integration.) For a 
circular pipe, A = nr2 and dA = 2nr dr. Therefore, 
where r is the radial distance measured outward from the center of 
the pipe. If the velocity distribution is given by 
-, 
to withstand an impact load. Use numerical integration to com- 
where ro is the total radius (in this case, 2 cm), compute Q using the pute the modulus of toughness for the stress-strain curve seen in 
Fig. P24.32b. multiple-application trapezoidal rule. Discuss the results. 
24.34 Using the following data, calculate the work done by stretching a 
spring that has a spring constant of k 3 x 10' N/m to x = 0.45 m. 
FIGURE 824.32 
[a) A rod under axial loading and [b) the resulting stress-strain curve where stress is in kips per 
square inch ( 1 O3 Ib/ln2) and strain is dimensionless. 
660 ENGINEERING APPLICATIONS: NUMERICAL INTEGRATION AND DIFFERENTIATION 
FIGURE P24.33 
24.35 A jet fighter's position on an aircraft carrier's runway was 
timed during landing: 
where x is the distance from the end of the carrier. Estimate 
(a) velocity (dxldt) and (b) acceleration (dv ld t ) using numerical 
differentiation. 
24.36 Employ the multiple-application trapezoidal rule to evaluate 
the vertical distance traveled by a rocket if the vertical velocity is 
given by 
24.37 The upward velocity of a rocket can be computed by the fol- 
lowing formula: 
where v =upward velocity, u = velocity at which fuel is expelled rel- 
ative to the rocket, mo = initial mass of the rocket at time t = 0, 
q = fuel consumption rate, and g = downward acceleration of grav- 
ity (assumed constant = 9.8 m/s2). If u = 2000 m/s, mo = 150,00Okg, 
and q = 2600 kg/s, use six-segment trapezoidal and Simpson's 113 
rule, six-point Gauss quadrature, and order of hx Romberg methods 
to determine how high the rocket will fly in 30 s. 
LOGUE: PART S 
PT6,4 TRADE-OFFS 
Table PT6.4 provides a summary of the trade-offs involved in numerical integrationor 
quadrature. Most of these methods are based on the simple physical interpretation of an in- 
tegral as the area under a curve. These techniques are designed to evaluate the integral of 
two different cases: (1) a mathematical function and (2) discrete data in tabular form. 
The Newton-Cotes formulas are the primary methods discussed in Chap. 21. They are 
applicable to both continuous and discrete functions. Both closed and open versions of 
these formulas are available, The open forms, which have integration limits that extend be- 
yond the range of the data, are rarely used for the evaluation of definite integrals. However, 
they have utility for the solution of ordinary differential equations and for evaluating im- 
proper integrals. 
The closed Newton-Cotes formulas are based on replacing a mathematical function or 
tabulated data by an interpolating polynomial that is easy to integrate. The simplest version 
is the trapezoidal rule, which is based on taking the area below a straight line joining adja- 
cent values of the function. One way to improve the accuracy of the trapezoidal rule is to 
divide the integration interval from a to b into a number of segments and apply the method 
to each segment. 
Aside from applying the trapezoidal rule with finer segmentation, another way to ob- 
tain a more accurate estimate of the integral is to use higher-order polynomials to connect 
the points. If a quadratic equation is employed, the result is Simpson's 1/3 rule. If a cubic 
TABLE PT6.4 Comparison of the characteristics of alternative methods For numerical integration. 
The comparisons are based on general experience and do not account for the 
behaviol of special functions. - 
Dab Points Data Points 
Required for Required For Truncation Programming 
Method One Application n Applications Error Application ENort Comments 
Trapezoidal rule 2 
Simpson's 113 rule 3 
Simpson's rule 3 o r 4 
( 1 / 3 and 3/81 
Higher-order ? 5 
Newton-Cotes 
Rornberg integrahon 3 
Gauss quadrature 1 2 
n+ l =h3i"[<) Wide 
2 n + 1 =h5fi4 ' [ t ) Wide 
2 3 h f ] Wide 
N / A =h7f161(t] Rare 
Requires f [x ] 
be known 
Requires f ix] 
be known 
Easy 
Easy 
Easy 
Easy 
Moderate Inappropriate for 
tabular data 
Easy Inappropriate for 
tabular data 
662 EPILOGUE: PART SIX 
equation is used, the result is Simpson's 318 rule. Because they are much more accurate 
than the trapezoidal rule, these formulas are usually preferred and multiple-application ver- 
sions are available. For situations with an even number of segments, the multiple applica- 
tion of the 113 rule is recommended. For an odd number of segments, the 3/8 rule can be 
applied to the last three segments and the 113 rule to the remaining segments. 
Higher-order Newton-Cotes formulas are also available. However, they are rarely used 
in practice. Where high accuracy is required, Romberg integration and Gauss quadrature 
formulas are available. It should be noted that both Romberg integration and Gauss quad- 
rature are of practical value only in cases where the function is available. These techniques 
are ill-suited for tabulated data. 
PT6.4 IMPBRTANT RELATIONSHIPS AND FORMUUS 
Table PT6.5 summarizes important information that was presented in Part Six. This table 
can be consulted to quickly access important relationships and formulas. 
PT6.6 ADVANCED METHODS AND ADDITIONAL REFERENCES 
Although we have reviewed a number of numerical integration techniques, there are other 
methods that have utility in engineering practice. For example, udaptive Simpson's inte- 
gration is based on dividing the integration interval into a series of subintervals of width h. 
Then Simpson's 113 rule is used to evaluate the integral of each subinterval by halving the 
step size in an iterative fashion, that is, with a step size of h, h/2, h /4 , h /8 , and so forth. 
The iterations are continued for each subinterval until an approximate error estimate falls 
below a prespecified stopping criterion E,. The total integral is then computed as the sum- 
mation of the integral estimates for the subintervals. This technique is especially valuable 
for complicated functions that have regions exhibiting both lower- and higher-order varia- 
tions. Discussions of adaptive integration may be found in Gerald and Wheatley (1984) and 
Rice (1983). In addition, adaptive schemes for solving ordinary differential equations can 
be used to evaluate complicated integrals, as was mentioned in PT6.1 and as will be dis- 
cussed in Chap. 25. 
Another method for obtaining integrals is to fit cubic splines to the data. The resulting 
cubic equations can be integrated easily (Forsythe et al., 1977). A similar approach is also 
sometimes used for differentiation. Finally, aside from the Gauss-Legendre formulas dis- 
cussed in Sec. 22.2, there are a variety of other quadrature formulas. Camahan, Luther, and 
Wilkes (1969) and Ralston and Rabinowitz (1978) summarize many of these approaches. 
In summary, the foregoing is intended to provide you with avenues for deeper explo- 
ration of the subject. Additionally, all the above references describe basic techniques cov- 
ered in Part Six. We urge you to consult these alternative sources to broaden your under- 
standing of numerical methods for integration. 
PT6.6 ADVANCED METHODS AND ADDITIONAL REFERENCES 663 
Method Formulation Gmphic Interpretations Error 
Trapezo~dal rule 
flai + flbi 
I - . ( b - a) ---- 
2 
Simpson's 
113 rule 
n- l 
ilxoi + 4 C fix.) + 2 5' flxll + flxol 
Multiple-application 1 -. lb - a) r = l 3 1=2 4 
Simpson's 113 rule 3n 
Simpson's 
318 r ~ l e 
Rornberg 
integration 
Gauss 
quad.ature 
a = x0 b = x,, x 
[b - a)' 
-- * 80 f'41(<1 
f i x ) 
(b - 015 f,4i -- 
1 8 0 1 ~ 
	Preface
	Contents
	Part One: Modeling, Computers and Error Analysis
	Modeling, Computers and Error Analysis
	1 Mathematical Modeling and Engineering Problem Solving
	2 Computers and Software
	3 Approximations and Round-Off Errors
	4 Truncation Errors and the Taylor Series
	Epilogue: Part One
	Part Two: Roots of Equations
	Roots of Equations
	5 Bracketing Methods
	6 Open Methods
	7 Roots of Polynomials
	8 Engineering Applications: Roots of Equations
	Epilogue: Part Two
	Part Three: Linear Algebraic Equations
	Linear Algebraic Equations
	9 Gauss Elimination
	10 LU Decomposition and Matrix Inversion
	11 Special Matrices and Gauss-Seidel
	12 Engineering Applications: Linear Algebraic Equations
	Epilogue: Part Three
	Part Four: Optimization
	Optimization
	13 One-Dimensional Unconstrained Optimization
	14 Multimensional Unconstrained Optimization
	15 Constrained Optimization
	16 Engineering Applications: Optimization
	Epilogue: Part Four
	Part Five: Curve Fitting
	Curve Fitting
	17 Least-Squares Regression
	18 Interpolation
	19 Fourier Approximation
	20 Engineering Applications: Curve Fitting
	Epilogue: Part Five
	Part Six: Numerical Differentation and Integration
	Numerical Differentation and Integration
	21 Newton-Cotes Integration Formulas
	22 Integration of Equations
	23 Numerical Differentiation
	24 Engineering Applications: Numerical Integration and Differentiation
	Epilogue: part Six
	Part Seven: Ordinary Differential Equations
	Ordinary Differential Equations
	25 Runge-Kutta Methods
	26 Stiffness and Multistep Methods
	27 Boundary-Value and Eingenvalue Problems
	28 Engineering Applications: Ordinary Differential Equations
	Epilogue: Part Seven
	Part Eight: Partial Differential Equations
	Partial Differential Equations
	29 Finite Difference: Elliptic Equations
	30 Finite Difference: Parabolic Equations
	31 Finite-Element Method
	32 Engineering Applications: Partial Differential Equations
	Epilogue: Part Eight
	Appendix A: The Fourier Series
	Appendix B: Getting Started with Mathcad
	Appendix C: Getting Started with Matlab
	Bibliography
	Index
<<
 /ASCII85EncodePages false
 /AllowTransparency false
 /AutoPositionEPSFiles true
 /AutoRotatePages /All
 /Binding /Left
 /CalGrayProfile (Dot Gain20%)
 /CalRGBProfile (sRGB IEC61966-2.1)
 /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2)
 /sRGBProfile (sRGB IEC61966-2.1)
 /CannotEmbedFontPolicy /Warning
 /CompatibilityLevel 1.4
 /CompressObjects /Tags
 /CompressPages true
 /ConvertImagesToIndexed true
 /PassThroughJPEGImages true
 /CreateJDFFile false
 /CreateJobTicket false
 /DefaultRenderingIntent /Default
 /DetectBlends true
 /DetectCurves 0.0000
 /ColorConversionStrategy /LeaveColorUnchanged
 /DoThumbnails false
 /EmbedAllFonts true
 /EmbedOpenType false
 /ParseICCProfilesInComments true
 /EmbedJobOptions true
 /DSCReportingLevel 0
 /EmitDSCWarnings false
 /EndPage -1
 /ImageMemory 1048576
 /LockDistillerParams false
 /MaxSubsetPct 100
 /Optimize true
 /OPM 1
 /ParseDSCComments true
 /ParseDSCCommentsForDocInfo true
 /PreserveCopyPage true
 /PreserveDICMYKValues true
 /PreserveEPSInfo true
 /PreserveFlatness true
 /PreserveHalftoneInfo false
 /PreserveOPIComments false
 /PreserveOverprintSettings true
 /StartPage 1
 /SubsetFonts true
 /TransferFunctionInfo /Apply
 /UCRandBGInfo /Preserve
 /UsePrologue false
 /ColorSettingsFile ()
 /AlwaysEmbed [ true
 ]
 /NeverEmbed [ true
 ]
 /AntiAliasColorImages false
 /CropColorImages true
 /ColorImageMinResolution 300
 /ColorImageMinResolutionPolicy /OK
 /DownsampleColorImages true
 /ColorImageDownsampleType /Bicubic
 /ColorImageResolution 300
 /ColorImageDepth -1
 /ColorImageMinDownsampleDepth 1
 /ColorImageDownsampleThreshold 1.50000
 /EncodeColorImages true
 /ColorImageFilter /DCTEncode
 /AutoFilterColorImages true
 /ColorImageAutoFilterStrategy /JPEG
 /ColorACSImageDict <<
 /QFactor 0.15
 /HSamples [1 1 1 1] /VSamples [1 1 1 1]
 >>
 /ColorImageDict <<
 /QFactor 0.15
 /HSamples [1 1 1 1] /VSamples [1 1 1 1]
 >>
 /JPEG2000ColorACSImageDict <<
 /TileWidth 256
 /TileHeight 256
 /Quality 30
 >>
 /JPEG2000ColorImageDict <<
 /TileWidth 256
 /TileHeight 256
 /Quality 30
 >>
 /AntiAliasGrayImages false
 /CropGrayImages true
 /GrayImageMinResolution 300
 /GrayImageMinResolutionPolicy /OK
 /DownsampleGrayImages true
 /GrayImageDownsampleType /Bicubic
 /GrayImageResolution 300
 /GrayImageDepth -1
 /GrayImageMinDownsampleDepth 2
 /GrayImageDownsampleThreshold 1.50000
 /EncodeGrayImages true
 /GrayImageFilter /DCTEncode
 /AutoFilterGrayImages true
 /GrayImageAutoFilterStrategy /JPEG
 /GrayACSImageDict <<
 /QFactor 0.15
 /HSamples [1 1 1 1] /VSamples [1 1 1 1]
 >>
 /GrayImageDict <<
 /QFactor 0.15
 /HSamples [1 1 1 1] /VSamples [1 1 1 1]
 >>
 /JPEG2000GrayACSImageDict <<
 /TileWidth 256
 /TileHeight 256
 /Quality 30
 >>
 /JPEG2000GrayImageDict <<
 /TileWidth 256
 /TileHeight 256
 /Quality 30
 >>
 /AntiAliasMonoImages false
 /CropMonoImages true
 /MonoImageMinResolution 1200
 /MonoImageMinResolutionPolicy /OK
 /DownsampleMonoImages true
 /MonoImageDownsampleType /Bicubic
 /MonoImageResolution 1200
 /MonoImageDepth -1
 /MonoImageDownsampleThreshold 1.50000
 /EncodeMonoImages true
 /MonoImageFilter /CCITTFaxEncode
 /MonoImageDict <<
 /K -1
 >>
 /AllowPSXObjects false
 /CheckCompliance [
 /None
 ]
 /PDFX1aCheck false
 /PDFX3Check false
 /PDFXCompliantPDFOnly false
 /PDFXNoTrimBoxError true
 /PDFXTrimBoxToMediaBoxOffset [
 0.00000
 0.00000
 0.00000
 0.00000
 ]
 /PDFXSetBleedBoxToMediaBox true
 /PDFXBleedBoxToTrimBoxOffset [
 0.00000
 0.00000
 0.00000
 0.00000
 ]
 /PDFXOutputIntentProfile ()
 /PDFXOutputConditionIdentifier ()
 /PDFXOutputCondition ()
 /PDFXRegistryName ()
 /PDFXTrapped /False
 /Description <<
 /CHS <FEFF4f7f75288fd94e9b8bbe5b9a521b5efa7684002000500044004600206587686353ef901a8fc7684c976262535370673a548c002000700072006f006f00660065007200208fdb884c9ad88d2891cf62535370300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c676562535f00521b5efa768400200050004400460020658768633002>
 /CHT <FEFF4f7f752890194e9b8a2d7f6e5efa7acb7684002000410064006f006200650020005000440046002065874ef653ef5728684c9762537088686a5f548c002000700072006f006f00660065007200204e0a73725f979ad854c18cea7684521753706548679c300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c4f86958b555f5df25efa7acb76840020005000440046002065874ef63002>
 /DAN <FEFF004200720075006700200069006e0064007300740069006c006c0069006e006700650072006e0065002000740069006c0020006100740020006f007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400650072002000740069006c0020006b00760061006c00690074006500740073007500640073006b007200690076006e0069006e006700200065006c006c006500720020006b006f007200720065006b007400750072006c00e60073006e0069006e0067002e0020004400650020006f007000720065007400740065006400650020005000440046002d0064006f006b0075006d0065006e0074006500720020006b0061006e002000e50062006e00650073002000690020004100630072006f00620061007400200065006c006c006500720020004100630072006f006200610074002000520065006100640065007200200035002e00300020006f00670020006e0079006500720065002e>
 /DEU <FEFF00560065007200770065006e00640065006e0020005300690065002000640069006500730065002000450069006e007300740065006c006c0075006e00670065006e0020007a0075006d002000450072007300740065006c006c0065006e00200076006f006e002000410064006f006200650020005000440046002d0044006f006b0075006d0065006e00740065006e002c00200076006f006e002000640065006e0065006e002000530069006500200068006f00630068007700650072007400690067006500200044007200750063006b006500200061007500660020004400650073006b0074006f0070002d0044007200750063006b00650072006e00200075006e0064002000500072006f006f0066002d00470065007200e400740065006e002000650072007a0065007500670065006e0020006d00f60063006800740065006e002e002000450072007300740065006c006c007400650020005000440046002d0044006f006b0075006d0065006e007400650020006b00f6006e006e0065006e0020006d006900740020004100630072006f00620061007400200075006e0064002000410064006f00620065002000520065006100640065007200200035002e00300020006f0064006500720020006800f600680065007200200067006500f600660066006e00650074002000770065007200640065006e002e>
 /ESP <FEFF005500740069006c0069006300650020006500730074006100200063006f006e0066006900670075007200610063006900f3006e0020007000610072006100200063007200650061007200200064006f00630075006d0065006e0074006f0073002000640065002000410064006f0062006500200050004400460020007000610072006100200063006f006e00730065006700750069007200200069006d0070007200650073006900f3006e002000640065002000630061006c006900640061006400200065006e00200069006d0070007200650073006f0072006100730020006400650020006500730063007200690074006f00720069006f00200079002000680065007200720061006d00690065006e00740061007300200064006500200063006f00720072006500630063006900f3006e002e002000530065002000700075006500640065006e00200061006200720069007200200064006f00630075006d0065006e0074006f00730020005000440046002000630072006500610064006f007300200063006f006e0020004100630072006f006200610074002c002000410064006f00620065002000520065006100640065007200200035002e003000200079002000760065007200730069006f006e0065007300200070006f00730074006500720069006f007200650073002e>
 /FRA <FEFF005500740069006c006900730065007a00200063006500730020006f007000740069006f006e00730020006100660069006e00200064006500200063007200e900650072002000640065007300200064006f00630075006d0065006e00740073002000410064006f00620065002000500044004600200070006f007500720020006400650073002000e90070007200650075007600650073002000650074002000640065007300200069006d007000720065007300730069006f006e00730020006400650020006800610075007400650020007100750061006c0069007400e90020007300750072002000640065007300200069006d007000720069006d0061006e0074006500730020006400650020006200750072006500610075002e0020004c0065007300200064006f00630075006d0065006e00740073002000500044004600200063007200e900e90073002000700065007500760065006e0074002000ea0074007200650020006f007500760065007200740073002000640061006e00730020004100630072006f006200610074002c002000610069006e00730069002000710075002700410064006f00620065002000520065006100640065007200200035002e0030002000650074002000760065007200730069006f006e007300200075006c007400e90072006900650075007200650073002e>/ITA <FEFF005500740069006c0069007a007a006100720065002000710075006500730074006500200069006d0070006f007300740061007a0069006f006e00690020007000650072002000630072006500610072006500200064006f00630075006d0065006e00740069002000410064006f006200650020005000440046002000700065007200200075006e00610020007300740061006d007000610020006400690020007100750061006c0069007400e00020007300750020007300740061006d00700061006e0074006900200065002000700072006f006f0066006500720020006400650073006b0074006f0070002e0020004900200064006f00630075006d0065006e007400690020005000440046002000630072006500610074006900200070006f00730073006f006e006f0020006500730073006500720065002000610070006500720074006900200063006f006e0020004100630072006f00620061007400200065002000410064006f00620065002000520065006100640065007200200035002e003000200065002000760065007200730069006f006e006900200073007500630063006500730073006900760065002e>
 /JPN <FEFF9ad854c18cea51fa529b7528002000410064006f0062006500200050004400460020658766f8306e4f5c6210306b4f7f75283057307e30593002537052376642306e753b8cea3092670059279650306b4fdd306430533068304c3067304d307e3059300230c730b930af30c830c330d730d730ea30f330bf3067306e53705237307e305f306f30d730eb30fc30d57528306b9069305730663044307e305930023053306e8a2d5b9a30674f5c62103055308c305f0020005000440046002030d530a130a430eb306f3001004100630072006f0062006100740020304a30883073002000410064006f00620065002000520065006100640065007200200035002e003000204ee5964d3067958b304f30533068304c3067304d307e30593002>
 /KOR <FEFFc7740020c124c815c7440020c0acc6a9d558c5ec0020b370c2a4d06cd0d10020d504b9b0d1300020bc0f0020ad50c815ae30c5d0c11c0020ace0d488c9c8b85c0020c778c1c4d560002000410064006f0062006500200050004400460020bb38c11cb97c0020c791c131d569b2c8b2e4002e0020c774b807ac8c0020c791c131b41c00200050004400460020bb38c11cb2940020004100630072006f0062006100740020bc0f002000410064006f00620065002000520065006100640065007200200035002e00300020c774c0c1c5d0c11c0020c5f40020c2180020c788c2b5b2c8b2e4002e>
 /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken voor kwaliteitsafdrukken op desktopprinters en proofers. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.)
 /NOR <FEFF004200720075006b00200064006900730073006500200069006e006e007300740069006c006c0069006e00670065006e0065002000740069006c002000e50020006f0070007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e00740065007200200066006f00720020007500740073006b00720069006600740020006100760020006800f800790020006b00760061006c00690074006500740020007000e500200062006f007200640073006b0072006900760065007200200065006c006c00650072002000700072006f006f006600650072002e0020005000440046002d0064006f006b0075006d0065006e00740065006e00650020006b0061006e002000e50070006e00650073002000690020004100630072006f00620061007400200065006c006c00650072002000410064006f00620065002000520065006100640065007200200035002e003000200065006c006c00650072002000730065006e006500720065002e>
 /PTB <FEFF005500740069006c0069007a006500200065007300730061007300200063006f006e00660069006700750072006100e700f50065007300200064006500200066006f0072006d00610020006100200063007200690061007200200064006f00630075006d0065006e0074006f0073002000410064006f0062006500200050004400460020007000610072006100200069006d0070007200650073007300f5006500730020006400650020007100750061006c0069006400610064006500200065006d00200069006d00700072006500730073006f0072006100730020006400650073006b0074006f00700020006500200064006900730070006f00730069007400690076006f0073002000640065002000700072006f00760061002e0020004f007300200064006f00630075006d0065006e0074006f00730020005000440046002000630072006900610064006f007300200070006f00640065006d0020007300650072002000610062006500720074006f007300200063006f006d0020006f0020004100630072006f006200610074002000650020006f002000410064006f00620065002000520065006100640065007200200035002e0030002000650020007600650072007300f50065007300200070006f00730074006500720069006f007200650073002e>
 /SUO <FEFF004b00e40079007400e40020006e00e40069007400e4002000610073006500740075006b007300690061002c0020006b0075006e0020006c0075006f0074002000410064006f0062006500200050004400460020002d0064006f006b0075006d0065006e007400740065006a00610020006c0061006100640075006b006100730074006100200074007900f6007000f60079007400e400740075006c006f0073007400750073007400610020006a00610020007600650064006f007300740075007300740061002000760061007200740065006e002e00200020004c0075006f0064007500740020005000440046002d0064006f006b0075006d0065006e00740069007400200076006f0069006400610061006e0020006100760061007400610020004100630072006f0062006100740069006c006c00610020006a0061002000410064006f00620065002000520065006100640065007200200035002e0030003a006c006c00610020006a006100200075007500640065006d006d0069006c006c0061002e>
 /SVE <FEFF0041006e007600e4006e00640020006400650020006800e4007200200069006e0073007400e4006c006c006e0069006e006700610072006e00610020006f006d002000640075002000760069006c006c00200073006b006100700061002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e00740020006600f600720020006b00760061006c00690074006500740073007500740073006b0072006900660074006500720020007000e5002000760061006e006c00690067006100200073006b0072006900760061007200650020006f006300680020006600f600720020006b006f007200720065006b007400750072002e002000200053006b006100700061006400650020005000440046002d0064006f006b0075006d0065006e00740020006b0061006e002000f600700070006e00610073002000690020004100630072006f0062006100740020006f00630068002000410064006f00620065002000520065006100640065007200200035002e00300020006f00630068002000730065006e006100720065002e>
 /ENU (Use these settings to create Adobe PDF documents for quality printing on desktop printers and proofers. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.)
 >>
 /Namespace [
 (Adobe)
 (Common)
 (1.0)
 ]
 /OtherNamespaces [
 <<
 /AsReaderSpreads false
 /CropImagesToFrames true
 /ErrorControl /WarnAndContinue
 /FlattenerIgnoreSpreadOverrides false
 /IncludeGuidesGrids false
 /IncludeNonPrinting false
 /IncludeSlug false
 /Namespace [
 (Adobe)
 (InDesign)
 (4.0)
 ]
 /OmitPlacedBitmaps false
 /OmitPlacedEPS false
 /OmitPlacedPDF false
 /SimulateOverprint /Legacy
 >>
 <<
 /AddBleedMarks false
 /AddColorBars false
 /AddCropMarks false
 /AddPageInfo false
 /AddRegMarks false
 /ConvertColors /NoConversion
 /DestinationProfileName ()
 /DestinationProfileSelector /NA
 /Downsample16BitImages true
 /FlattenerPreset <<
 /PresetSelector /MediumResolution
 >>
 /FormElements false
 /GenerateStructure true
 /IncludeBookmarks false
 /IncludeHyperlinks false
 /IncludeInteractive false
 /IncludeLayers false
 /IncludeProfiles true
 /MultimediaHandling /UseObjectSettings
 /Namespace [
 (Adobe)
 (CreativeSuite)
 (2.0)
 ]
 /PDFXOutputIntentProfileSelector /NA
 /PreserveEditing true
 /UntaggedCMYKHandling /LeaveUntagged
 /UntaggedRGBHandling /LeaveUntagged
 /UseDocumentBleed false
 >>
 ]
>> setdistillerparams
<<
 /HWResolution [2400 2400]
 /PageSize [612.000 792.000]
>> setpagedevice