Composite Simpson's rule — the same algorithm as scipy.integrate.simpson (SciPy ≥ 1.11).
x is either sample positions (same length as y, strictly monotonic — decreasing integrates with the
sign of the direction, as in SciPy) or a uniform spacing dx (default 1).
Odd number of samples (even number of intervals): every pair of intervals [x₀, x₁, x₂] is integrated
by the parabola through its three points, with the weights for unequal widths h₀, h₁:
(h₀+h₁)/6 · [ (2 − h₁/h₀)·y₀ + (h₀+h₁)²/(h₀h₁)·y₁ + (2 − h₀/h₁)·y₂ ]. Exact for cubics on any grid.
Even number of samples: Simpson over the first N − 1 samples plus Cartwright's correction for the last
interval (the parabola through the last three points integrated over the last interval only):
α·y_{N−1} + β·y_{N−2} − η·y_{N−3} with α = (2h₁² + 3h₀h₁)/(6(h₀+h₁)), β = (h₁² + 3h₀h₁)/(6h₀),
η = h₁³/(6h₀(h₀+h₁)). Exact for quadratics.
Two samples: the trapezoid (the only rule two points define). Fewer: 0.
Composite Simpson's rule — the same algorithm as
scipy.integrate.simpson(SciPy ≥ 1.11).xis either sample positions (same length asy, strictly monotonic — decreasing integrates with the sign of the direction, as in SciPy) or a uniform spacingdx(default 1).[x₀, x₁, x₂]is integrated by the parabola through its three points, with the weights for unequal widthsh₀, h₁:(h₀+h₁)/6 · [ (2 − h₁/h₀)·y₀ + (h₀+h₁)²/(h₀h₁)·y₁ + (2 − h₀/h₁)·y₂ ]. Exact for cubics on any grid.α·y_{N−1} + β·y_{N−2} − η·y_{N−3}withα = (2h₁² + 3h₀h₁)/(6(h₀+h₁)),β = (h₁² + 3h₀h₁)/(6h₀),η = h₁³/(6h₀(h₀+h₁)). Exact for quadratics.