Integrate is a small, lightweight Rust library for performing numerical integration of real-valued functions. It is designed to integrate functions, providing a simple and efficient way to approximate definite integrals using various numerical methods. Numerical integration is concerned with developing algorithms to approximate the integral of a function $f(x)$. The most commonly used algorithms are Newton-Cotes formulas, Romberg’s method, Gaussian quadrature, and to lesser extents Hermite’s formulas and certain adaptive techniques.

Caveats

All of the numerical integration techniques listed above assume that the integrand is continuous on the interval of integration. The error estimates generally require that the integrands are differentiable of whatever order is required so that the formula for the error estimate makes sense. Below is a check list which one should verify before using any of the numerical algorithms.

Integrand check list

• Are there any singularities of the integrand in the interval of integration?
  • If there are, then can they be removed?
  • A function which has jump discontinuities can be integrated by splitting the interval of integration into subintervals on which the function is continuous, and then numerically integrate over each subinterval and add the results.
  • Using integration by parts certain types of singular integrals can be expressed as the sum of a function and a non-singular integral.
  • In other cases, an approximation of the integral in a small neighborhood of the singularity can be obtained by hand and the interval of integration can be split into three intervals in which numerical integration can be used on two of them.
• Does the function oscillate over the region of integration? If it does, then make sure that the step size is chosen to be smaller than the wave length of the function. The interval of integration can also be split into subintervals in which each subinterval is half a wave length and the algorithm is applied to each subinterval.

Getting Started

Installation

Add integrate to your Cargo.toml:

[dependencies]
integrate = "0.1"

Your First Integral

Let’s compute \(\int_0^1 e^x \, dx = e - 1 \approx 1.71828\) using the trapezoidal rule:

use integrate::newton_cotes::trapezoidal_rule;

let result = trapezoidal_rule(|x: f64| x.exp(), 0.0_f64, 1.0_f64, 1_000_usize);
println!("Result:   {:.6}", result);
println!("Exact:    {:.6}", std::f64::consts::E - 1.0);

Using the Prelude

For convenience, you can import all methods at once with:

use integrate::prelude::*;

This exposes all 11 integration functions directly.

Choosing a Method

Generics: f32 and f64

All Newton-Cotes and Romberg methods are generic over float types. You can pass f32 functions:

use integrate::newton_cotes::simpson_rule;

let result = simpson_rule(|x: f32| x * x, 0.0_f32, 1.0_f32, 1_000_usize);
println!("{:.6}", result); // ≈ 0.333333

Newton-Cotes methods

Newton-Cotes methods approximate the integral of a function by summing a weighted combination of the function evaluated at equally-spaced points and nodes. If the endpoints of the interval of integration are excluded from the sum, the method is called an open Newton-Cotes method and if the endpoints are included the method is called a closed Newton-Cotes method.

• Rectangle Rule
• Trapezoidal Rule
• Simpson’s Rule
• Newton’s 3/8 Rule

Rectangle rule

Motivation

The rectangle rule (also called the midpoint rule) is the simplest Newton-Cotes formula for numerical integration. It approximates each subinterval’s contribution by the function value at its midpoint, making it easy to implement and exact for constant functions. Despite its low order of accuracy, it is a useful baseline and requires no endpoint evaluations.

Example

use integrate::newton_cotes::rectangle_rule;

let square = |x: f64| x * x;

let a = 0.0;
let b = 1.0;

let num_steps: usize = 1_000_000;

let integral = rectangle_rule(square, a, b, num_steps);
println!("{}", integral);

Understanding the Rectangle rule

Composite formula

The rectangle rule approximates the integral of \(f(x)\) on \([a, a+h]\) by the signed area of the rectangle with width \(h\) and height \(f\left(a + h/2\right)\).

The composite rule decomposes \([a, b]\) into \(n\) subintervals of equal length \(h = (b-a)/n\) and sums the result over each subinterval:

\[ R_h(f) = h \left[ f\left(a + \tfrac{h}{2}\right) + f\left(a + \tfrac{3h}{2}\right) + \cdots + f\left(b - \tfrac{h}{2}\right) \right] \]

Euler-Maclaurin expansion

The Euler-Maclaurin summation formula relates \(R_h(f)\) to the exact integral:

\[ \begin{align} R_h(f) &= \int_{a}^{b} f(x)\,dx - \frac{h^2}{24} \left[ f^\prime(b) - f^\prime(a) \right] + \frac{7h^4}{5760} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &+ \cdots + K h^{2p-2} \left[ f^{(2p-3)}(b) - f^{(2p-3)}(a) \right] + O(h^{2p}) \end{align} \]

The \(O(h^{2p})\) remainder means that for an infinitely differentiable function whose first \(2p-3\) derivatives vanish at both endpoints, the theorem guarantees only that

\[ \lim_{h \to 0} \left| \frac{R_h(f) - \int_{a}^{b} f(x),dx}{h^{2p}} \right| < M \]

for some \(M > 0\), not that the rule is exact.

The expansion also shows that \(n\) should be chosen so that \(h = (b-a)/n < 1\). For example, if \(h = 0.1\):

\[ \begin{split} R_{0.1}(f) &= \int_{a}^{b} f(x)\,dx - 0.00042 \left[ f^\prime(b) - f^\prime(a) \right] \\ &+ 0.00000012 \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

if \(h = 0.01\):

\[ \begin{split} R_{0.01}(f) &= \int_{a}^{b} f(x)\,dx - 0.0000042 \left[ f^\prime(b) - f^\prime(a) \right] \\ &+ 0.000000000012 \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

and if \(h = 10\):

\[ \begin{split} R_{10}(f) &= \int_{a}^{b} f(x)\,dx - 4.1667 \left[ f^\prime(b) - f^\prime(a) \right] \\ &+ 12.15 \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

Truncation error

If \(f \in C^2[a, b]\), applying the mean-value theorem to the leading term of the Euler-Maclaurin expansion gives the standard truncation error:

\[ R_h(f) - \int_{a}^{b} f(x)\,dx = -\frac{h^2}{24}(b-a)\,f^{\prime\prime}(c) \]

for some \(c \in [a, b]\). A corollary is that if \(f^{\prime\prime}(x) = 0\) on \([a, b]\) — i.e. if \(f\) is linear — then the rule is exact. If \(f\) is linear, \(n\) may be chosen to be 1.

Computation

The \(n\) nodes are the midpoints of the \(n\) subintervals:

\[ x_i = a + \left(i + \tfrac{1}{2}\right)h, \quad i = 0, 1, \ldots, n-1, \quad h = \frac{b-a}{n} \]

Each node carries weight \(h\), so the sum is:

\[ R_h(f) = h \sum_{i=0}^{n-1} f(x_i) \]

The implementation makes a single forward pass over the \(n\) midpoints, accumulating the weighted sum, then multiplies by \(h\). No endpoint evaluations are needed.

Limitations

The rectangle rule converges as \(O(h^2)\), making it the least accurate of the Newton-Cotes rules in this library. Prefer Simpson’s rule or Romberg’s method for smooth integrands. The rule is not suitable for functions with discontinuities or singularities inside \([a, b]\) without splitting the interval first.

Trapezoidal rule

Motivation

The trapezoidal rule approximates a definite integral by replacing the integrand with straight-line segments between adjacent sample points. It is a practical first choice for smooth functions and is the natural building block for Richardson extrapolation methods such as Romberg’s method.

Example

use integrate::newton_cotes::trapezoidal_rule;

let square = |x: f64| x * x;

let a = 0.0;
let b = 1.0;

let num_steps: usize = 1_000_000;

let integral = trapezoidal_rule(square, a, b, num_steps);
println!("{}", integral);

Understanding the Trapezoidal rule

Composite formula

The trapezoidal rule approximates the integral of \(f(x)\) on \([a, a+h]\) by the signed area of the trapezoid formed by the line connecting \((a, f(a))\) to \((a+h, f(a+h))\).

The composite rule decomposes \([a, b]\) into \(n\) subintervals of equal length \(h = (b-a)/n\) and sums the result over each subinterval:

\[ T_h(f) = h \left[ \frac{f(a)}{2} + f(a+h) + f(a+2h) + \cdots + f(b-h) + \frac{f(b)}{2} \right] \]

Euler-Maclaurin expansion

The Euler-Maclaurin summation formula relates \(T_h(f)\) to the exact integral:

\[ \begin{align} T_h(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^2}{12} \left[ f^\prime(b) - f^\prime(a) \right] - \frac{h^4}{720} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &+ \cdots + K h^{2p-2} \left[ f^{(2p-3)}(b) - f^{(2p-3)}(a) \right] + O(h^{2p}) \end{align} \]

The \(O(h^{2p})\) remainder means that for an infinitely differentiable function whose first \(2p-3\) derivatives vanish at both endpoints, the theorem guarantees only that

\[ \lim_{h \to 0} \left| \frac{T_h(f) - \int_{a}^{b} f(x),dx}{h^{2p}} \right| < M \]

for some \(M > 0\), not that the rule is exact.

The expansion also shows that \(n\) should be chosen so that \(h = (b-a)/n < 1\). For example, if \(h = 0.1\):

\[ \begin{split} T_{0.1}(f) &= \int_{a}^{b} f(x)\,dx + 0.00083 \left[ f^\prime(b) - f^\prime(a) \right] \\ &- 0.00000014 \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

if \(h = 0.01\):

\[ \begin{split} T_{0.01}(f) &= \int_{a}^{b} f(x)\,dx + 0.0000083 \left[ f^\prime(b) - f^\prime(a) \right] \\ &- 0.000000000014 \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

and if \(h = 10\):

\[ \begin{split} T_{10}(f) &= \int_{a}^{b} f(x)\,dx + 8.3333 \left[ f^\prime(b) - f^\prime(a) \right] \\ &- 13.89 \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

Truncation error

If \(f \in C^2[a, b]\), applying the mean-value theorem to the leading term of the Euler-Maclaurin expansion gives the standard truncation error:

\[ T_h(f) - \int_{a}^{b} f(x)\,dx = \frac{h^2}{12}(b-a)\,f^{\prime\prime}(c) \]

for some \(c \in [a, b]\). The positive sign reflects that the trapezoid overestimates the integral when \(f\) is convex (\(f^{\prime\prime} > 0\)) and underestimates when \(f\) is concave (\(f^{\prime\prime} < 0\)).

A corollary is that if \(f^{\prime\prime}(x) = 0\) on \([a, b]\) — i.e. if \(f\) is linear — then the rule is exact. If \(f\) is linear, \(n\) may be chosen to be 1.

Computation

The \(n+1\) nodes are the equally spaced points

\[ x_i = a + i,h, \quad i = 0, 1, \ldots, n, \quad h = \frac{b-a}{n} \]

The endpoint nodes carry weight \(h/2\) and all interior nodes carry weight \(h\), so the formula can be reorganised as

\[ T_h(f) = h,\left(\frac{f(a) + f(b)}{2} + \sum_{i=1}^{n-1} f(x_i)\right) \]

The implementation accumulates the interior sum in a single forward pass over \(i = 1, \ldots, n-1\), adds the half-weighted endpoints, then multiplies by \(h\). A total of \(n+1\) function evaluations are required.

Limitations

The trapezoidal rule converges as \(O(h^2)\) for general smooth functions. While more accurate than the rectangle rule, it is outperformed by Simpson’s rule for smooth integrands. It should not be used for functions with discontinuities or singularities in \([a, b]\) without splitting the interval first. For smooth periodic functions on \([a, b]\) the rule achieves spectral convergence (faster than any power of \(h\)) because the endpoint correction terms in the Euler-Maclaurin expansion all vanish.

Simpson’s rule

Motivation

Simpson’s rule fits a quadratic polynomial through three equally-spaced points per subinterval and integrates it exactly, yielding fourth-order convergence for smooth functions. It is one of the most widely used Newton-Cotes formulas and is exact for polynomials of degree three or less, despite using only second-degree interpolants — a consequence of the symmetric placement of nodes.

Example

use integrate::newton_cotes::simpson_rule;

let square = |x: f64| x * x;

let a = 0.0;
let b = 1.0;

let num_steps: usize = 1_000_000;

let integral = simpson_rule(square, a, b, num_steps);
println!("{}", integral);

Understanding Simpson’s rule

Composite formula

Simpson’s rule approximates the integral of \(f(x)\) on \([a, a+h]\) by integrating the quadratic through \((a, f(a))\), \((a+h/2, f(a+h/2))\), and \((a+h, f(a+h))\).

The composite rule decomposes \([a, b]\) into \(n\) subintervals of equal length \(h = (b-a)/n\) and sums over each subinterval:

\[ \begin{align} S_h(f) &= \frac{h}{6} \left[ f(a) + 4f\left(a+\tfrac{h}{2}\right) + 2f(a+h) \right. \\ & \left. + \cdots + 2f(b-h) + 4f\left(b-\tfrac{h}{2}\right) + f(b) \right] \end{align} \]

Euler-Maclaurin expansion

An immediate consequence of the Euler-Maclaurin summation formula relates \(S_h(f)\) to the exact integral:

\[ \begin{align} S_h(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^4}{2880} \left[ f^{(3)}(b) - f^{(3)}(a) \right] - \frac{h^6}{96768} \left[ f^{(5)}(b) - f^{(5)}(a) \right] \\ &+ \cdots + K h^{2p-2} \left[ f^{(2p-3)}(b) - f^{(2p-3)}(a) \right] + O(h^{2p}) \end{align} \]

The \(O(h^{2p})\) remainder means that for an infinitely differentiable function whose first \(2p-3\) derivatives vanish at both endpoints, the theorem guarantees only that

\[ \lim_{h \to 0} \left| \frac{S_h(f) - \int_{a}^{b} f(x),dx}{h^{2p}} \right| < M \]

for some \(M > 0\), not that the rule is exact.

The expansion also shows that \(n\) should be chosen so that \(h = (b-a)/n < 1\). For example, if \(h = 0.1\):

\[ \begin{align} S_{0.1}(f) &= \int_{a}^{b} f(x)\,dx + 3.5 \cdot 10^{-8} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- 1.033 \cdot 10^{-11} \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{align} \]

if \(h = 0.01\):

\[ \begin{align} S_{0.01}(f) &= \int_{a}^{b} f(x)\,dx + 3.5 \cdot 10^{-12} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- 1.033 \cdot 10^{-17} \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{align} \]

and if \(h = 10\):

\[ \begin{align} S_{10}(f) &= \int_{a}^{b} f(x)\,dx + 3.47 \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- 10.33 \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{align} \]

Truncation error

If \(f \in C^4[a, b]\), applying the mean-value theorem to the leading term of the Euler-Maclaurin expansion gives the standard truncation error:

\[ S_h(f) - \int_{a}^{b} f(x)\,dx = \frac{h^4}{2880}(b-a)\,f^{(4)}(c) \]

for some \(c \in [a, b]\). A corollary is that if \(f^{(4)}(x) = 0\) on \([a, b]\) — i.e. if \(f\) is a polynomial of degree at most 3 — then the rule is exact. If \(f\) is a cubic, \(n\) may be chosen to be 1.

Computation

The \(2n+1\) nodes are the endpoints, the \(n\) subinterval midpoints, and the \(n-1\) interior shared endpoints:

\[ x_i = a + \tfrac{i}{2}\,h, \quad i = 0, 1, \ldots, 2n, \quad h = \frac{b-a}{n} \]

The effective weights are:

• \(h/6\) at the endpoints \(x_0 = a\) and \(x_{2n} = b\)
• \(2h/3\) at each midpoint \(x_1, x_3, \ldots, x_{2n-1}\)
• \(h/3\) at each interior shared node \(x_2, x_4, \ldots, x_{2n-2}\)

The implementation makes a single forward pass over the \(n\) subintervals, accumulating three terms per step — \(f(a) + 4f(a+h/2)\) for the first, \(2f(a+ih) + 4f(a+ih+h/2)\) for interior subintervals, and \(f(b)\) for the last — then multiplies the total by \(h/6\). A total of \(2n+1\) function evaluations are required.

Limitations

Simpson’s rule converges as \(O(h^4)\) and requires at least four-times-differentiable integrands for its error estimate to apply. It is not suitable for functions with discontinuities or singularities without splitting the interval. The rule evaluates \(2n+1\) points per \(n\) subintervals — slightly more than the trapezoidal rule’s \(n+1\) — so Newton’s 3/8 rule may be preferred when the number of subintervals is a multiple of three and an equivalent accuracy is sought with more uniform node spacing.

Newton’s 3/8 rule

Motivation

Newton’s 3/8 rule fits a cubic polynomial through four equally-spaced points per subinterval and integrates it exactly. It achieves the same order of accuracy as Simpson’s rule (\(O(h^4)\)) but uses three equally-spaced subdivisions per panel, making it useful when the number of subintervals is a multiple of three or when more uniform node spacing is preferred.

Example

use integrate::newton_cotes::newton_rule;

let square = |x: f64| x * x;

let a = 0.0;
let b = 1.0;

let num_steps: usize = 1_000_000;

let integral = newton_rule(square, a, b, num_steps);
println!("{}", integral);

Understanding Newton’s 3/8 rule

Composite formula

Newton’s \(3/8\) rule approximates the integral of \(f(x)\) on \([a, a+h]\) by integrating the cubic through \((a, f(a))\), \((a+h/3, f(a+h/3))\), \((a+2h/3, f(a+2h/3))\), and \((a+h, f(a+h))\).

The composite rule decomposes \([a, b]\) into \(n\) subintervals of equal length \(h = (b-a)/n\) and sums over each subinterval:

\[ \begin{align} N_h(f) &= \frac{h}{8} \left[ f(a) + 3f\left(a+\tfrac{h}{3}\right) + 3f\left(a+\tfrac{2h}{3}\right) + 2f(a+h) \right. \\ & \left. + \cdots + 2f(b-h) + 3f\left(b-\tfrac{2h}{3}\right) + 3f\left(b-\tfrac{h}{3}\right) + f(b) \right] \end{align} \]

Euler-Maclaurin expansion

An immediate consequence of the Euler-Maclaurin summation formula relates \(N_h(f)\) to the exact integral:

\[ \begin{align} N_h(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^4}{6480} \left[ f^{(3)}(b) - f^{(3)}(a) \right] - \frac{h^6}{244944} \left[ f^{(5)}(b) - f^{(5)}(a) \right] \\ &+ \cdots + K h^{2p-2} \left[ f^{(2p-3)}(b) - f^{(2p-3)}(a) \right] + O(h^{2p}) \end{align} \]

The \(O(h^{2p})\) remainder means that for an infinitely differentiable function whose first \(2p-3\) derivatives vanish at both endpoints, the theorem guarantees only that

\[ \lim_{h \to 0} \left| \frac{N_h(f) - \int_{a}^{b} f(x),dx}{h^{2p}} \right| < M \]

for some \(M > 0\), not that the rule is exact.

The expansion also shows that \(n\) should be chosen so that \(h = (b-a)/n < 1\). For example, if \(h = 0.1\):

\[ \begin{align} N_{0.1}(f) &= \int_{a}^{b} f(x)\,dx + 1.5 \cdot 10^{-8} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- 4.1 \cdot 10^{-12} \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{align} \]

if \(h = 0.01\):

\[ \begin{align} N_{0.01}(f) &= \int_{a}^{b} f(x)\,dx + 1.5 \cdot 10^{-12} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- 4.1 \cdot 10^{-18} \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{align} \]

and if \(h = 10\):

\[ \begin{align} N_{10}(f) &= \int_{a}^{b} f(x)\,dx + 1.54 \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- 4.08 \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{align} \]

Truncation error

If \(f \in C^4[a, b]\), applying the mean-value theorem to the leading term of the Euler-Maclaurin expansion gives the standard truncation error:

\[ N_h(f) - \int_{a}^{b} f(x)\,dx = \frac{h^4}{6480}(b-a)\,f^{(4)}(c) \]

for some \(c \in [a, b]\). A corollary is that if \(f^{(4)}(x) = 0\) on \([a, b]\) — i.e. if \(f\) is a polynomial of degree at most 3 — then the rule is exact. If \(f\) is a cubic, \(n\) may be chosen to be 1.

Computation

The \(3n+1\) nodes are the endpoints and the two third-points of each subinterval:

\[ x_{3i+k} = a + \left(i + \tfrac{k}{3}\right)h, \quad i = 0, \ldots, n-1,\quad k = 0, 1, 2, \quad \text{plus } x_{3n} = b \]

The effective weights are:

• \(h/8\) at the endpoints \(x_0 = a\) and \(x_{3n} = b\)
• \(3h/8\) at each third-point \(x_1, x_2, x_4, x_5, \ldots\)
• \(h/4\) at each interior shared node \(x_3, x_6, \ldots, x_{3n-3}\)

The implementation makes a single forward pass over the \(n\) subintervals, accumulating four terms per step — \(f(a) + 3f(a+h/3) + 3f(a+2h/3)\) for the first, \(2f(a+ih) + 3f(a+ih+h/3) + 3f(a+ih+2h/3)\) for interior subintervals, and \(f(b)\) for the last — then multiplies the total by \(h/8\). A total of \(3n+1\) function evaluations are required.

Limitations

Newton’s 3/8 rule converges as \(O(h^4)\), the same order as Simpson’s rule, but requires \(3n+1\) evaluations versus Simpson’s \(2n+1\) for the same \(n\). For general use, Simpson’s rule is more efficient. Like all fixed-step Newton-Cotes formulas, it is not suitable for integrands with discontinuities or singularities inside \([a, b]\) without splitting the interval first.

Gaussian quadrature

Instead of evaluating the integrand at equally spaced nodes as in Newton-Cotes methods, Gaussian quadrature methods make a judicious choice of nodes so as to maximize the precision of the numerical integration relative to the number of integrand evaluations.

The common Gaussian quadrature methods are:

• Gauss-Legendre used to integrate a function \(f(x)\) over a closed and bounded interval \([a,b]\).
• Gauss-Laguerre used to integrate a function of the form \(f(x) e^{-x}\) over the positive x-axis \(\lbrace x \in \mathbb{R} : x > 0 \rbrace\).
• Gauss-Hermite used to integrate a function of the form \(f(x) e^{-x^2}\) over the entire x-axis, \(\lbrace x \in \mathbb{R} : -\infty < x < \infty \rbrace\).
• Gauss-Chebyshev First Kind used to integrate a function of the form \(\frac{f(x)}{\sqrt( 1-x^2 )}\) over the interval \([-1,1]\).
• Gauss-Chebyshev Second Kind used to integrate a function of the form \(f(x) * \sqrt{ 1-x^2 }\) over the interval \([-1,1]\).

Gauss-Legendre

Motivation

Instead of evaluating the integrand at equally spaced nodes as Newton-Cotes methods do, Gauss-Legendre quadrature places nodes at the zeros of Legendre polynomials — positions that are optimal in the sense that an \(n\)-point rule integrates all polynomials of degree up to \(2n-1\) exactly. This means it achieves the same accuracy as a \(2n\)-point trapezoidal rule with only half the function evaluations, making it the preferred method for smooth integrands on finite intervals.

Example

use integrate::gauss_quadrature::legendre_rule;

let square = |x: f64| x * x;

let a = 0.0;
let b = 1.0;

let n: usize = 100;

let integral = legendre_rule(square, a, b, n);
println!("{}", integral);

Understanding Gauss-Legendre rule

Legendre polynomials

With respect to the inner product

\[ \langle f, g \rangle = \int_{-1}^{1} f(x) \cdot g(x) \cdot w(x) \, dx \]

the Legendre polynomials, defined by Rodrigues’ formula

\[ P_n(x) = \frac{1}{2^n n!} \frac{\partial^{n} (x^2 - 1)^n}{\partial x^n} \quad \text{for} \quad n > 0 \]

and \(P_0(x) = 1\), form an orthogonal family with weight function \(w(x) = 1\) on \([-1, 1]\).

Quadrature formula on \([-1, 1]\)

The \(n\)-point Gauss-Legendre quadrature formula \(GL_n(f)\), for approximating \(\int_{-1}^{1} f(x)\,dx\), is given by

\[ GL_n(f) = A_1 f(x_1) + \cdots + A_n f(x_n) \]

where \(x_1, \ldots, x_n\) are the zeros of \(P_n\) and the weights are

\[ A_i = \frac{2}{1 - x_i^2} \cdot \frac{1}{n^2 P_{n-1}(x_i)^2} \quad \text{for} \quad i = 1, \ldots, n \]

All weights \(A_i\) are positive, which makes the rule numerically stable.

Change of variables to \([a, b]\)

To integrate over a general interval \([a, b]\), the substitution \(t = \dfrac{2(x - a)}{b - a} - 1\) gives

\[ \int_{a}^{b} f(x) \, dx = \frac{b - a}{2} \int_{-1}^{1} f\left( \frac{t(b-a) + (b+a)}{2} \right) dt \]

so the \(n\)-point Gauss-Legendre approximation on \([a, b]\) is

\[ GL_n(f, a, b) = A_1’ f(x_1’) + \cdots + A_n’ f(x_n’) \]

where the mapped nodes and weights are

\[ x_i’ = \frac{b-a}{2} \cdot x_i + \frac{b+a}{2}, \qquad A_i’ = \frac{b-a}{2} \cdot A_i \quad \text{for} \quad i = 1, \ldots, n \]

Truncation error

The truncation error for the rule on \([-1, 1]\) is

\[ \int_{-1}^{1} f(x) \, dx - GL_n(f) = K \cdot \frac{f^{(2n)}(c)}{(2n)!} \]

where \(c \in (-1, 1)\) is unknown and \(K\) is a constant determined by

\[ K = \int_{-1}^{1} x^{2n} \, dx - GL_n(x^{2n}) \]

The same form holds on \([a, b]\) with \(c \in (a, b)\). A corollary is that if \(f^{(2n)}(x) = 0\) for all \(x \in [a, b]\) — i.e. if \(f\) is a polynomial of degree at most \(2n - 1\) — then the rule is exact.

Node and weight computation

Most Gauss quadrature rules compute nodes and weights by building a tridiagonal Jacobi matrix and finding its eigenvalues (Golub-Welsch algorithm), which costs \(O(n^2)\) in total. The Legendre rule in this crate uses a different strategy that delivers each node and weight independently in \(O(1)\) time, regardless of \(n\):

• \(n \leq 100\): each \((x_i, A_i)\) pair is read directly from a precomputed table of \(\theta\)-coordinates and weights stored in the source. The lookup is a single array index with no arithmetic at all.

• \(n > 100\): Bogaert’s asymptotic formula expresses each node as \(x_i = \cos\theta_i\), where \(\theta_i\) is obtained from the \(k\)-th zero of the Bessel function \(J_0\) via a fixed-degree Chebyshev polynomial correction (6 terms for the node, 8 terms for the weight). The entire computation is a handful of fused-multiply-add operations — no Newton iteration, no eigenvalue solver.

The result is that the per-node cost is \(O(1)\) and the error on every node and weight is within a few ulps of the true value, even for very large \(n\).

Limitations

Gauss-Legendre quadrature is best suited for smooth functions on finite intervals. It is not appropriate for functions with singularities, discontinuities, or functions defined on unbounded domains. Rapidly oscillating functions may also require a very large number of nodes \(n\). For semi-infinite or doubly-infinite domains use the Gauss-Laguerre or Gauss-Hermite rules instead.

References

Node and weight computation based on the original C++ implementation by Ignace Bogaert, featuring \(O(1)\) per-node complexity and errors within a few ulps:

\[ [1]: \text{Bogaert I., } \textit{Iteration-free computation of Gauss-Legendre quadrature nodes and weights,} \\ \text{SIAM Journal on Scientific Computing, Volume 36, Number 3, 2014, pages A1008\text{–}A1026.} \]

For more details, here is a link to the article.

Gauss-Chebyshev Quadrature

Motivation

Chebyshev quadrature is designed for integrands that carry an algebraic singularity at the endpoints of \([-1, 1]\). Rather than struggling with the singularity numerically, the weight function absorbs it exactly into the orthogonality relation, leaving a smooth \(f(x)\) to be sampled.

Two weight functions are supported:

• First kind — \(w(x) = 1/\sqrt{1-x^2}\), which diverges at \(\pm 1\) and arises in arclength integrals, potential theory, and spectral methods.
• Second kind — \(w(x) = \sqrt{1-x^2}\), which vanishes at \(\pm 1\) and appears in problems with elliptic geometry or Jacobi-weight boundary conditions.

Because the zeros and weights of both families have exact trigonometric closed forms, node and weight computation requires only a single trigonometric evaluation per point — making this the most computationally efficient Gauss rule in the library.

Chebyshev Polynomials of the First Kind

Example

use integrate::gauss_quadrature::gauss_first_kind_chebyshev_rule;

let f = |x: f64| 1.0;

let n: usize = 100;

let integral = gauss_first_kind_chebyshev_rule(f, n);
println!("{}", integral);

Understanding Gauss-Chebyshev rule (first kind)

With respect to the inner product

\[ \langle f, g \rangle = \int_{-1}^{1} f(x) \cdot g(x) \cdot w(x) \, dx \]

the Chebyshev polynomials of the first kind, defined by

\[ T_n(x) = \cos(n \arccos(x)) \quad \text{for} \quad n > 0 \]

and \(T_0(x) = 1\), form an orthogonal family with weight function \(w(x) = 1/\sqrt{1-x^2}\) on \([-1, 1]\), satisfying the three-term recurrence

\[ T_{k+1}(x) = 2x\,T_k(x) - T_{k-1}(x) \]

The \(n\)-point Gauss-Chebyshev quadrature formula of the first kind, \(GC_n^{(1)}(f)\), approximating \(\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}}\,dx\), is

\[ GC_n^{(1)}(f) = A_1 f(x_1) + \cdots + A_n f(x_n) \]

where \(x_1, \ldots, x_n\) are the zeros of \(T_n\) and all weights are equal:

\[ A_i = \frac{\pi}{n} \quad \text{for} \quad i = 1, \ldots, n \]

Truncation error

\[ \begin{split} \int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}}\,dx - GC_n^{(1)}(f) \\ = \frac{\pi}{2^{2n-1}\,(2n)!}\, f^{(2n)}(\xi) \end{split} \]

for some \(\xi \in (-1, 1)\). The rule is exact for all polynomials of degree at most \(2n - 1\).

Node and weight computation

The zeros of \(T_n\) have the explicit closed form

\[ x_k = \cos\left(\frac{(2k-1),\pi}{2n}\right), \quad k = 1, \ldots, n \]

and all weights are the same constant \(\pi/n\). Each node requires a single cosine evaluation; no weight computation is needed beyond a single division. There is no Jacobi matrix, no eigenvalue solver, and no iterative refinement — the total cost is \(O(n)\) with a very small constant.

Chebyshev Polynomials of the Second Kind

Example

use integrate::gauss_quadrature::gauss_second_kind_chebyshev_rule;

let f = |x: f64| 1.0;

let n: usize = 100;

let integral = gauss_second_kind_chebyshev_rule(f, n);
println!("{}", integral);

Understanding Gauss-Chebyshev rule (second kind)

With respect to the same inner product structure, the Chebyshev polynomials of the second kind, defined by

\[ U_n(x) = \frac{\sin((n+1)\arccos(x))}{\sin(\arccos(x))} \quad \text{for} \quad n > 0 \]

and \(U_0(x) = 1\), form an orthogonal family with weight function \(w(x) = \sqrt{1-x^2}\) on \([-1, 1]\), satisfying the same three-term recurrence

\[ U_{k+1}(x) = 2x\,U_k(x) - U_{k-1}(x) \]

The \(n\)-point Gauss-Chebyshev quadrature formula of the second kind, \(GC_n^{(2)}(f)\), approximating \(\int_{-1}^{1} f(x)\sqrt{1-x^2}\,dx\), is

\[ GC_n^{(2)}(f) = A_1 f(x_1) + \cdots + A_n f(x_n) \]

where \(x_1, \ldots, x_n\) are the zeros of \(U_n\) and the weights are

\[ A_i = \frac{\pi}{n+1} \sin^2\left(\frac{i,\pi}{n+1}\right) \quad \text{for} \quad i = 1, \ldots, n \]

Truncation error

\[ \begin{split} \int_{-1}^{1} f(x)\sqrt{1-x^2}\,dx - GC_n^{(2)}(f) \\ = \frac{\pi}{2^{2n+1}\,(2n)!}\, f^{(2n)}(\xi) \end{split} \]

for some \(\xi \in (-1, 1)\). The rule is exact for all polynomials of degree at most \(2n - 1\).

Node and weight computation

The zeros of \(U_n\) have the explicit closed form

\[ x_k = \cos\left(\frac{k,\pi}{n+1}\right), \quad k = 1, \ldots, n \]

and each weight requires one additional sine evaluation:

\[ A_k = \frac{\pi}{n+1} \sin^2\left(\frac{k,\pi}{n+1}\right) \]

As with the first kind, there is no Jacobi matrix, no eigenvalue solver, and no iterative refinement. The total cost is \(O(n)\) with a very small constant.

The equal-angle spacing of \(k\pi/(n+1)\) is the same grid used in the discrete sine transform (DST), so the second-kind nodes and weights can also be interpreted as a DST-I evaluation.

Limitations

Gauss-Chebyshev quadrature is only suitable when the integrand can be expressed in the weighted form \(f(x)/\sqrt{1-x^2}\) (first kind) or \(f(x)\sqrt{1-x^2}\) (second kind) on \([-1, 1]\). Applying either rule to an unweighted integrand will produce inaccurate results, since the accuracy guarantee relies on orthogonality with respect to the corresponding weight function. Both rules are restricted to the interval \([-1, 1]\); for other domains or weight functions use a different Gauss rule.

Gauss-Laguerre

Motivation

Gauss-Laguerre quadrature is designed for integrals of the form \(\int_0^{+\infty} f(x)\,e^{-x}\,dx\), where the exponential decay naturally arises in probability distributions, physics, and engineering problems. By using the zeros and weights of Laguerre polynomials, the method achieves high accuracy with far fewer function evaluations than general-purpose rules on a truncated domain.

Example

use integrate::gauss_quadrature::gauss_laguerre_rule;

let f = |x: f64| 1.0;

let n: usize = 100;

let integral = gauss_laguerre_rule(f, n);
println!("{}", integral);

Understanding Gauss-Laguerre rule

Laguerre polynomials

With respect to the inner product

\[ \langle f, g \rangle = \int_{0}^{+\infty} f(x) \cdot g(x) \cdot w(x) \, dx \]

the Laguerre polynomials, defined by Rodrigues’ formula

\[ L_n(x) = e^x \frac{\partial^{n} (x^n e^{-x})}{\partial x^n} \quad \text{for} \quad n > 0 \]

and \(L_0(x) = 1\), form an orthogonal family with weight function \(w(x) = e^{-x}\) on the positive \(x\)-axis, satisfying the three-term recurrence

\[ L_{k+1}(x) = \frac{(2k + 1 - x)\,L_k(x) - k\,L_{k-1}(x)}{k + 1} \]

Quadrature formula

The \(n\)-point Gauss-Laguerre formula \(GL_n(f)\), approximating \(\int_0^{+\infty} f(x)\,e^{-x}\,dx\), is

\[ GL_n(f) = A_1 f(x_1) + \cdots + A_n f(x_n) \]

where \(x_1, \ldots, x_n\) are the zeros of \(L_n\) and the weights are

\[ A_i = \frac{x_i}{(n+1)^2 \, L_{n+1}(x_i)^2} \quad \text{for} \quad i = 1, \ldots, n \]

Truncation error

The truncation error for the \(n\)-point rule is

\[ \int_0^{+\infty} f(x)\,e^{-x}\,dx - GL_n(f) = \frac{(n!)^2}{(2n)!}\, f^{(2n)}(\xi) \]

for some \(\xi > 0\). A corollary is that if \(f\) is a polynomial of degree at most \(2n - 1\) then the rule is exact.

Node and weight computation

Nodes and weights are computed via the Golub-Welsch algorithm: the \(n\) zeros of \(L_n\) are the eigenvalues of the symmetric tridiagonal Jacobi matrix

\[ J = \begin{pmatrix} d_0 & e_1 & & \\ e_1 & d_1 & e_2 & \\ & \ddots & \ddots & e_{n-1} \\ & & e_{n-1} & d_{n-1} \end{pmatrix} \]

with diagonal entries \(d_k = 2k + 1\) and off-diagonal entries \(e_k = k\) for \(k = 0, 1, \ldots, n-1\).

Eigenvalues are located by bisection on Sturm sequences, costing \(O(n)\) per node and \(O(n^2)\) in total. Once each zero \(x_i\) is known, the three-term recurrence is used to evaluate \(L_{n+1}(x_i)\) and compute the weight \(A_i\) directly.

For large \(n\) the weights of the outermost nodes decrease extremely rapidly and may underflow to zero; a warning is printed to stderr when this occurs.

Limitations

Gauss-Laguerre quadrature is only appropriate for integrals of the form \(\int_0^{+\infty} f(x)\,e^{-x}\,dx\). If the integrand does not carry the exponential factor \(e^{-x}\), the method will produce inaccurate results. It is not suitable for integrands defined on finite intervals or on the full real line.

Gauss-Hermite

Motivation

Gauss-Hermite quadrature is designed for integrals of the form \(\int_{-\infty}^{+\infty} f(x)\,e^{-x^2}\,dx\), which appear naturally in quantum mechanics, probability theory (Gaussian distributions), and signal processing. By exploiting the orthogonality of Hermite polynomials, it provides exponential convergence for smooth integrands with Gaussian-type decay.

Example

use integrate::gauss_quadrature::gauss_hermite_rule;

let f = |x: f64| 1.0;

let n: usize = 100;

let integral = gauss_hermite_rule(f, n);
println!("{}", integral);

Understanding Gauss-Hermite rule

Hermite polynomials

With respect to the inner product

\[ \langle f, g \rangle = \int_{-\infty}^{+\infty} f(x) \cdot g(x) \cdot w(x) \, dx \]

the Hermite polynomials, defined by Rodrigues’ formula

\[ H_n(x) = (-1)^n \cdot e^{x^2} \cdot \frac{\partial^{n} e^{-x^2}}{\partial x^n} \quad \text{for} \quad n > 0 \]

and \(H_0(x) = 1\), form an orthogonal family with weight function \(w(x) = e^{-x^2}\) on the entire real line, satisfying the three-term recurrence

\[ H_{k+1}(x) = 2x\,H_k(x) - 2k\,H_{k-1}(x) \]

Because \(H_n\) is either even or odd, its zeros are symmetric about the origin: if \(x\) is a zero then so is \(-x\), and for odd \(n\) there is always a zero at \(x = 0\).

Quadrature formula

The \(n\)-point Gauss-Hermite formula \(GH_n(f)\), approximating \(\int_{-\infty}^{+\infty} f(x)\,e^{-x^2}\,dx\), is

\[ GH_n(f) = A_1 f(x_1) + \cdots + A_n f(x_n) \]

where \(x_1, \ldots, x_n\) are the zeros of \(H_n\) and the weights are

\[ A_i = \frac{2^{n+1} \cdot n! \cdot \sqrt{\pi}}{H_{n-1}(x_i)^2} \quad \text{for} \quad i = 1, \ldots, n \]

Truncation error

The truncation error for the \(n\)-point rule is

\[ \begin{split} \int_{-\infty}^{+\infty} f(x)\,e^{-x^2}\,dx - GH_n(f) \\ = \frac{n!\,\sqrt{\pi}}{2^n\,(2n)!}\, f^{(2n)}(\xi) \end{split} \]

for some \(\xi \in \mathbb{R}\). A corollary is that if \(f\) is a polynomial of degree at most \(2n - 1\) then the rule is exact.

Node and weight computation

Nodes and weights are computed via the Golub-Welsch algorithm: the \(n\) zeros of \(H_n\) are the eigenvalues of the symmetric tridiagonal Jacobi matrix with all-zero diagonal and off-diagonal entries

\[ e_k = \sqrt{\frac{k}{2}}, \quad k = 1, \ldots, n-1 \]

The all-zero diagonal directly reflects the symmetry of \(H_n\): eigenvalues come in \(\pm\) pairs (with \(0\) as an unpaired eigenvalue for odd \(n\)), halving the effective search space. Eigenvalues are located by bisection on Sturm sequences, costing \(O(n)\) per node and \(O(n^2)\) in total.

Once each zero \(x_i\) is known, \(H_{n-1}(x_i)\) is evaluated via the three-term recurrence and the weight \(A_i\) computed directly. For large \(n\) the factor \(2^{n+1} \cdot n!\) overflows f64, so the implementation switches to BigInt arithmetic for the numerator and denominator before converting back. A warning is printed to stderr if any weight underflows to zero after the conversion.

Limitations

Gauss-Hermite quadrature is only appropriate for integrals of the form \(\int_{-\infty}^{+\infty} f(x)\,e^{-x^2}\,dx\). If the integrand does not carry the Gaussian factor \(e^{-x^2}\), the method will produce inaccurate results. It is not suitable for integrands on finite intervals or those with non-Gaussian tail behavior. For large \(n\) the outermost weights are extremely small and may underflow, effectively reducing the usable order of the rule.

Gauss-Kronrod

Motivation

The \((2N+1)\)-point Gauss-Kronrod rule simultaneously produces a high-accuracy integral estimate and a reliable error bound in a single pass over \(2N+1\) function evaluations. It does so by augmenting an \(N\)-point Gauss-Legendre rule with \(N+1\) optimally placed Kronrod points: the original \(N\) Gauss nodes are reused, so no function evaluations are wasted.

This makes it the method of choice when both accuracy and a trustworthy error estimate are needed. It is the engine behind QUADPACK and SciPy’s scipy.integrate.quad.

Example

use integrate::gauss_kronrod::gauss_kronrod_rule;

let f = |x: f64| x.exp();

let a = 0.0_f64;
let b = 1.0_f64;
let n = 7_usize; // 15-point Gauss-Kronrod rule

match gauss_kronrod_rule(f, a, b, n) {
    Ok((integral, error)) => {
        println!("Integral:       {:.10}", integral);
        println!("Error estimate: {:.2e}", error);
    }
    Err(e) => println!("Error: {}", e),
}

Understanding Gauss-Kronrod rule

The Gauss-Kronrod rule is constructed in two stages.

Stage 1: Gauss-Legendre base rule

The \(N\)-point Gauss-Legendre rule approximates \(\int_{-1}^{1} f(x)\,dx\) by

\[ G_N(f) = \sum_{i=1}^{N} w_i^G f(x_i) \]

where \(x_1, \ldots, x_N\) are the zeros of the Legendre polynomial \(P_N\) and \(w_1^G, \ldots, w_N^G\) are positive weights chosen so that the formula is exact for all polynomials of degree at most \(2N - 1\).

Stage 2: Kronrod extension

The Stieltjes polynomial \(E_{N+1}\) of degree \(N+1\) is the unique monic polynomial satisfying

\[ \int_{-1}^{1} E_{N+1}(x) \cdot x^k \cdot P_N(x) \, dx = 0, \quad k = 0, 1, \ldots, N \]

Its \(N+1\) zeros \(t_1, \ldots, t_{N+1}\) interlace strictly with those of \(P_N\) and all lie in \((-1, 1)\). Adding them to the \(N\) existing Gauss-Legendre nodes gives \(2N+1\) distinct abscissas. Weights for all \(2N+1\) points can then be chosen so that the combined formula

\[ GK_{2N+1}(f) = \sum_{i=1}^{N} w_i^{GK} f(x_i) + \sum_{j=1}^{N+1} v_j^{GK} f(t_j) \]

is exact for all polynomials of degree at most \(3N + 1\) (when \(N\) is odd). Since the Gauss nodes are reused, the total cost is \(2N+1\) function evaluations.

Change of variables to \([a, b]\)

For an integral over a general interval \([a, b]\), the substitution \(x = \dfrac{(b-a)\,t + (b+a)}{2}\) gives

\[ \int_a^b f(x) \, dx = \frac{b-a}{2} \int_{-1}^{1} f\left(\frac{(b-a)\,t + (b+a)}{2}\right) dt \]

so the \((2N+1)\)-point Gauss-Kronrod approximation on \([a, b]\) is

\[ GK_{2N+1}(f, a, b) = \frac{b-a}{2} \cdot GK_{2N+1}(\tilde{f}), \quad \tilde{f}(t) = f\left(\frac{(b-a)\,t + (b+a)}{2}\right) \]

Error estimate

gauss_kronrod_rule returns both the integral and an error estimate:

• integral \(= GK_{2N+1}(f, a, b)\)
• error \(= \left\lvert GK_{2N+1}(f, a, b) - G_N(f, a, b) \right\rvert\)

The error is the absolute difference between the \((2N+1)\)-point Kronrod estimate and the \(N\)-point Gauss-Legendre estimate computed on the same set of function values. For smooth integrands this quantity is a reliable indicator of the absolute integration error.

Choosing \(N\)

The default \(N = 7\) (15-point rule) matches the default used by QUADPACK and SciPy on each subinterval and gives double-precision accuracy for most smooth integrands.

Limitations

The single-interval Gauss-Kronrod rule requires the integrand to be smooth on \([a, b]\). For integrands with singularities, discontinuities, or sharp peaks the error estimate may be unreliable and the accuracy may be poor. In such cases the interval should be split manually or an adaptive method should be used instead.

The rule is defined only on finite intervals \([a, b]\). For semi-infinite or doubly-infinite domains use the Gauss-Laguerre or Gauss-Hermite rules instead.

References

Node and weight computation ported from the original Fortran 77 by Robert Piessens and Maria Branders (C translation by John Burkardt, GNU LGPL):

\[ [1]: \text{Piessens R., Branders M., } \textit{A Note on the Optimal Addition of Abscissas to Quadrature Formulas of Gauss and Lobatto,} \\ \text{Mathematics of Computation, Vol. 28, No. 125, 1974, pp. 135\text{–}139.} \]

Adaptive Methods

Adaptive methods are used to estimate the integral of a function on a closed and bounded interval. The basic idea is to take smaller subintervals in regions where the function changes significantly and larger subintervals in regions where the function doesn’t change significantly. What is significant depends on the quadrature technique involved, for Simpson’s rule, if the function is nearly cubic (or quadratic, or linear, or constant) then the function doesn’t change significantly from being cubic and a large interval can be used.

Adaptive Simpson

Motivation

Fixed-step quadrature methods allocate function evaluations uniformly across the interval, which is wasteful when the integrand is smooth almost everywhere but has rapid variation or near-singularities in a small region. Adaptive quadrature automatically concentrates evaluations where the integrand is hardest to approximate, achieving a target accuracy with fewer total evaluations than any fixed-step rule of the same order.

The adaptive Simpson method combines the well-understood accuracy of Simpson’s rule with a local error indicator derived from Richardson extrapolation. Subintervals where the integrand is smooth are accepted quickly; subintervals where it varies rapidly are bisected and retried until the local error falls below a pro-rated fraction of the global tolerance.

Example

use integrate::adaptive_quadrature::adaptive_simpson_method;

let f = |x: f64| x.exp();

let a = 0.0;
let b = 1.0;

let tolerance = 10.0e-6;
let min_h = 10.0e-3;

let result = adaptive_simpson_method(f, a, b, min_h, tolerance);

match result {
    Ok(res)  => println!("{}", res),
    Err(err) => println!("{}", err),
};

Understanding the Adaptive Simpson method

Error indicator

On a subinterval \([l, r]\) of length \(H = r - l\), two Simpson estimates are computed from five evenly-spaced points \(l,\; l+H/4,\; m,\; r-H/4,\; r\) (where \(m = (l+r)/2\)):

• Single estimate \(S_1\) — standard Simpson using three points:

\[ S_1 = \frac{H}{6}\bigl[f(l) + 4f(m) + f(r)\bigr] \]

• Composite estimate \(S_2\) — Simpson applied to the two half-panels \([l, m]\) and \([m, r]\):

\[ S_2 = \frac{H}{12}\bigl[f(l) + 4f(l+H/4) + 2f(m) + 4f(r-H/4) + f(r)\bigr] \]

By Richardson extrapolation, \(S_2\) is approximately 16 times more accurate than \(S_1\) for smooth integrands (halving the step size reduces the \(O(h^5)\) per-panel error by \(2^4 = 16\)). Consequently

\[ S_1 - \int_{l}^{r} f\,dx \approx \frac{16}{15}(S_1 - S_2) \qquad \text{and} \qquad \left|\int_{l}^{r} f\,dx - S_2\right| \approx \frac{1}{15}|S_1 - S_2| \]

The difference \(|S_1 - S_2|\) is therefore a reliable local error indicator for \(S_2\), and the true error in the accepted estimate is roughly 15 times smaller than the acceptance threshold — making the method conservative.

Acceptance criterion

The global tolerance is distributed across subintervals in proportion to their length. A subinterval \([l, r]\) is accepted when

\[ |S_1 - S_2| < \varepsilon_{\text{local}} = 2 \cdot \text{tol} \cdot \frac{r - l}{b - a} \]

where \(\text{tol}\) is the user-specified global tolerance and \(b - a\) is the full integration interval. This proportional allocation ensures that the sum of all local errors is bounded by the global tolerance.

Algorithm

The implementation is iterative, using a linked stack of pending right-portions to avoid recursion-depth limits that afflict the classical recursive formulation. Starting from the full interval \([a, b]\):

1. Compute \(S_1\) and \(S_2\) on the current active subinterval.
2. If \(|S_1 - S_2| < \varepsilon_{\text{local}}\): accept — add \(S_2\) to the running total, pop the next pending subinterval from the stack, and reuse the already computed endpoint values \(f(r)\) and \(f(m_{\text{prev}})\) as the new left endpoint and midpoint.
3. If \(|S_1 - S_2| \geq \varepsilon_{\text{local}}\): reject — push the right half \([m, r]\) (with its cached function values) onto the stack, and set the active subinterval to the left half \([l, m]\).
4. Repeat from step 1. Terminate successfully when the right endpoint \(b\) is reached; terminate with an error if the active subinterval length falls below min_h.

Each accepted step reuses three of the five function values from the previous iteration (the three points on the new coarser subinterval coincide with points already evaluated); only two new evaluations (the quarter-points \(l+H/4\) and \(r-H/4\)) are needed per step.

Function evaluation budget

The number of evaluations depends on the integrand’s behaviour:

• Smooth integrands: few bisections are needed; the algorithm accepts most subintervals on the first or second attempt, using roughly 2–4 evaluations per accepted panel beyond the initial three.
• Rapidly varying integrands: many bisections concentrate evaluations in the difficult region, while smooth regions are traversed cheaply.
• Worst case: if every step requires bisection down to min_h, the number of evaluations is proportional to \((b-a)/\text{min_h}\).

Limitations

The adaptive Simpson method requires the integrand to be continuous on \([a, b]\). It may fail or produce unreliable results for functions with true singularities or jump discontinuities inside the interval. The error indicator \(|S_1 - S_2|\) is only reliable when the integrand is smooth on the current subinterval; it can be fooled by integrand features that happen to cancel at the five sample points.

If a subinterval is bisected until its length reaches min_h without the tolerance being satisfied, the method returns an error. This typically signals a singularity or extremely rapid variation; try splitting the interval manually at the problem point, reducing min_h, or increasing tolerance.

Romberg’s method

Motivation

Romberg’s method combines the composite trapezoidal rule with Richardson extrapolation to systematically cancel the leading error terms from the Euler-Maclaurin expansion. Each extrapolation step eliminates one error term, so with \(n\) columns the final estimate has accuracy \(O(h^{2n})\) — matching a Gauss rule of order \(2n\) but using only \(2^{n-1}+1\) function evaluations. It is an excellent choice when high precision is required for smooth integrands on finite intervals.

Example

use integrate::romberg::romberg_method;

fn square(x: f64) -> f64 {
    x.powi(2)
}

let a = 0.0;
let b = 1.0;

let num_steps: usize = 10;

let integral = romberg_method(square, a, b, num_steps);
println!("{}", integral);

Understanding Romberg’s method

Euler-Maclaurin expansion

The Euler-Maclaurin summation formula expresses the composite trapezoidal estimate \(T_h(f)\) with step size \(h = (b-a)/n\) as

\[ \begin{align} T_h(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^2}{12} \left[ f^\prime(b) - f^\prime(a) \right] - \frac{h^4}{720} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &+ \cdots + K h^{2p-2} \left[ f^{(2p-3)}(b) - f^{(2p-3)}(a) \right] + O(h^{2p}) \end{align} \]

Two features make this expansion ideal for Richardson extrapolation:

1. The error is a pure even-power series in \(h\): terms \(h^2, h^4, h^6, \ldots\) with coefficients that depend on the endpoints but not on \(h\).
2. When the step size is halved to \(h/2\), each coefficient is divided by the corresponding power of 4 — which allows exact cancellation.

Richardson extrapolation

Halving \(h\) gives a new estimate

\[ \begin{split} T_{h/2}(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^2}{4 \cdot 12} \left[ f^\prime(b) - f^\prime(a) \right] \\ &- \frac{h^4}{16 \cdot 720} \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

Forming the linear combination \(\frac{4 T_{h/2}(f) - T_h(f)}{3}\) cancels the \(h^2\) term exactly:

\[ \begin{split} \frac{4 T_{h/2}(f) - T_h(f)}{3} &= \int_{a}^{b} f(x)\,dx + \frac{h^4}{2880} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- \frac{h^6}{96768} \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{split} \]

This extrapolation step improves the order from \(O(h^2)\) to \(O(h^4)\). Applying it again to two such \(O(h^4)\) estimates cancels the \(h^4\) term and yields \(O(h^6)\), and so on.

The Romberg table

The full process is organised into a triangular table \(R[i, j]\), where \(i\) is the refinement level (step size \(h_i = (b-a)/2^i\)) and \(j\) is the extrapolation column:

\[ R[i, j] = \frac{4^j , R[i, j-1] - R[i-1, j-1]}{4^j - 1}, \qquad 1 \leq j \leq i \]

The first column \(R[i, 0]\) is the composite trapezoidal estimate at step size \(h_i\). Column \(j\) contains estimates that are \(O(h^{2(j+1)})\) accurate. The divisors \(4^j - 1\) take the values \(3, 15, 63, 255, \ldots\) as \(j\) increases.

For \(n\) columns the table has \(n(n+1)/2\) entries and the bottom-right entry \(R[n-1, n-1]\) is the final estimate, with error \(O(h^{2n})\). For example, with \(n = 10\) columns the error behaves as \(O(h^{20})\), which for smooth integrands gives double-precision accuracy with far fewer evaluations than any fixed Newton-Cotes rule of that order.

Computation

Midpoint reuse. Each first-column entry \(R[i, 0]\) is obtained from the previous one by adding only the \(2^{i-1}\) new midpoints — the \(2^{i-1}+1\) points already evaluated at level \(i-1\) lie exactly on the level-\(i\) grid:

\[ \begin{split} R[i, 0] &= \frac{1}{2}, R[i-1, 0] \\ &+ h_i \sum_{k=0}^{2^{i-1}-1} f\left(a + (2k+1),h_i\right) \end{split} \]

This means the total number of function evaluations across all \(n\) levels is \(2^{n-1}+1\), not \(O(n \cdot 2^{n-1})\) as would be required if each level recomputed all its points independently.

Rolling buffer. The Richardson extrapolation depends only on the previous row and the current row of the table. The implementation therefore keeps only two rows of length \(n\) in memory at any time, rather than the full \(n \times n\) table, keeping the memory footprint at \(O(n)\).

Summary of cost for \(n\) columns:

Limitations

Romberg’s method requires the integrand to be infinitely differentiable on \([a, b]\) for the Euler-Maclaurin expansion to hold and the cancellation to be effective. It is not appropriate for functions with singularities, discontinuities, or even kinks (non-differentiable points) in the integration interval — at such points the error expansion breaks down and the extrapolation accelerates to the wrong limit.

The number of columns \(n\) controls both the accuracy and the number of evaluations. For double-precision arithmetic, \(n \approx 10\) is typically sufficient for smooth integrands; using larger \(n\) beyond the point where \(R[n-1, n-1]\) has converged to machine precision adds evaluations without improving the result. For functions with near-singularities or rapid oscillation, prefer the adaptive Simpson method.