PDE Solvers

PDE solvers turn partial differential equations into finite numerical problems by combining a spatial discretization, boundary conditions, time integration or eigenvalue solution, and validation tests. In quantum mechanics the central examples are the time-dependent and time-independent Schrödinger equations in coordinate space.

This page is the numerical workflow home for PDEs. It does not replace Partial Differential Equations as the mathematical background, Finite Difference Methods as the stencil page, Spectral Methods as the global-basis page, or Time-Stepping Methods as the real-time dynamics page. It explains how those tools fit together in a complete computation.

Two Quantum PDE Tasks

The time-dependent Schrödinger equation is an evolution PDE:

$$i\hbar\frac{\partial\psi}{\partial t} = - \frac{\hbar^2}{2m}\nabla^2\psi + V(\mathbf r,t)\psi.$$

Given initial data and boundary conditions, a PDE solver must approximate $\psi(\mathbf r,t)$ over time.

The time-independent Schrödinger equation is a spatial eigenvalue PDE:

$$- \frac{\hbar^2}{2m}\nabla^2\varphi + V(\mathbf r)\varphi = E\varphi.$$

Given a domain, boundary conditions, and potential, a PDE solver must approximate selected eigenvalues and eigenfunctions.

These tasks are related but not identical. Time evolution cares about stability, phase accuracy, unitarity, and long-time drift. Stationary eigenproblems care about Hermiticity, spectral pollution, residuals, normalization, and convergence under spatial refinement.

Spatial Discretization First

A numerical PDE calculation starts by replacing functions on a continuum domain by finite data:

• grid samples on a line, rectangle, box, or mesh;
• basis coefficients in a finite expansion;
• finite-element degrees of freedom;
• spectral or pseudospectral coefficients;
• boundary unknowns in specialized integral-equation methods.

For a grid representation, a wavefunction becomes a vector

$$\psi_h = (\psi_1,\dots,\psi_N)^T.$$

The Laplacian becomes a matrix $L_h$, the potential becomes a multiplication matrix $V_h$, and the Hamiltonian becomes

$$H_h = - \frac{\hbar^2}{2m}L_h + V_h.$$

This finite matrix is the actual object passed to eigensolvers or time integrators. Its Hermiticity, sparsity, and boundary rows determine the numerical physics.

Method of Lines

For an evolution PDE, the most common workflow is the method of lines:

1. Discretize space.
2. Keep time continuous.
3. Solve the resulting large system of ODEs.

For the Schrödinger equation this gives

$$i\hbar \frac{d\psi_h}{dt} = H_h(t)\psi_h(t),$$

or

$$\frac{d\psi_h}{dt} = - \frac{i}{\hbar} H_h(t)\psi_h(t).$$

The ODE may be large, sparse, oscillatory, and complex-valued. Generic ODE Solvers can be useful, but closed-system quantum dynamics has extra structure: for Hermitian $H_h$, exact evolution is unitary. That is why Time-Stepping Methods and Matrix Exponentials Numerically are central.

Finite-Difference PDE Example

On a rectangular two-dimensional grid, a standard second-order Laplacian is

$$\begin{aligned} (\nabla_h^2\psi)_{j,k} &= \frac{\psi_{j+1,k}-2\psi_{j,k}+\psi_{j-1,k}}{h_x^2}\\ &\quad+ \frac{\psi_{j,k+1}-2\psi_{j,k}+\psi_{j,k-1}}{h_y^2}. \end{aligned}$$

With Dirichlet boundaries, the unknowns are often interior grid values only. The Hamiltonian matrix is sparse because each grid point couples only to nearby grid points.

For a stationary problem, one solves

$$H_h\varphi_h = E_h\varphi_h.$$

For a time-dependent problem, one solves

$$i\hbar \frac{d\psi_h}{dt} = H_h\psi_h.$$

The same spatial matrix can therefore support both eigenvalue and time-evolution calculations.

Spectral and Pseudospectral PDE Solvers

Spectral PDE solvers approximate the wavefunction using global basis functions. In a Fourier pseudospectral method:

1. Store $\psi_j$ on a periodic grid.
2. Use an FFT to obtain Fourier coefficients.
3. Multiply by $-k^2$ to apply the Laplacian.
4. Transform back to position space for potential multiplication.

This is efficient for smooth periodic or large-box wave packets. It also supports split-operator evolution, where kinetic-energy phases are applied in momentum space and potential phases in position space.

The method can fail quietly if periodicity, smoothness, or resolution assumptions are violated. Endpoint mismatch, rough potentials, and aliasing can all contaminate a calculation before the time integrator has a chance to be the limiting error.

Finite Elements and General Domains

Finite-element methods divide the domain into simple cells and approximate the wavefunction by local basis functions. They are useful for irregular geometry, nonuniform resolution, and complicated boundary conditions.

A finite-element stationary problem often has a generalized eigenvalue form

$$Hc = ESc,$$

where $S$ is the mass or overlap matrix. The inner product is therefore not the ordinary Euclidean dot product unless $S=I$. Norms, residuals, and orthogonality checks must use the correct discrete inner product.

Finite elements are not developed in detail here, but the same principles apply: preserve the operator domain, assemble Hermitian matrices when the continuum operator is self-adjoint, and refine the mesh to test convergence.

Time-Stepping Choices

Once the spatial Hamiltonian is finite, time evolution can be advanced in several ways.

For a time-independent Hamiltonian, exact finite-dimensional evolution is

$$\psi_h(t+\Delta t) = \exp\left( - \frac{i}{\hbar} H_h\Delta t \right)\psi_h(t).$$

For large systems, one usually computes the action of this exponential rather than the full dense matrix.

A Crank–Nicolson step is

$$\left( I+ \frac{i\Delta t}{2\hbar}H_h \right)\psi_h^{n+1} = \left( I- \frac{i\Delta t}{2\hbar}H_h \right)\psi_h^n.$$

For Hermitian $H_h$ and exact solves, this step is norm-preserving. The cost is a linear solve at every step.

For FFT-friendly Hamiltonians with $H=T+V$, a second-order split-operator step has the schematic form

$$\psi^{n+1} \approx e^{-iV\Delta t/(2\hbar)} e^{-iT\Delta t/\hbar} e^{-iV\Delta t/(2\hbar)} \psi^n.$$

This is efficient when $T$ is diagonal in momentum space and $V$ is diagonal in position space. Its error comes from noncommutativity of $T$ and $V$, time dependence, finite grid resolution, and FFT conventions.

Stability

Stability asks whether small numerical errors are amplified. For constant-coefficient grid problems, von Neumann analysis inserts a Fourier mode

$$\psi_j^n = G^n e^{ij\theta}$$

and studies the amplification factor $G(\theta)$.

For a mode that should remain bounded, a basic stability requirement is

$$\lvert G(\theta)\rvert \le 1$$

over the represented range of $\theta$. For Schrödinger evolution, boundedness is not enough: a good method should also approximate phases and preserve norm or inner products when the physics requires unitary evolution.

The finite Hamiltonian’s spectral radius often controls the allowed time step for explicit methods. Refining a grid increases the largest represented kinetic energy, so a stable time step may have to shrink even if the physical wave packet has low energy.

Consistency, Stability, and Convergence

A PDE method is consistent when its finite equations approach the continuum PDE as the grid spacing, time step, or basis truncation is refined. It is stable when perturbations do not grow in an uncontrolled way. It is convergent when the numerical solution approaches the true solution in the chosen norm.

For many linear well-posed initial-value problems, consistency plus stability implies convergence. This is the practical content behind the Lax equivalence theorem. The slogan is useful, but the assumptions matter: nonlinearities, boundaries, singular potentials, unbounded domains, and wrong norms can invalidate naive reasoning.

For a quantum grid calculation, a typical target estimate has the form

$$\lVert \psi_{h,\Delta t}(t) - \psi(t) \rVert \le C \left( h^p+\Delta t^q \right),$$

under smoothness, boundary, and stability assumptions. The exponents $p$ and $q$ must be observed by refinement, not merely quoted from the method description.

Boundary Conditions and Domains

Boundary conditions are part of the operator. A PDE solver must encode them in the finite problem:

• Dirichlet rows for hard walls;
• periodic wraparound couplings for rings or periodic boxes;
• Neumann or Robin derivative conditions;
• matching conditions at interfaces;
• absorbing layers or exterior complex scaling for outgoing waves;
• regularity conditions at coordinate singularities.

Boundary mistakes are especially dangerous because the interior stencil or basis can look correct. Always inspect the boundary implementation, not only the formula for the Laplacian.

Absorbing boundaries and complex potentials are useful numerical devices, but they change the Hamiltonian from Hermitian to non-Hermitian. Norm loss then represents designed absorption, not closed-system probability conservation.

Sparse Linear Algebra

Discretized PDEs in more than one dimension quickly become large. A grid with $N_xN_yN_z$ points has that many wavefunction unknowns, and dense matrices are usually impossible.

Local finite-difference and finite-element discretizations typically produce sparse matrices. Useful operations include:

• sparse matrix-vector products;
• sparse direct or iterative linear solves;
• Krylov exponential actions;
• selected sparse eigenvalue computations;
• preconditioned solves for implicit steps.

The storage and algorithmic background lives in Sparse Matrices and Sparse Eigensolvers.

Validation Workflow

Trustworthy PDE computations are built from checks:

• compare against a separable or exactly solvable case;
• refine the spatial grid, mesh, or basis size;
• refine the time step independently from the spatial resolution;
• monitor norm, energy, symmetry quantum numbers, and boundary flux;
• check eigenvalue residuals and orthogonality;
• compare finite-difference and spectral discretizations when possible;
• verify that the finite domain is large enough for localized states;
• inspect high-energy or high-wavenumber modes near the cutoff;
• report boundary conditions, grid conventions, tolerances, and units.

The point is not to make every computation expensive. The point is to identify which error source is currently limiting the result.

The general vocabulary for separating truncation, roundoff, algebraic solver, and statistical uncertainty is Error Estimates. The step-by-step refinement workflow is Convergence Tests.

Choosing a PDE Solver

Situation	Common first choice	Main caution
smooth periodic wave packet	Fourier pseudospectral plus FFT	aliasing and periodic wraparound
hard-wall box on a simple grid	finite differences	boundary rows and high-energy modes
irregular domain	finite elements	overlap matrix and mesh convergence
low-lying bound states	sparse eigensolver on $H_h$	residuals, box size, and spectral pollution
closed real-time dynamics	Crank–Nicolson or exponential action	unitarity and phase error
scattering with outgoing waves	absorbing layer or complex scaling	interpretation of norm loss
imaginary-time relaxation	diffusion-like time stepping	stiffness and normalization practices

The best method is the one whose failure modes you can test for the observable you care about.

Common Mistakes

• Treating the interior finite-difference stencil as the whole PDE solver.
• Refining the time step while leaving the spatial grid unconverged.
• Reporting a finite-box spectrum as if it were the continuum spectrum.
• Forgetting that spectral methods solve the periodic extension unless told otherwise.
• Using a non-Hermitian boundary or derivative implementation while expecting unitary evolution.
• Ignoring the overlap matrix in finite-element or nonorthogonal-basis calculations.
• Judging stability only from the physical low-energy modes and ignoring cutoff modes.
• Calling norm loss physical when it was created by an unintended non-Hermitian discretization.
• Comparing wavefunctions across grids without using a consistent norm and phase convention.

Cross-Links

• Partial Differential Equations
• Time-Dependent Schrödinger Equation
• Time-Independent Schrödinger Equation
• Discretization
• Finite Difference Methods
• Spectral Methods
• Fast Fourier Transform
• Time-Stepping Methods
• Matrix Exponentials Numerically
• Sparse Matrices
• Sparse Eigensolvers
• Conditioning and Stability
• Error Estimates
• Convergence Tests
• Benchmark Problems
• ODE Solvers

References

• R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007.
• J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed., SIAM, 2004.
• L. N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000.
• S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Springer, 2008.
• R. Kosloff, “Time-dependent quantum-mechanical methods for molecular dynamics”, Journal of Physical Chemistry 92, 2087-2100, 1988.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
• D. J. Tannor, Introduction to Quantum Mechanics: A Time-Dependent Perspective, University Science Books, 2007.

Exercises

1. Write the five-point finite-difference Laplacian on a two-dimensional uniform grid.

Solution

For spacings $h_x$ and $h_y$,

$$\begin{aligned} (\nabla_h^2\psi)_{j,k} &= \frac{\psi_{j+1,k}-2\psi_{j,k}+\psi_{j-1,k}}{h_x^2}\\ &\quad+ \frac{\psi_{j,k+1}-2\psi_{j,k}+\psi_{j,k-1}}{h_y^2}. \end{aligned}$$

Each grid point couples to its four nearest neighbors and itself, except where boundary conditions modify the formula.

2. Show how the spatially discretized Schrödinger equation becomes an ODE system.

Solution

After spatial discretization, the Hamiltonian becomes a finite matrix $H_h$. The wavefunction samples or coefficients form a vector $\psi_h(t)$. The PDE becomes

$$i\hbar \frac{d\psi_h}{dt} = H_h\psi_h.$$

Solving for the derivative gives

$$\frac{d\psi_h}{dt} = - \frac{i}{\hbar} H_h\psi_h,$$

which is a finite-dimensional linear ODE system.

3. For a Hermitian finite Hamiltonian, why is explicit Euler a poor real-time Schrödinger step?

Solution

For an energy eigenmode with $H_h\phi=E\phi$, the exact amplification factor over one step is $e^{-iE\Delta t/\hbar}$, which has modulus $1$. Explicit Euler gives

$$g = 1 - \frac{iE\Delta t}{\hbar}.$$

For real $E$,

$$\lvert g\rvert^2 = 1 + \left( \frac{E\Delta t}{\hbar} \right)^2 \gt 1$$

unless $E=0$. The method artificially grows the norm.

4. Show that the Crank–Nicolson matrix step is unitary for Hermitian $H_h$ in exact arithmetic.

Solution

Let

$$A = \frac{\Delta t}{2\hbar}H_h.$$

The step matrix is

$$U_{CN} = (I+iA)^{-1}(I-iA).$$

If $H_h$ is Hermitian, then $A$ is Hermitian. The factors $I+iA$ and $I-iA$ are adjoints of each other and commute because both are polynomials in $A$. Therefore

$$U_{CN}^\dagger U_{CN} = I.$$

The conclusion assumes exact linear solves; in floating-point arithmetic, solve tolerance also matters.

5. A second-order spatial method and a second-order time method are both used. If halving $\Delta t$ changes the answer but halving $h$ does not, what should be refined next?

Solution

The time-step error is currently more visible than the spatial error, so refine $\Delta t$ further while keeping the already-tested spatial grid fixed. After the time-step dependence is reduced, repeat a spatial refinement check because error sources can interact.