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:
 b1 b2 bn ∫ ∫ … ∫ (f1,f2, … ,fm)dxn … dx2dx1, a1 a2 an
$∫a1b1∫a2b2…∫anbn(f1,f2,…,fm)dxn…dx2dx1,$
where fi = fi(x1,x2,,xn)${f}_{i}={f}_{i}\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$, for i = 1,2,,m$i=1,2,\dots ,m$.
Upon entry, unless mincls has been set to a value less than or equal to 0$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 absestmax (absreq,finest × relreq)$‖{\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:     a(ndim) – double array
ndim, the dimension of the array, must satisfy the constraint ndim1${\mathbf{ndim}}\ge 1$.
The lower limits of integration, ai${a}_{i}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
2:     b(ndim) – double array
ndim, the dimension of the array, must satisfy the constraint ndim1${\mathbf{ndim}}\ge 1$.
The upper limits of integration, bi${b}_{i}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
3:     mincls – int64int32nag_int scalar
Must be set either to the minimum number of funsub calls to be allowed, in which case mincls0${\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 parameters 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.
Constraints:
• ${\mathbf{maxcls}}\ge {\mathbf{mincls}}$;
• maxclsr${\mathbf{maxcls}}\ge r$;
• where ​r = 2n + 2n2 + 2n + 1, if ​n < 11, or ​r = 1 + n(4n26n + 14) / 3, if ​n11$\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:     nfun – int64int32nag_int scalar
m$m$, the number of integrands.
Constraint: nfun1${\mathbf{nfun}}\ge 1$.
6:     funsub – function handle or string containing name of m-file
funsub must evaluate the integrands fi${f}_{i}$ at a given point.
[f] = funsub(ndim, z, nfun)

Input Parameters

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

Output Parameters

1:     f(nfun) – double array
The value of the i$i$th integrand at the given point.
7:     absreq – double scalar
The absolute accuracy required by you.
Constraint: absreq0.0${\mathbf{absreq}}\ge 0.0$.
8:     relreq – double scalar
The relative accuracy required by you.
Constraint: relreq0.0${\mathbf{relreq}}\ge 0.0$.
9:     wrkstr(lenwrk) – double array
lenwrk, the dimension of the array, must satisfy the constraint lenwrk8 × ndim + 11 × nfun + 3${\mathbf{lenwrk}}\ge 8×{\mathbf{ndim}}+11×{\mathbf{nfun}}+3$.
If mincls < 0${\mathbf{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$n$, the number of dimensions of the integrals.
Constraint: ndim1${\mathbf{ndim}}\ge 1$.
2:     lenwrk – int64int32nag_int scalar
Default: The dimension of the array wrkstr.
The dimension of the array wrkstr as declared in the (sub)program from which nag_quad_md_adapt_multi (d01ea) is called.
Suggested value: lenwrk6n + 9m + (n + m + 2)(1 + p / r)${\mathbf{lenwrk}}\ge 6n+9m+\left(n+m+2\right)\left(1+p/r\right)$, where p$p$ is the value of maxcls and r$r$ is defined under maxcls. If lenwrk is significantly smaller than this, the function will not work as efficiently and may even fail.
Constraint: lenwrk8 × ndim + 11 × nfun + 3${\mathbf{lenwrk}}\ge 8×{\mathbf{ndim}}+11×{\mathbf{nfun}}+3$.

None.

### 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${\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:     wrkstr(lenwrk) – double array
Contains information about the current subdivision which could be used in a continuation call.
3:     finest(nfun) – double array
finest(i)${\mathbf{finest}}\left(\mathit{i}\right)$ specifies the best estimate obtained from the i$\mathit{i}$th integral, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$.
4:     absest(nfun) – double array
absest(i)${\mathbf{absest}}\left(\mathit{i}\right)$ specifies the estimated absolute accuracy of finest(i)${\mathbf{finest}}\left(\mathit{i}\right)$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$.
5:     ifail – int64int32nag_int scalar
${\mathrm{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 ifail = 1${\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 ifail = 2${\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 ifail = 3${\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.
ifail = 4${\mathbf{ifail}}=4$
 On entry, ndim < 1${\mathbf{ndim}}<1$, or nfun < 1${\mathbf{nfun}}<1$, or ${\mathbf{maxcls}}<{\mathbf{mincls}}$, or maxcls < r${\mathbf{maxcls}} (see maxcls), or absreq < 0.0${\mathbf{absreq}}<0.0$, or relreq < 0.0${\mathbf{relreq}}<0.0$, or lenwrk < 8 × ndim + 11 × nfun + 3${\mathbf{lenwrk}}<8×{\mathbf{ndim}}+11×{\mathbf{nfun}}+3$.

## Accuracy

An absolute error estimate for each integrand is output in the array absest. The function exits with ${\mathbf{ifail}}={\mathbf{0}}$ if
 max (absest(i)) ≤ max (absreq,relreq × maxi |finest(i)|). i
$maxi(absesti)≤max(absreq,relreq×maxi|finesti|).$

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${\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 mincls2r${\mathbf{mincls}}\ge 2r$, (see specification of maxcls) and the results compared with those from the previous call.
If the function is called with mincls2r${\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 + log10(n3) $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$n$.

## Example

This example computes
 1 1 1 1 ∫ ∫ ∫ ∫ (f1,f2, … ,f10)dx4dx3dx2dx1, 0 0 0 0
$∫01 ∫01 ∫01 ∫01 (f1,f2,…,f10) dx4 dx3 dx2 dx1 ,$
where j = 1,2,,10$j=1,2,\dots ,10$, fj = ln(x1 + 2x2 + 3x3 + 4x4)sin(j + x1 + 2x2 + 3x3 + 4x4)${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$r$ must be changed if ndim > 10${\mathbf{ndim}}>10$ (see specification of maxcls).
```function nag_quad_md_adapt_multi_example
% Initialize variables and arrays.
nfun = 10;
ndim = 4;
absreq = 0;
relreq = 1e-3;
mincls = 0;
ircls = 2^ndim+2*ndim*ndim+2*ndim+1;
maxcls = ircls;
lenwrk = 6*ndim+9*nfun+(ndim+nfun+2)*(1+maxcls/ircls);
w = 1:nfun;
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
if (ndim < 10)
mulfac = 2^ndim;
else
mulfac = 2*(ndim^3);
end

% Turn off warnings from NAG routines, because we'll handle these ourselves.
% Then call the routine for the first time.  The initial parameter values
% are set up to demonstrate the continuation facility of the routine; it
% should exit with ifail = 1 (maxcls is too small).
warning('off', 'NAG:warning');
int64(maxcls), int64(nfun), @f, absreq, relreq, wrkstr);

% Check exit condition.  First, deal with problems associated with
% input parameters on this initial call.
if (ifail == 2)
% Length of the work array is too small.  Print message and exit.
error('Warning: nag_quad_md_adapt_multi returned with ifail = %1d (Lenwrk is too small to obtain required accuracy).',ifail);
elseif (ifail == 4)
% Incorrect parameter values.  Print message and exit.
end

% Continue loop until ifail isn't 1 or 3 (i.e. maxcls is big enough).
while (ifail == 1 || ifail == 3)
[mincls, wrkstr, finest, absest, ifail] = ...
@f, absreq, relreq, wrkstr);

% Is it big enough yet?  Output appropriate message.
if ifail ~= 0
fprintf('** Maxcls too small to obtain required accuracy\n');
fprintf('Maxcls = %7.0f    ifail = %1.0f \n\n',maxcls, ifail);
fprintf(['Results so far (%1.0f user function calls in last ', ...
else
fprintf(['\nFinal results (%1.0f user function calls in last ', ...
end

% Output results, and store them for plotting.
disp('  I     Integral   Estimated error')
for i = 1:nfun
fprintf(' %2.0f  ',w(i));
fprintf('%10.4f   %10.4f  \n',finest(i),absest(i))
fkeep(i,ncall) = finest(i);
akeep(i,ncall) = absest(i);
end

% Also, store the number of function calls.
mkeep(ncall,1) = mincls;

% Prepare for next call of routine.
mincls = -1;
maxcls = maxcls * mulfac;
ncall = ncall + 1;
if ifail ~= 0
fprintf('\n');
end
end
% Turn NAG warnings back on.
warning('on', 'NAG:warning');
% Plot results.
fig = figure('Number', 'off');
display_plot(fkeep, mkeep, 'Integrand Number');

fig = figure('Number', 'off');
display_plot(akeep, mkeep, 'Estimated Error');

function f = f(ndim, z, nfun)
% Evaluate the integrands.
f = zeros(nfun,1);
sm = 0;
for n=1:int32(ndim) % can't use int64 in loop range
sm = sm + double(n)*z(n);
end
for k=1:int32(nfun) % can't use int64 in loop range
f(k) = log(sm)*sin(double(k)+sm);
end
function plot(data, labels, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Use bars to plot data.
bar(data,'group');
% Label the x axis, and set its limits tightly.
xlabel('Integral Values', labFmt{:});
set(gca, 'XLim', [0.5 10.5]);
% Label the y axis.
ylabel(ylabelString, labFmt{:});
% Set the title and y axis based on y axis label.
if strncmp(ylabelString, 'Integrand', 9)
title({'Multi-dimensional Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YTick', [-0.5 -0.3 -0.1 0.1 0.3 0.5]);
else
title({'Errors in Computing Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YScale', 'log');
set(gca, 'YDir', 'reverse');
end
legend([num2str(labels(1)) ' calls'], ...
[num2str(labels(2)) ' calls'], ...
[num2str(labels(3)) ' calls'], 'Location', 'Best')
function display_plot(data, labels, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Use bars to plot data.
bar(data,'group');
% Label the x axis, and set its limits tightly.
xlabel('Integral Values', labFmt{:});
set(gca, 'XLim', [0.5 10.5]);
% Label the y axis.
ylabel(ylabelString, labFmt{:});
% Set the title and y axis based on y axis label.
if strncmp(ylabelString, 'Integrand', 9)
title({'Multi-dimensional Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YTick', [-0.5 -0.3 -0.1 0.1 0.3 0.5]);
else
title({'Errors in Computing Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YScale', 'log');
set(gca, 'YDir', 'reverse');
end
legend([num2str(labels(1)) ' calls'], ...
[num2str(labels(2)) ' calls'], ...
[num2str(labels(3)) ' calls'], 'Location', 'Best')
```
```

** Maxcls too small to obtain required accuracy
Maxcls =      57    ifail = 1

Results so far (57 user function calls in last call of nag_quad_md_adapt_multi)

I     Integral   Estimated error
1      0.0422       0.0086
2      0.3998       0.0038
3      0.3898       0.0127
4      0.0214       0.0099
5     -0.3666       0.0020
6     -0.4176       0.0120
7     -0.0846       0.0110
8      0.3261       0.0001
9      0.4371       0.0112
10      0.1461       0.0119

** Maxcls too small to obtain required accuracy
Maxcls =     912    ifail = 1

Results so far (798 user function calls in last call of nag_quad_md_adapt_multi)

I     Integral   Estimated error
1      0.0384       0.0006
2      0.4012       0.0006
3      0.3952       0.0006
4      0.0258       0.0006
5     -0.3673       0.0006
6     -0.4227       0.0006
7     -0.0895       0.0006
8      0.3260       0.0006
9      0.4417       0.0006
10      0.1514       0.0006

I     Integral   Estimated error
1      0.0384       0.0004
2      0.4012       0.0003
3      0.3952       0.0003
4      0.0258       0.0003
5     -0.3672       0.0003
6     -0.4227       0.0003
7     -0.0895       0.0003
8      0.3260       0.0003
9      0.4417       0.0003
10      0.1514       0.0003

```
```function d01ea_example
% Initialize variables and arrays.
nfun = 10;
ndim = 4;
absreq = 0;
relreq = 1e-3;
mincls = 0;
ircls = 2^ndim+2*ndim*ndim+2*ndim+1;
maxcls = ircls;
lenwrk = 6*ndim+9*nfun+(ndim+nfun+2)*(1+maxcls/ircls);
w = 1:nfun;
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

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

% Turn off warnings from NAG routines, because we'll handle these ourselves.
% Then call the routine for the first time.  The initial parameter values
% are set up to demonstrate the continuation facility of the routine; it
% should exit with ifail = 1 (maxcls is too small).
warning('off', 'NAG:warning');
[mincls, wrkstr, finest, absest, ifail] = d01ea(a, b, int64(mincls), ...
int64(maxcls), int64(nfun), @f, absreq, relreq, wrkstr);

% Check exit condition.  First, deal with problems associated with
% input parameters on this initial call.
if (ifail == 2)
% Length of the work array is too small.  Print message and exit.
error('Warning: d01ea returned with ifail = %1d\n(Lenwrk is too small to obtain required accuracy).',ifail);
elseif (ifail == 4)
% Incorrect parameter values.  Print message and exit.
error('Warning: d01ea returned with ifail = %1d ',ifail);
end

% Continue loop until ifail isn't 1 or 3 (i.e. maxcls is big enough).
while (ifail == 1 || ifail == 3)
[mincls, wrkstr, finest, absest, ifail] = ...
d01ea(a, b, int64(mincls), int64(maxcls), int64(nfun), ...
@f, absreq, relreq, wrkstr);

% Is it big enough yet?  Output appropriate message.
if ifail ~= 0
fprintf('** Maxcls too small to obtain required accuracy\n');
fprintf('Maxcls = %7.0f    ifail = %1.0f \n\n',maxcls, ifail);
fprintf(['Results so far (%1.0f user function calls in last ', ...
'call of d01ea)\n\n'], mincls);
else
fprintf(['\nFinal results (%1.0f user function calls in last ', ...
'call of d01ea)\n\n'], mincls);
end

% Output results, and store them for plotting.
disp('  I     Integral   Estimated error')
for i = 1:nfun
fprintf(' %2.0f  ',w(i));
fprintf('%10.4f   %10.4f  \n',finest(i),absest(i))
fkeep(i,ncall) = finest(i);
akeep(i,ncall) = absest(i);
end

% Also, store the number of function calls.
mkeep(ncall,1) = mincls;

% Prepare for next call of routine.
mincls = -1;
maxcls = maxcls * mulfac;
ncall = ncall + 1;
if ifail ~= 0
fprintf('\n');
end
end
% Turn NAG warnings back on.
warning('on', 'NAG:warning');
% Plot results.
fig = figure('Number', 'off');
display_plot(fkeep, mkeep, 'Integrand Number');

fig = figure('Number', 'off');
display_plot(akeep, mkeep, 'Estimated Error');

function f = f(ndim, z, nfun)
% Evaluate the integrands.
f = zeros(nfun,1);
sm = 0;
for n=1:int32(ndim) % can't use int64 in loop range
sm = sm + double(n)*z(n);
end
for k=1:int32(nfun) % can't use int64 in loop range
f(k) = log(sm)*sin(double(k)+sm);
end
function plot(data, labels, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Use bars to plot data.
bar(data,'group');
% Label the x axis, and set its limits tightly.
xlabel('Integral Values', labFmt{:});
set(gca, 'XLim', [0.5 10.5]);
% Label the y axis.
ylabel(ylabelString, labFmt{:});
% Set the title and y axis based on y axis label.
if strncmp(ylabelString, 'Integrand', 9)
title({'Multi-dimensional Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YTick', [-0.5 -0.3 -0.1 0.1 0.3 0.5]);
else
title({'Errors in Computing Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YScale', 'log');
set(gca, 'YDir', 'reverse');
end
legend([num2str(labels(1)) ' calls'], ...
[num2str(labels(2)) ' calls'], ...
[num2str(labels(3)) ' calls'], 'Location', 'Best')
function display_plot(data, labels, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Use bars to plot data.
bar(data,'group');
% Label the x axis, and set its limits tightly.
xlabel('Integral Values', labFmt{:});
set(gca, 'XLim', [0.5 10.5]);
% Label the y axis.
ylabel(ylabelString, labFmt{:});
% Set the title and y axis based on y axis label.
if strncmp(ylabelString, 'Integrand', 9)
title({'Multi-dimensional Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YTick', [-0.5 -0.3 -0.1 0.1 0.3 0.5]);
else
title({'Errors in Computing Integrals of ten Integrands with ', ...
'Various Numbers of Function Calls'}, titleFmt{:});
set(gca, 'YScale', 'log');
set(gca, 'YDir', 'reverse');
end
legend([num2str(labels(1)) ' calls'], ...
[num2str(labels(2)) ' calls'], ...
[num2str(labels(3)) ' calls'], 'Location', 'Best')
```
```
d01ea example program results

** Maxcls too small to obtain required accuracy
Maxcls =      57    ifail = 1

Results so far (57 user function calls in last call of d01ea)

I     Integral   Estimated error
1      0.0422       0.0086
2      0.3998       0.0038
3      0.3898       0.0127
4      0.0214       0.0099
5     -0.3666       0.0020
6     -0.4176       0.0120
7     -0.0846       0.0110
8      0.3261       0.0001
9      0.4371       0.0112
10      0.1461       0.0119

** Maxcls too small to obtain required accuracy
Maxcls =     912    ifail = 1

Results so far (798 user function calls in last call of d01ea)

I     Integral   Estimated error
1      0.0384       0.0006
2      0.4012       0.0006
3      0.3952       0.0006
4      0.0258       0.0006
5     -0.3673       0.0006
6     -0.4227       0.0006
7     -0.0895       0.0006
8      0.3260       0.0006
9      0.4417       0.0006
10      0.1514       0.0006

Final results (912 user function calls in last call of d01ea)

I     Integral   Estimated error
1      0.0384       0.0004
2      0.4012       0.0003
3      0.3952       0.0003
4      0.0258       0.0003
5     -0.3672       0.0003
6     -0.4227       0.0003
7     -0.0895       0.0003
8      0.3260       0.0003
9      0.4417       0.0003
10      0.1514       0.0003

```