NAG Library Routine Document
d01daf (dim2_fin)
1
Purpose
d01daf attempts to evaluate a double integral to a specified absolute accuracy by repeated applications of the method described by
Patterson (1968) and
Patterson (1969).
2
Specification
Fortran Interface
Integer, Intent (Inout)  ::  ifail  Integer, Intent (Out)  ::  npts  Real (Kind=nag_wp), External  ::  phi1, phi2, f  Real (Kind=nag_wp), Intent (In)  ::  ya, yb, absacc  Real (Kind=nag_wp), Intent (Out)  ::  ans 

C Header Interface
#include <nagmk26.h>
void 
d01daf_ (const double *ya, const double *yb, double (NAG_CALL *phi1)(const double *y), double (NAG_CALL *phi2)(const double *y), double (NAG_CALL *f)(const double *x, const double *y), const double *absacc, double *ans, Integer *npts, Integer *ifail) 

3
Description
d01daf attempts to evaluate a definite integral of the form
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
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.
4
References
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
5
Arguments
 1: $\mathbf{ya}$ – Real (Kind=nag_wp)Input

On entry: $a$, the lower limit of the integral.
 2: $\mathbf{yb}$ – Real (Kind=nag_wp)Input

On entry: $b$, the upper limit of the integral. It is not necessary that $a<b$.
 3: $\mathbf{phi1}$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

phi1 must return the lower limit of the inner integral for a given value of
$y$.
The specification of
phi1 is:
Fortran Interface
Real (Kind=nag_wp)  ::  phi1  Real (Kind=nag_wp), Intent (In)  ::  y 

C Header Interface
#include <nagmk26.h>
double 
phi1 (const double *y) 

 1: $\mathbf{y}$ – Real (Kind=nag_wp)Input

On entry: the value of $y$ for which the lower limit must be evaluated.
phi1 must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d01daf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: phi1 should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d01daf. If your code inadvertently
does return any NaNs or infinities,
d01daf is likely to produce unexpected results.
 4: $\mathbf{phi2}$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

phi2 must return the upper limit of the inner integral for a given value of
$y$.
The specification of
phi2 is:
Fortran Interface
Real (Kind=nag_wp)  ::  phi2  Real (Kind=nag_wp), Intent (In)  ::  y 

C Header Interface
#include <nagmk26.h>
double 
phi2 (const double *y) 

 1: $\mathbf{y}$ – Real (Kind=nag_wp)Input

On entry: the value of $y$ for which the upper limit must be evaluated.
phi2 must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d01daf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: phi2 should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d01daf. If your code inadvertently
does return any NaNs or infinities,
d01daf is likely to produce unexpected results.
 5: $\mathbf{f}$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

f must return the value of the integrand
$f$ at a given point.
The specification of
f is:
Fortran Interface
Real (Kind=nag_wp)  ::  f  Real (Kind=nag_wp), Intent (In)  ::  x, y 

C Header Interface
#include <nagmk26.h>
double 
f (const double *x, const double *y) 

 1: $\mathbf{x}$ – Real (Kind=nag_wp)Input
 2: $\mathbf{y}$ – Real (Kind=nag_wp)Input

On entry: the coordinates of the point $\left(x,y\right)$ at which the integrand $f$ must be evaluated.
f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d01daf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: f should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d01daf. If your code inadvertently
does return any NaNs or infinities,
d01daf is likely to produce unexpected results.
 6: $\mathbf{absacc}$ – Real (Kind=nag_wp)Input

On entry: the absolute accuracy requested.
 7: $\mathbf{ans}$ – Real (Kind=nag_wp)Output

On exit: the estimated value of the integral.
 8: $\mathbf{npts}$ – IntegerOutput

On exit: the total number of function evaluations.
 9: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if
${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit, the recommended value is
$1$.
When the value $\mathbf{1}\text{or}1$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Note: d01daf may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}\mathrm{mod}10\ne 0$

The outer integral has converged, but
$n$ of the inner integrals have not converged with
$255$ points:
$n=\u2329\mathit{\text{value}}\u232a$.
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=\u2329\mathit{\text{value}}\u232a$.
ans may still contain an approximate estimate of the integral, but its reliability will decrease as
$n$ increases.
 ${\mathbf{ifail}}=1$

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.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
The absolute accuracy is specified by the variable
absacc. If, on exit,
${\mathbf{ifail}}={\mathbf{0}}$ then the result is most likely correct to this accuracy. Even if
ifail is nonzero
on exit, it is still possible that the calculated result could differ from the true value by less than the given accuracy.
8
Parallelism and Performance
d01daf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
The time taken by d01daf 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 routine 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 routine is not wellsuited to nonsmooth integrands, i.e., integrands having some kind of analytic discontinuity (such as a discontinuous or infinite partial derivative of some loworder) 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$
This may be presented to
d01daf as
but here the outer integral has an induced squareroot singularity stemming from the way the region has been presented to
d01daf. This situation should be avoided by recasting the problem. For the example given, the use of polar coordinates would avoid the difficulty:
10
Example
This example evaluates the integral discussed in
Section 9, presenting it to
d01daf first as
and then as
Note the difference in the number of function evaluations.
10.1
Program Text
Program Text (d01dafe.f90)
10.2
Program Data
None.
10.3
Program Results
Program Results (d01dafe.r)