University Physics V · The Schrödinger Equation · 4.9
Matrix Hamiltonians & Numerical Eigenvalues
Almost no potential you will meet has a closed-form spectrum, so you will diagonalise. This lesson is about doing it honestly: naming which matrix you actually built, why it is symmetric, and how to quote an eigenvalue with an error bar instead of fifteen digits of borrowed confidence.
Build the model
Connect the measurement to the mechanism.
Every eigenvalue a computer returns belongs to a different operator from the one you wrote down, and the whole skill is knowing which way it differs. Sampling ψ on a grid replaces d²/dx² by a three-point difference, whose eigenvalue on the mode e(ikx) is not −k² but −(4/h²)sin²(kh/2); that is smaller in magnitude for every k, short by k⁴h²/12, so grid levels sit below the truth and the deficit grows as the fourth power of the wavenumber. Expanding instead in N basis functions and diagonalising ⟨φₘ|H|φₙ⟩ does the opposite: you are minimising the true quadratic form over a subspace of the true Hilbert space, so by the Rayleigh–Ritz theorem every root is an upper bound, falling monotonically as the basis grows.
Cutting the domain at finite walls is a third and independent lie, turning the continuum above the well into a ladder of invented levels and squeezing the genuine ones upward. Both matrices are real symmetric whenever V is real, which is why scipy.linalg.eigh rather than eig is the right call and why the eigenvectors come back orthogonal for free. What all this costs is a second run: the method of this subject is not "diagonalise" but "diagonalise again with a different h, a different basis size and a wider box, and quote only the digits that did not move."
- Simple definition
- A matrix Hamiltonian is the finite real symmetric matrix you get by sampling H on a grid or projecting it onto N basis states, so that H|ψ⟩ = E|ψ⟩ becomes an ordinary matrix eigenvalue problem a computer solves exactly.
- Example
- On [0, 1] with ħ = m = 1, V = 0 and four interior points (h = 0.2), H is the 4×4 tridiagonal matrix with 25 on the diagonal and −12.5 beside it; its lowest eigenvalue is 4.7746 against the exact π²/2 = 4.9348.
Real, symmetric and tridiagonal: hand it to eightridiagonal at O(N²) rather than dense eigh at O(N³).
xⱼ = xₘᵢₙ + jh for j = 1…N; h the grid spacing in m, V in J, entries in J
The grid's own kinetic spectrum. It saturates at 2ħ²/mh², so no eigenvalue above that ceiling can exist at all.
k in rad m⁻¹, h in m; the exact operator gives ħ²k²/2m, equal only as kh → 0
Halving h quarters the error, and every level is biased low, never high — the grid gives no upper bound.
kₙ = nπ/L for the box; the sign is fixed, the deficit is never negative
Two runs buy an h⁴ answer: 4.77458 and 4.89435 combine to 4.93427 against the exact 4.93480.
E(h) the eigenvalue computed on spacing h, in J; valid while the h² term dominates
Basis truncation errs upward while the grid errs downward, so running both brackets the exact level.
Eₖ(N) is the k-th root of the N×N matrix ⟨φₘ|H|φₙ⟩ over a nested basis
Pass S to scipy.linalg.eigh(H, S), or the roots are not variational bounds on anything.
S is dimensionless and positive definite; S = I only for an orthonormal basis
From a wavefunction to a column vector
Choose a window, place N interior points at xⱼ = xₘᵢₙ + jh with h = (xₘₐₓ − xₘᵢₙ)/(N + 1), and let ψⱼ = ψ(xⱼ). The three-point second difference ψ″(xⱼ) ≈ (ψ_(j−1) − 2ψⱼ + ψ_(j+1))/h² turns −ħ²ψ″/2m + Vψ = Eψ into a matrix equation with diagonal ħ²/(mh²) + V(xⱼ) and both off-diagonals −ħ²/(2mh²). The walls enter as ψ₀ = ψ_(N+1) = 0, which is why those terms simply drop out of the first and last rows: the Dirichlet condition is built into the matrix, not imposed on the answer afterwards. For a real V the matrix is real symmetric, so scipy.linalg.eigh returns real eigenvalues and an orthogonal eigenvector matrix by construction; scipy.linalg.eig returns the same numbers carrying imaginary parts of order 10⁻¹⁶ and no orthogonality guarantee. Because this H is tridiagonal, scipy.linalg.eigh_tridiagonal costs O(N²) against dense eigh's O(N³), which on N = 4000 is the difference between a second and a minute.
What the three-point Laplacian actually diagonalises
Feed the difference operator a plane wave ψⱼ = e(ikxⱼ). Then (ψ_(j−1) − 2ψⱼ + ψ_(j+1))/h² = ψⱼ(e(−ikh) − 2 + e(ikh))/h² = −ψⱼ (4/h²) sin²(kh/2). So the grid's kinetic energy is (2ħ²/mh²)sin²(kh/2), not ħ²k²/2m. Expanding, (4/h²)sin²(kh/2) = k² − k⁴h²/12 + …, so every mode is short by ħ²k⁴h²/24m: the bias is negative and grows as the fourth power of the wavenumber. Two consequences you should carry away as numbers. Demanding 1% accuracy needs k²h²/12 < 0.01, i.e. kh < 0.35, which is about 18 grid points per de Broglie wavelength — the answer to "how fine a grid?". And the kinetic energy saturates at 2ħ²/mh², a hard ceiling. For the box on [0, L] the result is exact rather than asymptotic, Eₙ(h) = (2ħ²/mh²)sin²(nπh/2L), so on N = 12 points the twelfth level comes back 53% low while the first is only 0.49% low.
Basis truncation errs the other way, and provably
The alternative is to keep the operator and shrink the space: pick N functions φ₁…φN obeying the boundary conditions, form Hₘₙ = ⟨φₘ|H|φₙ⟩ and Sₘₙ = ⟨φₘ|φₙ⟩, and solve Hc = ESc. This is Rayleigh–Ritz, and the Hylleraas–Undheim–MacDonald theorem says the k-th root is an upper bound on the k-th exact eigenvalue and decreases monotonically as functions are added, because you are minimising the same quadratic form over a subspace of the true Hilbert space. On [0, 1] with ħ = m = 1, the single trial function x(1 − x) gives ⟨H⟩/⟨1⟩ = (1/6)/(1/30) = 5.0000 against the exact 4.93480 — high by 1.3%, as it must be. Add x²(1 − x)² and the 2×2 generalised problem drops it to 4.93487. Two warnings. With a non-orthogonal basis you must pass S, or the roots bound nothing; and a nearly linearly dependent basis makes S ill-conditioned, at which point the extra digits are numerical noise rather than physics.
The walls invent a spectrum of their own
A grid needs a finite window, and that window is a third approximation entirely independent of h. Put walls a distance W apart around a well supporting, say, three bound states, and eigh will hand back several hundred eigenvalues. Only those below V(∞) are candidates for bound states; everything above is the continuum chopped into a ladder with spacing set by the box, roughly ħ²π²j²/(2mW²). The test is not to stare at the number but to move the wall: a genuine bound state's energy stops moving once W exceeds a few decay lengths κ⁻¹ = ħ/√(2m|E|), while an artefact level falls as 1/W². The walls also bias the genuine levels, and they bias them upward, because clipping the exponential tail is extra confinement. So the grid pushes levels down and the box pushes them up, and for one unlucky pair of h and W the two can cancel: a single run cannot tell you which error you are looking at.
Convergence is a two-parameter study, not just a fine grid
Report an eigenvalue the way you report a measurement. Halve h with the box fixed and watch the change: for the three-point Laplacian E(h) = E − Ch² + O(h⁴), so the error should fall by a factor of four, and if it does you have earned the right to Richardson-extrapolate, E ≈ [4E(h/2) − E(h)]/3. On the box that took the ground state from 4.89435 at h = 0.1 (error 4 × 10⁻²) to 4.93427 (error 5 × 10⁻⁴) for no extra matrix. Then widen the box with h fixed and watch again. The two studies must be run separately precisely because their errors have opposite signs. Quote only the digits that survive both. Three cheap assertions belong in the same script: check the matrix equals its own transpose to machine precision, count the nodes of eigenvector n and confirm you get n − 1, and evaluate the residual ‖Hψ − Eψ‖ instead of trusting the solver's return value.
Reading the eigenvector the solver hands back
eigh returns eigenvectors as columns of unit Euclidean norm, Σⱼ |ψⱼ|² = 1, which is not the quantum normalisation. The discrete inner product carries the measure: ⟨ψ|ψ⟩ = h Σⱼ |ψⱼ|², so divide the returned column by √h before computing ⟨x⟩ or overlaying an analytic ψ(x). The overall sign is arbitrary — eigh may equally hand you −ψ — so fix a convention, such as making the first nonzero component positive, before comparing two runs. Degeneracy needs the same care: inside a degenerate block the solver returns some orthonormal basis of the eigenspace, not the one your symmetry argument chose, so label states by a commuting observable rather than by whatever order came back. And a pair split by less than the grid error is not resolved at all, however many digits the printout shows.
Change one variable at a time
Make the relationship visible.
Drag N down to 12 and step n up: the grid tracks the exact parabola while there are more than about 18 points per wavelength, then peels away, the top level landing at under half its true energy. Push N to 60 and that same level climbs back to within a few per cent.
EXACT Eₙ493.48 ħ²/mL²
GRID Eₙ408.04 ħ²/mL²
GRID DEFICIT17.31 % low
POINTS PER WAVELENGTH4.2
Live interpretationEXACT Eₙ: 493.48 ħ²/mL². GRID Eₙ: 408.04 ħ²/mL². GRID DEFICIT: 17.31 % low. POINTS PER WAVELENGTH: 4.2
Catch the common trap
Explain before calculating.
The infinite well on [0, L] is solved by diagonalising a three-point finite-difference Hamiltonian on a uniform grid of N interior points, spacing h = L/(N + 1). Compared with the exact Eₙ = n²π²ħ²/(2mL²), where do the computed levels sit, and which of them are worst?
Choose an answer to test the model.
Practice & worked examples
Reason from the model, then test the result.
EasyTake the infinite well on [0, 1] with ħ = m = 1 and V = 0 inside, and put four interior grid points on it. Write down the 4×4 finite-difference Hamiltonian, then find its lowest eigenvalue using the standard result that a tridiagonal matrix of size N with diagonal d and off-diagonals −c has eigenvalues d − 2c cos(jπ/(N + 1)). Compare with the exact ground state.
- Grid: h = L/(N + 1) = 1/5 = 0.2, with interior points at 0.2, 0.4, 0.6 and 0.8. The walls give ψ₀ = ψ₅ = 0, so those terms drop out of the first and last rows and the matrix is 4×4, not 6×6 — the Dirichlet condition is built into H, not applied to the answer later.
- Entries: d = ħ²/(mh²) + V = 1/0.04 = 25 and c = ħ²/(2mh²) = 12.5, so H = tridiag(−12.5, 25, −12.5). It is real and symmetric, so eigh applies and the four eigenvalues are real by construction.
- Lowest eigenvalue: j = 1 gives E₁ = 25 − 25 cos(π/5) = 25(1 − 0.80902) = 4.7746.
- Exact: E₁ = π²ħ²/(2mL²) = π²/2 = 4.9348. The grid is low by 0.1602, or 3.2%, and it is low rather than high because (1/h²)(1 − cos πh) = (2/h²)sin²(πh/2) < π²/2 for every h > 0.
AnswerH = tridiag(−12.5, 25, −12.5); E₁ = 4.7746 against the exact π²/2 = 4.9348 — 3.2% low, and low because the three-point Laplacian underestimates k².
MediumKeep the same box but use nine interior points, h = 0.1, with the closed form Eₙ(h) = (2ħ²/mh²)sin²(nπh/2L) for the grid eigenvalues. Show that the error in E₁ scales as h², combine the two grids into a better value, then work out what the same grid does to the top level, n = 9.
- At h = 0.1: E₁ = 200 sin²(0.05π) = 200 × 0.0244717 = 4.89435.
- Errors: at h = 0.2 the deficit was 4.93480 − 4.77458 = 0.16023; at h = 0.1 it is 4.93480 − 4.89435 = 0.04045. The ratio is 3.96, so halving h has quartered the error — the h² law that ħ²k⁴h²/24m predicts (0.16235 and 0.04059 from the formula).
- Richardson: with E(h) = E − Ch², the combination [4E(h/2) − E(h)]/3 kills the h² term. (4 × 4.89435 − 4.77458)/3 = 14.80282/3 = 4.93427, which is 5 × 10⁻⁴ from the exact 4.93480 instead of 4 × 10⁻².
- The top of the spectrum is another matter. For n = 9, nπh/2 = 1.41372 rad, so E₉ = 200 sin²(1.41372) = 200 × 0.97553 = 195.11 against the exact 81π²/2 = 399.72 — 51% low. There are 2(N + 1)/n = 20/9 = 2.2 grid points per wavelength there, and the h² expansion failed long before.
AnswerE₁(0.1) = 4.89435, error 0.0405 — a quarter of the h = 0.2 error; Richardson gives 4.93427 against the exact 4.93480. But E₉ = 195.1 against 399.7, 51% low, at 2.2 points per wavelength.
HardBracket the ground state of the same box from both sides, with ħ = m = L = 1. (a) Take the h = 0.1 grid value as the low estimate. (b) Compute the Rayleigh quotient of the trial function φ(x) = x(1 − x), which satisfies both boundary conditions. (c) Explain why the two numbers must straddle the exact energy, verify that they do, and then tighten the upper bound with a second trial function φ₂(x) = x²(1 − x)².
- Lower side: E₁(h = 0.1) = 200 sin²(0.05π) = 4.89435. Because (4/h²)sin²(kh/2) ≤ k² on every Fourier mode, and the box eigenfunctions are exactly those modes, this genuinely lies below the exact answer — a statement that holds here because V = 0, not one to assume for a general potential.
- Upper side: E[φ] = ⟨φ|H|φ⟩/⟨φ|φ⟩ with H = −½ d²/dx². Since φ(0) = φ(1) = 0, integrating by parts gives ⟨φ|H|φ⟩ = ½∫₀¹(φ′)² dx with φ′ = 1 − 2x, so ½∫₀¹(1 − 4x + 4x²) dx = ½(1 − 2 + 4/3) = 1/6.
- Norm: ∫₀¹ x²(1 − x)² dx = 1/3 − 1/2 + 1/5 = 1/30. Hence E[φ] = (1/6)/(1/30) = 5.0000 exactly.
- Why they straddle: the Rayleigh quotient minimises the true H over a one-dimensional subspace of the true space, so it can only exceed the true minimum, while the grid replaces H by a strictly softer operator, so it can only fall short. Check: 4.89435 < π²/2 = 4.93480 < 5.00000, a bracket 0.10565 wide with the exact answer 38% of the way across it.
- Tighten it: add φ₂ = x²(1 − x)². Then H = [[1/6, 1/30], [1/30, 1/105]] and S = [[1/30, 1/140], [1/140, 1/630]], and det(H − ES) = 0 reduces to E² − 56E + 252 = 0. Its lower root, (56 − √2128)/2 = 4.93487, is still above π²/2 = 4.93480, but now by only 1.5 parts in 10⁵.
Answer4.89435 (grid, biased low) and 5.00000 (one-term Rayleigh-Ritz, rigorously high) bracket π²/2 = 4.93480; a second trial function tightens the upper bound to 4.93487.