# NAG Library Function Document

## 1Purpose

nag_quad_2d_fin (d01dac) attempts to evaluate a double integral to a specified absolute accuracy by repeated applications of the method described by Patterson (1968) and Patterson (1969).

## 2Specification

 #include #include
void  nag_quad_2d_fin (double ya, double yb,
 double (*phi1)(double y, Nag_Comm *comm),
 double (*phi2)(double y, Nag_Comm *comm),
 double (*f)(double x, double y, Nag_Comm *comm),
double absacc, double *ans, Integer *npts, Nag_Comm *comm, NagError *fail)

## 3Description

nag_quad_2d_fin (d01dac) attempts to evaluate a definite integral of the form
 $I= ∫ab ∫ ϕ1y ϕ2y fx,y dx dy$
where $a$ and $b$ are constants and ${\varphi }_{1}\left(y\right)$ and ${\varphi }_{2}\left(y\right)$ are functions of the variable $y$.
The integral is evaluated by expressing it as
 $I=∫abFydy, where Fy= ∫ ϕ1y ϕ2y fx,ydx.$
Both the outer integral $I$ and the inner integrals $F\left(y\right)$ are evaluated by the method, described by Patterson (1968) and Patterson (1969), of the optimum addition of points to Gauss quadrature formulae.
This method uses a family of interlacing common point formulae. Beginning with the $3$-point Gauss rule, formulae using $7$, $15$, $31$, $63$, $127$ and finally $255$ points are derived. Each new formula contains all the points of the earlier formulae so that no function evaluations are wasted. Each integral is evaluated by applying these formulae successively until two results are obtained which differ by less than the specified absolute accuracy.

## 4References

Patterson T N L (1968) On some Gauss and Lobatto based integration formulae Math. Comput. 22 877–881
Patterson T N L (1969) The optimum addition of points to quadrature formulae, errata Math. Comput. 23 892

## 5Arguments

1:    $\mathbf{ya}$doubleInput
On entry: $a$, the lower limit of the integral.
2:    $\mathbf{yb}$doubleInput
On entry: $b$, the upper limit of the integral. It is not necessary that $a.
3:    $\mathbf{phi1}$function, supplied by the userExternal Function
phi1 must return the lower limit of the inner integral for a given value of $y$.
The specification of phi1 is:
 double phi1 (double y, Nag_Comm *comm)
1:    $\mathbf{y}$doubleInput
On entry: the value of $y$ for which the lower limit must be evaluated.
2:    $\mathbf{comm}$Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to phi1.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling nag_quad_2d_fin (d01dac) you may allocate memory and initialize these pointers with various quantities for use by phi1 when called from nag_quad_2d_fin (d01dac) (see Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: phi1 should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_quad_2d_fin (d01dac). If your code inadvertently does return any NaNs or infinities, nag_quad_2d_fin (d01dac) is likely to produce unexpected results.
4:    $\mathbf{phi2}$function, supplied by the userExternal Function
phi2 must return the upper limit of the inner integral for a given value of $y$.
The specification of phi2 is:
 double phi2 (double y, Nag_Comm *comm)
1:    $\mathbf{y}$doubleInput
On entry: the value of $y$ for which the upper limit must be evaluated.
2:    $\mathbf{comm}$Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to phi2.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling nag_quad_2d_fin (d01dac) you may allocate memory and initialize these pointers with various quantities for use by phi2 when called from nag_quad_2d_fin (d01dac) (see Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: phi2 should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_quad_2d_fin (d01dac). If your code inadvertently does return any NaNs or infinities, nag_quad_2d_fin (d01dac) is likely to produce unexpected results.
5:    $\mathbf{f}$function, supplied by the userExternal Function
f must return the value of the integrand $f$ at a given point.
The specification of f is:
 double f (double x, double y, Nag_Comm *comm)
1:    $\mathbf{x}$doubleInput
2:    $\mathbf{y}$doubleInput
On entry: the coordinates of the point $\left(x,y\right)$ at which the integrand $f$ must be evaluated.
3:    $\mathbf{comm}$Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to f.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling nag_quad_2d_fin (d01dac) you may allocate memory and initialize these pointers with various quantities for use by f when called from nag_quad_2d_fin (d01dac) (see Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: f should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_quad_2d_fin (d01dac). If your code inadvertently does return any NaNs or infinities, nag_quad_2d_fin (d01dac) is likely to produce unexpected results.
6:    $\mathbf{absacc}$doubleInput
On entry: the absolute accuracy requested.
7:    $\mathbf{ans}$double *Output
On exit: the estimated value of the integral.
8:    $\mathbf{npts}$Integer *Output
On exit: the total number of function evaluations.
9:    $\mathbf{comm}$Nag_Comm *
The NAG communication argument (see Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
10:  $\mathbf{fail}$NagError *Input/Output
The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

## 6Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 2.3.1.2 in How to Use the NAG Library and its Documentation for further information.
On entry, argument $〈\mathit{\text{value}}〉$ had an illegal value.
NE_CONVERGENCE
The outer integral has converged, but $n$ of the inner integrals have not converged with $255$ points: $n=〈\mathit{\text{value}}〉$. ans may still contain an approximate estimate of the integral, but its reliability will decrease as $n$ increases.
The outer integral has not converged, and $n$ of the inner integrals have not converged with $255$ points: $n=〈\mathit{\text{value}}〉$. ans may still contain an approximate estimate of the integral, but its reliability will decrease as $n$ increases.
The outer integral has not converged with $255$ points. However, all the inner integrals have converged, and ans may still contain an approximate estimate of the integral.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.

## 7Accuracy

The absolute accuracy is specified by the variable absacc. If, on exit, ${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_NOERROR then the result is most likely correct to this accuracy. Even if ${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_CONVERGENCE on exit, it is still possible that the calculated result could differ from the true value by less than the given accuracy.

## 8Parallelism and Performance

Please consult the x06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

The time taken by nag_quad_2d_fin (d01dac) depends upon the complexity of the integrand and the accuracy requested.
With Patterson's method accidental convergence may occasionally occur, when two estimates of an integral agree to within the requested accuracy, but both estimates differ considerably from the true result. This could occur in either the outer integral or in one or more of the inner integrals.
If it occurs in the outer integral then apparent convergence is likely to be obtained with considerably fewer integrand evaluations than may be expected. If it occurs in an inner integral, the incorrect value could make the function $F\left(y\right)$ appear to be badly behaved, in which case a very large number of pivots may be needed for the overall evaluation of the integral. Thus both unexpectedly small and unexpectedly large numbers of integrand evaluations should be considered as indicating possible trouble. If accidental convergence is suspected, the integral may be recomputed, requesting better accuracy; if the new request is more stringent than the degree of accidental agreement (which is of course unknown), improved results should be obtained. This is only possible when the accidental agreement is not better than machine accuracy. It should be noted that the function requests the same accuracy for the inner integrals as for the outer integral. In practice it has been found that in the vast majority of cases this has proved to be adequate for the overall result of the double integral to be accurate to within the specified value.
The function is not well-suited to non-smooth integrands, i.e., integrands having some kind of analytic discontinuity (such as a discontinuous or infinite partial derivative of some low-order) in, on the boundary of, or near, the region of integration. Warning: such singularities may be induced by incautiously presenting an apparently smooth interval over the positive quadrant of the unit circle, $R$
 $I=∫Rx+ydx dy.$
This may be presented to nag_quad_2d_fin (d01dac) as
 $I=∫01 dy ∫01-y2 x+ydx=∫01 121-y2+y⁢1-y2 dy$
but here the outer integral has an induced square-root singularity stemming from the way the region has been presented to nag_quad_2d_fin (d01dac). This situation should be avoided by re-casting the problem. For the example given, the use of polar coordinates would avoid the difficulty:
 $I=∫01dr∫0π2r2cos⁡υ+sin⁡υdυ.$

## 10Example

This example evaluates the integral discussed in Section 9, presenting it to nag_quad_2d_fin (d01dac) first as
 $∫01 ∫01-y2x+y dx dy$
and then as
 $∫01∫0π2r2cos⁡υ+sin⁡υdυ dr.$
Note the difference in the number of function evaluations.

### 10.1Program Text

Program Text (d01dace.c)

None.

### 10.3Program Results

Program Results (d01dace.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017