We would like to approximate the integral of $f$ over an interval $[a,b]$ with positive measure $\mu$. To this end, we choose $n$ points $x_1, x_2, \ldots, x_n$ inside the open interval $(a,b)$ together with the end points $a$ and $b$ and we choose $n+2$ weights $\lambda_a, \lambda_b, \lambda_{x_1}, \lambda_{x_2}, \ldots, \lambda_{x_n}$ in such a way that we hope that the formula (called the Gauss-Lobatta quadrature formula) $\displaystyle \lambda_af(a) + \lambda_bf(b) + \sum_{k=1}^n \lambda_{x_k}f(x_k)$ is a good approximation to $\displaystyle \int_a^b f(x) d\mu(x)$. If this formula is exact whenever $f$ is a polynomial of degree at most $m$, we say that our formula has degree of accuracy $m$.
Polynomials $p_0, p_1, \ldots, p_n, \ldots$ (where $p_n$ has degree $n$) are called orthogonal polynomials for the measure $\mu$ if $\displaystyle \int_a^b p_i(x)p_j(x) d\mu(x) = 0$ for $i \neq j$.
I would like to prove the following: our formula has degree of accuracy $2n+1$ if, and only if, $x_1, x_2, \ldots, x_n$ are the zeros of the orthogonal polynomial of degree $n$ for the modified measure $(x-a)(x-b)d\mu(x)$.
Any help would be much appreciated!
Let $f$ be a polynomial of degree at most $2n+1$ and let $L$ be the Lagrange interpolant for $f$ at the interpolation points $a, x_1, \ldots, x_n, b$, where the $x_i$ are the zeros of the orthogonal polynomial $p_n$ for the measure $\mu_0$ with $d\mu_0(x) = (x-a)(x-b)d\mu(x)$. We can write
$$L = \displaystyle f(a)\ell_a(x) + f(b)\ell_b(x) + \sum_{j=1}^n f(x_j)\ell_{x_j}(x)$$
where $\ell_a, \ell_b$ and $\ell_{x_j}$ are the fundamental polynomials of Lagrange interpolation Then $f - L$ is a polynomial of degree at most $2n+1$ with zeros at the interpolation points, so we can write it as $p_n(x)(x-a)(x-b)q(x)$, where $q(x)$ is a polynomial of degree at most $n-1$. Integrating this and using the orthogonality relation we obtain:
$$\int_a^b f(x) - L(x) d\mu(x) = \int_a^b p_n(x)(x-a)(x-b)q(x) d\mu(x) \\ = \int_a^b p_n(x)q(x) d\mu_0(x) = 0$$
And we conclude that with $\lambda_a = \int_a^b \ell_a(x) d\mu(x)$ with similar expressions for the other weights, our formula has degree of accuracy $2n+1$.
While on the other hand, if the $x_i$ are not the zeros of $p_n(x)$, then for $p(x) = (x-x_1)\cdots(x-x_n)$ there exists a polynomial $q$ with degree at most $n-1$ such that $\int_a^b p(x)q(x) d\mu_0(x) \neq 0$.
Note that degree of accuracy $2n+2$ cannot be achieved if the support of $\mu$ contains more than $n$ points, since for $q(x) = (x-x_1)\cdots(x-x_n)$ (with the $x_i$ the zeros of $p_n(x)$) we have that $q_0(x) = q(x)^2(x-a)(x-b)$ is a polynomial of degree $2n+2$ with
$$\displaystyle \int_a^b q_0(x)\mu(x) = \int_a^b q(x)^2d\mu_0(x) > 0$$
while $\displaystyle \lambda_aq_0(a) + \lambda_bq_0(b) + \sum_{k=1}^n \lambda_{x_k}q_0(x_k) = 0$, since we evaluate $q_0$ at its zeros in that sum.