The Miracle of Newton-Cotes 
*

*Background and motivation:


Suppose I want to approximate the integral $\int_{a}^{b}f(x)dx$ using evenly-spaced sampled values of $f$: $f(a + \frac{i}{n}(b-a)), i=0,\cdots,n$. By the linearity of the integral, it is natural to try a linear approximation:
$$ \int_{a}^{b}f(x)dx \sim \sum_{i=0}^{n} w_i f(a + \frac{i}{n}(b-a))$$
This is the idea of the Newton-Cotes formulas. The weights $w_i$ are chosen in such a way  that the formula is accurate for all polynomials of degree up to $n$ (and then the error depends on the $(n+1)$'th derivative of $f$). 
In other words, just solve the following $n+1$ equations:
$$ \int_{a}^{b} x^{k} dx = \frac{b^{k+1}-a^{k+1}}{k+1}= \sum_{i=0}^{n} w_i (a + \frac{i}{n}(b-a))^k, k=0,\cdots,n$$
The relevant matrix is a Vandermonde matrix, so there is indeed a solution.


*

*The even miracle and my question:


According to this wikipedia page, it seems that for even $n$, this formula is also accurate for $k=n+1$, and in general the error term depends on the $(n+2)$'th derivative of the $f$. This is a non-trivial linear algebra statement.
So my question is: why is true?
 A: It's not a miracle, just a nice property, and there are much more impressive methods out there, such as Gauss, Clenshaw-Curtis, double exponential rule. See also mathworld.
What happens is easier to see on an interval $[-h,h]$ centered around $0$. If you have an analytic function $f$, then
$$\int_{-h}^h f(x)\,dx = \sum_{m\geq0} \frac{f^{(m)}(0)}{(m+1)!} 2h^{m+1}[m\text{ even}], $$
where the sum involves only terms with even $m$, because $\int_{-h}^h x^m$ vanishes when $m$ is odd.
With an equally spaced choice of points $x_j$, and $w_j = w_{n-j}$ (which happens for the Newton-Cotes rule), the sum
$$ \sum_{j=0}^{n} w_j x_j^m $$
will vanish when $m$ is odd, because we will have $x_j^m = -x_{n-j}^m$. To understand why $w_j = w_{n-j}$, note that
$$ w_j = \int_{-h}^h L_j(x)\,dx, $$
where $\{L_j(x)\}$ are the basis for the Lagrange interpolating polynomial of $f(x)$ on the points $x_j$.
So the rule correctly integrates all odd terms in the power series to zero, leaving only even terms in the error.
