$$\int_{-1}^{1}f(x)dx \approx c_1f(x_1)+c_2f(x_2)+c_3f(x_3)\qquad(*)$$ for three nodes
We have three weights, $c_1,c_2,c_3$ and three nodes ,$x_1,x_2,x_3$ to find. Now choose $$f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5$$
$(2\times3-1)=5$ degree polynomials. Integrate $(*)$ from $-1\text{ to }1$ we found
$$\int_{-1}^1f(x)dx=a_0(2)+a_1(0)+a_2\left(\frac{2}{3}\right)+a_3(0)+a_4\left(\frac{2}{3}\right)+a_5(0)\qquad(1)$$
Now, $$c_1f(x_1)+c_2f(x_2)+c_3f(x_3)=c_1(a_0+a_1x_1+a_2x^2_1+a_3x^3_1+a_4x^4_1+a_5x^5_1)+c_2(a_0+a_1x_2+a_2x^2_2+a_3x^3_2+a_4x^4_2+a_5x^5_2)+c_3(a_0+a_1x_3+a_2x^2_3+a_3x^3_3+a_4x^4_3+a_5x^5_3)\qquad(2)$$
One last work equating $(1)\text{ and }(2)$ and collect those coefficients.
$$\begin{eqnarray}
c_1+c_2+c_3 &=& 2 \\
c_1x_1+c_2x_2+c_3x_3 &=& 0 \\
c_1x^2_1+c_2x^2_2+c_3x^2_3 &=&\frac{2}{3} \\
c_1x^3_1+c_2x^3_2+c_3x^3_3 &=& 0 \\
c_1x^4_1+c_2x^4_2+c_3x^4_3 &=& \frac{2}{5} \\
c_1x^5_1+c_2x^5_2+c_3x^5_3 &=& 0
\end{eqnarray}$$
I'm not sure how to get the rest because these are non linear equation then how to get $c_1,c_2,c_3,x_1,x_2,x_3$. According to Wikipedia: $$\int_{-1}^{1}f(x)dx \approx \frac{8}{9}f(0)+\frac{5}{9}f\left(\sqrt{\frac{3}{5}}\right)+\frac{5}{9}f\left(-\sqrt{\frac{3}{5}}\right)$$
One more question: What's the geometrical intuition of Gaussian quadrature$?$ Why it works like this$?$ Because without knowing it seem I am doing thing not learning. And if there is any better way to solve nonlinear equations like this please provide also. Sorry for asking too many questions in one.
Thanks in Advance.
Best Answer
To know why Gauss quadrature works, you should look at the proof.
As for how to do it, you need to do Gram-Schmidt on the standard polynomial basis to get a degree three orthogonal polynomial. The roots of this polynomial will constitute the nodes. Then, test your rule on the polynomials $1,x$, and $x^2$ to get a linear system of equations to find the weights. The proof explains why this is what you want to do, and this is the algorithm that you always follow for such a problem.
EDIT: How it's done-
First, we apply Gram-Schmidt to the standard basis $v_1=1,v_2=x,v_3=x^2,v_4=x^3$ with respect to the inner product $\langle f,g\rangle =\int\limits_{-1}^1 f(x)g(x)\, dx$. We go up to a cubic because the roots will be the nodes, and we want three nodes. In particular, we want to compute \begin{align*} u_1&=v_1\\ v_2&=u_2-\frac{\langle u_2,v_1\rangle}{\langle v_1,v_1\rangle}v_1\\ v_3&=u_3-\frac{\langle u_3,v_1\rangle}{\langle v_1,v_1\rangle}v_1-\frac{\langle u_3,v_2\rangle}{\langle v_2,v_2\rangle}v_2\\ v_4&=u_4-\frac{\langle u_4,v_1\rangle}{\langle v_1,v_1\rangle}v_1-\frac{\langle u_4,v_2\rangle}{\langle v_2,v_2\rangle}v_2-\frac{\langle u_4,v_3\rangle}{\langle v_3,v_3\rangle}v_3 \end{align*} You can find that \begin{align*} u_1&=1\\ u_2&=x\\ u_3&=x^2-\frac{1}{3}\\ u_4&=x^3-\frac{3}{5}x \end{align*} The roots of $u_4$ are $x_{1,2,3}=0,\pm\sqrt{\frac{3}{5}},$ so these are the nodes.
So, our rule is $$\int\limits_{-1}^1f(x)\, dx=c_1f(0)+c_2f\left(\sqrt{\frac{3}{5}}\right)+c_3f\left(-\sqrt{\frac{3}{5}}\right),$$ and all that we have to do is find $c_1,c_2,c_3.$ We impose exactness on $f(x)=1,x,x^2.$ That is, we require that the rule is an equality for these $f$. If we compute what this means, we get that \begin{align*} \text{ if }f(x)=1&:\ \ \int\limits_{-1}^1 \, dx=2=c_1+c_2+c_3\\ \text{ if }f(x)=x&:\ \ \int\limits_{-1}^1x \, dx=0=c_2\sqrt{\frac{3}{5}}-c_3\sqrt{\frac{3}{5}}\\ \text{ if }f(x)=x^2&:\ \ \int\limits_{-1}^1x^2 \, dx=\frac{2}{3}=c_2\frac{3}{5}+c_3\frac{3}{5} \end{align*} This is a system of linear equations, and you can solve it and find that
\begin{align*} c_1&=\frac{8}{9}\\ c_2&=\frac{5}{9}\\ c_3&=\frac{5}{9}\\ \end{align*} This is consistent with the rule that you were looking for.
You can find a proof of the general method of Gaussian quadrature here: https://www.johndcook.com/OrthogonalPolynomials.pdf