Suppose I estimate $\int_0^1 f(x)dx$ as $\frac12(f(a)+f(b))$, with $a,\,b$ chosen to achieve the lowest-order possible error. We assume $f$ equals its Maclaurin series and $\int$ commutes with $\sum$, viz.$$\int_0^1f(x)dx=\sum_{n\ge0}\frac{f^{(n)}(0)}{(n+1)!},\,\frac12(f(a)+f(b))=\sum_{n\ge0}\frac{f^{(n)}(0)(a^n+b^n)}{2\cdot n!}.$$We thus want as many small $n$ as possible to satisfy$$a^n+b^n=\frac{2}{n+1}.$$The case $n=0$ is trivial; the case $n=1$ is $a+b=1$, which we expected; the case $n=2$ is $a^2+b^2=\frac23$, and we can now show $a,\,b$ must be $\frac12\pm\frac{1}{\sqrt{12}}$. Beautifully, the case $n=3$ is now also correct, because$$a^3+b^3=\tfrac12(a+b)(3(a^2+b^2)-(a+b)^2)=\frac12\cdot 1\cdot(3\cdot\tfrac23-1^2)=\frac12,$$which is what we needed. Unfortunately, that's as far as it works, since$$a,\,b=\frac12\pm\frac{1}{\sqrt{12}}\implies ab=\frac13\implies a^4+b^4=(a^2+b^2)^2-2(ab)^2=\frac49-\frac29=\frac29\ne\frac25.$$But the success up to $n=3$ is a big deal, since any other $a,\,b$ fail at least two terms earlier. In particular, this buys us an extra $O(h^2)$ factor in the error term when we cut an integral into $N\gg1$ strips of width $h\propto 1/N$.
None of this is original to me. I've seen it before, but I can't remember what using these values of $a,\,b$ is called. I've checked a number of algorithms with Wikipedia articles, but none of them are this one.
