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).

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}d\mu (x)$, $n=0,1,2,\mathrm{\dots}$, can these be used to find a unique
$\mu (x)$, a non-decreasing, real, function of $x$, 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).

In what follows we consider only the simple, illustrative, case that $\mu (x)$ is continuously differentiable so that $d\mu (x)=w(x)dx$, with $w(x)$ real, positive, and continuous on a real interval $[a,b].$ 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 $w(x)$, 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 $N$ be a positive integer and define

18.40.1 | ${P}_{1,1}$ | $=1,{P}_{m,1}=0,$ | ||

$m=2,\mathrm{\dots},2N+4$, | ||||

${P}_{m,2}$ | $={(-1)}^{m-1}{\mu}_{m-1},$ | |||

$m=1,\mathrm{\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\ge 3$, $m\ge 1$, $m+n\le 2N+5$, | ||||

use the first row of this $P$-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,\mathrm{\dots},2N+3$, | |||

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,\mathrm{\dots},N$. | |||

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 $w(x)$.

We have from (18.2.38) that

18.40.4 | $$\underset{N\to \mathrm{\infty}}{lim}{F}_{N}(z)=F(z)\equiv \frac{1}{{\mu}_{0}}{\int}_{a}^{b}\frac{w(x)dx}{z-x},$$ | ||

$z\in \u2102\backslash [a,b]$, | |||

in which

18.40.5 | $${F}_{N}(z)=\frac{1}{{\mu}_{0}}\sum _{n=1}^{N}\frac{{w}_{n}}{z-{x}_{n}}.$$ | ||

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

18.40.6 | $$\underset{\epsilon \to 0+}{lim}{\int}_{a}^{b}\frac{w(x)dx}{{x}^{\prime}+\mathrm{i}\epsilon -x}dx={\u2a0d}_{a}^{b}\frac{w(x)dx}{{x}^{\prime}-x}-\mathrm{i}\pi w({x}^{\prime}),$$ | ||

The question is then: how is this possible given only ${F}_{N}(z)$, rather than $F(z)$ itself?
${F}_{N}(z)$ often converges to smooth results for $z$ off the real axis for $\mathrm{\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 $w(x)$ are easily obtained for $N\sim 10$ to $20$. Gautschi (2004, p. 119–120) has explored the
$\epsilon \to {0}^{+}$ limit via the Wynn $\epsilon $-algorithm, (3.9.11) to accelerate convergence, finding four to eight digits of precision in $w(x)$,
depending smoothly on ${x}^{\prime}$, for $N\approx 4000$, for an example involving *first numerator* Legendre OP’s.

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)$, | |||

$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}d\mu (x)$, as shown in Figure 18.40.1 for $N=12$ and $120$, shown here for the repulsive Coulomb–Pollaczek OP’s of Figure 18.39.2, with the parameters as listed therein.

The bottom and top of the steps at the ${x}_{i}$ are lower and upper bounds to ${\int}_{a}^{{x}_{i}}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 $x$ yields a Stieltjes–Perron inversion to obtain ${w}^{\mathrm{RCP}}(x)$
to a precision of $\sim 4$ decimal digits for $N=120$. 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).

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}}{{\frac{dx(j,N)}{dj}|}_{j=i}}.$$ | ||

This allows Stieltjes–Perron inversion for the $w({x}_{i,N})$, given the quadrature weights and points. Here $x(t,N)$ is an interpolation of the abscissas ${x}_{i,N},i=1,2,\mathrm{\dots},N$, that is, $x(i,N)={x}_{i,N}$, allowing differentiation by $i$. 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)=\frac{{x}_{1,N}}{1+}\frac{{a}_{1}(t-1)}{1+}\frac{{a}_{2}(t-2)}{1+}\mathrm{\cdots}\frac{{a}_{N-1}(t-(N-1))}{1},$$ | ||

$t\in (0,\mathrm{\infty})$, | |||

where the coefficients are defined recursively via ${a}_{1}=\frac{{x}_{1,N}}{{x}_{2,N}}-1$, and

18.40.10 | $${a}_{\mathrm{\ell}}=-1-\frac{2{a}_{\mathrm{\ell}-1}}{1+}\frac{3{a}_{\mathrm{\ell}-2}}{1+}\mathrm{\cdots}\frac{\mathrm{\ell}{a}_{1}}{1-\frac{{x}_{1,N}}{{x}_{\mathrm{\ell}+1,N}}},$$ | ||

$\mathrm{\ell}=2,\mathrm{\dots},N-1$. | |||

The PWCF $x(t,N)$ is a minimally oscillatory algebraic interpolation of the abscissas ${x}_{i,N},i=1,2,\mathrm{\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 $[-1,1]$ has an essential singularity at $x=-1$.

Further, *exponential* convergence in $N$, 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 $w(x)$ for these OP systems
on $x\in [-1,1]$ and $(-\mathrm{\infty},\mathrm{\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.