This blog post shows how to numerically integrate a one-dimensional function by using the QUAD subroutine in SAS/IML software. The name "quad" is short for quadrature, which means numerical integration.
You can use the QUAD subroutine to numerically find the definite integral of a function on a finite, semi-infinite, or infinte domain. In order to use the QUAD subroutine, you must be able to write down an expression for the function that you want to integrate. The QUAD subroutine evaluates the expression at certain locations (chosen by the algorithm) so that the integral is computed to a high level of accuracy with relatively few function evaluations.
In order to use the QUAD subroutine, you must supply at least two input arguments and one output argument. The simplest syntax for the QUAD subroutine is
call quad(result, "FcnName", limits);where the arguments are as follows:
- result is an output argument that will contain the value of the integral
- "FcnName" is the name of a SAS/IML function module that evaluates the integrand at an arbitrary point within the limits of integration
- limits defines the limits of integration
A Simple Integral
For a simple example, suppose that you want to integrate the cosine function over the interval [a, b]. The following SAS/IML module, MyFunc, evaluates the integrand at an arbitrary point:
proc iml; start MyFunc(x); return( cos(x) ); finish;
In order to integrate the cosine function over, say, the interval [0, π/2], use the following statements:
/** integrate cos(x) on [a,b] **/ a = 0; b = constant("PI") / 2; call quad(R, "MyFunc", a||b);
The value of the integral is stored in the 1x1 scalar R. For this example, you can compare the numerical result to the exact value of the integral by using the indefinite integral (which is the sine function):
/** compare with exact answer */ exactInt = sin(b) - sin(a); diff = exactInt - R; print R exactInt diff;
The numerical result is within 10-9 of the exact value. By default, SAS/IML prints values by using a BEST9. format, so the value of R is displayed as 1. To see more digits of precision, you can explicitly specify a w.d format in the PRINT statement:
How Does the QUAD Subroutine Work?
In a broad sense, there are two type of numerical integration routines: those that integrate data and those that integrate functions. The trapezoidal rule and Simpson's rule are examples of techniques that can approximate the integral when given data in the form of (x, y) pairs. These techniques are not adaptive; they use the data points that they are given. I will discuss these methods in a future blog post.
In contrast, the numerical method used by the QUAD subroutine requires a function that can be evaluated at any point in the domain of integration.
The QUAD routine implements a Romberg-type integration scheme. In its simplest form, Romberg integration is a numerical integration method that uses extrapolation of trapezoidal sums to approximate an integral over a domain. Roughly speaking, Romberg's method computes a sequence of trapezoidal approximations with decreasing step sizes, and extrapolates the result to the limiting case of zero step size. The QUAD subroutine is actually more sophisticated than this and implements additional features for dealing with singularities, functions with large derivatives, and infinite domains.
The SAS/IML documentation has further details on the QUAD subroutine and includes examples of integration over infinite and semi-infinite domains. You can also use the QUAD routine to evaluate a double (iterated) integral.