FQRUL
Computes a Fejér quadrature rule with various classical weight functions.
Required Arguments
N — Number of quadrature points. (Input)
A — Lower limit of integration. (Input)
B — Upper limit of integration. (Input)
B must be greater than A.
QX — Array of length N containing quadrature points. (Output)
QW — Array of length N containing quadrature weights. (Output)
Optional Arguments
IWEIGH — Index of the weight function. (Input)
Default: IWEIGH = 1.
IWEIGH | WT(X) |
1 | 1 |
2 | 1/(X − ALPHA) |
3 | (B − X)α (X − A)β |
4 | (B − X)α (X − A)β log(X − A) |
5 | (B − X)α (X − A)β log(B − X) |
ALPHA — Parameter used in the weight function (except if IWEIGH = 1, it is ignored). (Input)
If IWEIGH = 2, then it must satisfy A.LT.ALPHA.LT.B. If IWEIGH = 3, 4, or 5, then ALPHA must be greater than −1.
Default: ALPHA= 0.0.
BETAW — Parameter used in the weight function (ignored if IWEIGH = 1 or 2). (Input)
BETAW must be greater than −1.0.
Default: BETAW= 0.0.
FORTRAN 90 Interface
Generic: CALL FQRUL (N, A, B, QX, QW [, …])
Specific: The specific interface names are S_FQRUL and D_FQRUL.
FORTRAN 77 Interface
Single: CALL FQRUL (N, A, B, IWEIGH, ALPHA, BETAW, QX, QW)
Double: The double precision name is DFQRUL.
Description
The routine FQRUL produces the weights and points for the Fejér quadrature rule. Since this computation is based on a quarter-wave cosine transform, the computations are most efficient when N, the number of points, is a product of small primes. These quadrature formulas may be an intermediate step in a more complicated situation, see for instance Gautschi and Milovanofic (1985).
The Fejér quadrature rules are based on polynomial interpolation. First, choose classical abscissas (in our case, the Gauss points for the Chebyshev weight function (1 − x2)-1/2), then derive the quadrature rule for a different weight. In order to keep the presentation simple, we will describe the case where the interval of integration is [−1, 1] even though FQRUL allows rescaling to an arbitrary interval [a, b].
We are looking for quadrature rules of the form
where the
are the zeros of the N-th Chebyshev polynomial (of the first kind) TN (x) = cos(N arccos x). The weights in the quadrature rule Q are chosen so that, for all polynomials p of degree less than N,
for some weight function w. In FQRUL, the user has the option of choosing w from five families of functions with various algebraic and logarithmic endpoint singularities.
These Fejér rules are important because they can be computed using specialized FFT quarter-wave transform routines. This means that rules with a large number of abscissas may be computed efficiently. If we insert Tl for p in the above formula, we obtain
for l = 0, …, N − 1. This is a system of linear equations for the unknown weights wj that can be simplified by noting that
and hence,
The last expression is the cosine quarter-wave forward transform for the sequence
that is implemented in
Chapter 6, “Transforms” under the name
QCOSF. More importantly,
QCOSF has an inverse
QCOSB. It follows that if the integrals on the left in the last expression can be computed, then the Fejér rule can be derived efficiently for highly composite integers
N utilizing
QCOSB. For more information on this topic, consult Davis and Rabinowitz (1984, pages 84
−86) and Gautschi (1968, page 259).
Comments
1. Workspace may be explicitly provided, if desired, by use of F2RUL/DF2RUL. The reference is:
CALL F2RUL (N, A, B, IWEIGH, ALPHA, BETAW, QX, QW, WK)
The additional argument is:
WK — Work array of length 3 * N + 15.
2. If IWEIGH specifies the weight WT(X) and the interval (A, B), then approximately
3. The routine FQRUL uses an fft, so it is most efficient when N is the product of small primes.
Example
Here, we obtain the Fejér quadrature rules using 10, 100, and 200 points. With these rules, we get successively better approximations to the integral
USE FQRUL_INT
USE UMACH_INT
USE CONST_INT
IMPLICIT NONE
INTEGER NMAX
PARAMETER (NMAX=200)
INTEGER I, K, N, NOUT
REAL A, ANSWER, B, F, QW(NMAX), &
QX(NMAX), SIN, SUM, X, PI, ERROR
INTRINSIC SIN, ABS
!
F(X) = X*SIN(41.0*PI*X**2)
! Get output unit number
CALL UMACH (2, NOUT)
!
PI = CONST('PI')
DO 20 K=1, 3
IF (K .EQ. 1) N = 10
IF (K .EQ. 2) N = 100
IF (K .EQ. 3) N = 200
A = 0.0
B = 1.0
! Get points and weights from FQRUL
CALL FQRUL (N, A, B, QX, QW)
! Evaluate the integral from these
! points and weights
SUM = 0.0
DO 10 I=1, N
SUM = SUM + F(QX(I))*QW(I)
10 CONTINUE
ANSWER = SUM
ERROR = ABS(ANSWER - 1.0/(41.0*PI))
WRITE (NOUT,99999) N, ANSWER, ERROR
20 CONTINUE
!
99999 FORMAT (/, 1X, 'When N = ', I3, ', the quadrature result making ' &
, 'use of these points ', /, ' and weights is ', 1PE11.4, &
', with error ', 1PE9.2, '.')
END
Output
When N = 10, the quadrature result making use of these points and weights is -1.6523E-01, with error 1.73E-01.
When N = 100, the quadrature result making use of these points and weights is 7.7637E-03, with error 2.79E-08.
When N = 200, the quadrature result making use of these points and weights is 7.7636E-03, with error 1.40E-08.
