7.3 Numerical integration

When the integral of a function $\int f(x) \mathrm{d}x$ cannot be found analytically, one must resort to ``quadrature,'' numerically evaluating the integral. Some seemingly simple functions, do not have analytical integrals. For example, $f(x)=e^{-x^2}$ 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 $f(x)$ 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 $\int_a^b f(x) \mathrm{d}x$ can be interpreted as the area under the curve between the points $a$ and $b$. If $f(x)$ 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 $a$ and $b$ into multiple ``panels'' and summing the area of the panels. Each panel has a width $h$ and a height given by the value of $f(x)$ 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),

$$\lim_{h\rightarrow0//N\rightarrow\infty}\sum_{i=1}^N f(x_i)h$$

Numerically, the value of the integral can be determined to arbitrary precision by successively increasing the number of panels between the limits $a$ and $b$. The following script illustrates this approach to numerical integration for the simple case of

$$\int_0^2 x^2 \mathrm{d}x=x^3/3\vert^2_0=8/3$$

% 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 = 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 $s_\mathrm{exact}$ and $s_\mathrm{approx}$ are the analytic (exact) and approximated answers respectively, then the relative error $\varepsilon$ is expressed as $\varepsilon=\vert\frac{s_\mathrm{exact}-s_\mathrm{approx}}{s_\mathrm{exact}}\vert$ The script above plots the relative error as a function of $1/h$ where $h$ is the step size. Since $1/h$ 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 $-1$ 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, $\varepsilon\propto(1/h)$. 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 $f(x)$ 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 $h$ times the average of its sides $(f(x_i)+f(x_{i+1}))/2$. Fig. 7.8 compares the convergence rates of the trapezoidal method and the rectangular method. For $\sim2000$ 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 $-2$. In other words, the method has ``quadratic convergence'' ( $\epsilon\propto(1/h)^2$). That is, doubling the number of panels (halving the step size $h$) 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 $[a,b]$ are chosen and the function $f(x)$ is interpolated linearly between each point. (Refer again to Fig. 7.7.) For sufficiently closely spaced points, this interpolation approximates $f(x)$ 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)