18 Orthogonal Polynomials Computation 18.39 Applications in the Physical Sciences 18.41 Tables

§18.40 Methods of Computation

ⓘ

Referenced by:: Erratum (V1.2.0) Section 18.40, Version 1.2.0 (March 27, 2024)
Permalink:: http://dlmf.nist.gov/18.40
See also:: Annotations for Ch.18

§18.40(i) Computation of Polynomials
§18.40(ii) The Classical Moment Problem

§18.40(i) Computation of Polynomials

ⓘ

Keywords:: Clenshaw’s algorithm, classical orthogonal polynomials, differential equations, in series of Chebyshev polynomials, solutions
Referenced by:: Erratum (V1.2.0) Section 18.40, Version 1.2.0 (March 27, 2024)
Permalink:: http://dlmf.nist.gov/18.40.i
Addition (effective with 1.2.0):: This subsection was added.
See also:: Annotations for §18.40 and Ch.18

Orthogonal polynomials can be computed from their explicit polynomial form by Horner’s scheme (§1.11(i)). Usually, however, other methods are more efficient, especially the numerical solution of difference equations (§3.6) and the application of uniform asymptotic expansions (when available) for OP’s of large degree. For applications in which the OP’s appear only as terms in series expansions (compare §18.18(i)) the need to compute them can be avoided altogether by use instead of Clenshaw’s algorithm (§3.11(ii)) and its straightforward generalization to OP’s other than Chebyshev. For further information see Clenshaw (1955), Gautschi (2004, §§2.1, 8.1), and Mason and Handscomb (2003, §2.4).

§18.40(ii) The Classical Moment Problem

ⓘ

Keywords:: Stieltjes–Perron inversion
Referenced by:: §18.2(ix), §18.2(x), Figure 18.39.2, Figure 18.39.2, §18.39(iv), Erratum (V1.2.0) Section 18.40, Version 1.2.0 (March 27, 2024)
Permalink:: http://dlmf.nist.gov/18.40.ii
Addition (effective with 1.2.0):: This subsection was added.
See also:: Annotations for §18.40 and Ch.18

The problem of moments is simply stated and the early work of Stieltjes, Markov, and Chebyshev on this problem was the origin of the understanding of the importance of both continued fractions and OP’s in many areas of analysis. Given the power moments, $\mu_{n}=\int_{a}^{b}x^{n}\,\mathrm{d}\mu(x)$ , $n=0,1,2,\dots$ , can these be used to find a unique $\mu(x)$ , a non-decreasing, real, function of , in the case that the moment problem is determined? Should a unique solution not exist the moment problem is then indeterminant. The theory behind these remarks is in Shohat and Tamarkin (1970), Akhiezer (2021), Chihara (1978).

A numerical approach to the recursion coefficients and quadrature abscissas and weights

ⓘ

See also:: Annotations for §18.40(ii), §18.40 and Ch.18

In what follows we consider only the simple, illustrative, case that $\mu(x)$ is continuously differentiable so that $\,\mathrm{d}\mu(x)=w(x)\,\mathrm{d}x$ , with real, positive, and continuous on a real interval The strategy will be to: 1) use the moments to determine the recursion coefficients $\alpha_{n},\beta_{n}$ of equations (18.2.11_5) and (18.2.11_8); then, 2) to construct the quadrature abscissas $x_{i}$ and weights (or Christoffel numbers) $w_{i}$ from the J-matrix of §3.5(vi), equations (3.5.31) and(3.5.32). These quadrature weights and abscissas will then allow construction of a convergent sequence of approximations to , as will be considered in the following paragraphs.

There are many ways to implement these first two steps, noting that the expressions for $\alpha_{n}$ and $\beta_{n}$ of equation (18.2.30) are of little practical numerical value, see Gautschi (2004) and Golub and Meurant (2010). A simple set of choices is spelled out in Gordon (1968) which gives a numerically stable algorithm for direct computation of the recursion coefficients in terms of the moments, followed by construction of the J-matrix and quadrature weights and abscissas, and we will follow this approach: Let be a positive integer and define

18.40.1

$P_{1,1}=1,~P_{m,1}=0,$ $m=2,\dots,2N+4$ ,

$P_{m,2}=(-1)^{m-1}\mu_{m-1},$ $m=1,\dots,2N+3$ ,

$P_{m,n}=P_{1,n-1}P_{m+1,n-2}-P_{1,n-2}P_{m+1,n-1},$ $n\geq 3$ , $m\geq 1$ , $m+n\leq 2N+5$ ,

ⓘ

Symbols:: : positive integer, : nonnegative integer and : nonnegative integer
Referenced by:: Erratum (V1.2.0) Section 18.40
Permalink:: http://dlmf.nist.gov/18.40.E1
Encodings:: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

use the first row of this -matrix for

18.40.2 $a_{1}=\mu_{0},~a_{n}=\frac{P_{1,n+1}}{P_{1,n}P_{1,n-1}},$ $n=2,\dots,2N+3$ ,

ⓘ

Symbols:: : positive integer and : nonnegative integer
Permalink:: http://dlmf.nist.gov/18.40.E2
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

and these can be used for the recursion coefficients

18.40.3 $\alpha_{0}=a_{2},~\alpha_{n}=a_{2n+1}+a_{2n+2},~\beta_{n}=a_{2n}a_{2n+1},$ $n=1,\dots,N$ .

ⓘ

Symbols:: : positive integer and : nonnegative integer
Permalink:: http://dlmf.nist.gov/18.40.E3
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

The quadrature abscissas $x_{n}$ and weights $w_{n}$ then follow from the discussion of §3.5(vi). See Gautschi (1983) for examples of numerically stable and unstable use of the above recursion relations, and how one can then usefully differentiate between numerical results of low and high precision, as produced thereby.

Having now directly connected computation of the quadrature abscissas and weights to the moments, what follows uses these for a Stieltjes–Perron inversion to regain .

Stieltjes Inversion via (approximate) Analytic Continuation

ⓘ

Keywords:: Stieltjes–Perron inversion, analytic continuation
See also:: Annotations for §18.40(ii), §18.40 and Ch.18

We have from (18.2.38) that

18.40.4 $\lim_{N\to\infty}F_{N}(z)=F(z)\equiv\frac{1}{\mu_{0}}\int_{a}^{b}\frac{w(x)\,% \mathrm{d}x}{z-x},$ $z\in\mathbb{C}\backslash[a,b]$ ,

ⓘ

Symbols:: $[\NVar{a},\NVar{b}]$ : closed interval, $\mathbb{C}$ : complex plane, $\,\mathrm{d}\NVar{x}$ : differential of , $\in$ : element of, $\equiv$ : equals by definition, $\int$ : integral, : positive integer, : weight function, : complex variable and : real variable
Proved:: Ismail (2009, p. 36)(proved)
Permalink:: http://dlmf.nist.gov/18.40.E4
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

in which

18.40.5 $F_{N}(z)=\frac{1}{\mu_{0}}\sum_{n=1}^{N}\frac{w_{n}}{z-x_{n}}.$

ⓘ

Symbols:: : positive integer, $w_{x}$ : weights, : complex variable, : nonnegative integer and : real variable
Permalink:: http://dlmf.nist.gov/18.40.E5
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

Let $x^{\prime}\in(a,b)$ . It is now necessary to take the limit $\varepsilon\to 0{+}$ of $F(x^{\prime}+\mathrm{i}\varepsilon)$ , and the imaginary part is the required Stieltjes–Perron inversion:

18.40.6 $\lim_{\varepsilon\to 0{+}}\int_{a}^{b}\frac{w(x)\,\mathrm{d}x}{x^{\prime}+% \mathrm{i}\varepsilon-x}\,\mathrm{d}x=\pvint_{a}^{b}\frac{w(x)\,\mathrm{d}x}{x% ^{\prime}-x}-\mathrm{i}\pi w(x^{\prime}),$

ⓘ

Symbols:: $\,\mathrm{d}\NVar{x}$ : differential of , $\mathrm{i}$ : imaginary unit, $\int$ : integral, $\pvint_{\NVar{a}}^{\NVar{b}}$ : Cauchy principal value, : weight function and : real variable
Permalink:: http://dlmf.nist.gov/18.40.E6
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

The question is then: how is this possible given only $F_{N}(z)$ , rather than itself? $F_{N}(z)$ often converges to smooth results for off the real axis for $\Im{z}$ at a distance greater than the pole spacing of the $x_{n}$ , this may then be followed by approximate numerical analytic continuation via fitting to lower order continued fractions (either Padé, see §3.11(iv), or pointwise continued fraction approximants, see Schlessinger (1968, Appendix)), to $F_{N}(z)$ and evaluating these on the real axis in regions of higher pole density that those of the approximating function. Results of low (2 to 3 decimal digits) precision for are easily obtained for $N\sim 10$ to 20. Gautschi (2004, p. 119–120) has explored the $\varepsilon\to 0^{+}$ limit via the Wynn $\varepsilon$ -algorithm, (3.9.11) to accelerate convergence, finding four to eight digits of precision in , depending smoothly on $x^{\prime}$ , for $N\approx 4000$ , for an example involving first numerator Legendre OP’s.

Histogram Approach

ⓘ

Keywords:: Stieltjes–Perron inversion, graphic, histogram construction
See also:: Annotations for §18.40(ii), §18.40 and Ch.18

The quadrature points and weights can be put to a more direct and efficient use. Define

18.40.7 $\mu_{N}(x)=\sum_{n=1}^{N}w_{n}H\left(x-x_{n}\right),$ $x\in(a,b)$ ,

ⓘ

Symbols:: $H\left(\NVar{x}\right)$ : Heaviside function, $\in$ : element of, $(\NVar{a},\NVar{b})$ : open interval, : positive integer, $w_{x}$ : weights, : nonnegative integer and : real variable
Source:: Shohat and Tamarkin (1970, pp. 42-43)
Referenced by:: §18.40(ii)
Permalink:: http://dlmf.nist.gov/18.40.E7
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

$H\left(x\right)$ being the Heaviside step-function, see (1.16.13).

Equation (18.40.7) provides step-histogram approximations to $\int_{a}^{x}\,\mathrm{d}\mu(x)$ , as shown in Figure 18.40.1 for and 120, shown here for the repulsive Coulomb–Pollaczek OP’s of Figure 18.39.2, with the parameters as listed therein.

Figure 18.40.1: Histogram approximations to the Repulsive Coulomb–Pollaczek, RCP, weight function integrated over

, see Figure 18.39.2 for an exact result, for

, shown for

and

. Magnify

ⓘ

Symbols:: $[\NVar{a},\NVar{b})$ : half-closed interval, : positive integer and : real variable
Referenced by:: §18.40(ii), Erratum (V1.2.0) Section 18.40, Version 1.2.0 (March 27, 2024)
Permalink:: http://dlmf.nist.gov/18.40.F1
Encodings:: pdf, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

The bottom and top of the steps at the $x_{i}$ are lower and upper bounds to $\int_{a}^{x_{i}}\,\mathrm{d}\mu(x)$ as made explicit via the Chebyshev inequalities discussed by Shohat and Tamarkin (1970, pp. 42–43). Interpolation of the midpoints of the jumps followed by differentiation with respect to yields a Stieltjes–Perron inversion to obtain $w^{\mathrm{RCP}}(x)$ to a precision of $\sim 4$ decimal digits for . Convergence is $\sim O\left(N^{-2}\right)$ . Results similar to these appear in Langhoff et al. (1976) in methods developed for physics applications, and which includes treatments of systems with discontinuities in $\mu(x)$ , using what is referred to as the Stieltjes derivative which may be traced back to Stieltjes, as discussed by Deltour (1968, Eq. 12).

Derivative Rule Approach

ⓘ

Keywords:: Stieltjes–Perron inversion, derivative rule, graphic, histogram construction
See also:: Annotations for §18.40(ii), §18.40 and Ch.18

An alternate, and highly efficient, approach follows from the derivative rule conjecture, see Yamani and Reinhardt (1975), and references therein, namely that

18.40.8 $w(x_{i,N})\approx\frac{w_{i,N}}{\left.\frac{\mathrm{d}x(j,N)}{\mathrm{d}j}% \right|_{j=i}}.$

ⓘ

Symbols:: $\frac{\mathrm{d}\NVar{f}}{\mathrm{d}\NVar{x}}$ : derivative of with respect to , $(\NVar{a},\NVar{b})$ : open interval, : positive integer, : weight function, $w_{x}$ : weights and : real variable
Referenced by:: Figure 18.40.2, Figure 18.40.2
Permalink:: http://dlmf.nist.gov/18.40.E8
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

This allows Stieltjes–Perron inversion for the $w(x_{i,N})$ , given the quadrature weights and points. Here is an interpolation of the abscissas $x_{i,N},i=1,2,\dots,N$ , that is, $x(i,N)=x_{i,N}$ , allowing differentiation by . In what follows this is accomplished in two ways: i) via the Lagrange interpolation of §3.3(i) ; and ii) by constructing a pointwise continued fraction, or PWCF, as follows:

18.40.9 $x(t,N)=\cfrac{x_{1,N}}{1+\cfrac{a_{1}(t-1)}{1+\cfrac{a_{2}(t-2)}{1+\cdots}}}% \frac{a_{N-1}(t-(N-1))}{1},$ $t\in(0,\infty)$ ,

ⓘ

Symbols:: $\in$ : element of, $(\NVar{a},\NVar{b})$ : open interval, : positive integer, : real variable and : real variable
Source:: Schlessinger (1968)
Permalink:: http://dlmf.nist.gov/18.40.E9
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

where the coefficients are defined recursively via $a_{1}=\frac{x_{1,N}}{x_{2,N}}-1$ , and

18.40.10 $a_{\ell}=-1-\cfrac{2a_{\ell-1}}{1+\cfrac{3a_{\ell-2}}{1+\cdots}}\frac{\ell a_{% 1}}{1-\frac{x_{1,N}}{x_{\ell+1,N}}},$ $\ell=2,\dots,N-1$ .

ⓘ

Symbols:: : positive integer, $\ell$ : nonnegative integer and : real variable
Referenced by:: Erratum (V1.2.0) Section 18.40
Permalink:: http://dlmf.nist.gov/18.40.E10
Encodings:: TeX, pMML, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

The PWCF is a minimally oscillatory algebraic interpolation of the abscissas $x_{i,N},i=1,2,\dots,N$ .

Comparisons of the precisions of Lagrange and PWCF interpolations to obtain the derivatives, are shown in Figure 18.40.2. The example chosen is inversion from the $\alpha_{n},\beta_{n}$ for the weight function for the repulsive Coulomb–Pollaczek, RCP, polynomials of (18.39.50). This is a challenging case as the desired $w^{\mathrm{RCP}}(x)$ on has an essential singularity at .

Figure 18.40.2: Derivative Rule inversions for $w^{\mathrm{RCP}}(x)$ carried out via Lagrange and PWCF interpolations. Shown are the absolute errors of approximation (18.40.8) at the points $x_{i,N}$ , $i=1,2,\dots,N$ for

. For the derivative rule Lagrange interpolation (red points) gives $\sim 15$ digits in the central region, while PWCF interpolation (blue points) gives $\sim 25$ . Magnify

ⓘ

Symbols:: : positive integer, : weight function and : real variable
Referenced by:: §18.40(ii), §18.40(ii), Erratum (V1.2.0) Section 18.40, Version 1.2.0 (March 27, 2024)
Permalink:: http://dlmf.nist.gov/18.40.F2
Encodings:: pdf, png
See also:: Annotations for §18.40(ii), §18.40(ii), §18.40 and Ch.18

Further, exponential convergence in , via the Derivative Rule, rather than the power-law convergence of the histogram methods, is found for the inversion of Gegenbauer, Attractive, as well as Repulsive, Coulomb–Pollaczek, and Hermite weights and zeros to approximate for these OP systems on $x\in[-1,1]$ and $(-\infty,\infty)$ respectively, Reinhardt (2018), and Reinhardt (2021b), Reinhardt (2021a). Achieving precisions at this level shown above requires higher than normal computational precision, see Gautschi (2009). In Figure 18.40.2 the approximations were carried out with a precision of 50 decimal digits.

18.39 Applications in the Physical Sciences 18.41 Tables

Site Privacy Accessibility Privacy Program Copyrights Vulnerability Disclosure No Fear Act Policy FOIA Environmental Policy Scientific Integrity Information Quality Standards Commerce.gov Science.gov USA.gov