Two-point Gaussian Quadrature Rule with weight function $w(x)=x$ I am trying to construct a formula of the form $\int_0^1 xf(x)dx=A_0f(x_0)+A_1f(x_1)$ with degree of precision 3. 
A Gaussian Quadrature is the only interpolatory quadrature with degree of precision 2n+1 with $x_0,x_1,...,x_n$ being the nodes chosen that are also the $n+1$ zeros of the $(n+1)$th orthogonal polynomial. 
Looking through similar questions, most recommend using the method of undetermined coefficients to find $x_0,x_1,A_0,$ and $A_1$ by using $f(x) = 1,x,x^2,$ and $x^3$. Does using the method of undetermined coefficients inherently give you nodes that are the $n+1$ zeros of the $(n+1)$th orthogonal polynomial? 
For $w(x)=1$, we can use the Legendre orthogonal polynomials to get the nodes, but from what I could find online there doesn't exist such a thing for $w(x)=x$. 
 A: Yes, the method of undetermined coefficients does give you the the nodes. We want
$$\int_0^1x\cdot x^ndx=A_0x_0^n+A_1x_1^n=\frac1{n+2}$$
for $x\in\{0,1,2,3\}$' We can write this out as a matrix equation
$$\begin{bmatrix}1&1\\x_0&x_1\\x_0^2&x_1^2\\x_0^3&x_1^3\end{bmatrix}\begin{bmatrix}A_0\\A_1\end{bmatrix}=\begin{bmatrix}\frac12\\\frac13\\\frac14\\\frac15\end{bmatrix}$$
To solve this vanderMonde system we can do Gaussian elimination from the bottom up, subtracting $x_0$ times the third row from the fourth, then $x_0$ times the second row from the third and then $x_0$ times the first row from the second to get
$$\begin{bmatrix}1&1\\0&x_1-x_0\\0&x_1(x_1-x_0)\\0&x_1^2(x_1-x_0)\end{bmatrix}\begin{bmatrix}A_0\\A_1\end{bmatrix}=\begin{bmatrix}\frac12\\\frac13-\frac12x_0\\\frac14-\frac13x_0\\\frac15-\frac14x_0\end{bmatrix}$$
Another round yields
$$\begin{bmatrix}1&1\\0&x_1-x_0\\0&0\\0&0\end{bmatrix}\begin{bmatrix}A_0\\A_1\end{bmatrix}=\begin{bmatrix}\frac12\\\frac13-\frac12x_0\\\frac14-\frac13(x_0+x_1)+\frac12x_0x_1\\\frac15-\frac14(x_0+x_1)+\frac13x_0x_1\end{bmatrix}$$
The last two rows are linear equations for the coefficients of the orthogonal polynomial
$$\begin{bmatrix}\frac13&-\frac12\\\frac14&-\frac13\end{bmatrix}\begin{bmatrix}x_0+x_1\\x_0x_1\end{bmatrix}=\begin{bmatrix}\frac14\\\frac15\end{bmatrix}$$
With solutions $x_0+x_1=6/5$, $x_0x_1=3/10$. So the equation for the nodes is
$$x^2-(x_0+x_1)x+x_0x_1=x^2-\frac65x+\frac3{10}=0$$
With solutions
$$x_0=\frac{6-\sqrt6}{10}$$
$$x_1=\frac{6+\sqrt6}{10}$$
Then from the second row we have
$$A_1=\frac{\frac13-\frac12x_0}{x_1-x_0}=\frac{3\sqrt6+2}{12\sqrt6}$$
And from the first row
$$A_0=\frac12-A_1=\frac{3\sqrt6-2}{12\sqrt6}$$
There is in fact a family of orthogonal polynomials appropriate to this problem, the Jacobi polynomials with $\alpha=0$ and $\beta=1$. You could have worked them out yourself by assuming $\pi(x)=(x-x_0)(x-x_1)=x^2+ax+b$ and used orthogonality to $f(x)=1$ and $f(x)=x$ to get $a$ and $b$.
A: Given $x_0, x_1,$ with $0 < x_0 < x_1 < 1,$ we can write any linear
polynomial as
$$
f(x) = \frac{x - x_1}{x_0 - x_1}f(x_0) +
\frac{x - x_0}{x_1 - x_0}f(x_1),
$$
where $f(x_0)$ and $f(x_1)$ are arbitrary, so we must have
\begin{align*}
A_0 & = \int_0^1x\frac{x - x_1}{x_0 - x_1}\,dx
= \frac{2 - 3x_1}{6(x_0 - x_1)}, \\
A_1 & = \int_0^1x\frac{x - x_0}{x_1 - x_0}\,dx
= \frac{2 - 3x_0}{6(x_1 - x_0)}.
\end{align*}
Lagrange interpolation similarly gives the $A_i$ in terms of the
$x_i$ for any weight function $w(x),$ and any number of points.
There is also a general construction of the $x_i$ as the zeros of
one of a series of polynomials $\phi_0, \phi_1, \phi_2, \ldots$
orthogonal with respect to an inner product and norm determined by
$w(x).$ In this case,
\begin{gather*}
(f, g) = \int_0^1xf(x)g(x)\,dx, \\
\|f\| = \sqrt{(f, f)} = \sqrt{\int_0^1xf(x)^2\,dx}.
\end{gather*}
The general proof, for weight function $w(x)$ and $k+1$ points
$x_0, x_1, \ldots, x_k,$ is short enough to give here. (See Theorem
12.3 in M. J. D. Powell, Approximation theory and methods,
where the calculation for the case $w(x) = x,$ $k = 1$ is set as
Exercise 12.2.) Any polynomial $f$ of degree $2k+1$ or less may be
written as $p(x)\phi_{k+1}(x) + q(x),$ where $p, q$ are of degree
$k$ or less.  By orthogonality (functions are defined on $[a, b]$),
$$
\int_a^bw(x)f(x)\,dx = \int_a^bw(x)q(x)\,dx.
$$
If $x_0, x_1, \ldots, x_k$ are [the] zeros of $\phi_{k+1},$ then,
for any weights $A_i,$
$$
\sum_{i=0}^kA_if(x_i) = \sum_{i=0}^kA_iq(x_i).
$$
But by Lagrange interpolation,
$$
\int_a^bw(x)q(x)\,dx = \sum_{i=0}^kA_iq(x_i),
$$
where the weights $A_i$ are now given by
$$
A_i =
\int_a^bw(x)\bigg(\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\bigg)\,dx.
$$
Therefore
$$
\int_a^bw(x)f(x)\,dx = \sum_{i=0}^kA_if(x_i),
$$
as required. $\square$
The problem now reduces to obtaining
$\phi_0, \phi_1, \phi_2, \ldots$ by a Gram-Schmidt process, i.e.
iteratively taking a polynomial of degree $j+1$
($j = 0, 1, 2, \ldots$) and subtracting from it its projection onto
the space of polynomials of degree $j$ or less.  Instead of doing
this with the polynomial $x^{j+1},$ we recursively use $x\phi_j(x).$
The proof is a little too long to give here, but the result is the
recurrence relation
\begin{gather*}
\phi_0(x) = 1, \\
\phi_1(x) = x - \alpha_0, \\
\phi_{j+1}(x) = (x - \alpha_j)\phi_j(x) - \beta_j\phi_{j-1}(x)
\quad (j \geqslant 1),
\end{gather*}
where
\begin{gather*}
\alpha_j = \frac{(\phi_j, x\phi_j)}{\|\phi_j\|^2}
\quad (j \geqslant 0), \\
\beta_j = \frac{\|\phi_j\|^2}{\|\phi_{j-1}\|^2}
\quad (j \geqslant 1).
\end{gather*}
In the present case, with $w(x) = x,$ and $k = 1,$ we wish to find
$\phi_2.$ Here is a summary of the calculation:
\begin{gather*}
\|\phi_0\|^2 = \int_0^1x\,dx = \frac12, \\
(\phi_0, x\phi_0) = \int_0^1x^2\,dx = \frac13, \\
\therefore \alpha_0 = \frac23, \ \phi_1(x) = x - \frac23; \\
\|\phi_1\|^2 = \int_0^1x\left(x - \frac23\right)^2\,dx
= \frac1{36}, \\
(\phi_1, x\phi_1) = \int_0^1x^2\left(x - \frac23\right)^2\,dx
= \frac2{135}, \\
\therefore\ \alpha_1 = \frac8{15}, \ \beta_1 = \frac1{18}, \\
\phi_2(x) =
\left(x - \frac8{15}\right)\left(x - \frac23\right) - \frac1{18}
= \boxed{x^2 - \frac{6x}5 + \frac3{10}},
\end{gather*}
giving
\begin{align*}
x_0 = \frac{6 - \sqrt6}{10}, \\
x_1 = \frac{6 + \sqrt6}{10},
\end{align*}
and
\begin{align*}
A_0 & = \frac{2 - 3x_1}{6(x_0 - x_1)} = \frac{9 - \sqrt6}{36}, \\
A_1 & = \frac{2 - 3x_0}{6(x_1 - x_0)} = \frac{9 + \sqrt6}{36},
\end{align*}
in agreement with user5713492's answer.
