Next: 7.4 Numerical derivatives
Up: 7. Workhorses of Computational
Previous: 7.2.2 Newton's Method and
Contents
When the integral of a function
cannot be found
analytically, one must resort to ``quadrature,'' numerically evaluating the
integral. Some seemingly simple functions, do not have analytical
integrals. For example, occurs frequently in physics
problems but cannot be integrated analytically. Recall that the
geometric interpretation of a definite integral is the area under the curve
between the limits of the integral.
Figure 7.6:
Illustration of a numerical approximation to integrating a
function over a definite interval.

For example, consider the example illustrated in Fig. 7.6.
The integral
can be interpreted as the area
under the curve between the points and . If cannot be
integrated analytically, then the integral must be evaluated numerically.
Conceptually, the area under the curve can be approximated by dividing the
interval between and into multiple ``panels'' and summing the area
of the panels. Each panel has a width and a height given by the value
of at the left edge of the panel. You learned in your calculus class
that the result of the definite integral is determined by taking the limit
of the sum as the width of each panel goes to zero (and thus the number of
panels goes to infinity),
Numerically, the value of the integral can be determined to arbitrary
precision by successively increasing the number of panels between the
limits and . The following script illustrates this approach to
numerical integration for the simple case of
% Simple numerical integration example
clc; clear all; close all
f=inline('x^2','x'); % Define f(x) as an inline function
fint = inline('x^3/3','x'); % Define the analytic integral of f(x)
a = 0; b = 2; % Define the limits of the integral
N = 2; % Number of panels to start with
step = 0; % Keep track of the number of different h sizes used
relerr = 1; % Initialize relative error so the while loop doesn't exit immediately
eps = 1e4; % Specify the desired accuracy (maximum relative error)
ans = fint(b)fint(a); % Compute the analytic answer so we can calculate the error
fprintf('step N width answer relative error\n')
while relerr>eps % Continue until we've reached the desired accuracy
step = step + 1; % Update number of h sizes
x = linspace(a,b,N+1); % Divide the interval into N panels (N+1 points on the interval)
sum = 0; % Initialize the sum
h = (ba)/N; % Calculate the width of panels
for i = 1:N % Loop over each panel
sum = sum + f(x(i))*h; % Calculate area of i_th panel
end
relerr(step) = abs((anssum)/ans); % Calculate relative error
fprintf('%4d %6d %9.5f %11.7f %12.3d\n', step, N, h, sum, relerr(step))
stepsize(step) = h; %Store the step sizes to be plotted
N = 2*N; % Double the number of panels and go again
if N>5e4; break; end % Exit the while loop if we're starting to use too many panels
end
loglog(1./stepsize,relerr,'*') % Note use of ./ instead of /
xlabel('step size^{1} (proportional to the number of panels)')
ylabel('Relative error')
title('Convergence speed of the rectangular integration method')
Note that we
have introduced the idea of a relative error. The relative error is
the error as a fraction of the actual answer. If
and
are the analytic (exact) and approximated answers
respectively, then the relative error is expressed as
The script above plots the relative error as a function of where
is the step size. Since is proportional the number of panels, the
graph indicates how rapidly the relative error decreases as the number of
panels is increased. The slope of the resulting graph is approximately
which indicates that halving the step size (doubling the number of panels)
results in a halving of the relative error. This is called ``linear
convergence'' because the error is linearly proportional to the number of
panels,
.
While this approach to numerical integration is intuitive and follows the
definition of definite integrals that you learned in your calculus class,
it is inefficient. It's convergence is poor. A very simple improvement is
illustrated in Fig. 7.7. In this case, instead of
rectangles, we use trapezoidshaped panels. This method does not
require any more computer time because we were already calculating
for each panel but obviously is far more accurate than the rectangular
method.
Figure 7.7:
Illustration of the trapezoidal method of quadrature (numerical integration).

The area of each panel is merely its width times the average of its
sides
. Fig. 7.8 compares the
convergence rates of the trapezoidal method and the rectangular method. For
panels, the error in the trapezoidal method is nearly one
million times smaller.
Figure 7.8:
Comparison of the convergence rates of the rectangular
integration method versus the trapezoidal method.

The slope of the graph for the trapezoidal method is about . In other
words, the method has ``quadratic convergence''
(
). That is, doubling the number of panels (halving
the step size ) reduces the relative error by a factor of four.
We can think of the trapezoidal method as a kind of interpolation scheme. A
number of evenly spaced points on the interval are chosen and the
function is interpolated linearly between each point. (Refer
again to Fig. 7.7.) For sufficiently closely spaced points,
this interpolation approximates with fair accuracy. This resulting
piecewiselinear function is then integrated exactly. Since the
trapezoidal method worked so much better than the rectangular method, you
might speculate that a method that interpolated across each panel using
parabolas instead of just straight lines would yield a method even
better than the trapezoidal method. In fact, you are correct and the method
is called Simpson's method.
Simpson's method should always be your first choice for integrating
numerically and often you won't need anything fancier. (discuss method in detail)
Next: 7.4 Numerical derivatives
Up: 7. Workhorses of Computational
Previous: 7.2.2 Newton's Method and
Contents
Gus Hart
20050128