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.

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.

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.

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.

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.