Finite-Difference Schrödinger Equation

A finite-difference Hamiltonian replaces spatial derivatives by local combinations of wavefunction samples. For a constant-mass one-dimensional problem, the centered second derivative produces a tridiagonal Hermitian eigenproblem whose low-energy spectrum converges quadratically under suitable smoothness and boundary assumptions. Solving that matrix accurately is only the first part of approximating the continuum problem.

Required background. Boundary conditions in computation specifies endpoint rows; array representations distinguishes samples from Euclidean coefficients.

The centered kinetic stencil

Let $\alpha=\hbar^2/(2m)$ and $H=-\alpha\partial_x^2+V(x)$ on $(a,b)$ with Dirichlet endpoints. Use $x_j=a+jh$, $j=1,\ldots,N$, and $h=(b-a)/(N+1)$. Taylor expansion gives

$$\frac{\psi(x+h)-2\psi(x)+\psi(x-h)}{h^2} =\psi''(x)+\frac{h^2}{12}\psi^{(4)}(x)+O(h^4).$$

Hence

$$H_{jj}=\frac{2\alpha}{h^2}+V(x_j), \qquad H_{j,j+1}=H_{j+1,j}=-\frac{\alpha}{h^2}.$$

The expansion assumes enough smoothness locally. Discontinuous potentials and singular endpoints can reduce the observed order or introduce alignment effects. Treat interfaces explicitly or test refinement with care; the formal stencil order alone does not determine the error of every observable.

Use a tridiagonal or sparse solver for large $N$. Dense diagonalization is useful as a small transparent reference calculation. A matrix-free implementation must use the same endpoint rule and should agree with the assembled matrix on test vectors.

Normalization and observables

If an eigensolver returns $c^\dagger c=1$, the continuum samples are approximately $\psi_j=c_j/\sqrt h$. Then $h\sum_j|\psi_j|^2=1$. For a multiplication observable $f(x)$,

$$\langle f\rangle\approx h\sum_j f(x_j)|\psi_j|^2 =\sum_j f(x_j)|c_j|^2.$$

Do not multiply the last expression by $h$ again. Derivative observables need a compatible discrete operator and boundary interpretation; pointwise differentiating a poorly resolved eigenvector can produce much larger errors than its energy suggests.

Exact discrete square-well benchmark

For $V=0$ and $L=b-a$, try $c_j=\sin(j\theta)$ with $\theta=n\pi/(N+1)$, $n=1,\ldots,N$. These samples vanish at both endpoints. The identity $c_{j+1}+c_{j-1}=2\cos\theta\,c_j$ makes the kinetic stencil a multiple of the same vector:

$$(Hc)_j=\frac{2\alpha}{h^2}(1-\cos\theta)c_j.$$

Thus the exact discrete eigenvalues are

$$E_n^{(N)}=\frac{4\alpha}{h^2} \sin^2\!\left(\frac{n\pi}{2(N+1)}\right).$$

For fixed $n$ and $h\to0$,

$$E_n^{(N)} =\alpha\left(\frac{n\pi}{L}\right)^2 -\frac{\alpha h^2}{12}\left(\frac{n\pi}{L}\right)^4+O(h^4).$$

The eigenvalues approach the continuum values from below. This finite-difference construction is not automatically a conforming Rayleigh–Ritz approximation with upper bounds. Its relative dispersion error grows approximately as $(n\pi h/L)^2/12$, explaining why a grid that resolves the ground state can fail badly near its highest modes.

A short NumPy check in units $\alpha=L=1$ is:

import numpy as np
N = 80
h = 1.0 / (N + 1)
H = (np.diag(np.full(N, 2.0))
     - np.diag(np.ones(N - 1), 1)
     - np.diag(np.ones(N - 1), -1)) / h**2
E, C = np.linalg.eigh(H)
E_discrete = 4.0 / h**2 * np.sin(np.pi * h / 2.0)**2
print(E[0], abs(E[0] - E_discrete), abs(E[0] - np.pi**2))
print(np.linalg.norm(H @ C[:, 0] - E[0] * C[:, 0]))

The discrete target is about $9.868367239$; its continuum error is about $0.001237162$. The matrix-solver discrepancy and residual should be much smaller in binary64. The program separates these errors instead of mistaking a small residual for a continuum accuracy certificate.

Potentials and the independent box limit

For a harmonic oscillator in dimensionless units, insert $V_j=x_j^2/2$ and use $\alpha=1/2$. On a finite symmetric box, first refine $h$ and then enlarge the box at fixed sufficiently small $h$. The target ground energy is $1/2$; tail truncation and stencil error can have opposite signs and accidentally cancel.

A finite square well provides a different challenge: moving the mesh relative to a discontinuous edge can change the sampled potential abruptly. Align the interface deliberately, integrate its effect in another discretization, or document an interface-aware refinement sequence. A smooth-potential order test does not automatically validate this case.

Diagnostics and limitations

Check Hermiticity, endpoint behavior, normalization, node count, and residuals. Compare low states with analytic benchmarks and run an observable-specific convergence study. Increase the box independently of the number of interior points.

Higher-order stencils require matching higher-order endpoint treatment. They can also widen the matrix bandwidth. For smooth problems, matrix Numerov or spectral methods may achieve a target accuracy with fewer points; near interfaces, finite elements offer useful flexibility.

Exercises

At fixed mode index $n$, approximately how does halving $h$ change the leading eigenvalue error? What if $n$ is doubled while $h$ stays fixed and the mode remains well resolved?

Solution

Halving $h$ reduces the leading absolute error by four. Doubling $n$ increases the leading absolute error by sixteen and the relative error by four. The expansion requires $n\pi h/L\ll1$; it is not valid near the grid cutoff.

References

• LeVeque, R. J., Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007, author resources — stencil consistency, stability, and convergence.
• Trefethen, L. N., and D. Bau III, Numerical Linear Algebra, SIAM, 1997 — structured eigenproblems and residual checks.