\chapter{Numerical Integration: The Trapezoid Rule}


Suppose we want to find an integral $\int_a^bf(x)dx$, i.e., the area under
the curve $y=f(x)$ between $x=a$ and $x=b$, and suppose that we are unable
to evaluate this integral by the antiderivative method.  The purpose of
this section is to describe an approximation method which, while not giving
us the exact value of the integral, can be used to obtain a value which is
accurate enough for practical purposes.

\section[Riemann sums]{Riemann sums}

The most basic way to approximate the area under a curve is by {\bf Riemann
sums}.  A \index{Riemann sum} Riemann sum is the total area of a number of thin rectangles that
we draw so as to fill out roughly the area under the curve.  More
precisely, we first describe the {\bf n-th left endpoint Riemann sum}.  We
choose a number $n$, and divide the interval $[a,b]$ into $n$ equal
subintervals, each of length $(b-a)/n$.  We let $\Delta x=(b-a)/n$ denote
the length of a subinterval.  If $n$ is large, then we have very small
subintervals, and in that case it turns out that the Riemann sum we later
get will be close to the true value for the area under the curve.

\smallskip
Next, we draw a vertical line up to the curve at each endpoint of the
subintervals.  The first line has length $f(a)$, the second line has length
$f(a+\Delta x)$, the third line has length $f(a+2\Delta x)$, and so on,
until we get to the last line, which has length $f(a+n\Delta x)=f(b)$.  Now
form a rectangle by drawing a little horizontal line of length $\Delta x$
from left to right at the top of each of the vertical lines (except for the
last one at $b$).

\begin{rightindent}{2.5 truein}
 Each thin rectangle has height equal to the value of $f(x)$ at the 
{\it left} endpoint of the little subinterval that forms its base; and each
rectangle has width equal to $\Delta x$.  Thus, the area of the first thin
rectangle is $\Delta x\cdot f(a)$, the area of the second rectangle is
$\Delta x\cdot f(a+\Delta x)$, the area of the third
is $\Delta x\cdot
f(a+2\Delta x)$, and so on; finally, the area of the $n$-th thin rectangle
is 
\[
     \Delta x\cdot f(a+(n-1)\Delta x)\,.
\] 
\end{rightindent}
\hfill
%-figure EA -----------------------------
\begin{psfigure}{2.5in}{1.5in}{figEA.ps}
%\putcrosshair{0,0}
\puttext{50, 37}{\Small$y = f(x)$}
\puttext{-37., -50}{\Small$a$}
\puttext{20, -50}{\Small$b$}
\puttext{-30, -40}{\Small$\Delta x$}
\puttext{-15, -40}{\Small$\Delta x$}
\puttext{-0.6, -40}{\Small$\Delta x$}
\puttext{15, -40}{\Small$\Delta x$}
\puttext{0,-60}{\Small Left endpoint Riemann sum}
\end{psfigure}
%------------------------------

\medskip

\noindent
Hence, the sum of the areas of the
$n$ rectangles is
\[
\Delta x\Bigl(f(a)+f(a+\Delta x)+f(a+2\Delta
x)+\cdots+f(a+(n-1)\Delta x)
\Bigr).
\]  
This ``left endpoint Riemann sum'' is approximately equal to
the area under the curve.  Using the summation notation, we can write:
\[
\int_a^bf(x)dx\approx\Delta x\sum_{i=0}^{n-1}f(a+i\Delta x),
\text{ where }\Delta x=\frac{b-a}{ n}.
\]
(To interpret the summation notation, let $i$ start at the value written
below the $\sum$, namely 0, and run through the numbers until you get to
the top limit of the $\sum$, namely $n-1$; plug each value of $i$ into
$f(a+i\Delta x)$ and add the resulting terms.)

\smallskip

Alternately, we could use the {\it \index{right endpoint Riemann sum} right endpoint Riemann sum}, where we
take rectangles that go to the opposite side of the curve from the
rectangles in the left endpoint Riemann sum.  (In other words, if a
rectangle in the left endpoint sum is a little 
\begin{rightindent}{2.5in}
\noindent
below the curve --- which happens when $y=f(x)$ is an increasing function
--- then the corresponding right endpoint rectangle will be a little above
the curve; and if a left endpoint rectangle is a little above the curve ---
which happens when you have a decreasing function --- then the
corresponding right endpoint rectangle will be a little below the curve.)
In the formula for the right endpoint Riemann sum, the $f$-values start at
$f(a+\Delta x)$ and go through $f(b)$:

\end{rightindent}
\hfill
%-figure EAb -----------------------------
\begin{psfigure}{2.5in}{1.5in}{figEAb.ps}
%\putcrosshair{0,0}
\puttext{50, 37}{\Small$y = f(x)$}
\puttext{-37., -50}{\Small$a$}
\puttext{20, -50}{\Small$b$}
\puttext{-30, -40}{\Small$\Delta x$}
\puttext{-15, -40}{\Small$\Delta x$}
\puttext{-0.6, -40}{\Small$\Delta x$}
\puttext{15, -40}{\Small$\Delta x$}
\puttext{0,-60}{\Small Right endpoint Riemann sum}
\end{psfigure}
%------------------------------
\medskip

\noindent
\[
\Delta x\Bigl(f(a+\Delta x)+f(a+2\Delta
x)+\cdots+f(a+(n-1)\Delta x) +f(b)\Bigr)\,.
\]
(Compare this with the formula for the left endpoint sum.)

\smallskip

The accuracy of an approximation is greater if we choose $n$ large.
Generally, neither the left endpoint nor the right endpoint Riemann sum
method is very accurate.  That is, to get a value close to the true area
under the curve, you need to use very large $n$.


\section[Trapezoid Rule]{The Trapezoid Rule}

A better method can be obtained by {\it averaging} the
left and right endpoint Riemann sums.  This is equivalent to drawing a
slanted line across the top (see the figure below at the right).


\medskip

\begin{rightindent}{2.5 truein}

The resulting sum is called the {\bf \index{trapezoid rule} trapezoid rule}.  
The trapezoid rule is more accurate than the left or right endpoint Riemann
sum.  It is easy to find a formula for the trapezoid rule by adding
the left endpoint formula and the right endpoint formula and then
dividing by 2.  That is, we take the sum

\end{rightindent}
\hfill
%-figure EB -----------------------------
\begin{psfigure}{2.5in}{1.5in}{figEB.ps}
%\putcrosshair{0,0}
\puttext{50, 37}{\Small$y = f(x)$}
\puttext{-37., -50}{\Small$a$}
\puttext{20, -50}{\Small$b$}
\puttext{-30, -40}{\Small$\Delta x$}
\puttext{-15, -40}{\Small$\Delta x$}
\puttext{-0.6, -40}{\Small$\Delta x$}
\puttext{15, -40}{\Small$\Delta x$}
\puttext{0,-60}{\Small The Trapezoidal Rule}
\end{psfigure}
%------------------------------

\medskip

\begin{align*}
\Delta x\Bigl(f(a)+&f(a+\Delta x)+f(a+2\Delta x)+\cdots+
f(a+(n-1)\Delta x)\Bigr)+\\
+\Delta x\Bigl(&f(a+\Delta x)+f(a+2\Delta x)+\cdots+f(a+(n-1)\Delta x)+
f(b)\Bigr)
\end{align*}
and divide by 2.  The result is the trapezoid rule:

\smallskip

\begin{align*} 
\int_a^bf(x)dx\approx\Delta x&\Bigl(\frac{1}{ 2}f(a)+f(a+\Delta x)+
f(a+2\Delta x)+f(a+3\Delta x)+\\
&+\cdots+ f(a+(n-1)\Delta x)+\frac{1}{ 2}f(b)\Bigr).
\end{align*} 
Notice the difference with the Riemann sum formula: both 
$f(a)$ and $f(b)$ appear at the two extremes of the sum, each with coefficient
$\frac{1}{ 2}$.

\begin{example}
%{\bf Example 1.} 
Compare the approximations to $\int_1^2 x^{-2}dx$ using
the \index{left endpoint Riemann sum} left endpoint Riemann sum, the right endpoint Riemann sum, and the
trapezoid rule --- all with $n=4$ subdivisions --- with
the true value of the integral, which is $-x^{-1}\bigm|_1^2=1/2$.

\smallskip

Here $f(x)=1/x^2$, $\Delta x=1/4$, and the 4 subintervals have endpoints
$1$, $5/4$, $3/2$, $7/4$ and $2$.  We have

\medskip

\begin{centering}
\begin{tabular}{rcl}
$\text{left Riemann sum}=$
&
$\disp\frac{1}{ 4}\left( 1+\frac{1}{ (5/4)^2}
   +\frac{1}{(3/2)^2}+\frac{1}{(7/4)^2}\right)$
&$=0.603$;\\
$\text{right Riemann sum}=$
&
$\disp\frac{1}{ 4}\left(\frac{1}{(5/4)^2}
   +\frac{1}{(3/2)^2}+\frac{1}{(7/4)^2}+\frac{1}{
2^2}\right)$
&$=0.415$;\\ 
\\
$\text{trapezoid rule}=$
&
$\disp\frac{0.603+0.415}{2}$
&$=0.509$.
\end{tabular}\\
\end{centering}

\medskip

\noindent
Note that the left (or right) Riemann sum is pretty inaccurate, while
the trapezoid rule is a lot closer, even though the number $n=4$ of
subdivisions is rather small.  We would get an even closer approximation
to the true value of the integral if we used a larger value of $n$.

\end{example}


\begin{example}
%{\bf Example 2.} 
Use the trapezoid rule with $n=6$ subdivisions to find the
length of one arch of the sine curve.

\smallskip

In Example 5 of the last section, we saw that an arch of the sine curve
$y=f(x)=\sin x$ has length equal to $\int_0^{\pi}\sqrt{1+(f'(x))^2}dx=
\int_0^{\pi}\sqrt{1+\cos^2x}\,dx$.  We now apply the trapezoid rule to this
integral.  When using the trapezoid rule, be careful about what function you
are substituting for the $f(x)$.  In this problem, because we are taking
the arclength rather than the area under the sine curve, the function to
substitute in the trapezoid rule is \underline{not} our original
function $\sin x$.  The $f(x)$ in the trapezoid rule is $\sqrt{1+\cos^2x}$.
That is, it is a {\it different} $f(x)$ from the $f(x)$ we substituted into
the arclength formula $\int_a^b\sqrt{1+(f'(x))^2}dx$.  Be careful not to
get confused: the $f(x)$ in the arclength formula is not the $f(x)$ in
the trapezoid rule.

\smallskip

Taking $n=6$, we have $\Delta x=(\pi-0)/6=\pi/6$.  Thus,
\begin{align*}
 \int_0^{\pi}&\sqrt{1+\cos^2x}\,dx\approx\frac{\pi}{ 6}\Bigl(\frac{1}{ 2}
\sqrt{1+\cos^20}+\sqrt{1+\cos^2(\pi/6)}+\sqrt{1+\cos^2(\pi/3)}\\
&+\sqrt{1+ \cos^2(\pi/2)} 
+ \sqrt{1+\cos^2(2\pi/3)}+\sqrt{1+\cos^2(5\pi/6)}+\frac{1}{ 2}
\sqrt{1+\cos^2(\pi)}\Bigr)\\
&=\frac{\pi}{ 6}\Bigl(\frac{1}{ 2}
\sqrt{2}+\sqrt{7/4}+\sqrt{5/4}+1+\sqrt{5/4}+\sqrt{7/4}+\frac{1}{ 2}\sqrt{2}
\Bigr)=3.8202.\end{align*}

Notice that the straight-line distance from the beginning to the end of the
arch is $\pi=3.1416$; so our answer of $3.8202$ for the arclength is
reasonable.  In doing arclength problems, it is a good idea to compare your
answer with a rough estimate for the length in question: if the answer
makes no sense, then go over your work again looking for careless errors.

\end{example}




\section[Interpolation]{Functions Obtained by Interpolating Experimental Data} 

There are two basic uses of numerical integration: (1) finding the
definite integral of a function given by a formula for which you cannot
find the antiderivative; and (2) finding the definite integral when the
function is not given by a formula at all. The latter situation occurs
frequently in practice. Suppose, for example, that you have a function
$f(t)$ whose values come from instrument readings at regular time
intervals from $t=a$ to $t=b$. You can graph $f(t)$ by tabulating and
plotting the points corresponding to the instrument readings, and then
drawing a smooth curve connecting these points. In this situation you
can use the table of values along with the trapezoid rule to compute
$\int_a^b f(t)dt$, or perhaps the integral of some other function made
up from $f(t)$. It is generally best to choose the number $n$ of
subdivisions to equal the number of time intervals over which you took
instrument readings.  \index{interpolation}

\begin{example}
%{\bf Example 3.} 
The table below shows the speedometer reading in
miles/hour when a car tried to brake before an accident.  Use the trapezoid
rule to determine how far the car traveled in this time.

\begin{center}
\begin{tabular}{rlllllllll}
$t$ (in seconds) &0  &1 &2 &3 &4 &5 &6 &7 &8\\
speed~(mph)      & 70&65&59&52&44&36&27&18&10
\end{tabular}\\
\end{center}

First, we note that to convert miles/hour to feet/sec we must multiply by
5280 and divide by 3600.  Thus, the distance in feet (using the trapezoid rule
with $n=8$ and $\Delta t=1$) is 
\[
\frac{5280}{ 3600}\cdot \Bigl(\frac{1}{ 2}\cdot 70+65+59+52+44+36+
27+18+\frac{1}{ 2}\cdot 10\Bigr)=500.13\>\ft.
\] 
Since your table of speeds is
given with only two significant figures of accuracy, in practical terms it
is misleading to give the distance traveled to five \index{significant figures} significant figures.
We should simply say that the car traveled 500 ft.

\end{example}

\begin{example} 
A wine cork exactly 3 cm high has circular bottom of
radius 0.542 cm and circular top of radius 0.859 cm.  The following table
shows the radius in centimeters of the cork's \index{circular cross-section} circular cross-section at
intervals of 0.5 cm from bottom to top: 

\begin{centering}
\begin{tabular}{rrrrrrr}
0.542 & 0.608 & 0.668 & 0.723 & 0.774 & 0.820 & 0.859
\end{tabular}\\
\end{centering}
\noindent
Find the volume of the cork.

\smallskip

\begin{rightindent}{3.5 truein}
The table shows the distance $g(y)$ from the \index{axis of symmetry} axis of symmetry (which we
take as the $y$-axis) to the curve that is rotated around the $y$-axis to
get the cork, taken at the endpoints of the 6 intervals of length $\Delta
y=0.5$ from $y=0$ to $y=3$.  Thus, it makes sense to apply the trapezoid
rule with $n=6$ subdivisions.  The function being integrated is $g(y)^2$:
\end{rightindent}
\hfill
%-figure ED -----------------------------
\begin{psfigure}{0.75in}{1.0in}{figED.ps}
%\putcrosshair{0,0}
%\puttext{0,-1}{}
\end{psfigure}
%------------------------------

\begin{align*} 
\text{Vol}&=\pi\int_0^3g(y)^2\,dy\\
&\approx
\pi\Delta y\Bigl(\frac{1}{ 2}g(0)^2+g(0.5)^2+g(1)^2+g(1.5)^2+g(2)^2+
g(2.5)^2+\frac{1}{ 2}g(3)^2\Bigr)\\ 
&=\frac{\pi}{ 2}\Bigl(\frac{1}{ 2}\cdot 0.542^2+ 0.608^2+
0.668^2+ 0.723^2+ 0.774^2+ 0.820^2+\frac{1}{ 2}\cdot 0.859^2\Bigr)\\
&=4.91\>\text{cm}^3.
\end{align*}

\end{example}

\begin{example}
%{\bf Example 5.}
(Measuring \index{cardiac output} cardiac output.\footnote{Taken from {\it
Measuring Cardiac Output}, by B.~Horelick and S.~Koont, Project UMAP,
Education Development Center, 1978.})
The rate $R$ at which a person's heart pumps blood is called the person's
{\bf cardiac output}.  The normal rate is about 5 liters/min, but the
cardiac output can be quite different if the person has certain illnesses
or is engaged in strenuous activity.

\smallskip

One way to measure $R$ is to inject a known amount $D$ of dye in a vein
near the heart.  The dyed blood then circulates once through the system in
a little less than half a minute.  Meanwhile, the concentration of dye in
the blood, measured in milligrams per liter, can be monitored in the aorta,
where the first of the dye will appear a few seconds after injection.  The
table below at the left shows concentration readings in mg/liter taken at
1-second intervals until all of the dye has passed the monitoring
instrument.  Here $t=0$ is the instant when $D=5.0$ mg of dye was injected
into the blood.  The graph below on the right is obtained by plotting the
points in the table and \index{interpolating} interpolating a smooth curve between those points.

\medskip

\begin{centering}
\begin{tabular}{rr||rr}
$t$    &$C(t)$   &$t$  &  $C(t)$  \\ \hline      
0      &0.0      &13  & 4.0\\     
1      &0.0      &14   &3.8\\      
2      &0.0      &15   &3.7\\      
3      &0.1      &16   &2.9\\      
4      &0.6      &17   &2.2\\      
5      &0.9      &18   &1.5\\      
6      &1.4      &19   &1.1\\      
7      &1.9      &20   &0.9\\      
8      &2.7      &21   &0.8\\      
9      &3.0      &22   &0.5\\     
10     &3.7      &23   &0.3\\      
11     &4.0      &24   &0.0\\        
12     &4.1          &&
\end{tabular} 
% NOTE:
%   The values in the table below are slightly off from
%   the values Duchamp got from the function whose graph is shown
%
\hspace{0.2 truein}
\begin{minipage}[c]{3in}
%-figure EE -----------------------------
\begin{psfigure}{3in}{2.0in}{figEE.ps}
%\putcrosshair{0,0}
\end{psfigure}
\end{minipage}\\
%------------------------------
\end{centering}

\medskip

We now derive a relationship between the cardiac output $R$, the amount of
dye $D$, and the integral under the $C(t)$ curve.  To do this, we define a
function $D(t)$ to be the amount of dye that has passed by the monitoring
point within the first $t$ seconds.  Then $dD/dt$ is the rate at which dye
is passing.  We have 
\[
\frac{dD}{ dt}=R\cdot C(t),
\] 
i.e., {\it the number of
mg per sec of dye is equal to the number of liters per second of blood
times the number of mg of dye per liter of blood.} This is an important
practical principle, which we shall encounter again when we treat mixing
problems: {\it the rate (in mass per unit time) at which an impurity (such
as dye) is flowing past a point is equal to the product of the rate (in
liters per unit time) that the liquid is flowing past and the concentration
(mass per liter).} For example, if 2 liters per second are going by, and
each liter contains 3 mg of impurity, then 6 mg per second of impurity are
going by.

\smallskip

Returning to the dye, we now have an expression for the derivative of our
function $D(t)$.  Then the total dye that passes --- which is the amount
$D$ that we originally injected --- is equal to the definite integral of
that expression, i.e., 
\[
\frac{dD}{
dt}=R\>C(t)
\text{ implies~that }
D=D(24)-D(0)=\int_0^{24}R \,C(t)\,dt=R \,\int_0^{24}C(t)\,dt,
\] 
where we took
the constant $R$ outside of the integral.  Suppose we find the integral of
$C(t)$.  If we knew $R$, then this formula would tell us what $D$ is.  As
it is, we do not know $R$, but we do know $D$.  So we want to use this
formula to solve for $R$ in terms of $D$:
\[
R=\frac{D}{\int_0^{24}C(t)\,dt}=\frac{5.0}{\int_0^{24}C(t)\,dt}.
\] 
It remains
to find the integral of $C(t)$, which we do using the trapezoid rule with
$n=24$ subdivisions (so that $\Delta t=1$ sec): 
\[
\int_0^{24}C(t)\,dt \approx
\frac{1}{ 2}f(0)+f(1)+f(2)+f(3)+
\cdots+f(23)+\frac{1}{ 2}f(24)=44.1\,.
\]  
Finally, $R=5.0/44.1=0.113$ liters/sec $=6.8$ liters/min.  Notice that
if $f(a)=0$ and $f(b)=0$, as is the case in this example, then there is
no difference between the left (or right) endpoint Riemann sum and the
trapezoid rule.

\end{example}

