NAG Library Function Document
nag_1d_withdraw_quad_gauss_1 (d01tac)
1
Purpose
nag_1d_withdraw_quad_gauss_1 (d01tac) computes an estimate of the definite integral of a function of known analytical form, using a Gaussian quadrature formula with a specified number of abscissae. Formulae are provided for a finite interval (Gauss–Legendre), a semiinfinite interval (Gauss–Laguerre, rational Gauss), and an infinite interval (Gauss–Hermite).
2
Specification
#include <nag.h> 
#include <nagd01.h> 
double 
nag_1d_withdraw_quad_gauss_1 (Nag_GaussFormulae quadrule,
double 
(*f)(double x,
Nag_User *comm),


double a,
double b,
Integer npts,
Nag_User *comm,
NagError *fail) 

3
Description
3.1
General
nag_1d_withdraw_quad_gauss_1 (d01tac) evaluates an estimate of the definite integral of a function
$f\left(x\right)$, over a finite or infinite interval, by
$n$point Gaussian quadrature (see
Davis and Rabinowitz (1967),
Fröberg (1970),
Ward (1975) or
Stroud and Secrest (1966)). The integral is approximated by a summation
where
${\omega}_{i}$ are called the weights, and the
${x}_{i}$ the abscissae. A selection of values of
$n$ is available. (See
Section 5.)
3.2
Both Limits Finite
The Gauss–Legendre weights and abscissae are used, and the formula is exact for any function of the form:
The formula is appropriate for functions which can be well approximated by such a polynomial over
$\left[a,b\right]$. It is inappropriate for functions with algebraic singularities at one or both ends of the interval, such as
$\frac{1}{\sqrt{\left(1+x\right)}}$ on
$\left[1,1\right]$.
3.3
One Limit Infinite
Two quadrature formulae are available for these integrals:
(a) 
The Gauss–Laguerre formula is exact for any function of the form:

(b) 
This formula is appropriate for functions decaying exponentially at infinity; the argument $b$ should be chosen if possible to match the decay rate of the function. 
(c) 
The rational Gauss formula is exact for any function of the form:

(d) 
This formula is likely to be more accurate for functions having only an inverse power rate of decay for large $x$. Here the choice of a suitable value of $b$ may be more difficult; unfortunately a poor choice of $b$ can make a large difference to the accuracy of the computed integral. 
3.4
Both Limits Infinite
The Gauss–Hermite weights and abscissae are used, and the formula is exact for any function of the form:
Again, for general functions not of this exact form, the argument
$b$ should be chosen to match if possible the decay rate at
$\pm \infty $.
4
References
Davis P J and Rabinowitz P (1967) Numerical Integration 33–52 Blaisdell Publishing Company
Fröberg C E (1970) Introduction to Numerical Analysis Addison–Wesley
Stroud A H and Secrest D (1966) Gaussian Quadrature Formulas Prentice–Hall
Ward R C (1975) The combination shift $QZ$ algorithm SIAM J. Numer. Anal. 12 835–853
5
Arguments
 1:
$\mathbf{quadrule}$ – Nag_GaussFormulaeInput

On entry: indicates the quadrature formula:
 ${\mathbf{quadrule}}=\mathrm{Nag\_Legendre}$, for Gauss–Legendre quadrature on a finite interval;
 ${\mathbf{quadrule}}=\mathrm{Nag\_Rational}$, for rational Gauss quadrature on a semiinfinite interval;
 ${\mathbf{quadrule}}=\mathrm{Nag\_Laguerre}$, for Gauss–Laguerre quadrature on a semiinfinite interval;
 ${\mathbf{quadrule}}=\mathrm{Nag\_Hermite}$, for Gauss–Hermite quadrature on an infinite interval.
Constraint:
${\mathbf{quadrule}}=\mathrm{Nag\_Legendre}$, $\mathrm{Nag\_Rational}$, $\mathrm{Nag\_Laguerre}$ or $\mathrm{Nag\_Hermite}$.
 2:
$\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,
Nag_User *comm)


 1:
$\mathbf{x}$ – doubleInput

On entry: the point at which the integrand $f$ must be evaluated.
 2:
$\mathbf{comm}$ – Nag_User *

Pointer to a structure of type Nag_User with the following member:
 p – Pointer

On entry/exit: the pointer
$\mathbf{comm}\mathbf{\to}\mathbf{p}$ should be cast to the required type, e.g.,
struct user *s = (struct user *)comm → p, to obtain the original object's address with appropriate type. (See the argument
comm below.)
Note: f should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
nag_1d_withdraw_quad_gauss_1 (d01tac). If your code inadvertently
does return any NaNs or infinities,
nag_1d_withdraw_quad_gauss_1 (d01tac) is likely to produce unexpected results.
Some points to bear in mind when coding
f are mentioned in
Section 9.
 3:
$\mathbf{a}$ – doubleInput
 4:
$\mathbf{b}$ – doubleInput

On entry: the arguments
$a$ and
$b$ which occur in the integration formulae:
Gauss–Legendre: $a$ is the lower limit and $b$ is the upper limit of the integral. It is not necessary that $a<b$.
Rational Gauss:
$b$ must be chosen so as to make the integrand match as closely as possible the exact form given in
Section 3. The interval of integration is
$\left[a,\infty \right)$ if
$a+b>0$, and
$\left(\infty ,a\right]$ if
$a+b<0$.
Gauss–Laguerre:
$b$ must be chosen so as to make the integrand match as closely as possible the exact form given in
Section 3. The interval of integration is
$\left[a,\infty \right)$ if
$b>0$, and
$\left(\infty ,a\right]$ if
$b<0$.
Gauss–Hermite:
$a$ and
$b$ must be chosen so as to make the integrand match as closely as possible the exact form given in
Section 3.
Constraints:
 Rational Gauss: ${\mathbf{a}}+{\mathbf{b}}\ne 0.0$;
 Gauss–Laguerre: ${\mathbf{b}}\ne 0.0$;
 Gauss–Hermite: ${\mathbf{b}}>0.0$.
 5:
$\mathbf{npts}$ – IntegerInput

On entry: the number of abscissae to be used, $n$.
Constraint:
${\mathbf{npts}}=1$, $2$, $3$, $4$, $5$, $6$, $8$, $10$, $12$, $14$, $16$, $20$, $24$, $32$, $48$ or $64$.
 6:
$\mathbf{comm}$ – Nag_User *

Pointer to a structure of type Nag_User with the following member:
 p – Pointer

On entry/exit: the pointer
$\mathbf{comm}\mathbf{\to}\mathbf{p}$, of type Pointer, allows you to communicate information to and from
f(). An object of the required type should be declared, e.g., a structure, and its address assigned to the pointer
$\mathbf{comm}\mathbf{\to}\mathbf{p}$ by means of a cast to Pointer in the calling program, e.g.,
comm.p = (Pointer)&s. The type pointer will be
void * with a C compiler that defines
void * and
char * otherwise.
 7:
$\mathbf{fail}$ – NagError *Input/Output

The NAG error argument (see
Section 3.7 in How to Use the NAG Library and its Documentation).
6
Error Indicators and Warnings
 NE_BAD_PARAM

On entry, argument
quadrule had an illegal value.
 NE_QUAD_GAUSS_CONS

Gauss–Hermite input is invalid with ${\mathbf{b}}\le 0.0$.
Constraint: ${\mathbf{b}}>0.0$.
Gauss–Laguerre input is invalid with ${\mathbf{b}}=0.0$.
Constraint: ${\mathbf{b}}\ne 0.0$.
Rational Gauss input is invalid with ${\mathbf{a}}+{\mathbf{b}}=0.0$.
Constraint: ${\mathbf{a}}+{\mathbf{b}}\ne 0.0$.
The answer is returned as zero.
 NE_QUAD_GAUSS_NPTS_RULE

The
$N$point rule is not among those stored.
The answer is evaluated for
$\u2329\mathit{\text{value}}\u232a$, the largest possible value of
npts less than the requested value,
$\u2329\mathit{\text{value}}\u232a$.
7
Accuracy
The accuracy depends on the behaviour of the integrand, and on the number of abscissae used. No tests are carried out in nag_1d_withdraw_quad_gauss_1 (d01tac) to estimate the accuracy of the result. If such an estimate is required, the function may be called more than once, with a different number of abscissae each time, and the answers compared. It is to be expected that for sufficiently smooth functions a larger number of abscissae will give improved accuracy.
Alternatively, the interval of integration may be subdivided, the integral estimated separately for each subinterval, and the sum of these estimates compared with the estimate over the whole interval.
The coding of the function
f may also have a bearing on the accuracy. For example, if a highorder Gauss–Laguerre formula is used, and the integrand is of the form
it is possible that the exponential term may underflow for some large abscissae. Depending on the machine, this may produce an error, or simply be assumed to be zero. In any case, it would be better to evaluate the expression as:
Another situation requiring care is exemplified by
The integrand here assumes very large values; for example, for
$m=63$, the peak value exceeds
$3\times {10}^{33}$. Now, if the machine holds floatingpoint numbers to an accuracy of
$k$ significant decimal digits, we could not expect such terms to cancel in the summation leaving an answer of much less than
${10}^{33k}$ (the weights being of order unity); that is instead of zero, we obtain a rather large answer through rounding error. Fortunately, such situations are characterised by great variability in the answers returned by formulae with different values of
$n$. In general, you should be aware of the order of magnitude of the integrand, and should judge the answer in that light.
8
Parallelism and Performance
nag_1d_withdraw_quad_gauss_1 (d01tac) is not threaded in any implementation.
The time taken by nag_1d_withdraw_quad_gauss_1 (d01tac) depends on the complexity of the expression for the integrand and on the number of abscissae required.
10
Example
This example evaluates the integrals
by Gauss–Legendre quadrature;
by rational Gauss quadrature with
$b=0$;
by Gauss–Laguerre quadrature with
$b=1$; and
by Gauss–Hermite quadrature with
$a=1$ and
$b=3$.
The formulae with $n=4,8,16$ are used in each case.
10.1
Program Text
Program Text (d01tace.c)
10.2
Program Data
None.
10.3
Program Results
Program Results (d01tace.r)