Recurrence Relations for Special Functions

Special functions almost never live alone. The Bessel functions $J_\nu(x)$ , the Legendre polynomials $P_n(x)$ , the Hermite and Laguerre polynomials, and the Gamma function itself all come in indexed families, and the members of a family are welded together by recurrence relations — algebraic identities that link neighbouring orders, degrees, or arguments with a handful of simple coefficients. These relations are not decorative. They are the most efficient route to numerical values, the cleanest tool for proving identities, and the algebraic skeleton from which the defining differential equation can be reconstructed. A physicist who knows a family's recurrences can compute it, differentiate it, integrate it, and factor its operator without ever looking up a series. Three ideas organize the whole subject. First, there is a single universal structure: orthogonality with respect to a positive weight forces every family of orthogonal polynomials into a three-term recurrence $p_{n+1}=(a_n x+b_n)p_n-c_n p_{n-1}$ , and Favard's theorem makes this an if-and-only-if statement. Bessel and Legendre obey close cousins of the same form. Second, every pure recurrence has a partner derivative relation; combine the two and you eliminate a neighbour and recover the governing second-order ODE, while reading them as raising and lowering operators factors that ODE into a ladder — the same algebra that solves the quantum oscillator. Third, and most practically, a three-term recurrence has two independent solutions, one growing and one decaying with $n$ , so it is numerically trustworthy in one direction only. This is why $J_n(x)$ must be computed downward by Miller's algorithm while $Y_n(x)$ marches upward — a lesson that, fittingly, recurs throughout computational physics.

Key takeaways

A recurrence relation ties a special function $f_n$ to its neighbours of nearby order, degree, or argument through a finite linear formula with elementary coefficients.
The three-term recurrence $p_{n+1}=(a_n x+b_n)p_n-c_n p_{n-1}$ is universal for orthogonal polynomials; Bessel obeys $J_{\nu-1}+J_{\nu+1}=\tfrac{2\nu}{x}J_\nu$ and Legendre obeys Bonnet's $(n+1)P_{n+1}=(2n+1)xP_n-nP_{n-1}$ .
Each recurrence pairs with a derivative relation; together they regenerate the ODE, and read as ladder operators they factor it — the classical shadow of $a^\dagger a$ and of $L_\pm$ .
The three-term recurrence is a symmetric tridiagonal (Jacobi) matrix: its eigenvalues are the Gauss quadrature nodes (Golub–Welsch).
A recurrence is useful only in its numerically stable direction: $J_n$ downward via Miller, $Y_n$ upward, $I_n$ downward, $K_n$ upward.
Key idea: indexed special functions obey finite linear relations among neighbouring orders; the three-term recurrence is universal for orthogonal polynomials, pairs with a derivative relation to regenerate the ODE, and must be iterated only in its numerically stable direction.

The core insight

Indexed special functions obey finite linear relations among neighbouring orders; the three-term recurrence is universal for orthogonal polynomials, pairs with a derivative relation to regenerate the governing ODE, and must be iterated only in its numerically stable direction.

Where this lives today

Recurrence relations are the silent infrastructure of scientific computing. Every numerical library that evaluates Bessel, Legendre, Coulomb, or spheroidal functions — from GNU GSL to the routines inside MATLAB, SciPy, and Boost — does so through recurrences, with the upward-versus-downward decision made by the stability analysis above; a wrong direction is a classic silent bug that returns plausible garbage. In cosmology the analysis of Cosmic Microwave Background maps demands associated Legendre functions $P_\ell^m$ up to $\ell\sim3000$ , computed by carefully stabilized recurrences in both $\ell$ and $m$ , and the same relations reappear as the classical shadow of the quantum angular-momentum ladder $L_\pm$ acting on $|\ell,m\rangle$ . Gaussian quadrature — the default integration rule in finite-element solvers, lattice QCD, and quantum-chemistry codes — is generated by the Golub–Welsch eigenvalue algorithm straight from recurrence coefficients. And in symbolic computation, Zeilberger's algorithm turns the search for hypergeometric-sum identities into the discovery of a recurrence and its certificate, letting a computer prove closed forms that once took a career. The humble three-term relation is, quietly, one of the most-executed pieces of mathematics in physics.

The physical picture

Think of a recurrence as a discrete dynamical system in the order $n$ : hand it two starting members of the family and it marches out the entire sequence, exactly as a two-step integrator marches an ODE forward in time. The two independent solutions of the recurrence are the discrete analogues of the two solutions of the underlying second-order equation — one that grows with $n$ and one that decays. The first figure shows this marching concretely for the Legendre polynomials, where Bonnet's recursion turns the seeds $P_0=1$ and $P_1=x$ into the whole tower of polynomials with nothing but multiplication by $x$ and a subtraction.

Figure

Bonnet's recursion marches along the order n: two neighbours Pₙ₋₁ and Pₙ generate the next Legendre polynomial through (n+1)Pₙ₊₁ = (2n+1)x Pₙ − n Pₙ₋₁, so the seeds P₀ = 1 and P₁ = x build P₂ = ½(3x²−1), then P₃ = ½(5x³−3x), then P₄ = ⅛(35x⁴−30x²+3)

The coefficient structure is not arbitrary; it encodes the geometry of the weight function. For orthogonal polynomials the recurrence coefficients are moments of the weight, and the whole recurrence assembles into a symmetric tridiagonal Jacobi matrix whose eigenvalues are the Gauss quadrature nodes — a three-term recurrence literally is a tridiagonal operator. Reading the same relations differently, a raising step $f_n\to f_{n+1}$ and a lowering step $f_n\to f_{n-1}$ multiply to a second-order operator that annihilates $f_n$ , the classical mirror of $a^\dagger a$ for the oscillator. But the sharpest practical fact is stability, and the second figure is its punchline: because $J_n(x)$ is the decaying solution of its recurrence, iterating toward it upward is catastrophic, while iterating away from it downward is exact.

Figure

Numerical stability of the Bessel order recurrence at x = 1: Jₙ is the minimal solution and falls from J₀ = 0.7652 to J₉ ≈ 5.2×10⁻⁹, but upward iteration seeded with a 10⁻¹⁰ error peels off near n = 7 and is wrong by five orders of magnitude by n = 9, whereas Miller downward recursion recovers Jₙ to eight digits

The core equations

The single most important formula in the theory is the universal three-term recurrence for orthogonal polynomials. It says that if a family $\{p_n\}$ is orthogonal with respect to any positive weight, then $p_{n+1}$ is a linear-in- $x$ multiple of $p_n$ minus a constant multiple of $p_{n-1}$ — no other members ever appear. Favard's theorem supplies the converse: any such relation with $c_n>0$ generates an orthogonal family.

p_{n+1}(x)=(a_n x+b_n)\,p_n(x)-c_n\,p_{n-1}(x) \qquad (1)

where: $p_n$ — polynomial of degree $n$ , orthogonal w.r.t. weight $w(x)$ ; $x$ — real argument; $a_n,b_n,c_n$ — coefficients fixed by the weight and normalization, with $a_n,c_n>0$ . Reading: given two consecutive members and the coefficients, every higher member follows by arithmetic; the sign structure $-c_n$ is what makes the family oscillate and its roots interlace. Bessel functions of the first kind satisfy the direct analogue, with the order $\nu$ playing the role of the degree and the argument entering only through the $1/x$ coefficient. It is the workhorse identity of cylinder-function computation.

J_{\nu-1}(x)+J_{\nu+1}(x)=\frac{2\nu}{x}\,J_\nu(x) \qquad (2)

where: $J_\nu$ — Bessel function of the first kind, order $\nu$ ; $x$ — argument; $\nu$ — order (integer or real). Reading: three consecutive orders at a fixed $x$ are linearly dependent; the same relation holds for $Y_\nu$ , so any cylinder function is fixed by two starting orders. Every pure recurrence has a derivative partner. For Bessel the two forms below are the raising and lowering relations in disguise; adding and subtracting them exchanges the pure recurrence (2) for the derivative and back.

J_\nu'(x)=\tfrac{1}{2}\big(J_{\nu-1}(x)-J_{\nu+1}(x)\big),\qquad J_\nu'(x)=J_{\nu-1}(x)-\tfrac{\nu}{x}J_\nu(x) \qquad (3)

where: $J_\nu'$ — derivative of $J_\nu$ with respect to $x$ . Reading: the derivative of one order is a combination of neighbouring orders — the "differential-difference" property that lets a recurrence carry the full content of the ODE. For Legendre polynomials the three-term recurrence has a name, Bonnet's recursion, and explicit integer coefficients, which is what makes it the standard way to tabulate $P_n$ and to build spherical harmonics.

(n+1)P_{n+1}(x)=(2n+1)\,x\,P_n(x)-n\,P_{n-1}(x) \qquad (4)

where: $P_n$ — Legendre polynomial of degree $n$ , with $P_0=1,\ P_1=x$ ; $x$ — argument, $-1\le x\le 1$ in physics; $n$ — degree. Reading: this is equation (1) with $a_n=\tfrac{2n+1}{n+1}$ , $b_n=0$ , $c_n=\tfrac{n}{n+1}$ , the coefficients dictated by the weight $w=1$ on $[-1,1]$ . The Hermite polynomials that solve the quantum oscillator carry both a pure recurrence and a clean lowering relation, and the two together are the ladder operators of the harmonic oscillator written for polynomials.

H_{n+1}(x)=2x\,H_n(x)-2n\,H_{n-1}(x),\qquad H_n'(x)=2n\,H_{n-1}(x) \qquad (5)

where: $H_n$ — physicists' Hermite polynomial, $H_0=1,\ H_1=2x$ . Reading: $H_n'=2nH_{n-1}$ is exactly the lowering operator $a$ ; multiplying by the raising step reproduces Hermite's equation and the oscillator ladder. The Gamma function shows that a recurrence need not step a discrete index at all — its functional equation is a one-step recurrence in a continuous parameter, and it is the engine of $\Gamma$ 's analytic continuation.

\Gamma(z+1)=z\,\Gamma(z),\qquad \Gamma(n+1)=n! \qquad (6)

where: $\Gamma$ — Euler Gamma function; $z$ — complex argument, $z\neq 0,-1,-2,\dots$ . Reading: iterating downward pushes any $z$ into a tabulated unit strip; iterating the relation across $z=0$ is what defines $\Gamma$ at negative non-integers and exposes its poles. The same pattern threads through the rest of the special-function zoo. The modified Bessel functions and the associated Legendre functions each obey their own three-term relations, which drive their libraries.

I_{\nu-1}(x)-I_{\nu+1}(x)=\frac{2\nu}{x}I_\nu(x),\qquad (2\ell+1)\,x\,P_\ell^m(x)=(\ell-m+1)P_{\ell+1}^m(x)+(\ell+m)P_{\ell-1}^m(x) \qquad (7)

where: $I_\nu$ — modified Bessel function of the first kind; $P_\ell^m$ — associated Legendre function of degree $\ell$ , order $m$ . Reading: the sign flip in the $I_\nu$ relation (a difference, not a sum) is precisely what reverses which solution is minimal, so $I_n$ is computed downward while $K_n$ is computed upward.

Derivation AOrthogonality forces a three-term recurrence

The claim is that orthogonality alone — not any special feature of a particular weight — forces every family of orthogonal polynomials into a relation with exactly three terms. The engine is one structural fact: multiplication by $x$ is self-adjoint under the weighted inner product, so an overlap of $p_n$ with $x\,p_k$ can be shifted onto $p_n$ itself, and a degree count then annihilates almost every term in the expansion of $x\,p_n$ . We expand, apply self-adjointness, and count degrees, and the three surviving terms are the recurrence.

Step 1. Expand $x\,p_n$ in the orthogonal basis, $x\,p_n(x)=\sum_{k=0}^{n+1}\alpha_k\,p_k(x)$ . Why: $x\,p_n$ has degree $n+1$ , and $\{p_0,\dots,p_{n+1}\}$ is a basis for all polynomials of degree $\le n+1$ , so the sum is finite and exact.

Step 2. Project onto $p_k$ using the inner product $\langle f,g\rangle=\int f\,g\,w\,\mathrm{d}x$ , giving $\alpha_k=\dfrac{\langle x\,p_n,p_k\rangle}{\langle p_k,p_k\rangle}=\dfrac{\langle p_n,\,x\,p_k\rangle}{\langle p_k,p_k\rangle}$ . Why: multiplication by $x$ is self-adjoint, $\langle x f,g\rangle=\langle f,x g\rangle$ , so the $x$ can be moved off $p_n$ and onto $p_k$ .

Step 3. Note that $x\,p_k$ has degree $k+1$ ; whenever $k+1<n$ , i.e. $k<n-1$ , it lies in the span of $\{p_0,\dots,p_{n-1}\}$ and is orthogonal to $p_n$ , so $\alpha_k=0$ . Why: only $k=n-1,\,n,\,n+1$ can survive, which is exactly three terms — no matter what the weight is.

Step 4. Rename the three survivors and solve for the top member to get $p_{n+1}=(a_n x+b_n)p_n-c_n p_{n-1}$ , with $a_n,b_n,c_n$ built from leading coefficients and norms; Favard's theorem is the converse, that any such relation with $c_n>0$ defines an orthogonal family. Why: this closes the loop — the recurrence and orthogonality are logically equivalent, not merely compatible.

Go deeper · Derivation B · 1 minFrom ladder relations back to Bessel's equation

Derivation BFrom ladder relations back to Bessel's equation

The point of this derivation is that a recurrence is not a weaker fact than the differential equation — it is an equivalent one. Starting from the two first-order ladder relations for $J_\nu$ , we differentiate one of them and use the other to eliminate every neighbouring order, and what remains is Bessel's second-order equation. This is the general mechanism by which a recurrence plus a derivative relation regenerates the governing ODE for any of the classical families.

Step 1. Take the raising relation $x\,J_\nu'(x)=\nu J_\nu(x)-x\,J_{\nu+1}(x)$ , equivalent to $\tfrac{\mathrm{d}}{\mathrm{d}x}(x^{-\nu}J_\nu)=-x^{-\nu}J_{\nu+1}$ evaluated at order $\nu$ . Why: this expresses the derivative of $J_\nu$ using the single higher neighbour $J_{\nu+1}$ , giving us something to differentiate.

Step 2. Record the lowering relation $x\,J_\nu'(x)=-\nu J_\nu(x)+x\,J_{\nu-1}(x)$ , from $\tfrac{\mathrm{d}}{\mathrm{d}x}(x^{\nu}J_\nu)=x^{\nu}J_{\nu-1}$ . Why: adding the two ladder relations recovers the derivative form (3), and subtracting them recovers pure recurrence (2) — we now hold both.

Step 3. Differentiate the raising relation with the product rule, $x\,J_\nu''+J_\nu'=\nu J_\nu'-J_{\nu+1}-x\,J_{\nu+1}'$ . Why: introducing $J_\nu''$ is the only new ingredient a second-order equation can need; everything else must cancel.

Step 4. Substitute $J_{\nu+1}=\tfrac{\nu}{x}J_\nu-J_\nu'$ from the raising relation, and $J_{\nu+1}'$ from the lowering relation one order up, so every neighbour is written in $J_\nu$ and $J_\nu'$ alone. Why: eliminating $J_{\nu\pm1}$ is what collapses the difference structure into a pure differential one.

Step 5. Collect terms; the neighbours cancel identically and leave $x^2 J_\nu''+x J_\nu'+(x^2-\nu^2)J_\nu=0$ , Bessel's equation. Why: the recurrence carried the full dynamical content all along — the ODE was hiding inside two first-order relations.

Go deeper · Derivation C · 1 minThe three-term recurrence as a Jacobi matrix

Derivation CThe three-term recurrence as a Jacobi matrix

Here the payoff is structural and computational at once. Written for the orthonormal family, the three-term recurrence becomes symmetric, and stacking the polynomials into a vector turns it into the action of a symmetric tridiagonal matrix — the Jacobi matrix. Its eigenvalues are exactly the roots of the polynomials, hence the nodes of Gaussian quadrature, and its eigenvectors give the weights. This is the Golub–Welsch algorithm, and it is why one never root-finds a high-degree polynomial directly.

Step 1. For the orthonormal polynomials $\tilde p_n$ the recurrence is symmetric, $x\,\tilde p_n=\beta_n\tilde p_{n+1}+\alpha_n\tilde p_n+\beta_{n-1}\tilde p_{n-1}$ . Why: rescaling to unit norm makes the sub- and super-diagonal coefficients equal, which is what lets a real symmetric matrix appear.

Step 2. Stack the first $N$ members into $\mathbf p(x)=(\tilde p_0,\dots,\tilde p_{N-1})^{T}$ , so the recurrence reads $x\,\mathbf p(x)=J_N\,\mathbf p(x)+\beta_{N-1}\tilde p_N(x)\,\mathbf e_N$ , where $J_N$ is the tridiagonal Jacobi matrix with diagonal $\alpha_n$ and off-diagonal $\beta_n$ . Why: the recurrence for every index except the last is captured exactly by $J_N$ ; the leftover term is the single boundary defect.

Step 3. At a root $x_j$ of $\tilde p_N$ the defect vanishes, $\tilde p_N(x_j)=0$ , so $J_N\,\mathbf p(x_j)=x_j\,\mathbf p(x_j)$ — the quadrature nodes are precisely the eigenvalues of $J_N$ . Why: this converts high-degree root-finding into a stable symmetric eigenproblem.

Step 4. The weight for node $x_j$ is $w_j=\mu_0\,(v_{j,0})^2$ , the squared first component of the normalized eigenvector times $\mu_0=\int w\,\mathrm{d}x$ . Why: this completes Golub–Welsch — recurrence coefficients in, Gauss nodes and weights out, with no explicit polynomials evaluated.

Limiting cases

Which solution of a recurrence dominates, and therefore which direction is safe to iterate, depends entirely on the regime of order versus argument. The table reads off the asymptotics that decide the computation direction.

Case	Result	Meaning
Bessel, $n\gg x$	$J_n(x)\sim\dfrac{1}{\sqrt{2\pi n}}\Big(\dfrac{ex}{2n}\Big)^{n}$	$J_n$ is the minimal (rapidly decaying) solution; upward recursion is unstable, so use downward Miller recursion.
Bessel, $Y_n$	$Y_n(x)\sim-\dfrac{(n-1)!}{\pi}\Big(\dfrac{2}{x}\Big)^{n}$	$Y_n$ is the dominant (growing) solution; upward recursion is stable and standard.
Legendre, large $n$	$\vert P_n(\cos\theta)\vert\lesssim\sqrt{\dfrac{2}{\pi n\sin\theta}}$	Both solutions of Bonnet's recurrence stay bounded on $(-1,1)$ , so upward recursion is stable there.
Modified Bessel	$I_n$ decays, $K_n$ grows with $n$	The sign flip in (7) swaps roles: $I_n$ downward (minimal), $K_n$ upward (dominant).
Gamma recurrence	$\Gamma(z+1)=z\,\Gamma(z)$	Iterated to shift $z$ into a tabulated strip; continues $\Gamma$ analytically to all $z\neq0,-1,-2,\dots$

Common trapA recurrence is not safe to iterate in either direction

The most expensive mistake in special-function computation is assuming a recurrence can be run whichever way is convenient. A three-term recurrence supports two independent solutions, one that grows with $n$ (the dominant solution) and one that decays (the minimal solution). Iterating toward the minimal solution is catastrophically unstable: rounding error seeds a trace of the dominant solution, and it grows geometrically until it swamps the answer within a handful of steps. For $J_n(x)$ with $n\gg x$ the minimal solution is the one you want, so upward iteration destroys it and you must recur downward (Miller's algorithm), fixing the overall scale at the end with a normalization sum. The rule generalizes: always iterate away from the minimal solution, never toward it.

Worked examples

1 · Computing Bessel values by Miller's backward recursion (fully worked)

Compute $J_0(1),J_1(1),\dots,J_5(1)$ using only the recurrence $J_{n-1}=\tfrac{2n}{x}J_n-J_{n+1}$ at $x=1$ , starting from the top with arbitrary seeds — and show why the naive upward march fails.

Step 1. Seed at a high order well above $x=1$ : take $N=8$ and set the trial values $f_9=0,\ f_8=1$ (any small seed works; the scale is fixed later). At $x=1$ the recurrence is $f_{n-1}=2n\,f_n-f_{n+1}$ , pure integer arithmetic.

Step 2. Iterate downward, $n=8,7,\dots,1$ : the chain comes out $f_8=1,\ f_7=16,\ f_6=223,\ f_5=2660,\ f_4=26377,\ f_3=208356,\ f_2=1223759,\ f_1=4686680,\ f_0=8149601$ . Each step damps the unwanted dominant component relative to the wanted $J_n$ .

Step 3. Fix the scale with the closure identity $1=J_0(x)+2\sum_{k\ge1}J_{2k}(x)$ , i.e. divide by $S=f_0+2(f_2+f_4+f_6+f_8)$ . Here $S=8149601+2(1223759+26377+223+1)=\mathbf{10\,650\,321}$ .

Step 4. Normalize: $J_n(1)=f_n/S$ . This gives $J_0=8149601/10650321=\mathbf{0.7651977}$ , $J_1=\mathbf{0.4400506}$ , $J_2=\mathbf{0.1149035}$ , $J_3=\mathbf{0.0195634}$ , $J_4=\mathbf{0.0024766}$ , $J_5=\mathbf{0.00024976}$ — every one correct to the digits shown, from nothing but the recurrence and a normalization.

Step 5. Contrast the upward march. Start from the exact $J_0(1)=0.7651977$ and $J_1(1)=0.4400506$ but with a $10^{-10}$ error, and run $J_{n+1}=2n\,J_n-J_{n-1}$ upward: the values track truth to $n=6$ , then peel away — by $n=9$ the recurrence returns $\sim8\times10^{-4}$ against the true $5.2\times10^{-9}$ , wrong by five orders of magnitude. Same relation, opposite direction, opposite fate. Downward is exact to 8 digits; upward blows up.

2 · Building Legendre polynomials with Bonnet's recursion

Use Bonnet's recursion (4) with the seeds $P_0=1,\ P_1=x$ to generate $P_2,P_3,P_4$ , and evaluate $P_3$ at $x=0.5$ .

Step 1. $n=1$ : $(1+1)P_2=(2\cdot1+1)xP_1-1\cdot P_0=3x\cdot x-1$ , so $P_2=\tfrac12(3x^2-1)$ — the standard quadrupole polynomial.

Step 2. $n=2$ : $3P_3=5x\,P_2-2P_1=5x\cdot\tfrac12(3x^2-1)-2x=\tfrac12(15x^3-5x)-2x=\tfrac12(15x^3-9x)$ , hence $P_3=\tfrac12(5x^3-3x)$ .

Step 3. $n=3$ : $4P_4=7x\,P_3-3P_2=7x\cdot\tfrac12(5x^3-3x)-3\cdot\tfrac12(3x^2-1)=\tfrac12(35x^4-21x^2-9x^2+3)$ , giving $P_4=\tfrac18(35x^4-30x^2+3)$ .

Step 4. Numerical check at $x=0.5$ : $P_3(0.5)=\tfrac12(5\cdot0.125-3\cdot0.5)=\tfrac12(0.625-1.5)=\mathbf{-0.4375}$ , and re-running Bonnet numerically ( $P_2(0.5)=-0.125$ , then $3P_3=5(0.5)(-0.125)-2(0.5)=-1.3125$ ) gives the same $P_3(0.5)=-0.4375$ . Recurrence and closed form agree exactly.

Physics of the portrait

Recurrence relations sit at the crossroads of three centuries of analysis: an empirical shortcut for computing tables, a theorem tying orthogonality to a tridiagonal operator, and the discrete face of the ladder algebra that quantizes the oscillator and angular momentum. The portrait here is the structure itself.

Origins · 1 minthree-term recurrence and how this was discovered

three-term recurrence

the identity pₙ₊₁ = (aₙx + bₙ)pₙ − cₙpₙ₋₁ linking neighbouring orders · orthogonal polynomials and special-function computation

The three-term recurrence is one of those objects that was used for a century before it was understood. Legendre and Bessel wrote down the relations for their own functions in the early 1800s as computational conveniences, with no sense that a single theorem governed them all. Pafnuty Chebyshev, working on continued fractions and least-squares approximation in the 1850s, was the first to treat orthogonal polynomials as a general class and to see the recurrence as their defining feature. Thomas Stieltjes made the connection exact in the 1890s: he showed that a three-term recurrence is the same data as a continued fraction and as a moment problem, so the coefficients $a_n,b_n,c_n$ encode the weight completely. The converse — that any recurrence with $c_n>0$ produces an orthogonal family — was proved by Jean Favard in 1935 (and independently by Shohat and Natanson), closing the equivalence and giving the theorem its name. The computational side matured later: J.C.P. Miller's 1952 Bessel tables introduced backward recursion as the cure for instability, and Golub and Welsch in 1969 recognized the recurrence as a symmetric tridiagonal matrix whose eigenvalues are the Gauss quadrature nodes. What Legendre used as arithmetic is now known to be a self-adjoint operator, a continued fraction, and a ladder algebra at once.

History

1782

— Legendre introduces the polynomials Pₙ(x) in the theory of the gravitational potential

1824

— Bessel systematizes Jν(x) and their order recurrence in planetary perturbation theory

1935

— Favard proves the converse, that the recurrence characterizes orthogonal families

1952

— Miller's backward-recurrence algorithm stabilizes the numerical computation of Jₙ

1969

— Golub and Welsch read Gauss quadrature nodes off the Jacobi matrix

Going further

References. Stone & Goldbart, Mathematics for Physics (CUP, 2009), chapters on special functions and orthogonal polynomials · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., chs. 14–15 (Bessel, Legendre) · Hassani, Mathematical Physics, 2nd ed., chs. on Sturm–Liouville systems and classical orthogonal polynomials · Byron & Fuller, Mathematics of Classical and Quantum Physics, ch. 6 · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., ch. 18 · Abramowitz & Stegun, Handbook of Mathematical Functions, secs. 8.5, 9.1, 9.12, 22.7 · NIST DLMF, secs. 8.5, 10.6, 10.74, 14.10, 18.9. The Miller worked example was computed directly: seeds $f_9=0,\ f_8=1$ at $x=1$ give the integer chain $f_0=8\,149\,601,\dots$ , normalization $S=10\,650\,321$ , and $J_0(1)=0.7651977$ (true value $0.76519769$ ), with all six orders matching the known values to the digits shown; the upward march seeded with a $10^{-10}$ error reaches $\sim8\times10^{-4}$ at $n=9$ against the true $5.2\times10^{-9}$ . Bonnet gives $P_2=\tfrac12(3x^2-1)$ , $P_3=\tfrac12(5x^3-3x)$ with $P_3(0.5)=-0.4375$ , and $P_4=\tfrac18(35x^4-30x^2+3)$ ; Hermite gives $H_2=4x^2-2$ ; and $\Gamma(z+1)=z\Gamma(z)$ yields $\Gamma(n+1)=n!$ .

Builds on: Gamma Function · Orthogonal Polynomials: General Theory · Generating Functions for Special Functions. Leads to: Bessel Functions of the First and Second Kind · Legendre Polynomials · Associated Legendre Functions.

Ask the tutor: Ask me to run Miller's algorithm on a Bessel value of your choice, or to build any special-function family up from its recurrence step by step.

What is the difference between a recurrence relation and a derivative relation?

A pure recurrence links a function to its neighbours of different order at the same argument — for Bessel, $J_{\nu-1}+J_{\nu+1}=\tfrac{2\nu}{x}J_\nu$ connects three orders with no derivatives. A derivative (or differential–difference) relation instead expresses $f_n'$ in terms of neighbouring orders, such as $J_\nu'=\tfrac12(J_{\nu-1}-J_{\nu+1})$ . They are partners: adding and subtracting the two ladder forms converts one into the other, and combining a pure recurrence with its derivative relation eliminates a neighbour and regenerates the governing second-order ODE. Together they carry exactly the information of the differential equation, which is why you can compute, differentiate, and factor a special function knowing only its recurrences.

Why is Favard's theorem important if we already have the recurrence?

The forward direction — orthogonality implies a three-term recurrence — is what lets you compute a known family. Favard's theorem runs the implication backward: any three-term recurrence of the form $p_{n+1}=(a_nx+b_n)p_n-c_np_{n-1}$ with $c_n>0$ automatically generates a family that is orthogonal with respect to some positive measure. That converse is what makes the recurrence a genuine characterization rather than a mere consequence. It underlies the moment problem, the theory of continued fractions, and the construction of custom orthogonal families for a chosen weight — you can invent the coefficients and be guaranteed an orthogonality behind them, without ever writing down the weight explicitly.

How do I know which direction to iterate a recurrence?

Look at the asymptotics of the two independent solutions as the index grows. One solution dominates (grows) and one is minimal (decays); you must iterate away from the minimal solution, never toward it. Concretely: compute or recall how each solution behaves for large $n$ . For Bessel with $n\gg x$ , $J_n$ decays super-exponentially while $Y_n$ grows, so $J_n$ is minimal and needs downward recursion, whereas $Y_n$ is dominant and is safely built upward. The modified functions flip: $I_n$ is minimal (downward), $K_n$ dominant (upward). If both solutions stay comparable — as for Legendre on $(-1,1)$ — either direction is stable. When in doubt, test a case with known values and watch for blow-up.

What does Miller's algorithm actually do, and why the strange normalization?

Miller's algorithm computes the minimal solution (like $J_n$ ) by exploiting the very instability that dooms upward recursion. You start at an order $N$ comfortably above both $x$ and the highest order you need, seed arbitrary values $f_{N+1}=0,\ f_N=1$ , and iterate downward. Because downward iteration drives toward the minimal solution, each step suppresses the unwanted dominant component, and after enough steps the ratios $f_n/f_m$ converge to the true ratios $J_n/J_m$ regardless of the arbitrary seed. What you have lost is the overall scale — the seed was arbitrary — so you recover it at the end with a known closure identity such as $1=J_0(x)+2\sum_{k\ge1}J_{2k}(x)$ , which fixes the single missing constant. It is one of the most elegant tricks in numerical analysis: turning instability into precision.

Are recurrence relations connected to quantum ladder operators?

Directly. The raising and lowering relations of a special function are the classical, polynomial version of the creation and annihilation operators of quantum mechanics. For Hermite polynomials the lowering relation $H_n'=2nH_{n-1}$ is literally the action of the oscillator's annihilation operator $a$ on the number states, and the pure recurrence $H_{n+1}=2xH_n-2nH_{n-1}$ encodes $x=\tfrac{1}{\sqrt2}(a+a^\dagger)$ . The same holds for angular momentum: the recurrences that step the order $m$ of the associated Legendre functions are the shadow of $L_\pm$ moving between $|\ell,m\rangle$ states. So the two first-order relations that factor a special-function ODE are the same factorization that Dirac used to solve the oscillator algebraically — recurrences and ladders are one idea in two languages.

Do recurrence relations only apply to orthogonal polynomials?

No — they are far more general. Any indexed family defined by a linear structure tends to satisfy one. The Gamma function obeys $\Gamma(z+1)=z\Gamma(z)$ , a recurrence in a continuous parameter with no polynomial in sight. Bessel and Hankel functions, Coulomb wavefunctions, spherical harmonics, hypergeometric functions (through Gauss's contiguous relations), and Clebsch–Gordan coefficients all satisfy recurrences. What orthogonal polynomials add is the guarantee of a three-term relation and its equivalence to a self-adjoint tridiagonal operator. Outside that class the number of terms can be larger and the coefficients messier, but the underlying idea — that neighbouring members of a family are linearly dependent in a computable way — is one of the most pervasive structures in mathematical physics.