GQRCF

Computes a Gauss, Gauss-Radau or Gauss-Lobatto quadrature rule given the recurrence coefficients for the monic polynomials orthogonal with respect to the weight function.

Required Arguments

N — Number of quadrature points. (Input)

B — Array of length N containing the recurrence coefficients. (Input)

See Comments for definitions.

See Comments for definitions.

C — Array of length N containing the recurrence coefficients. (Input)

See Comments for definitions.

See Comments for definitions.

QX — Array of length N containing quadrature points. (Output)

QW — Array of length N containing quadrature weights. (Output)

Optional Arguments

NFIX — Number of fixed quadrature points. (Input)

NFIX = 0, 1 or 2. For the usual Gauss quadrature rules NFIX = 0.

Default: NFIX = 0.

NFIX = 0, 1 or 2. For the usual Gauss quadrature rules NFIX = 0.

Default: NFIX = 0.

QXFIX — Array of length NFIX (ignored if NFIX = 0) containing the preset quadrature point(s). (Input)

FORTRAN 90 Interface

Generic: CALL GQRCF (N, B, C, QX, QW [, …])

Specific: The specific interface names are S_GQRCF and D_GQRCF.

FORTRAN 77 Interface

Single: CALL GQRCF (N, B, C, NFIX, QXFIX, QX, QW)

Double: The double precision name is DGQRCF.

Description

The routine GQRCF produces the points and weights for the Gauss, Gauss-Radau, or Gauss-Lobatto quadrature formulas given the three-term recurrence relation for the orthogonal polynomials. In particular, it is assumed that the orthogonal polynomials are monic, and hence, the three-term recursion may be written as

where p0 = 1 and p-1 = 0. It is obvious from this representation that the degree of pi is i and that pi is monic. In order for the recurrence to give rise to a sequence of orthogonal polynomials (with respect to a nonnegative measure), it is necessary and sufficient that ci > 0. This routine is a modification of the subroutine GAUSSQUADRULE (Golub and Welsch 1969). In the simple case when NFIX = 0, the routine returns points in x = QX and weights in w = QW so that

for all functions f that are polynomials of degree less than 2N. Here, w is any weight function for which the above recurrence produces the orthogonal polynomials pi on the interval [a, b] and w is normalized by

If NFIX = 1, then one of the above xi equals the first component of QXFIX. Similarly, if NFIX = 2, then two of the components of x will equal the first two components of QXFIX. In general, the accuracy of the above quadrature formula degrades when NFIX increases. The quadrature rule will integrate all functions f that are polynomials of degree less than 2N − NFIX.

Comments

1. Workspace may be explicitly provided, if desired, by use of G2RCF/DG2RCF. The reference is:

CALL G2RCF (N, B, C, NFIX, QXFIX, QX, QW, WK)

The additional argument is:

WK — Work array of length N.

2. Informational error

Type | Code | Description |
---|---|---|

4 | 1 | No convergence in 100 iterations. |

3. The recurrence coefficients B(I) and C(I) define the monic polynomials via the relation P(I) = (X − B(I + 1)) * P(I − 1) − C(I + 1) * P(I − 2). C(1) contains the zero-th moment

of the weight function. Each element of C must be greater than zero.

4. If WT(X) is the weight specified by the coefficients and the interval is (a, b), then approximately

5. Gaussian quadrature is always the method of choice when the function F(X) behaves like a polynomial. Gaussian quadrature is also useful on infinite intervals (with appropriate weight functions) because other techniques often fail.

Example

We compute the Gauss quadrature rule (with N = 6) for the Chebyshev weight, (1 + x2) (-1/2), from the recurrence coefficients. These coefficients are obtained by a call to the IMSL routine RECCF.

USE GQRCF_INT

USE UMACH_INT

USE RECCF_INT

IMPLICIT NONE

INTEGER N

PARAMETER (N=6)

INTEGER I, NFIX, NOUT

REAL B(N), C(N), QW(N), QX(N), QXFIX(2)

! Get output unit number

CALL UMACH (2, NOUT)

! Recursion coefficients will come from

! routine RECCF.

! The call to RECCF finds recurrence

! coefficients for Chebyshev

! polynomials of the 1st kind.

CALL RECCF (N, B, C)

!

! The call to GQRCF will compute the

! quadrature rule from the recurrence

! coefficients determined above.

CALL GQRCF (N, B, C, QX, QW)

WRITE (NOUT,99999) (I,QX(I),I,QW(I),I=1,N)

99999 FORMAT (6(6X,'QX(',I1,') = ',F8.4,7X,'QW(',I1,') = ',F8.5,/))

!

END

Output

QX(1) = -0.9325 QW(1) = 0.17132

QX(2) = -0.6612 QW(2) = 0.36076

QX(3) = -0.2386 QW(3) = 0.46791

QX(4) = 0.2386 QW(4) = 0.46791

QX(5) = 0.6612 QW(5) = 0.36076

QX(6) = 0.9325 QW(6) = 0.17132