Numerical Mathematics

A trustworthy numerical result is a mathematical argument supported by computation. It identifies the intended continuum or finite-dimensional problem, explains how that problem was represented, chooses an algorithm appropriate to its structure, quantifies error, and survives independent checks. A solver returning a number is only the beginning of that argument.

This chapter develops method-level foundations. Software interfaces, implementation environments, and project workflows are outside its canonical scope; Math Needed for Computational QM connects these tools to computation-facing study routes.

The numerical trust chain

Most quantum calculations pass through several distinct problems:

$$\begin{gathered} \text{physical model} \longrightarrow \text{mathematical problem} \longrightarrow \text{finite representation}\\ \longrightarrow \text{algebraic algorithm} \longrightarrow \text{floating-point result} \longrightarrow \text{reported observable}. \end{gathered}$$

Each arrow introduces different assumptions and failure modes.

Stage	Typical choice	Essential evidence
mathematical problem	operator domain, boundary data, initial state	well-posed equations and physical units
finite representation	grid, basis, box, cutoff, quadrature	refinement and boundary checks
algorithm	eigensolver, linear solve, time step, transform	stability and residual diagnostics
arithmetic	precision, scaling, summation order	roundoff and conditioning analysis
reported quantity	energy, norm, phase, rate, expectation value	propagated uncertainty and benchmark comparison

Validation should target the observable being claimed. A small eigenpair residual does not by itself prove that the continuum eigenvalue is converged; it only says that the computed vector nearly solves the finite matrix problem.

Arithmetic, conditioning, and stability

Standard floating-point analysis models a basic rounded operation by

$$\operatorname{fl}(x\mathbin{\circ}y) = (x\mathbin{\circ}y)(1+\delta), \qquad \lvert\delta\rvert\le u,$$

when the exact result is in the normal range and no exceptional case intervenes. Here $u$ is the unit roundoff. Repeated operations can amplify these perturbations, especially through cancellation, poor scaling, or ill-conditioned transformations.

Conditioning belongs to the mathematical problem. For an invertible matrix,

$$\kappa(A) = \lVert A\rVert \lVert A^{-1}\rVert$$

measures worst-case sensitivity in the chosen norm. Stability belongs to the algorithm. A backward-stable method returns the exact answer to a nearby problem. Their roles combine schematically as

$$\text{forward error} \lesssim \text{condition number} \times \text{backward error}.$$

Floating-Point Arithmetic owns representation, roundoff, cancellation, scaling, and reproducibility. Conditioning and Stability separates sensitive problems from unstable algorithms and explains why residuals must be interpreted in scale-aware norms.

Discretizing continuum quantum mechanics

A discretization replaces an infinite-dimensional problem by finite data. That replacement includes more than a grid spacing or basis size:

• the spatial domain and its truncation;
• boundary conditions;
• grid points or basis functions;
• the discrete inner product and quadrature weights;
• operator matrices;
• ultraviolet and infrared resolution;
• symmetry and Hermiticity properties.

For example, the centered second derivative is

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

This interior order does not guarantee second-order convergence if boundary rows are lower order or the solution lacks the required smoothness. A nontrivial quadrature matrix $W$ may also change the discrete adjoint condition from $H^\dagger=H$ to

$$H^\dagger W = WH.$$

Discretization owns the continuum-to-finite transition. Finite Difference Methods develops local stencils and Laplacian matrices. Spectral Methods develops global basis and pseudospectral approximations, whose rapid convergence depends on smoothness and boundary matching. Numerical Quadrature supplies the weighted sums used for normalization, overlaps, expectation values, and matrix elements.

Dense and sparse spectral problems

After discretization, a stationary problem often becomes

$$Hv=\lambda v$$

or, in a nonorthogonal basis,

$$Hv=\lambda Sv.$$

Dense Hermitian diagonalization is appropriate when the matrix fits comfortably in memory and many eigenpairs are needed. Large local Hamiltonians are often sparse: only a small fraction of matrix entries are nonzero, and the important primitive is then the matrix-vector product rather than explicit factorization.

For a normalized approximate eigenvector $\tilde v$, the residual

$$r = H\tilde v - \tilde\lambda\tilde v$$

tests the finite algebraic problem. A useful scale-aware form is

$$\eta = \frac{\lVert r\rVert} {\lVert H\rVert\lVert\tilde v\rVert +\lvert\tilde\lambda\rvert\lVert\tilde v\rVert}.$$

Residuals should be accompanied by orthogonality checks, symmetry labels, and refinement in the underlying grid or basis. Near degeneracy, the invariant subspace is usually better conditioned than an individual basis of eigenvectors.

Matrix Diagonalization covers dense Hermitian eigensolvers and eigenpair validation. Sparse Matrices covers storage, matrix-free actions, and weighted Hermiticity. Sparse Eigensolvers covers Lanczos and Arnoldi methods, targeting, restarts, and convergence diagnostics.

Time evolution and differential equations

Spatial discretization turns the time-dependent Schrödinger equation into

$$i\hbar\frac{d\psi}{dt} = H(t)\psi.$$

For a time-independent Hermitian Hamiltonian, exact evolution over one step is

$$\psi_{n+1} = \exp\left( - \frac{i\Delta t}{\hbar}H \right) \psi_n.$$

Different methods approximate this structure in different ways. Crank–Nicolson uses

$$\left( I+\frac{i\Delta t}{2\hbar}H \right)\psi_{n+1} = \left( I-\frac{i\Delta t}{2\hbar}H \right)\psi_n$$

and is unitary for a time-independent Hermitian $H$ in exact arithmetic. Exponential-action methods approximate $e^A v$ without necessarily forming $e^A$. Split-operator methods alternate kinetic and potential evolution, often using a fast Fourier transform to move between position and momentum grids.

General ODE and PDE solvers add concerns such as stiffness, adaptive local error control, boundary-value shooting, multidimensional domains, and sparse linear solves. Norm conservation alone is not enough: a propagator can preserve norm while accumulating unacceptable phase or observable error.

Time-Stepping Methods compares explicit, implicit, and structure-aware propagators. Matrix Exponentials Numerically covers diagonalization, scaling and squaring, and Krylov exponential actions. Fast Fourier Transform fixes discrete transform and momentum-grid conventions. ODE Solvers and PDE Solvers own the broader initial-value, boundary-value, stability, and convergence frameworks.

Error budgets and convergence evidence

A numerical error budget may include:

• model or approximation error;
• finite-domain error;
• spatial discretization or basis-truncation error;
• time-step error;
• algebraic solver error;
• quadrature error;
• floating-point error;
• statistical uncertainty and bias.

These terms need not be independent, so adding all estimates in quadrature is not automatically justified. The dominant contribution should be identified by controlled variation.

Suppose a scalar result has asymptotic behavior

$$Q(h) = Q_\star + Ch^p + o(h^p).$$

Three refinements estimate the observed order:

$$p_{\mathrm{obs}} = \log_2 \left( \frac{ \lvert Q(h)-Q(h/2)\rvert }{ \lvert Q(h/2)-Q(h/4)\rvert } \right).$$

Agreement with the formal order is meaningful only after entering the asymptotic regime and controlling other errors such as box size, solver tolerance, and roundoff. Refining every parameter simultaneously can hide which source controls the result.

Error Estimates gives the error taxonomy and reporting standards. Convergence Tests gives the refinement workflow. Benchmark Problems provides exact and controlled tests for spectra, tunneling, radial equations, two-level dynamics, and periodic grids.

Recommended routes

• First stationary spectrum: Floating-Point Arithmetic → Discretization → Finite Difference Methods → Matrix Diagonalization → Convergence Tests.
• Large sparse Hamiltonian: Sparse Matrices → Sparse Eigensolvers → Conditioning and Stability → Error Estimates.
• Real-time wave packet: Discretization → Time-Stepping Methods → Matrix Exponentials Numerically → Fast Fourier Transform → Benchmark Problems.
• Radial boundary problem: ODE Solvers → Numerical Quadrature → Convergence Tests.
• Continuum PDE: Discretization → Finite Difference Methods or Spectral Methods → PDE Solvers → Error Estimates.

Page map

Page	Central question
Floating-Point Arithmetic	How does finite representation of real numbers affect quantum calculations?
Conditioning and Stability	Is sensitivity intrinsic to the problem or introduced by the algorithm?
Discretization	Which finite problem is replacing the continuum one?
Finite Difference Methods	How do local derivative stencils become operator matrices?
Spectral Methods	When do global basis expansions converge rapidly?
Numerical Quadrature	How are normalization and matrix-element integrals approximated reliably?
Matrix Diagonalization	How are dense Hermitian eigenpairs computed and checked?
Sparse Matrices	How should large local Hamiltonians be stored or applied?
Sparse Eigensolvers	How are selected eigenpairs extracted without dense diagonalization?
Time-Stepping Methods	Which time steps control stability, phase error, and unitarity?
Matrix Exponentials Numerically	When should one compute an exponential or only its action?
Fast Fourier Transform	How do discrete Fourier grids support derivatives and split evolution?
ODE Solvers	How are initial-value and shooting problems integrated and validated?
PDE Solvers	How are multidimensional quantum equations discretized and evolved?
Error Estimates	Which uncertainties limit the claimed observable?
Convergence Tests	Does the answer approach a stable limit under controlled refinement?
Benchmark Problems	Which exact or controlled cases expose implementation failures?

Common mistakes

Mistake	Correction
Treating displayed digits as accuracy	report only digits supported by an error estimate
Confusing conditioning with algorithmic stability	diagnose the problem map and the algorithm separately
Refining grid spacing while holding a too-small box fixed	vary ultraviolet and infrared controls independently
Assuming a Hermitian continuum operator yields a Hermitian array automatically	include boundary rows and the discrete inner product
Accepting an eigenvalue because the solver converged	check residuals, symmetries, subspaces, and discretization refinement
Choosing a time step only from norm drift	check phase, energy, observables, and spectral stability
Comparing wavefunctions without aligning phase and normalization	compare invariant observables or phase-aligned states
Using one benchmark that shares the production method’s failure mode	triangulate with exact limits and independent algorithms

Exercises

1. Centered second derivative

Use Taylor expansions of $\psi(x+h)$ and $\psi(x-h)$ through fifth order to derive the centered second-derivative formula and its leading error.

Solution

The expansions are

$$\begin{aligned} \psi(x+h) &=\psi+h\psi' +\frac{h^2}{2}\psi'' +\frac{h^3}{6}\psi''' +\frac{h^4}{24}\psi^{(4)} +\frac{h^5}{120}\psi^{(5)} +O(h^6),\\ \psi(x-h) &=\psi-h\psi' +\frac{h^2}{2}\psi'' -\frac{h^3}{6}\psi''' +\frac{h^4}{24}\psi^{(4)} -\frac{h^5}{120}\psi^{(5)} +O(h^6). \end{aligned}$$

Adding and solving for $\psi''$ 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).$$

Thus the approximation has leading truncation error $h^2\psi^{(4)}/12$ and is second order when the required derivatives exist.

2. What an eigenpair residual guarantees

Let $H$ be Hermitian, $\lVert\tilde v\rVert=1$, and $r=H\tilde v-\tilde\lambda\tilde v$. Show that at least one exact eigenvalue $\lambda_j$ satisfies

$$\lvert\lambda_j-\tilde\lambda\rvert \le \lVert r\rVert.$$

Solution

Expand $\tilde v=\sum_j c_jv_j$ in an orthonormal eigenbasis of $H$. Then

$$\lVert r\rVert^2 = \sum_j \lvert c_j\rvert^2 \lvert\lambda_j-\tilde\lambda\rvert^2.$$

If every eigenvalue were farther than $\lVert r\rVert$ from $\tilde\lambda$, the right-hand side would be strictly greater than

$$\lVert r\rVert^2 \sum_j\lvert c_j\rvert^2 = \lVert r\rVert^2,$$

a contradiction. The result concerns the finite Hermitian matrix; it does not bound continuum discretization error.

3. Crank–Nicolson unitarity

For time-independent Hermitian $H$, define $A=\Delta t\,H/(2\hbar)$ and

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

Show that $C^\dagger C=I$.

Solution

Because $A^\dagger=A$,

$$C^\dagger = (I+iA)(I-iA)^{-1}.$$

Every factor is a function of the same matrix $A$, so the factors commute. Therefore

$$\begin{aligned} C^\dagger C &=(I+iA)(I-iA)^{-1} (I+iA)^{-1}(I-iA)\\ &=I. \end{aligned}$$

Finite linear-solver tolerances and roundoff can still introduce small norm errors in an implementation.

4. Observed convergence order

A quantity is computed on three successively halved grids:

$$Q(h)=1.0400, \qquad Q(h/2)=1.0100, \qquad Q(h/4)=1.0025.$$

Estimate the observed order.

Solution

The consecutive differences are $0.0300$ and $0.0075$, whose ratio is $4$. Hence

$$p_{\mathrm{obs}} = \log_2 4 =2.$$

This is consistent with second-order asymptotic convergence, but more refinements and independent control of box size and solver tolerance are still needed before reporting an error bar.

References

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002.
• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997.
• R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007.
• L. N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000.
• E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration, 2nd ed., Springer, 2006.
• Y. Saad, Numerical Methods for Large Eigenvalue Problems, revised ed., SIAM, 2011.
• 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.