Orthogonal Polynomials: General Theory

Orthogonal polynomials are the polynomial eigenbasis that quietly organises an astonishing amount of mathematical physics: separation of variables in the Laplacian, the harmonic-oscillator and hydrogen wavefunctions, multipole expansions, and high-accuracy numerical integration all reduce to one of a small number of classical families. The unifying datum is a single object — a weight function $w(x)\ge 0$ on an interval $[a,b]$ with finite moments $\mu_k=\int_a^b x^k w\,\mathrm{d}x$ . Feed the monomials $1,x,x^2,\dots$ through Gram-Schmidt in the inner product $\langle f,g\rangle=\int_a^b fg\,w\,\mathrm{d}x$ and out comes a unique (up to scale) sequence $p_0,p_1,p_2,\dots$ with $\deg p_n=n$ , pairwise orthogonal by construction. From that one construction three apparently unrelated structures emerge and turn out to be equivalent. First, a universal three-term recurrence: because $x p_n$ has degree $n+1$ and orthogonality annihilates every overlap except with $p_{n+1},p_n,p_{n-1}$ , every family obeys a tridiagonal recurrence — orthogonal polynomials are the eigen-directions of a Jacobi matrix. Second, a Sturm-Liouville origin: the classical families are eigenfunctions of a self-adjoint second-order operator whose weight is forced by the Pearson equation $(ws)'=wt$ . Third, a Rodrigues formula, a closed form that produces $p_n$ by $n$ -fold differentiation of $w s^n$ . Threaded through all three are the zeros: real, simple, confined to $(a,b)$ , interlacing from one degree to the next — and, remarkably, the optimal nodes of Gaussian quadrature. This article develops the general theory; Legendre, associated Legendre, and spherical harmonics are treated in their own articles as worked instances of it.

Key takeaways

A family $\{p_n\}_{n\ge0}$ with $\deg p_n=n$ is fixed by one hypothesis: orthogonality $\int_a^b p_m p_n\,w\,\mathrm{d}x=h_n\delta_{mn}$ against a positive weight $w$ on $[a,b]$ .
Every such family obeys a universal three-term recurrence $p_{n+1}=(A_n x+B_n)p_n-C_n p_{n-1}$ , so an orthogonal family is exactly a symmetric tridiagonal (Jacobi) operator in disguise.
The classical families — Hermite, Laguerre, Jacobi (with Legendre, Chebyshev, Gegenbauer as special cases) — are eigenfunctions of a singular Sturm-Liouville operator and carry a Rodrigues formula $p_n\propto w^{-1}(\mathrm{d}/\mathrm{d}x)^n[w\,s^n]$ .
The $n$ zeros of $p_n$ are real, simple, interior to $(a,b)$ , and interlace those of $p_{n+1}$ ; they are precisely the nodes of the $n$ -point Gaussian rule, which is exact to degree $2n-1$ .
Favard's theorem is the converse: any three-term recurrence with positive coefficients reconstructs a positive weight, so "orthogonal family" and "tridiagonal recurrence" name the same object.
Key idea: Orthogonality against a positive weight $w$ on $[a,b]$ is the single hypothesis; the three-term recurrence, the real interlacing zeros, Gaussian quadrature, and (for classical weights) a Sturm-Liouville eigenproblem and Rodrigues formula are all forced consequences — and Favard's converse closes the loop by rebuilding the weight from the recurrence.

The core insight

Orthogonality against a positive weight $w$ on $[a,b]$ is the single hypothesis; everything else — the three-term recurrence, real interlacing zeros, Gaussian quadrature, and (for classical weights) a Sturm-Liouville eigenvalue problem and a Rodrigues formula — is a forced consequence, and Favard's theorem runs the implication backward, so an orthogonal family and a positive-coefficient tridiagonal recurrence are two faces of one object.

Where this lives today

Orthogonal polynomials are not a museum piece; they are load-bearing in modern computation and physics. Every serious numerical-integration library computes Gauss nodes and weights on the fly by the Golub-Welsch algorithm, which finds them as eigenvalues and first-eigenvector components of the symmetric tridiagonal Jacobi matrix built from the recurrence coefficients — the three-term recurrence turned into linear algebra. Chebyshev polynomials are the workhorse of spectral methods for PDEs, where their minimax property and FFT-friendly node placement give exponentially convergent solvers used in climate models, fluid dynamics, and quantum chemistry. In quantum many-body physics the kernel-polynomial method expands spectral densities and Green's functions in Chebyshev series to reach system sizes direct diagonalisation cannot touch, and the Lanczos algorithm — the backbone of large-scale eigensolvers and of matrix-product-state simulations — is exactly the construction of orthogonal polynomials with respect to a spectral measure. Meanwhile the theory keeps growing: multiple orthogonal polynomials and matrix-valued families now drive results in random-matrix theory and integrable probability, where the distribution of the largest eigenvalue of a random matrix is governed by the very same Christoffel-Darboux kernel that appears in equation (5).

The physical picture

Think of the weighted inner product $\langle f,g\rangle=\int_a^b fg\,w\,\mathrm{d}x$ as turning the space of polynomials into a Hilbert space, and of Gram-Schmidt as the tool that straightens the natural but oblique basis $1,x,x^2,\dots$ into a rigid orthogonal frame. The weight $w$ is the whole story: it says where on $[a,b]$ agreement is expensive. Hermite's $e^{-x^2}$ concentrates measure near the origin, Chebyshev's $(1-x^2)^{-1/2}$ piles it against the endpoints, Legendre's flat $w=1$ treats every point alike. Change the weight and you change the family — but not the skeleton of theorems every family shares.

Figure

The first Legendre polynomials P₁ blue, P₂ amber, P₃ green on the interval −1 to 1: Gram-Schmidt straightens the monomials 1, x, x², … into an orthogonal basis obeying ∫ Pₘ Pₙ dx = 2/(2n+1) · δₘₙ, and Pₙ has exactly n simple zeros interior to the interval

The three-term recurrence is the first non-obvious payoff, and it is no accident. Multiplying $p_n$ by $x$ raises the degree by one, and because multiplication by $x$ is symmetric — $\langle x p_n,p_k\rangle=\langle p_n,x p_k\rangle$ — the overlap with $p_k$ vanishes the moment $\deg(x p_k)<n$ , i.e. for every $k<n-1$ . Only the immediate neighbours $p_{n+1},p_n,p_{n-1}$ can appear. Written as a matrix, "multiply by $x$ " is therefore symmetric and tridiagonal — a Jacobi matrix — and the zeros of $p_n$ are literally its eigenvalues. That single fact makes orthogonal polynomials the discrete-spectral backbone of everything from Gaussian quadrature to the Lanczos algorithm. Those zeros are where the theory becomes visual. All $n$ of them lie strictly inside $(a,b)$ , all are simple, and the zeros of consecutive degrees interlace like the teeth of a zip. The second figure lays the zeros of the first five Legendre polynomials on a shared axis so the interlacing is unmistakable — and marks the payoff: the two zeros of $P_2$ are exactly the nodes of the two-point Gauss-Legendre rule.

Figure

Zeros of the Legendre polynomials P₁ through P₅ interlace on the interval −1 to 1; the two zeros of P₂ at ±0.5774 are precisely the nodes of the 2-point Gauss-Legendre rule, which integrates every polynomial of degree ≤ 3, that is 2n−1, exactly

The core equations

The whole theory is generated by one line. Fix a weight and demand that the polynomials be mutually orthogonal in the inner product it defines; every result below is a theorem about the consequences of this single demand.

\langle p_m,p_n\rangle=\int_a^b p_m(x)\,p_n(x)\,w(x)\,\mathrm{d}x=h_n\,\delta_{mn} \qquad (1)

where: $w(x)$ — weight function, $w\ge0$ on $[a,b]$ with finite moments; $h_n=\|p_n\|^2>0$ — the squared norm; $\delta_{mn}$ — Kronecker delta. Setting $m\ne n$ says "different degrees don't overlap"; setting $m=n$ merely names the norm. Whether such a family exists at all is decided by the moments of the weight, packaged into a Hankel determinant. Positivity of these determinants is exactly the statement that the weight is a genuine positive measure.

\mu_k=\int_a^b x^k w(x)\,\mathrm{d}x,\qquad \Delta_n=\det\big[\mu_{i+j}\big]_{i,j=0}^{n}>0 \qquad (2)

where: $\mu_k$ — the $k$ -th moment; $\Delta_n$ — the Hankel (moment) determinant. The Gram matrix of the monomials in this inner product is precisely $[\mu_{i+j}]$ , so $\Delta_n>0$ for all $n$ is positive-definiteness in disguise. When the Hankel determinants are positive the monic family is not merely guaranteed to exist but written in closed form as a ratio of moment determinants — an explicit, if unwieldy, definition.

p_n(x)=\frac{1}{\Delta_{n-1}}\det\!\begin{pmatrix}\mu_0&\mu_1&\cdots&\mu_n\\ \vdots&\vdots&&\vdots\\ \mu_{n-1}&\mu_n&\cdots&\mu_{2n-1}\\ 1&x&\cdots&x^n\end{pmatrix} \qquad (3)

where: each of the top $n$ rows reproduces a Gram constraint, forcing $\langle p_n,x^j\rangle=0$ for $j<n$ ; $\Delta_{n-1}\ne0$ makes the monic polynomial well-defined and unique. Reading it: the determinant is orthogonality made into an object. Everything computational, though, runs off the recurrence, not the determinant. It holds for every orthogonal family without exception.

p_{n+1}(x)=(A_n x+B_n)\,p_n(x)-C_n\,p_{n-1}(x),\qquad C_n=\frac{A_n}{A_{n-1}}\frac{h_n}{h_{n-1}}>0 \qquad (4)

where: $A_n,B_n,C_n$ — coefficients fixed by $w$ and the normalisation; $C_n>0$ always. In monic form $A_n=1$ , $B_n=-b_n$ , and it reads $p_{n+1}=(x-b_n)p_n-c_n p_{n-1}$ with $c_n=h_n/h_{n-1}>0$ . The recurrence collapses the reproducing kernel — the sum of $p_k(x)p_k(y)/h_k$ that projects onto degree- $n$ polynomials — into a single closed expression, the Christoffel-Darboux identity.

\sum_{k=0}^{n}\frac{p_k(x)\,p_k(y)}{h_k}=\frac{1}{h_n A_n}\,\frac{p_{n+1}(x)\,p_n(y)-p_n(x)\,p_{n+1}(y)}{x-y} \qquad (5)

where: the left side is the kernel $K_n(x,y)$ ; taking $y\to x$ (l'Hôpital) yields the confluent form that delivers the Gaussian weights. This identity is the engine behind completeness and convergence proofs. For the classical families only, the weight is not free: it must satisfy the Pearson equation, and then the polynomials solve a self-adjoint Sturm-Liouville eigenproblem with a purely algebraic spectrum.

(w s)'=w\,t,\qquad \frac{1}{w}\frac{\mathrm{d}}{\mathrm{d}x}\!\left(w\,s\,\frac{\mathrm{d}p_n}{\mathrm{d}x}\right)=\lambda_n\,p_n,\qquad \lambda_n=n\!\left(t'+\tfrac12(n-1)\,s''\right) \qquad (6)

where: $s(x)$ — polynomial of degree $\le2$ ; $t(x)$ — polynomial of degree $1$ ; $\lambda_n$ — eigenvalue, quadratic in $n$ . Reading it on cases: Hermite ( $s=1$ , $t=-2x$ ) gives $\lambda_n=-2n$ ; Legendre ( $s=1-x^2$ , $t=-2x$ ) gives $\lambda_n=-n(n+1)$ . The same $s$ and $w$ generate $p_n$ explicitly by repeated differentiation — the Rodrigues formula, the discrete shadow of that eigenproblem.

p_n(x)=\frac{1}{K_n\,w(x)}\frac{\mathrm{d}^n}{\mathrm{d}x^n}\!\left[w(x)\,s(x)^n\right] \qquad (7)

where: $K_n$ — normalisation constant fixing the convention; $s^n$ — the $n$ -th power of the Pearson polynomial (e.g. $(x^2-1)^n$ for Legendre, giving $P_n=\tfrac{1}{2^n n!}\tfrac{\mathrm{d}^n}{\mathrm{d}x^n}(x^2-1)^n$ ). Finally the zeros pay off in numerics: put the nodes at the zeros of $p_n$ and choose the Christoffel weights, and the quadrature is exact to a degree no other $n$ -node rule can reach.

\int_a^b f(x)\,w(x)\,\mathrm{d}x=\sum_{i=1}^{n}w_i\,f(x_i)\ \ \text{exact for }\deg f\le 2n-1,\qquad w_i=\left(\sum_{k=0}^{n-1}\frac{p_k(x_i)^2}{h_k}\right)^{-1}>0 \qquad (8)

where: $x_i$ — the $n$ zeros of $p_n$ ; $w_i$ — Christoffel numbers, manifestly positive as reciprocals of a sum of squares. Equivalently $w_i=-A_n h_n/[A_{n-1}\,p_{n+1}(x_i)\,p_n'(x_i)]$ , which is positive because $p_n'(x_i)$ and $p_{n+1}(x_i)$ carry opposite signs at a simple zero.

Derivation AWhy every orthogonal family obeys a three-term recurrence

The goal is to prove that orthogonality alone — with no assumption about the weight beyond positivity — forces the entire family into a tridiagonal recurrence, so that an orthogonal polynomial system is nothing more nor less than a symmetric tridiagonal operator acting on a graded basis. We show that when $x p_n$ is expanded in the orthogonal basis, all but three coefficients are annihilated by the symmetry of multiplication by $x$ , and we read off the recurrence and the sign of its trailing coefficient.

Step 1. Since $\deg(x p_n)=n+1$ , expand it in the complete basis of lower-degree orthogonal polynomials: $x p_n(x)=\sum_{k=0}^{n+1} c_k\,p_k(x)$ , with coefficients $c_k h_k=\langle x p_n,p_k\rangle$ obtained by taking the inner product with $p_k$ and using orthogonality. Why: an orthogonal basis makes every coefficient an explicit projection, so the problem reduces to deciding which projections vanish.

Step 2. Multiplication by $x$ is symmetric in the weighted inner product, $\langle x p_n,p_k\rangle=\langle p_n,x p_k\rangle$ , because $x$ is a real scalar function and can be moved across the integrand freely. Why: this self-adjointness is the whole mechanism — it lets us trade the degree-raising on $p_n$ for degree-raising on $p_k$ , where it can be defeated.

Step 3. For $k<n-1$ the polynomial $x p_k$ has degree $k+1<n$ , so it lies in the span of $p_0,\dots,p_{n-1}$ , all orthogonal to $p_n$ ; hence $\langle p_n,x p_k\rangle=0$ and $c_k=0$ . Only $k=n-1,n,n+1$ survive. Why: degree-counting plus orthogonality kills the tail, leaving a strictly three-diagonal band.

Step 4. Rearranging $x p_n=c_{n+1}p_{n+1}+c_n p_n+c_{n-1}p_{n-1}$ for $p_{n+1}$ gives $p_{n+1}=(A_n x+B_n)p_n-C_n p_{n-1}$ . Taking $\langle\,\cdot\,,p_{n-1}\rangle$ shows $C_n h_{n-1}=\langle x p_n,p_{n-1}\rangle=\langle p_n,x p_{n-1}\rangle=(A_{n-1}/A_n)^{-1}h_n$ up to the leading-coefficient bookkeeping, yielding $C_n=(A_n/A_{n-1})(h_n/h_{n-1})>0$ . Why: the trailing coefficient inherits the ratio of consecutive norms, so its positivity is automatic — the fact Favard's converse later exploits. *In monic form the recurrence reads $p_{n+1}=(x-b_n)p_n-c_n p_{n-1}$ with $c_n=h_n/h_{n-1}$ .*

Go deeper · Derivation B · 1 minZeros are real, simple, and interlace inside the interval

Derivation BZeros are real, simple, and interlace inside the interval

The goal is to establish the geometric backbone of the whole subject: that positivity of the weight alone confines all $n$ roots of $p_n$ to the open interval $(a,b)$ , forces them to be simple, and makes consecutive degrees interlace — the property that turns the zeros into quadrature nodes and eigenvalues. We argue by a sign-change contradiction against orthogonality, then read interlacing off the Christoffel-Darboux identity.

Step 1. Let $x_1,\dots,x_m$ be the interior points where $p_n$ actually changes sign (roots of odd multiplicity) and suppose for contradiction that $m<n$ . Build $q(x)=\prod_{j=1}^m (x-x_j)$ , a polynomial of degree $m<n$ whose sign flips exactly where $p_n$ 's does. Why: matching $q$ 's sign changes to $p_n$ 's is the device that manufactures a fixed-sign integrand.

Step 2. The product $p_n(x)\,q(x)$ never changes sign on $(a,b)$ : at each $x_j$ both factors flip together, so the product keeps one sign throughout, and it is not identically zero. Why: a fixed-sign, non-trivial function integrated against a strictly positive weight cannot integrate to zero.

Step 3. Therefore $\int_a^b p_n q\,w\,\mathrm{d}x\ne 0$ . But $\deg q=m<n$ , so orthogonality demands $\langle p_n,q\rangle=0$ — a contradiction. Hence $m=n$ : $p_n$ has $n$ interior sign changes, and being degree $n$ these are all its zeros, each simple. Why: the contradiction pins the count at exactly $n$ and rules out repeated roots in one stroke.

Step 4. Interlacing follows from Christoffel-Darboux at $y=x$ : $\sum_{k=0}^n p_k(x)^2/h_k = [p_{n+1}'(x)p_n(x)-p_n'(x)p_{n+1}(x)]/(h_n A_n)>0$ . At a zero $x_i$ of $p_n$ this forces $p_{n+1}(x_i)\,p_n'(x_i)<0$ , so $p_{n+1}$ takes opposite signs at consecutive zeros of $p_n$ and must vanish once between each pair. Why: a strictly positive kernel converts "consecutive polynomials share no zeros" into the sharper statement that their zeros strictly alternate. *The Legendre zeros illustrating the interlacing are: $P_1$ : $0$ ; $P_2$ : $\pm0.5774$ ; $P_3$ : $0,\,\pm0.7746$ ; $P_4$ : $\pm0.3400,\,\pm0.8611$ ; $P_5$ : $0,\,\pm0.5385,\,\pm0.9062$ .*

Go deeper · Derivation C · 1 minGaussian quadrature is exact to degree 2n minus 1

Derivation CGaussian quadrature is exact to degree 2n minus 1

The goal is to prove the central numerical fact of the theory: choosing the $n$ nodes to be the zeros of $p_n$ and the weights to be the Christoffel numbers integrates every polynomial of degree up to $2n-1$ exactly — twice the reach of any rule with prescribed nodes — and that $2n-1$ is the best possible, no $n$ -node rule doing better. We use polynomial division by $p_n$ and the vanishing of $p_n$ at the nodes.

Step 1. Take any polynomial $f$ with $\deg f\le 2n-1$ and divide by $p_n$ : $f=q\,p_n+r$ with $\deg q,\deg r\le n-1$ , since $\deg p_n=n$ . Why: division isolates the part of $f$ that orthogonality can kill from the part the nodes can see.

Step 2. Integrate against the weight. The quotient term vanishes, $\int_a^b q\,p_n\,w\,\mathrm{d}x=\langle p_n,q\rangle=0$ , because $\deg q\le n-1$ and $p_n$ is orthogonal to every lower-degree polynomial, leaving $\int_a^b f\,w=\int_a^b r\,w$ . Why: this is where orthogonality does the work — half of $f$ is invisible to the weighted integral.

Step 3. At each node $x_i$ , a zero of $p_n$ , we have $f(x_i)=q(x_i)p_n(x_i)+r(x_i)=r(x_i)$ ; the quadrature sum therefore reads $\sum_i w_i f(x_i)=\sum_i w_i r(x_i)$ . Why: the nodes are chosen precisely so the quotient term is invisible to the sum as well, matching what the integral already discarded.

Step 4. An interpolatory $n$ -node rule is exact for every $r$ with $\deg r\le n-1$ , so $\int_a^b r\,w=\sum_i w_i r(x_i)=\sum_i w_i f(x_i)=\int_a^b f\,w$ . Optimality: the degree- $2n$ polynomial $p_n(x)^2$ has positive weighted integral yet the sum gives zero (every node is a root), so $2n$ is unreachable. Why: exhibiting one failing degree- $2n$ integrand proves $2n-1$ is sharp, not merely sufficient. *As a concrete check: the two-point Gauss-Legendre rule reproduces $\int_{-1}^1 x^2\,\mathrm{d}x=2/3$ and $\int_{-1}^1 x^3\,\mathrm{d}x=0$ exactly, and fails at $x^4$ — giving $2/9$ instead of the correct $2/5$ .*

Limiting cases

The three classical families are distinguished entirely by the support of $w$ and the degree of the Pearson polynomial $s$ — a complete taxonomy of the "nice" orthogonal polynomials with a Rodrigues formula.

Case	Result	Meaning
Hermite (infinite line)	$w=e^{-x^2},\ x\in(-\infty,\infty),\ s=1$	Quadratic-exponential weight; eigenfunctions of the quantum harmonic oscillator. Recurrence $H_{n+1}=2xH_n-2nH_{n-1}$ , eigenvalue $\lambda_n=-2n$ .
Laguerre (half-line)	$w=x^{\alpha}e^{-x},\ x\in[0,\infty),\ s=x$	Generalised Laguerre $L_n^{(\alpha)}$ ; the radial hydrogen wavefunctions. Needs $\alpha>-1$ for integrability.
Jacobi (finite interval)	$w=(1-x)^{\alpha}(1+x)^{\beta},\ x\in[-1,1],\ s=1-x^2$	Contains Legendre ( $\alpha=\beta=0$ ), Chebyshev ( $\alpha=\beta=\mp\tfrac12$ ), Gegenbauer ( $\alpha=\beta$ ). Needs $\alpha,\beta>-1$ .
Discrete weight (mass points)	$w=\sum_i m_i\,\delta(x-x_i)$ , $N$ atoms	Only $N$ orthogonal polynomials exist ( $p_0,\dots,p_{N-1}$ ); the family truncates — Charlier, Krawtchouk, Hahn.
Degenerate weight ( $h_n\to0$ )	recurrence breaks, $C_n\to0$	A non-positive "weight" makes a Hankel determinant vanish; the monic polynomial is undefined and the recurrence degenerates.

Common trapNormalisation conventions silently rewrite every constant

The single most common error is quoting a recurrence or norm from one book while normalising as in another. The same family appears as monic ( $A_n=1$ ), $L^2$ -orthonormal ( $h_n=1$ ), or in a "standard" scaling — $P_n(1)=1$ for Legendre, leading coefficient $2^n$ for Hermite, $T_n(\cos\theta)=\cos n\theta$ for Chebyshev — and each convention changes $A_n,B_n,C_n$ and $h_n$ throughout. Worse, "the Hermite weight" is ambiguous: the physicists' $H_n$ are orthogonal in $e^{-x^2}$ , the probabilists' $He_n$ in $e^{-x^2/2}$ , and the two differ by a rescaling of $x$ that alters every coefficient. Fix the convention before writing a single formula, and never mix a recurrence from one source with a norm from another.

Worked examples

1 · Building the first Legendre polynomials from moments (confidence builder)

Construct the monic orthogonal polynomials on $[-1,1]$ with weight $w=1$ directly from the moments, and read off the zeros of $p_2$ .

Step 1. Compute the moments $\mu_k=\int_{-1}^1 x^k\,\mathrm{d}x$ : these vanish for odd $k$ by symmetry, and for even $k$ give $\mu_0=2$ , $\mu_2=2/3$ , $\mu_4=2/5$ . Positivity check: $\Delta_0=\mu_0=2>0$ and $\Delta_1=\mu_0\mu_2-\mu_1^2=4/3>0$ , so a family exists.

Step 2. Start with $p_0=1$ . Gram-Schmidt for $p_1$ : $p_1=x-\dfrac{\langle x,1\rangle}{\langle1,1\rangle}=x-\dfrac{\mu_1}{\mu_0}=\mathbf{x}$ , since $\mu_1=0$ .

Step 3. For $p_2=x^2-\alpha\,p_1-\beta\,p_0$ impose orthogonality. Against $p_0$ : $\beta=\dfrac{\langle x^2,1\rangle}{\langle1,1\rangle}=\dfrac{\mu_2}{\mu_0}=\dfrac{2/3}{2}=\dfrac13$ . Against $p_1$ : $\alpha=\dfrac{\langle x^2,x\rangle}{\langle x,x\rangle}=\dfrac{\mu_3}{\mu_2}=0$ . Hence $p_2=\mathbf{x^2-\tfrac13}$ .

Step 4. This is the monic Legendre polynomial; multiplying by $3/2$ recovers the standard $P_2=\tfrac12(3x^2-1)$ . Its zeros solve $x^2=\tfrac13$ , i.e. $x=\pm1/\sqrt3=\mathbf{\pm0.5774}$ — real, simple, interior to $(-1,1)$ , exactly as the separation theorem promises, and (as example 2 shows) exactly the two-point Gauss-Legendre nodes.

2 · Two-point Gauss-Legendre quadrature (fully worked)

Derive the two-point Gauss-Legendre rule on $[-1,1]$ , verify its exactness to degree $3$ , and use it to estimate $\int_{-1}^1 e^x\,\mathrm{d}x$ .

Step 1. The nodes are the zeros of $P_2$ , i.e. $x_{1,2}=\pm1/\sqrt3$ from example 1. The Christoffel weights follow from $w_i=2/[(1-x_i^2)\,P_2'(x_i)^2]$ with $P_2'=3x$ : at $x_i^2=\tfrac13$ this is $w_i=\dfrac{2}{(2/3)(3)}=\mathbf{1}$ , so $w_1=w_2=1$ and $\sum w_i=2=\int_{-1}^1\mathrm{d}x$ as a sanity check.

Step 2. Verify exactness. For $f=x^2$ : rule gives $1\cdot\tfrac13+1\cdot\tfrac13=\tfrac23$ , matching $\int_{-1}^1 x^2=\tfrac23$ . For $f=x^3$ (odd): rule gives $(1/\sqrt3)^3-(1/\sqrt3)^3=0=\int_{-1}^1 x^3$ . So the rule is exact through degree $2n-1=\mathbf{3}$ .

Step 3. Confirm it is not exact at degree $4$ : for $f=x^4$ the rule gives $2\cdot(1/3)^2=\tfrac29\approx0.2222$ , but $\int_{-1}^1 x^4=\tfrac25=0.4$ . The failure is expected — $x^4$ exceeds $2n-1$ .

Step 4. Apply to $\int_{-1}^1 e^x\,\mathrm{d}x$ . The rule gives $e^{0.5774}+e^{-0.5774}=1.78129+0.56139=\mathbf{2.34270}$ , against the exact $e-e^{-1}=2.35040$ . With just two function evaluations the error is $0.0077$ , a relative error of only $\mathbf{0.33\%}$ — the payoff of placing nodes at orthogonal-polynomial zeros rather than on an equal grid.

Physics of the portrait

The general theory has no single author, but it has a decisive founder. Legendre, Hermite, Laguerre, and Chebyshev each gave the world a family; it was Chebyshev who first asked what all such families have in common, and answered with the machinery — continued fractions, the moment problem, and mechanical quadrature — that underlies this entire article.

Origins · 1 minPafnuty Chebyshev and how this was discovered

Pafnuty Chebyshev

1821–1894, St Petersburg school · founder of the general theory of orthogonal polynomials and Gaussian quadrature

Pafnuty Lvovich Chebyshev built the general theory of orthogonal polynomials out of a problem in approximation. Studying how best to fit data and how to interpolate with least error, he was led in the 1850s to continued-fraction expansions of the integral $\int w(x)/(z-x)\,\mathrm{d}x$ , whose convergents are ratios of exactly the polynomials orthogonal in the weight $w$ — the first genuinely general construction, not tied to any one family. From it he extracted the three-term recurrence, the interior real zeros, and the mechanical-quadrature rule that Gauss had found for the Legendre weight, now generalised to any weight. His own polynomials $T_n$ , minimising the maximum deviation from zero, emerged as the extremal case. Chebyshev's students carried the programme forward: Andrey Markov sharpened the moment inequalities, and Thomas Stieltjes, working in the same tradition, solved the moment problem outright in 1894, tying a positive measure to its continued fraction and completing the circle of ideas — weight, moments, recurrence, quadrature — that Chebyshev had opened. The modern subject, codified in Szegő's 1939 treatise, is recognisably his.

History

1782

— Legendre introduces his polynomials in the theory of gravitational attraction

1814

— Gauss publishes mechanical quadrature using the Legendre nodes

1859

— Chebyshev founds the general theory via continued fractions and least-squares approximation

1894

— Stieltjes solves the moment problem, linking positive measures to continued fractions

1935

— Favard proves the converse: a positive-coefficient three-term recurrence defines an orthogonal family

Going further

Texts. Szegő, Orthogonal Polynomials, AMS Colloquium Publications vol. 23 (the canonical reference) · Chihara, An Introduction to Orthogonal Polynomials (Favard's theorem, moment problem) · Stone & Goldbart, Mathematics for Physics, ch. 8 (Sturm-Liouville and special functions) · Arfken & Weber, Mathematical Methods for Physicists, 7th ed., ch. 15–18 · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., ch. 18 · Abramowitz & Stegun, Handbook of Mathematical Functions, ch. 22. The Sturm-Liouville eigenvalue is stated in the operator-equals- $\lambda_n$ convention as $\lambda_n=n(t'+\tfrac12(n-1)s'')$ , which gives $\lambda_n=-2n$ for Hermite and $\lambda_n=-n(n+1)$ for Legendre — note this differs by a sign from sources that write the ODE as $sy''+ty'+\lambda y=0$ . Key anchor values: moments $\mu_0=2$ , $\mu_2=2/3$ , $\mu_4=2/5$ ; monic Legendre $p_2=x^2-\tfrac13$ with zeros $\pm1/\sqrt3$ ; two-point Gauss weights $w_1=w_2=1$ and the estimate $\int_{-1}^1 e^x\approx2.34270$ against the exact $2.35040$ (0.33% error); and $C_n>0$ in the recurrence throughout. Conventions are flagged in the pitfall box above.

Prerequisites. Gamma Function · Beta Function · Recurrence Relations for Special Functions · Generating Functions for Special Functions. Leads to. Legendre Polynomials · Associated Legendre Functions · Spherical Harmonics.

Ask the tutor: Ask me to derive the recurrence, Rodrigues formula, or Gaussian nodes for any specific family — Hermite, Laguerre, Legendre, or Chebyshev — and watch the general theory specialise to it.

Is orthogonality really the only assumption, or is a recurrence assumed too?

Orthogonality against a positive weight is the only assumption; the three-term recurrence is a theorem, not a hypothesis. Given the family is orthogonal, the symmetry of multiplication by $x$ forces $x p_n$ to expand in only $p_{n+1},p_n,p_{n-1}$ , and that is the recurrence. The striking part is Favard's converse: start instead from any recurrence $p_{n+1}=(x-b_n)p_n-c_n p_{n-1}$ with every $c_n>0$ , and a positive measure making the $p_n$ orthogonal must exist. So the two properties are logically equivalent — you can take either as the definition and derive the other.

What makes a family "classical" as opposed to merely orthogonal?

Any positive weight generates an orthogonal family, but only three (up to affine change of variable) also solve a second-order differential equation and carry a Rodrigues formula: Hermite, Laguerre, and Jacobi. The distinguishing feature is that their weight satisfies the Pearson equation $(ws)'=wt$ with $s$ of degree $\le2$ and $t$ of degree $1$ . That algebraic constraint is exactly what makes the Sturm-Liouville operator self-adjoint with polynomial eigenfunctions. A generic weight — say $w=e^{-x^4}$ — still yields a perfectly good orthogonal family with a three-term recurrence, but no differential equation of finite order and no Rodrigues formula.

Why are the zeros guaranteed to be simple and real?

Both facts come from a single contradiction argument. If $p_n$ had fewer than $n$ interior sign changes, you could build a low-degree polynomial $q$ that makes $p_n q$ one-signed on the interval; then $\int p_n q\,w$ would be nonzero, yet orthogonality forces it to vanish since $\deg q<n$ . The only escape is exactly $n$ interior sign changes, which for a degree- $n$ polynomial means $n$ distinct real roots inside $(a,b)$ — no repeated roots, none on or outside the boundary. Positivity of the weight is essential: it is what guarantees a one-signed integrand cannot integrate to zero.

Why is Gaussian quadrature exact to degree 2n−1 and not simply n−1?

An $n$ -node interpolatory rule is exact to degree $n-1$ for any node placement — that is just polynomial interpolation. Gauss doubles this by choosing the nodes cleverly: put them at the zeros of $p_n$ . Then for a polynomial $f$ of degree up to $2n-1$ , write $f=q\,p_n+r$ ; the $q\,p_n$ part integrates to zero by orthogonality, and it also contributes nothing to the sum because $p_n$ vanishes at every node. Only the degree- $\le n-1$ remainder survives, and the rule handles that exactly. You buy $n$ extra degrees of exactness by spending your freedom in node location — which fixed-node rules like Newton-Cotes throw away.

How does the Jacobi matrix connect to the recurrence and the zeros?

Write the recurrence in orthonormal form and it reads $x\,\hat p_n=a_{n+1}\hat p_{n+1}+b_n\hat p_n+a_n\hat p_{n-1}$ — multiplication by $x$ acting as a symmetric tridiagonal matrix $J$ with the $b_n$ on the diagonal and the $a_n$ off it. Truncate $J$ to its top-left $n\times n$ block and its eigenvalues are exactly the zeros of $p_n$ , while the squared first components of its eigenvectors are the Gauss weights. This is the Golub-Welsch algorithm, and it is why computing quadrature nodes is really an eigenvalue problem: the abstract statement "orthogonal polynomials are eigen-directions of a Jacobi operator" is a literal, computable identity.

Do these families form a complete basis, so I can expand any function?

For the classical weights on their natural intervals, yes: the orthonormal polynomials are complete in the weighted space $L^2(w)$ , so any square-integrable function has a mean-square-convergent expansion $f=\sum_n \langle f,\hat p_n\rangle\,\hat p_n$ , the analogue of a Fourier series. Completeness follows from the Christoffel-Darboux kernel acting as an approximate identity, together with density of polynomials — which for unbounded intervals needs the weight to decay fast enough (a Carleman-type moment condition). Pointwise and uniform convergence are more delicate and depend on the smoothness of $f$ near the endpoints, exactly as with Fourier series, but the $L^2$ completeness that physics relies on is secure.