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}.$$

  • $\begingroup$ 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. $\endgroup$ – Amir Sagiv Nov 22 '16 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.