d01ra:: Quadrature (NAG Toolbox)

nag_quad_1d_gen_vec_multi_rcomm (d01ra) is a general purpose adaptive integrator which calculates an approximation to a vector of definite integrals

F

over a finite range

[a, b]

, given the vector of integrands

f (x)

.

F = \int_{a}^{b} f (x) d x

If the same subdivisions of the range are equally good for functions

f_{1} (x)

and

f_{2} (x)

, because

f_{1} (x)

and

f_{2} (x)

have common areas of the range where they vary slowly and where they vary quickly, then we say that

f_{1} (x)

and

f_{2} (x)

are ‘similar’. nag_quad_1d_gen_vec_multi_rcomm (d01ra) is particularly effective for the integration of a vector of similar functions.

nag_quad_1d_gen_vec_multi_rcomm (d01ra) is an extension to various QUADPACK routines, including QAG, QAGS and QAGP. 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 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).

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]

. This will have been divided into

s_{pri}

level

1

segments, where

s_{pri}

is the number of Primary Divisions, and will use any provided break-points if

Primary Division Mode = MANUAL

.

Once a complete integral estimate over

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

and

ε_{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

\max_{sdiv}

subdivisions are allowed if sufficient memory is provided, where

\max_{sdiv}

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}

can be set using Extrapolation Safeguard, and the extrapolated solution will only be considered if

ε_{safe} ε_{q} \leq ε_{ex}

, where

ε_{q}

and

ε_{ex}

are the estimated error directly from the quadrature and from the extrapolation respectively. If extrapolation is successful for the computation of integral

j

, the extrapolated solution will be returned in

dinest (j)

on completion of nag_quad_1d_gen_vec_multi_rcomm (d01ra). Otherwise the direct solution will be returned in

dinest (j)

. This is indicated by the value of

needi (j)

on completion.

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

e_{m} (S_{n})

transformation Math. Tables Aids Comput. 10 91–96

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 argument irevcm. Between intermediate exits and re-entries, all arguments other than irevcm, needi and fm must remain unchanged.

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:

nag_quad_1d_gen_vec_multi_rcomm (d01ra) cannot guarantee, but in practice usually achieves, the following accuracy for each integral

F_{j}

:

|F_{j} - dinest (j)| \leq tol

where

tol = \max (ε_{a}, ε_{r} \times |F_{j}|)

ε_{a}

and

ε_{r}

are the error tolerances Absolute Tolerance and Relative Tolerance respectively. Moreover, it returns errest, the entries of which in normal circumstances satisfy,

|F_{j} - dinest (j)| \leq errest (j) \leq tol .

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

f_{j}

.

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,

\int_{0}^{1} (\begin{matrix} \log (x) \\ x^{- \frac{1}{2}} \\ x^{2} \end{matrix}) d x

will be quite efficient, as the irregular behaviour of the first two integrands is at

x = 0

. On the contrary, the evaluation of the integrals,

\int_{0}^{1} (\begin{matrix} \log (x) \\ \log (1 - x) \end{matrix}) d x

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

\bar{k}

with bounds

[a_{\bar{k}}, b_{\bar{k}}]

satisfying

|b_{\bar{k}} - a_{\bar{k}}| < \max (δ_{a}, δ_{r} \times |b - a|)

has a local error estimate greater than the requested tolerance for at least one integral. The values

δ_{a}

and

δ_{r}

can be set through the optional parameters Absolute Interval Minimum and Relative Interval Minimum respectively.

This section is recommended for expert users only. It describes the contents of the arrays com and icom upon exit from nag_quad_1d_gen_vec_multi_rcomm (d01ra) with

ifail = 0

,

1

,

2

or

3

, and provided at least one iteration completed, failure due to insufficient licom or lcom.

The arrays icom and com contain details of the integration, including various scalars, one-dimensional arrays, and (effectively) two-dimensional arrays. The dimensions of these arrays vary depending on the arguments and options used and the progress of the algorithm. Here we describe some of these details, including how and where they are stored in icom and com.

Scalar quantities:

The indices in icom including the following scalars are available via query only options, see Diagnostic Options. For example,

I_{ldi}

is the integer value returned by the option Index LDI.

$ldi$	The leading dimension of the two-dimensional integer arrays stored in icom detailed below. $ldi = icom (I_{ldi})$ .
$ldr$	The leading dimension of the two-dimensional real arrays stored in com detailed below. $ldr = icom (I_{ldr})$ .
$nsdiv$	The number of segments that have been subdivided during the adaptive process. $nsdiv = icom (I_{nsdiv})$ .
$nseg$	The total number of segments formed. $nseg = 2 nsdiv + s_{pri}$ . $nseg = icom (I_{nseg})$ .
$dsp$	The reference of the first element of the array $ds$ stored in com. $dsp = icom (I_{dsp})$ .
$esp$	The reference of the first element of the array $es$ stored in com. $esp = icom (I_{esp})$ .
$evalsp$	The reference of the first element of the array $evals$ stored in icom. $evalsp = icom (I_{evalsp})$ .
$fcp$	The reference of the first element of the array $fcount$ stored in icom. $fcp = icom (I_{fcp})$ .
$sinforp$	The reference of the first element of the array $sinfor$ stored in com. $sinforp = icom (I_{sinforp})$ .
$sinfoip$	The reference of the first element of the array $sinfoi$ stored in icom. $sinfoip = icom (I_{sinfoip})$ .

One-dimensional arrays:

$fcount (n_{i})$: $fcount (1) = icom (fcp)$ .
$fcount (j)$ contains the number of different approximations of integral $j$ calculated, for $j = 1, 2, \dots, n_{i}$ .

Two-dimensional arrays:

$sinfoi (5, nseg)$: $sinfoi (1, 1) = icom (sinfoip)$ .
$sinfoi$ contains information about the hierarchy of splitting.
$sinfoi (1, k)$ contains the split identifier for segment $k$ , for $k = 1, 2, \dots, nseg$ .

$sinfoi (2, k)$ contains the parent segment number of segment $k$ (i.e., the segment was split to create segment $k$ ), for $k = 1, 2, \dots, nseg$ .

$sinfoi (3, k)$ and $sinfoi (4, k)$ contain the segment numbers of the two child segments formed from segment $k$ , if segment $k$ has been split. If segment $k$ has not been split, these will be negative.

$sinfoi (5, k)$ contains the level at which the segment exists, corresponding to $n_{a} + 1$ , where $n_{a}$ is the number of ancestor segments of segment $k$ , for $k = 1, 2, \dots, nseg$ . A negative level indicates that segment $k$ will not be split further, the level is then given by the absolute value of $sinfoi (5, k)$ .
$sinfor (2, nseg)$: $sinfor (1, 1) = com (sinforp)$ .
$sinfor$ contains the bounds of each segment.
$sinfor (1, k)$ contains the lower bound of segment $k$ , for $k = 1, 2, \dots, nseg$ .

$sinfor (2, k)$ contains the upper bound of segment $k$ , for $k = 1, 2, \dots, nseg$ .
$evals (n_{i}, nseg)$: $evals (1, 1) = icom (evalsp)$ .
$evals$ contains information to indicate whether an estimate of the integral $j$ has been obtained over segment $k$ , and if so whether this evaluation still contributes to the direct estimate of $F_{j}$ , for $j = 1, 2, \dots, n_{i}$ and $k = 1, 2, \dots, nseg$ .

$evals (j, k) = 0$ indicates that integral $j$ has not been evaluated over segment $k$ .

$evals (j, k) = 1$ indicates that integral $j$ has been evaluated over segment $k$ , and that this evaluation contributes to the direct estimate of $F_{j}$ .

$evals (j, k) = 2$ indicates that integral $j$ has been evaluated over segment $k$ , that this evaluation contributes to the direct estimate of $F_{j}$ , and that you have requested no further evaluation of this integral at this segment by setting $needi (j) < 0$ .

$evals (j, k) = 3$ indicates that integral $j$ has been evaluated over segment $k$ , and this evaluation no longer contributes to the direct estimate of $F_{j}$ .

$evals (j, k) = 4$ indicates that integral $j$ has been evaluated over segment $k$ , that this evaluation contributes to the direct estimate of $F_{j}$ , and that this segment is too small for any further splitting to be performed. Integral $j$ also has a local error estimate over this segment above the requested tolerance. Such segments cause nag_quad_1d_gen_vec_multi_rcomm (d01ra) to return $ifail = 2$ or $3$ , indicating extremely bad behaviour.

$evals (j, k) = 5$ indicates that integral $j$ has been evaluated over segment $k$ , that this evaluation contributes to the direct estimate of $F_{j}$ , and that this segment is too small for any further splitting to be performed. The local error estimate is however below the requested tolerance.
$ds (n_{i}, nseg)$: $ds (1, 1) = com (dsp)$ .
$ds (j, k)$ contains the definite integral estimate of the $j$ th integral over the $k$ th segment, ${ds}_{j, k}$ , provided it has been evaluated, for $j = 1, 2, \dots, n_{i}$ and $k = 1, 2, \dots, nseg$ .
$es (n_{i}, nseg)$: $es (1, 1) = com (esp)$ .
$es (j, k)$ contains the definite integral error estimate of the $j$ th integral over the $k$ th segment, ${es}_{j, k}$ , provided it has been evaluated, for $j = 1, 2, \dots, n_{i}$ and $k = 1, 2, \dots, nseg$ .

For each integral

j

, the direct approximation

D_{j}

of

F_{j}

, and its error estimate

E_{j}

, may be constructed as,

\begin{matrix} F_{j} \approx D_{j} = \sum_{K_{j}} {ds}_{j, k}, \\ |F_{j} - D_{j}| \approx E_{j} = \sum_{K_{j}} {es}_{j, k}, \end{matrix}

where

K_{j}

is the set of all contributing segments,

K_{j} = \{k ∣ evals (j, k) = 1, ​ 2, ​ 4 ​ or ​ 5, ​ 1 \leq k \leq nseg\}

.

D_{j}

will have been returned in

dinest (j)

, unless extrapolation was successful, as indicated by

needi (j)

.

Similarly,

E_{j}

will have been returned in

errest (j)

unless extrapolation was successful, in which case the error estimate from the extrapolation will have been returned. If for a given integral

j

one or more contributing segments have unacceptable error estimates, it may be possible to improve the direct approximation by replacing the contributions from these segments with more accurate estimates should these be calculable by some means. Indeed for any segment

\bar{k} \in k

, with lower bound

a_{\bar{k}} = sinfor (1, \bar{k})

and upper bound

b_{\bar{k}} = sinfor (2, \bar{k})

, one may alter the direct approximation

D_{j}

by the following,

\begin{matrix} {ds}_{j, \bar{k}}^{new} \approx \int_{a_{\bar{k}}}^{b_{\bar{k}}} f_{j} (x) ​ ​ dx \\ D_{j} = \sum_{K_{j}} {ds}_{j, k} - {ds}_{j, \bar{k}} + {ds}_{j, \bar{k}}^{new} . \end{matrix}

The error estimate

E_{j}

may be altered similarly.

This example integrates

F = \int_{0}^{π} (\begin{matrix} x \sin (2 x) \cos (15 x) \\ x^{2} \sin (2 x) \cos (50 x) \end{matrix}) d x .

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 Description of the s.

The following optional parameters, see Diagnostic Options, may be utilized by expert users in conjunction with the information provided in Details of the Computation.

Index LDI
Index LDR
Index NSDIV
Index NSEG
Index FCP
Index EVALSP
Index DSP
Index ESP
Index SINFOIP
Index SINFORP

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:

$ε$ represents the machine precision (see nag_machine_precision (x02aj)).

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.

Unsetable options will return the appropriate value when calling nag_quad_opt_get (d01zl). They will have no effect if passed to nag_quad_opt_set (d01zk).

For nag_quad_1d_gen_vec_multi_rcomm (d01ra) the maximum length of the argument cvalue used by nag_quad_opt_get (d01zl) is

15

.

r = δ_{a}

, the absolute lower limit for a segment to be considered for subdivision. See also Relative Interval Minimum and Further Comments.

Constraint:

r \geq 128 ε

.

r = ε_{a}

, the absolute tolerance required. See also Relative Tolerance and Description.

Constraint:

r \geq 0.0

.

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.

$ON$: Use extrapolation.
$OFF$: Disable extrapolation.

r = ε_{safe}

. If

ε_{q}

is the estimated error from the quadrature evaluation alone, and

ε_{ex}

is the error estimate determined using extrapolation, then the extrapolated solution will only be accepted if

ε_{safe} ε_{q} \leq ε_{ex}

.

i = \max_{sdiv}

, the maximum number of subdivisions the algorithm may use in the adaptive phase, forming at most an additional

(2 \times \max_{sdiv})

segments.

i = s_{pri}

, the number of initial segments of the domain

[a, b]

. By default the initial segment is the entire domain.

Constraint:

0 < i < 1000000

.

Determines how the initial set of

s_{pri}

segments will be generated.

Note: an absolute bound on the size of an initial segment of

10.0 ε

is automatically applied in all cases, and will result in fewer initial subdivisions being generated if automatically generated or supplied break-points result in segments smaller than this..

Indicates how new subdivisions of segments sustaining unacceptable local errors for integrals should be prioritized.

$LEVEL$: 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 argument sid), and to more integrals being satisfactorily estimated if computational exhaustion occurs.
$MAXERR$: 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 argument sid), and potentially fewer integrals satisfying the requested tolerances if computational exhaustion occurs.

The basic quadrature rule to be used during the integration. Currently

6

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 $7$ Gauss points and $15$ Kronrod points.
GK21: The Gauss–Kronrod rule based on $10$ Gauss points and $21$ Kronrod points. This is the rule used by nag_quad_1d_fin_bad_vec (d01at)
GK31: The Gauss–Kronrod rule based on $15$ Gauss points and $31$ Kronrod points.
GK41: The Gauss–Kronrod rule based on $20$ Gauss points and $41$ Kronrod points.
GK51: The Gauss–Kronrod rule based on $25$ Gauss points and $51$ Kronrod points.
GK61: The Gauss–Kronrod rule based on $30$ Gauss points and $61$ Kronrod points. This is the highest order rule, most suitable for highly oscilliatory integrals.

r = δ_{r}

, the relative factor in the lower limit,

δ_{r} |b - a|

, for a segment to be considered for subdivision. See also Absolute Interval Minimum and Further Comments.

Constraint:

r \geq 0.0

.

r = ε_{r}

, the required relative tolerance. See also Absolute Tolerance and Description.

Constraint:

r \geq 0.0

.

Note: setting both

ε_{r} = ε_{a} = 0.0

is possible, although it will most likely result in an excessive amount of computational effort.

These options are provided for expert users who wish to examine and modify the precise details of the computation. They should only be used after nag_quad_1d_gen_vec_multi_rcomm (d01ra) returns, as opposed to the options listed in Description of the s which must be used before the first call to nag_quad_1d_gen_vec_multi_rcomm (d01ra).