OPQ: A MATLAB SUITE OF PROGRAMS FOR GENERATING ORTHOGONAL POLYNOMIALS<BR>AND RELATED QUADRATURE RULES

OPQ: A MATLAB SUITE OF PROGRAMS FOR GENERATING ORTHOGONAL POLYNOMIALS
AND RELATED QUADRATURE RULES

Walter Gautschi

This set of Matlab codes is a companion piece to the book ``Orthogonal Polynomials: Computation and Approximation'', Clarendon Press, Oxford, 2004. The routines, among others, implement all computational procedures discussed therein and provide code for the examples, tables, and figures. The book is referenced below as ``OPCA''.

ORTHOGONAL POLYNOMIALS

1. Three-Term Recurrence Relation

This section contains M-files for generating the coefficients alpha_k and beta_k in the three-term recurrence relation

p_{k+1}(t)=(t-alpha_k) p_k(t) - beta_k p_{k-1}(t)

for monic orthogonal polynomials. Many of these are classical, others are not. The next two sections provide additional routines for computing these coefficients, based on the modified Chebyshev algorithm, Lanczos's algorithm, and a multiple-component discretization method. Many of the routines are transcriptions of Fortran programs in W. Gautschi, ``Algorithm 726: ORTHPOL -- a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules'', ACM Trans. Math. Software 20 (1994), 21-62. See OPCA, §§1.3, 1.5, 2.1.9, 2.2.5.

quadgj.m quadjp.m r_Bspline_dis.m r_Emfin.m r_Enumod.m r_OPhrbimod.m r_algexp.m r_algexpabslog.m r_alglog1.m r_alglog2.m r_assoc_legendre.m r_blockhahn.m r_cBspline_dis.m r_cBspline_dis_old.m r_chebleg.m r_dblexp.m r_dblexp_freud.m r_ellcheb.m r_elliptic.m r_explog.m r_explog01.m r_gjacobi.m r_hahn.m r_hatfcn_dis.m r_hatfcn_mcheb.m r_hermite.m r_hrhermite.m r_jaclog.m r_jacobi.m r_jacobi01.m r_jacplusdis.m r_jacplusdis_cheb.m r_lagplusdis.m r_laguerre.m r_logist.m r_logistic.m r_lower_subrange_jacobi.m r_lsrbinet.m r_marchenko_pastur_cheb.m r_marchenko_pastur_cheb_old.m r_marchenko_pastur_dis.m r_marchenko_pastur_dis_old.m r_marchenko_pastur_mod.m r_maxwell.m r_meixner_pollaczek.m r_modbess.m r_morse.m r_mpplusdis.m r_radtrans_dis.m r_symm_subrange_jacobi.m

2. Modified Chebyshev Algorithm. See OPCA, §2.1.7.

chebyshev.m

3. Discretization Methods. See OPCA, §§2.2.2--2.2.6.

fejer.m quadgp.m stieltjes.m lanczos.m mcdis.m mcchebyshev.m

4. Computing Cauchy Integrals

The M-files in this section compute Cauchy integrals and related functions, the kernels of the remainder term in Gauss quadrature for analytic functions, and provide routines for the Stieltjes-Perron inversion formula and Hilbert transforms. There are also two routines that apply Cauchy integrals to the problem of modifying a measure (and its recurrence coefficients) by dividing it by a polynomial. All M-files, except for the last three, are transcriptions of Fortran routines in W. Gautschi, ``Algorithm 726: ORTHPOL -- a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules'', ACM Trans. Math. Software 20(1994), 21-62. See OPCA, §2.3

cauchy.m nu0jac.m nu0lag.m nu0her.m kernel.m gchri1.m gchri2.m SPHT.m SPjac.m HTjac.m

5. Modification Algorithms

Given a weight function (or measure) w having orthogonal polynomials whose three-term recurrence relation is known, the problem here is to obtain the three-term recurrence relation of the polynomials orthogonal relative to the weight function rw, where r is a rational function. The M-files in this section provide the solution to this problem when r is a linear or quadratic factor or divisor. All except chri7.m are new. The routine chri7.m is a transcription of the routine chri, iopt=7, of W. Gautschi, ``Algorithm 726: ORTHPOL -- a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules'', ACM Trans. Math. Software 20(1994), 21-62. See OPCA, §2.4.

chri1.m chri2.m chri3.m chri4.m chri5.m chri6.m chri7.m chri8.m indop.m

6. Sobolev Orthogonal Polynomials

These are polynomials pi_k, k=0,1,2,..., of exact degrees k orthogonal with respect to an inner product involving derivatives,

(p,q)_S=int [p(t)q(t)d_0(t)+p'(t)q'(t)d_1(t)+... +p^{(s)}(t)q^{(s)}(t)d_s(t)].

Here, s>=1 is a given integer, and d_j(t), j=0,1,...s, are given positive measures, some of which may be discrete, others absolutely continuous with possibly a discrete component superimposed. If the polynomials pi_k are assumed monic, there is an infinite upper Hessenberg matrix B[1:Inf,1:Inf] with B(i+1,i)=1, i=1,2,..., such that

pi_k(t)=B(k+1,k)*t*pi_{k-1}(t)-B(1,k)*pi_{k-1}(t)- B(2,k)*pi_{k-2}(t)-...-B(k,k)*pi_0(t), k=1,2,3,... .

The coefficients involved (from left to right) in this recurrence relation are those in the k-th column of B, moving (upward) from the element 1 in position (k+1,k) to the element B(k,k) in position (1,k). The zeros of pi_k are the eigenvalues of the k-th order leading principal minor matrix B_k=B[1:k,1:k] of B. The first two M-files in this section provide routines for generating the upper triangular part of the matrix B_N, whereas the third M-file generates the zeros of pi_N. The last routine implements a Clenshaw-like algorithm for evaluating finite sums and their derivatives in Sobolev orthogonal polynomials. See OPCA, §2.5.

chebyshev_sob.m stieltjes_sob.m althammer.m sobzeros.m clenshaw_sob.m

QUADRATURE

1. Gauss-Type Quadrature

This section provides M-files for generating Gauss, Gauss-Radau, Gauss-Lobatto, generalized Gauss-Radau, and generalized Gauss-Lobatto quadrature formulae from the recurrence coefficients of the underlying weight function (or measure). See OPCA, §3.1.1.

gauss.m radau.m radau_jacobi.m radau_laguerre.m lobatto.m lobatto_jacobi.m gradau.m globatto.m

2. Gauss-Kronrod Quadrature

The two M-files in this section generate respectively the Jacobi-Kronrod matrix for a given weight function and the Gauss-Kronrod quadrature rule (if it exists). See OPCA, §3.1.2.

r_kronrod.m kronrod.m

3. Gauss-Turán Quadrature

The M-file below generates the Gauss-Turán quadrature rule. See OPCA, §3.1.3.

turan.m

4. Polynomial/Rational Gauss-Type Quadrature

Functions having poles are integrated more profitably by quadrature rules that are exact not only for polynomials, but also for elementary rational functions having the same poles, or at least those closest to the interval of integration. TheM-files in this section are designed to help generate such quadrature rules that in some sense have maximum degree of exactness. A detailed description of how they are constructed is given in W. Gautschi, ``Algorithm 793: GQRAT -- Gauss quadrature for rational functions'', ACM Trans. Math. Software 25 (1999), 213-239. See OPCA, §3.1.4.

r_mod.m quadrat.m gauss_rational.m radau_rational.m lobatto_rational.m kronrod_rational.m

5. Fermi-Dirac and Bose-Einstein Integrals

These integrals, of interest in solid state physics, are amenable to polynomial/rational quadrature methods. The M-files below implement these methods, also in cases where the integrand has poles very close to the interval of integration. See OPCA, §3.1.4.3.

fermi_dirac.m fex30.m bose_einstein.m fex31.m bose_einstein_diffpole.m fex32.m

6. Cauchy Principal Value Integrals

Two versions of Gauss quadrature are used to compute Cauchy principal value integrals. See OPCA, §3.1.5.

cauchyPVI.m

7. Polynomials Orthogonal on Several Intervals

There are Stieltjes and modified Chebyshev algorithms designed for measures on multiple domains. Both use matrix formulations of Gauss quadrature sums. See OPCA, §3.1.6.

r_multidomain_sti.m r_multidomain_cheb.m

8. Least Squares Approximation

The first routine produces discrete polynomial least squares approximants and associated Fourier coefficients. The next two routines generate special Fourier coefficients based on Gauss and Gauss-Lobatto abscissae and weights. The last routine produces discrete Sobolev least-squares approximants and related Sobolev-Fourier coefficients. See OPCA, §§3.2.1, 3.2.3.

least_squares.m fourier_gauss.m fourier_lobatto.m least_squares_sob.m

9. Weighted Newton-Cotes Quadrature

This routine generates the weights of a weighted n-point Newton-Cotes quadrature formula with given nodes and weight function w.

NewtCotes.m

UTILITY ROUTINES

A number of auxiliary Matlab routines and files are colleced here, which are useful in various contexts.

1. Mollifier Function. See OPCA, §2.1.2.

moll.m

2. Inverse of a Confluent Vandermonde Matrix. See OPCA, §2.1.4.

Tinv.m

3. Modified Moments

The first M-file below generates Legendre moments for a logarithmic weight function; it is used in r_jaclog.m. The second M-file computes Legendre moments for a squared logarithmic weight function; it is used in r_jaclogsq.m. The third M-file computes Chebyshev moments for an elliptic weight function; it is used in r_elliptic.m. See OPCA, §2.1.9.

mm_log.m mm_logsq.m mm_ell.m

4. Condition Numbers of Moment Maps. See OPCA, §§2.1.4--2.1.5.

acondGlow.m gcondG.m acondG.m rcondG.m g_n.m

5. Recurrence Coefficients. See OPCA, Examples 2.30, 2.41, 2.42.

abhrhermite   half-range Hermite

abeinstein1   Einstein

abeinstein2   Einstein squared

abfermi1   Fermi

abfermi2   Fermi squared

absqp1einstein1   Einstein modified by square root factor

absqp1einstein2   Einstein squared modified by square root factor

absqm1einstein1   Einstein modified by square root divisor

absqm1einstein2   Einstein squared modified by square root divisor

absqp1fermi1   Fermi modified by square root factor

absqp1fermi2   Fermi squared modified by square root factor

absqm1fermi1   Fermi modified by square root divisor

absqm1fermi2   Fermi squared modified by square root divisor

abjaclog   Modified logarithm

absqm1log1   Logarithm modified by square root divisor

absqm1log2   Logarithm squared modified by square root divisor