§18.40 Methods of Computation

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

use the first row of this $P$-matrix for

and these can be used for the recursion coefficients

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

in which

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:

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

$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

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:

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

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.