Finite Difference Methods

Finite difference methods approximate derivatives by differences of nearby grid values. In quantum mechanics, they turn differential Hamiltonians into sparse matrices that can be diagonalized, propagated in time by Time-Stepping Methods, or used in boundary-value calculations. The full PDE workflow is organized in PDE Solvers.

The core idea is local: if a wavefunction is smooth on the scale of the grid spacing $h$, Taylor expansions relate derivatives at $x_j$ to samples at neighboring points. The numerical danger is equally local: if the grid does not resolve the wavefunction or the boundary condition is implemented incorrectly, the resulting matrix solves a different problem.

Why Quantum Mechanics Needs Them

Finite differences appear in:

• one-dimensional bound-state calculations;
• finite-box approximations to continuum problems;
• tunneling through smooth barriers;
• radial Schrödinger equations after coordinate reduction;
• time-dependent wave-packet propagation on grids;
• finite-volume or finite-element methods in more advanced settings;
• benchmark calculations for validating more elaborate codes.

They are not always the most accurate method for smooth periodic problems, where Spectral Methods may converge faster. But finite differences are transparent, sparse, local, and easy to inspect.

Uniform Grid

On a uniform one-dimensional grid,

$$x_j=x_0+jh, \qquad \psi_j\approx\psi(x_j).$$

Taylor expansions give

$$\psi(x_j+h) = \psi(x_j) + h\psi'(x_j) + \frac{h^2}{2}\psi''(x_j) + \frac{h^3}{6}\psi'''(x_j) + O(h^4),$$

and

$$\psi(x_j-h) = \psi(x_j) - h\psi'(x_j) + \frac{h^2}{2}\psi''(x_j) - \frac{h^3}{6}\psi'''(x_j) + O(h^4).$$

Adding or subtracting these expansions produces derivative formulas.

First Derivatives

The centered first derivative is

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

One-sided formulas are useful near boundaries. The forward difference

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

is only first-order accurate. Higher-order one-sided formulas exist, but boundary accuracy and boundary conditions should be handled deliberately rather than patched after the fact.

In quantum mechanics, first derivatives appear in momentum operators, current densities, radial transformations, and gauge-coupled Hamiltonians. A careless first-derivative discretization can break Hermiticity or time-reversal properties.

Second Derivatives

The centered second derivative is

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

For the one-dimensional Hamiltonian

$$H = - \frac{\hbar^2}{2m} \frac{d^2}{dx^2} + V(x),$$

this gives the grid action

$$(H_h\psi)_j = - \frac{\hbar^2}{2m} \frac{ \psi_{j+1}-2\psi_j+\psi_{j-1} }{h^2} + V_j\psi_j,$$

where $V_j=V(x_j)$.

This is the standard starting point for finite-difference bound-state calculations.

Laplacian Matrix

For $N$ interior grid points with Dirichlet boundary conditions at the endpoints, the second-difference matrix is

$$D_2 = \frac{1}{h^2} \begin{pmatrix} -2 & 1 & 0 & \cdots & 0\\ 1 & -2 & 1 & \cdots & 0\\ 0 & 1 & -2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & 1\\ 0 & 0 & 0 & 1 & -2 \end{pmatrix}.$$

The Hamiltonian matrix is

$$H_h = - \frac{\hbar^2}{2m}D_2 + \operatorname{diag}(V_1,\dots,V_N).$$

For real $V_j$, this matrix is real symmetric. With the usual uniform-grid inner product, it is Hermitian and can be passed to a Hermitian eigensolver.

The resulting matrix is sparse: each row couples only neighboring grid points. This sparsity is one reason finite differences scale to much larger grids than dense methods. The storage and matrix-vector-product language is developed in Sparse Matrices.

Boundary Conditions

Boundary conditions change the first and last rows of the matrix.

For Dirichlet boundaries, one often stores only interior points and sets the exterior boundary values to zero. For periodic boundaries, the first and last interior points are coupled:

$$\psi_{-1}=\psi_{N-1}, \qquad \psi_N=\psi_0.$$

For Neumann boundaries, one approximates a derivative condition such as $\psi'(x_{\min})=0$. A simple ghost-point implementation sets

$$\psi_{-1}=\psi_1,$$

which enforces the centered derivative to vanish at the boundary point.

Boundary rows are common sources of errors because they are not the same as the interior stencil. Always state how the boundary condition was implemented.

Higher-Order Stencils

The basic centered second derivative is second-order accurate. A fourth-order version is

$$\psi''(x_j) \approx \frac{ -\psi_{j+2} + 16\psi_{j+1} - 30\psi_j + 16\psi_{j-1} - \psi_{j-2} }{12h^2}.$$

Higher-order stencils can improve accuracy for smooth wavefunctions, but they use wider neighborhoods. Wider stencils complicate boundaries, increase matrix bandwidth, and may behave poorly near nonsmooth potentials or discontinuities.

Order is not the only criterion. A lower-order method with correct boundary treatment can outperform a higher-order stencil with inconsistent boundary rows.

Multidimensional Laplacians

On a rectangular two-dimensional grid with spacings $h_x$ and $h_y$, a common Laplacian approximation is

$$\nabla_h^2 = D_x\otimes I_y + I_x\otimes D_y,$$

where $D_x$ and $D_y$ are one-dimensional second-difference matrices and $\otimes$ is the matrix Kronecker product.

This tensor-product structure is valuable. It gives sparse Hamiltonians and makes separable benchmark problems easy to check. On irregular domains or curvilinear coordinates, additional metric factors and boundary geometry must be handled carefully.

Accuracy and Convergence

For a second-order stencil, the local derivative error is $O(h^2)$ when the wavefunction is sufficiently smooth. For eigenvalues, the observed convergence can depend on the state, boundary conditions, potential smoothness, and box size.

A standard grid-refinement test compares results at $h$, $h/2$, and $h/4$. If the leading error is proportional to $h^2$, then the error should shrink by about a factor of $4$ when $h$ is halved.

This is only meaningful after domain truncation error is under control. For bound states in a finite box, increase the box size and refine the grid separately.

For a systematic refinement workflow, see Convergence Tests.

Roundoff and Conditioning

Derivative matrices grow as $h$ shrinks. The second-difference matrix has entries of size $1/h^2$, so roundoff and conditioning can become more visible on very fine grids.

For a second derivative, a schematic error balance is

$$\text{error} \sim C_1h^2 + C_2\frac{u}{h^2},$$

where $u$ is the unit roundoff. The first term decreases with refinement; the second can increase. See Floating-Point Arithmetic and Conditioning and Stability for the numerical background.

Physical Validation

Finite-difference Hamiltonians should be tested on known cases:

• infinite square well energies;
• harmonic oscillator low-lying levels;
• free-particle dispersion on a periodic grid;
• symmetry of even and odd states in symmetric potentials;
• normalization and orthogonality under the discrete inner product.

For the infinite well, the exact continuum energies are proportional to $n^2$. A finite-difference grid should reproduce low-$n$ levels increasingly well under refinement, while high-$n$ levels near the grid cutoff converge poorly.

Common Mistakes

• Using the interior stencil at a boundary without implementing the boundary condition.
• Forgetting the factor $1/h^2$ in the Laplacian.
• Assuming high-energy eigenvalues converge as quickly as low-energy eigenvalues.
• Treating a finite-box spectrum as a continuum spectrum.
• Refining $h$ while ignoring roundoff or box-size effects.
• Breaking Hermiticity with an asymmetric derivative stencil.
• Comparing grid wavefunctions without using the correct discrete inner product.
• Using a high-order stencil near a discontinuity and expecting high-order convergence.

Cross-Links

• Discretization
• Spectral Methods
• Floating-Point Arithmetic
• Conditioning and Stability
• Matrix Diagonalization
• Sparse Matrices
• Time-Stepping Methods
• PDE Solvers
• Convergence Tests
• Boundary Conditions
• Eigenvalue Problems
• Matrices as Linear Maps
• Infinite Square Well
• Quantum Harmonic Oscillator
• Imaginary-Time Projection Notebook

References

• R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
• T. Pang, An Introduction to Computational Physics, 2nd ed., Cambridge University Press, 2006.
• W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007.
• G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd ed., Oxford University Press, 1985.

Exercises

1. Derive the centered second-derivative formula by adding Taylor expansions at $x+h$ and $x-h$.

Solution

Taylor expansion gives

$$\psi(x+h) = \psi + h\psi' + \frac{h^2}{2}\psi'' + \frac{h^3}{6}\psi''' + \frac{h^4}{24}\psi^{(4)} + O(h^5),$$

and

$$\psi(x-h) = \psi - h\psi' + \frac{h^2}{2}\psi'' - \frac{h^3}{6}\psi''' + \frac{h^4}{24}\psi^{(4)} + O(h^5).$$

Adding and solving for $\psi''$ gives

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

2. Write the second-difference matrix for three interior points with Dirichlet boundary conditions.

Solution

With three interior points, the matrix is

$$D_2 = \frac{1}{h^2} \begin{pmatrix} -2 & 1 & 0\\ 1 & -2 & 1\\ 0 & 1 & -2 \end{pmatrix}.$$

The missing boundary values are fixed to zero by the Dirichlet condition.

3. For real $V_j$, explain why the standard finite-difference Hamiltonian with Dirichlet boundaries is Hermitian.

Solution

The second-difference matrix is real symmetric, so $D_2^\dagger=D_2$. The potential matrix $\operatorname{diag}(V_j)$ is also real diagonal, hence Hermitian. A real linear combination of Hermitian matrices is Hermitian, so

$$H_h = - \frac{\hbar^2}{2m}D_2 + \operatorname{diag}(V_j)$$

is Hermitian.

4. A second-order finite-difference eigenvalue error is dominated by $Ch^2$. If halving $h$ does not reduce the error by about a factor of $4$, name two possible explanations.

Solution

One possibility is that another error source dominates, such as finite-box error, roundoff, or eigensolver tolerance. Another is that the assumptions behind second-order convergence are not met, for example because the potential or wavefunction is not sufficiently smooth or because the boundary stencil is only first-order accurate.

5. How does a periodic boundary condition modify the one-dimensional second-difference matrix?

Solution

The first and last grid points become neighbors. In addition to the usual tridiagonal entries, the matrix has corner entries connecting the first row to the last column and the last row to the first column. These entries implement

$$\psi_{-1}=\psi_{N-1}, \qquad \psi_N=\psi_0.$$