Numerical Integration
Approximate definite integrals using quadrature rules — visualize rectangles, trapezoids, and parabolas, compare convergence rates, and compute subdivisions needed for a given tolerance.
The Numerical Integration Problem
Given a function \(f\) on \([a,b]\), numerical quadrature approximates the definite integral \(\int_a^b f(x)\,dx\) using a weighted sum of function values. Every quadrature rule has the form:
where \(x_i\) are the nodes (sample points) and \(w_i\) are the weights. The choice of nodes and weights determines the degree of exactness — the highest-degree polynomial that the rule integrates exactly.
“Degree” in this tool means degree of exactness (also called the degree of precision) — the largest \(d\) such that the rule is exact for every polynomial of degree \(\le d\). The badge “degree 3” on a method means it integrates cubics exactly.
This is not the same as the order of convergence \(p\) in the error \(O(h^{p})\), which describes how fast the composite error shrinks as the step \(h\to 0\). The two are related: a composite rule of degree \(d\) converges at order \(p=d+1\) when \(d\) is odd and \(p=d+2\) when \(d\) is even (symmetry buys one extra order). For example, Simpson’s rule has degree 3 but converges at order 4.
Families of Quadrature Rules
The rules in this tool fall into a few families, ordered roughly by how cleverly they place their nodes. Each rule’s formula and a detailed explanation live in the Methods panel below.
- Riemann sums sample one corner of each subinterval — the crudest rules, error \(O(h)\).
- Newton–Cotes rules integrate the polynomial through equispaced nodes. Closed rules use the panel endpoints (trapezoidal, Simpson, Boole…); open rules use only interior nodes, handy when \(f\) misbehaves at an endpoint.
- Gaussian rules choose both nodes and weights to maximise the degree of exactness (Gauss–Legendre), optionally pinning one or both endpoints (Radau, Lobatto).
- Spectral & double-exponential rules (Clenshaw–Curtis, tanh–sinh) are single global formulas that converge faster than any fixed power of \(h\) for smooth integrands.
- Weighted Gaussian rules (Chebyshev, Jacobi) integrate \(\int_a^b w(x)f(x)\,dx\) with a singular weight \(w\) built into the nodes — exact for the weight, so they nail integrable endpoint singularities.
- Adaptive & extrapolation methods (adaptive Simpson, Romberg) refine or extrapolate automatically until a tolerance is met.
Error Estimates
For the composite rules on \([a,b]\) with step size \(h=(b-a)/n\):
- Left / Right: \(|E| \le \dfrac{(b-a)^2}{2n}\,\max|f'|\) — \(O(h)\)
- Midpoint: \(|E| \le \dfrac{(b-a)^3}{24n^2}\,\max|f''|\) — \(O(h^2)\)
- Trapezoidal: \(|E| \le \dfrac{(b-a)^3}{12n^2}\,\max|f''|\) — \(O(h^2)\)
- Simpson: \(|E| \le \dfrac{(b-a)^5}{180n^4}\,\max|f^{(4)}|\) — \(O(h^4)\)
- Composite degree-\(p\) rule: error \(O(h^{p+1})\) for odd \(p\), \(O(h^{p+2})\) for even \(p\) — one extra order from symmetry.
Extended-precision error. A double-precision computation caps the resolvable error near machine epsilon (\(\approx 10^{-16}\)): once a high-order or spectral rule reaches that floor, its true error is buried in round-off. To see past it, nearly every rule — the Riemann, Newton–Cotes, Gauss–Legendre, Gauss–Lobatto/Radau, Romberg, Clenshaw–Curtis and tanh–sinh rules, plus the weighted Gaussian (Chebyshev, Jacobi), infinite-domain (Hermite, Laguerre) and double-exponential (sinh–sinh, exp–sinh) families — re-evaluates both the approximation \(Q_n\) and the true value \(I\) with ~40-digit arithmetic, so the absolute error \(|Q_n-I|\) is resolved 10+ orders past the \(10^{-16}\) double-precision floor (and shown in concise scientific notation — the exponent already carries the magnitude). This uses extended-precision nodes and weights computed on the fly: Gauss–Legendre/Hermite/Laguerre/Jacobi abscissae come from the eigenvalues of the recurrence’s Jacobi matrix (Golub–Welsch, via a high-precision QL iteration, with a Spouge \(\Gamma\) function for the Jacobi weights), and the Chebyshev and tanh–sinh/double-exponential nodes from their closed forms. The value cards show the first 16 digits. The same extended precision drives the convergence plot (whose curves keep falling past the \(10^{-16}\) floor instead of flattening) and the tolerance calculator (which can now solve for accuracies well below machine epsilon). The adaptive-Simpson and Gauss–Kronrod rules keep the standard double-precision estimate (the latter reports its own embedded \(|I_{15}-I_{7}|\) estimator).
This Tool
- Pick a rule from the grouped method menu above — families are organised as Riemann sums, closed/open Newton–Cotes, Gauss–Legendre/Lobatto/Radau, spectral & double-exponential, weighted Gaussian (Chebyshev / Jacobi), and adaptive & extrapolation, each with its degree of precision.
- Choose a built-in example or enter a custom integrand, interval, and number of subintervals. For Gauss–Jacobi, also set the weight exponents \(\alpha,\beta\) (the fields appear when that rule is selected).
- Click Compute to see the quadrature formula, the numerical approximation, exact value (when known), absolute error, and the integration graph.
- Convergence plot: A log-log plot of absolute error vs. number of subintervals, confirming each rule’s theoretical rate; spectral rules visibly fall faster than any straight reference line.
- Tolerance calculator: Enter a desired accuracy \(\varepsilon\) and the tool finds the smallest number of subintervals \(n\) for which the measured error \(|Q_n - I|\) drops below \(\varepsilon\), via doubling-then-bisection.
Riemann sums
Left Riemann Sumdegree 0
Approximate each panel by a rectangle whose height is the value at its left endpoint.
Exact only for constants; the error is \(O(h)\) because the rectangle ignores the slope across the panel.
Right Riemann Sumdegree 0
The same idea using the right endpoint of each panel.
Also \(O(h)\); left and right sums bracket the integral for a monotone \(f(x)\).
Closed Newton–Cotes (equispaced, endpoints included)
Trapezoidal Ruledegree 1
Connect consecutive points by straight lines — integrate the piecewise-linear interpolant.
Exact for linears; error \(O(h^2)\). The seed for Romberg.
Simpson’s 1/3 Ruledegree 3
Fit a parabola through each pair of panels (needs an even \(n\)).
Exact through cubics (a bonus order from symmetry), error \(O(h^4)\).
Simpson’s 3/8 Ruledegree 3
Fit a cubic through each group of three panels (\(n\) a multiple of 3).
Same degree 3 / \(O(h^4)\) as the 1/3 rule; handy when \(n\) is a multiple of 3 rather than 2.
Boole’s Ruledegree 5
Fit a quartic through each group of four panels.
Exact through quintics, error \(O(h^6)\).
Closed 6-point Ruledegree 5
The five-panel (six-node) closed Newton–Cotes formula.
Degree 5, error \(O(h^6)\). Higher closed rules start to develop negative weights and lose reliability.
Open Newton–Cotes (interior nodes only)
Midpoint Ruledegree 1
One rectangle per panel at its midpoint — the simplest open rule, never touching the endpoints.
Exact for linears, error \(O(h^2)\); usually a touch more accurate than the trapezoidal rule.
Open 2-point Ruledegree 1
Two interior nodes per panel of three; endpoints are never sampled.
Degree 1, error \(O(h^2)\). Useful when \(f(x)\) is singular at \(a\) or \(b\).
Milne’s Rule (open 3-point)degree 3
Three interior nodes per panel of four.
Degree 3, error \(O(h^4)\) — the open counterpart of Simpson. Note the negative middle weight.
Open 4-point Ruledegree 3
Four interior nodes per panel of five.
Degree 3, error \(O(h^4)\).
Gauss–Legendre (optimal free nodes)
Gauss–Legendre 2-pointdegree 3
Choosing both nodes and weights, an \(m\)-point rule is exact for degree \(\le 2m-1\). Nodes are roots of the Legendre polynomial on \([-1,1]\), mapped to each subinterval.
With \(m=2\) the nodes are the roots of \(P_2(x)=\tfrac12(3x^2-1)\):
Degree 3, composite error \(O(h^4)\) at just two evaluations per panel.
Gauss–Legendre 3-pointdegree 5
Three Legendre nodes per subinterval, with nodes and weights again chosen together to maximise exactness; applied compositely over the \(n\) panels with the same mapping as the 2-point rule.
Exact for degree \(\le 2m-1=5\), composite error \(O(h^6)\).
Gauss–Legendre 5-pointdegree 9
Five Legendre nodes per subinterval — the highest-order plain Gauss rule here. As with every Gauss–Legendre rule, both the nodes (roots of the degree-5 Legendre polynomial on \([-1,1]\)) and the weights are chosen to integrate the highest possible polynomial degree, then mapped to each panel.
So accurate that a handful of panels already reach the round-off floor on smooth integrands. Also used internally as the high-accuracy reference fallback.
Gauss–Legendre 7-pointdegree 13
Seven Legendre nodes per subinterval — the same 7 nodes that anchor the Gauss–Kronrod pair below. Exact for degree \(\le 2m-1=13\).
Two more nodes than the 5-point rule, lifting the exactness from degree 9 to 13; on smooth integrands it reaches round-off with just one or two panels.
Gauss–Kronrod 7–15error est.
The 15-point Kronrod rule \(K_{15}\) extends the 7-point Gauss–Legendre rule \(G_7\): it keeps all 7 Gauss nodes and adds 8 optimally-placed ones, integrating polynomials up to degree 22 (versus 13 for \(G_7\)).
Because the two rules share the 7 Gauss nodes, the difference \(\bigl|I_{15}-I_{7}\bigr|\) costs no extra evaluations yet gives a reliable a posteriori error estimate — the engine behind QUADPACK’s adaptive integrators. Here the pair is applied compositely over \(n\) panels and the panel estimates are summed; the result card reports the embedded estimate next to the true error.
Gauss–Lobatto & Gauss–Radau (constrained nodes)
Gauss–Lobatto 3-pointdegree 3
Lobatto rules fix both endpoints \(\pm1\) and optimise the interior nodes; an \(m\)-point rule is exact for degree \(\le 2m-3\).
This 3-point rule is exactly Simpson’s 1/3. Fixing endpoints lets adjacent panels share evaluations.
Gauss–Lobatto 4-pointdegree 5
Four nodes per subinterval with both endpoints \(\pm1\) fixed and the two interior nodes optimised; exact for degree \(\le 2m-3=5\).
Degree 5, composite error \(O(h^6)\).
Gauss–Lobatto 5-pointdegree 7
Five nodes per subinterval pinning both endpoints \(\pm1\), the three interior nodes optimised; exact for degree \(\le 2m-3=7\).
Degree 7, composite error \(O(h^8)\).
Gauss–Radau 2-pointdegree 2
Radau rules fix one endpoint (here \(-1\)); an \(m\)-point rule is exact for degree \(\le 2m-2\). The rule is asymmetric.
Degree 2.
Gauss–Radau 3-pointdegree 4
Three nodes per subinterval with the left endpoint \(-1\) fixed and the other two optimised; exact for degree \(\le 2m-2=4\). Asymmetric, like every Radau rule.
Degree 4, composite error \(O(h^5)\).
Spectral & double-exponential (global rules)
Clenshaw–Curtisspectral
A single global rule at the Chebyshev extreme points; integrates the cosine interpolant of \(f(x)\) exactly.
Converges spectrally (faster than any \(O(h^p)\)) for smooth \(f(x)\); \(n\) sets \(N\), giving \(N+1\) nodes clustered toward the endpoints.
Tanh–sinh (double exponential)double-exp
Substitute \(x=\tanh\!\big(\tfrac{\pi}{2}\sinh t\big)\) and apply the trapezoidal rule in \(t\).
The endpoints are crushed double-exponentially, so convergence is almost exponential and integrable endpoint singularities (e.g. \(\int_0^1 x^{-1/2}\,dx\)) are handled gracefully. Nodes cluster toward the centre.
Weighted Gaussian (singular weight built in)
Gauss–Chebyshev, 1st kindweighted
Targets the weighted integral with \(w(x)=1/\sqrt{(b-x)(x-a)}\) (Jacobi \(\alpha=\beta=-\tfrac12\)); nodes and weights are closed-form.
All weights equal \(\pi/n\). Exact for \(f(x)\) of degree \(\le 2n-1\); e.g. \(\int_{-1}^{1}\tfrac{dx}{\sqrt{1-x^2}}=\pi\) with \(f(x)\equiv1\).
Gauss–Chebyshev, 2nd kindweighted
Weight \(w(x)=\sqrt{(b-x)(x-a)}\) (Jacobi \(\alpha=\beta=+\tfrac12\)).
Example: \(\int_{-1}^{1}\sqrt{1-x^2}\,dx=\tfrac{\pi}{2}\).
Gauss–Jacobi (\(\alpha,\beta\))weighted
General weight \(w(x)=(b-x)^{\alpha}(x-a)^{\beta}\), \(\alpha,\beta>-1\). Nodes/weights come from the symmetric tridiagonal Jacobi matrix of the recurrence coefficients (Golub–Welsch).
The nodes \(x_k\) are the roots of the Jacobi polynomial \(P_n^{(\alpha,\beta)}\) on \([-1,1]\), and \(\Gamma\) is the gamma function. With \(\alpha=0,\ \beta=-\tfrac12\) on \([0,1]\) the weight is \(x^{-1/2}\), so \(\int_0^1 x^{-1/2}\,dx=2\) comes out exact with \(f(x)\equiv1\). Set \(\alpha,\beta\) in the fields that appear when this rule is selected.
Infinite-domain
Gauss–Hermiteweight \(e^{-x^2}\) on \(\mathbb{R}\)
Integrates over the whole real line against the Gaussian weight \(w(x)=e^{-x^2}\); the \([a,b]\) inputs do not apply. Nodes are the roots of the Hermite polynomial \(H_n\) and the weights follow from Golub–Welsch (the Jacobi matrix has zero diagonal and off-diagonals \(\sqrt{k/2}\)).
Exact when \(f(x)\) is a polynomial of degree \(\le 2n-1\); e.g. \(\int_{-\infty}^{\infty} e^{-x^2}\,dx=\sqrt{\pi}\) with \(f(x)\equiv1\). The nodes are symmetric about 0.
Gauss–Laguerreweight \(e^{-x}\) on \([0,\infty)\)
Integrates over the half-line \([0,\infty)\) against \(w(x)=e^{-x}\). Nodes are the roots of the Laguerre polynomial \(L_n\); the Jacobi matrix has diagonal \(2k+1\) and off-diagonals \(k\) (Golub–Welsch).
Exact for \(f(x)\) of degree \(\le 2n-1\); e.g. \(\int_0^{\infty} e^{-x}\cos x\,dx=\tfrac12\). Ideal for integrands that already carry an exponential decay.
Sinh–sinh\(\mathbb{R}\), unweighted
An unweighted double-exponential rule for the whole line: it integrates \(f\) itself (no built-in weight), unlike Gauss–Hermite. The map \(x=\sinh\!\big(\tfrac{\pi}{2}\sinh t\big)\) sends \(t\in\mathbb{R}\) to \(x\in\mathbb{R}\) and crushes both tails double-exponentially, so a plain trapezoidal rule in \(t\) converges almost exponentially.
Use it for integrands with algebraic or exponential decay that aren’t of the form \(e^{-x^2}\!\times\)(smooth), e.g. \(\int_{-\infty}^{\infty}\tfrac{dx}{1+x^2}=\pi\).
Exp–sinh\([0,\infty)\), unweighted
The half-line counterpart: \(x=\exp\!\big(\tfrac{\pi}{2}\sinh t\big)\) maps \(t\in\mathbb{R}\) to \(x\in(0,\infty)\), again with a trapezoidal rule in \(t\). Unweighted, so it integrates \(f\) directly over \([0,\infty)\) without the \(e^{-x}\) weight that Gauss–Laguerre assumes.
Handles endpoint singularities and slow decay gracefully, e.g. \(\int_0^{\infty}\tfrac{dx}{1+x^2}=\tfrac{\pi}{2}\) or \(\int_0^{\infty} e^{-x}/\sqrt{x}\,dx=\sqrt{\pi}\).
Highly oscillatory
Ordinary rules need several nodes per wavelength, so their cost and error grow with the frequency \(\omega\). These rules instead integrate the oscillation analytically, staying accurate as \(\omega\to\infty\). Enter the smooth amplitude \(f(x)\) in the integrand box; the oscillator \(\sin(\omega x)\) or \(\cos(\omega x)\) and the frequency \(\omega\) are set with the extra controls.
Filon (Simpson)\(\int_a^b f\,e^{i\omega x}\)
On each pair of panels the amplitude \(f\) is approximated by a quadratic and multiplied by \(\{\sin,\cos\}(\omega x)\), which is then integrated exactly. The panel weights \(\alpha,\beta,\gamma\) depend on \(\theta=\omega h\) and absorb the oscillation, so the rule is exact for the oscillator times any quadratic and its error is bounded uniformly in \(\omega\) — where composite Simpson explodes once \(\omega h\) is no longer \(\ll 1\).
Finite interval \([a,b]\). Try \(\int_0^1 x\sin(40x)\,dx\) or \(\int_{-1}^{1} e^{x}\cos(30x)\,dx\) — both nailed with a handful of panels. (Filon 1928; Abramowitz & Stegun 25.4.47.)
Ooura–Mori DE\(\int_0^\infty f\,e^{i\omega x}\)
A double-exponential rule for Fourier-type integrals on \([0,\infty)\). The change of variable \(x=\tfrac{M}{\omega}\phi(t)\), \(M=\pi/h\), is chosen so the far nodes fall asymptotically on the zeros of \(\sin/\cos(\omega x)\): the slowly-decaying oscillatory tail telescopes to zero instead of being sampled blindly, succeeding where sinh–sinh stalls on the same integrand.
Half-line \([0,\infty)\). Try \(\int_0^{\infty}\tfrac{\sin x}{x}\,dx=\tfrac{\pi}{2}\) or \(\int_0^{\infty}\tfrac{\cos 2x}{1+x^2}\,dx=\tfrac{\pi}{2e^{2}}\). (Ooura & Mori, J. Comput. Appl. Math. 112 (1999) 229–241.)
Ooura–Mori DE, full line\(\int_{-\infty}^\infty f\,e^{i\omega x}\)
The whole-line Fourier integral, obtained by folding the negative half onto \([0,\infty)\) and running the same double-exponential rule on the symmetrised amplitude. Because \(\sin\) is odd and \(\cos\) even,
so a single half-line Ooura–Mori sweep on \(g(x)=f(x)\mp f(-x)\) captures both tails, keeping the far nodes on the zeros of \(\sin/\cos(\omega x)\). Full line \((-\infty,\infty)\). Try \(\int_{-\infty}^{\infty}\tfrac{\cos 2x}{1+x^2}\,dx=\tfrac{\pi}{e^{2}}\) or \(\int_{-\infty}^{\infty}\tfrac{x\sin x}{1+x^2}\,dx=\tfrac{\pi}{e}\).
Adaptive & extrapolation
Adaptive Simpsonadaptive
Recursively bisect a panel only where it is needed, comparing a whole-panel Simpson estimate \(S\) with its split.
Concentrates nodes where the integrand is hardest instead of spreading them uniformly.
Rombergextrapolation
Richardson-extrapolate the composite trapezoidal rule at successively halved steps.
Each column cancels the next error term, so the diagonal \(R_{m,m}\) converges very rapidly for smooth \(f(x)\).
Examples
Cite this tool
Kapita, S. (2026). Numerical Integration. Math Tools. https://doi.org/10.5281/zenodo.20981237
Kapita, Shelvean. "Numerical Integration." Math Tools, 2026, doi.org/10.5281/zenodo.20981237.
@online{kapita2026integration,
author = {Shelvean Kapita},
title = {{Numerical Integration}},
year = {2026},
organization = {Math Tools},
doi = {10.5281/zenodo.20981237},
url = {https://doi.org/10.5281/zenodo.20981237}
}