University Physics V · The Quantum Wavefunction · 3.9
Discretised States & Matrix Operators
Nothing in this unit is computable until |ψ⟩ becomes a column of numbers. Sample it on a uniform grid and the integrals turn into dot products, the operators into matrices, and numpy.linalg.eigh into a legitimate method. The price is a box, a spacing, and a boundary condition you now have to defend.
Build the model
Connect the measurement to the mechanism.
Discretisation is a change of Hilbert space, not a numerical shortcut. Choose a box and N sample points a distance Δx apart, replace the vector |ψ⟩ by the column ψⱼ = ψ(xⱼ), and the L2 inner product ∫φ*ψ dx becomes the midpoint sum Δx φ†ψ: CN under that weight is the space you are now working in. Every operator inherits a matrix.
X is diagonal because multiplying by x acts pointwise, and its eigenvectors — the grid unit vectors — are this space's stand-in for the |x⟩ that never belonged to L2. P has no pointwise action, so it becomes a stencil: the centred difference gives a Hermitian tridiagonal matrix, the FFT gives exact multiplication by ħk in the transformed basis. Expectation values collapse to ⟨A⟩ = Δx ψ†Aψ and eigenvalue problems to numpy.linalg.eigh.
What it costs is everything the continuum kept in its boundary conditions. Truncating at N rows silently imposes ψ = 0 outside the box, or periodicity if you used the FFT; the spacing caps representable momentum at πħ/Δx and biases every kinetic eigenvalue low by (kΔx)²/12; and since Tr[A, B] = 0 for any two N×N matrices while Tr(iħI) = iħN, the canonical commutator cannot hold exactly at any N. In finite dimensions Hermitian and self-adjoint coincide, so the grid does not answer the domain question that made −iħ d/dx delicate on a half-line — it decides it for you, and you have to know which answer you wanted.
- Simple definition
- Discretisation replaces the state |ψ⟩ by its samples on a uniform grid, turning the Hilbert space into CN with inner product Δx ψ†φ, so every observable becomes an N×N Hermitian matrix and every expectation value a vector-matrix-vector product.
- Example
- An electron in a 1.00 nm box with N = 4 interior points has Δx = 0.20 nm, X = diag(0.20, 0.40, 0.60, 0.80) nm, and kinetic matrix (ħ²/2mΔx²)⋅tridiag(−1, 2, −1) whose prefactor is 0.952 eV.
Turns every integral in this unit into one weighted dot product.
ψⱼ = ψ(xⱼ) in m(−1/2), Δx in m; normalisation is Δx ψ†ψ = 1
Its eigenvectors eⱼ are the grid's normalisable stand-in for the δ-normalised |x⟩.
entries in m and real, so X = X† with no work; X² = diag(xⱼ²)
A forward difference leaves ‖P − P†‖ ~ ħ/Δx and returns an imaginary ⟨p⟩ on a real ψ.
entries in kg m s⁻¹; P = P† only because the stencil is centred
Exactly Hermitian and spectrally accurate, but it wires in periodic walls, not a box.
F the unitary DFT matrix; n from −N/2 to N/2 − 1; ħkₙ in kg m s⁻¹
Expands to (ħ²k²/2m)[1 − (kΔx)²/12], so a grid always biases levels low.
for an electron at Δx = 0.20 nm, ħ²/2mΔx² = 0.952 eV; band top 4× that
⟨[X, P]⟩ → iħ[1 − Δx²⟨p²⟩/2ħ²] on smooth states, but the operator identity never holds.
S± shift by one site; eigenvalues of [X, P] are iħcos(nπ/(N+1))
The grid is a quadrature rule, and it changes the space
Fix a box of length L and N interior points xⱼ = jΔx with Δx = L/(N+1). Sampling replaces |ψ⟩ by the column (ψ₁ … ψN), and ∫φ*ψ dx by the midpoint sum Δx Σⱼ φⱼ*ψⱼ = Δx φ†ψ. That weight is not decoration: ψ carries units of m(−1/2), so ψ†ψ has units of m⁻¹ and only Δx ψ†ψ is dimensionless. Normalisation becomes Δx Σⱼ |ψⱼ|² = 1, and Born's rule for an interval becomes P(a < x < b) = Δx Σ over the samples inside it. Two independent errors have already entered. The quadrature error in the inner product is O(Δx²) for a smooth ψ, and far smaller than that for one vanishing at both ends. The error in the operators, which is what the rest of this topic is about, is a separate account entirely. For L = 1.00 nm and N = 4, Δx = 0.20 nm and the samples sit at 0.20, 0.40, 0.60 and 0.80 nm: four complex numbers standing in for a function.
X is diagonal, and that is all the position basis ever meant
Multiplication by x acts pointwise, so on the grid it is X = diag(x₁ … xN) and ⟨x⟩ = Δx Σⱼ xⱼ|ψⱼ|². Nothing is approximated: X is the position operator restricted to the sample points. Its eigenvectors are the unit columns eⱼ, and once rescaled by 1/√Δx they are this space's stand-in for |x⟩ — which is a genuine gain, because the continuum's position eigenkets are δ-normalised distributions outside L2 while eⱼ is an honest vector in CN with an honest eigenvalue xⱼ. It is also a genuine loss. The grid has quietly declared the resolution of a position measurement to be Δx: structure in |ψ|² finer than 0.20 nm is simply not represented, and X² = diag(xⱼ²) inherits the same blindness, so ⟨x²⟩ is a four-term sum rather than an integral.
P is a stencil, and only the centred one is Hermitian
p̂ = −iħ d/dx has no pointwise action, so it must be built from neighbours. The centred difference gives (Pψ)ⱼ = −iħ(ψ_{j+1} − ψ_(j−1))/2Δx, so P_{j, j+1} = −iħ/2Δx and P_(j, j−1) = +iħ/2Δx. Check the adjoint: (P†)_{j+1, j} = (P_{j, j+1})* = +iħ/2Δx = P_{j+1, j}. Hermitian — and it is the centring that does it. Build P from the forward difference (ψ_{j+1} − ψⱼ)/Δx instead and the matrix carries entries above the diagonal and nothing below, so ‖P − P†‖ is of order ħ/Δx and grows as the grid is refined. On a real ψ that matrix returns a pure-imaginary ⟨p⟩, which is the fastest way to find the bug. The FFT route differs in kind: P = F†diag(ħkₙ)F with kₙ = 2πn/(NΔx) is exactly Hermitian and spectrally accurate, but it imposes periodic walls, so it solves a ring rather than a box.
Every expectation is a sandwich; P⋅P is not the kinetic matrix
With matrices in hand, every question in this unit is linear algebra: ⟨A⟩ = Δx ψ†Aψ, Var(A) = Δx ψ†A²ψ − (Δx ψ†Aψ)², and H|n⟩ = Eₙ|n⟩ is one call to numpy.linalg.eigh. One trap deserves naming. The product P@P is not the three-point kinetic matrix. Applying the centred difference twice reaches out to 2Δx, and its plane-wave symbol is (ħ sin(kΔx)/Δx)², while the three-point Laplacian gives (2ħ sin(kΔx/2)/Δx)². Expand both: the three-point form is low by (kΔx)²/12, P⋅P by (kΔx)²/3, four times worse. At Δx = 0.20 nm and k = 4.0 nm⁻¹, the exact ħ²k²/2m is 0.610 eV, the three-point matrix returns 0.578 eV and P⋅P returns 0.490 eV. Worse still, sin(kΔx) vanishes at kΔx = π, so P⋅P assigns zero kinetic energy to the alternating vector (+1, −1, +1, −1, …) — the odd-even decoupling that makes the three-point second derivative the right choice inside a Hamiltonian.
The wall and the spacing are two errors pointing opposite ways
Cutting the matrix off at N rows is the statement ψ₀ = ψ_{N+1} = 0: a Dirichlet wall you never wrote down, which raises the energy of any state with amplitude beyond the box. The spacing pushes the other way. The eigenvalues of (ħ²/2mΔx²)tridiag(−1,2,−1) are exactly (2ħ²/mΔx²)sin²(nπΔx/2L), always below ħ²kₙ²/2m by the factor 1 − (kₙΔx)²/12. For an electron in a 1.00 nm box with N = 4, the prefactor ħ²/2mΔx² is 0.952 eV, and the ground level comes out at 0.364 eV against the exact 0.376 eV — 3.2% low, matching the predicted −(πΔx/L)²/12 = −3.3%. The fourth level is 3.45 eV against 6.02 eV, 43% low. That is the general shape: the bottom of a grid spectrum is trustworthy and the top is not, so keep kΔx well under 1 for every level you intend to quote. Beyond the Nyquist cut |k| ≤ π/Δx — 15.7 nm⁻¹, or 9.40 eV of kinetic energy, at Δx = 0.20 nm — the grid cannot represent the state at all.
Finite dimensions hide the domain and cannot hold [x, p] = iħ
In CN, Hermitian and self-adjoint are one condition: a matrix equal to its conjugate transpose, no domain to state, no boundary term to check. That reads as a simplification and is really a concealment. The continuum question — which self-adjoint extension of −(ħ²/2m)d²/dx² you meant, whether −iħ d/dx has one at all on a half-line — has not been answered; the discretisation answered it for you when you chose a stencil and a truncation. The sharpest symptom is the canonical commutator. Since P is nonzero only on the first off-diagonals, [X, P]_{jk} = (xⱼ − xₖ)P_{jk} gives [X, P] = (iħ/2)(S₊ + S₋), a matrix with no diagonal at all, eigenvalues iħcos(nπ/(N+1)) and trace zero. It cannot equal iħI at any N, because Tr[A, B] = 0 for any two finite matrices while Tr(iħIN) = iħN. What does converge is the expectation on states the grid resolves: ⟨[X, P]⟩ = iħ[1 − Δx²⟨p²⟩/2ħ²].
Change one variable at a time
Make the relationship visible.
Hold k at 4 nm⁻¹ and drag h from 0.10 to 0.38 nm: the exact parabola never moves while both grid bands sag away from it, the three-point band by (kh)²/12 and the P⋅P band by (kh)²/3. Watch the faint P⋅P curve dive back towards zero as kh approaches π.
EXACT ħ²k²/2m0.610 eV
THREE-POINT LAPLACIAN0.578 eV
P⋅P, CENTRED STENCIL0.490 eV
NYQUIST k = π/h15.7 nm⁻¹
Live interpretationEXACT ħ²k²/2m: 0.610 eV. THREE-POINT LAPLACIAN: 0.578 eV. P⋅P, CENTRED STENCIL: 0.490 eV. NYQUIST k = π/h: 15.7 nm⁻¹
Catch the common trap
Explain before calculating.
On an N-point grid with X = diag(xⱼ) and P built from the centred difference, you form the matrix [X, P] and then refine the spacing Δx. What happens?
Choose an answer to test the model.
Practice & worked examples
Reason from the model, then test the result.
EasyAn electron sits in a 1.00 nm infinite box. Discretise it with N = 4 interior points: write down X and the kinetic matrix, then get the lowest grid eigenvalue and compare it with the exact ground state.
- N = 4 interior points on [0, L] means Δx = L/(N+1) = 1.00/5 = 0.20 nm, samples at 0.20, 0.40, 0.60, 0.80 nm, so X = diag(0.20, 0.40, 0.60, 0.80) nm. Truncating at four rows is the Dirichlet statement ψ₀ = ψ₅ = 0.
- Kinetic prefactor: ħ²/2mΔx² = (1.0546×10⁻³⁴)²/(2 × 9.109×10⁻³¹ × (2.0×10⁻¹⁰)²) = 1.526×10⁻¹⁹ J = 0.952 eV, so T = 0.952 eV × tridiag(−1, 2, −1).
- The eigenvalues of tridiag(−1,2,−1) are 2 − 2cos(nπ/5) = 4sin²(nπ/10). For n = 1: 4(0.30902)² = 0.38197, so E₁ = 0.952 × 0.38197 = 0.364 eV.
- Exact: E₁ = π²ħ²/2mL² = 0.376 eV. The grid is 3.2% low, matching the predicted −(k₁Δx)²/12 = −(π × 0.20/1.00)²/12 = −3.3%.
AnswerX = diag(0.20, 0.40, 0.60, 0.80) nm; E₁ = 0.364 eV on the grid against 0.376 eV exactly — 3.2% low, and always low, because the three-point Laplacian underestimates k².
MediumOn a grid of spacing Δx = 0.20 nm an electron state is close to a plane wave with k = 4.0 nm⁻¹. Find the kinetic energy returned by the exact operator, by the three-point Laplacian, and by P⋅P built from the centred difference. Rank them and say why.
- Set the scale: ħ²/2mΔx² = 0.952 eV, and kΔx = 4.0 nm⁻¹ × 0.20 nm = 0.80 rad, comfortably inside the Nyquist limit kΔx ≤ π.
- Exact: ħ²k²/2m = 0.952 eV × (kΔx)² = 0.952 × 0.64 = 0.610 eV.
- Three-point Laplacian: (ħ²/2mΔx²) · 4sin²(kΔx/2) = 0.952 × 4sin²(0.40) = 0.952 × 0.6066 = 0.578 eV, which is 5.2% low and matches −(kΔx)²/12 = −5.3%.
- P⋅P from the centred stencil: (ħ²/2mΔx²) · sin²(kΔx) = 0.952 × 0.5146 = 0.490 eV, which is 19.6% low and matches −(kΔx)²/3 — four times the leading error, because two centred differences reach out to 2Δx.
- The ranking is forced, not accidental: sin²θ = 4sin²(θ/2)cos²(θ/2) ≤ 4sin²(θ/2) ≤ θ², so P⋅P ≤ three-point ≤ exact at every k. At kΔx = π, sin(kΔx) = 0 and P⋅P gives the alternating vector zero kinetic energy, while the three-point matrix gives it 4 × 0.952 = 3.81 eV.
AnswerExact 0.610 eV; three-point 0.578 eV (5.2% low); P⋅P 0.490 eV (19.6% low). Use the three-point Laplacian for T, never P@P.
HardAn electron is sampled at spacing Δx = 0.20 nm in a minimum-uncertainty Gaussian of position width σ = 0.50 nm. (a) Evaluate ⟨[X, P]⟩ for the centred-difference matrices. (b) What spacing would cut the shortfall to 0.1%? (c) Show that no spacing removes it from the operator itself.
- Matrix elements: [X, P]_{jk} = (xⱼ − xₖ)P_{jk}, and P lives only on the first off-diagonals, so [X, P] = (iħ/2)(S₊ + S₋) — no diagonal, hence nothing like iħI.
- Expectation: ⟨[X, P]⟩ = iħ Δx Σⱼ ψⱼ*(ψ_{j+1} + ψ_(j−1))/2. Taylor-expand: (ψ_{j+1} + ψ_(j−1))/2 = ψⱼ + (Δx²/2)ψ″ⱼ + O(Δx⁴), and ∫ψ*ψ″ dx = −∫|ψ′|² dx = −⟨p²⟩/ħ², so ⟨[X, P]⟩ = iħ[1 − Δx²⟨p²⟩/2ħ²].
- For a minimum-uncertainty Gaussian, ⟨p²⟩ = ħ²/4σ², so the bracket is 1 − Δx²/8σ² = 1 − (0.20)²/(8 × 0.50²) = 1 − 0.0400/2.00 = 0.980. Hence ⟨[X, P]⟩ = 0.980 iħ, a 2.0% shortfall.
- For 0.1%: Δx²/8σ² = 1.0×10⁻³ gives Δx = σ√(8×10⁻³) = 0.50 × 0.0894 = 0.045 nm — 4.5 times finer, so 4.5 times the matrix dimension over the same box, for an error falling only as Δx².
- But the operator never converges. Tr[X, P] = Tr(XP) − Tr(PX) = 0 for any two N×N matrices, while Tr(iħIN) = iħN. Concretely the eigenvalues of [X, P] are iħcos(nπ/(N+1)): near +iħ on the smoothest grid vector, near −iħ on the most oscillatory one, averaging to zero.
Answer⟨[X, P]⟩ = 0.980 iħ at Δx = 0.20 nm; Δx ≈ 0.045 nm for a 0.1% shortfall; and [X, P] = iħI at no N, since Tr[X, P] = 0 while Tr(iħIN) = iħN.