Finite differences for stationary states

The potentials that can be solved exactly in the main text can be counted on one hand. But with a computer, the bound states of any one-dimensional potential can be solved in under a second — the method is to turn the differential equation into a matrix eigenvalue problem and hand it to a mature diagonalisation routine. Section 4.4 and the custom-potential solver of section 4.5 run on exactly this algorithm; here we give the full derivation, code and checks.

The idea: functions become vectors, operators become matrices

The stationary equation (natural units $\hbar=m=1$, used consistently in the code):

$$-\frac{1}{2}\psi''(x)+V(x)\psi(x)=E\psi(x)\tag{C.1.1}$$

Discretisation: chop the continuous interval into $N$ grid points $x_j$ with spacing $h$. The wavefunction $\psi(x)$ becomes an $N$-dimensional vector $(\psi_1,\dots,\psi_N)$ — the literal realisation of Chapter 3’s “a function is an infinite-dimensional vector”, except the dimension is now finite.

The remaining question: how do you express the second derivative $\psi''$ using grid values?

The three-point scheme and its error orderbasic~5 minexpand

Taylor-expand $\psi(x_j\pm h)$ about $x_j$:

$$\psi_{j\pm1}=\psi_j\pm h\psi_j'+\frac{h^2}{2}\psi_j''\pm\frac{h^3}{6}\psi_j'''+\frac{h^4}{24}\psi_j^{(4)}+\cdots\tag{C.1.2}$$

Add the two: the odd-order terms ($\psi'$, $\psi'''$) cancel in pairs,

$$\psi_{j+1}+\psi_{j-1}=2\psi_j+h^2\psi_j''+\frac{h^4}{12}\psi_j^{(4)}+\cdots\tag{C.1.3}$$

Rearranging gives the three-point difference scheme:

$$\psi_j''=\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{h^2}-\frac{h^2}{12}\psi_j^{(4)}+\cdots\tag{C.1.4}$$

The leading error term is $\propto h^2$: halve the spacing and the error drops to a quarter. This is called second-order accuracy, and in a moment we will verify that “factor of 4” directly by numerical experiment.

Substituting into the stationary equation, row $j$ becomes an algebraic equation connecting only neighbouring grid points:

$$-\frac{1}{2h^2}\psi_{j-1}+\left(\frac{1}{h^2}+V_j\right)\psi_j-\frac{1}{2h^2}\psi_{j+1}=E\psi_j\tag{C.1.5}$$

In matrix form, $H$ is a tridiagonal matrix: main diagonal $1/h^2+V_j$, both off-diagonals uniformly $-1/2h^2$. Boundary conditions: take a box large enough and demand $\psi_0=\psi_{N+1}=0$ (Dirichlet boundaries — equivalent to wrapping everything in an infinite well). As long as the box is much larger than the range over which the wavefunction decays, this artificial wall has no effect.

Complete code

Solve the harmonic oscillator $V=x^2/2$ and compare with the exact levels $E_n=n+\tfrac12$:

import numpy as np
from scipy.linalg import eigh_tridiagonal

# ---- Parameters (natural units hbar = m = 1) ----
L = 16.0          # solution interval [-L/2, L/2]
N = 2000          # number of interior grid points
x = np.linspace(-L/2, L/2, N + 2)[1:-1]   # drop the two boundary points
h = x[1] - x[0]

# ---- Potential: harmonic oscillator V = x^2 / 2 (change only this line for another potential) ----
V = 0.5 * x**2

# ---- Hamiltonian matrix: tridiagonal ----
diag = 1.0 / h**2 + V                     # main diagonal
off  = -0.5 / h**2 * np.ones(N - 1)       # off-diagonals

# ---- Diagonalise (only the lowest few eigenvalues are needed) ----
E, psi = eigh_tridiagonal(diag, off, select='i', select_range=(0, 5))

# ---- Normalisation: eigh returns vectors with sum(psi^2)=1; convert to integral normalisation ----
psi /= np.sqrt(h)

print("  n   numerical E_n   exact E_n   error")
for n in range(6):
    exact = n + 0.5
    print(f"  {n}   {E[n]:.8f}      {exact:.1f}         {abs(E[n]-exact):.2e}")

Output (measured):

  n   numerical E_n   exact E_n   error
  0   0.49999800      0.5         2.00e-06
  1   1.49999001      1.5         9.99e-06
  2   2.49997403      2.5         2.60e-05
  3   3.49995005      3.5         5.00e-05
  4   4.49991808      4.5         8.19e-05
  5   5.49987812      5.5         1.22e-04

Six significant figures, in under a second. The key lines, one by one:

• np.linspace(...)[1:-1]: generate the $N+2$ points including both ends, then throw the ends away. On the boundary $\psi=0$ is known — it should not be an unknown; the matrix then naturally excludes those points. That is the entire implementation of Dirichlet boundaries.
• eigh_tridiagonal: a diagonalisation routine that eats tridiagonal matrices exclusively, at far lower cost than brute-forcing the full matrix (select_range=(0, 5) asks for only the six lowest eigenpairs, indices 0 to 5 — faster still).
• psi /= np.sqrt(h): the numerical routine normalises by the vector inner product $\sum_j\psi_j^2=1$, while physics wants $\int\psi^2\dd x\approx h\sum_j\psi_j^2=1$ — a factor of $\sqrt h$ apart.
• Note that the error grows with the level: higher excited states oscillate faster, and the same grid spacing strains harder to resolve them.

Convergence check: seeing the h² law with your own eyes

Every numerical result deserves this step: double the grid density and watch by what factor the error shrinks.

import numpy as np
from scipy.linalg import eigh_tridiagonal

def ground_energy(N, L=16.0):
    x = np.linspace(-L/2, L/2, N + 2)[1:-1]
    h = x[1] - x[0]
    diag = 1.0 / h**2 + 0.5 * x**2
    off = -0.5 / h**2 * np.ones(N - 1)
    E, _ = eigh_tridiagonal(diag, off, select='i', select_range=(0, 0))
    return E[0], h

print("   N        h        |E0 - 0.5|    ratio")
prev = None
for N in [125, 250, 500, 1000, 2000]:
    E0, h = ground_energy(N)
    err = abs(E0 - 0.5)
    ratio = f"{prev/err:.2f}" if prev else "  --"
    print(f"{N:5d}   {h:.4f}   {err:.3e}   {ratio}")
    prev = err

Measured output:

   N        h        |E0 - 0.5|    ratio
  125   0.1270   5.044e-04     --
  250   0.0637   1.270e-04   3.97
  500   0.0319   3.187e-05   3.98
 1000   0.0160   7.984e-06   3.99
 2000   0.0080   1.998e-06   4.00

The ratio converges cleanly to 4 — exactly the $O(h^2)$ the derivation predicted. If your program yields a ratio that is not 4 (say 2, or something erratic), another error source is interfering — most often one of the traps below.

Try another potential

Swap out the single line V = 0.5 * x**2 in the code. A few worth trying:

Potential	Code	What to look for
Quartic oscillator $x^4$	V = x**4	No exact solution! Level spacing widens with $n$ (contrast the evenly spaced oscillator)
Double well	V = (x**2 - 4)**2 / 8	The lowest two levels are nearly degenerate — tunnelling splitting, connecting to 11.4 Instantons
Finite well	V = np.where(abs(x) < 2, -1.0, 0.0)	Finitely many bound states; compare section 2.8
Linear potential $\lvert x\rvert$	V = np.abs(x)	A toy model of quark confinement; levels given by zeros of the Airy function

∎Key formulas

• Three-point difference

  $$\psi_j''=\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{h^2}+O(h^2)$$

  Add the Taylor expansions; odd orders cancel automatically

• Discrete Hamiltonian

  $$H_{jj}=\frac{1}{h^2}+V_j,\qquad H_{j,j\pm1}=-\frac{1}{2h^2}$$

  Tridiagonal matrix; Dirichlet boundary = simply drop the boundary unknowns

• Convergence criterion

  $$\frac{\text{err}(h)}{\text{err}(h/2)}\to 4$$

  The fingerprint of a second-order method; a wrong ratio means another error source

?Self-check3 questions

1. 1.

   Shrink the grid spacing h to 1/3 of its value. The discretisation error of the three-point scheme becomes roughly:

   • 1/9
   • 1/3
   • 1/27
   • unchanged

   Submit
2. 2.

   Solving the oscillator numerically, you find the grid is already very fine, yet the ground-state energy still deviates from 0.5 and no longer improves with N. The most likely cause is:

   • The box length L is too small — the wavefunction is squeezed by the artificial wall
   • scipy's diagonalisation routine lacks precision
   • The harmonic oscillator has no exact solution anyway
   • The wavefunction was never normalised

   Submit
3. 3.

   Diagonalising the N = 2000 Hamiltonian yields 2000 eigenvalues. The trustworthy ones are:

   • The low-lying bound states whose wavelength is much larger than h and which decay to zero before the boundary
   • All of them — they are exact eigenvalues of the matrix
   • Only the ground state
   • The high-lying ones, since high-energy states care less about the details of the potential

   Submit