For the closed Newton-Cotes quadrature over $[x_1, x_n]$, the coefficients $H_{n,i}$ for$$\int_{x_1}^{x_n} f(x)\:\text{d}x = h \sum_{i=1}^n H_{n,i} \; f(x_i)$$are given explicitly by$$H_{n,r+1} =\frac{(-1)^{n-r}}{r!(n-r)!}\int_0^n \frac{\prod_{k=0}^n (t-k)}{t-r}\:\text{d}t;$$see http://mathworld.wolfram.com/Newton-CotesFormulas.html.To compute those values numerically, one could go ahead and evaluate the integral, but perhaps that's not the most efficient thing to do. SciPy does something else, but I don't quite get it.Any hints?
computing Newton-Cotes weights
quadrature
I can explain what the scipy code is doing (https://github.com/scipy/scipy/blob/v0.18.0/scipy/integrate/quadrature.py#L833-L842). I haven't checked this directly, but this should be correct in principle.The line ti = 2 * yi - 1 remaps the interval $(0,1)$ to $(-1,1)$, so now the weights will be the integrals of Lagrange elementary interpolants over $(-1,1)$ with nodes $(-1:2/N:1)$.The line C = ti ** nvec[:, np.newaxis] constructs the (transpose of) the Vandermonde matrix, with entries $C_{ij} = t_i ^ j$. This matrix has the property that when multiplied on the left by a vector $\alpha$, it produces$$\alpha^\top C = (p_\alpha(t_0), \ldots, p_\alpha(t_N))$$which is the vector of $p_\alpha$ (polynomial with coefficients being the entries of $\alpha$) evaluated at each $t_0,\ldots,t_N$.The inverse of this, $C^{-1}$, will then consist of vectors $\alpha_k$ that give the coefficients of Lagrange interpolants, because $\alpha_k^\top C = e_k$, and having $p_{\alpha_k}(t_l) = \delta_{kl}$ is the definition of Lagrange interpolants.The Vandermonde matrix is very ill-conditioned, so scipy applies twice the Newton-Raphson iteration for the (matrix) reciprocal: Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv).Once the coefficients of Lagrange interpolants, $C^{-1}$, are known, it is straightforward to integrate them, by multiplying $C^{-1}$ on the right with $(2/1, 0, 2/3, 0, 2/5, \ldots)^\top$ (the integrals of $x^m$ for odd $m$ are zerothe integral is over $(-1,1)$ because of ti above), which gives the Newton-Cotes weights: ai = Cinv[:, ::2].dot(vec) * (N / 2.)