Prerequisite definite integral
Gauss quadrature can be approximated by summation
\begin{equation}
\int f(x) \,\mathrm{d}{x} \approx \sum_i w_i f(x_i)
\end{equation}
For a set of $x_i, w_i$ in a certain interval, if it is of degree $N$, then when $f(x)$ is a polynomial of degree less than or equal to $N$, the above formula is exactly true.
Gauss-Lobatto points
In Gauss-Lobatto integration, let the number of sampling points (including two end points) be $N$. If the integrand is a polynomial of order $2N-3$ ($f(x) = x^{2n-3} + \dots$), there is no error in the integration.
Note that the Gauss-Lobatto integral is symmetric
\begin{equation}
x_j = -x_{N-j+1} \qquad w_{0j} = w_{0,N-j+1}
\end{equation}
Orthogonal normalized basis
Each base is a polynomial of order $N-1$. Since the order is the same as the number of zeros, the polynomial can be uniquely determined, that is, the Lagrange interpolation polynomial
\begin{equation}
\begin{aligned}
p_n(x) &= \prod_{i=1}^{n-1} \frac{x-x_i}{x_n-x_i} \prod_{i=n+1}^{N} \frac{x-x_i}{x_n-x_i}\\
&= \frac{x-x_1}{x_n-x_1} \times\dots\times\frac{x-x_{n-1}}{x_n-x_{n-1}}\frac{x-x-{n+1}}{x_n-x_{n+1}} \dots \frac{x-x_N}{x_n-x_N}
\end{aligned} \end{equation}
\begin{equation}
p_n(x_{n'}) = \delta_{n, n'}
\end{equation}
Since any two base products are only polynomials of order $2N-2$, the sum can be used to replace their integrals.
\begin{equation}
\int_{-1}^1 p_i(x) p_j(x) \,\mathrm{d}{x} = \sum_k w_k p_i(x_k) p_j(x_k) = w_i \delta_{ij}
\end{equation}
So $N$ orthogonal normalization basis is
\begin{equation}
f_n(x) = \frac{1}{\sqrt{w_n}} p_n(x)
\end{equation}
Satisfy
\begin{equation}
f_i(x_j) = \frac{1}{\sqrt{w_i}} \delta_{ij}
\end{equation}
Another equivalent representation method is to use the polynomial corresponding to the Gauss-Lobatto numerical integration of order $N$, the derivative of Legendre polynomial of order $N-1$, $P'_{N-1}(y)$, to construct a polynomial that satisfies the conditions. According to the definition, its $N-2$ zeros are respectively $y_2, y_3\dots y_{N-1}$. In order to add the two zeros of $y_1=-1$ and $y_n=1$, they are transformed into a $N$ degree polynomial
\begin{equation}
(1-y^2)P'_{N-1}(y)
\end{equation}
However,
eq. 2 requires $p_n(y_n)=1$, so we divide
eq. 8 by its own tangent at $y_n$ and form the limit type $0/0=1$ at $y=y_n$ to get the order polynomial $u_n(y)$.
\begin{equation}
u_n(y) = \frac{(1-y^2)P'_{N-1}(y)}{[(1-y^2)P'_{N-1}(y)]'_y} = y_n (y-y_n)
\end{equation}
This formula is actually the same polynomial as
eq. 3 , because all polynomials of order $N-1$ with $N-1$ zeros can be factored into
eq. 3 multiplied by an undetermined constant. Use
eq. 9 to quickly expand polynomials (because the coefficients of Legendre polynomials can be directly calculated by formulas).