Sturm-Liouville Theory and Eigenfunction Expansions

The theory is the silent backbone of computational physics. Every finite-element and spectral-method PDE solver diagonalizes a discretized self-adjoint operator and expands the solution in its eigenmodes — the descendant of $f = \sum c_n y_n$. In quantum mechanics the Schrödinger equation is a singular S–L problem, and bound-state expansions, orthonormal eigenstates and variational energy bounds (the Rayleigh–Ritz method behind quantum chemistry and density-functional theory) are S–L theory in physics notation; an electronic-structure code computing a molecule's ground state is Example 2 at scale. Structural engineering reads the eigenvalues as vibration frequencies and buckling loads — the modal analysis certifying a bridge or turbine blade; in photonics and acoustics they are the guided modes of a fiber or waveguide, with the S–L weight as the mode-overlap metric. Even the open frontier is S–L-shaped: "can you hear the shape of a drum?", $\mathcal{PT}$-symmetric operators with real spectra, and quantum graphs are all live programs about the same operator, its weight, and its boundary conditions.

Hold the whole theory in one image: the Sturm–Liouville operator is an infinite-dimensional Hermitian matrix. A real symmetric $M$ has real eigenvalues and orthogonal eigenvectors spanning the space; the self-adjoint $L$ has real eigenvalues $\lambda_n$ and orthogonal eigenfunctions $y_n$ spanning $L^2_w$. Separating variables is just diagonalizing $L$ and re-expressing the data in its eigenbasis — rotating a vector into the principal axes of a quadratic form. The first figure shows that eigenbasis for the cleanest case, a string clamped at both ends, and carries the quantization rule Sturm's theorem makes visual: mode $n$ has exactly $n-1$ interior zeros, so counting nodes is ordering the spectrum.

The second image is the payoff of completeness. Because $\{y_n\}$ spans $L^2_w$, any square-integrable datum can be written $f = \sum_n c_n y_n$ with $c_n$ the projection onto mode $n$. The picture rebuilds a blunt target, $f = 1$, from a few sine modes: one term is a lone arch, five odd terms track the flat top, and the overshoot near each end is the Gibbs phenomenon — the signature that $L^2$ (mean-square) convergence is not convergence at every point. Completeness promises the mean-square error vanishes, not a pointwise match at a jump.

The starting point is the equation itself, in divergence (self-adjoint) form where the top term appears as the derivative of $p\,y'$ — the packaging that later lets an integration by parts land cleanly.

where: $p(x)$ — leading coefficient, $p \in C^1$ with $p > 0$ on $(a,b)$ for the regular case; $q(x)$ — real continuous potential-like term; $w(x)$ — weight (density), $w > 0$, defining the inner product; $\lambda$ — eigenvalue (separation constant); $y(x)$ — eigenfunction. Reading it: a stiffness $p$, a restoring term $q$, and an eigenvalue-scaled mass density $w$ set the balance a mode must strike. The same equation rearranges into an operator eigenvalue problem $L[y] = \lambda y$, the form that makes the analogy with a Hermitian matrix exact.

where: $L$ — the self-adjoint Sturm–Liouville operator on the weighted space $L^2_w[a,b]$. Reading it: its eigenvalues are the physical spectrum (frequencies, decay rates, energies) and its eigenfunctions are the modes. The theory is universal because every second-order linear ODE can be dragged into form (1) by an integrating factor — which is why Bessel, Legendre, Hermite and Laguerre all inherit the S–L guarantees.

where: $a_2, a_1, a_0$ — the original coefficients, $a_2 \neq 0$. Reading it: the weight $w = p/a_2$ is forced by the ODE, and computing coefficients without it is the most common source of wrong answers. Orthogonality only makes sense once the inner product is fixed, and for S–L it is the weighted one — $w$ is the metric of the function space ($w = 1$ for Fourier, $w = r$ for Bessel, $w = e^{-x^2}$ for Hermite).

where: $\overline{f}$ — complex conjugate (eigenfunctions may be chosen real in the regular case). Reading it: "orthogonal" always means $\langle y_m, y_n\rangle_w = 0$, never the unweighted integral. The heart of the machinery is Lagrange's identity, which converts the difference of two operator actions into a pure boundary term. Self-adjointness — and every spectral guarantee — is exactly the statement that this boundary term vanishes.

where: $u, v$ — twice-differentiable functions on $[a,b]$. Reading it: Dirichlet, Neumann, Robin and periodic conditions all zero the right-hand side, making $\langle u, Lv\rangle_w = \langle Lu, v\rangle_w$. Once $L$ is self-adjoint, expansions become projection: any datum in $L^2_w$ is the sum of its eigenmode components, each coefficient a single inner product.

where: $c_n$ — the generalized Fourier coefficient (projection onto mode $n$); $f$ — any function in $L^2_w[a,b]$. Reading it: a coefficient is the data dotted against a normalized mode. Because $L$ is bounded below, its lowest eigenvalue is the minimum of a Rayleigh quotient — yielding rigorous upper bounds from any trial function.

where: $\mathcal{R}[y]$ — the Rayleigh quotient; its infimum over admissible $y$ is $\lambda_1$. Reading it: any admissible shape plugged in can only overestimate the fundamental — the basis of the Rayleigh–Ritz method. Finally, completeness is a closure relation: the eigenfunctions resolve the identity, and the same sum weighted by $1/(\lambda_n - \lambda)$ builds the Green's function whose poles are the spectrum.

where: $\delta$ — the Dirac delta; $G$ — the resolvent kernel of $L - \lambda$. Reading it: completeness is the delta identity; the Green's function is the same eigenbasis with energy denominators, tying boundary-value problems to the spectrum.

The character of the spectrum — discrete or continuous, real or complex, complete or not — is set by the regularity of $p$, $w$ and the interval, and by whether the boundary conditions are self-adjoint.

1 · Heat flow on a rod, solved by eigenfunction expansion (fully worked)

An insulated rod of length $L = 1$ (units where the diffusivity $\kappa = 1$) is held at zero temperature at both ends and starts with the parabolic profile $u(x, 0) = x(1 - x)$. Find the exact modal solution, the leading coefficients, and the times over which each mode dies.

Step 1. Separate $u(x,t) = X(x)T(t)$ in $u_t = u_{xx}$: the spatial factor solves $X'' + \lambda X = 0$ with $X(0) = X(1) = 0$ — Derivation 3 with $L = 1$, giving $\lambda_n = n^2\pi^2$ and $X_n = \sqrt{2}\,\sin(n\pi x)$; the time factor decays as $e^{-\lambda_n t}$, so $u(x,t) = \sum_n b_n \sin(n\pi x)\,e^{-n^2\pi^2 t}$.

Step 2. Project the data onto the modes (equation 6, $w = 1$): $b_n = 2\int_0^1 x(1-x)\sin(n\pi x)\,\mathrm{d}x$. This is $4/(n^3\pi^3)$ for odd $n$ and vanishes for even $n$ (the data is symmetric about $x = 1/2$, killing the antisymmetric modes), so $\mathbf{b_n = 8/(n^3\pi^3)}$ for odd $n$.

Step 3. In numbers: $b_1 = 8/\pi^3 = \mathbf{0.2580}$, $b_3 = 8/(27\pi^3) = \mathbf{0.009556}$, $b_5 = \mathbf{0.002064}$. The coefficients fall off like $n^{-3}$, so the parabola is nearly a single sine arch, as the two-term check $b_1 - b_3 = 0.2485$ against the exact peak $u(0.5,0) = 0.25$ confirms.

Step 4. Decay rates: mode $n$ dies with e-folding time $\tau_n = 1/(n^2\pi^2)$, so $\tau_1 = \mathbf{0.1013}$ while $\tau_3 = \mathbf{0.01126}$ — the third mode fades nine times faster. After a time of order $\tau_1 \approx 0.10$ the profile is, to better than 1%, the pure fundamental $0.2580\,\sin(\pi x)\,e^{-\pi^2 t}$ — why smooth initial data relaxes to its lowest mode.

2 · A Rayleigh–Ritz bound on the fundamental (fully worked)

Without solving the eigenproblem, estimate the lowest eigenvalue of $y'' + \lambda y = 0$ on $[0,1]$ with $y(0) = y(1) = 0$ from a single trial function, and compare with the exact value.

Step 1. Choose the simplest admissible shape: $\phi(x) = x(1 - x)$, vanishing at both ends. With $p = w = 1$, $q = 0$ and $\phi$ meeting the Dirichlet conditions, the boundary term in (7) drops and $\mathcal{R}[\phi] = \int_0^1 \phi'^2\,\mathrm{d}x \big/ \int_0^1 \phi^2\,\mathrm{d}x$.

Step 2. Compute the two integrals. Numerator: $\phi' = 1 - 2x$, so $\int_0^1 (1-2x)^2\,\mathrm{d}x = 1/3$. Denominator: $\int_0^1 x^2(1-x)^2\,\mathrm{d}x = 1/3 - 1/2 + 1/5 = 1/30$.

Step 3. Divide: $\mathcal{R}[\phi] = (1/3)/(1/30) = \mathbf{10}$. The exact fundamental is $\lambda_1 = \pi^2 = \mathbf{9.8696}$, so the crude one-parameter trial overestimates by only 1.3% — and by the variational principle it must be an upper bound. A parabola guessed in one line pins the ground-state eigenvalue to better than 2%, which is why Rayleigh–Ritz is the workhorse for drums, quantum wells and molecular energies where exact solutions are unavailable.

3 · Reading off the weight from a singular equation

Bessel's equation $x^2 y'' + x y' + (x^2 - \nu^2)y = 0$ governs cylindrical modes. Cast the eigenvalue version into S–L form and identify the weight.

Step 1. Write the parametric form on $[0, a]$: $x^2 y'' + x y' + (\lambda x^2 - \nu^2)y = 0$, with $\lambda$ the eigenvalue. Divide by $x$ to expose the structure: $x y'' + y' + (\lambda x - \nu^2/x)y = 0$.

Step 2. Recognize $x y'' + y' = (x y')'$ exactly — already divergence form with $p = x$. Matching to (1) gives $\mathbf{p = x}$, $q = -\nu^2/x$, and $\mathbf{w = x}$.

Step 3. Because $p(0) = 0$, the endpoint $x = 0$ is singular and boundedness of $y$ there replaces an explicit boundary value. Orthogonality of the Bessel modes $J_\nu(\alpha_{\nu k} x / a)$ is therefore weighted by $w = x$: $\int_0^a J_\nu(\alpha_{\nu j}x/a)\,J_\nu(\alpha_{\nu k}x/a)\,x\,\mathrm{d}x = 0$ for $j \neq k$ — the "$r\,\mathrm{d}r$" every student meets in cylindrical integrals is the S–L weight, not an accident of the area element.

Sturm–Liouville theory is where two nineteenth-century friends built, almost simultaneously, the first infinite-dimensional spectral theory — decades before Hilbert spaces existed to name it. The result reads today like linear algebra done in a space of functions.

Texts. Stone & Goldbart, Mathematics for Physics, Ch. 8 (Linear Differential Operators) · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 8 · Hassani, Mathematical Physics, Ch. 19–20 (Sturm–Liouville systems) · Byron & Fuller, Mathematics of Classical and Quantum Physics, Ch. 5 · Courant & Hilbert, Methods of Mathematical Physics, Vol. I, Ch. V–VI · Coddington & Levinson, Theory of Ordinary Differential Equations, Ch. 7–8.

Provenance. All results are standard Sturm–Liouville theory. A verifier should recompute: (1) the integrating-factor reduction $p = \exp\int(a_1/a_2)\,\mathrm{d}x$, $w = p/a_2$; (2) the Lagrange boundary-term cancellation for separated and periodic conditions, hence reality and orthogonality; (3) the string spectrum $\lambda_n = n^2\pi^2/L^2$ with $\int_0^L \sin^2(n\pi x/L)\,\mathrm{d}x = L/2$; (4) the heat-rod coefficients $b_n = 8/(n^3\pi^3)$ odd ($b_1 = 0.2580$) and decay times $\tau_n = 1/(n^2\pi^2)$ ($\tau_1 = 0.1013$); (5) the Rayleigh bound $\mathcal{R}[x(1-x)] = 10$ versus exact $\pi^2 = 9.8696$; (6) the Bessel weight $w = x$. No nonstandard constants are asserted.

Prerequisites. Boundary Conditions and Well-Posedness · Separation of Variables in Cartesian Coordinates · Classification of Second-Order Linear PDEs. Leads to. Separation of Variables in Spherical Coordinates · Separation of Variables in Cylindrical Coordinates · Green's Functions for PDEs · The Helmholtz Equation.

Ask the tutor: ask me to reduce any specific ODE — Bessel, Legendre, Hermite, or one from your own problem set — to self-adjoint form and read off its weight, spectrum, and orthogonality relation step by step.