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.
C — Array of length N containing the recurrence coefficients. (Input)
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.
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
