About the Project
18 Orthogonal PolynomialsComputation

§18.40 Methods of Computation

Contents
  1. §18.40(i) Computation of Polynomials
  2. §18.40(ii) The Classical Moment Problem

§18.40(i) Computation of Polynomials

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

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, μn=abxndμ(x), n=0,1,2,, can these be used to find a unique μ(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).

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

In what follows we consider only the simple, illustrative, case that μ(x) is continuously differentiable so that dμ(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 αn,βn of equations (18.2.11_5) and (18.2.11_8); then, 2) to construct the quadrature abscissas xi and weights (or Christoffel numbers) wi 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 αn and β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 P1,1 =1,Pm,1=0,
m=2,,2N+4,
Pm,2 =(1)m1μm1,
m=1,,2N+3,
Pm,n =P1,n1Pm+1,n2P1,n2Pm+1,n1,
n3, m1, m+n2N+5,

use the first row of this P-matrix for

18.40.2 a1=μ0,an=P1,n+1P1,nP1,n1,
n=2,,2N+3,

and these can be used for the recursion coefficients

18.40.3 α0=a2,αn=a2n+1+a2n+2,βn=a2na2n+1,
n=1,,N.

The quadrature abscissas xn and weights wn 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).

Stieltjes Inversion via (approximate) Analytic Continuation

We have from (18.2.38) that

18.40.4 limNFN(z)=F(z)1μ0abw(x)dxzx,
z\[a,b],

in which

18.40.5 FN(z)=1μ0n=1Nwnzxn.

Let x(a,b). It is now necessary to take the limit ε0+ of F(x+iε), and the imaginary part is the required Stieltjes–Perron inversion:

18.40.6 limε0+abw(x)dxx+iεxdx=abw(x)dxxxiπw(x),

The question is then: how is this possible given only FN(z), rather than F(z) itself? FN(z) often converges to smooth results for z off the real axis for z at a distance greater than the pole spacing of the xn, 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 FN(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 N10 to 20. Gautschi (2004, p. 119–120) has explored the ε0+ limit via the Wynn ε-algorithm, (3.9.11) to accelerate convergence, finding four to eight digits of precision in w(x), depending smoothly on x, for N4000, for an example involving first numerator Legendre OP’s.

Histogram Approach

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

18.40.7 μN(x)=n=1NwnH(xxn),
x(a,b),

H(x) being the Heaviside step-function, see (1.16.13).

Equation (18.40.7) provides step-histogram approximations to axdμ(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.

See accompanying text
Figure 18.40.1: Histogram approximations to the Repulsive Coulomb–Pollaczek, RCP, weight function integrated over [1,x), see Figure 18.39.2 for an exact result, for Z=+1, shown for N=12 and N=120. Magnify

The bottom and top of the steps at the xi are lower and upper bounds to axidμ(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 wRCP(x) to a precision of 4 decimal digits for N=120. Convergence is O(N2). 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 μ(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

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(xi,N)wi,Ndx(j,N)dj|j=i.

This allows Stieltjes–Perron inversion for the w(xi,N), given the quadrature weights and points. Here x(t,N) is an interpolation of the abscissas xi,N,i=1,2,,N, that is, x(i,N)=xi,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)=x1,N1+a1(t1)1+a2(t2)1+aN1(t(N1))1,
t(0,),

where the coefficients are defined recursively via a1=x1,Nx2,N1, and

18.40.10 a=12a11+3a21+a11x1,Nx+1,N,
=2,,N1.

The PWCF x(t,N) is a minimally oscillatory algebraic interpolation of the abscissas xi,N,i=1,2,,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 αn,βn for the weight function for the repulsive Coulomb–Pollaczek, RCP, polynomials of (18.39.50). This is a challenging case as the desired wRCP(x) on [1,1] has an essential singularity at x=1.

See accompanying text
Figure 18.40.2: Derivative Rule inversions for wRCP(x) carried out via Lagrange and PWCF interpolations. Shown are the absolute errors of approximation (18.40.8) at the points xi,N, i=1,2,,N for N=40. For the derivative rule Lagrange interpolation (red points) gives 15 digits in the central region, while PWCF interpolation (blue points) gives 25. Magnify

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[1,1] and (,) 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.