Next: 7.4 Numerical derivatives Up: 7. Workhorses of Computational Previous: 7.2.2 Newton's Method and   Contents

7.3 Numerical integration

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.
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
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 = 1e-4;   % 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 = (b-a)/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((ans-sum)/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 trapezoid-shaped 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.
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.
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 piecewise-linear 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 2005-01-28