If $\left( p\_n \right)\_{n=0}^{\infty}$ is a family of orthogonal polynoamials with respect to a measure $\mu$ on $[-1,1]$, and $\left( x\_j, w\_j \right)$ are the quadrature points and weights for the respective Gaussian quadrature rule, we can easily prove that 
$$ f\_N (x) : \, = \sum\limits\_{n=0}^{N-1} \sum\limits\_{j=1}^{N} p\_n (x\_j) f(x\_j) w\_j p\_n (x) \, ,$$
is the interpolation polynomial of degree $N-1$ for $f$ at the quadrature points $x\_1, \ldots x\_N$.


**My question:** While I could prove it, I couldn't find a reference for this proof in any textbook. Could you help with that?


This is cross-posted from [this](https://math.stackexchange.com/q/2595490/161591) post in MSE.



Proof Reference - Polynomial interpolation at quadrature points
Here, three possible references for the formula:



> 
> P. J. Davis and P. Rabinowitz, Methods of Numerical Integration,
>  Computer Science and Applied Mathematics. Academic Press, New York,
>  1984 (see p.88)
> 
> 
> J. C. Mason and D. C. Handscomb, Chebyshev Polynomials, CRC Press, New
>  York, 2003 (see section 8.3.2 for the case of Chebyshev Polynomials)
> 
> 
> J. Shen, T. Tang, L-L. Wang, Spectral methods. Algorithms, analysis
>  and applications. Springer Series in Computational Mathematics, 41.
>  Springer, Heidelberg, 2011 (see Thm 3.9 p.64)
> 
> 
> 


