nag_quad_1d_gauss_vec (d01uac) (PDF version)
d01 Chapter Contents
d01 Chapter Introduction
NAG Library Manual

NAG Library Function Document

nag_quad_1d_gauss_vec (d01uac)

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

nag_quad_1d_gauss_vec (d01uac) 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 semi-infinite interval (Gauss–Laguerre, rational Gauss), and an infinite interval (Gauss–Hermite).

2  Specification

#include <nag.h>
#include <nagd01.h>
void  nag_quad_1d_gauss_vec (Nag_QuadType quad_type, double a, double b, Integer n,
void (*f)(const double x[], Integer nx, double fv[], Integer *iflag, Nag_Comm *comm),
double *dinest, Nag_Comm *comm, NagError *fail)

3  Description

3.1  General

nag_quad_1d_gauss_vec (d01uac) evaluates an estimate of the definite integral of a function fx, over a finite or infinite range, by n-point Gaussian quadrature (see Davis and Rabinowitz (1975)Fröberg (1970)Ralston (1965) or Stroud and Secrest (1966)). The integral is approximated by a summation of the product of a set of weights and a set of function evaluations at a corresponding set of abscissae xi. For adjusted weights, the function values correspond to the values of the integrand f, and hence the sum will be
i=1nwifxi
where the wi are called the weights, and the xi the abscissae. A selection of values of n is available. (See Section 5.)
Where applicable, normal weights may instead be used, in which case the corresponding weight function ω is factored out of the integrand as fx=ωxgx and hence the sum will be
i=-1 n w-i gx ,
where the normal weights w-i=wiωxi are computed internally.
nag_quad_1d_gauss_vec (d01uac) uses a vectorized f to evaluate the integrand or normalized integrand at a set of abscissae, xi, for i=1,2,,nx. If adjusted weights are used, the integrand fxi must be evaluated otherwise the normalized integrand gxi must be evaluated.

3.2  Both Limits Finite

abfxdx.
The Gauss–Legendre weights and abscissae are used, and the formula is exact for any function of the form:
fx=i= 0 2n- 1 ci xi.
The formula is appropriate for functions which can be well approximated by such a polynomial over a,b. It is inappropriate for functions with algebraic singularities at one or both ends of the interval, such as 1+x-1/2 on -1,1.

3.3  One Limit Infinite

afxdx  or  -afxdx.
Two quadrature formulae are available for these integrals.
(a) The Gauss–Laguerre formula is exact for any function of the form:
fx=e-bx i= 0 2n- 1 ci xi.
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.
If the adjusted weights are selected, the complete integrand fx should be provided through f.
If the normal form is selected, the contribution of e-bx is accounted for internally, and f should only return gx, where fx=e-bxgx.
If b<0 is supplied, the interval of integration will be a,. Otherwise if b>0 is supplied, the interval of integration will be -,a.
(b) The rational Gauss formula is exact for any function of the form:
fx=i=2 2n+1cix+bi=i=0 2n-1c2n+1-ix+bix+b2n+1.
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.
Only the adjusted form of the rational Gauss formula is available, and as such, the complete integrand fx must be supplied in f.
If a+b<0, the interval of integration will be a,. Otherwise if a+b>0, the interval of integration will be -,a.

3.4  Both Limits Infinite

- +fxdx.
The Gauss–Hermite weights and abscissae are used, and the formula is exact for any function of the form:
fx = e - b x-a 2 i=0 2n- 1 ci xi ,
where b>0. Again, for general functions not of this exact form, the argument b should be chosen to match if possible the decay rate at ±.
If the adjusted weights are selected, the complete integrand fx should be provided through f.
If the normal form is selected, the contribution of e-bx-a2 is accounted for internally, and f should only return gx, where fx=e-bx-a2gx.

4  References

Davis P J and Rabinowitz P (1975) Methods of Numerical Integration Academic Press
Fröberg C E (1970) Introduction to Numerical Analysis Addison–Wesley
Ralston A (1965) A First Course in Numerical Analysis pp. 87–90 McGraw–Hill
Stroud A H and Secrest D (1966) Gaussian Quadrature Formulas Prentice–Hall

5  Arguments

1:     quad_typeNag_QuadTypeInput
On entry: indicates the quadrature formula.
quad_type=Nag_Quad_Gauss_Legendre
Gauss–Legendre quadrature on a finite interval, using normal weights.
quad_type=Nag_Quad_Gauss_Laguerre
Gauss–Laguerre quadrature on a semi-infinite interval, using normal weights.
quad_type=Nag_Quad_Gauss_Laguerre_Adjusted
Gauss–Laguerre quadrature on a semi-infinite interval, using adjusted weights.
quad_type=Nag_Quad_Gauss_Hermite
Gauss–Hermite quadrature on an infinite interval, using normal weights.
quad_type=Nag_Quad_Gauss_Hermite_Adjusted
Gauss–Hermite quadrature on an infinite interval, using adjusted weights.
quad_type=Nag_Quad_Gauss_Rational_Adjusted
Rational Gauss quadrature on a semi-infinite interval, using adjusted weights.
Constraint: quad_type=Nag_Quad_Gauss_Legendre, Nag_Quad_Gauss_Laguerre, Nag_Quad_Gauss_Laguerre_Adjusted, Nag_Quad_Gauss_Hermite, Nag_Quad_Gauss_Hermite_Adjusted or Nag_Quad_Gauss_Rational_Adjusted.
2:     adoubleInput
3:     bdoubleInput
On entry: the quantities a and b as described in the appropriate subsection of Section 3.
Constraints:
  • Rational Gauss: a+b0.0;
  • Gauss–Laguerre: b0.0;
  • Gauss–Hermite: b>0.
4:     nIntegerInput
On entry: n, the number of abscissae to be used.
Constraint: n=1, 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 20, 24, 32, 48 or 64.
If the soft fail option is used, the answer is evaluated for the largest valid value of n less than the requested value.
5:     ffunction, supplied by the userExternal Function
f must return the value of the integrand f, or the normalized integrand g, at a specified point.
The specification of f is:
void  f (const double x[], Integer nx, double fv[], Integer *iflag, Nag_Comm *comm)
1:     x[nx]const doubleInput
On entry: the abscissae, xi, for i=1,2,,nx at which function values are required.
2:     nxIntegerInput
On entry: nx, the number of abscissae.
3:     fv[nx]doubleOutput
On exit: if adjusted weights are used, the values of the integrand f. fv[i-1]=fxi, for i=1,2,,nx.
Otherwise the values of the normalized integrand g. fv[i-1]=gxi, for i=1,2,,nx.
4:     iflagInteger *Input/Output
On entry: iflag=0.
On exit: set iflag<0 if you wish to force an immediate exit from nag_quad_1d_gauss_vec (d01uac) with fail.code= NE_USER_STOP.
5:     commNag_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_1d_gauss_vec (d01uac) you may allocate memory and initialize these pointers with various quantities for use by f when called from nag_quad_1d_gauss_vec (d01uac) (see Section 3.2.1.1 in the Essential Introduction).
Some points to bear in mind when coding f are mentioned in Section 7.
6:     dinestdouble *Output
On exit: the estimate of the definite integral.
7:     commNag_Comm *Communication Structure
The NAG communication argument (see Section 3.2.1.1 in the Essential Introduction).
8:     failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

6  Error Indicators and Warnings

The value of a and/or b is invalid for the chosen quad_type. Either:
NE_BAD_PARAM
On entry, argument value had an illegal value.
The value of a and/or b is invalid.
On entry, quad_type=value.
On entry, a=value and b=value.
Constraint: a+b>0.0.
The value of a and/or b is invalid.
On entry, quad_type=value.
On entry, a=value and b=value.
Constraint: b>0.0.
The value of a and/or b is invalid.
On entry, quad_type=value.
On entry, a=value and b=value.
Constraint: b>0.0.
NE_INT
On entry, n=value.
Constraint: n>0.
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.
NE_QUAD_GAUSS_NPTS_RULE
The n-point rule is not among those stored.
On entry: n=value.
n-point rule used: n=value.
NE_TOO_SMALL
Underflow occurred in calculation of normal weights.
Reduce n or use adjusted weights: n=value.
NE_USER_STOP
Exit requested from f with iflag=value.
NE_WEIGHT_ZERO
No nonzero weights were generated for the provided parameters.

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_quad_1d_gauss_vec (d01uac) 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 range of integration may be subdivided, the integral estimated separately for each sub-interval, and the sum of these estimates compared with the estimate over the whole range.
The coding of f may also have a bearing on the accuracy. For example, if a high-order Gauss–Laguerre formula is used, and the integrand is of the form
fx=e-bxgx
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 with
fx = sgn gx × exp -bx+ lngx
Another situation requiring care is exemplified by
- +e-x2xmdx=0,  m​ odd.
The integrand here assumes very large values; for example, when m=63, the peak value exceeds 3×1033. Now, if the machine holds floating-point 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 1033-k (the weights being of order unity); that is, instead of zero we obtain a rather large answer through rounding error. 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_quad_1d_gauss_vec (d01uac) is currently neither directly nor indirectly threaded. In particular, the user-supplied argument f is not called from within a parallel region initialized inside nag_quad_1d_gauss_vec (d01uac).
The user-supplied argument f uses a vectorized interface, allowing for the required vector of function values to be evaluated in parallel; for example by placing appropriate OpenMP directives in the code for the user-supplied argument f.

9  Further Comments

The time taken by nag_quad_1d_gauss_vec (d01uac) depends on the complexity of the expression for the integrand and on the number of abscissae required.

10  Example

This example evaluates the integrals
0141+x2 dx=π
by Gauss–Legendre quadrature;
2 1x2 lnx dx =0.378671
by rational Gauss quadrature with b=0;
2e-xxdx=0.048901
by Gauss–Laguerre quadrature with b=1; and
- +e-3x2-4x-1dx=- +e-3 x+1 2e2x+2dx=1.428167
by Gauss–Hermite quadrature with a=-1 and b=3.
The formulae with n=2,4,8,16,32 and 64 are used in each case. Both adjusted and normal weights are used for Gauss–Laguerre and Gauss–Hermite quadrature.

10.1  Program Text

Program Text (d01uace.c)

10.2  Program Data

None.

10.3  Program Results

Program Results (d01uace.r)


nag_quad_1d_gauss_vec (d01uac) (PDF version)
d01 Chapter Contents
d01 Chapter Introduction
NAG Library Manual

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