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_1d_gen_vec_multi_rcomm (d01ra)

Purpose

nag_quad_1d_gen_vec_multi_rcomm (d01ra) is a general purpose adaptive integrator which calculates an approximation to a vector of definite integrals FF over a finite range [a,b] [a,b] , given the vector of integrands f(x)f(x).
F = ab f(x) d x
F = a b f(x) d x
If the same subdivisions of the range are equally good for functions f1(x)f1(x) and f2(x)f2(x), because f1(x)f1(x) and f2(x)f2(x) have common areas of the range where they vary slowly and where they vary quickly, then we say that f1(x)f1(x) and f2(x)f2(x) are ‘similar’. nag_quad_1d_gen_vec_multi_rcomm (d01ra) is particularly effective for the integration of a vector of similar functions.

Syntax

[irevcm, sid, needi, x, nx, dinest, errest, icom, com, ifail] = d01ra(irevcm, a, b, needi, x, nx, fm, dinest, errest, iopts, opts, icom, com, 'ni', ni, 'lenx', lenx)
[irevcm, sid, needi, x, nx, dinest, errest, icom, com, ifail] = nag_quad_1d_gen_vec_multi_rcomm(irevcm, a, b, needi, x, nx, fm, dinest, errest, iopts, opts, icom, com, 'ni', ni, 'lenx', lenx)

Description

nag_quad_1d_gen_vec_multi_rcomm (d01ra) is an extension to various QUADPACK routines, including QAG, QAGS and QAWSE. The extensions made allow multiple integrands to be evaluated simultaneously, using a vectorized interface and reverse communication.
The quadrature scheme employed by nag_quad_1d_gen_vec_multi_rcomm (d01ra) can be chosen by you. Six Gauss–Kronrod schemes are available. The algorithm incorporates a global acceptance criterion (as defined by Malcolm and Simpson (1976)), optionally together with the ε-algorithm (see Wynn (1956)) to perform extrapolation. The local error estimation is described in Piessens et al. (1983).
nag_quad_1d_gen_vec_multi_rcomm (d01ra) is the integration function in the suite of functions nag_quad_1d_gen_vec_multi_rcomm (d01ra) and nag_quad_1d_gen_vec_multi_dimreq (d01rc). It also uses optional parameters, which can be set and queried using the functions nag_quad_opt_set (d01zk) and nag_quad_opt_get (d01zl) respectively. The options available are described in Section [Optional Parameters].
First, the option arrays iopts and opts must be initialized using nag_quad_opt_set (d01zk). Thereafter any required options must be set before calling nag_quad_1d_gen_vec_multi_rcomm (d01ra), or the function nag_quad_1d_gen_vec_multi_dimreq (d01rc). The options available are described in Section [Optional Parameters].
A typical usage of this suite of functions is (in pseudo-code for clarity),
Setup phase
iopts = zeros(100, 1, nag_int_name);
opts  = zeros(100, 1);

% initialize option arrays
[iopts, opts, ifail] = d01zk('Initialize = d01ra', iopts, opts);

% set any non-default options required
[iopts, opts, ifail] = d01zk('option = value', iopts, opts);
...

% determine maximum required array lengths
[lenxrq, ldfmrq, sdfmrq, licmin, licmax, lcmin, lcmax, ifail] = ...
      d01rc(ni, iopts, opts);

% allocate remaining arrays
needi  = zeros(ni, 1, nag_int_name);
com    = zeros(lcmax, 1);
icom   = zeros(licmax, 1, nag_int_name);
fm     = zeros(ldfmrq, sdfmrq);
dinest = zeros(ni, 1); 
errest = zeros(ni, 1); 
x      = zeros(1, lenxrq);
Solve phase
irevcm = nag_int(1);

while irevcm ~= 0
  [irevcm, sid, needi, x, nx, dinest, errest, icomm, comm, ifail] = ...
    d01ra(irevcm, a, b, needi, x, nx, fm, dinest, errest, ...
          iopts, opts, icom, com);

  switch irevcm
    case 11
        Initial solve phase 
        evaluate fm(1:ni,1:nx)
    case(12)
        Adaptive solve phase 
        evaluate fm(needi(1:ni)=1,1:nx)
    case(0)
        Final return
  end
Diagnostic phase
[ivalue, rvalue, cvalue, optype, ifail] = d01zl('option', iopts, opts);
During the initial solve phase, the first estimation of the definite integral and error estimate is constructed over the interval [a,b] [a,b] . This will have been divided into spri spri  level 11 segments, where spri spri  is the number of Primary Divisions, and will use any provided breakpoints if Primary Division Mode = MANUALPrimary Division Mode=MANUAL.
Once a complete integral estimate over [a,b] [a,b]  is available, i.e., after all the estimates for the level 1 segments have been evaluated, the function enters the adaptive phase. The estimated errors are tested against the requested tolerances εa εa  and εr εr , corresponding to the Absolute Tolerance and Relative Tolerance respectively. Should this test fail, and additional subdivision be allowed, a segment is selected for subdivision, and is subsequently replaced by two new segments at the next level of refinement. How this segment is chosen may be altered by setting Prioritize Error to either favour the segment with the maximum error, or the segment with the lowest level supporting an unacceptable (although potentially non-maximal) error. Up to maxsdiv maxsdiv  subdivisions are allowed provided sufficient memory, where maxsdiv maxsdiv  is the value of Maximum Subdivisions.
Once a sufficient number of error estimates have been constructed for a particular integral, the function may optionally use Extrapolation to attempt to accelerate convergence. This may significantly lower the amount of work required for a given integration. To minimize the risk of premature convergence from extrapolation, a safeguard εsafe εsafe  can be set using Extrapolation Safeguard, and the extrapolated solution will only be considered if εsafe εq εex εsafe εq εex , where εq εq  and εex εex  are the estimated error directly from the quadrature and from the extrapolation respectively.

References

Malcolm M A and Simpson R B (1976) Local versus global strategies for adaptive quadrature ACM Trans. Math. Software 1 129–146
Piessens R (1973) An algorithm for automatic integration Angew. Inf. 15 399–401
Piessens R, de Doncker–Kapenga E, Überhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration Springer–Verlag
Wynn P (1956) On a device for computing the em(Sn)em(Sn) transformation Math. Tables Aids Comput. 10 91–96

Parameters

Note: this function uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the parameter irevcm. Between intermediate exits and re-entries, all parameters other than irevcm, needi and fm must remain unchanged.

Compulsory Input Parameters

1:     irevcm – int64int32nag_int scalar
On initial entry: irevcm = 1irevcm=1.
irevcm = 1irevcm=1
Set up data structures in icom and com and start a new integration.
Constraint: irevcm = 1irevcm=1 on initial entry.
On intermediate re-entry: irevcm should normally be left unchanged. However, if irevcm is set to a negative value, nag_quad_1d_gen_vec_multi_rcomm (d01ra) will terminate, (see irevcm = 11irevcm=11 and irevcm = 12irevcm=12 above).
2:     a – double scalar
aa, the lower bound of the domain.
3:     b – double scalar
bb, the upper bound of the domain.
If |ba| < 10ε|b-a|<10ε, where εε is the machine precision (see nag_machine_precision (x02aj)), then nag_quad_1d_gen_vec_multi_rcomm (d01ra) will return dinest(j) = errest(j) = 0.0dinestj=errestj=0.0, for j = 1,2,,nij=1,2,,ni.
4:     needi(ni) – int64int32nag_int array
On initial entry: need not be set.
On intermediate re-entry: needi(j)needij may be used to indicate what action you have taken for integral jj.
needi(j) = 1needij=1
You have provided the values fj(xi)fj(xi) in fm(j,i)fmji, for i = 1,2,,nxi=1,2,,nx.
needi(j) < 0needij<0
You are abandoning the evaluation of integral jj. The current values of dinest(j)dinestj and errest(j)errestj will be returned on final completion.
Otherwise you have not provided the value fj(xi)fj(xi).
5:     x(lenx) – double array
lenx, the dimension of the array, must satisfy the constraint lenxlenxrqlenxlenxrq, where lenxrqlenxrq is dependent upon the options currently set (see Section [Optional Parameters]). lenxrqlenxrq is returned as lenxrq from nag_quad_1d_gen_vec_multi_dimreq (d01rc).
On initial entry: if Primary Division Mode = AUTOMATICPrimary Division Mode=AUTOMATIC, x need not be set. This is the default behaviour.
If Primary Division Mode = MANUALPrimary Division Mode=MANUAL, x may be used to supply a set of initial ‘break points’ inside the domain of integration. Specifically, x(i)xi must contain a break point xi0xi0, for i = 1,2,,(spri1)i=1,2,,(spri-1), where sprispri is the number of Primary Divisions.
Constraint: if break points are supplied, xi0(a,b)xi0(a,b), |xi0a| > 10.0ε|xi0-a|>10.0ε, |xi0b| > 10.0ε|xi0-b|>10.0ε, for i = 1,2,,(spri1)i=1,2,,(spri-1).
6:     nx – int64int32nag_int scalar
On initial entry: need not be set.
On intermediate re-entry: must not be changed.
7:     fm(ldfm, : :) – double array
The first dimension of the array fm must be at least ldfmrqldfmrq, where ldfmrqldfmrq is dependent upon nini and the options currently set. ldfmrqldfmrq is returned as ldfmrq from nag_quad_1d_gen_vec_multi_dimreq (d01rc). If default options are chosen, ldfmrq = nildfmrq=ni, implying ldfmnildfmni
The second dimension of the array must be at least sdfmrqsdfmrq, where sdfmrqsdfmrq is dependent upon nini and the options currently set. sdfmrqsdfmrq is returned as sdfmrq from nag_quad_1d_gen_vec_multi_dimreq (d01rc). If default options are chosen, sdfmrq = lenxrqsdfmrq=lenxrq
On initial entry: need not be set.
On intermediate re-entry: if indicated by needi(j)needij you must supply the values fj(xi)fj(xi) in fm(j,i)fmji, for i = 1,2,,nxi=1,2,,nx and j = 1,2,,nij=1,2,,ni.
8:     dinest(ni) – double array
dinest(j)dinestj contains the current estimate of the definite integral FjFj.
On initial entry: need not be set.
On intermediate re-entry: must not be altered.
9:     errest(ni) – double array
errest(j)errestj contains the current error estimate of the definite integral FjFj.
On initial entry: need not be set.
On intermediate re-entry: must not be altered.
10:   iopts(100100) – int64int32nag_int array
Integer option array generated by nag_quad_opt_set (d01zk).
Constraint: iopts must not be altered between calls to nag_quad_1d_gen_vec_multi_rcomm (d01ra), nag_quad_1d_gen_vec_multi_dimreq (d01rc), nag_quad_opt_set (d01zk) and nag_quad_opt_get (d01zl).
11:   opts(100100) – double array
Real option array generated by nag_quad_opt_set (d01zk).
Constraint: opts must not be altered between calls to nag_quad_1d_gen_vec_multi_rcomm (d01ra), nag_quad_1d_gen_vec_multi_dimreq (d01rc), nag_quad_opt_set (d01zk) and nag_quad_opt_get (d01zl).
12:   icom(licom) – int64int32nag_int array
licom, the dimension of the array, must satisfy the constraint licomlicminlicomlicmin, where licminlicmin is dependent upon ni and the current options set. licminlicmin is returned as licmin from nag_quad_1d_gen_vec_multi_dimreq (d01rc). If the default options are set, then licmin = 55 + 6 × nilicmin=55+6×ni. Larger values than licminlicmin are recommended if you anticipate that any integrals will require the domain to be further subdivided..
The contents of this array must not be changed directly once the integration process has started. icom contains details of the integration procedure, including information on the integration of the nini integrals over individual segments. This data is stored sequentially in the order that segments are created. You are free to examine the contents of icom using the diagnostic function at any time during, or after, the integration process.
13:   com(lcom) – double array
lcom, the dimension of the array, must satisfy the constraint lcom > lcminlcom>lcmin, where lcminlcmin is dependent upon ni, sprispri and the current options set. lcminlcmin is returned as lcmin from nag_quad_1d_gen_vec_multi_dimreq (d01rc). If default options are set, then lcmin = 96 + 12 × ni lcmin =96+12×ni. Larger values are recommended if you anticipate that any integrals will require the domain to be further subdivided.
Given the current options and parameters, the maximum value, lcmaxlcmax, of lcom that may be required, is returned as lcmax from nag_quad_1d_gen_vec_multi_dimreq (d01rc). If default options are chosen, lcmax = 94 + 9 × ni + ni / 2 + (spri + 100) × (2 + ni / 2 + 2 × ni)lcmax=94+9×ni+ni/2 +(spri+100)×(2+ni/2+2×ni).
.
The contents of this array must not be changed directly once the integration process has started. com contains details of the integration procedure, including information on the integration of the nini integrals over individual segments. This data is stored sequentially in the order that segments are created. You are free to examine the contents of com using the diagnostic function at any time during, or after, the solution phase.

Optional Input Parameters

1:     ni – int64int32nag_int scalar
Default: The dimension of the arrays needi, dinest, errest and the first dimension of the array fm. (An error is raised if these dimensions are not equal.)
nini, number of integrands.
2:     lenx – int64int32nag_int scalar
Default: The dimension of the array x.
The dimension of the array x as declared in the (sub)program from which nag_quad_1d_gen_vec_multi_rcomm (d01ra) is called. Currently lenx = max (122,spri1)lenx=max(122,spri-1) will be sufficient for all cases.
Constraint: lenxlenxrqlenxlenxrq, where lenxrqlenxrq is dependent upon the options currently set (see Section [Optional Parameters]). lenxrqlenxrq is returned as lenxrq from nag_quad_1d_gen_vec_multi_dimreq (d01rc).

Input Parameters Omitted from the MATLAB Interface

ldfm licom lcom

Output Parameters

1:     irevcm – int64int32nag_int scalar
On intermediate exit: irevcm = 11irevcm=11 or 1212.
irevcm requests the integrands fj(xi)fj(xi) be evaluated for all required j1,,nij1,,ni as indicated by needi, and at all the points xixi, for i = 1,2,,nxi=1,2,,nx. Abscissae xixi are provided in x(i)xi and fj(xi)fj(xi) must be returned in fm(j,i)fmji.
During the initial solve phase:
irevcm = 11irevcm=11
Function values are required to construct the initial estimates of the definite integrals.
If needi(j) = 1needij=1, fj(xi)fj(xi) must be supplied in fm(j,i)fmji. This will be the case unless you have abandoned the evaluation of specific integrals on a previous call.
If needi(j) = 0needij=0, you have previously abandoned the evaluation of integral jj, and hence should not supply the value of fj(xi)fj(xi).
dinest and errest contain incomplete information during this phase. As such you should not abandon the evaluation of any integrals during this phase unless you do not require their estimate.
If irevcm is set to a negative value during this phase, needi(j)needij, for j = 1,2,,nij=1,2,,ni, will be set to this negative value and ifail = 1ifail=-1 will be returned.
During the adaptive solve phase:
irevcm = 12irevcm=12
Function values are required to improve the estimates of the definite integrals.
If needi(j) = 0needij=0, any evaluation of fj(xi)fj(xi) will be discarded, so there is no need to provide them.
If needi(j) = 1needij=1, fj(xi)fj(xi) must be returned.
If needi(j) = 2needij=2, 33 or 44, the current error estimate of integral jj does not require integrand jj to be evaluated. Should you choose to, integrand jj can be evaluated in which case needi(j)needij must be set to 11.
dinest and errest contain complete information during this phase.
If irevcm is set to a negative value during this phase either ifail = 1ifail=1, 22 or 33 will be returned and the elements of needi will reflect the current state of of the adaptive process.
On final exit: irevcm = 0irevcm=0.
irevcm = 0irevcm=0
Indicates the algorithm has completed.
2:     sid – int64int32nag_int scalar
On intermediate exit: sid identifies a specific set of abscissae, xx, returned during the integration process. When a new set of abscissae are generated the value of sid is incremented by 11. Advanced users may store calculations required for an identified set xx, and reuse them should nag_quad_1d_gen_vec_multi_rcomm (d01ra) return the same value of sid, i.e., the same set of abscissae was used.
3:     needi(ni) – int64int32nag_int array
On intermediate exit: needi(j)needij indicates what action must be taken for integral j = 1,2,nij=1,2,ni (see irevcm).
needi(j) = 0needij=0
Do not provide fj(xi)fj(xi). Any provided values will be ignored.
needi(j) = 1needij=1
The values fj(xi)fj(xi) must be provided in fm(j,i)fmji, for i = 1,2,,nxi=1,2,,nx.
needi(j) = 2needij=2
The values fj(xi)fj(xi) are not definitely required, however the error estimate for integral jj is still above the requested tolerance. If you wish to provide values for the evaluation of integral jj, set needi(j) = 1needij=1, and supply fj(xi)fj(xi) in fm(j,i)fmji, for i = 1,2,,nxi=1,2,,nx.
needi(j) = 3needij=3
The error estimate for integral jj cannot be improved to below the requested tolerance directly, either because no more new splits may be performed due to exhaustion, or due to the detection of extremely bad integrand behaviour. However, providing the values fj(xi)fj(xi) may still lead to some improvement, and may lead to an acceptable error estimate indirectly using Wynn's epsilon algorithm. If you wish to provide values for the evaluation of integral jj, set needi(j) = 1needij=1, and supply fj(xi)fj(xi) in fm(j,i)fmji, for i = 1,2,,nxi=1,2,,nx.
needi(j) = 4needij=4
The error estimate of integral jj is below the requested tolerance. If you believe this to be false, if for example the result in dinest(j)dinestj is greatly different to what you may expect, you may force the algorithm to re-evaluate this conclusion by including the evaluations of integrand jj at xixi, for i = 1,2,,nxi=1,2,,nx, and setting needi(j) = 1needij=1. Integral and error estimation will be performed again during the next iteration.
On final exit: needi(j)needij indicates the final state of integral jj.
needi(j) = 0needij=0
The error estimate for FjFj is below the requested tolerance.
needi(j) = 1needij=1
The error estimate for FjFj is below the requested tolerance after extrapolation.
needi(j) = 2needij=2
The error estimate for FjFj is above the requested tolerance.
needi(j) = 3needij=3
The error estimate for FjFj is above the requested tolerance, and extremely bad behaviour of integral jj has been detected.
needi(j) < 0needij<0
You prohibited further evaluation of integral jj.
4:     x(lenx) – double array
On intermediate exit: x(i)xi is the abscissa xixi, for i = 1,2,,nxi=1,2,,nx, at which the appropriate integrals must be evaluated.
5:     nx – int64int32nag_int scalar
On intermediate exit: nxnx, the number of abscissae at which integrands are required.
6:     dinest(ni) – double array
dinest(j)dinestj contains the current estimate of the definite integral FjFj.
Contains the current estimates of the ni integrals. If irevcm = 0irevcm=0, this will be the final solution.
7:     errest(ni) – double array
errest(j)errestj contains the current error estimate of the definite integral FjFj.
Contains the current error estimates for the ni integrals. If irevcm = 0irevcm=0, errest contains the final error estimates of the ni integrals.
8:     icom(licom) – int64int32nag_int array
9:     com(lcom) – double array
10:   ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_quad_1d_gen_vec_multi_rcomm (d01ra) may return useful information for one or more of the following detected errors or 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 = 1ifail=1
At least one error estimate exceeded the requested tolerances.
W ifail = 2ifail=2
Extremely bad behaviour was detected for at least one integral.
W ifail = 3ifail=3
Extremely bad behaviour was detected for at least one integral. At least one other integral error estimate was above the requested tolerance.
  ifail = 11ifail=11
irevcm had an illegal value.
  ifail = 21ifail=21
Constraint: ni1ni1.
  ifail = 71ifail=71
On entry, Primary Division Mode = MANUALPrimary Division Mode=MANUAL and at least one supplied breakpoint in x is outside of the domain of integration.
  ifail = 81ifail=81
  ifail = 111ifail=111
W ifail = 171ifail=171
W ifail = 191ifail=191
  ifail = 1001ifail=1001
Either the option arrays iopts and opts have not been initialized for nag_quad_1d_gen_vec_multi_rcomm (d01ra), or they have become corrupted.
  ifail = 1101ifail=1101
On entry, one of icom and com has become corrupted.
W ifail = 1ifail=-1
Evaluation of all integrals has been stopped during the initial phase.

Accuracy

nag_quad_1d_gen_vec_multi_rcomm (d01ra) cannot guarantee, but in practice usually achieves, the following accuracy for each integral FjFj:
|Fjdinest(j)| tol
| Fj - dinestj | tol
where
tol = max (εa, εr × |Fj| )
tol = max(εa, εr × |Fj| )
εaεa and εrεr are the error tolerances Absolute Tolerance and Relative Tolerance respectively. Moreover, it returns errest, the entries of which in normal circumstances satisfy,
|Fjdinest(j)| errest(j) tol .
| Fj - dinestj | errestj tol .

Further Comments

The time required by nag_quad_1d_gen_vec_multi_rcomm (d01ra) is usually dominated by the time required to evaluate the values of the integrands fjfj.
nag_quad_1d_gen_vec_multi_rcomm (d01ra) will be most efficient if any badly behaved integrands provided have irregularities over similar subsections of the domain. For example, evaluation of the integrals,
∫ 01
  log(x) x − (1/2) x2  
dx
0 1 log(x) x-12 x2 dx
will be quite efficient, as the irregular behaviour of the first two integrands is at x = 0x=0. On the contrary, the evaluation of the integrals,
01
(log(x))
log(1x)
dx
0 1 log(x) log(1-x) dx
will be less efficient, as the two integrands have singularities at opposite ends of the domain, which will result in subdivisions which are only of use to one integrand. In such cases, it will be more efficient to use two sets of calls to nag_quad_1d_gen_vec_multi_rcomm (d01ra).
nag_quad_1d_gen_vec_multi_rcomm (d01ra) will flag extremely bad behaviour if a sub-interval k k-  with bounds [ak,bk] [ak-,bk-]  satisfying |bkak| < max (δa, δr × |ba| ) | bk- - ak- | < max(δa, δr × |b-a| )  has a local error estimate greater than the requested tolerance for at least one integral. The values δaδa and δrδr can be set through the optional parameters Absolute Interval Minimum and Relative Interval Minimum respectively. You may try using nag_quad_1d_gen_vec_multi_rcomm (d01ra) again over such sub-intervals with only the integrands that exhibited bad behaviour and a new set of communication arrays, or you may use another quadrature routine designed to deal with any known specific singular behaviour over this sub-interval. If these provide sufficient accuracy over the specific sub-intervals, you may correct the definite integral evaluations over [a,b][a,b] by subtracting the inadequate sub-interval estimations from dinest and errest and adding the acceptable sub-interval estimations generated by the direct integration of the sub-interval where bad behaviour was detected.

Example

function nag_quad_1d_gen_vec_multi_rcomm_example

  % Setup phase.

  % set problem parameters
  ni = int64(2);
  nx = int64(0);
  % lower (a) and upper (b) bounds
  a = 0;
  b = pi;
  iopts = zeros(100, 1, 'int64');
  opts  = zeros(100, 1);

  % initialize option arrays
  [iopts, opts, ifail] = nag_quad_opt_set('Initialize = nag_quad_1d_gen_vec_multi_rcomm', iopts, opts);

  % set any non-default options required
  [iopts, opts, ifail] = nag_quad_opt_set('Quadrature Rule = gk41', iopts, opts);
  [iopts, opts, ifail] = nag_quad_opt_set('Absolute Tolerance = 1.0e-7', iopts, opts);
  [iopts, opts, ifail] = nag_quad_opt_set('Relative Tolerance = 1.0e-7', iopts, opts);

  % determine maximum required array lengths
  [lenxrq, ldfmrq, sdfmrq, licmin, licmax, lcmin, lcmax, ifail] = ...
        nag_quad_1d_gen_vec_multi_dimreq(ni, iopts, opts);

  % allocate remaining arrays
  needi  = zeros(ni, 1, 'int64');
  comm   = zeros(lcmax, 1);
  icomm  = zeros(licmax, 1, 'int64');
  fm     = zeros(ldfmrq, sdfmrq);
  dinest = zeros(ni, 1);
  errest = zeros(ni, 1);
  x      = zeros(1, lenxrq);

  % Solve phase.

  % Use nag_quad_1d_gen_vec_multi_rcomm to evaluate the definate integrals of:
  %   f_1 = (x*sin(2*x))*cos(15*x)
  %   f_2 = (x*sin(2*x))*(x*cos(50*x))

  % set initial irevcm
  irevcm = int64(1);

  while irevcm ~= 0
    [irevcm, sid, needi, x, nx, dinest, errest, icomm, comm, ifail] = ...
      nag_quad_1d_gen_vec_multi_rcomm(irevcm, a, b, needi, x, nx, fm, dinest, errest, ...
            iopts, opts, icomm, comm);

    switch irevcm
      case 11
        % Initial returns.
        % These will occur during the non-adaptive phase.
        % All values must be supplied.
        % dinest and errest do not contain approximations
        % over the complete interval at this stage.

        % Calculate x*sin(2*x), storing the result in fm(2,1:nx) for re-use.
        fm(2, :) = x.*sin(2*x);

        % Calculate f_1
        fm(1, :) = fm(2, :).*cos(15*x);

        % Calculate f_2
        fm(2, :) = fm(2, :).*x.*cos(50*x);
      case 12
        % Intermediate returns.
        % These will occur during the adaptive phase.
        % All requested values must be supplied.
        % dinest and errest do not contain approximations
        % over the complete interval at this stage.

        % Calculate x*sin(2*x).
        fm(2, :) = x.*sin(2*x);

        % Calculate f_1 if required
        if needi(1) == 1
          fm(1, :) = fm(2, :).*cos(15*x);
        end

        % Complete f_2 calculation if required.
        if needi(2) == 1
          fm(2, :) = fm(2, :).*x.*cos(50*x);
        end
      case 0
        % Final return
    end
  end

  % query some currently set options and statistics.
  [ivalue, rvalue, cvalue, optype, ifail] = nag_quad_opt_get('Quadrature rule', iopts, opts);
  display_option('Quadrature rule',optype,ivalue,rvalue,cvalue);
  [ivalue, rvalue, cvalue, optype, ifail] = nag_quad_opt_get('Maximum Subdivisions', iopts, opts);
  display_option('Maximum Subdivisions',optype,ivalue,rvalue,cvalue);
  [ivalue, rvalue, cvalue, optype, ifail] = nag_quad_opt_get('Extrapolation', iopts, opts);
  display_option('Extrapolation',optype,ivalue,rvalue,cvalue);
  [ivalue, rvalue, cvalue, optype, ifail] = nag_quad_opt_get('Extrapolation Safeguard', iopts, opts);
  display_option('Extrapolation Safeguard',optype,ivalue,rvalue,cvalue);

  % print solution
  fprintf('\nIntegral |  needi  |   dinest   |   errest   \n');
  for j=1:ni
    fprintf('%9d %9d %12.4e %12.4e\n', j, needi(j), dinest(j), errest(j));
  end


function [dinest, errest, user] = monit(ni, ns, dinest, errest, fcount, ...
                                        sinfoi, evals, ldi, sinfor, fs, ...
                                        es, ldr, user)
  % Display information on individual segments
  fprintf('\nInformation on splitting and evaluations over subregions.\n');
  for k=1:ns
    sid = sinfoi(1,k);
    parent = sinfoi(2,k);
    child1 = sinfoi(3,k);
    child2 = sinfoi(4,k);
    level = sinfoi(5,k);
    lbnd = sinfor(1,k);
    ubnd = sinfor(2,k);
    fprintf('\nSegment %3d Sid = %3d Parent = %3d Level = %3d.\n', k, sid, parent, level);
    if (child1>0)
      fprintf('Children = (%3d, %3d)\n', child1, child2);
    end
    fprintf('Bounds (%11.4e, %11.4e)\n', lbnd, ubnd);
    for j = 1:ni
      if (evals(j,k) ~= 0)
        fprintf('Integral %2d approximation %11.4e\n', j, fs(j,k));
        fprintf('Integral %2d error estimate %11.4e\n', j, es(j,k));
        if (evals(j,k) ~= 1)
          fprintf('Integral %2d evaluation has been superseded by descendants.\n', j);
        end
      end
    end
  end
function display_option(optstr,optype,ivalue,rvalue,cvalue)
  % Query optype and print the appropriate option values

  switch optype
    case 1
      fprintf('%30s: %13d\n', optstr, ivalue);
    case 2
      fprintf('%30s: %13.4e\n', optstr, rvalue);
    case 3
      fprintf('%30s: %16s\n', optstr, cvalue);
    case 4
      fprintf('%30s: %3d  %16s\n', optstr, ivalue, cvalue);
    case 5
      fprintf('%30s: %14.4e  %16s\n', optstr, rvalue, cvalue);
  end
 
               Quadrature rule: GK41                            
          Maximum Subdivisions:            50
                 Extrapolation: ON                              
       Extrapolation Safeguard:    1.0000e-12

Integral |  needi  |   dinest   |   errest   
        1         0  -2.8431e-02   1.1234e-14
        2         0   7.9083e-03   2.6600e-09

function d01ra_example

  % Setup phase.

  % set problem parameters
  ni = int64(2);
  nx = int64(0);
  % lower (a) and upper (b) bounds
  a = 0;
  b = pi;
  iopts = zeros(100, 1, 'int64');
  opts  = zeros(100, 1);

  % initialize option arrays
  [iopts, opts, ifail] = d01zk('Initialize = d01ra', iopts, opts);

  % set any non-default options required
  [iopts, opts, ifail] = d01zk('Quadrature Rule = gk41', iopts, opts);
  [iopts, opts, ifail] = d01zk('Absolute Tolerance = 1.0e-7', iopts, opts);
  [iopts, opts, ifail] = d01zk('Relative Tolerance = 1.0e-7', iopts, opts);

  % determine maximum required array lengths
  [lenxrq, ldfmrq, sdfmrq, licmin, licmax, lcmin, lcmax, ifail] = ...
        d01rc(ni, iopts, opts);

  % allocate remaining arrays
  needi  = zeros(ni, 1, 'int64');
  comm   = zeros(lcmax, 1);
  icomm  = zeros(licmax, 1, 'int64');
  fm     = zeros(ldfmrq, sdfmrq);
  dinest = zeros(ni, 1);
  errest = zeros(ni, 1);
  x      = zeros(1, lenxrq);

  % Solve phase.

  % Use d01ra to evaluate the definate integrals of:
  %   f_1 = (x*sin(2*x))*cos(15*x)
  %   f_2 = (x*sin(2*x))*(x*cos(50*x))

  % set initial irevcm
  irevcm = int64(1);

  while irevcm ~= 0
    [irevcm, sid, needi, x, nx, dinest, errest, icomm, comm, ifail] = ...
      d01ra(irevcm, a, b, needi, x, nx, fm, dinest, errest, ...
            iopts, opts, icomm, comm);

    switch irevcm
      case 11
        % Initial returns.
        % These will occur during the non-adaptive phase.
        % All values must be supplied.
        % dinest and errest do not contain approximations
        % over the complete interval at this stage.

        % Calculate x*sin(2*x), storing the result in fm(2,1:nx) for re-use.
        fm(2, :) = x.*sin(2*x);

        % Calculate f_1
        fm(1, :) = fm(2, :).*cos(15*x);

        % Calculate f_2
        fm(2, :) = fm(2, :).*x.*cos(50*x);
      case 12
        % Intermediate returns.
        % These will occur during the adaptive phase.
        % All requested values must be supplied.
        % dinest and errest do not contain approximations
        % over the complete interval at this stage.

        % Calculate x*sin(2*x).
        fm(2, :) = x.*sin(2*x);

        % Calculate f_1 if required
        if needi(1) == 1
          fm(1, :) = fm(2, :).*cos(15*x);
        end

        % Complete f_2 calculation if required.
        if needi(2) == 1
          fm(2, :) = fm(2, :).*x.*cos(50*x);
        end
      case 0
        % Final return
    end
  end

  % query some currently set options and statistics.
  [ivalue, rvalue, cvalue, optype, ifail] = d01zl('Quadrature rule', iopts, opts);
  display_option('Quadrature rule',optype,ivalue,rvalue,cvalue);
  [ivalue, rvalue, cvalue, optype, ifail] = d01zl('Maximum Subdivisions', iopts, opts);
  display_option('Maximum Subdivisions',optype,ivalue,rvalue,cvalue);
  [ivalue, rvalue, cvalue, optype, ifail] = d01zl('Extrapolation', iopts, opts);
  display_option('Extrapolation',optype,ivalue,rvalue,cvalue);
  [ivalue, rvalue, cvalue, optype, ifail] = d01zl('Extrapolation Safeguard', iopts, opts);
  display_option('Extrapolation Safeguard',optype,ivalue,rvalue,cvalue);

  % print solution
  fprintf('\nIntegral |  needi  |   dinest   |   errest   \n');
  for j=1:ni
    fprintf('%9d %9d %12.4e %12.4e\n', j, needi(j), dinest(j), errest(j));
  end


function [dinest, errest, user] = monit(ni, ns, dinest, errest, fcount, ...
                                        sinfoi, evals, ldi, sinfor, fs, ...
                                        es, ldr, user)
  % Display information on individual segments
  fprintf('\nInformation on splitting and evaluations over subregions.\n');
  for k=1:ns
    sid = sinfoi(1,k);
    parent = sinfoi(2,k);
    child1 = sinfoi(3,k);
    child2 = sinfoi(4,k);
    level = sinfoi(5,k);
    lbnd = sinfor(1,k);
    ubnd = sinfor(2,k);
    fprintf('\nSegment %3d Sid = %3d Parent = %3d Level = %3d.\n', k, sid, parent, level);
    if (child1>0)
      fprintf('Children = (%3d, %3d)\n', child1, child2);
    end
    fprintf('Bounds (%11.4e, %11.4e)\n', lbnd, ubnd);
    for j = 1:ni
      if (evals(j,k) ~= 0)
        fprintf('Integral %2d approximation %11.4e\n', j, fs(j,k));
        fprintf('Integral %2d error estimate %11.4e\n', j, es(j,k));
        if (evals(j,k) ~= 1)
          fprintf('Integral %2d evaluation has been superseded by descendants.\n', j);
        end
      end
    end
  end
function display_option(optstr,optype,ivalue,rvalue,cvalue)
  % Query optype and print the appropriate option values

  switch optype
    case 1
      fprintf('%30s: %13d\n', optstr, ivalue);
    case 2
      fprintf('%30s: %13.4e\n', optstr, rvalue);
    case 3
      fprintf('%30s: %16s\n', optstr, cvalue);
    case 4
      fprintf('%30s: %3d  %16s\n', optstr, ivalue, cvalue);
    case 5
      fprintf('%30s: %14.4e  %16s\n', optstr, rvalue, cvalue);
  end
 
               Quadrature rule: GK41                            
          Maximum Subdivisions:            50
                 Extrapolation: ON                              
       Extrapolation Safeguard:    1.0000e-12

Integral |  needi  |   dinest   |   errest   
        1         0  -2.8431e-02   1.1234e-14
        2         0   7.9083e-03   2.6600e-09

Optional Parameters

This section can be skipped if you wish to use the default values for all optional parameters, otherwise, the following is a list of the optional parameters available. A full description of each optional parameter is provided in Section [Description of the optional parameters].

Description of the Optional Parameters

For each option, we give a summary line, a description of the optional parameter and details of constraints.
The summary line contains:
The following symbol represents various machine constants:
All options accept the value ‘DEFAULT’ in order to return single options to their default states.
Keywords and character values are case insensitive, however they must be separated by at least one space.
For nag_quad_1d_gen_vec_multi_rcomm (d01ra) the maximum length of the parameter cvalue used by nag_quad_opt_get (d01zl) is 1515.
Absolute Interval Minimum  rr
Default = 128.0ε=128.0ε
r = δar=δa, the absolute lower limit for a segment to be considered for subdivision. See also Relative Interval Minimum and Section [Further Comments].
Constraint: r128εr128ε.
Absolute Tolerance  rr
Default = 1024ε=1024ε
r = εar=εa, the absolute tolerance required. See also Relative Tolerance and Section [Description].
Constraint: r0.0r0.0.
Extrapolation  aa
Default = ON=ON
Activate or deactivate the use of the εε algorithm (Wynn (1956)). Extrapolation often reduces the number of iterations required to achieve the desired solution, but it can occasionally lead to premature convergence towards an incorrect answer.
ONON
Use extrapolation.
OFFOFF
Disable extrapolation.
Extrapolation Safeguard  rr
Default = 1.0e−12=1.0e−12
r = εsafer=εsafe. If εqεq is the estimated error from the quadrature evaluation alone, and εexεex is the error estimate determined using extrapolation, then the extrapolated solution will only be accepted if εsafe εq εex εsafe εq εex.
Maximum Subdivisions  ii
Default = 50=50
i = maxsdivi=maxsdiv, the maximum number of subdivisions the algorithm may use in the adaptive phase, forming at most an additional (2 × maxsdiv) (2×maxsdiv)  segments.
Primary Divisions  ii
Default = 1=1
i = sprii=spri, the number of initial segments of the domain [a,b][a,b]. By default the initial segment is the entire domain.
Constraint: 0 < i < 10000000<i<1000000.
Primary Division Mode  aa
Default = AUTOMATIC=AUTOMATIC
Determines how the initial set of sprispri segments will be generated.
AUTOMATIC
nag_quad_1d_gen_vec_multi_rcomm (d01ra) will automatically generate sprispri segments of equal size covering the interval [a,b][a,b].
MANUAL
nag_quad_1d_gen_vec_multi_rcomm (d01ra) will use the break points xi0xi0, for i = 1,2,,spri1i=1,2,,spri-1, supplied in x on initial entry to generate the initial segments covering [a,b][a,b]. These may be supplied in any order, however it will be more efficient to supply them in ascending (or descending if a > ba>b) order. Repeated break points are allowed, although this will generate fewer initial segments.
Note: an absolute bound on the size of an initial segment of 10.0ε10.0ε is automatically applied in all cases, and will result in fewer initial subdivisions being generated if automatically generated or supplied breakpoints result in segments smaller than this..
Prioritize Error  aa
Default = LEVEL=LEVEL
Indicates how new subdivisions of segments sustaining unacceptable local errors for integrals should be prioritized.
LEVELLEVEL
Segments with lower level with unsatisfactory error estimates will be chosen over segments with greater error on higher levels. This will probably lead to more integrals being improved in earlier iterations of the algorithm, and hence will probably lead to fewer repeated returns (see parameter sid), and to more integrals being satisfactorily estimated if computational exhaustion occurs.
MAXERRMAXERR
The segment with the worst overall error will be split, regardless of level. This will more rapidly improve the worst integral estimates, although it will probably result in the fewest integrals being improved in earlier iterations, and may hence lead to more repeated returns (see parameter sid), and potentially fewer integrals satisfying the requested tolerances if computational exhaustion occurs.
Quadrature Rule  aa
Default = GK15=GK15
The basic quadrature rule to be used during the integration. Currently 66 Gauss–Kronrod rules are available, all identifiable by the letters GK followed by the number of points required by the Kronrod rule. Higher order rules generally provide higher accuracy with fewer subdivisons. However, for integrands with sharp singularities, lower order rules may be more efficient, particularly if the integrand away from the singularity is well behaved. With higher order rules, you may need to increase the Absolute Interval Minimum and the Relative Interval Minimum to maintain numerical difference between the abscissae and the segment bounds.
GK15
The Gauss–Kronrod rule based on 77 Gauss points and 1515 Kronrod points.
GK21
The Gauss–Kronrod rule based on 1010 Gauss points and 2121 Kronrod points. This is the rule used by nag_quad_1d_fin_bad_vec (d01at)
GK31
The Gauss–Kronrod rule based on 1515 Gauss points and 3131 Kronrod points.
GK41
The Gauss–Kronrod rule based on 2020 Gauss points and 4141 Kronrod points.
GK51
The Gauss–Kronrod rule based on 2525 Gauss points and 5151 Kronrod points.
GK61
The Gauss–Kronrod rule based on 3030 Gauss points and 6161 Kronrod points. This is the highest order rule, most suitable for highly oscilliatory integrals.
Relative Interval Minimum  rr
Default = 1.0e−6=1.0e−6
r = δrr=δr, the relative factor in the lower limit, δr|ba|δr|b-a|, for a segment to be considered for subdivision. See also Absolute Interval Minimum and Section [Further Comments].
Constraint: r0.0r0.0.
Relative Tolerance  rr
Default = sqrt(ε)=ε 
r = εrr=εr, the required relative tolerance. See also Absolute Tolerance and Section [Description].
Constraint: r0.0r0.0.
Note: setting both εr = εa = 0.0εr=εa=0.0 is possible, although it will most likely result in an excessive amount of computational effort.

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–2013