Companion Website

Computational Physics

From Equations to Simulation

Chapter 7

Numerical Integration

The Chapter 7 companion programs implement the main numerical integration methods discussed in the chapter in MATLAB, Python, C++, and Fortran. The resources below connect these implementations with widely used numerical libraries and provide benchmark problems for testing and validating results.

Chapter 7 codes Chapter 7 exercises & projects Chapter 7 updates & errata Chapter 7 overview

Software and library resources

The table below gives representative routines and libraries for the numerical tasks discussed in Chapter 7. Different libraries may use different algorithms, transformations, and error estimates, so check the linked documentation before choosing a routine for a particular problem.

On a narrow screen, scroll horizontally to view all language columns.

Numerical task	MATLAB	Python	C++	Fortran
Composite Newton–Cotes rules Basic trapezoidal and Simpson rules for integrating tabulated data or a function sampled on a regular grid.	`trapz` for trapezoidal integration; companion Simpson implementation	`scipy.integrate.trapezoid` `scipy.integrate.simpson`	Companion trapezoidal and Simpson implementations; these simple rules are also straightforward to code directly.	Companion trapezoidal and Simpson implementations; these simple rules are also straightforward to code directly.
General adaptive integration A good first choice for many one-dimensional integrals when automatic error control is useful.	`integral` `quadgk`	`scipy.integrate.quad`	`Boost.Math gauss_kronrod` `GSL qag / qags` (C library, callable from C++)	`Netlib QUADPACK DQAG / DQAGS` `Modern Fortran QUADPACK`
Integrable endpoint singularities Cases such as 1/√x, where the integral exists even though the integrand is singular at an endpoint.	`integral` `quadgk`	`scipy.integrate.quad`	`Boost.Math tanh_sinh` `GSL qags` (C library, callable from C++)	`QUADPACK DQAGS`
Infinite intervals Improper integrals with one or both limits extending to infinity.	`integral` `quadgk`	`quad with infinite bounds`	`Boost.Math exp_sinh` or `sinh_sinh`; `GSL qagi / qagiu / qagil` (C library, callable from C++)	`QUADPACK DQAGI`
Principal-value integrals Integrals containing a simple pole inside the interval and interpreted as a Cauchy principal value.	Analytical subtraction with `integral` or `quadgk`	`quad` `weight='cauchy'`	`GSL qawc`	`QUADPACK DQAWC`
Oscillatory integrals Integrands that change sign rapidly or contain explicit sine or cosine factors.	`quadgk` for moderate oscillation; companion Filon implementation for strongly oscillatory sine or cosine factors	`quad` `weight='sin'` or `'cos'`	`Boost.Math Ooura Fourier sine/cosine` `GSL qawo / qawf` (C library, callable from C++)	`QUADPACK DQAWO / DQAWF`
Fixed Gaussian quadrature High-order quadrature for smooth integrands when the nodes and weights are chosen in advance.	Companion Gauss–Legendre implementation or user-supplied nodes and weights	`scipy.integrate.fixed_quad` `numpy.polynomial.legendre.leggauss` `scipy.special.roots_legendre`	`Boost.Math gauss` `GSL fixed quadrature routines` (C library, callable from C++)	Library-specific Gaussian rules; custom Gauss–Legendre implementations are also common.
Multidimensional integration Nested or product quadrature in low dimensions, with other approaches often preferred as dimension increases.	`integral2` `integral3`	`nquad`, `dblquad`, `tplquad`; `cubature` (SciPy 1.15 or later); `qmc_quad` or the `vegas` package for higher dimensions	Nested one-dimensional quadrature; GSL VEGAS / MISER (C library) or the Cuba library for higher dimensions.	Nested QUADPACK calls for low dimensions; Cuba or other Monte Carlo methods for higher dimensions.

Benchmark integrals and reference values

These examples can be used to test a new implementation, compare numerical methods, or check changes to the companion programs. When an exact or high-accuracy reference value is available, compare it with the numerical result. If the method also provides an error estimate, check how well that estimate reflects the actual error.

Example	Reference value	Why it is useful
∫₀^π sin(x) dx	`2`	Basic check for composite rules and Gaussian quadrature.
∫₀¹ dx / √x	`2`	Tests an integrable endpoint singularity without requiring a difficult reference value.
∫₀¹ x⁷ dx using four-point Gauss–Legendre quadrature	`1/8`	A four-point Gauss–Legendre rule is exact through degree 7, making this a sharp test of generated nodes, weights, and interval mapping.
∫₀¹ dx / [1 + 400(x - 0.2)²]	`[atan(16) + atan(4)] / 20` `≈ 0.1417097590`	Shows whether an adaptive method can find and resolve a narrow localized feature.
∫₀^∞ e^-x sin(x) dx	`1/2`	Tests an improper integral on an infinite interval.
PV ∫_-1¹ cos(x)/(x - 0.2) dx	`≈ -0.5912784964342435`	Tests correct treatment of an interior principal-value singularity.
∫₀¹ e^-x cos(50x) dx	`[1 + e^-1(50 sin 50 - cos 50)] / 2501` `≈ -1.6717740468665e-3`	Shows how rapid oscillations can challenge a general-purpose adaptive routine.
∫₀^2π e^cos(x) dx	`2π I₀(1)` `≈ 7.954926521012845`	For an analytic periodic integrand, the periodic trapezoidal rule converges exponentially and can be unexpectedly competitive with higher-order polynomial rules.
∫₀¹∫₀¹ 1/(1+x²+y²) dy dx	`≈ 0.6395103518703110`	Tests product or nested quadrature in two dimensions.
∫₀¹∫₀^sin(x) x²/(y²+2) dy dx	`≈ 0.1034464976429315`	Tests nested integration when the inner limit depends on the outer variable.

The Lorentzian benchmark uses the same coefficient, 400, as the companion example. For a harder resolution test, replace 400 by A = 10⁴. The exact value is [atan(0.8√A) + atan(0.2√A)] / √A, so the reference remains available in closed form.

Choosing an integration method

Character of the problem	Method to consider	Practical note
Smooth function on a finite interval	Gauss–Legendre or a general adaptive library routine	If the integrand is smooth, a fixed Gaussian rule can be very efficient; a library adaptive routine is often the easiest first choice.
Automatic error control is needed	Adaptive Gauss–Kronrod	Choose tolerances that match the accuracy needed in the physical problem; requesting many more digits is usually unnecessary.
Localized peak or rapidly varying region	Adaptive Simpson or adaptive Gauss–Kronrod	Plot the integrand first. If the location of the peak is known, splitting the interval there can also help.
Integrable endpoint singularity	Tanh–sinh quadrature, another suitable variable transformation, or an adaptive routine designed for difficult integrands	Avoid evaluating an undefined endpoint value directly. The double-exponential tanh–sinh transformation is especially effective for many endpoint singularities.
Infinite interval	Variable transformation or a library routine supporting infinite limits	Check how rapidly the tail decays and inspect the transformed integrand if a change of variable is used.
Simple pole inside the interval	Analytical subtraction or a dedicated principal-value routine	Ordinary quadrature across the pole is not appropriate; use a principal-value treatment.
Rapid oscillations	Filon-type, weighted, or other oscillatory-specific quadrature	General adaptive methods may spend many function evaluations trying to follow the oscillations.
Periodic analytic integrand over a full period	Periodic trapezoidal rule	The error can decrease exponentially with the number of equally spaced points, making this simple rule exceptionally effective.
Tabulated data	Composite trapezoidal or Simpson rule using the supplied coordinates	The trapezoidal rule works naturally with irregular spacing; check the spacing requirements of any Simpson implementation.
Smooth integrand in two dimensions	Product Gaussian quadrature or adaptive multidimensional routine	Tensor-product rules work well in low dimensions, but the number of points grows quickly as dimension increases.
Non-rectangular domain (variable inner limits)	Nested integration	The accuracy of the inner integral affects the outer result, so the two tolerances should be chosen consistently.
High-dimensional integral	Monte Carlo or quasi-Monte Carlo integration	For finite-variance integrands, ordinary Monte Carlo error decreases approximately as `N^-1/2`, with an exponent independent of dimension; the error constant can still worsen as dimension increases.

Selected references and documentation

R. Piessens, E. de Doncker-Kapenga, C. W. Überhuber, and D. Kahaner, QUADPACK: A Subroutine Package for Automatic Integration, Springer-Verlag, 1983. The source of the adaptive routines underlying many of the library integrators listed above.
G. E. Forsythe, M. A. Malcolm, and C. B. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, 1977. The QUANC8 routine used as the basis for the companion implementation is described in this work.
P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, 1984. A comprehensive reference on classical quadrature formulas, error analysis, and special-purpose methods.
G. H. Golub and J. H. Welsch, “Calculation of Gauss Quadrature Rules,” Mathematics of Computation 23 (1969), 221–230. A classic method for computing Gaussian quadrature nodes and weights.
L. F. Shampine, “Vectorized Adaptive Quadrature in MATLAB,” Journal of Computational and Applied Mathematics 211 (2008), 131–140. Describes the adaptive strategy behind MATLAB's modern vectorized quadrature approach.
L. N. G. Filon, “On a Quadrature Formula for Trigonometric Integrals,” Proceedings of the Royal Society of Edinburgh 49 (1928), 38–47. The original basis for the Filon-type oscillatory method used in the companion program.
G. P. Lepage, “A New Algorithm for Adaptive Multidimensional Integration,” Journal of Computational Physics 27 (1978), 192–203. doi:10.1016/0021-9991(78)90004-9 Introduces the VEGAS adaptive Monte Carlo algorithm.
W. Gander and W. Gautschi, “Adaptive Quadrature—Revisited,” BIT Numerical Mathematics 40 (2000), 84–101. doi:10.1023/A:1022318402393 Compares and develops practical adaptive quadrature strategies.
L. N. Trefethen and J. A. C. Weideman, “The Exponentially Convergent Trapezoidal Rule,” SIAM Review 56 (2014), 385–458. doi:10.1137/130932132 Explains the exceptional performance of the trapezoidal rule for periodic analytic functions and related problems.
T. Hahn, “Cuba—A Library for Multidimensional Numerical Integration,” Computer Physics Communications 168 (2005), 78–95. doi:10.1016/j.cpc.2005.01.010 Describes the Vegas, Suave, Divonne, and Cuhre algorithms provided by the Cuba library.
NIST Digital Library of Mathematical Functions: Quadrature (§3.5)
MATLAB integration and differentiation documentation
SciPy integration documentation
GNU Scientific Library: Numerical Integration
Netlib QUADPACK archive

← Back to Chapter 7