University Physics V · Mathematical Foundations of Quantum Physics · 1.9

Differential equations as eigenvalue problems

Every bound-state problem you will meet is this one in costume: a second-order ODE, two homogeneous boundary conditions, and the question of which λ let anything survive. Learn to spot the Sturm–Liouville shape, let the boundary conditions do the quantising, and trade the operator for a matrix when the solvable shortlist runs out.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

Ask what −u″ = λu means on an interval and the honest answer is: nothing yet. For every λ the equation has two independent solutions, so the differential expression alone has no spectrum. Impose u(0) = u(L) = 0 and almost every λ dies; only λₙ = (nπ/L)² survives, one sine per integer.

The eigenvalue problem is therefore the pair — operator and boundary conditions — and its general shape is Sturm–Liouville: −(pu′)′ + qu = λwu. Integrate by parts twice and the operator's symmetry comes down to a single boundary term, [p(u*′v − u*v′)]; boundary conditions that kill it make the problem self-adjoint, and self-adjointness buys the whole toolkit this unit has assembled — real eigenvalues, eigenfunctions orthogonal under the weight w, and a complete basis to expand on. The cost is that closed-form solutions stop at a shortlist of famous equations.

Beyond it you trade the ODE for a matrix: replace u″ by its central difference on N interior points and the operator becomes a real symmetric tridiagonal matrix whose numpy.linalg.eigh spectrum approximates the true one — always from below, with a deficit growing as (kπh/L)²/12, because a grid of N points holds exactly N modes and bends the stiff ones flat.

Simple definition: A differential eigenvalue problem asks for the values λ at which an ODE plus homogeneous boundary conditions admits a nonzero solution; the boundary conditions, not the equation, select the discrete spectrum.
Example: −u″ = λu with u(0) = u(1) = 0 has nonzero solutions only at λₙ = n²π² = 9.87, 39.5, 88.8, … (in m⁻² for x in metres), one eigenfunction √2 sin(nπx) for each integer n.

Sturm–Liouville form−(p(x) u′)′ + q(x) u = λ w(x) u

Any a₀u″ + a₁u′ + a₂u = λu reaches this shape via an integrating factor — and the form hands you the weight for orthogonality.

p(x) > 0 and w(x) > 0 on [a, b]; for −u″ = λu on a metre interval, λ carries m⁻²

Dirichlet spectrum of −d²/dx²λₙ = (nπ/L)², uₙ(x) = √(2/L) sin(nπx/L)

The model problem. The integer n is manufactured by the boundary conditions, not by the ODE.

n = 1, 2, 3, …; ⟨uₘ|uₙ⟩ = δₘₙ on [0, L]

Lagrange identity boundary term⟨u|Lv⟩ − ⟨Lu|v⟩ = [p(u*′v − u*v′)]ₐᵇ

Self-adjointness is decided here: Dirichlet, Neumann, Robin and periodic conditions all make the bracket vanish.

for L = −(pu′)′ + qu with real p and q, evaluated between the endpoints a and b for any admissible pair u, v

Central-difference matrix(Hu)ⱼ = (−uⱼ₋₁ + 2uⱼ − uⱼ₊₁)/h²

A real symmetric tridiagonal matrix — exactly what numpy.linalg.eigh is built for.

grid xⱼ = jh, j = 1 … N, h = L/(N+1); Dirichlet built in because the endpoint values u(0) = u(L) = 0 are never unknowns

Spectrum of the difference matrixλₖʰ = (4/h²) sin²(kπh/2L) ≤ (kπ/L)²

Every grid eigenvalue sits low, and the error grows as k² — trust the bottom of an eigh spectrum, never its top.

k = 1 … N; fractional deficit ≈ (kπh/L)²/12, second order in h

Read the ODE as an operator, and λ as its spectrum

Solve −u″ = λu and you cannot: for every λ ≠ 0 there are two independent solutions, sin(√λ x) and cos(√λ x), and at λ = 0 the pair 1 and x, so nothing has been asked yet. The real question arrives with homogeneous boundary conditions. Demand u(0) = 0 and the cosine dies; demand u(L) = 0 as well and sin(√λ L) = 0 forces √λ L = nπ. Out of a continuum of candidates the interval keeps λₙ = (nπ/L)² and discards the rest — for L = 1 m that is 9.87, 39.48, 88.83 m⁻², the harmonics of a string and, dressed with ħ²/2m, the levels of the infinite well. This is the template for every eigenvalue problem in the course: a linear operator L, a space of functions the boundary conditions carve out, and the discrete set of λ at which L − λ fails to be invertible. Quantisation is not a postulate bolted onto wave mechanics; it is what homogeneous boundary conditions do to a differential operator.

The integrating factor that reveals Sturm–Liouville form

Most equations do not arrive as −(pu′)′ + qu = λwu; they arrive as a₀(x)u″ + a₁(x)u′ + a₂(x)u = λu, with the first-derivative term spoiling the symmetry. One integrating factor repairs it: multiply through by (1/a₀)exp(∫a₁/a₀ dx) and the first two terms collapse into the single derivative (pu′)′ with p = exp(∫a₁/a₀ dx). Hermite's equation u″ − 2xu′ = −2nu becomes −(e^(−x²) u′)′ = 2n e^(−x²) u, so p = e^(−x²) and — crucially — the same exponential appears as the weight w. That is the point of the exercise: the weight names the inner product the problem lives in, ⟨f|g⟩ = ∫ f*g w dx, and it is in that inner product, not the bare L² one, that eigenfunctions of different eigenvalues are orthogonal. Miss the weight and the Hermite polynomials look non-orthogonal; keep it and ∫ Hₘ Hₙ e^(−x²) dx = 0 for m ≠ n drops out for free.

Self-adjointness lives in the boundary term

Integrate ⟨u|Lv⟩ by parts twice for L = −(d/dx)(p d/dx) + q and every interior term cancels, leaving one surface term: ⟨u|Lv⟩ − ⟨Lu|v⟩ = [p(u*′v − u*v′)]ₐᵇ. The differential expression is symmetric by construction; self-adjointness is decided by whether the boundary conditions kill this bracket for every admissible pair. Dirichlet u = 0, Neumann u′ = 0, and any real Robin combination αu + βu′ = 0 kill it end by end; periodic conditions u(a) = u(b), u′(a) = u′(b) kill it jointly, provided p(a) = p(b). Change the conditions and the spectrum changes with them: on [0, 1] the same −d²/dx² has eigenvalues (nπ)² under Dirichlet, gains λ = 0 under Neumann, and under periodic conditions keeps λ = 0 but turns every (2nπ)² doubly degenerate. And when no choice of conditions can kill the bracket on the given domain — the momentum operator −iħ d/dx on the half-line x ≥ 0 is the standard case — the operator is symmetric but not self-adjoint, and the spectral theorem's promise of real eigenvalues and a complete eigenbasis is off the table.

Discretise: the operator becomes a matrix for eigh

When the equation leaves the solvable shortlist, put it on a grid. Take xⱼ = jh with h = L/(N + 1) and keep only the N interior values uⱼ, j = 1 … N: the Dirichlet conditions are not applied afterwards but built in, because the endpoint values u(0) and u(L) never appear as unknowns. The central difference u″(xⱼ) ≈ (uⱼ₊₁ − 2uⱼ + uⱼ₋₁)/h² then turns −d²/dx² into the tridiagonal matrix with 2/h² down the diagonal and −1/h² beside it — real, symmetric, and exactly what numpy.linalg.eigh exists for: real eigenvalues in ascending order, orthonormal eigenvectors as columns. This model problem can even be diagonalised by hand: the trial vector uⱼ = sin(kπjh/L) and the identity sin(A(j+1)) + sin(A(j−1)) = 2 cos A · sin(Aj) give λₖʰ = (2/h²)(1 − cos(kπh/L)) = (4/h²) sin²(kπh/2L). On L = 1 with N = 4 (h = 0.2) the lowest eigenvalue is 100 sin²(0.1π) = 9.549 against the true π² = 9.870 — the machinery works, and its bias is already visible.

Trust the bottom of the spectrum, never the top

The grid is a basis with N vectors, so the matrix owns exactly N modes, and mode k — which draws k half-waves across the interval — is sampled by only (N+1)/k points per half-wave. Comparing spectra shows what the sampling costs: the true (kπ/L)² against the computed (4/h²) sin²(kπh/2L). Since sin θ < θ, every computed eigenvalue sits below its target, and expanding the sine gives the deficit: λₖʰ ≈ (kπ/L)²[1 − (kπh/L)²/12]. Two lessons are folded in. First, the error is second order in h: halve the spacing and each fixed mode's deficit quarters, which is the slope of 2 that a log-log convergence plot should show. Second, the error grows as k², so the top of the spectrum is junk however fine the grid — at k = N the computed value falls 40–60% short. A working rule: ten points per half-wave gives a deficit of (π/10)²/12 ≈ 0.8%. Demand that many for the highest mode you intend to quote, then double N and check the quoted ones stop moving.

This is the Schrödinger equation's shape

The reason this sits in week one: the time-independent Schrödinger equation −(ħ²/2m)ψ″ + V(x)ψ = Eψ is a Sturm–Liouville problem with p = ħ²/2m, q = V and weight w = 1, and bound states are the eigenpairs selected by ψ → 0 at the walls. Its discretisation is the same tridiagonal kinetic matrix plus V evaluated on the grid down the diagonal — one extra line of NumPy. Every guarantee in this lesson transfers: Hermitian, so real energies; boundary term killed, so orthogonal states you may expand in; a discrete spectrum because the boundary conditions demand one. And every warning transfers too: a finite box stands in for infinity, so states that leak to the wall are corrupted along with the high-k modes the grid cannot bend. When the box of Unit 5, the finite wells of Unit 6 and the oscillator of Unit 7 arrive, the solver is already on the table; all that changes is the diagonal.

Change one variable at a time

Make the relationship visible.

Interactive model

Interior grid points N12

Mode number k12

Start at k = N = 12, where the computed eigenvalue sits 53% low, then raise N at fixed k = 12 and watch the filled marker climb toward the parabola — 18% low by N = 24 — while sliding k down shows the bottom of the spectrum is already sound.

TRUE λ = (kπ)²1421

GRID λ FROM eigh666

DEFICIT53.1 %

POINTS PER HALF-WAVE1.1

Live interpretationTRUE λ = (kπ)²: 1421. GRID λ FROM eigh: 666. DEFICIT: 53.1 %. POINTS PER HALF-WAVE: 1.1

Catch the common trap

Explain before calculating.

Misconception check

Students treat boundary conditions as bookkeeping — solve the equation, then use the conditions to fix the constants — as if the spectrum belonged to the differential expression alone. It does not. For every λ the ODE has a full set of solutions, so the expression by itself selects nothing: no boundary conditions, no eigenvalue problem at all. Nor do the conditions merely trim a spectrum that was already there; they define the operator. Dirichlet, Neumann and periodic conditions on the one interval [0, 1] turn the one formula −d²/dx² into three different self-adjoint operators with three different spectra — the lowest eigenvalue is π² for the first and 0 for the other two, and only the third has degenerate levels. The boundary conditions that quantise the spectrum are the same ones that decide, through the bracket [p(u*′v − u*v′)], whether the operator is self-adjoint at all, which is why the physics convention writes them next to the operator, never after the solution.

You discretise −u″ = λu on [0, 1] with u(0) = u(1) = 0 on N = 49 interior points (h = 0.02) and call `numpy.linalg.eigh` on the tridiagonal matrix. What comes back for mode k = 40, whose true eigenvalue is (40π)² ≈ 15 791?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyFind the three lowest eigenvalues and normalised eigenfunctions of −u″ = λu on [0, 1 m] with u(0) = u(1) = 0, and state exactly which step of the calculation created the integer n.

General solution for λ > 0: u = A sin(√λ x) + B cos(√λ x). (λ ≤ 0 gives lines and real exponentials, which cannot vanish at both ends without vanishing everywhere.)
u(0) = 0 forces B = 0 — the first boundary condition halves the solution space but quantises nothing.
u(1) = 0 with A ≠ 0 demands sin √λ = 0, so √λ = nπ: the second boundary condition is the quantiser.
λₙ = n²π² = 9.87, 39.48, 88.83 m⁻²; normalising ∫₀¹ A² sin²(nπx) dx = A²/2 = 1 gives A = √2, so uₙ = √2 sin(nπx).

Answerλ₁ = π² ≈ 9.87 m⁻², λ₂ = 4π² ≈ 39.48 m⁻², λ₃ = 9π² ≈ 88.83 m⁻², with uₙ = √2 sin(nπx). The integer exists because sin √λ must vanish at x = 1 — quantisation enters at the second boundary condition.

MediumLaguerre's equation xu″ + (1 − x)u′ + ν = 0 is the simplest member of the family whose associated form governs the radial part of hydrogen. Bring it to Sturm–Liouville form, identify p, q, w and λ, and state the orthogonality relation the form guarantees on [0, ∞).

Match to a₀u″ + a₁u′ + a₂u = 0: a₀ = x, a₁ = 1 − x. The first-derivative term blocks symmetry, so build the integrating factor from a₁/a₀.
∫(1 − x)/x dx = ln x − x, so p = exp(∫a₁/a₀ dx) = x e⁻ˣ; multiply the equation by p/a₀ = e⁻ˣ to get x e⁻ˣ u″ + (1 − x) e⁻ˣ u′ + n e⁻ˣ u = 0.
The first two terms are now an exact derivative, since (x e⁻ˣ u′)′ = x e⁻ˣ u″ + (1 − x) e⁻ˣ u′. The equation reads −(x e⁻ˣ u′)′ = n e⁻ˣ u.
Read off the Sturm–Liouville data: p = x e⁻ˣ, q = 0, w = e⁻ˣ, λ = n. Note p(0) = 0 and p, w → 0 at infinity: the problem is singular, and the endpoint conditions are boundedness at 0 and no faster than polynomial growth at infinity, not Dirichlet — that is enough to make p(u*′v − u*v′) vanish at both ends.
Self-adjointness under the weight w delivers ∫₀^∞ Lₘ(x) Lₙ(x) e⁻ˣ dx = 0 for m ≠ n — the measure e⁻ˣ dx comes out of the form, not from a separate calculation.

Answer−(x e⁻ˣ u′)′ = n e⁻ˣ u, so p = x e⁻ˣ, q = 0, w = e⁻ˣ, λ = n, and the Laguerre polynomials obey ∫₀^∞ Lₘ Lₙ e⁻ˣ dx = 0 for m ≠ n.

HardDiscretise −u″ = λu on [0, 1] with Dirichlet conditions on N = 9 interior points. Using λₖʰ = (4/h²) sin²(kπh/2), evaluate the lowest computed eigenvalue and its fractional error, verify the (kπh)²/12 deficit rule, and find the N that brings λ₁'s error below 0.01%.

h = 1/(N + 1) = 0.1, so the matrix is tridiag(−100, 200, −100) and its exact spectrum is λₖʰ = (4/h²) sin²(kπh/2) = 400 sin²(kπ/20).
k = 1: λ₁ʰ = 400 sin²(π/20) = 400 × (0.15643)² = 400 × 0.024472 = 9.7887, against the true π² = 9.8696.
Fractional deficit: (9.8696 − 9.7887)/9.8696 = 0.0082, i.e. 0.82% low. The rule predicts (kπh)²/12 = (0.1π)²/12 = 0.098696/12 = 0.0082 — agreement in the second significant figure.
Second-order convergence: the deficit scales as h², so doubling the grid to h = 0.05 (N = 19) should quarter it. Check: 1600 sin²(π/40) = 9.8493, which is 0.205% low — a quarter of 0.82%.
Set (πh)²/12 = 10⁻⁴: h² = 12 × 10⁻⁴/π² = 1.216 × 10⁻⁴, h = 0.0110, so N + 1 = 1/h ≈ 90.7. Take N = 90 (h = 1/91): deficit = (π/91)²/12 = 9.9 × 10⁻⁵ = 0.0099%, whereas N = 89 would leave 0.0102%.

Answerλ₁ʰ = 9.7887, sitting 0.82% below π² — exactly the (πh)²/12 the expansion predicts. Halving h quarters the deficit, and N = 90 interior points push the error to 0.0099%, just under the 0.01% target.