computing Newton-Cotes weights

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?

• Is your question why SciPy does something different? Are they just using a different equation that, however, evaluates to the same number? – Wolfgang Bangerth Oct 28 '16 at 15:22
• They are the same numbers, but I don't understand what they are doing (any why). No reference given in the source code. – Nico Schlömer Oct 28 '16 at 15:25

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 zero—the integral is over $(-1,1)$ because of ti above), which gives the Newton-Cotes weights: ai = Cinv[:, ::2].dot(vec) * (N / 2.)

• You should propose a pull request that adds your explanation to the documentation of that function. – Wolfgang Bangerth Oct 28 '16 at 20:37

In the pass I did something similarly. At the time I applied Gauss-Legendre quadrature to evaluate $H_{n,r+1}$ intergals; this calculation adopted (1) Gauss-Legendre weight and points calculated in advance (then, change the variable $t$ in the integrand to $\xi \in [-1,1]$), and (2) multi-precision floating point library fmlib.

I guess that you obtain $H_{n,i}$ from Lagrange Interpolation Formulae. Let me review some basic thing so that you may recognize the same derivation. From the Lagrange Interpolation, we have \begin{align} f(x) = \sum_{i=1}^{n} l_i (x) f_i + R_n (x) \end{align} where \begin{align} l_i (x) &= \frac{ \Pi_n( x) }{ (x-x_i) } . \frac{ 1 }{ \Pi _n ' (x_i) } = \frac{ ( x - x_1 ) \ldots ( x - x_{i-1} ) ( x - x_{i+1} ) \ldots ( x - x_{n} ) }{ ( x_i - x_1 ) \ldots ( x_i - x_{i-1} ) ( x_i - x_{i+1} ) \ldots ( x_i - x_{n} ) } , \\ R_n (x) &= \frac{ f^{(n)} (\xi) }{ n! } \cdot \Pi _n (x) , \end{align} for some $\xi = \xi(x) \in (x_1, x_n )$. Note that we defined $\Pi_n (x) = \prod_{i=1}^{n} (x-x_i)$, and we also have \begin{align} \left| R_n (x) \right| &\le \frac{ ( x_n - x_1 )^{n} }{ n! } \max_{ x_1 \le x \le x_n } \left| f^{(n)} (x) \right| \end{align} which can be proved using the Role’s theorem.

Now we take \begin{align} \int_{x_1}^{x_n} f(x) d x &= \sum_{i=1}^{n} \underbrace{ \left( \int_{x_1}^{x_n} l_i (x) d x \right) }_{ A_ i } f_i + \underbrace{ \int_{x_1}^{x_n} R_n (x) d x }_{ B } \end{align} Define \begin{align} A_ i &= \int_{x_1}^{x_n} l_i (x) d x ~~~(\equiv H_{n,i} ? ) \\ B &= \int_{x_1}^{x_n} R_n (x) d x \end{align} We are going to estimate those integrals. Change variable, \begin{align} h &= \frac{ x_n - x_1 }{ n - 1 } ,~~~~~ x = x_1 + (t-1) h ,~~ \forall t \in [1,n] \end{align} Using the latter one, $x - x_j = x_1 + (t-1)h - x_j$ and $x_j = x_1 + (j-1) h$, we have $$x - x_j = x_1 + (t-1)h - \left[ x_1 + (j-1)h \right] = (t-j)h$$ Thus \begin{align} A_ i &= \int_{x_1}^{x_n} l_i (x) d x = \int_{x_1}^{x_n} \frac{ ( x - x_1 ) \ldots ( x - x_{i-1} ) ( x - x_{i+1} ) \ldots ( x - x_{n} ) }{ ( x_i - x_1 ) \ldots ( x_i - x_{i-1} ) ( x_i - x_{i+1} ) \ldots ( x_i - x_{n} ) } d x \\ &= \int_{1}^{n} \frac{ h ( t - 1 ) \cdot h ( t - 2 ) \cdots h ( t-i+1 ) \cdot h (t-i-1) \cdots h(t-n) }{ h ( i - 1 ) \cdot h ( i - 2 ) \cdots h ( i-i+1 ) \cdot h (i-i-1) \cdots h(i-n) } \cdot h d t \end{align} Or \begin{align} A_ i &= h \int_{1}^{n} \frac{ ( t - 1 ) \cdot ( t - 2 ) \cdots ( t-i+1 ) \cdot (t-i-1) \cdots (t-n) }{ ( i - 1 ) \cdot ( i - 2 ) \cdots ( i-i+1 ) \cdot (i-i-1) \cdots (i-n) } d t \end{align} And \begin{align} B &= \int_{x_1}^{x_n} R_n (x) d x = \int_{x_1}^{x_n} ( x - x_1 ) \ldots ( x - x_{i} ) \ldots ( x - x_{n} ) \cdot \frac{ f^{(n)} (\xi (x) ) }{ n! } d x \\ &= \int_{1}^{n} h ( t - 1 ) \cdot h ( t - 2 ) \cdots h ( t - i) \cdots h(t-n) \cdot \frac{ f^{(n)} (\xi (x(t)) ) }{ n! } \cdot h d t \end{align} We estimate \begin{align} | B | &\le M h^{n+1} \end{align} where $M = |C_n | \sup_{ \xi \in [x_1,x_n] } \left| f^{(n)} (\xi) \right|$ and \begin{align} C_n &= \frac{1}{ n! } \int_{1}^{n} ( t - 1 ) \cdots ( t - i) \cdots ( t - n ) d t \label{c} \end{align} Here integrals $A_i$ and $C_n$ can be estimated EXACTLY using Gauss-Legendre formulas with $n/2$ points.

As I could see in numerical result at the time, we have (a) $A_i$ are increased to very large values as $n$ increasing, for example $\max A_i \sim 10^{0}, 10^{2}, 10^{4}, 10^{7}, 10^{10}$ and $10^{24}$ when $n=10,20,30,40,50$, and $100$, respectively; (b) $A_i$ is symmetry with repect to $i$, i.e. $A_{i} = A_{n-i+1}$; (c) $\sum_{i=1}^{n} A_{i} = n-1$; and (c) $C_n = 0$ for $n$ odd, and $C_n$ is bounded by a small value when $n$ even.