Differential equations as eigenvalue problems

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.

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.

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.

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.

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.

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.