nag_quad_1d_gauss (d01ba) 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).

Note: this function is scheduled to be withdrawn, please see d01ba in Advice on Replacement Calls for Withdrawn/Superseded Routines..

Syntax

[result, ifail] = d01ba(d01xxx, a, b, n, fun)

[result, ifail] = nag_quad_withdraw_1d_gauss(d01xxx, a, b, n, fun)

Description

General

nag_quad_1d_gauss (d01ba) evaluates an estimate of the definite integral of a function

f (x)

, 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

\sum_{i = 1}^{n} w_{i} f (x_{i})

where the

w_{i}

are called the weights, and the

x_{i}

the abscissae. A selection of values of

n

is available. (See Arguments.)

Both Limits Finite

\int_{a}^{b} f (x) d x .

The Gauss–Legendre weights and abscissae are used, and the formula is exact for any function of the form:

f (x) = \sum_{i = 0}^{2 n - 1} c_{i} x^{i} .

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}

[- 1, 1]

One Limit Infinite

\int_{a}^{\infty} f (x) d x or \int_{- \infty}^{a} f (x) d x .

Two quadrature formulae are available for these integrals.

(a)

The Gauss–Laguerre formula is exact for any function of the form:

f (x) = e^{- b x} \sum_{i = 0}^{2 n - 1} c_{i} x^{i} .

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.

(b)

The rational Gauss formula is exact for any function of the form:

f (x) = \sum_{i = 2}^{2 n + 1} \frac{c_{i}}{{(x + b)}^{i}} = \frac{\sum_{i = 0}^{2 n - 1} c_{2 n + 1 - i} {(x + b)}^{i}}{{(x + b)}^{2 n + 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.

Both Limits Infinite

\int_{- \infty}^{+ \infty} f (x) d x .

The Gauss–Hermite weights and abscissae are used, and the formula is exact for any function of the form:

f (x) = e^{- b {(x - a)}^{2}} \sum_{i = 0}^{2 n - 1} c_{i} x^{i} .

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

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: – string

String specifying the quadrature formula to be used:

'd01baz', for Gauss–Legendre quadrature on a finite interval;
'd01bay', for rational Gauss quadrature on a semi-infinite interval;
'd01bax', for Gauss–Laguerre quadrature on a semi-infinite interval;
'd01baw', for Gauss–Hermite quadrature on an infinite interval.

2: $a$ – double scalar

3: $b$ – double scalar

The parameters

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 One Limit Infinite(b). The range of integration is $[a \infty)$ if $a + b > 0$ , and $(- \infty a]$ 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 One Limit Infinite(a). The range of integration is $[a \infty)$ if $b > 0$ , and $(- \infty a]$ is $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 Both Limits Infinite.

Constraints:

Rational Gauss: $a + b \neq 0.0$ ;
Gauss–Laguerre: $b \neq 0.0$ ;
Gauss–Hermite: $b > 0$ .

4: $n$ – int64int32nag_int scalar

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

64

5: $fun$ – function handle or string containing name of m-file

fun must return the value of the integrand

f

at a specified point.

[result] = fun(x)

Input Parameters

1: $x$ – double scalar: The point at which the integrand $f$ must be evaluated.

Output Parameters

1: $result$ – double scalar: The value of $f (x)$ evaluated at x.

Some points to bear in mind when coding fun are mentioned in Accuracy.

Optional Input Parameters

None.

Output Parameters

1: $result$ – double scalar: The result of the function.
2: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Note: nag_quad_1d_gauss (d01ba) 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 = 1$: The N-point rule is not among those stored. If the soft fail option is used, the answer is evaluated for the largest valid value of n less than the requested value.

$ifail = 2$: The value of a and/or b is invalid.

Rational Gauss: $a + b = 0.0$ .

Gauss–Laguerre: $b = 0.0$ .

Gauss–Hermite: $b \leq 0.0$ .

If the soft fail option is used, the answer is returned as zero.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

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 (d01ba) 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 fun 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) = e^{- b x} g (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 as:

f (x) = \exp (- b x + \ln g (x))

Another situation requiring care is exemplified by

\int_{- \infty}^{+ \infty} e^{- x^{2}} x^{m} d x = 0, m ​ odd.

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 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

10^{33 - k}

(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.

Further Comments

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

Example

This example evaluates the integrals

\int_{0}^{1} \frac{4}{1 + x^{2}} d x = π

by Gauss–Legendre quadrature;

\int_{2}^{\infty} \frac{1}{x^{2} \ln x} d x = 0.378671

by rational Gauss quadrature with

b = 0

;

\int_{2}^{\infty} \frac{e^{- x}}{x} d x = 0.048901

by Gauss–Laguerre quadrature with

b = 1

; and

\int_{- \infty}^{+ \infty} e^{- 3 x^{2} - 4 x - 1} d x = \int_{- \infty}^{+ \infty} e^{- 3 {(x + 1)}^{2}} e^{2 x + 2} d x = 1.428167

by Gauss–Hermite quadrature with

a = - 1

and

b = 3

The formulae with

n = 4, 8, 16

are used in each case.

Open in the MATLAB editor: d01ba_example

function d01ba_example


fprintf('d01ba example results\n\n');

global funid;

for funid=1:4
  switch funid
    case 1
      fprintf('\nGauss-Legendre example\n');
      a = 0.0;
      b = 1.0;
      key = 'd01baz';
    case 2
      fprintf('\nRational Gauss example\n');
      a = 2.0;
      b = 0.0;
      key = 'd01bay';
    case 3
      fprintf('\nGauss-Laguerre example\n');
      a = 2.0;
      b = 1.0;
      key = 'd01bax';
    case 4
      fprintf('\nGauss-Hermite example\n');
      a = -1.0;
      b = 3.0;
      key = 'd01baw';
  end

  for i=1:3
    n = int64(2*2^i);
    [dinest, ifail] = d01ba(key, a, b, n, @f);

    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)

  global funid;

  switch funid
    case 1
      fv = 4./(1+x.*x);
    case 2
      fv = 1./(x.*x.*log(x));
    case 3
      fv = exp(-x)./x;
    case 4
      fv = exp(-3.*x.*x-4.*x-1);
    otherwise
      fv = zeros(nx, 1);
      iflag = -1;
  end

d01ba example results


Gauss-Legendre example
 4 Points   Answer =    3.14161
 8 Points   Answer =    3.14159
16 Points   Answer =    3.14159

Rational Gauss example
 4 Points   Answer =    0.37910
 8 Points   Answer =    0.37876
16 Points   Answer =    0.37869

Gauss-Laguerre example
 4 Points   Answer =    0.04887
 8 Points   Answer =    0.04890
16 Points   Answer =    0.04890

Gauss-Hermite example
 4 Points   Answer =    1.42803
 8 Points   Answer =    1.42817
16 Points   Answer =    1.42817

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox