naginterfaces.library.quad.dim1_​gen_​vec_​multi_​rcomm

naginterfaces.library.quad.dim1_gen_vec_multi_rcomm(irevcm, a, b, needi, x, nx, fm, dinest, errest, comm)

dim1_gen_vec_multi_rcomm 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\left(x\right)$.

$$F={\int }_a^{b}f\left(x\right)dx$$

If the same subdivisions of the range are equally good for functions $f_1\left(x\right)$ and $f_2\left(x\right)$, because $f_1\left(x\right)$ and $f_2\left(x\right)$ have common areas of the range where they vary slowly and where they vary quickly, then we say that $f_1\left(x\right)$ and $f_2\left(x\right)$ are ‘similar’. dim1_gen_vec_multi_rcomm is particularly effective for the integration of a vector of similar functions.

For full information please refer to the NAG Library document for d01ra

https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/d01/d01raf.html

Parameters

irevcmint

On initial entry: $irevcm=1$.

$irevcm=1$

Sets up data structures in $comm$[‘icom’] and $comm$[‘com’] and starts a new integration.

On intermediate entry: $irevcm$ should normally be left unchanged. However, if $irevcm$ is set to a negative value, dim1_gen_vec_multi_rcomm will terminate, (see $irevcm=11$ and $irevcm=12$ above).

afloat

$a$, the lower bound of the domain.

bfloat

$b$, the upper bound of the domain.

If $|b-a|<10ϵ$, where $ϵ$ is the machine precision (see machine.precision), dim1_gen_vec_multi_rcomm will return $dinest[\text{j}-1]=errest[\text{j}-1]=0.0$, for $\text{j}=1,2,\dots ,n_i$.

neediint, ndarray, shape $\left(\text{ni}\right)$, modified in place

On initial entry: need not be set.

On intermediate exit: $needi[j-1]$ indicates what action must be taken for integral $j=1,2,\dots n_i$ (see $irevcm$).

$needi[j-1]=0$

Do not provide $f_j\left(x_i\right)$. Any provided values will be ignored.

$needi[j-1]=1$

The values $f_j\left(x_{\text{i}}\right)$ must be provided in $fm[j-1,\text{i}-1]$, for $\text{i}=1,2,\dots ,n_x$.

$needi[j-1]=2$

The values $f_j\left(x_i\right)$ are not required, however the error estimate for integral $j$ is still above the requested tolerance. If you wish to provide values for the evaluation of integral $j$, set $needi[j-1]=1$, and supply $f_j\left(x_{\text{i}}\right)$ in $fm[j-1,\text{i}-1]$, for $\text{i}=1,2,\dots ,n_x$.

$needi[j-1]=3$

The error estimate for integral $j$ 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 $f_j\left(x_i\right)$ 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 $j$, set $needi[j-1]=1$, and supply $f_j\left(x_{\text{i}}\right)$ in $fm[j-1,\text{i}-1]$, for $\text{i}=1,2,\dots ,n_x$.

$needi[j-1]=4$

The error estimate of integral $j$ is below the requested tolerance. If you believe this to be false, if for example the result in $dinest[j-1]$ is greatly different to what you may expect, you may force the algorithm to re-evaluate this conclusion by including the evaluations of integrand $j$ at $x_{\text{i}}$, for $\text{i}=1,2,\dots ,n_x$, and setting $needi[j-1]=1$. Integral and error estimation will be performed again during the next iteration.

On intermediate entry: $needi[j-1]$ may be used to indicate what action you have taken for integral $j$.

$needi[j-1]=1$

You have provided the values $f_j\left(x_{\text{i}}\right)$ in $fm[j-1,\text{i}-1]$, for $\text{i}=1,2,\dots ,n_x$.

$needi[j-1]<0$

You are abandoning the evaluation of integral $j$. The current values of $dinest[j-1]$ and $errest[j-1]$ will be returned on final completion.

Otherwise you have not provided the value $f_j\left(x_i\right)$.

On final exit: $needi[j-1]$ indicates the final state of integral $j$.

$needi[j-1]=0$

The error estimate for $F_j$ is below the requested tolerance.

$needi[j-1]=1$

The error estimate for $F_j$ is below the requested tolerance after extrapolation.

$needi[j-1]=2$

The error estimate for $F_j$ is above the requested tolerance.

$needi[j-1]=3$

The error estimate for $F_j$ is above the requested tolerance, and extremely bad behaviour of integral $j$ has been detected.

$needi[j-1]<0$

You prohibited further evaluation of integral $j$.

xfloat, ndarray, shape $\left(\text{lenxrq}\right)$, modified in place

On initial entry: if $\text{‘Primary Division Mode'}=\text{'AUTOMATIC'}$, $x$ need not be set. This is the default behaviour.

If $\text{‘Primary Division Mode'}=\text{'MANUAL'}$, $x$ is used to supply a set of initial ‘break-points’ inside the domain of integration.

Specifically, $x[i-1]$ must contain a break-point $x_{\text{i}}^0$, for $\text{i}=1,2,\dots ,(s_{\text{pri}}-1)$, where $s_{\text{pri}}$ is the number of ‘Primary Divisions’.

On intermediate exit: $x[\text{i}-1]$ is the abscissa $x_{\text{i}}$, for $\text{i}=1,2,\dots ,n_x$, at which the appropriate integrals must be evaluated.

nxint

On initial entry: need not be set.

On intermediate entry: must not be changed.

fmfloat, array-like, shape $(\text{ldfmrq},\text{sdfmrq})$

On initial entry: need not be set.

On intermediate entry: if indicated by $needi[j-1]$ you must supply the values $f_j\left(x_i\right)$ in $fm[\text{j}-1,\text{i}-1]$, for $\text{j}=1,2,\dots ,n_i$, for $\text{i}=1,2,\dots ,n_x$.

dinestfloat, ndarray, shape $\left(\text{ni}\right)$, modified in place

$dinest[j-1]$ contains the current estimate of the definite integral $F_j$.

On initial entry: need not be set.

On intermediate entry: must not be altered.

On exit: contains the current estimates of the $\text{ni}$ integrals. If $irevcm=0$, this will be the final solution.

errestfloat, ndarray, shape $\left(\text{ni}\right)$, modified in place

$errest[j-1]$ contains the current error estimate of the definite integral $F_j$.

On initial entry: need not be set.

On intermediate entry: must not be altered.

On exit: contains the current error estimates for the $\text{ni}$ integrals. If $irevcm=0$, $errest$ contains the final error estimates of the $\text{ni}$ integrals.

commdict, communication object, modified in place

Communication structure.

This argument must have been initialized by a prior call to opt_set().

Returns

irevcmint

On intermediate exit: $irevcm=11$ or $12$.

$irevcm$ requests the integrands $f_j\left(x_i\right)$ be evaluated for all required $j\in 1,\dots ,n_i$ as indicated by $needi$, and at all the points $x_{\text{i}}$, for $\text{i}=1,2,\dots ,n_x$.

Abscissae $x_i$ are provided in $x[i-1]$ and $f_j\left(x_i\right)$ must be returned in $fm[j-1,i-1]$.

During the initial solve phase:

$irevcm=11$

Function values are required to construct the initial estimates of the definite integrals.

If $needi[j-1]=1$, $f_j\left(x_i\right)$ must be supplied in $fm[j-1,i-1]$.

This will be the case unless you have abandoned the evaluation of specific integrals on a previous call.

If $needi[j-1]=0$, you have previously abandoned the evaluation of integral $j$, and hence should not supply the value of $f_j\left(x_i\right)$.

$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[\text{j}-1]$, for $\text{j}=1,2,\dots ,n_i$, will be set to this negative value and $errno$ = -1 will be returned.

During the adaptive solve phase:

$irevcm=12$

Function values are required to improve the estimates of the definite integrals.

If $needi[j-1]=0$, any evaluation of $f_j\left(x_i\right)$ will be discarded, so there is no need to provide them.

If $needi[j-1]=1$, $f_j\left(x_i\right)$ must be provided in $fm[j-1,i-1]$.

If $needi[j-1]=2$, $3$ or $4$, the current error estimate of integral $j$ does not require integrand $j$ to be evaluated and provided in $fm[j-1,i-1]$.

Should you choose to, integrand $j$ can be evaluated in which case $needi[j-1]$ must be set to $1$.

$dinest$ and $errest$ contain complete information during this phase.

If $irevcm$ is set to a negative value during this phase $errno$ = 1, 2 or 3 will be returned and the elements of $needi$ will reflect the current state of the adaptive process.

On final exit: $irevcm=0$.

$irevcm=0$

Indicates the algorithm has completed.

sidint

For advanced users.

On intermediate exit: $sid$ identifies a specific set of abscissae, $x$, returned during the integration process. When a new set of abscissae are generated the value of $sid$ is incremented by $1$. Advanced users may store calculations required for an identified set $x$, and reuse them should dim1_gen_vec_multi_rcomm return the same value of $sid$, i.e., the same set of abscissae was used.

nxint

On intermediate exit: $n_x$, the number of abscissae at which integrands are required.

Other Parameters

‘Absolute Interval Minimum’float

Default $=128.0ϵ$

$r={\delta }_a$, the absolute lower limit for a segment to be considered for subdivision. See also ‘Relative Interval Minimum’ and Further Comments.

Constraint: $r\ge 128ϵ$.

‘Absolute Tolerance’float

Default $=1024ϵ$

$r={ϵ}_a$, the absolute tolerance required. See also ‘Relative Tolerance’ and Notes.

Constraint: $r\ge 0.0$.

‘Extrapolation’str

Default $=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.

$ON$

Use extrapolation.

$OFF$

Disable extrapolation.

‘Extrapolation Safeguard’float

Default $=1.0e-12$

$r={ϵ}_{\text{safe}}$. If ${ϵ}_q$ is the estimated error from the quadrature evaluation alone, and ${ϵ}_{\text{ex}}$ is the error estimate determined using extrapolation, then the extrapolated solution will only be accepted if ${ϵ}_{\text{safe}}{ϵ}_q\le {ϵ}_{\text{ex}}$.

‘Maximum Subdivisions’int

Default $=50$

$i={max}_{\text{sdiv}}$, the maximum number of subdivisions the algorithm may use in the adaptive phase, forming at most an additional $(2\times {max}_{\text{sdiv}})$ segments.

‘Primary Divisions’int

Default $=1$

$i=s_{\text{pri}}$, the number of initial segments of the domain $[a,b]$. By default the initial segment is the entire domain.

Constraint: $0<i<1000000$.

‘Primary Division Mode’str

Default $=\text{'AUTOMATIC'}$

Determines how the initial set of $s_{\text{pri}}$ segments will be generated.

AUTOMATIC

dim1_gen_vec_multi_rcomm will automatically generate $s_{\text{pri}}$ segments of equal size covering the interval $[a,b]$.

MANUAL

dim1_gen_vec_multi_rcomm will use the break-points $x_{\text{i}}^0$, for $\text{i}=1,2,\dots ,s_{\text{pri}}-1$, supplied in $x$ on initial entry to generate the initial segments covering $[a,b]$. These may be supplied in any order, however it will be more efficient to supply them in ascending (or descending if $a>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ϵ$ 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.

‘Prioritize Error’str

Default $=\text{'LEVEL'}$

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

$\text{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.

$\text{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.

‘Quadrature Rule’str

Default $=\text{GK15}$

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 dim1_fin_general().

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.

‘Relative Interval Minimum’float

Default $=1.0e-6$

$r={\delta }_r$, the relative factor in the lower limit, ${\delta }_r|b-a|$, for a segment to be considered for subdivision. See also ‘Absolute Interval Minimum’ and Further Comments.

Constraint: $r\ge 0.0$.

‘Relative Tolerance’float

Default $=\sqrt{ϵ}$

$r={ϵ}_r$, the required relative tolerance. See also ‘Absolute Tolerance’ and Notes.

Constraint: $r\ge 0.0$.

Note: setting both ${ϵ}_r={ϵ}_a=0.0$ is possible, although it will most likely result in an excessive amount of computational effort.

‘Index LDI’int

query only

$I_{\text{ldi}}$, the index of $comm$[‘icom’] required for obtaining $\text{ldi}$. See Details of the Computation.

‘Index LDR’int

query only

$I_{\text{ldr}}$, the index of $comm$[‘icom’] required for obtaining $\text{ldr}$. See Details of the Computation.

‘Index NSDIV’int

query only

$I_{\text{nsdiv}}$, the index of $comm$[‘icom’] required for obtaining $\text{nsdiv}$. See Details of the Computation.

‘Index NSEG’int

query only

$I_{\text{nseg}}$, the index of $comm$[‘icom’] required for obtaining $\text{nseg}$. See Details of the Computation.

‘Index FCP’int

query only

$I_{\text{fcp}}$, the index of $comm$[‘icom’] required for obtaining $\text{fcp}$. See Details of the Computation.

‘Index EVALSP’int

query only

$I_{\text{evalsp}}$, the index of $comm$[‘icom’] required for obtaining $\text{evalsp}$. See Details of the Computation.

‘Index DSP’int

query only

$I_{\text{dsp}}$, the index of $comm$[‘icom’] required for obtaining $\text{dsp}$. See Details of the Computation.

‘Index ESP’int

query only

$I_{\text{esp}}$, the index of $comm$[‘icom’] required for obtaining $\text{esp}$. See Details of the Computation.

‘Index SINFOIP’int

query only

$I_{\text{sinfoip}}$, the index of $comm$[‘icom’] required for obtaining $\text{sinfoip}$. See Details of the Computation.

‘Index SINFORP’int

query only

$I_{\text{sinforp}}$, the index of $comm$[‘icom’] required for obtaining $\text{sinforp}$. See Details of the Computation.

Raises

NagValueError

(errno $11$)

$irevcm$ had an illegal value.

On entry, $irevcm=⟨value⟩$.

(errno $21$)

On entry, $\text{ni}=⟨value⟩$.

Constraint: $\text{ni}\ge 1$.

(errno $71$)

On entry, $\text{‘Primary Division Mode'}=\text{'MANUAL'}$ and at least one supplied break-point in $x$ is outside of the domain of integration.

(errno $81$)

$\text{lenx}$ is insufficient for the chosen options.

On entry, $\text{lenx}=⟨value⟩$.

Constraint: $\text{lenx}\ge ⟨value⟩$.

(errno $1001$)

Either the option arrays $comm$[‘iopts’] and $comm$[‘opts’] have not been initialized for dim1_gen_vec_multi_rcomm, or they have become corrupted.

(errno $1101$)

On entry, one of $comm$[‘icom’] and $comm$[‘com’] has become corrupted.

Warns

NagAlgorithmicWarning

(errno $-1$)

Evaluation of all integrals has been stopped during the initial phase.

(errno $1$)

At least one error estimate exceeded the requested tolerances.

(errno $2$)

Extremely bad behaviour was detected for at least one integral.

(errno $3$)

Extremely bad behaviour was detected for at least one integral. At least one other integral error estimate was above the requested tolerance.

Notes

dim1_gen_vec_multi_rcomm 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 dim1_gen_vec_multi_rcomm 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).

dim1_gen_vec_multi_rcomm is the integration function in the suite of functions dim1_gen_vec_multi_rcomm and dim1_gen_vec_multi_dimreq(). It also uses options, which can be set and queried using the functions opt_set() and opt_get() respectively. The options available are described in Other Parameters.

The algorithm used by dim1_gen_vec_multi_rcomm is divided into an initial solve phase and an adaptive solve phase.

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_{\text{pri}}$ level $1$ segments, where $s_{\text{pri}}$ is the number of ‘Primary Divisions’, and will use any provided break-points if $\text{‘Primary Division Mode'}=\text{'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}_{\text{sdiv}}$ subdivisions are allowed if sufficient memory is provided, where ${max}_{\text{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 ${ϵ}_{\text{safe}}$ can be set using ‘Extrapolation Safeguard’, and the extrapolated solution will only be considered if ${ϵ}_{\text{safe}}{ϵ}_q\le {ϵ}_{\text{ex}}$, where ${ϵ}_q$ and ${ϵ}_{\text{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-1]$ on completion of dim1_gen_vec_multi_rcomm. Otherwise the direct solution will be returned in $dinest[j-1]$. This is indicated by the value of $needi[j-1]$ on completion.

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 $e_m\left(S_n\right)$ transformation, Math. Tables Aids Comput. (10), 91–96