Orthogonal Polynomials

The weighted inner product

The starting point is a weight function w(x) > 0 on an interval (a, b), with ∫ w < ∞ on enough subintervals. Two functions f and g are declared orthogonal under w if:

⟨f, g⟩_w = ∫_a^b f(x) g(x) w(x) dx = 0

A family {p_0, p_1, p_2, …} with deg p_n = n is an orthogonal polynomial family if ⟨p_n, p_m⟩_w = 0 whenever n ≠ m. Such a family always exists — the moments μ_k = ∫ xᵏ w(x) dx must all be finite, and Gram–Schmidt applied to {1, x, x², …} produces it. The result is unique up to scalar normalisation per polynomial.

Gram–Schmidt construction

Start with monomials {1, x, x², …} and subtract their projections onto previous orthogonal polynomials:

p_0(x) = 1
p_1(x) = x − ⟨x, p_0⟩ / ⟨p_0, p_0⟩ · p_0(x)
p_2(x) = x² − ⟨x², p_0⟩ / ⟨p_0, p_0⟩ · p_0 − ⟨x², p_1⟩ / ⟨p_1, p_1⟩ · p_1
...
p_{n+1}(x) = x · p_n(x) − ∑_{k=0}^{n} (⟨x p_n, p_k⟩ / ⟨p_k, p_k⟩) · p_k(x)

A remarkable simplification: due to the fact that multiplication by x shifts degree by 1, only the last two projections (k = n and k = n − 1) are nonzero. So the construction collapses to a three-term recurrence:

p_{n+1}(x) = (α_n x + β_n) p_n(x) − γ_n p_{n−1}(x)

with constants α_n, β_n, γ_n determined by ⟨x p_n, p_n⟩, ⟨p_n, p_n⟩, ⟨p_{n−1}, p_{n−1}⟩. This is the most consequential structural theorem in the subject — Favard's theorem says the converse holds, too: any sequence satisfying such a recurrence with appropriate coefficient signs is the orthogonal polynomial family of some positive weight.

The four classical families

Family	Interval	Weight w(x)	Recurrence (monic)	Where it appears
Legendre P_n	[−1, 1]	1	(n+1) P_{n+1} = (2n+1) x P_n − n P_{n−1}	Multipole expansion, sphere harmonics, plain-old integration
Chebyshev T_n	[−1, 1]	1/√(1 − x²)	T_{n+1} = 2x T_n − T_{n−1}	Minimax approximation, equiripple filters, libm polynomials
Hermite H_n	(−∞, ∞)	e^{−x²}	H_{n+1} = 2x H_n − 2n H_{n−1}	QM harmonic oscillator, Gaussian quadrature, Edgeworth series
Laguerre L_n	[0, ∞)	e^{−x}	(n+1) L_{n+1} = (2n+1 − x) L_n − n L_{n−1}	Hydrogen atom radial wavefunctions, queueing theory, light scattering
Jacobi P_n^(α,β)	[−1, 1]	(1−x)^α (1+x)^β	3-term, coefficient formulas more involved	General family containing Legendre (α = β = 0), Chebyshev (α = β = ±½), Gegenbauer
Gegenbauer C_n^(λ)	[−1, 1]	(1 − x²)^(λ − ½)	3-term	Higher-dimensional spherical harmonics, ultraspherical methods

Jacobi and Gegenbauer are generalisations containing the others as special cases. In practice "orthogonal polynomials" usually means one of the first four, and each has tabulated roots, weights, and software support (NumPy, SciPy, Boost, GSL).

Worked example — first three Hermite polynomials

Using the physicists' convention with weight e^{−x²} on (−∞, ∞):

H_0(x) = 1
H_1(x) = 2x
H_2(x) = 4x² − 2
H_3(x) = 8x³ − 12x
H_4(x) = 16x⁴ − 48x² + 12

Verify orthogonality of H_1 and H_2:

∫_{−∞}^{∞} H_1(x) H_2(x) e^{−x²} dx = ∫_{−∞}^{∞} (2x)(4x² − 2) e^{−x²} dx
                                    = ∫_{−∞}^{∞} (8x³ − 4x) e^{−x²} dx
                                    = 0  (odd integrand on symmetric interval)

The integrand is odd in x, the interval is symmetric — the integral vanishes. Orthogonality verified. (For H_1 and H_3, both odd, the integrand is even but cancels via the specific coefficient structure — the heavy lifting happens for same-parity pairs.)

Roots of orthogonal polynomials

Three universal facts about the roots of p_n in an orthogonal family:

1. All n roots are real. No complex roots can appear (this follows from the three-term recurrence; equivalently from the fact that p_n is the characteristic polynomial of a symmetric tridiagonal matrix).
2. All n roots lie strictly inside the support (a, b). No root sits on the boundary or outside it.
3. The roots of p_n and p_{n+1} strictly interlace. Between any two consecutive roots of p_n there is exactly one root of p_{n+1}.

These roots are exactly the optimal nodes for Gauss quadrature: integrating ∫ f(x) w(x) dx ≈ Σ w_i f(x_i) where x_i are roots of p_n is exact for polynomials f of degree ≤ 2n − 1.

Series expansion in an orthogonal basis

Any function f with ∫ f² w < ∞ can be expanded in the orthogonal basis:

f(x) = ∑_{n=0}^∞ c_n p_n(x)
c_n = ⟨f, p_n⟩_w / ⟨p_n, p_n⟩_w
    = (∫ f(x) p_n(x) w(x) dx) / (∫ p_n(x)² w(x) dx)

This is the orthogonal-polynomial analogue of a Fourier series — and for the right weight, often equivalent. Chebyshev expansion on [−1, 1] with substitution x = cos θ is literally a Fourier cosine series. Hermite expansion is the natural basis for Gaussian-weighted functions and underlies the Edgeworth correction in statistics.

Hermite polynomials and the quantum harmonic oscillator

The 1D quantum harmonic oscillator with frequency ω has Schrödinger equation:

−ℏ²/(2m) · ψ'' + ½ m ω² x² ψ = E ψ

After rescaling ξ = x √(mω/ℏ) and ε = 2E/(ℏω), this becomes ψ'' + (ε − ξ²) ψ = 0. The square-integrable solutions are:

ψ_n(ξ) = (mω/(πℏ))^(1/4) · (1/√(2ⁿ n!)) · H_n(ξ) · e^{−ξ²/2}
E_n = ℏω (n + ½)         for n = 0, 1, 2, ...

The Gaussian e^{−ξ²/2} envelope and the Hermite polynomial together form the eigenfunction ψ_n. The orthogonality of {H_n} with weight e^{−x²} on ℝ becomes orthonormality of {ψ_n} under the standard L² inner product on ℝ — energy eigenstates form an orthonormal basis for the Hilbert space, the cornerstone of quantum mechanics.

JavaScript — Legendre, Hermite, Chebyshev, Laguerre by recurrence

// All four classical families via three-term recurrence
function legendre(n, x) {
  if (n === 0) return 1;
  if (n === 1) return x;
  let p0 = 1, p1 = x;
  for (let k = 1; k < n; k++) {
    const pk = ((2*k + 1) * x * p1 - k * p0) / (k + 1);
    p0 = p1; p1 = pk;
  }
  return p1;
}

function hermite(n, x) {           // physicists' convention
  if (n === 0) return 1;
  if (n === 1) return 2 * x;
  let h0 = 1, h1 = 2 * x;
  for (let k = 1; k < n; k++) {
    const hk = 2 * x * h1 - 2 * k * h0;
    h0 = h1; h1 = hk;
  }
  return h1;
}

function chebyshev(n, x) {
  if (n === 0) return 1;
  if (n === 1) return x;
  let t0 = 1, t1 = x;
  for (let k = 1; k < n; k++) {
    const tk = 2 * x * t1 - t0;
    t0 = t1; t1 = tk;
  }
  return t1;
}

function laguerre(n, x) {
  if (n === 0) return 1;
  if (n === 1) return 1 - x;
  let l0 = 1, l1 = 1 - x;
  for (let k = 1; k < n; k++) {
    const lk = ((2*k + 1 - x) * l1 - k * l0) / (k + 1);
    l0 = l1; l1 = lk;
  }
  return l1;
}

console.log(legendre(3, 0.5));   // P_3(0.5) = -0.4375
console.log(hermite(4, 1));      // H_4(1) = -20
console.log(chebyshev(5, 0.5));  // T_5(0.5) = 0.5 (cos(5*60°) = cos(300°) = 0.5)
console.log(laguerre(3, 1));     // L_3(1) ≈ -0.1666...

Where orthogonal polynomials matter

• Gauss quadrature. Roots of the orthogonal polynomial for weight w are the optimal nodes for ∫ f w dx ≈ Σ w_i f(x_i). Gauss–Legendre, Gauss–Hermite, Gauss–Chebyshev, Gauss–Laguerre — each tailored to its weight.
• Spectral methods for PDEs. Solutions expanded in Chebyshev or Legendre polynomials converge spectrally (exponentially fast for smooth solutions). Used in climate models, plasma physics, computational fluid dynamics.
• Quantum mechanics. Hermite for the harmonic oscillator, Laguerre for the hydrogen atom's radial part, Legendre for the angular part (spherical harmonics use associated Legendre).
• Polynomial chaos in uncertainty quantification. Random variables expanded in Hermite (for Gaussian inputs) or Legendre (for uniform) — used in aerospace and finance to propagate uncertainty efficiently.
• Probability and statistics. Edgeworth and Gram–Charlier expansions correct normal approximations using Hermite polynomials. Generalized Laguerre is the orthogonal family for the gamma distribution.
• Random matrix theory. Distributions of eigenvalues of Gaussian random matrices are described by Hermite polynomials; of Wishart matrices by Laguerre polynomials.
• Signal processing. Hermite polynomials × Gaussian = Hermite functions = optimal time-frequency basis. Used in pulse shaping, EEG analysis, optical tweezers.

Common mistakes

• Forgetting the weight. "Orthogonal polynomials" alone doesn't pin down a family — you must specify the weight. Legendre and Chebyshev both live on [−1, 1] but are different families because their weights differ.
• Mixing physicists' and probabilists' Hermite. H_n with weight e^{−x²} (physicists') differs from He_n with weight e^{−x²/2}/√(2π) (probabilists') by a scaling. The recurrence and orthogonality constants differ. Check which convention your library uses.
• Computing polynomials from explicit coefficients. The coefficients alternate sign and grow combinatorially. Use the three-term recurrence in practice.
• Using Gauss–Legendre quadrature on a problem with a Gaussian weight. The right answer for ∫ f(x) e^{−x²} dx uses Gauss–Hermite. Each weight has its own quadrature rule.
• Assuming the orthogonal polynomials are normalized. Most textbooks use a fixed leading coefficient (e.g. T_n has leading coefficient 2^(n−1)) instead of unit L²-norm. Your library may use yet another normalisation. Always look up the convention before computing coefficients from explicit formulas.
• Confusing orthogonal polynomials with orthonormal functions. Orthogonal means ⟨p_n, p_m⟩ = 0 for n ≠ m; orthonormal additionally requires ‖p_n‖ = 1. Most classical families are not orthonormal as standardly defined — you must divide by ‖p_n‖ to make them so.