Let ${\displaystyle f}$ be a continuous function on ${\displaystyle [a,b]}$. Let ${\displaystyle n}$ be an integer, ${\displaystyle h=(b-a)/n}$, ${\displaystyle x_{i}=a+ih}$, ${\displaystyle f_{i}=f(x_{i})}$. Then we can estimate ${\displaystyle \int _{a}^{b}f(x)\,dx}$ as the linear combination ${\displaystyle \sum c_{i}f_{i}}$ for some choice of coefficients ${\displaystyle c_{i}}$.

In the case ${\displaystyle n=6}$ we have, for example ${\displaystyle \sum c_{i}f_{i}={\frac {h}{2}}(f_{0}+2f_{1}+2f_{2}+2f_{3}+2f_{4}+2f_{5}+f_{6}){\mbox{ (trapezoidal rule)}}}$ ${\displaystyle \sum c_{i}f_{i}={\frac {h}{3}}(f_{0}+4f_{1}+2f_{2}+4f_{3}+2f_{4}+4f_{5}+f_{6}){\mbox{ (Simpson's rule)}}}$ Our aim is to explain where the magic coefficients come from.

Suppose the points ${\displaystyle \{x_{i}\}}$ are divided into ${\displaystyle p}$ contiguous panes, each with ${\displaystyle (m+1)}$ points, so that ${\displaystyle n=pm}$. For the trapezoidal rule, we have ${\displaystyle p=6}$, ${\displaystyle m=1}$, and for Simpson's rule we have ${\displaystyle p=3}$, ${\displaystyle m=2}$. Consider a single pane with ${\displaystyle (m+1)}$ points, and suppose ${\displaystyle a=0}$, ${\displaystyle b=m}$, ${\displaystyle h=1}$. The general result will follow from the linearity properties of the integral. Let ${\displaystyle L_{i}(x)}$ be the Lagrange basis polynomials of degree ${\displaystyle m}$ satisfying ${\displaystyle L_{i}(j)=\delta _{ij}}$ for ${\displaystyle 0\leq i,j\leq m}$. Then ${\displaystyle f(x)\simeq \sum L_{i}(x)f_{i}}$ ${\displaystyle \int _{a}^{b}f(x)\,dx\simeq \sum \left(\int _{a}^{b}L_{i}(x)\,dx\right)f_{i}=\sum c_{i}f_{i}}$

require 'numeric'
V=:^/~ @: i.@:#
L=:clean@(%. V)"1
basis=:L @: =@i.@>:
int=:(0&p.. p. <:@:#)"1
c=:int@basis
c 1
0.5 0.5
c 2
0.333333 1.33333 0.333333


For a single pane, these give

${\displaystyle \int _{x_{0}}^{x_{1}}f(x)\,dx\simeq {\frac {h}{2}}(f_{0}+f_{1}){\mbox{ (trapezoidal rule)}}}$ ${\displaystyle \int _{x_{0}}^{x_{2}}f(x)\,dx\simeq {\frac {h}{3}}(f_{0}+4f_{1}+f_{2}){\mbox{ (Simpson's rule)}}}$ The result for multiple panes is simple additivity of the single pane result:

tz=:#~ -.@(*./\.)@(0&=)  NB. remove trailing zeros
ppr=:+//.@(*/)           NB. polynomial product
pprtz=:tz@:ppr
coeff=:(c@:]) pprtz 0 = ] | i.@:*
6 coeff 1
0.5 1 1 1 1 1 0.5
3 coeff 2
0.333333 1.33333 0.666667 1.33333 0.666667 1.33333 0.333333


The same method works for generating formulas with more interpolating points (although there are good reasons for stopping at Simpson's rule). If we take ${\displaystyle m=3}$, we get:

   2 coeff 3
0.375 1.125 1.125 0.75 1.125 1.125 0.375


${\displaystyle \sum c_{i}f_{i}={\frac {3h}{8}}(f_{0}+3f_{1}+3f_{2}+2f_{3}+3f_{4}+3f_{5}+f_{6})\ \ (Simpson's\ {\frac {3}{8}}\ rule)}$