hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_quad_md_adapt_multi (d01ea)


    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example


nag_quad_md_adapt_multi (d01ea) computes approximations to the integrals of a vector of similar functions, each defined over the same multidimensional hyper-rectangular region. The function uses an adaptive subdivision strategy, and also computes absolute error estimates.


[mincls, wrkstr, finest, absest, ifail] = d01ea(a, b, mincls, maxcls, nfun, funsub, absreq, relreq, wrkstr, 'ndim', ndim, 'lenwrk', lenwrk)
[mincls, wrkstr, finest, absest, ifail] = nag_quad_md_adapt_multi(a, b, mincls, maxcls, nfun, funsub, absreq, relreq, wrkstr, 'ndim', ndim, 'lenwrk', lenwrk)


nag_quad_md_adapt_multi (d01ea) uses a globally adaptive method based on the algorithm described by van Dooren and de Ridder (1976) and Genz and Malik (1980). It is implemented for integrals in the form:
where fi=fix1,x2,,xn, for i=1,2,,m.
Upon entry, unless mincls has been set to a value less than or equal to 0, nag_quad_md_adapt_multi (d01ea) divides the integration region into a number of subregions with randomly selected volumes. Inside each subregion the integrals and their errors are estimated. The initial number of subregions is chosen to be as large as possible without using more than mincls calls to funsub. The results are stored in a partially ordered list (a heap). The function then proceeds in stages. At each stage the subregion with the largest error (measured using the maximum norm) is halved along the coordinate axis where the integrands have largest absolute fourth differences. The basic rule is applied to each half of this subregion and the results are stored in the list. The results from the two halves are used to update the global integral and error estimates (finest and absest) and the function continues unless absestmaxabsreq,finest×relreq where the norm . is the maximum norm, or further subdivision would use more than maxcls calls to funsub. If at some stage there is insufficient working storage to keep the results for the next subdivision, the function switches to a less efficient mode; only if this mode of operation breaks down is insufficient storage reported.


Genz A C and Malik A A (1980) An adaptive algorithm for numerical integration over an N-dimensional rectangular region J. Comput. Appl. Math. 6 295–302
van Dooren P and de Ridder L (1976) An adaptive algorithm for numerical integration over an N-dimensional cube J. Comput. Appl. Math. 2 207–217


Compulsory Input Parameters

1:     andim – double array
The lower limits of integration, ai, for i=1,2,,n.
2:     bndim – double array
The upper limits of integration, bi, for i=1,2,,n.
3:     mincls int64int32nag_int scalar
Must be set either to the minimum number of funsub calls to be allowed, in which case mincls0 or to a negative value. In this case, the function continues the calculation started in a previous call with the same integrands and integration limits: no arguments other than mincls, maxcls, absreq, relreq or ifail must be changed between the calls.
4:     maxcls int64int32nag_int scalar
The maximum number of funsub calls to be allowed. In the continuation case this is the number of new funsub calls to be allowed.
  • maxclsmincls;
  • maxclsr;
  • where ​r=2n+2n2+2n+1, if ​n<11, or ​r=1+n4n2-6n+14/3, if ​n11.
5:     nfun int64int32nag_int scalar
m, the number of integrands.
Constraint: nfun1.
6:     funsub – function handle or string containing name of m-file
funsub must evaluate the integrands fi at a given point.
[f] = funsub(ndim, z, nfun)

Input Parameters

1:     ndim int64int32nag_int scalar
n, the number of dimensions of the integrals.
2:     zndim – double array
The coordinates of the point at which the integrands must be evaluated.
3:     nfun int64int32nag_int scalar
m, the number of integrands.

Output Parameters

1:     fnfun – double array
The value of the ith integrand at the given point.
7:     absreq – double scalar
The absolute accuracy required by you.
Constraint: absreq0.0.
8:     relreq – double scalar
The relative accuracy required by you.
Constraint: relreq0.0.
9:     wrkstrlenwrk – double array
If mincls<0, wrkstr must be unchanged from the previous call of nag_quad_md_adapt_multi (d01ea).

Optional Input Parameters

1:     ndim int64int32nag_int scalar
Default: the dimension of the arrays a, b. (An error is raised if these dimensions are not equal.)
n, the number of dimensions of the integrals.
Constraint: ndim1.
2:     lenwrk int64int32nag_int scalar
Suggested value: lenwrk6n+9m+n+m+21+p/r, where p is the value of maxcls and r is defined under maxcls. If lenwrk is significantly smaller than this, the function will not work as efficiently and may even fail.
Default: the dimension of the array wrkstr.
The dimension of the array wrkstr.
Constraint: lenwrk8×ndim+11×nfun+3.

Output Parameters

1:     mincls int64int32nag_int scalar
Gives the number of funsub calls actually used by nag_quad_md_adapt_multi (d01ea). For the continuation case (mincls<0 on entry) this is the number of new funsub calls on the current call to nag_quad_md_adapt_multi (d01ea).
2:     wrkstrlenwrk – double array
Contains information about the current subdivision which could be used in a continuation call.
3:     finestnfun – double array
finesti specifies the best estimate obtained from the ith integral, for i=1,2,,m.
4:     absestnfun – double array
absesti specifies the estimated absolute accuracy of finesti, for i=1,2,,m.
5:     ifail int64int32nag_int scalar
ifail=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W  ifail=1
maxcls was too small for nag_quad_md_adapt_multi (d01ea) to obtain the required accuracy. The arrays finest and absest respectively contain current estimates for the integrals and errors.
W  ifail=2
lenwrk is too small for the function to continue. The arrays finest and absest respectively contain current estimates for the integrals and errors.
W  ifail=3
On a continuation call, maxcls was set too small to make any progress. Increase maxcls before calling nag_quad_md_adapt_multi (d01ea) again.
On entry,ndim<1,
ormaxcls<r (see maxcls),
An unexpected error has been triggered by this routine. Please contact NAG.
Your licence key may have expired or may not have been installed correctly.
Dynamic memory allocation failed.


An absolute error estimate for each integrand is output in the array absest. The function exits with ifail=0 if

Further Comments

Usually the running time for nag_quad_md_adapt_multi (d01ea) will be dominated by the time in funsub, so the maximum time that could be used by nag_quad_md_adapt_multi (d01ea) will be proportional to maxcls multiplied by the cost of a call to funsub.
On a normal call, you should set mincls=0 on entry.
For some integrands, particularly those that are poorly behaved in a small part of the integration region, nag_quad_md_adapt_multi (d01ea) may terminate prematurely with values of absest that are significantly smaller than the actual absolute errors. This behaviour should be suspected if the returned value of mincls is small relative to the expected difficulty of the integrals. When this occurs nag_quad_md_adapt_multi (d01ea) should be called again, but with an entry value of mincls2r, (see specification of maxcls) and the results compared with those from the previous call.
If the function is called with mincls2r, the exact values of finest and absest on return will depend (within statistical limits) on the sequence of random numbers generated internally within nag_quad_md_adapt_multi (d01ea) by calls to nag_rand_dist_uniform01 (g05sa). Separate runs will produce identical answers unless the part of the program executed prior to calling nag_quad_md_adapt_multi (d01ea) also calls (directly or indirectly) functions from Chapter G05, and, in addition, the series of such calls differs between runs.
Because of moderate instability in the application of the basic integration rule, approximately the last 1+ log10n3  decimal digits may be inaccurate when using nag_quad_md_adapt_multi (d01ea) for large values of n.


This example computes
01 01 01 01 f1,f2,,f10 dx4 dx3 dx2 dx1 ,  
where j=1,2,,10, fj=lnx1+2x2+3x3+4x4sinj+x1+2x2+3x3+4x4. The program is intended to show how to exploit the continuation facility provided with nag_quad_md_adapt_multi (d01ea): the function exits with ifail=1 (printing an explanatory error message) and is re-entered with maxcls reset to a larger value. The program can be used with any values of ndim and nfun, except that the expression for r must be changed if ndim>10 (see specification of maxcls).
function d01ea_example

fprintf('d01ea example results\n\n');

% Initialize variables and arrays.
nfun = int64(10);
ndim = 4;
absreq = 0;
relreq = 1e-3;
mincls = int64(0);
ircls = 2^ndim+2*ndim*ndim+2*ndim+1;
maxcls = int64(ircls);
lenwrk = 6*ndim+9*nfun+(ndim+nfun+2)*(1+maxcls/ircls);
ncall = 1;
a = zeros(1, ndim);
b = ones(1, ndim);
wrkstr = zeros(lenwrk, 1);
mkeep = zeros(1,1);
fkeep = zeros(1,1);
akeep = zeros(1,1);

% Factor for scaling the maxcls parameter (which controls the maximum
% number of times f is to be called by d01ea).
if (ndim < 10)
    mulfac = 2^ndim;
    mulfac = 2*(ndim^3);

% Turn off warnings from NAG routines, catching them instead.
wstat = warning();
% First call: the initial parameter values are set up to demonstrate the
% continuation facility of the routine; it should exit with ifail = 1
% indicating that maxcls is too small.
% Loop until ifail is not 1 or 3 (i.e. maxcls is big enough).
ifail = int64(1);
while (ifail == 1 || ifail == 3)
    [mincls, wrkstr, finest, absest, ifail] = ...
        d01ea(a, b, mincls, maxcls, nfun, @f, absreq, relreq, wrkstr);

    % Is it big enough yet?  Output appropriate message.
    if ifail ~= 0
        fprintf('maxcls = %7d mincls = %7d  ifail = %d \n\n', maxcls, ...
                mincls, ifail);
        fprintf('Success: %7d user function calls in last call\n\n', mincls);

    % Store results for plotting.
    fkeep(1:nfun,ncall) = finest(1:nfun);
    akeep(1:nfun,ncall) = absest(1:nfun);
    mkeep(ncall,1) = mincls;

    % Prepare for next call of routine.
    mincls = int64(-1);
    maxcls = maxcls * mulfac;
    ncall = ncall + 1;
    if ifail ~= 0
        fprintf('  .. increasing maxcls to %7d\n\n', maxcls);

% Turn NAG warnings back to origianl setting.
% Plot results.
fig1 = figure;
display_plot(fkeep, mkeep, 'Integrand Number');
fig2 = figure;
display_plot(akeep, mkeep, 'Estimated Error');

function f = f(ndim, z, nfun)
% Evaluate the integrands.
f = zeros(nfun,1);
sm = 0;
for n=1:ndim
    sm = sm + double(n)*z(n);
for k=1:nfun
    f(k) = log(sm)*sin(double(k)+sm);

function display_plot(data, labels, ylabelString)
% Formatting for title and axis labels.
% Use bars to plot intregral data and points on logscale for errors.
if strncmp(ylabelString, 'Integrand', 9)
    title({'Multi-dimensional Integrals of ten Integrands with ', ...
        'Various Numbers of Function Calls'});
    set(gca, 'YTick', [-0.5 -0.3 -0.1 0.1 0.3 0.5]);
    title({'Errors in Computing Integrals of ten Integrands with ', ...
        'Various Numbers of Function Calls'});
    set(gca, 'YScale', 'log');
    set(gca, 'XLim', [0.5 10.5]);
xlabel('Integral Index');
% Add a legend.
legend([num2str(labels(1)) ' calls'], ...
    [num2str(labels(2)) ' calls'], ...
    [num2str(labels(3)) ' calls'], 'Location', 'Best')
d01ea example results

maxcls =      57 mincls =      57  ifail = 1 

  .. increasing maxcls to     912

maxcls =     912 mincls =     798  ifail = 1 

  .. increasing maxcls to   14592

Success:     912 user function calls in last call


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015