Gauss–Legendre quadrature

Set $\omega(x)=\prod_{i=1}^n(x-x_i)$. The degree-$2n$ polynomial $\omega^2$ vanishes at every node, so its quadrature value is zero. Its integral is strictly positive because it is nonnegative and not identically zero. Thus the rule fails on this polynomial, regardless of the signs of its weights.

For necessity, $\omega q$ has degree at most $2n-1$ when $\deg q\leq n-1$. It vanishes at all nodes, so exactness forces its integral to vanish.

For sufficiency, divide any $p$ of degree at most $2n-1$ as $p=\omega q+r$, with $\deg q\leq n-1$ and $\deg r\leq n-1$. Orthogonality makes the integral of $\omega q$ zero, while its node values make its quadrature zero. Interpolatory weights integrate $r$ exactly. Hence they integrate $p$ exactly.

Rodrigues’ formula has degree $n$ and the displayed leading coefficient. Integrate its product with a lower-degree polynomial $q$ by parts $n$ times. The first $n-1$ derivatives of $(x^2-1)^n$ vanish at both endpoints, so all boundary terms vanish, and $q^{(n)}=0$ proves orthogonality. At $x=1$, differentiating $(x-1)^n(x+1)^n$ exactly $n$ times leaves the single nonzero term $n!2^n$, proving the normalization.

For roots, multiply $P_n$ by the product $q$ of the linear factors corresponding to all its sign-changing roots inside $(-1,1)$. Then $P_nq$ has one sign throughout the interval and a nonzero integral. If there were fewer than $n$ such roots, $\deg q<n$ would contradict orthogonality. There are therefore $n$ interior roots; the degree forces all to be simple. Two monic orthogonal degree-$n$ polynomials have a lower-degree difference orthogonal to itself, hence zero. Normalization gives uniqueness of $P_n$.

For the norm, integrate Rodrigues’ formula against $P_n$ by parts $n$ times:

Integrating $[x(1-x^2)^n]'$ gives $(2n+1)J_n=2nJ_{n-1}$. Starting from $J_0=2$, this yields

Substitution of $\sigma_n$ gives the norm.

Expand $xP_n$ in the orthogonal polynomial basis. For $k\leq n-2$, $\langle xP_n,P_k\rangle=\langle P_n,xP_k\rangle=0$. Rodrigues’ formula gives parity, so $\langle xP_n,P_n\rangle=0$. Only $P_{n+1}$ and $P_{n-1}$ remain. Comparing leading coefficients gives $(n+1)/(2n+1)$ for the first; evaluating at one gives $n/(2n+1)$ for the second. Rearrange.

Construct $H_{2n-1}$ matching the function and derivative at every node. Existence and uniqueness follow from finite-dimensional linear algebra: homogeneous conditions force a polynomial of degree at most $2n-1$ to be divisible by $\omega^2$, hence zero; the data map between two dimension-$2n$ spaces is invertible.

Gauss is exact on $H_{2n-1}$ and its node values equal those of $f$, so $Q_n(f)=\int H_{2n-1}$. At any non-node $x$, form

The $n$ double node zeros and the additional zero $x$ give $2n+1$ zeros. Repeated Rolle’s theorem gives $C=f^{(2n)}(\xi_x)/(2n)!$. At nodes the error is zero.

Bound the continuous derivative between its minimum and maximum, multiply by nonnegative $\omega^2$, and integrate. Divide by $\int\omega^2>0$ and apply the intermediate value theorem to obtain one common $\xi$. No continuity of $\xi_x$ is needed. Finally, $\omega=P_n/\sigma_n$, so

Substitute to conclude.

Multiply the reference-interval error coefficient by $(h/2)^{2n+1}$ for each panel. The local bound is $B_{2n}h^{2n+1}(n!)^4/((2n+1)[(2n)!]^3)$. Summing over $N=(b-a)/h$ panels proves the bound. Thus $n=1$ recovers second-order composite midpoint, and $n=2$ gives fourth-order global accuracy.

Multiply and subtract the three-term recurrences at $x$ and $y$, then sum over their indices. Adjacent terms telescope, giving the Christoffel–Darboux identity

In detail, multiply the index-$k$ recurrences by $P_k(y)$ and $P_k(x)$ respectively and subtract; adjacent terms have coefficients $k+1$ and $k$, so summation leaves only the last term.

Set $y=x_i$. The right side becomes $nP_{n-1}(x_i)P_n(x)/(2(x-x_i))$. Divide by its limit at $x_i$, namely $nP_{n-1}(x_i)P_n'(x_i)/2$, to get $\ell_i(x)$. The integral of the left side is one, since only $k=0$ is not orthogonal to constants. Hence

To prove the derivative identity, let $G=(1-x^2)P_n'+nxP_n$. Its highest-degree term cancels, so its degree is at most $n-1$. For every $q$ of degree at most $n-2$, integration by parts gives

Thus $G$ is a multiple of $P_{n-1}$. Evaluating at one gives multiplier $n$, so $G=nP_{n-1}$; the case $n=1$ is immediate directly. At a root, $(1-x_i^2)P_n'(x_i)=nP_{n-1}(x_i)$. Substitute into the weight formula.