GAMING

The D01 chapter - Quadrature

The numerical evaluation of definite integrals in one or more dimensions is a commonly encountered task in analysis. Routines from the D01 chapter provide a variety of algorithms for use in this area; the applicability of individual routines depends on the form of the integrand. If its functional form is known analytically, then d01aj is most generally suited, and can be used when the integrand contains singularities, especially of an algebraic or logarithmic type. Other, more specialized, routines include d01ak if the integrand is oscillatory, and d01al if it exhibits discontinuities at known points.

The algorithm implemented by d01aj is adaptive - that is, it divides the interval over which the function is to be integrated into a set of subintervals, which are in turn subdivided until some accuracy condition (specified by the variables epsabs and epsrel in the fragment below) is met.

Here is the code which calls the routine to calculate the integral.

epsabs = 0.0;
epsrel = 0.000001;

% lw & liw are the dimensions of the w and iw output arrays.  lw 
% should be between 800 & 2000; liw should be a quarter of that.
lw = int32(1000);
liw = lw/4;
iw = zeros(liw, 'int32');
w = zeros(lw, 'double');

% Call the NAG routine, and check for errors.
[result, abserr, w, iw, ifail] = d01aj(func, a, b, epsabs, epsrel);
if ifail ~= 0
    disp([' Non-zero ifail after d01aj call ', num2str(ifail)]);
    return;
end;

In this fragment, the variable func, passed to d01aj, is a string containing the name of a MATLAB m-file; this file must contain a function which returns the value of the integrand at a given point. Here are the contents of an example file of this type:

function [result] = myfunc(x)
    result = x^2*abs(sin(5*x));

(this is the integrand used in figure 2, below). In addition to the approximation of the integral value (result), d01aj's output contains the specification of the the final set of subintervals, along with the error associated with each one (some manipulation is required to get these arrays into a form that MATLAB's plotting routines can use). Figure 2 shows results for the contributions from each subinterval to the integral, and the associated errors. It can be seen that, for this integrand, the algorithm has used narrower subintervals in regions where the function is changing rapidly; the subintervals associated with (relatively) large errors are those whose width is very small indeed (near where the function goes to zero).

Figure 2: Calculating an approximation to the integral of a function.

The MATLAB script for this demo is available as the file NAGToolboxDemos/Quadrature/d01aj_demo.m, distributed in this archive.