Numerical Quadrature

Numerical quadrature approximates integrals by weighted sums. In quantum mechanics, quadrature appears whenever a continuous wavefunction is normalized, an expectation value is evaluated, a matrix element is assembled, or a probability is integrated over a region.

The basic form is

$$\int_a^b f(x)\,dx \approx \sum_{j=0}^{N} w_j f(x_j),$$

where $x_j$ are nodes and $w_j$ are weights.

The weights are part of the numerical inner product. Dropping them can change normalization, Hermiticity, and expectation values.

Why Quantum Mechanics Needs It

Quadrature is used to compute:

• normalization integrals $\int \lvert\psi(x)\rvert^2\,dx$;
• probabilities over finite regions;
• expectation values $\langle\psi\lvert A\rvert\psi\rangle$;
• overlap matrices $\langle\phi_m,\phi_n\rangle$;
• Hamiltonian matrix elements $\langle\phi_m,H\phi_n\rangle$;
• Fourier-type oscillatory amplitudes;
• radial integrals in central potentials;
• thermal traces and density-of-states approximations.

Even when the final calculation is a matrix problem, quadrature may be hidden inside how the matrix was built.

Rectangle and Trapezoidal Rules

On a uniform grid

$$x_j=a+jh, \qquad h=\frac{b-a}{N},$$

the composite trapezoidal rule is

$$\int_a^b f(x)\,dx \approx h \left[ \frac12 f(x_0) + \sum_{j=1}^{N-1}f(x_j) + \frac12 f(x_N) \right].$$

For smooth nonperiodic functions, its error is typically $O(h^2)$. For smooth periodic functions sampled over a full period, the trapezoidal rule can converge much faster because the endpoint behavior matches.

For a grid wavefunction, the same rule gives a discrete norm:

$$\lVert\psi\rVert_h^2 = h \left[ \frac12\lvert\psi_0\rvert^2 + \sum_{j=1}^{N-1}\lvert\psi_j\rvert^2 + \frac12\lvert\psi_N\rvert^2 \right].$$

This is why normalization on a grid is a quadrature question, not just a vector-length question.

Simpson’s Rule

For an even number of subintervals, Simpson’s rule is

$$\int_a^b f(x)\,dx \approx \frac{h}{3} \left[ f(x_0) + f(x_N) + 4\sum_{\substack{j=1\\ j\ \mathrm{odd}}}^{N-1} f(x_j) + 2\sum_{\substack{j=2\\ j\ \mathrm{even}}}^{N-2} f(x_j) \right].$$

For sufficiently smooth functions, Simpson’s rule has error $O(h^4)$. It is a good default for many one-dimensional smooth integrals, but it assumes equally spaced points and an even number of subintervals.

Higher formal order does not guarantee a better answer if the integrand has endpoint singularities, discontinuities, unresolved oscillations, or roundoff-dominated cancellation.

Gaussian Quadrature

Gaussian quadrature chooses nodes and weights so that polynomials up to high degree are integrated exactly. In a standard weighted form,

$$\int_a^b f(x)w(x)\,dx \approx \sum_{j=1}^{N} W_j f(x_j).$$

The nodes are roots of orthogonal polynomials associated with the weight $w(x)$. An $N$-point Gaussian rule is exact for polynomials of degree up to $2N-1$ when the integrand has the form polynomial times the weight.

Common cases include:

• Gauss–Legendre quadrature on finite intervals;
• Gauss–Hermite quadrature for weights involving $e^{-x^2}$;
• Gauss–Laguerre quadrature for weights involving $e^{-x}$ on $[0,\infty)$;
• Gauss-Jacobi quadrature for endpoint power-law weights.

This connects directly to Orthogonal Polynomials.

Matrix Elements and Overlaps

For basis functions $\phi_m$ and $\phi_n$, a matrix element such as

$$V_{mn} = \int \phi_m^\ast(x)V(x)\phi_n(x)\,dx$$

is often evaluated by quadrature:

$$V_{mn} \approx \sum_j w_j \phi_m^\ast(x_j)V(x_j)\phi_n(x_j).$$

If the same quadrature rule is used consistently, Hermitian matrix elements can remain Hermitian up to roundoff when $V$ is real. If inconsistent quadrature is used for different entries, artificial non-Hermiticity can appear.

For nonorthogonal bases, the overlap matrix

$$S_{mn} = \int \phi_m^\ast(x)\phi_n(x)\,dx$$

is also a quadrature object. Poor quadrature can make a good basis look ill conditioned or hide a genuinely ill-conditioned basis.

Infinite Domains

Many quantum integrals are over $\mathbb R$ or $[0,\infty)$. Common strategies include:

• truncate the domain and check tail error;
• map the infinite interval to a finite interval;
• use a Gaussian rule matched to the weight, such as Hermite or Laguerre;
• split the integral into regions with different behavior;
• subtract or factor known asymptotic behavior.

For a localized bound state, domain truncation may be harmless once the tail is negligible. For scattering states or slowly decaying functions, the same truncation can be misleading.

Oscillatory Integrals

Oscillatory integrals such as

$$\int a(x)e^{iS(x)/\hbar}\,dx$$

are delicate because positive and negative phase contributions cancel. A rule that samples too coarsely can return a plausible but wrong small number.

Good diagnostics include:

• refine until several points resolve each local wavelength;
• compare with analytic stationary-phase estimates when available;
• split the integral at stationary points or singularities;
• avoid subtracting two large noisy partial sums to obtain a tiny amplitude.

For asymptotic background, see Asymptotic Analysis.

Multidimensional Quadrature

Tensor-product quadrature in $d$ dimensions grows rapidly:

$$N_{\mathrm{points}} = N^d.$$

This curse of dimensionality is one reason many-body quantum integrals are hard. Low-dimensional tensor grids are useful for benchmark problems, but high-dimensional problems often require separability, sparse grids, Monte Carlo, variational structure, or problem-specific approximations.

For statistical sampling methods, see Monte Carlo Basics.

Error Checks

Useful checks include:

• refine the grid or increase the quadrature order;
• compare trapezoidal, Simpson, and Gaussian rules on the same integral;
• test against an exactly known normalization or expectation value;
• check symmetry, realness, and positivity properties;
• monitor cancellation by comparing the result with the sum of magnitudes;
• verify that the quadrature weights match the grid and coordinate measure.

For radial integrals, remember the measure. In three dimensions,

$$d^3x = r^2\sin\theta\,dr\,d\theta\,d\phi.$$

Forgetting the measure factor is not a quadrature error; it is the wrong integral.

For combining quadrature error with roundoff, solver residuals, and statistical uncertainty, see Error Estimates.

Common Mistakes

• Normalizing a grid wavefunction without the quadrature weights.
• Using Simpson’s rule with an odd number of subintervals.
• Applying a high-order rule to a nonsmooth integrand and expecting high-order convergence.
• Forgetting Jacobian factors after a coordinate change.
• Treating an oscillatory integral as resolved without checking local wavelength.
• Using inconsistent quadrature for matrix elements that should form a Hermitian matrix.
• Confusing quadrature error with eigensolver error after the matrix has been built.
• Assuming an infinite-domain integral is safe because the plotted tail looks small.

Cross-Links

• Discretization
• Floating-Point Arithmetic
• Spectral Methods
• Error Estimates
• Inner Products
• Expectation Values
• Probability Densities
• Orthogonal Polynomials
• Asymptotic Analysis
• Monte Carlo Basics

References

• P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, 1984.
• 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.
• J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 3rd ed., Springer, 2002.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
• M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions, Dover, 1965.

Exercises

1. Write the trapezoidal approximation to the normalization integral of a wavefunction sampled at $N+1$ equally spaced points on $[a,b]$.

Solution

With $h=(b-a)/N$ and samples $\psi_j=\psi(x_j)$, the trapezoidal normalization is

$$\int_a^b \lvert\psi(x)\rvert^2\,dx \approx h \left[ \frac12\lvert\psi_0\rvert^2 + \sum_{j=1}^{N-1}\lvert\psi_j\rvert^2 + \frac12\lvert\psi_N\rvert^2 \right].$$

2. Why does Simpson’s rule require an even number of subintervals?

Solution

Composite Simpson’s rule fits a quadratic polynomial across pairs of adjacent subintervals. Therefore the interval must be divided into an even number of subintervals so they can be grouped into pairs. Equivalently, the number of grid points must be odd.

3. What does it mean that an $N$-point Gaussian quadrature rule is exact for polynomials up to degree $2N-1$?

Solution

It means that, for any polynomial $p(x)$ of degree at most $2N-1$, the quadrature sum gives the exact weighted integral:

$$\int_a^b p(x)w(x)\,dx = \sum_{j=1}^N W_jp(x_j).$$

The nodes and weights are chosen to satisfy this exactness property.

4. Why are Gauss–Hermite rules natural for integrals involving oscillator wavefunctions?

Solution

Harmonic-oscillator wavefunctions contain Gaussian factors times Hermite polynomials. Gauss–Hermite quadrature is designed for integrals with Gaussian weights on the real line, so it matches the natural weight structure of oscillator matrix elements and overlaps.

5. Give one reason an oscillatory integral can appear converged while still being wrong.

Solution

If the grid does not resolve the local oscillation wavelength, samples may miss cancellations or stationary-phase regions. The resulting sum can settle to a small-looking number that is controlled by aliasing or sampling error rather than by the true integral. Refining the grid and comparing with phase-based estimates are essential checks.