Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

## Purpose

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.

## Syntax

[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)

## Description

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:
 $∫a1b1∫a2b2…∫anbnf1,f2,…,fmdxn…dx2dx1,$
where ${f}_{i}={f}_{i}\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$, for $i=1,2,\dots ,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 $‖{\mathbf{absest}}‖\le \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{absreq}},‖{\mathbf{finest}}‖×{\mathbf{relreq}}\right)$ 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.

## References

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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{a}\left({\mathbf{ndim}}\right)$ – double array
The lower limits of integration, ${a}_{i}$, for $\mathit{i}=1,2,\dots ,n$.
2:     $\mathrm{b}\left({\mathbf{ndim}}\right)$ – double array
The upper limits of integration, ${b}_{i}$, for $\mathit{i}=1,2,\dots ,n$.
3:     $\mathrm{mincls}$int64int32nag_int scalar
Must be set either to the minimum number of funsub calls to be allowed, in which case ${\mathbf{mincls}}\ge 0$ 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:     $\mathrm{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.
Constraints:
• ${\mathbf{maxcls}}\ge {\mathbf{mincls}}$;
• ${\mathbf{maxcls}}\ge r$;
• $\text{where ​}r={2}^{n}+2{n}^{2}+2n+1\text{, if ​}n<11\text{, or ​}r=1+n\left(4{n}^{2}-6n+14\right)/3\text{, if ​}n\ge 11$.
5:     $\mathrm{nfun}$int64int32nag_int scalar
$m$, the number of integrands.
Constraint: ${\mathbf{nfun}}\ge 1$.
6:     $\mathrm{funsub}$ – function handle or string containing name of m-file
funsub must evaluate the integrands ${f}_{i}$ at a given point.
[f] = funsub(ndim, z, nfun)

Input Parameters

1:     $\mathrm{ndim}$int64int32nag_int scalar
$n$, the number of dimensions of the integrals.
2:     $\mathrm{z}\left({\mathbf{ndim}}\right)$ – double array
The coordinates of the point at which the integrands must be evaluated.
3:     $\mathrm{nfun}$int64int32nag_int scalar
$m$, the number of integrands.

Output Parameters

1:     $\mathrm{f}\left({\mathbf{nfun}}\right)$ – double array
The value of the $i$th integrand at the given point.
7:     $\mathrm{absreq}$ – double scalar
The absolute accuracy required by you.
Constraint: ${\mathbf{absreq}}\ge 0.0$.
8:     $\mathrm{relreq}$ – double scalar
The relative accuracy required by you.
Constraint: ${\mathbf{relreq}}\ge 0.0$.
9:     $\mathrm{wrkstr}\left({\mathbf{lenwrk}}\right)$ – double array
If ${\mathbf{mincls}}<0$, wrkstr must be unchanged from the previous call of nag_quad_md_adapt_multi (d01ea).

### Optional Input Parameters

1:     $\mathrm{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: ${\mathbf{ndim}}\ge 1$.
2:     $\mathrm{lenwrk}$int64int32nag_int scalar
Suggested value: ${\mathbf{lenwrk}}\ge 6n+9m+\left(n+m+2\right)\left(1+p/r\right)$, 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: ${\mathbf{lenwrk}}\ge 8×{\mathbf{ndim}}+11×{\mathbf{nfun}}+3$.

### Output Parameters

1:     $\mathrm{mincls}$int64int32nag_int scalar
Gives the number of funsub calls actually used by nag_quad_md_adapt_multi (d01ea). For the continuation case (${\mathbf{mincls}}<0$ on entry) this is the number of new funsub calls on the current call to nag_quad_md_adapt_multi (d01ea).
2:     $\mathrm{wrkstr}\left({\mathbf{lenwrk}}\right)$ – double array
Contains information about the current subdivision which could be used in a continuation call.
3:     $\mathrm{finest}\left({\mathbf{nfun}}\right)$ – double array
${\mathbf{finest}}\left(\mathit{i}\right)$ specifies the best estimate obtained from the $\mathit{i}$th integral, for $\mathit{i}=1,2,\dots ,m$.
4:     $\mathrm{absest}\left({\mathbf{nfun}}\right)$ – double array
${\mathbf{absest}}\left(\mathit{i}\right)$ specifies the estimated absolute accuracy of ${\mathbf{finest}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$.
5:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{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  ${\mathbf{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  ${\mathbf{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  ${\mathbf{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.
${\mathbf{ifail}}=4$
 On entry, ${\mathbf{ndim}}<1$, or ${\mathbf{nfun}}<1$, or ${\mathbf{maxcls}}<{\mathbf{mincls}}$, or ${\mathbf{maxcls}} (see maxcls), or ${\mathbf{absreq}}<0.0$, or ${\mathbf{relreq}}<0.0$, or ${\mathbf{lenwrk}}<8×{\mathbf{ndim}}+11×{\mathbf{nfun}}+3$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

An absolute error estimate for each integrand is output in the array absest. The function exits with ${\mathbf{ifail}}={\mathbf{0}}$ if
 $maxiabsesti≤maxabsreq,relreq×maxifinesti.$

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 ${\mathbf{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 ${\mathbf{mincls}}\ge 2r$, (see specification of maxcls) and the results compared with those from the previous call.
If the function is called with ${\mathbf{mincls}}\ge 2r$, 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+{\mathrm{log}}_{10}\left({n}^{3}\right)$ decimal digits may be inaccurate when using nag_quad_md_adapt_multi (d01ea) for large values of $n$.

## Example

This example computes
 $∫01 ∫01 ∫01 ∫01 f1,f2,…,f10 dx4 dx3 dx2 dx1 ,$
where $j=1,2,\dots ,10$, ${f}_{j}=\mathrm{ln}\left({x}_{1}+2{x}_{2}+3{x}_{3}+4{x}_{4}\right)\mathrm{sin}\left(j+{x}_{1}+2{x}_{2}+3{x}_{3}+4{x}_{4}\right)$. The program is intended to show how to exploit the continuation facility provided with nag_quad_md_adapt_multi (d01ea): the function exits with ${\mathbf{ifail}}={\mathbf{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 ${\mathbf{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;
else
mulfac = 2*(ndim^3);
end

% Turn off warnings from NAG routines, catching them instead.
wstat = warning();
warning('OFF');
% 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);
else
fprintf('Success: %7d user function calls in last call\n\n', mincls);
end

% 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);
end
end
warning(wstat);

% 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);
end
for k=1:nfun
f(k) = log(sm)*sin(double(k)+sm);
end

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)
bar(data,'group');
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]);
else
plot(data,'*');
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]);
end
xlabel('Integral Index');
ylabel(ylabelString);
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

```