hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_quad_1d_gauss_vec (d01ua)

Purpose

nag_quad_1d_gauss_vec (d01ua) 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).

Syntax

[dinest, user, ifail] = d01ua(key, a, b, n, f, 'user', user)
[dinest, user, ifail] = nag_quad_1d_gauss_vec(key, a, b, n, f, 'user', user)

Description

General

nag_quad_1d_gauss_vec (d01ua) evaluates an estimate of the definite integral of a function f(x)f(x), over a finite or infinite range, by nn-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 xixi. For adjusted weights, the function values correspond to the values of the integrand ff, and hence the sum will be
n
wif(xi)
i = 1
i=1nwif(xi)
where the wiwi are called the weights, and the xixi the abscissae. A selection of values of nn is available. (See Section [Parameters].)
Where applicable, normal weights may instead be used, in which case the corresponding weight function ωω is factored out of the integrand as f(x) = ω(x)g(x)f(x)=ω(x)g(x) and hence the sum will be
n
wig(x),
i = 1
i=-1 n w-i g(x) ,
where the normal weights wi = wiω(xi)w-i=wiω(xi) are computed internally.
nag_quad_1d_gauss_vec (d01ua) uses a vectorized f to evaluate the integrand or normalized integrand at a set of abscissae, xixi, for i = 1,2,,nxi=1,2,,nx. If adjusted weights are used, the integrand f(xi)f(xi) must be evaluated otherwise the normalized integrand g(xi)g(xi) must be evaluated.

Both Limits Finite

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

One Limit Infinite

a
f(x)dx  or  f(x)dx.
a
af(x)dx  or  -af(x)dx.
Two quadrature formulae are available for these integrals.
(a) The Gauss–Laguerre formula is exact for any function of the form:
2n 1
f(x) = ebxcixi.
i = 0
f(x)=e-bx i= 0 2n- 1 ci xi.
This formula is appropriate for functions decaying exponentially at infinity; the parameter bb should be chosen if possible to match the decay rate of the function.
If the adjusted weights are selected, the complete integrand f(x)f(x) should be provided through f.
If the normal form is selected, the contribution of ebxe-bx is accounted for internally, and f should only return g(x)g(x), where f(x) = ebxg(x)f(x)=e-bxg(x).
If b < 0b<0 is supplied, the interval of integration will be [a,)[a,). Otherwise if b > 0b>0 is supplied, the interval of integration will be (,a](-,a].
(b) The rational Gauss formula is exact for any function of the form:
2n + 1
f(x) = (ci)/((x + b)i) = (i = 02n1c2n + 1i(x + b)i)/((x + b)2n + 1).
i = 2
f(x)=i=2 2n+1ci(x+b)i=i=0 2n-1c2n+1-i(x+b)i(x+b)2n+1.
This formula is likely to be more accurate for functions having only an inverse power rate of decay for large xx. Here the choice of a suitable value of bb may be more difficult; unfortunately a poor choice of bb 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 f(x)f(x) must be supplied in f.
If a + b < 0a+b<0, the interval of integration will be [a,)[a,). Otherwise if a + b > 0a+b>0, the interval of integration will be (,a](-,a].

Both Limits Infinite

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

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

Parameters

Compulsory Input Parameters

1:     key – int64int32nag_int scalar
Indicates the quadrature formula.
key = 0key=0
Gauss–Legendre quadrature on a finite interval, using normal weights.
key = 3key=3
Gauss–Laguerre quadrature on a semi-infinite interval, using normal weights.
key = -3key=-3
Gauss–Laguerre quadrature on a semi-infinite interval, using adjusted weights.
key = 4key=4
Gauss–Hermite quadrature on an infinite interval, using normal weights.
key = -4key=-4
Gauss–Hermite quadrature on an infinite interval, using adjusted weights.
key = -5key=-5
Rational Gauss quadrature on a semi-infinite interval, using adjusted weights.
Constraint: key = 0key=0, 33, -3-3, 44, -4-4 or -5-5.
2:     a – double scalar
3:     b – double scalar
The quantities aa and bb as described in the appropriate subsection of Section [Description].
Constraints:
  • Rational Gauss: a + b0.0a+b0.0;
  • Gauss–Laguerre: b0.0b0.0;
  • Gauss–Hermite: b > 0b>0.
4:     n – int64int32nag_int scalar
nn, the number of abscissae to be used.
Constraint: n = 1n=1, 22, 33, 44, 55, 66, 88, 1010, 1212, 1414, 1616, 2020, 2424, 3232, 4848 or 6464.
If the soft fail option is used, the answer is evaluated for the largest valid value of n less than the requested value.
5:     f – function handle or string containing name of m-file
f must return the value of the integrand ff, or the normalized integrand gg, at a specified point.
[fv, iflag, user] = f(x, nx, iflag, user)

Input Parameters

1:     x(nx) – double array
The abscissae, xixi, for i = 1,2,,nxi=1,2,,nx at which function values are required.
2:     nx – int64int32nag_int scalar
nxnx, the number of abscissae.
3:     iflag – int64int32nag_int scalar
iflag = 0iflag=0.
4:     user – Any MATLAB object
f is called from nag_quad_1d_gauss_vec (d01ua) with the object supplied to nag_quad_1d_gauss_vec (d01ua).

Output Parameters

1:     fv(nx) – double array
If adjusted weights are used, the values of the integrand ff. fv(i) = f(xi)fvi=f(xi), for i = 1,2,,nxi=1,2,,nx.
Otherwise the values of the normalized integrand gg. fv(i) = g(xi)fvi=g(xi), for i = 1,2,,nxi=1,2,,nx.
2:     iflag – int64int32nag_int scalar
Set iflag < 0iflag<0 if you wish to force an immediate exit from nag_quad_1d_gauss_vec (d01ua) with ifail = 1ifail=-1.
3:     user – Any MATLAB object
Some points to bear in mind when coding f are mentioned in Section [Accuracy].

Optional Input Parameters

1:     user – Any MATLAB object
user is not used by nag_quad_1d_gauss_vec (d01ua), but is passed to f. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Input Parameters Omitted from the MATLAB Interface

iuser ruser

Output Parameters

1:     dinest – double scalar
The estimate of the definite integral.
2:     user – Any MATLAB object
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_quad_1d_gauss_vec (d01ua) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W ifail = 1ifail=1
The nn-point rule is not among those stored.
W ifail = 2ifail=2
Underflow occurred in calculation of normal weights.
W ifail = 3ifail=3
No nonzero weights were generated for the provided parameters.
  ifail = 11ifail=11
On entry, key = _key=_.
Constraint: key = 0key=0, 33, -3-3, 44, -4-4 or -5-5.
  ifail = 12ifail=12
The value of a and/or b is invalid for the chosen key. Either:
  • Constraint: |a + b| > 0.0|a+b|>0.0.
  • Constraint: |b| > 0.0|b|>0.0.
  • Constraint: b > 0.0b>0.0.
  ifail = 14ifail=14
Constraint: n > 0n>0.
W ifail = 1ifail=-1
Exit requested from f with .

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 (d01ua) 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
f(x) = ebxg(x)
f(x)=e-bxg(x)
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
f(x) = sgn (g(x)) × exp(bx + ln|g(x)|)
f(x) = sgn (g(x)) × exp( -bx+ ln|g(x)| )
Another situation requiring care is exemplified by
+
ex2xmdx = 0,  m​ odd.
- +e-x2xmdx=0,  m​ odd.
The integrand here assumes very large values; for example, when m = 63m=63, the peak value exceeds 3 × 10333×1033. Now, if the machine holds floating point numbers to an accuracy of kk significant decimal digits, we could not expect such terms to cancel in the summation leaving an answer of much less than 1033k1033-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 nn.
In general, you should be aware of the order of magnitude of the integrand, and should judge the answer in that light.

Further Comments

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

Example

function nag_quad_1d_gauss_vec_example
  for funid=1:6
    switch funid
      case 1
        fprintf('\nGauss-Legendre example\n');
        a = 0.0;
        b = 1.0;
        key = int64(0);
      case 2
        fprintf('\nRational Gauss example\n');
        a = 2.0;
        b = 0.0;
        key = int64(-5);
      case 3
        fprintf('\nGauss-Laguerre example (adjusted weights)\n');
        a = 2.0;
        b = 1.0;
        key = int64(-3);
      case 4
        fprintf('\nGauss-Laguerre example (normal weights)\n');
        a = 2.0;
        b = 1.0;
        key = int64(3);
      case 5
        fprintf('\nGauss-Hermite example (adjusted weights)\n');
        a = -1.0;
        b = 3.0;
        key = int64(-4);
      case 6
        fprintf('\nGauss-Hermite example (normal weights)\n');
        a = -1.0;
        b = 3.0;
        key = int64(4);
    end

    for i=1:6
      n = int64(2^i);
      [dinest, user, ifail] = nag_quad_1d_gauss_vec(key, a, b, n, @f, 'user', funid);

      if ifail == 0 || ifail == 1
        fprintf('%2d Points   Answer = %10.5f\n', n, dinest);
      end
    end
  end


function [fv, iflag, user] = f(x, nx, iflag, user)
  switch user
    case 1
      fv = 4./(1+x.*x);
    case 2
      fv = 1./(x.*x.*log(x));
    case 3
      fv = exp(-x)./x;
    case 4
      fv = 1./x;
    case 5
      fv = exp(-3.*x.*x-4.*x-1);
    case 6
      fv = exp(2*x+2);
    otherwise
      fv = zeros(nx, 1);
      iflag = -1;
  end
 

Gauss-Legendre example
 2 Points   Answer =    3.14754
 4 Points   Answer =    3.14161
 8 Points   Answer =    3.14159
16 Points   Answer =    3.14159
32 Points   Answer =    3.14159
64 Points   Answer =    3.14159

Rational Gauss example
 2 Points   Answer =    0.37989
 4 Points   Answer =    0.37910
 8 Points   Answer =    0.37876
16 Points   Answer =    0.37869
32 Points   Answer =    0.37867
64 Points   Answer =    0.37867

Gauss-Laguerre example (adjusted weights)
 2 Points   Answer =    0.04833
 4 Points   Answer =    0.04887
 8 Points   Answer =    0.04890
16 Points   Answer =    0.04890
32 Points   Answer =    0.04890
64 Points   Answer =    0.04890

Gauss-Laguerre example (normal weights)
 2 Points   Answer =    0.04833
 4 Points   Answer =    0.04887
 8 Points   Answer =    0.04890
16 Points   Answer =    0.04890
32 Points   Answer =    0.04890
64 Points   Answer =    0.04890

Gauss-Hermite example (adjusted weights)
 2 Points   Answer =    1.38381
 4 Points   Answer =    1.42803
 8 Points   Answer =    1.42817
16 Points   Answer =    1.42817
32 Points   Answer =    1.42817
64 Points   Answer =    1.42817

Gauss-Hermite example (normal weights)
 2 Points   Answer =    1.38381
 4 Points   Answer =    1.42803
 8 Points   Answer =    1.42817
16 Points   Answer =    1.42817
32 Points   Answer =    1.42817
64 Points   Answer =    1.42817

function d01ua_example
  for funid=1:6
    switch funid
      case 1
        fprintf('\nGauss-Legendre example\n');
        a = 0.0;
        b = 1.0;
        key = int64(0);
      case 2
        fprintf('\nRational Gauss example\n');
        a = 2.0;
        b = 0.0;
        key = int64(-5);
      case 3
        fprintf('\nGauss-Laguerre example (adjusted weights)\n');
        a = 2.0;
        b = 1.0;
        key = int64(-3);
      case 4
        fprintf('\nGauss-Laguerre example (normal weights)\n');
        a = 2.0;
        b = 1.0;
        key = int64(3);
      case 5
        fprintf('\nGauss-Hermite example (adjusted weights)\n');
        a = -1.0;
        b = 3.0;
        key = int64(-4);
      case 6
        fprintf('\nGauss-Hermite example (normal weights)\n');
        a = -1.0;
        b = 3.0;
        key = int64(4);
    end

    for i=1:6
      n = int64(2^i);
      [dinest, user, ifail] = d01ua(key, a, b, n, @f, 'user', funid);

      if ifail == 0 || ifail == 1
        fprintf('%2d Points   Answer = %10.5f\n', n, dinest);
      end
    end
  end


function [fv, iflag, user] = f(x, nx, iflag, user)
  switch user
    case 1
      fv = 4./(1+x.*x);
    case 2
      fv = 1./(x.*x.*log(x));
    case 3
      fv = exp(-x)./x;
    case 4
      fv = 1./x;
    case 5
      fv = exp(-3.*x.*x-4.*x-1);
    case 6
      fv = exp(2*x+2);
    otherwise
      fv = zeros(nx, 1);
      iflag = -1;
  end
 

Gauss-Legendre example
 2 Points   Answer =    3.14754
 4 Points   Answer =    3.14161
 8 Points   Answer =    3.14159
16 Points   Answer =    3.14159
32 Points   Answer =    3.14159
64 Points   Answer =    3.14159

Rational Gauss example
 2 Points   Answer =    0.37989
 4 Points   Answer =    0.37910
 8 Points   Answer =    0.37876
16 Points   Answer =    0.37869
32 Points   Answer =    0.37867
64 Points   Answer =    0.37867

Gauss-Laguerre example (adjusted weights)
 2 Points   Answer =    0.04833
 4 Points   Answer =    0.04887
 8 Points   Answer =    0.04890
16 Points   Answer =    0.04890
32 Points   Answer =    0.04890
64 Points   Answer =    0.04890

Gauss-Laguerre example (normal weights)
 2 Points   Answer =    0.04833
 4 Points   Answer =    0.04887
 8 Points   Answer =    0.04890
16 Points   Answer =    0.04890
32 Points   Answer =    0.04890
64 Points   Answer =    0.04890

Gauss-Hermite example (adjusted weights)
 2 Points   Answer =    1.38381
 4 Points   Answer =    1.42803
 8 Points   Answer =    1.42817
16 Points   Answer =    1.42817
32 Points   Answer =    1.42817
64 Points   Answer =    1.42817

Gauss-Hermite example (normal weights)
 2 Points   Answer =    1.38381
 4 Points   Answer =    1.42803
 8 Points   Answer =    1.42817
16 Points   Answer =    1.42817
32 Points   Answer =    1.42817
64 Points   Answer =    1.42817


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013