RECQR

Computes recurrence coefficients for monic polynomials given a quadrature rule.

Required Arguments

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

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

B — Array of length NTERM containing recurrence coefficients. (Output)

C — Array of length NTERM containing recurrence coefficients. (Output)

Optional Arguments

N — Number of quadrature points. (Input)

Default: N = size (QX,1).

Default: N = size (QX,1).

NTERM — Number of recurrence coefficients. (Input)

NTERM must be less than or equal to N.

Default: NTERM = size (B,1).

NTERM must be less than or equal to N.

Default: NTERM = size (B,1).

FORTRAN 90 Interface

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

Specific: The specific interface names are S_RECQR and D_RECQR.

FORTRAN 77 Interface

Single: CALL RECQR (N, QX, QW, NTERM, B, C)

Double: The double precision name is DRECQR.

Description

The routine RECQR produces the recurrence coefficients for the orthogonal polynomials given the points and weights for the Gauss quadrature formula. It is assumed that the orthogonal polynomials are monic; hence the three-term recursion may be written

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 an inverse routine to GQRCF. Given the recurrence coefficients, the routine GQRCF produces the corresponding Gauss quadrature formula, whereas the routine RECQR produces the recurrence coefficients given the quadrature formula.

Comments

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

CALL R2CQR (N, QX, QW, NTERM, B, C, WK)

The additional argument is:

WK — Work array of length 2 * N.

2. 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). The zero-th moment

of the weight function is returned in C(1).

Example

To illustrate the use of RECQR, we will input a simple choice of recurrence coefficients, call GQRCF for the quadrature formula, put this information into RECQR, and recover the recurrence coefficients.

USE RECQR_INT

USE UMACH_INT

USE GQRCF_INT

IMPLICIT NONE

INTEGER N

PARAMETER (N=5)

INTEGER I, J, NFIX, NOUT, NTERM

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

INTRINSIC FLOAT

! Get output unit number

CALL UMACH (2, NOUT)

NFIX = 0

! Set arrays B and C of recurrence

! coefficients

DO 10 J=1, N

B(J) = FLOAT(J)

C(J) = FLOAT(J)/2.0

10 CONTINUE

WRITE (NOUT,99995)

99995 FORMAT (1X, 'Original recurrence coefficients')

WRITE (NOUT,99996) (I,B(I),I,C(I),I=1,N)

99996 FORMAT (5(6X,'B(',I1,') = ',F8.4,7X,'C(',I1,') = ',F8.5,/))

!

! The call to GQRCF will compute the

! quadrature rule from the recurrence

! coefficients given above.

!

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

WRITE (NOUT,99997)

99997 FORMAT (/, 1X, 'Quadrature rule from the recurrence coefficients' &

)

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

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

!

! Call RECQR to recover the original

! recurrence coefficients

NTERM = N

CALL RECQR (QX, QW, B, C)

WRITE (NOUT,99999)

99999 FORMAT (/, 1X, 'Recurrence coefficients determined by RECQR')

WRITE (NOUT,99996) (I,B(I),I,C(I),I=1,N)

!

END

Output

Original recurrence coefficients

B(1) = 1.0000 C(1) = 0.50000

B(2) = 2.0000 C(2) = 1.00000

B(3) = 3.0000 C(3) = 1.50000

B(4) = 4.0000 C(4) = 2.00000

B(5) = 5.0000 C(5) = 2.50000

Quadrature rule from the recurrence coefficients

QX(1) = 0.1525 QW(1) = 0.25328

QX(2) = 1.4237 QW(2) = 0.17172

QX(3) = 2.7211 QW(3) = 0.06698

QX(4) = 4.2856 QW(4) = 0.00790

QX(5) = 6.4171 QW(5) = 0.00012

Recurrence coefficients determined by RECQR

B(1) = 1.0000 C(1) = 0.50000

B(2) = 2.0000 C(2) = 1.00000

B(3) = 3.0000 C(3) = 1.50000

B(4) = 4.0000 C(4) = 2.00000

B(5) = 5.0000 C(5) = 2.50000