Matrix Hamiltonians and numerical eigenvalues

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 as a complex array, with nothing to guarantee that the imaginary parts vanish or that the eigenvectors are orthogonal. 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.

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.

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/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.

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 set by the box, with rungs roughly ħ²π²j²/(2mW²) above V(∞). 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(V(∞) − 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.

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.

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.