# Numerical Computation of Orthogonal Polynomials Recurrence Relations

Background and notations: Given an interval $I\subseteq \mathbb{R}$ and a continuous finite measure $d\mu = w(x)dx$, and denote $p_n(x)$ the orthogonal polynomials with respect to $d\mu$. We have the following recurrence relation $$p_j (x) = (a_j x +b_j)p_{j-1}(x) +c_jp_{j-2}(x) \, ,\quad \forall j\geq 1\, \quad p_0(x) \equiv 1,\,p_{-1}(x)\equiv 0 \, ,$$ where the constants are determined by the measure.

In the famous Golub-Welsch paper, section 4, they give a numerical method to calculate the constants. However, it requires the numerical value of $\int\limits_{I}x^\ell w(x)\,dx$ for all non-negative integers $\ell$.

The problem is that for numerical integration we usually need a quadrature formula, for which weights we need the recurrence relation (see the other sections of the same paper, for example).

Question 1: Is there a way to compute $\int\limits_{I}x^\ell w(x)\,dx$ without any quadrature formula?

Question 2: Is there a way to compute the recurrence constants without evaluating integrals?

Remark: We can always evaluate $\int\limits_{I}x^\ell w(x)\,dx$ using the Gauss Legendre quadrature. This means that we need only the well-known recurrence for the Lebesgue measure to compute these integrals for all other continuous measures. I'm looking for something else, though.

Computation of the coefficients of the recurrence relations of orthogonal polynomials is studied in details in the standart reference:

Gautschi, W., Orthogonal polynomials: computation and approximation. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2004

In particular, discretization methods are described in Section 2.2. They consist in the following :

1) approach the given measure $d\mu$ by a discrete $N$-point measure $d\mu_{N}$,

2) compute the recurrence coefficients of $d\mu_{N}$ and let $N$ go to infinity.

The three main issues, namely, appropriate choice of discretization, computation of recurrence coefficients of discrete measures, and convergence as $N$ tends to infinity, are discussed in Sections 2.2.1--2.2.4. Nontrivial examples are also given.

The following recent book should also be of interest :

Gautschi, W., Orthogonal polynomials in MATLAB. Exercises and solutions. Software, Environments, and Tools, 26. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2016.

I'm not sure if this is what you are looking for, but here we go. From orthogonality, you get $$\int_Ip_n(x)p_m(x)\,d\mu(x)=h_n\delta_{m,n}.$$ Suppose now that the leading coefficient of $p_n(x)$ is $k_n$. Then, $$a_n=\frac{k_{n+1}}{k_n} \qquad \text{and} \qquad -c_{n+1}=\frac{a_{n+1}h_{n+1}}{a_nh_n}.$$

• Thanks. A question and a remark: (1) What is $A_n$? (2) The problem with this method is that it is numerically unstable. $k_n$ are exponentially big in $n$, and so it is not stable compute their ration. Inasmuch as I gather, this is why the Golub-Welsch eigenvalue kind of methods were initially devised. – Amir Sagiv Nov 22 '16 at 14:52