Why are quadrature points given by the zeros of orthogonal polynomials? We know there exists unique Gaussian quadrature formula. Its quadrature points are given by the zeros of the orthogonal polynomial. Why do we use only the zeros of the orthogonal polynomials in Gaussian quadrature formula? Is it possible to use any points instead of the zeros of the orthogonal polynomials?  
 A: You can use any points you wish in a quadrature formula; it just will not be very precise if you don't choose them wisely. For example, 
$$
\int_{-1}^1 f(x)\,dx \approx  7 f(-1) - 10 f(1/4) + 5 f(1)
$$
is a quadrature formula with points $-1,1/4,1$. It's not a very precise formula, because I just made up the numbers. If I put more effort into the coefficients, I would make this formula exact  for polynomials of degrees up to $2$. But it still wouldn't work as well as the Gaussian quadrature on three points ($0$ and $\pm \sqrt{3/5}$), which is accurate for degrees up to $5$.
Here is why the zeros of orthogonal polynomials (in particular, Legendre polynomials) make a good choice of sample points. Given a general polynomial $p$ of degree $2n-1$, we can use long division to obtain $p=qL_n+r$ where $L_n$ is the Legendre polynomial of degree $n$, and both $q$ and $r$ have degree $<n$. Then
$$
\int_{-1}^1 p(x)\,dx = \int_{-1}^1 q(x)L_n(x)\,dx+\int_{-1}^1 r(x)\,dx
= \int_{-1}^1 r(x)\,dx
$$
because $L_n$ is orthogonal to any polynomial of degree less than $n$. Since the term $qL_n$ contributes zero to the integral, it should also contribute zero to our quadrature formula. Choosing the sample points to be the zeros of $L_n$ achieves this goal. 
Thus, one only needs to find the weights that give the exact value of $\int_{-1}^1 r(x)\,dx$ (where $\deg r\le n-1$) to obtain a quadrature formula that is exact for all degrees up to $2n-1$.
