Simpson's rule (${\int }_a^{a+2h}f\approx \frac{h}{3}[f(a)+4f(a+h)+f(a+2h)]$)

Simpson's rule approximates a definite integral by integrating the quadratic polynomial obtained by polynomial interpolation at two endpoints and their midpoint. For an interval $[a,b]$ with $a<b$, write $m=(a+b)/2$ and $h=(b-a)/2$. The three quadrature nodes are $a,m,b$, and the quadrature rule isThus its quadrature weights are $h/3,4h/3,h/3$. Their sum is $2h=b-a$, as required for exactness on constant functions. This is the three-node Newton-Cotes closed quadrature; its degree of precision of a quadrature rule is three, even though the interpolating polynomial has polynomial degree at most two.^[ref]

For an original worked example, take the cubic polynomial $f(x)=x^3-2x+4$ on the interval $[-1,2]$. Here $m=1/2$ and $h=3/2$, so the function evaluations giveSimpson's rule givesAn antiderivative supplies an independent exact definite integral:The quadratic polynomial through the same quadrature nodes is $p_2(x)=\frac{3}{2}{x}^2-\frac{1}{2}x+3$. It differs from $f$ between the quadrature nodes, but its definite integral agrees exactly. The picture shows why: the signed polynomial interpolation error has two lobes with opposite integrals.

Three function evaluations therefore capture more information about an integral than about the whole function. When an antiderivative is unavailable, composite Simpson's rule repeats this inexpensive local construction over a partition of an interval. For a sufficiently smooth function, halving the node spacing reduces its truncation error by approximately a factor of sixteen.

Polynomial interpolation makes Simpson's rule exact for every polynomial of polynomial degree at most two. To understand the extra polynomial degree, centre the coordinate at the midpoint: $y=x-m$, with endpoints $-h,h$. For any cubic polynomial $p$, let $p_2$ be its interpolating quadratic polynomial. The polynomial $p-p_2$ vanishes at all three quadrature nodes, sowhere $A$ is the coefficient of $x^3$ in $p$. As a function of $y$, this is an odd function, whose integral over the symmetric interval $[-h,h]$ is zero. Hence integrating $p_2$ also integrates $p$ exactly. This is the symmetry mechanism in the nodal-polynomial criterion for quadrature exactness.

In the worked example, the polynomial difference factors asIts definite integrals on $[-1,1/2]$ and $[1/2,2]$ are respectively $81/64$ and $-81/64$. Their cancellation concerns signed integrals, rather than pointwise accuracy of polynomial interpolation.

The degree of precision of a quadrature rule stops at three: on the same interval, the monomial $x^4$ gives $S(x^4)=69/8$, whereas ${\int }_{-1}^2{x}^{4}dx=33/5$. Their difference, exact integral minus approximation, is $-81/40$.

Composite Simpson's rule applies Simpson's rule to consecutive pairs of equal subintervals. Choose a positive even integer $n\ge 2$, put $h=(b-a)/n$, and take quadrature nodes $x_j=a+jh$ for $j=0,\dots ,n$. There are $n$ subintervals, $n/2$ two-subinterval panels, and $n+1$ quadrature nodes. ThenThe second sum is empty when $n=2$, recovering one panel. The quadrature weights, before multiplication by $h/3$, are $1,4,2,4,\dots ,2,4,1$. Every interior even-indexed quadrature node is an endpoint of two panels, so its two endpoint contributions add to two; every odd-indexed quadrature node is a midpoint with weight four. An odd $n$ cannot be covered by these pairs. The displayed formula requires equal spacing.^[ref]

For the exponential function on the interval $[0,1]$, the exact integral is $I=e-1$. With $n=2$, the formula gives $S_2\approx 1.718861151876593$, an overestimate of about $5.793234175\times 10^{-4}$. With $n=128$, it gives $S_{128}\approx 1.718281828494607$, an overestimate of about $3.556140064\times 10^{-11}$. The figure shows the intermediate refinements and compares every computed absolute error with the error bound for composite Simpson's rule.

For this exponential function, a geometric series also gives a compact independent formula. With $q=e^h$ and $h=1/n$, summing the $n/2$ panel contributions yieldsA Taylor series at $h=0$ givesThe Big O notation describes the remaining terms as $h$ tends to zero. Consequently $(S_n-I)/(S_{2n}-I)$ tends to sixteen: the computed ratios increase from about $15.65169$ to $15.99965$. Fourth-order truncation error is also obtained by Richardson extrapolation of the composite trapezoidal rule.^[ref]

Suppose $f$ has a fourth derivative that is a continuous function on the interval $[a,b]$, and $|f^{(4)}(x)|\le M$ there. With $I={\int }_a^{b}f(x)dx$, the signed quadrature error here is $I-S_n$. Then, for some point $\xi \in (a,b)$,In particular, a strictly positive fourth derivative makes $I-S_n<0$: composite Simpson's rule overestimates the integral.^[ref]

The bound's scale follows by adding the panel errors. A panel of width $2h$ has signed quadrature error $-h^5{f}^{(4)}({\xi }_j)/90$ and absolute error at most $M{h}^5/90$. There are $(b-a)/(2h)$ panels, so the triangle inequality gives $M(b-a)h^4/180$. The intermediate value theorem for the continuous fourth derivative replaces the average of the panel values $f^{(4)}({\xi }_j)$ by a single value $f^{(4)}(\xi )$, giving the signed formula.

The underlying single-panel calculation can be seen through the Peano kernel theorem. On the interval $[0,2]$, define the raw Peano kernel $K(t)=L[(x-t{)}_{+}^3]$, where $L$ is integral minus Simpson's rule and $(x-t{)}_{+}=\max (x-t,0)$ is the positive part. Then $K(t)=-t^3(4-3t)/12$ for $0\le t\le 1$ and $K(t)=-(2-t{)}^3(3t-2)/12$ for $1\le t\le 2$. It is strictly negative inside the interval and has integral $-1/15$. ThusThe uniform norm $\Vert f^{(4)}{\Vert }_{\infty }$ bounds the fourth derivative. The polynomial $f(x)=x^4/24$ attains equality. Scaling to endpoint spacing $h$ multiplies the coefficient by $h^5$, giving the sharp single-panel bound $h^5\Vert f^{(4)}{\Vert }_{\infty }/90$. The normalized Peano kernel is $K/3!$, so the external factorial must be omitted when using that normalization.

For the exponential function on $[0,1]$, $f^{(4)}(x)=e^x$, so $M=e$ and the bound is $e/(180n^4)$. Its leading actual absolute error is $(e-1)/(180n^4)$, smaller than the bound because a uniform upper bound on the fourth derivative does not use its variation across the interval. The theorem bounds truncation error with exact function evaluations; sample errors and roundoff error are additional contributions, as described by noise sensitivity of a quadrature rule.