- (Note that we do ''not'' use any of the original DCUHRE code by Genz, which is not under a free/open-source license.) Our code is based in part on code borrowed from the [http://mint.sbg.ac.at/HIntLib/ HIntLib numeric-integration library] by Rudolf Schürer and (for 1d integrals) from code for Gauss-Kronrod quadrature from the [http://www.gnu.org/software/gsl/ GNU Scientific Library], both of which are free software under the GNU GPL. + (Note that we do ''not'' use any of the original DCUHRE code by Genz, which is not under a free/open-source license.) Our code is based in part on code borrowed from the [http://mint.sbg.ac.at/HIntLib/ HIntLib numeric-integration library] by Rudolf Schürer and (for 1d integrals) from code for Gauss-Kronrod quadrature from the [http://www.gnu.org/software/gsl/ GNU Scientific Library], both of which are free software under the GNU GPL. (Another free-software multi-dimensional integration library, unrelated to our code here but also implementing the Genz–Malik algorithm, is [http://www.feynarts.de/cuba/ Cuba].)
I am also grateful to Dmitry Turbiner (dturbiner&nbsp;&alpha;&tau;&nbsp;alum.mit.edu), who implemented an initial prototype of the "vectorized" functionality (see below) for evaluating multiple points in a single call (as opposed to multiple functions in a single call), which facilitates parallelization of the integrand evaluation. I am also grateful to Dmitry Turbiner (dturbiner&nbsp;&alpha;&tau;&nbsp;alum.mit.edu), who implemented an initial prototype of the "vectorized" functionality (see below) for evaluating multiple points in a single call (as opposed to multiple functions in a single call), which facilitates parallelization of the integrand evaluation.
Steven G. Johnson has written a simple C subroutine for adaptive multidimensional integration (cubature) of vector-valued integrands over hypercubes, i.e. to compute integrals of the form:
(Of course, it can handle scalar integrands as the special case where is a one-dimensional vector: the dimensionalities of and are independent.) The code, which is distributed as free software under the terms of the GNU General Public License (v2 or later), is based on the algorithms described in:
(Note that we do not use any of the original DCUHRE code by Genz, which is not under a free/open-source license.) Our code is based in part on code borrowed from the HIntLib numeric-integration library by Rudolf Schürer and (for 1d integrals) from code for Gauss-Kronrod quadrature from the GNU Scientific Library, both of which are free software under the GNU GPL. (Another free-software multi-dimensional integration library, unrelated to our code here but also implementing the Genz–Malik algorithm, is Cuba.)
I am also grateful to Dmitry Turbiner (dturbiner ατ alum.mit.edu), who implemented an initial prototype of the "vectorized" functionality (see below) for evaluating multiple points in a single call (as opposed to multiple functions in a single call), which facilitates parallelization of the integrand evaluation.
a gzipped tar file. This unpacks to a directory containing a README file with instructions and a stand-alone cubature.c file that you can compile and link into your program, and a header file cubature.h that you #include. (The cubature.c code is small enough that it doesn't seem worthwhile to go to the trouble of installing it as a library.) The cubature.c also contains a little test program which is produced if you compile that file with -DTEST_INTEGRATOR (and is commented out otherwise), as described below.
B. Narasimhan wrote a GNU R interface for these routines, which can be downloaded here: http://cran.r-project.org/web/packages/cubature/index.html.
VAL and ERR are arrays of length FDIM, which upon return are the computed integral values and estimated errors, respectively. (The estimated errors are based on an embedded cubature rule of lower order; for smooth functions, this estimate is usually conservative.)
The return value of adapt_integrate is 0 on success and nonzero if there was an error (currently, only out-of-memory situations).
Here, the input is an array X of length NDIM (the point to be evaluated), the output is an array FVAL of length FDIM (the vector of function values at the point X).
The FDATA argument of F is equal to the FDATA argument passed to adapt_integrate—this can be used by the caller to pass any additional information through to F as needed (rather than using global variables, which are not re-entrant). If F does not need any additional data, you can just pass FDATA = NULL and ignore the FDATA argument to F.
These integration algorithms evaluate the integrand in "batches" of several points at a time. It is often useful to have access to this information so that your integrand function is not called for one point at a time, but rather for a whole "vector" of many points at once. For example, you may want to evaluate the integrand in parallel at different points. This functionality is available by calling:
The size of NPTS will vary with the dimensionality of the problem; higher-dimensional problems will have (exponentially) larger NPTS, allowing for the possibility of more parallelism. (However, for sufficiently high dimensions, e.g. >8, a different algorithm like Monte-Carlo integration or sparse grids is likely to be more efficient.) Currently, NPTS is 30 in 1d, 34 in 2d, and 66 in 3d, although the first call to F will have half that many points.
As a simple example, consider the Gaussian integral of the scalar function over the hypercube [ − 2,2]3 in 3 dimensions. You could compute this integral via code that looks like:
Note that the estimated relative error is 0.00136919/13.69609043 = 9.9969×10–5, within our requested tolerance of 10 − 4. The actual error in the integral value, as can be determined e.g. by running the integration with a much lower tolerance, is much smaller: the integral is too small by about 0.00002, for an actual relative error of about 1.4×10–6. As mentioned above, for smooth integrands the estimated error is almost always conservative (which means, unfortunately, that the integrator usually does more function evaluations than it needs to).
Integrals over infinite or semi-infinite intervals is possible by a change of variables. This is best illustrated in one dimension.
In multiple dimensions, one simply performs this change of variables on each dimension separately, as desired, multiplying the integrand by the corresponding Jacobian factor for each dimension being transformed.
The Jacobian factors diverge as the endpoints are approached. However, in order for your integral to have been finite in the first place it must be true that f(x) goes to zero at infinity faster than 1/x, hence the limit of the integrand multiplied by the Jacobian is zero at the endpoints. In any case, the quadrature/cubature rules currently employed in cubature.c do not evaluate the integrand at the endpoints, so you need not implement special handling for |t|=1.
7: a test integrand by Morokoff and Caflisch, a simple product of dim-th roots of the coordinates (weakly singular at the boundary)
Notice that it finds the integral after 82203 function evaluations with an estimated error of about 10–5, but the true error (compared to the exact result) is much smaller (2.5×10–8): the error estimation is typically conservative when applied to smooth functions like this.