Legendre Polynomials

Definition

The Legendre polynomials P_n(x), for n = 0, 1, 2, …, are polynomials of degree n on the interval [−1, 1], normalized by P_n(1) = 1, and orthogonal to each other:

∫_{−1}^{1} P_n(x) P_m(x) dx = 0     for n ≠ m
∫_{−1}^{1} P_n(x)² dx = 2 / (2n + 1)

Three equivalent ways to define them:

1. Gram–Schmidt on 1, x, x², x³, … with the inner product ⟨f, g⟩ = ∫_{−1}^{1} fg dx and the normalization P_n(1) = 1.
2. Rodrigues' formula: P_n(x) = (1/(2ⁿ n!)) dⁿ/dxⁿ [(x² − 1)ⁿ].
3. Legendre's ODE: bounded solutions of ((1 − x²) y')' + n(n+1) y = 0 on [−1, 1].

The first few polynomials

P₀(x) = 1
P₁(x) = x
P₂(x) = (3x² − 1) / 2
P₃(x) = (5x³ − 3x) / 2
P₄(x) = (35x⁴ − 30x² + 3) / 8
P₅(x) = (63x⁵ − 70x³ + 15x) / 8

A few patterns visible at a glance:

• P_n has the same parity as n: P_{2k} is even, P_{2k+1} is odd.
• P_n(1) = 1 and P_n(−1) = (−1)ⁿ for all n.
• P_n has exactly n distinct real roots, all strictly inside (−1, 1).
• The maximum of |P_n(x)| on [−1, 1] is 1, achieved at the endpoints.

Bonnet's recurrence

Given P_0 = 1 and P_1 = x, every P_n is generated by

(n + 1) P_{n+1}(x) = (2n + 1) x P_n(x) − n P_{n−1}(x)

This three-term recurrence is numerically stable for evaluation. Try n = 1: 2 P₂ = 3 x P₁ − P₀ = 3x² − 1, so P₂ = (3x² − 1)/2. ✓ Library implementations evaluate P_n in O(n) operations using this recurrence directly.

Rodrigues' formula

A single closed-form definition:

P_n(x) = (1 / (2ⁿ n!)) · dⁿ/dxⁿ [(x² − 1)ⁿ]

Check n = 2: d²/dx² [(x² − 1)²] = d²/dx² [x⁴ − 2x² + 1] = 12x² − 4. Divide by 4·2! = 8: P₂ = (12x² − 4)/8 = (3x² − 1)/2. ✓

Rodrigues lets you derive structural properties almost for free. Differentiating (x² − 1)ⁿ and integrating by parts inside ⟨P_n, P_m⟩ gives the orthogonality relations in a few lines.

The Legendre differential equation

P_n is the unique polynomial solution of

((1 − x²) y')' + n(n + 1) y = 0

(1 − x²) y'' − 2x y' + n(n + 1) y = 0

This is the form (d/dx)(p(x) dy/dx) + λ ρ(x) y = 0 of a Sturm–Liouville problem with p(x) = 1 − x², ρ(x) = 1, and λ = n(n + 1). The endpoints x = ±1 are regular singular points; demanding bounded solutions there forces λ = n(n+1) for non-negative integer n. The Sturm–Liouville machinery then guarantees orthogonality and completeness.

Connection to spherical harmonics

Separating Laplace's equation ∇²u = 0 in spherical coordinates (r, θ, φ), assuming u = R(r) Θ(θ) Φ(φ), reduces the angular part to:

(1/sin θ) (sin θ Θ')' + [n(n+1) − m²/sin² θ] Θ = 0

For φ-independent solutions (m = 0), substitute u = cos θ to get exactly Legendre's equation with x = cos θ. The regular solutions are P_n(cos θ). The m = 0 spherical harmonic is

Y_n^0(θ, φ) = √((2n + 1) / (4π)) · P_n(cos θ)

Legendre polynomials are the rotationally-symmetric-around-z spherical harmonics. The full set of Y_n^m includes "associated Legendre functions" P_n^m(cos θ) for m ≠ 0 (a closely-related family).

Multipole expansion of 1 / |x − x′|

The Coulomb / Newton kernel has a beautiful expansion in Legendre polynomials:

1 / |x − x′| = ∑_{n=0}^∞ (r<^n / r>^{n+1}) P_n(cos γ)

where r< = min(|x|, |x′|), r> = max, and γ is the angle between x and x′. This is the generating function for P_n:

1 / √(1 − 2 x t + t²) = ∑_{n=0}^∞ P_n(x) tⁿ

Pluggable into any potential theory: the gravity or Coulomb potential of a localized source, evaluated outside the source, expands as monopole (P₀) + dipole (P₁) + quadrupole (P₂) + …, with each term carrying angular pattern P_n(cos γ) and radial decay 1/r^{n+1}.

Gauss–Legendre quadrature

An n-point quadrature rule is a formula

∫_{−1}^{1} f(x) dx ≈ ∑_{i=1}^{n} w_i f(x_i)

that takes a function value at n points and outputs an integral. Gauss–Legendre chooses the nodes x_i and weights w_i optimally: x_i are the n roots of P_n(x), and

w_i = 2 / [(1 − x_i²) (P_n'(x_i))²]

The rule is exact for any polynomial of degree ≤ 2n − 1. Twice the polynomial degree of equally-spaced Newton–Cotes rules (trapezoidal: degree 1, Simpson: degree 3). For smooth integrands, the error decays super-algebraically with n.

A 5-point Gauss–Legendre rule exactly integrates any polynomial of degree 9. To 9 decimal places:

n=5 nodes/weights:
 ±0.906180        w = 0.236927
 ±0.538469        w = 0.478629
  0               w = 0.568889

This 5-point rule beats a 100-point Simpson rule for many smooth integrands. Modern numerical-integration libraries (SciPy, NumPy, Boost, GSL) ship Gauss–Legendre as a workhorse.

Worked example — expand cos(x) in Legendre series on [−1, 1]

The expansion is f(x) = Σ_n a_n P_n(x) with a_n = ((2n + 1)/2) ∫_{−1}^{1} f(x) P_n(x) dx.

For f(x) = cos(x), the integrals can be done in closed form (giving Bessel-function expressions); numerically:

a_0 = 0.8414710     (avg value of cos x on [-1,1])
a_1 = 0              (cos is even)
a_2 = -0.4317712
a_3 = 0
a_4 = 0.0235076
a_5 = 0
a_6 = -0.000529
...

Truncating at n = 4 gives a maximum error ≈ 4 × 10⁻⁴ on [−1, 1] — better than a degree-4 Taylor series at the origin. Legendre series converge optimally in the L² sense; for smooth functions they also converge pointwise and uniformly.

Legendre vs other classical orthogonal polynomials

	Legendre P_n	Chebyshev T_n	Hermite H_n	Laguerre L_n	Jacobi P_n^{α,β}
Interval	[−1, 1]	[−1, 1]	(−∞, ∞)	[0, ∞)	[−1, 1]
Weight w(x)	1	1/√(1 − x²)	e^{−x²}	e^{−x}	(1−x)^α (1+x)^β
P₂ example	(3x² − 1)/2	2x² − 1	4x² − 2	x² − 4x + 2 (×1/2)	family
Sturm–Liouville	((1−x²)y')' + λy = 0	((1−x²)y'/√(1−x²))' + …	y'' − 2xy' + λy = 0	(xy')' + (λ − x)y = 0	generalized
Used for	Sphere, gravity	Polynomial fitting	QHO, Gauss-Hermite	Hydrogen radial	Many — interpolation
Quadrature	Gauss-Legendre	Clenshaw-Curtis	Gauss-Hermite	Gauss-Laguerre	Gauss-Jacobi

Where Legendre polynomials appear

• Multipole expansion in electrostatics and gravity. 1/|x − x′| = Σ (r<^n / r>^{n+1}) P_n(cos γ). The basis of multipole moments and the angular structure of every static potential.
• Spherical harmonics (m = 0). Y_n^0(θ, φ) = √((2n+1)/(4π)) P_n(cos θ). The rotationally-symmetric-around-z spherical harmonics.
• Gaussian quadrature. Gauss–Legendre is the default high-precision integration rule for smooth integrands on [−1, 1].
• Geodesy. Earth's gravity field expansion in spherical harmonics — the J_n coefficients of the m = 0 terms are exactly Legendre-polynomial moments of the mass distribution.
• Radiative transfer. Phase functions for light scattering in atmospheres expand as Legendre series in the cosine of the scattering angle.
• Neutron transport. The angular distribution of neutron flux in reactor physics is conventionally expanded in Legendre polynomials of the cosine of the angle.
• Polynomial regression and orthogonal regression. Legendre polynomials are an orthogonal basis on [−1, 1] — fitting them avoids the conditioning problems of plain monomials.

Common mistakes

• Confusing P_n with associated Legendre P_n^m. P_n is the "plain" Legendre polynomial (m = 0). P_n^m (for m ≠ 0) is the associated Legendre function — different family, used in spherical harmonics with m ≠ 0.
• Forgetting the 2/(2n+1) normalization. ∫ P_n² dx = 2/(2n+1), not 1. Many references and bugs come from sloppily-defined "orthonormal" Legendre.
• Using equispaced points to integrate against P_n. The whole point of Gauss–Legendre is that the optimal nodes are at the zeros of P_n. Trapezoidal/Simpson at equispaced points throws away exponential accuracy.
• Expecting convergence near singularities. If f has a jump or a singularity in [−1, 1], the Legendre series converges in L² but with Gibbs-style oscillation near the singularity. Same caveat as Fourier series.
• Mixing variable conventions. P_n(x) on [−1, 1] is the standard; P_n(cos θ) on θ ∈ [0, π] is the "physics" form. Get the substitution right or signs go wrong.