492 nuMerIcAL InTeGrATIOn FOrMuLAs

Integrals are also employed by engineers and scientists 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 × 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 ΔVi:
n
Mass = ∑ ci ΔVi
i=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) dx dy dz

or

Mass = ∫∫ ∫ c(V ) d V
V

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 and science. For ex-
ample, 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

Flux = ∫∫ flux d A
A
which is referred to as an areal integral, where A = area.
These are just a few of the applications of integration that you might face regularly
in the pursuit of your profession. When the functions to be analyzed are simple, you will
normally choose to evaluate them analytically. However, it is often difficult or impossible
when the function is complicated, as is typically the case in more realistic examples. In
addition, the underlying function is often unknown and defined only by measurement at
discrete points. For both these cases, you must have the ability to obtain approximate val-
ues for integrals using numerical techniques as described next.

19.2 NEWTON-COTES FORMULAS

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 a poly-
nomial that is easy to integrate:
b b
I= ∫a f (x) dx ≅ ∫a fn (x) dx (19.8)
19.2 newTOn-cOTes FOrMuLAs 493

where fn(x) = a polynomial of the form

fn(x) = a0 + a1x + · · · + an−1xn−1 + an xn (19.9)
where n is the order of the polynomial. For example, in Fig. 19.4a, a first-order polynomial
(a straight line) is used as an approximation. In Fig. 19.4b, 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. 19.5, three

FIGURE 19.4
The approximation of an integral by the area under (a) a straight line and (b) a parabola.

f (x) f (x)

a b x a b x
(a) (b)

FIGURE 19.5
The approximation of an integral by the area under three straight-line segments.

f (x)

a b x
494 nuMerIcAL InTeGrATIOn FOrMuLAs

f (x) f (x)

a b x a b x
(a) (b)

FIGURE 19.6
The difference between (a) closed and (b) open integration formulas.

straight-line segments are used to approximate the integral. Higher-order polynomials can
be utilized for the same purpose.
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. 19.6a). The open forms have integration limits that extend beyond the range
of the data (Fig. 19.6b). This chapter emphasizes the closed forms. However, material on
open Newton-Cotes formulas is briefly introduced in Sec. 19.7.

19.3 THE TRAPEZOIDAL RULE

The trapezoidal rule is the first of the Newton-Cotes closed integration formulas. It cor-
responds to the case where the polynomial in Eq. (19.8) is first-order:

b
f (b) − f (a)
[ f (a) + (x − a)] dx
I= ∫a _________ (19.10)
b−a

The result of the integration is

f (a) + f (b)
I = (b − a) _________ (19.11)
2
which is called the trapezoidal rule.
Geometrically, the trapezoidal rule is equivalent to approximating the area of the trap-
ezoid under the straight line connecting f (a) and f (b) in Fig. 19.7. Recall from geometry
that the formula for computing the area of a trapezoid is the height times the average of the
bases. In our case, the concept is the same but the trapezoid is on its side. Therefore, the
integral estimate can be represented as

I = width × average height (19.12)

 19.3 The TrApeZOIdAL ruLe 495

f (x)

f (b)

f (a)

a b x

FIGURE 19.7
Graphical depiction of the trapezoidal rule.

or

I = (b − a) × average height (19.13)

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. (19.13). That is, they differ only with respect to the formulation of the average height.

19.3.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. 19.8). An esti-
mate for the local truncation error of a single application of the trapezoidal rule is
1 f ″ (ξ )(b − a)3
Et = −___ (19.14)
12
where ξ lies somewhere in the interval from a to b. Equation (19.14) indicates that if the
function being integrated is linear, the trapezoidal rule will be exact because the second
derivative of a straight line is zero. Otherwise, for functions with second- and higher-order
derivatives (i.e., with curvature), some error can occur.

eXAMpLe 19.1 single Application of the Trapezoidal rule

problem statement. Use Eq. (19.11) to numerically integrate

f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5

from a = 0 to b = 0.8. Note that the exact value of the integral can be determined analytically
to be 1.640533.
496 nuMerIcAL InTeGrATIOn FOrMuLAs

f (x)

2.0

Error

Integral estimate
0 0.8 x

FIGURE 19.8
Graphical depiction of the use of a single application of the trapezoidal rule to approximate
the integral of f (x) = 0.2 + 25x − 200x 2 + 675x 3 − 900x 4 + 400x5 from x = 0 to 0.8.

solution. The function values f (0) = 0.2 and f (0.8) = 0.232 can be substituted into
Eq. (19.11) to yield

0.2 + 0.232 = 0.1728

I = (0.8 − 0) __________
2
which represents an error of Et = 1.640533 − 0.1728 = 1.467733, which corresponds to
a percent relative error of εt = 89.5%. The reason for this large error is evident from the
graphical depiction in Fig. 19.8. Notice that the area under the straight line neglects a sig-
nificant 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
derivative over the interval can be computed by differentiating the original function twice
to give

f ″(x) = −400 + 4,050x − 10,800x 2 + 8,000x 3

The average value of the second derivative can be computed as [Eq. (19.7)]
0.8
_ ∫0 (−400 + 4,050x − 10,800x2 + 8,000x3) dx
f ″(x) = ______________________________________ = −60
0.8 − 0
which can be substituted into Eq. (19.14) to yield
1 (−60)(0.8)3 = 2.56
Ea = −___
12
19.3 The TrApeZOIdAL ruLe 497

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 is not necessarily an accurate approximation of f ″(ξ). Thus, we denote that the
error is approximate by using the notation Ea, rather than exact by using Et.

19.3.2 The Composite Trapezoidal Rule

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 (Fig. 19.9).
The areas of individual segments can then be added to yield the integral for the entire
interval. The resulting equations are called composite, or multiple-segment, integration
formulas.
Figure 19.9 shows the general format and nomenclature we will use to characterize
composite integrals. There are n + 1 equally spaced base points (x0, x1, x2, . . . , xn ). Conse-
quently, there are n segments of equal width:

h = _____
b−a
n (19.15)

If a and b are designated as x0, and xn, respectively, the total integral can be repre-
sented as
x1 x2 xn
I= ∫x
0
f (x) dx + ∫x
1
f (x) dx + · · · + ∫x
n−1
f (x) dx

FIGURE 19.9
Composite trapezoidal rule.

f(x)

x0 x1 x2 x3 x4 x5
x0 = a h= b− a xn = b
n
498 nuMerIcAL InTeGrATIOn FOrMuLAs

Substituting the trapezoidal rule for each integral yields

f (x0) + f (x1) f (x1) + f (x2) f (xn−1) + f (xn)

I = h ___________ + h __________ + · · · + h ____________ (19.16)
2 2 2
or, grouping terms:

2[ 0 n ]
n−1
I = __
h f (x ) + 2 ∑ f (x ) + f (x )
i (19.17)
i=1

or, using Eq. (19.15) to express Eq. (19.17) in the general form of Eq. (19.13):
n−1
f (x0) + 2 ∑ f (xi ) + f (xn )
I = (b − a) _____________________
i=1
2n
⏟
Width  (19.18)
Average height

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. (19.18), the interior points are given twice the weight of the two end points f (x 0)
and f (xn ).
An error for the composite trapezoidal rule can be obtained by summing the individual
errors for each segment to give

(b − a)3 n
Et = − _______ ∑ f ″(ξi ) (19.19)
12n3 i=1

where f ″(ξi ) is the second derivative at a point ξi 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
n
_
∑ f ″(ξi )
f″ ≅ ________
i=1
(19.20)
n
_
Therefore ∑ f ″(ξi ) ≅ n f ″ and Eq. (19.19) can be rewritten as
(b − a)3 _
Ea = − _______ f″ (19.21)
12n2
Thus, if the number of segments is doubled, the truncation error will be quartered. Note
that Eq. (19.21) is an approximate error because of the approximate nature of Eq. (19.20).

eXAMpLe 19.2 composite Application of the Trapezoidal rule

problem statement. Use the two-segment trapezoidal rule to estimate the integral of

f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5

from a = 0 to b = 0.8. Employ Eq. (19.21) to estimate the error. Recall that the exact value
of the integral is 1.640533.
19.3 The TrApeZOIdAL ruLe 499

solution. For n = 2 (h = 0.4):

f (0) = 0.2 f (0.4) = 2.456 f (0.8) = 0.232
0.2 + 2(2.456) + 0.232
I = 0.8 ______________________ = 1.0688
4
Et = 1.640533 − 1.0688 = 0.57173 εt = 34.9%

0.83 (−60) = 0.64

Ea = − ______
12(2)2
where −60 is the average second derivative determined previously in Example 19.1.

The results of the previous example, along with three- through ten-segment
applications of the trapezoidal rule, are summarized in Table 19.1. Notice how the
error decreases as the number of segments increases. However, also notice that the rate
of decrease is gradual. This is because the error is inversely related to the square of n
[Eq. (19.21)]. Therefore, doubling the number of segments quarters the error. In subse-
quent sections we develop higher-order formulas that are more accurate and that con-
verge more quickly on the true integral as the segments are increased. However, before
investigating these formulas, we will first discuss how MATLAB can be used to imple-
ment the trapezoidal rule.

19.3.3 MATLAB M-file: trap

A simple algorithm to implement the composite trapezoidal rule can be written as in
Fig. 19.10. The function to be integrated is passed into the M-file along with the limits of
integration and the number of segments. A loop is then employed to generate the integral
following Eq. (19.18).

TABLE 19.1 results for the composite trapezoidal rule to

estimate the integral of f (x) = 0.2 + 25x −
200x2 + 675x3 − 900x4 + 400x5 from
x = 0 to 0.8. The exact value is 1.640533.

n h I εt (%)

2 0.4 1.0688 34.9

3 0.2667 1.3695 16.5
4 0.2 1.4848 9.5
5 0.16 1.5399 6.1
6 0.1333 1.5703 4.3
7 0.1143 1.5887 3.2
8 0.1 1.6008 2.4
9 0.0889 1.6091 1.9
10 0.08 1.6150 1.6
500 nuMerIcAL InTeGrATIOn FOrMuLAs

function I = trap(func,a,b,n,varargin)
% trap: composite trapezoidal rule quadrature
% I = trap(func,a,b,n,pl,p2,...):
% composite trapezoidal rule
% input:
% func = name of function to be integrated
% a, b = integration limits
% n = number of segments (default = 100)
% pl,p2,... = additional parameters used by func
% output:
% I = integral estimate

if nargin<3,error('at least 3 input arguments required'),end

if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4|isempty(n),n = 100;end
x = a; h = (b − a)/n;
s=func(a,varargin{:});
for i = 1 : n − 1
x = x + h;
s = s + 2*func(x,varargin{:});
end
s = s + func(b,varargin{:});
I = (b − a) * s/(2*n);

FIGURE 19.10
M-file to implement the composite trapezoidal rule.

An application of the M-file can be developed to determine the distance fallen by

the free-falling bungee jumper in the first 3 s by evaluating the integral of Eq. (19.3). For
this example, assume the following parameter values: g = 9.81 m/s2, m = 68.1 kg, and
cd = 0.25 kg/m. Note that the exact value of the integral can be computed with Eq. (19.4)
as 41.94805.
The function to be integrated can be developed as an M-file or with an anonymous
function,

>> v=@(t) sqrt(9.81*68.1/0.25)*tanh(sqrt(9.81*0.25/68.1)*t)

v =
@(t) sqrt(9.81*68.1/0.25)*tanh(sqrt(9.81*0.25/68.1)*t)

First, let’s evaluate the integral with a crude five-segment approximation:

>> format long

>> trap(v,0,3,5)

ans =
41.86992959072735